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.

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.

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)
```

- 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)
```

- 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:

- 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)
```

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

In []:

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
```

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
```

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.

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 []:

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.

In []:

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

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

In []:

- How are the results different from the SIM case?

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