<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Modeling Gene Networks Using Ordinary Differential Equations
Author
Step2:  Visualizing the activating Hill function
Now we'll setup a plot to visualize what this function looks like.
Step3:  <h2> <font color='firebrick'> Question Set 1 (3 pts) </font> </h2>
What happens to the shape of the Hill function curve when you vary $n$ between 1 and 8?
What happens to the curve when you vary $\beta$ between 2 and 20?
What happens to the curve when you vary $K$ between 2 and 20?
Modeling transcriptional repression
If rather than stimulating the production of $Y$, $X$ "represses" $Y$, we can write the corresponding Hill function as
Step4:  <h2> <font color='firebrick'> Question Set 2 (3 pts) </font> </h2>
As before change the values of the simulation to answer these questions
What happens to the shape of the Hill repressive function curve when you vary $n$ between 1 and 8?
What happens to the curve when you vary $B$ between 2 and 20?
What happens to the curve when you vary $K$ between 2 and 20?
Simplifying Models using Logic Approximations
To simplify analysis it's often convenient to approximate step-like sigmoidal functions like those produced by the Hill equation with functions using logic approximations. 
We'll assume that when the transcription factor, $X$, is above a threshold, $K$, then  gene $Y$ is transcribed at a rate, $\beta$.  When $X$ is below the threshold, $K$, gene $Y$ is not be transcribed.  To represent this situation,  we can rewrite the formula for $Y$ as
Step5:  And now we can generate some plot to compare the logic approximation to the Hill function, for the activating case
Step6:  <h2> <font color='firebrick'> Question Set 3 (3 pts) </font> </h2>
As before change the values of the simulation to answer these questions
How does the parameter K in the Hill function relate to the location of the "step" of the logic approximation?
How does the parameter $\beta$ relate to the rate of production of $Y$?
For waht values of $n$ is the logic approximation most like the hill function?
Logic approximations for multi-input functions
What if a gene needs two or more activator proteins to be transcribed?  We can describe the amount of $Z$ transcribed as a function of active forms of $X$ and $Y$ with a function like
Step7:  <h2> <font color='firebrick'> Question Set 4 (3 pts) </font> </h2>
The concentration of $Y$ eventually reaches a steady state, $Y_{st}$.  How does $Y_{st}$ relate to $\beta$ and $\alpha$?
The response time of a dynamical system, $T_{1/2}$is defined as the time it takes for it to go half-way between it's initial and final value.  
a) How does the response time change as you vary $\beta$?
b) How does the response time change as you vary $\alpha$?
Estimate the response time as you vary the parameter  $\alpha$ and see if you can create a plot (in Python, or R, or Excel) showing the relationship between $\alpha$ and response time, for $0.01 \leq \alpha \leq 0.1$.
Toggling the activator X
In the proceeding example the activator $X$ was on at the beginning of the simulation and just stayed on.  Let's see what happens when $X$ has pulsatile dynamics.  This would be akin to  toggling $X$ on then off, and asking what happens to $Y$.
Step8:  Feed Forward Loops
We're now going to use some of these tools to look at a class of network motifs (small network topologies), called Feed Forward Loops (FFLs), found in signaling and regulatory networks. FFLs involve interactions between three components, with the basic topology illustrated below. Depending on the signs of the edges (whether activating or repressing) we can classify FFLs as "coherent" or "incoherent." We'll take a look at an example of each class.
A Coherent FFL
The most common type of coherent FFL is illustrated in the figure below.  In this system $X$ is an activator of $Y$ and both $X$ and $Y$ regulate the production of $Z$ with AND logic (i.e. both $X$ and $Y$ must be above particular thresholds in order to trigger the production of $Z$). 
Using our logic approximation framework we will model the coherent FFL network illustrated above as follows.
Gene Y
Step9:  Define Python functions for dY/dt and dZ/dt
Recall from above that
\begin{eqnarray}
\frac{dY}{dt} & = & \beta_y\ \Theta(X > K_{xy}) - \alpha_{y}Y \ \
\frac{dZ}{dt} & = & \beta_z\ \Theta(X > K_{xz})\Theta(Y > K_{yz}) - \alpha_{z}Z
\end{eqnarray}
Step10:  <h2> <font color='firebrick'>Question Set 5 (2 pts)</font> </h2>
How do the dynamics of $Y$ and $Z$ differ in the simulation above?
Try varying the length of the first (short) pulse.  How does changing the length of the pulse affect the dynamics of $Y$ and $Z$?
Performance of the Coherent FFL under noisy inputs
Let's further explore the behavior of the coherent FFL defined given noisy inputs.  As before we're going to define an input signal, $X$, that has a short and long pulse, but now we're going to pollute $X$ with random noise.
Step11:  <h2> <font color='firebrick'>Question Set 6 (2 pts) </font> </h2>
In the code cell above, try changing the duration of the first pulse and the scale of the noise (see comments in code) to get a sense of how good a filter the FFL is.  Is there a bias to the filtering with respect to Z turning on versus Z turning off?
Dynamics of Y and Z in the Coherent FFL
As before we can solve for Y as a function of time and calculate what its steady state value will be
Step12:  Type 1 Coherent FFLs can act as a Sign-Sensitive Delays
As discussed in the article by Shen-Orr et al. a feed forward loop of the type we've just discussed can act as a type of filter -- a sign-sensitive delay that keeps $Z$ from firing in response to transient noisy signals from $X$, but shuts down $Z$ immediately once the signal from $X$ is removed. 
An Incoherent FFL
Consider the FFL illustrated in the figure below.  
In this incoherent FFL, the logic function that regulates $Z$ is "X AND NOT Y".  That is $Z$ turns on once $X$ is above a given threshold, but only stays on fully as long as $Y$ is below another threshold. Again for simplicity we assume $K_{xy} = K_{yz}$. 
Dynamics of Y
As before, the dynamics of $Y$ are described by
<ASSISTANT_TASK:>
Python Code:
# import statements to make numeric and plotting functions available
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
def hill_activating(X, B, K, n):
     Hill function for an activator
    return (B * X**n)/(K**n + X**n)
