Exercise 5: Simulation of dynamic models of gene expression and regulation (DL 13 October)

Lecturer: Antti Honkela

Your name:

Name of your pair:

In this exercise you will implement simulators to simulate biochemical reaction systems similar to ones that could be used to simulate gene expression.

Task 1: Discrete event simulation I: Simple systems (2 points)

In this task your aim is to perform stochastic simulation of simple discrete reaction systems with a single species X in different scenarios. Each scenario comes with a number of questions. Please answer these by having your code print out the answers.

a) Simulate the decay of a population of 100 individuals with decay rate 1.

The corresponding reaction system is: $$ X \rightarrow \emptyset $$

  • Implement a function that runs the simulation.
  • Repeat the simulation 100 times. What is the average time until extinction? What are the minimum and maximum times?
  • What is the mean and standard deviation of the number of individuals left alive after 1 unit of time?
  • Plot the simulated numbers of individuals as function of time together in a single plot.

Hints:

  • numpy.random (imported as npr) has many useful functions for simulating random variables with specific distributions, e.g. npr.exponential and npr.choice (for subtasks b and c)
In [ ]:
%matplotlib inline
# Import some useful modules
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as npr

# Suggested structure: simulate the decay of a population of N0 individuals at rate 'rate' for 'numsteps' steps
def simulate_decay(N0, rate, numsteps):
    # Current population size
    N = N0
    # Current time
    t_cur = 0
    # History of population size across the simulation
    Nhist = np.zeros(numsteps+1)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps+1)
    Thist[0] = t_cur
    # ...
    # ADD CODE FOR PERFORMING THE SIMULATION HERE
    # ...
    return (Thist, Nhist)

b) Simulate a birth-death process of a population of 100 individuals with decay rate 1 and birth rates [0.8, 1, 1.2] for 10000 steps.

The corresponding reaction system is: \begin{align} X &\rightarrow 2X \ X &\rightarrow \emptyset \ \end{align}

  • Implement a function that runs the simulation.
  • Repeat each simulation 100 times. How many simulations lead to extinction in each scenario? What are the statistics of the final populations in each scenario (mean, min, max)?
In [ ]:
def simulate_birthdeath(N0, birthrate, deathrate, numsteps):
    # Current population size
    N = N0
    # Current time
    t_cur = 0
    # History of population size across the simulation
    Nhist = np.zeros(numsteps+1)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps+1)
    Thist[0] = t_cur
    # ...
    # ADD CODE FOR PERFORMING THE SIMULATION HERE
    # ...
    return (Thist, Nhist)

c) Simulate a migration-death process of a population of 100 individuals with decay rate 1 and migration rates [1.0, 10.0, 100.0] for 10000 steps.

The corresponding reaction system is: \begin{align} \emptyset &\rightarrow X \ X &\rightarrow \emptyset \ \end{align}

  • Implement a function that runs the simulation.
  • Repeat each simulation 100 times. How many simulations lead to extinction in each scenario? What are the statistics of the final populations in each scenario (mean, min, max)?
In [ ]:
def simulate_migrationdeath(N0, migrationrate, deathrate, numsteps):
    # Current population size
    N = N0
    # Current time
    t_cur = 0
    # History of population size across the simulation
    Nhist = np.zeros(numsteps+1)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps+1)
    Thist[0] = t_cur
    # ...
    # ADD CODE FOR PERFORMING THE SIMULATION HERE
    # ...
    return (Thist, Nhist)

Briefly summarise how each model behaves qualitatively in the text area below:

Task 2: Discrete event simulation II: The general case (3 points)

a) Implement the general Stochastic Simulation Algorithm of Gillespie (1977).

  • Implement a function that runs the simulator.
  • Test your simulator by re-implementing exercises 1a and 1b to make sure your results are similar.
In [ ]:
# Suggested skeleton for the simulation algorithm
# Suggested arguments are NumPy arrays with sizes as follows:
# inputs: [reactions x state]
# outputs: [reactions x state]
# rates: [reactions]
def stochastic_simulation_algorithm(initstate, inputs, outputs, rates, numsteps):
    # Current time and state
    t_cur = 0
    state = 1.0 * initstate
    # State history
    Nhist = np.zeros((numsteps+1, len(state)))
    Nhist[0,:] = state
    # Time history
    Thist = np.zeros(numsteps+1)
    Thist[0] = t_cur
    # ...
    # ADD CODE FOR PERFORMING THE SIMULATION HERE
    # ...
    return (Thist, Nhist)

