This HTML version of the exercises is only intended to give you a quick view of what the exercises are about. Use the .ipynb file to solve the exercises!

Exercise 3: Simulation of gene expression and regulation (DL 20 November)

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. Each task is worth 2 points.

Task 1: Discrete event simulation I: Simple systems¶

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.

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)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps)
    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:

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)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps)
    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:

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)
    Nhist[0] = N
    # Times of the corresponding population sizes
    Thist = np.zeros(numsteps)
    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¶

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, len(state)))
    Nhist[0,:] = state
    # Time history
    Thist = np.zeros(numsteps)
    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:

Try to find rates for the different equations that produce interesting behaviour.
Plot examples of the behaviour under different rates.

Task 3: Continuous-state simulation¶

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 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 repression¶

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:

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.
Repeat the stochastic simulations 100 times and compute the average. Visualise a number of the individual simulations.

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.

Task 5: Simulating network motifs¶

a) Use the simulators to simulate from the SIM (single input motif) network motif.

Using a model with regulator \(TF\) and targets \(X_i\) the system can be expressed through reactions:

Starting with relatively high TF abundance and fairly low decay will give interesting profiles.

You can experiment with rate constants to obtain interesting results. One set of rate constants that yield reasonable behaviour is:

target_transcription_rates = np.logspace(-1, 1, N) target_decay_rates = np.logspace(-1, 2, N) TF_decay_rate = 10 TF_initial_abundance = 100

Simulate the system and visualise the results.

Briefly describe what kind of behaviours you observe and how the rates affect the outcome.

b) Use the simulators to simulate from the feed-forward loop network motif.

Similar as in a) above, construct a set of target genes following the feed-forward loop motif, simulate them and visualise the results.

How are the results different from the SIM case?

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\).