Explanation: Modeling Gene Networks Using Ordinary Differential Equations
Author: Paul M. Magwene 
<br> 
Date: April 2016
Background readings
Shen-Orr, S. S., et al. 2002. Network motifs in the transcriptional regulation network of Escherichia coli. Nat Genet 31(1):64-8. http://dx.doi.org/10.1038/ng881
Alon, U. 2007. Network motifs: theory and experimental approaches. Nat Rev Genet, 8(6):450–61. http://dx.doi.org/10.1038/nrg2102.
To gain some intuition for how systems biologists build mathematical models of gene networks we're going to use computer simulations to explore the dynamical behavior of simple transcriptional networks.
In each of our simulations we will keep track of the the concentration of a different genes of interest as they change over time.  The basic approach we will use to calculate changes in the quantity of different molecules are  differential equations, which are simply a way of describing the instanteous change in a quantity of interest.  
All of our differential equations will be of this form:
\begin{eqnarray}
\frac{dY}{dt} = \mbox{rate of production} - \mbox{rate of decay}
\end{eqnarray}
To state this in words --  the amount of gene $Y$  changes over time is a function of two things: 1) a growth term which represents the rate at which the gene is being transcribed and translated; and 2) a decay term which gives the rate at which $Y$ trascsripts and protein are being degraded.  
In general we will assume that the "rate of production" is a function of the concentration of the genes that regulate $Y$(i.e. it's inputs in the transcriptional network), while the "rate of decay" is a proportional to the amount of $Y$ that is present. So the above formula will take the following structure:
$$
\frac{dY}{dt} =  f(X_1, X_2, \ldots) - \alpha Y
$$
The $f(X_1, X_2, \ldots)$ term represents the growth term and is a function of the transcription factors that regulate $Y$. The term, $\alpha Y$ represents the rate at which $Y$ is being broken down or diluted.  Notice that the decay rate is a proportional to the amount of $Y$ that is present. 
If $\frac{dy}{dt}$ is positive than the concentration of gene $Y$ is increasing, if $\frac{dy}{dt}$ is negative the concentration of $Y$ is decreasing, and if $\frac{dy}{dt} = 0$ than $Y$ is at steady state.
Modeling the rate of production term with the Hill Function
An appropriate approach for modeling the rate of production of a protein, $Y$, as a function of it's inputs, $X_1, X_2,..$, is a with the "Hill Function".  The Hill Function for a single transcriptional activator is:
$$
f(X) = \frac{\beta X^n}{K^n + X^n}
$$
$X$ represents the concentration of a transcriptional activator and $f(X)$  represents the the combined transcription and translation of the gene $Y$ that is regulated by $X$.
Modeling transcriptional activation
Here a Python function to represent transcriptional activation based on the Hill function given above:
End of explanation
## generate a plot using the hill_activating function defined above
# setup paramters for our simulation
# CHANGE THESE VALUES TO EXPLORE THE EFFECT OF THESE PARAMTERS
# (see questions below)
n = 1
B = 5
K = 10
# generate a range of x values, representing a range of concentrations of our
# transcription factor X
x = linspace(0,30,200)  # generate 200 evenly spaced points between 0 and 30
# calculating corresponding rates of production of Y
y = hill_activating(x, B, K, n) 
# plot the hill fxn with the user set parameters
plot(x, y, label='B = {}, K = {}, n={}\n(user specified)'.format(B, K, n))
# plot the hill fxn with a set of reference parameters to facilitate comparison
plot(x, hill_activating(x, 5, 10, 1), label='B = 5, K = 10, n=1\n(reference)', alpha=0.75)
xlabel('Concentration of X')
ylabel('Rate of production of Y')
legend(loc='best')
pass # suppress further output
Explanation: Visualizing the activating Hill function
Now we'll setup a plot to visualize what this function looks like.
End of explanation
## Python implementation of repressive Hill function 
def hill_repressing(X, B, K, n):
    return B/(1.0 + (X/K)**n)
## generate a plot using the hill_repressing function defined above
# CHANGE THESE VALUES TO EXPLORE THE EFFECT OF THESE PARAMTERS
# (see questions below)
n = 1
B = 5
K = 10
# generate a range of x values, representing a range of concentrations of our
# transcription factor X
x = linspace(0,30,200)  # generate 200 evenly spaced points between 0 and 30
# calculating corresponding rates of production of Y
y = hill_repressing(x, B, K, n) 
# plot the hill fxn with the user set parameters
plot(x, y, label='B = {}, K = {}, n={}\n(user specified)'.format(B, K, n))
# plot the hill fxn with a set of reference parameters to facilitate comparison
plot(x, hill_repressing(x, 5, 10, 1), label='B = 5, K = 10, n=1\n(reference)')
xlabel('Concentration of X')
ylabel('Rate of production of Y')
legend(loc='best')
pass # suppress further output
Explanation: <h2> <font color='firebrick'> Question Set 1 (3 pts) </font> </h2>
What happens to the shape of the Hill function curve when you vary $n$ between 1 and 8?
What happens to the curve when you vary $\beta$ between 2 and 20?
What happens to the curve when you vary $K$ between 2 and 20?
Modeling transcriptional repression
If rather than stimulating the production of $Y$, $X$ "represses" $Y$, we can write the corresponding Hill function as:
$$
f(X) = \frac{\beta}{1 + (X/K)^n}
$$
Remember that both of these Hill functions (activating and repressing)  describe the production of $Y$ as a function of the levels of $X$, not the temporal dynamics of $Y$ which we'll look at after developing a few more ideas.
End of explanation
def logic_activating(X, B, K):
    if X > K:
        theta = 1
    else:
        theta = 0
    return B*theta 
def logic_repressing(X, B, K):
    if X < K:
        theta = 1
    else:
        theta = 0
    return B*theta
Explanation: <h2> <font color='firebrick'> Question Set 2 (3 pts) </font> </h2>
As before change the values of the simulation to answer these questions
What happens to the shape of the Hill repressive function curve when you vary $n$ between 1 and 8?
What happens to the curve when you vary $B$ between 2 and 20?
What happens to the curve when you vary $K$ between 2 and 20?
Simplifying Models using Logic Approximations
To simplify analysis it's often convenient to approximate step-like sigmoidal functions like those produced by the Hill equation with functions using logic approximations. 
We'll assume that when the transcription factor, $X$, is above a threshold, $K$, then  gene $Y$ is transcribed at a rate, $\beta$.  When $X$ is below the threshold, $K$, gene $Y$ is not be transcribed.  To represent this situation,  we can rewrite the formula for $Y$ as:
$$
f(X) = \beta\ \Theta(X > K)
$$
where the function $\Theta = 0$ if the statement inside the parentheses is false and $\Theta = 1$ if the statement is true. An alternate way to write this is:
$$
f(X) = 
\begin{cases}
    0,      &\text{if $X > K$;} \
    \beta,  &\text{otherwise.}
\end{cases}
$$
When $X$ is a repressor we can write:
$$
f(X) = \beta\ \Theta(X < K)
$$
Python functions for the logic approximation
We can write Python functions to represent the logic approximations for activation and repression as follows:
End of explanation
## generate plots using your hill_activating and logic_activating functions defined above
## For X values range from 0 to 30
# CHANGE THESE VALUES TO EXPLORE THE EFFECT OF THESE PARAMTERS
n = 4
B = 5
K = 10
x = linspace(0, 30, 200)
plot(x, hill_activating(x, B, K, n), 
     label='B = {}, K = {}, n={}'.format(B, K, n))
logicx = [logic_activating(i, B, K) for i in x]
plot(x, logicx, label='logic approximation')
xlabel('Concentration of X')
ylabel('Rate of production of Y')
legend(loc='best')
ylim(-0.5, B*1.1)
pass
Explanation: And now we can generate some plot to compare the logic approximation to the Hill function, for the activating case:
End of explanation
## write a function to represent the simple differential equation above
def dYdt(B, K, a, X, Y):
    production = logic_activating(X, B, K)
    decay = a*Y
    return production - decay
## generate a plot of conc of Y over time using your dY function defined above
## Evaluated over 200 time units
B = 5
K = 10
X = K + 1
Y = [0] # initial value of Y
a = 0.05
nsteps = 200
for i in range(nsteps):
    yprevious = Y[-1]
    deltay = dYdt(B, K, a, X, yprevious)
    ynew = Y[-1] + deltay
    Y.append(ynew)
    
plot(Y)
xlabel('Time units')
ylabel('Concentration of Y')
ylim(0, (B/a)*1.1)
pass
Explanation: <h2> <font color='firebrick'> Question Set 3 (3 pts) </font> </h2>
As before change the values of the simulation to answer these questions
How does the parameter K in the Hill function relate to the location of the "step" of the logic approximation?
How does the parameter $\beta$ relate to the rate of production of $Y$?
For waht values of $n$ is the logic approximation most like the hill function?
Logic approximations for multi-input functions
What if a gene needs two or more activator proteins to be transcribed?  We can describe the amount of $Z$ transcribed as a function of active forms of $X$ and $Y$ with a function like:
$$
    f(X,Y) = \beta\ \Theta(X > K_x \land Y > K_y)
$$
The above equation describes "AND" logic (i.e. both X and Y have to be above their threshold levels, $K_x$ and $K_y$, for Z to be transcribed). In a similar manner we can define "OR" logic:
$$
f(X,Y) = \beta\ \Theta(X > K_x \lor Y > K_y)
$$
A SUM function would be defined like this:
$$
f(X,Y) = \beta_x \Theta(X > K_x) + \beta_y \Theta (Y > K_y)
$$
Modeling changes in network components over time
Up until this point we've been considering how the rate of production of a protein $Y$ changes with the concentration of a transcriptional activator/repressor that regulates $Y$.  Now we want to turn to the question of how the absolute amount of $Y$ changes over time.
As we discussed at the beginning of this notebook, how the amount of $Y$  changes over time is a function of two things: 1) a growth term which represents the rate of production of $Y$; and 2) a decay term which gives the rate at which $Y$ is degraded.  A  differential equation describing this as follows:
$$
\frac{dY}{dt} =  f(X_1, X_2, \ldots) - \alpha Y
$$
The $f(X_1, X_2, \ldots)$ term represents the growth term and is a function of the transcription factors that regulate $Y$. We've already seen a couple of ways to model the rate of producting -- using the Hill function or its logic approximation. For the sake of simplicity we'll use the logic approximation to model the growth term.  For example, in the case $Y$ is regulated by a single input we might use $f(X) = \beta \theta(X > K_1)$.  For the equivalent function where $Y$ was regulated by two transcription factor, $X_1$ and $X_2$, and both are required to be above the respective threshold, we could use the function $f(X_1, X_2) = \beta \theta (X_1 > K_1 \land X_2 > K_2)$. 
The second term, $\alpha Y$ represents the rate at which $Y$ is being broken down or diluted.  Notice that the decay rate is a proportional to the amount of $Y$ that is present.
Change in  concentration under constant activation
Now let's explore a simple model of regulation for the two gene network, $X \longrightarrow Y$. Here we assume that at time 0 the activator, $X$, rises above the threshold, $K$, necessary to induce transcription of $Y$ at the rate $\beta$.  $X$ remains above this threshold for the entire simulation. Therefore, we can write $dY/dt$ as:
$$
\frac{dY}{dt} =  \beta - \alpha Y
$$
Write a Python function to represent the change in $Y$ in a given time increment, under this assumption of constant activation:
End of explanation
B = 5
K = 10
a = 0.05
# setup pulse of X
# off (0) for first 50 steps, on for next 100 steps, off again for last 100 steps
X = [0]*50 + [3*K]*100 + [0]*100
Y = [0] # initial value of Y
nsteps = 250
for i in range(1, nsteps):
    xnow = X[i]
    yprevious = Y[-1]
    deltay = dYdt(B, K, a, xnow, yprevious)
    ynew = yprevious + deltay
    Y.append(ynew)
    