b) Use your simulator to simulate the Lotka-Volterra predator-prey system:

\begin{align} X_1 + X_2 &\rightarrow 2 X_1 \quad &\text{# Predator-prey interaction} \\ X_2 &\rightarrow 2 X_2 &\text{# Prey reproduction} \\ X_1 &\rightarrow \emptyset &\text{# Predator death} \\ \end{align}
  • Try to find rates for the different equations that produce interesting behaviour.
  • Plot examples of the behaviour under different rates.
In [ ]:

Task 3: Continuous-state simulation (3 points)

a) Implement a deterministic ordinary differential equation simulator using the same input as the stochastic simulator above.

Your simulator should follow Eq. (11) in Gillespie (2000) paper. In order to discretise the time, please select a finite time step \tau and apply the Euler discretisation of approximating $$ \frac{dX(t)}{dt} = f(X(t)) $$ by $$ X(t + \tau) \approx X(t) + \tau \cdot f(X(t)). $$

  • Repeat exercises 1a, 1b and 2b using your simulator. How do the results differ from the original?
In [ ]:
# Suggested skeleton with same arguments as above, except tau is the time step
def simulate_reaction_ode(initstate, inputs, outputs, rates, numsteps, tau):
    pass

b) Implement a stochastic continuous-state simulator using the same input as above.

Your simulator should follow Eq. (21) of the Gillespie (2000) paper.

  • Repeat exercises 1a, 1b and 2b using your simulator. How do the results differ from the original and from 3a?
In [ ]:
# Suggested skeleton with same arguments as above, except tau is the time step
def simulate_chemical_langevin_eq(initstate, inputs, outputs, rates, numsteps, tau):
    pass

Task 4: Simulating a system of repression of expression (2 points)

In this task you will implement a simplified model of repressive regulation, a simplified "repressilator". This is more complicated than the predator-prey system in 2b because preventing preduction is not the same as enhancing degradation. (In the predator analogue, predators would only target the young and therefore decrease the birth rate of prey but leave the existing (adult) population intact.

The repressilator model has 3 genes, each of which represses the next in a cycle A -> B -> C -> A. A realistic model would include 6 species to model both mRNAs and protein levels of each gene. In this task we simplify things a bit by combining the mRNA and protein for each gene.

In order to simulate repression, we augment the model by introducing additional states $aX$ and $rX$ that represent "active" and "repressed" promoter states of gene $X \in \{A, B, C\}$. Production of gene $X$ only happens when the promoter is in the active state and the repressor moves the promoter into repressed state. The system can be expressed through reactions: \begin{align} aX &\rightarrow aX + 10 X \quad &\text{# transcription} \ Y + aX &\rightarrow rX &\text{# repression, } Y \text{ is the repressor}\ rX &\rightarrow aX &\text{# recovery from repression} \ X &\rightarrow \emptyset &\text{# decay} \end{align} These reactions are repeated for all genes, resulting in total 12 reactions for the 9 species in the system.

  • Implement the reactions and simulate them using the simulators you implemented in Tasks 2 and/or 3. Use rate 0.1 for the decay equation and rate 1.0 for the others.
  • Remember that aX and rX represent different states of the single promoter of a gene, there are not supposed to be many copies of them.
  • Repeat the stochastic simulations 100 times and compute the average. Visualise the gene counts from a number of the individual simulations.
In [ ]:

Bonus task: implement a full repressilator model including both mRNAs and proteins. Try to assign realistic rates for all the reactions. Simulate your model using the simulators you implemented in Tasks 2 and/or 3 and visualise the results.

In [ ]:

General hints and warnings

  • Remember that Python passes function arguments by reference. This means that if you change the argument in a function, it will also change in the caller. If you have a simulator that changes its initial state, you may inadvertedly change its value in the caller too, which may affect future calls if you are running the simulator in a loop. Directly assigning the value to a new variable is not enough because assignments are also by reference.
  • The continuous-state simulators can be sensitive to the chosen $\tau$ value. If the simulation is unstable, try decreasing $\tau$.