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.
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$.
In [2]:
# import statements to make numeric and plotting functions available
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
In [3]:
## define your function in this cell
def hill_activating(X, B, K, n):
Xn = X**n
return (B * Xn)/(K**n + Xn)
In [5]:
## generate a plot using your hill_activating function defined above
# setup paramters for our simulation
B = 5
K = 10
x = linspace(0,30,200) # generate 200 evenly spaced points between 0 and 30
y = hill_activating(x, B, K, 1) # hill fxn with n = 1
plot(x, y, label='n=1')
xlabel('Concentration of X')
ylabel('Promoter activity')
legend(loc='best')
ylim(0, 6)
pass
Following the example above, generate plots for the Hill function where $n = {1, 2, 4, 8}$. Note that you can generate multiple curves in the same figure by repeatedly calling the plot function.
In [12]:
## generate curves for different n here
# setup paramters for our simulation
B = 5
K = 10
x = linspace(0,30,200) # generate 200 evenly spaced points between 0 and 30
y1 = hill_activating(x, B, K, 1) # hill fxn with n = 1
y2 = hill_activating(x, B, K, 2) # hill fxn with n = 2
y4 = hill_activating(x, B, K, 4) # hill fxn with n = 4
y8 = hill_activating(x, B, K, 8) # hill fxn with n = 8
plot(x, y1, label='n=1')
plot(x, y2, label='n=2')
plot(x, y4, label='n=4')
plot(x, y8, label='n=8')
xlabel('Concentration of X')
ylabel('Rate of production of Y')
legend(loc='best')
ylim(0, 6)
pass
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.
In [8]:
## define your repressive hill function in this cell
def hill_repressing(X, B, K, n):
return B/(1.0 + (X/K)**n)
In [11]:
## generate a plot using your hill_activating function defined above
## For X values range from 0 to 30
B = 5
K = 10
x = linspace(0,30,200)
plot(x, hill_repressing(x, B, K, 1), label='n=1')
plot(x, hill_repressing(x, B, K, 2), label='n=2')
plot(x, hill_repressing(x, B, K, 4), label='n=4')
plot(x, hill_repressing(x, B, K, 8), label='n=8')
xlabel('Conc. of X')
ylabel('Rate of production of Y')
legend(loc='best')
ylim(0, 6)
pass
Download the file hill-fxn.py and pnsim.py from the course website. Run this application from your terminal with the command: python hill-fxn.py.
Run the hill-fxn.py script from your terminal by typing python hill-fxn.py.
There are three sliders at the bottom of the application window. You can drag the blue regions of these sliders left or right to change the indicated parameter values. The exact values of each parameter are shown to the right of the sliders. As you drag the sliders the plot will update to show you what the Hill function looks like for the combination of parameters you have currently specified.
Also note there is a dashed vertical line in the plot window. When you move your mouse over the plot window this line will follow your position. As you do so, x- and y-plot values in the lower left of the application window will update to show you the exact position your mouse is pointing to in the plot. The dashed line and the plot readout are useful for reading values off the plot.
Vary the parameter $n$ over the range 1 to 10.
Vary the parameter $\beta$. How does changing $\beta$ change:
Vary the parameter $K$. How does changing $K$ change:
Download and run the script hill-fxn-wlogic.py -- This is like the previous hill-fxn.py script except it now include a set of buttons for toggling the logic approximation on and off. As before vary the parameters $n$, $\beta$ and $K$.
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$ is zero if the statement inside the parentheses is false or one 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) $$
In [16]:
## write your logic approximation functions here
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
Generate functions comparing the logic approximation to the Hill function, for the activating case:
In [17]:
## generate plots using your hill_activating and logic_activating functions defined above
## For X values range from 0 to 30
B = 5
K = 10
n = 4
x = linspace(0, 30, 200)
plot(x, hill_activating(x, B, K, n), label='n=8')
logicx = [logic_activating(i, B, K) for i in x]
plot(x, logicx, label='logic approximation')
xlabel('Concentration of X')
ylabel('Promoter activity')
ylim(-0.1, 5.5)
legend(loc='best')
pass
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) $$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.
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:
In [28]:
## write a function to represent the simple differential equation above
def dYdt(B,a,Y):
return B - a*Y # write your code here
In [24]:
## generate a plot using your dY function defined above
## Evaluated over 200 time units
Y = [0] # initial value of Y
B = 0.2
a = 0.05
nsteps = 200
for i in range(nsteps):
deltay = dYdt(B, a, Y[-1])
ynew = Y[-1] + deltay
Y.append(ynew)
plot(Y)
ylim(0, 4.5)
xlabel('Time units')
ylabel('Concentration of Y')
pass
In [1]:
# setup pulse of X
# off (0) for first 50 steps, on for next 100 steps, off again for last 100 steps
X = [0]*50 + [1]*100 + [0]*100
Y = [0]
B = 0.2
K = 0.5
a = 0.05
nsteps = 250
for i in range(1, nsteps):
xnow = X[i]
growth = logic_activating(xnow, B, K)
decay = a*Y[-1]
deltay = growth - decay
ynew = Y[-1] + deltay
Y.append(ynew)
plot(X, color='red', linestyle='dashed', label="X")
plot(Y, color='blue', label="Y")
ylim(0, 4.5)
xlabel('Time units')
ylabel('Concentration')
legend(loc="best")
pass
The simple-activation.py script simulates the case of constant transcriptional activation as described above. Run the simple-activation.py script from the terminal, using the command python simple-activation.py. This script includes three sliders:
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.
Estimate the response time from the simple-activation.py simulation and see if you can create a plot showing the relationship between $\alpha$ and response time, for $0.1 \leq \alpha \leq 2.0$.