plot(X, color='red', linestyle='dashed', label="X")
plot(Y, color='blue', label="Y")
ylim(0, max(max(X)*1.1, (B/a)*1.1))
xlabel('Time units')
ylabel('Concentration')
legend(loc="best")
pass
Explanation: <h2> <font color='firebrick'> Question Set 4 (3 pts) </font> </h2>
The concentration of $Y$ eventually reaches a steady state, $Y_{st}$.  How does $Y_{st}$ relate to $\beta$ and $\alpha$?
The response time of a dynamical system, $T_{1/2}$is defined as the time it takes for it to go half-way between it's initial and final value.  
a) How does the response time change as you vary $\beta$?
b) How does the response time change as you vary $\alpha$?
Estimate the response time as you vary the parameter  $\alpha$ and see if you can create a plot (in Python, or R, or Excel) showing the relationship between $\alpha$ and response time, for $0.01 \leq \alpha \leq 0.1$.
Toggling the activator X
In the proceeding example the activator $X$ was on at the beginning of the simulation and just stayed on.  Let's see what happens when $X$ has pulsatile dynamics.  This would be akin to  toggling $X$ on then off, and asking what happens to $Y$.
End of explanation
## We'll specify the behavior of X as a series of pulse of different length
## so we'll define a function to generate pulses
def pulse(ontime, offtime,  ntimes, onval=1):
    if ontime >= offtime:
        raise Exception("Invalid on/off times.")
    signal = np.zeros(ntimes)
    signal[ontime:offtime] = onval
    return signal
