In this assignment, we will use ideas from non-linear dynamics to investigate genetic switches and oscillators.

The assignment should be handed in as a report with copies of the code you used. Make sure you include a copy of all your code, including that for generating the figures, and that the code is commented with the appropriate question number.

A Jupyter notebook is the recommended way to carry out the assignment. Use

  import numpy as np
  import matplotlib.pyplot as plt
  %matplotlib inline

at the start of your notebook to import the modules needed and to have figures appearing in the notebook rather than in a separate window.

To print from Jupyter, you can create a Print Preview and then use your web browser to print this preview. Don't worry if extra space is added to the printed version.

You can work in pairs, but not in groups larger than two. Within a pair, you can have the same computer code, but you must have different write-ups. Please mark on your script the examination number of the person with whom you are collaborating. If there are more than two submitted assignments that have sections of the same computer code then we will have to give zero marks to those assignments.

Marks for each section are given in bold.

You only need to use odeint, and not fsolve, in the assignment because we are studying only the stable steady states.

A genetic toggle switch

Gardner and colleagues constructed a genetic switch using two repressors, each of which repressed the other. To be able to switch the genetic network from one state to the other, they used inducers – small molecules that induce gene expression. The inducers inhibit the action of the repressors because a repressor bound by an inducer is unable to bind to DNA. If the two repressors are $u$ and $v$ and $I_u$ is the inducer that binds to $u$ and $I_v$ is the inducer that binds to $v$, then we can model the system as

and

where $a_u$ is the maximum rate of expression of $u$ and $a_v$ is the maximum rate of expression of $v$. The binding affinity of an inducer molecule to a repressor is given by $K$, and $b$ and $n$ are Hill numbers: $b$ describes the cooperativity of repression and $n$ describes the cooperativity of inducer binding to repressors. The degradation rates of $u$ and $v$ are assumed equal and set the unit of time (and so are both one).

Question 1.1: A true switch should exhibit history-dependent behaviour. Gardner and colleagues demonstrated such hysteresis by applying a series of inducers to their genetic construct. We will repeat this experiment with the model. Write some code to simulate the dynamics of $u$ and $v$ using odeint. You will need to create a function that returns $du\/dt$ and $dv\/dt$ given $u$ and $v$ and $t$. As we will vary the concentration of inducers, the function should also take the concentrations of the inducers as inputs. Let $a_u= 10$, $a_v= 9$, $K= 3$, and both Hill numbers equal 2. We will use this code in the next question. [2]

Question 1.2: Starting with $u=v=0$, simulate the application of $I_u$ for 50 time units (set $I_u=100$ and $I_v=0$), then remove $I_u$ for another 50 time units. Following this removal, apply $I_v$ for 50 time units (set $I_u=0$ and $I_v=100$) and then remove $I_v$ for another 50 time units. You will need to call odeint four times and use the final values of $u$ and $v$ from the last call to odeint as the inital values for the current call to odeint. Plot your results (remembering to use show() to view the figure). [4] Do you see hysteresis (history-dependent behaviour)? Why? Change both Hill numbers to 1. Rerun the experiment and explain what you observe. [3]

Question 1.3: Determine the nullclines of the model system. To investigate the stability of the fixed points, plot the nullclines on the same graph with $u$ on the $x$-axis and $v$ on the $y$-axis and for $I_u= I_v= 0$. Picking 10 initial conditions for $u$ and $v$, use odeint to determine their time evolution and plot the trajectories on the same graph as your nullclines. You can use the alpha option in plt.plot to increase the transparency of the plotted trajectory. Mark the start and the end points of each trajectory with a point. Explain which fixed points are stable and, from the trajectories you plotted, why. [5]

Question 1.4: Determining the bifurcation diagram for a bistable system is another way to demonstrate that the system exhibits hysteresis. Assume that the gene coding $u$ is positively controlled by another transcription factor. Increasing the concentration of this transcription factor will increase the parameter $a_u$ in the equation for $du\/dt$. We will assume that the concentration of transcription factor changes only slowly and that the system reaches steady-state before the concentration of transcription factor changes again. Assuming $I_u= I_v= 0$, write code to increase $a_u$ from 0 to 35, and find the steady-state value of $u$ for each value of $a_u$ using the steady-state values of $u$ and $v$ for the previous value of $a_u$ as the initial conditions for the current value of $a_u$. Now decrease $a_u$ from 35 to 0 using, as before, the steady-state values for the previous $a_u$ as the initial values for the new value of $a_u$. You can use the np.linspace command to set the values of $a_u$. Plot on the same graph the steady-state value of $u$ against $a_u$ when $a_u$ is increased and when $a_u$ is decreased. Explain what you see. Is there a range of values for $a_u$ for which the switch once thrown on cannot be reset? [6]

A genetic oscillator

Elowitz and Leibler were the first people to construct a synthetic genetic oscillator. They added three genes to E. coli: the first gene repressed the second; the second gene repressed the third; and the third gene repressed the first.

Question 2.1: Explain why this system has negative feedback. [2]

In the model used by Elowitz and Leibler to design their system, transcription for each gene is described as

where $m_i$ is the mRNA of the $i$'th gene and $p_i$ is the repressor encoded by that mRNA. For example, transcription of the second gene is repressed by the protein transcribed from the first gene with this repression described by a Hill function with a Hill number of 2. Time is given in units of the mRNA lifetime, and $a$ is the maximum number of mRNAs transcribed per unit time. Both the mRNA lifetime and $a$ are assumed to be the same for all three genes, and note that $p_0$ is understood to mean $p_3$. Translation and protein degradation are modelled as

where $b$ is the ratio of mRNA to protein lifetimes.

Question 2.2: Fix $b=10$ and by varying $a$, produce six plots of the concentration of one of the proteins as the system goes through a bifurcation from a steady-state to a limit cycle (hint: start your investigation around the value of 14). Use odeint to perform the simulations and plt.subplot to draw all six plots in the same figure window. You will need to create a function that returns $dm_i\/dt$ and $dp_i\/dt$ for all three genes ($i=1$, 2 and 3). Make sure that at least one of the repressors has a positive initial value. If there are no repressors and no mRNAs, the system stays at this unstable steady-state. Label each subplot with its corresponding value of $a$. [4]

Question 2.3: An inefficient but effective way to determine the bifurcation diagram for this system is to to form a grid of $a$ and $b$ values, to simulate the system for each point on the grid, and then classify its behaviour as either oscillatory or steady-state. Write some Python code to automate this procedure. Let $a$ run between 1 and 10³ and $b$ run between 0.1 and 10³ (use np.logspace). Loop through each point on the grid (each particular value of $a$ and $b$) and use odeint to simulate the dynamics with these parameter values. Select a portion of the time-series of one of the proteins and use np.diff to determine if the protein is either at steady-state or is oscillating. To plot your bifurcation diagram, use the plt.contourf function. To set the axes to a log-scale, you can use, for example,

plt.xscale('log')

and 'linear' to change the x-axis back. [4] Explain why your bifurcation diagram has the shape it does. [4]