nsteps = 150
short_pulse = pulse(20, 23, nsteps)  # 5 sec pulse
long_pulse = pulse(50, 100, nsteps)  # 50 sec pulse
X = short_pulse + long_pulse  # we can then add the pulses to create
                              # a single time trace
plot(X, color='black')
xlabel('Time units')
ylabel('Amount of Gene Product')
ylim(0, 1.5)
pass
Explanation: Feed Forward Loops
We're now going to use some of these tools to look at a class of network motifs (small network topologies), called Feed Forward Loops (FFLs), found in signaling and regulatory networks. FFLs involve interactions between three components, with the basic topology illustrated below. Depending on the signs of the edges (whether activating or repressing) we can classify FFLs as "coherent" or "incoherent." We'll take a look at an example of each class.
A Coherent FFL
The most common type of coherent FFL is illustrated in the figure below.  In this system $X$ is an activator of $Y$ and both $X$ and $Y$ regulate the production of $Z$ with AND logic (i.e. both $X$ and $Y$ must be above particular thresholds in order to trigger the production of $Z$). 
Using our logic approximation framework we will model the coherent FFL network illustrated above as follows.
Gene Y:
\begin{eqnarray}
Y = f(X) = \beta_y\ \Theta(X > K_{xy})
\
\
\frac{dY}{dt} = \beta_y\ \Theta(X > K_{xy}) - \alpha_{y}Y
\end{eqnarray}    
Gene Z:
\begin{eqnarray}
Z = g(X,Y) = \beta_z\ \Theta(X > K_{xz})\Theta(Y > K_{yz})
\
\
\frac{dZ}{dt} = \beta_z\ \Theta(X > K_{xz})\Theta(Y > K_{yz}) - \alpha_{z}Z
\end{eqnarray}
End of explanation
def dYdt(B, K, a, X, Y):
    if X > K:
        theta = 1
    else:
        theta = 0
    return B * theta - a * Y
def dZdt(B, Kx, Ky, a, X, Y, Z):
    theta = 0
    if (X > Kx) and (Y > Ky):
        theta = 1
    return B * theta - a * Z
## Plot X, Y, and Z on the same time scale
nsteps = 150
short_pulse = pulse(20, 23, nsteps)  # 5 sec pulse
long_pulse = pulse(50, 100, nsteps)  # 50 sec pulse
X = short_pulse + long_pulse
# setup parameters for Y and Z
Y = [0]
betay, alphay = 0.2, 0.1
Kxy = 0.5
Z = [0]
betaz, alphaz = 0.2, 0.1
Kxz = 0.5
Kyz = 1
for i in range(nsteps):
    xnow = X[i]
    ynow, znow = Y[-1], Z[-1]
    
    ynew = ynow + dYdt(betay, Kxy, alphay, xnow, ynow)
    znew = znow + dZdt(betaz, Kxz, Kyz, alphaz, xnow, ynow, znow)
    
    Y.append(ynew)
    Z.append(znew)
    
plot(X, 'k--', label='X', linewidth=1.5)
plot(Y, 'b', label='Y')
plot(Z, 'r', label='Z')
ylim(-0.1, 2.5)
xlabel("Time")
ylabel("Concentration")
legend()
pass
Explanation: Define Python functions for dY/dt and dZ/dt
Recall from above that
\begin{eqnarray}
\frac{dY}{dt} & = & \beta_y\ \Theta(X > K_{xy}) - \alpha_{y}Y \ \
\frac{dZ}{dt} & = & \beta_z\ \Theta(X > K_{xz})\Theta(Y > K_{yz}) - \alpha_{z}Z
\end{eqnarray}
End of explanation
nsteps = 150
p1start = 10
p1duration = 5
p2start = 50
p2duration = 50
short_pulse = pulse(p1start, p1start + p1duration, nsteps)  # short pulse
long_pulse = pulse(p2start, p2start + p2duration, nsteps)  # long pulse
X = short_pulse + long_pulse 
# change this `scale` argument to increase/decrease noise
noise = np.random.normal(loc=0, scale=0.2, size=nsteps) # mean=0, sd=0.2
X = X + noise
# setup parameters for Y and Z
Y = [0]
betay, alphay = 0.2, 0.1
Kxy = 0.5
Z = [0]
betaz, alphaz = 0.2, 0.1
Kxz = 0.5
Kyz = 1
for i in range(nsteps):
    xnow = X[i]
    ynow, znow = Y[-1], Z[-1]
    ynew = ynow + dYdt(betay, Kxy, alphay, xnow, ynow)
    znew = znow + dZdt(betaz, Kxz, Kyz, alphaz, xnow, ynow, znow)
    Y.append(ynew)
    Z.append(znew)
# draw each trace as a subfigure
# subfigures stacked in a vertical grid
subplot2grid((3,1),(0,0))
plot(X, 'k', label='X', linewidth=1)
legend()
subplot2grid((3,1),(1,0))
plot(Y, 'b', label='Y', linewidth=2)
legend()
subplot2grid((3,1),(2,0))
plot(Z, 'r', label='Z', linewidth=2)
vlines(p1start, min(Z),max(Z)*1.1,color='black',linestyle='dashed')
annotate("pulse 1 on", xy=(p1start,1),xytext=(40,20),
         textcoords='offset points',
         horizontalalignment="center",
         verticalalignment="bottom",
         arrowprops=dict(arrowstyle="->",color='black',
                         connectionstyle='arc3,rad=0.5',
                         linewidth=1))
vlines(p2start, min(Z),max(Z)*1.1,color='black',linestyle='dashed')
annotate("pulse 2 on", xy=(p2start,1),xytext=(-40,0),
         textcoords='offset points',
         horizontalalignment="center",
         verticalalignment="bottom",
         arrowprops=dict(arrowstyle="->",color='black',
                         connectionstyle='arc3,rad=0.5',
                         linewidth=1))
legend()
pass
Explanation: <h2> <font color='firebrick'>Question Set 5 (2 pts)</font> </h2>
How do the dynamics of $Y$ and $Z$ differ in the simulation above?
Try varying the length of the first (short) pulse.  How does changing the length of the pulse affect the dynamics of $Y$ and $Z$?
Performance of the Coherent FFL under noisy inputs
Let's further explore the behavior of the coherent FFL defined given noisy inputs.  As before we're going to define an input signal, $X$, that has a short and long pulse, but now we're going to pollute $X$ with random noise.
End of explanation
def Ton(alpha, KYratio):
    return (1.0/alpha) * log(1.0/(1.0-KYratio))
## Create a contour plot for a range of alpha and Kyz/Yst
x = alpha = linspace(0.01, 0.2, 100)
y = KYratio = linspace(0.01, 0.99, 100)
X,Y = meshgrid(x, y)
Z = Ton(X,Y)
levels = MaxNLocator(nbins=20).tick_values(Z.min(), Z.max())
im = contourf(X,Y,Z, cmap=cm.afmhot_r, levels=levels)
contour(X, Y, Z, levels,
                  colors=('k',),
                  linewidths=(0.5,))
colorbar(im)
xlabel('alpha')
ylabel("Kyz/Yst")
pass
Explanation: <h2> <font color='firebrick'>Question Set 6 (2 pts) </font> </h2>
In the code cell above, try changing the duration of the first pulse and the scale of the noise (see comments in code) to get a sense of how good a filter the FFL is.  Is there a bias to the filtering with respect to Z turning on versus Z turning off?
Dynamics of Y and Z in the Coherent FFL
As before we can solve for Y as a function of time and calculate what its steady state value will be:
$$
Y(t) = Y_{st}(1-e^{-\alpha_{y}t})
$$
and 
$$
Y_{st}=\frac{\beta_y}{\alpha_y}
$$
How about $Z$?
Since $Z$ is governed by an AND function it needs both $X$ and $Y$ to be above their respective thresholds, $K_{xz}$ and $K_{yz}$. For the sake of simplicity let's assume that both $Y$ and $Z$ have the same threshold with respect to $X$, i.e. $K_{xy} = K_{xz}$. This allows us just to consider how long it takes for $Y$ to reach the threshold value $K_{yz}$. Given this we can calculate the delay before $Z$ turns on, $T_{\mathrm{on}}$ as follows.
$$
Y(T_{\mathrm{on}}) = Y_{st}(1-e^{-\alpha_y T_{\mathrm{on}}}) = K_{yz}
$$
and solving for $T_{\mathrm{on}}$ we find:
$$
T_{\mathrm{on}} = \frac{1}{\alpha_y} \log\left[\frac{1}{(1-K_{yz}/Y_{st})}\right]
$$
Thus we see that the delay before $Z$ turns on is a function of the degradation rate of $Y$ and the ratio between $Y_{st}$ and $K_{yz}$.  
Exploring the Parameter space of $Z$'s turn-on time
From the above formula, we see that there are two parameters that affect the turn-on time of $Z$ -- $\alpha_y$ (the scaling factor for the decay rate of $Y$)  and the compound parameter $K_{yz}/Y_{st}$ (the threshold concentration where $Y$ activate $Z$ relative to the steady state of $Y$).  To explore the two-dimensional parameter space of $Z's$ $T_on$ we can create a contour plot.
End of explanation
## A Python function that represents dZ/dt for the Incoherent FFL
## our dY function previously defined stays the same
def dZ_incoh(B1,B2,Kx,Ky,a,X,Y,Z):
    pass # define the function here
def dZ_incoh(B1,B2,Kx,Ky,a,X,Y,Z):
    theta = 0
    B = 0
    if (X > Kx) and (Y < Ky):
        theta = 1
        B = B1
    elif (X > Kx) and (Y >= Ky):
        theta = 1
        B = B2
    return B * theta - a * Z
nsteps = 150
short_pulse = pulse(20, 25, nsteps)  # 5 sec pulse
long_pulse = pulse(50, 100, nsteps)  # 50 sec pulse
X = short_pulse + long_pulse
# setup parameters for Y and Z
Y = [0]
betay, alphay = 0.2, 0.1
Kxy = 0.5
Z = [0]
betaz1, betaz2 = 0.2, 0.001
alphaz = 0.1
Kxz = 0.5
Kyz = 0.5
for i in range(nsteps):
    xnow = X[i]
    ynow, znow = Y[-1], Z[-1]
    
    ynew = ynow + dYdt(betay, Kxy, alphay, xnow, ynow)
    znew = znow + dZ_incoh(betaz1, betaz2, Kxz, Kyz, alphaz, xnow, ynow, znow)
    
    Y.append(ynew)
    Z.append(znew)
    
# draw each trace as a subfigure
# subfigures stacked in a vertical grid
subplot2grid((3,1),(0,0))
plot(X, 'k', label='X', linewidth=1)
legend()
ylim(0,1.1)
subplot2grid((3,1),(1,0))
plot(Y, 'b', label='Y', linewidth=2)
legend()
ylim(0,2.1)
subplot2grid((3,1),(2,0))
plot(Z, 'r', label='Z', linewidth=2)
legend()
ylim(0,0.7)
pass
Explanation: Type 1 Coherent FFLs can act as a Sign-Sensitive Delays
As discussed in the article by Shen-Orr et al. a feed forward loop of the type we've just discussed can act as a type of filter -- a sign-sensitive delay that keeps $Z$ from firing in response to transient noisy signals from $X$, but shuts down $Z$ immediately once the signal from $X$ is removed. 
An Incoherent FFL
Consider the FFL illustrated in the figure below.  
In this incoherent FFL, the logic function that regulates $Z$ is "X AND NOT Y".  That is $Z$ turns on once $X$ is above a given threshold, but only stays on fully as long as $Y$ is below another threshold. Again for simplicity we assume $K_{xy} = K_{yz}$. 
Dynamics of Y
As before, the dynamics of $Y$ are described by:
$$
\frac{dY}{dt} = \beta_y\ \Theta(X > K_{xy}) - \alpha_{y}Y
$$
and   
$$
Y(t) = Y_{st}(1-e^{-\alpha_{y}t})
$$
Dynamics of Z
To describe $Z$ we consider two phases - 1) while $Y < K_{yz}$ and 2) while $Y > K_{yz}$. 
Z, Phase 1
For the first phase:
$$
\frac{dZ}{dt} = \beta_z\ \Theta(X > K_{xz}) - \alpha_{z}Z
$$
and
$$
Z(t) = Z_{m}(1-e^{-\alpha_{z}t})
$$
As we did in the case of the coherent FFL, we can calculate the time until $Y$ reaches the treshold $K_{yz}$. We'll call this $T_{\mathrm{rep}}$ and it is the same formula we found for $T_{\mathrm{on}}$ previously.
$$
T_{\mathrm{rep}} = \frac{1}{\alpha_y \log[\frac{1}{1-K_{yz}/Y_{st}}]}
$$
Z, Phase 2
After a delay, $T_{\mathrm{rep}}$, $Y$ starts to repress the transcription of $Z$ and $Z$ decays to a new lower steady state, $Z_{st} = \beta_{z}^{'}/\alpha$. The value of  $\beta_{z}^{'}$ depends on how leaky the repression of $Z$ is by $Y$.  
The dynamics of $Z$ in Phase 2 is given by:
$$
Z(t) = Z_{st} + (Z_0 - Z_{st})e^{-\alpha_{z}(t-T_{\mathrm{rep}})}
$$
where
$$
Z_0 = Z_{m}(1-e^{-\alpha_{z}T_{\mathrm{rep}}})
$$
Combining the two phases of Z
We can combine the two phases of $Z$ into a single function:
$$
f(X,Y) = \beta_z\Theta(X > K_{xz} \land Y < K_{yz}) + \beta_{z}^{'}\Theta(Y \geq K_{yz}) - \alpha_z Z
$$
End of explanation
<END_TASK>