BENG0091 Coursework 1
To be submitted on Moodle by 5
th Mar 2021
“Chemical clocks” are reactions in which the concentration of a reactant or product species, exhibits
oscillations over time. When such reactions are performed under well-mixed conditions in large
volumes, the resulting oscillations appear quite regular; however, if they occur in small volumes (e.g.
the interior of a biological cell), they are subject to noise. This happens because reaction occurrence
is a random process: the precise timing of a reaction event cannot be known with certainty, and
therefore the molecular populations exhibit stochastic fluctuations. Due to the “law of large
numbers” these fluctuations become smaller (and eventually negligible) as the system size increases.
In order to study the fluctuations that are prevalent on small systems, an algorithm proposed by
Gillespie [1] can be used. The algorithm keeps track of the populations of the molecular species (i.e.
the numbers of molecules) that exist in the system at all times. These populations change as
reactions occur, for instance in a dimerisation reaction like 2A B, two molecules of species A are
consumed, and one molecule of species B is produced. The timing of the reaction events is
calculated by the “waiting times” for reaction events to occur, and these waiting times are modelled
by exponentially distributed random variables. Moreover, if there is more than one possible
reaction, the algorithm uses a discrete random variable to select the reaction to occur; each reaction
has probability of occurrence that is proportional its “propensity”. The latter can be calculated easily
from the numbers of molecules and the kinetic constant of the reaction (which will be known). To
learn more about Gillespie’s algorithm you can consult Ref. [1]; for the purposes of this coursework,
the pseudocode given at the Appendix (in the end of this document) is sufficient.
Thus, for this coursework you are asked to do the following:
1. Code (in Matlab or similar) Gillespie’s algorithm that simulates the following two decay reactions
of molecular species X and Y:
Of course mass has to be conserved, so in our reactions, symbol
denotes molecular species
that we are not interested in. Thus, your algorithm will keep track of the populations of species
X and Y over time. These populations (numbers of molecules) can take the values X = 0, 1, 2, …
(similarly for Y). Use the following values for the kinetic constants cd1 and cd2, and the initial
populations, X0 and Y0:
cdX = 0.01
cdY = 1.5
X0 = 300
Y0 = 350
The propensities are given as:
(a) In the same graph, plot 10 realisations of the time evolution of the system as propagated
by your stochastic algorithm, for times t = 0 to 4 (for each realisation you will have to use
a different seed for your random number generator). [10]
Always in the same plot, add the corresponding deterministic solution for the time
evolution of the system, as given by the following equations:
The above equations were obtained by solving ordinary differential equations (ODEs)
corresponding to your stochastic process. Comment on how well these ODE solutions
capture the behaviour of the system. [5]
(b) Now we will increase the size of the system by a factor of 10. Thus, repeat the procedures
of question (1.a) for the following parameters:
cdX = 0.001
cdY = 1.5
X0 = 3000
Y0 = 3500
Make a graph that shows 10 realisations of the stochastic process along with the
deterministic solutions (again for time t = 0 to 4). Compare and contrast these plots with
those of question (1.a) and comment on how well the deterministic solutions capture the
behaviour of this “larger” system. [15]
2. Now, let us model a more complicated system, which was devised by I. Prigogine and co-workers
at Université Libre de Bruxelles [2] (and was named “Brusselator” after the city of Brussels,
Belgium), as a prototype chemical oscillator. The following three reactions are possible in this
system (note that we have simplified the notation, compared to the original work of Prigogine et al.):
denotes other species that we are not interested in monitoring. Your algorithm
will therefore only focus on the population of two species: X and Y. We will denote the size of
the system with S, and thus, the values of the kinetic parameters are given as:
The propensities are given as:
(a) In 3 separate graphs, plot 3 realisations of the time evolution of the system for the
following values of size: S = 0.1, 1 and 10. The trajectories of both species should appear
in each plot. For all simulations: start from an initial population of X0 = Y0 = 0, and go up to
a maximum time tmax = 10, set the maximum number of events nmax = 108 and the
sampling time tsample = 0.005. Comment on the behaviour of the system, as captured by
these simulations. [20]BENG0091编程代写、代做R编程
(b) For the system with size parameter S = 0.1, we would like to estimate the probability
distribution of the population of molecular species X. To this end, simulate a longer
realisation setting tmax = 100, and take samples of the process over constant time
intervals of tsample = 0.001 time units. Plot the probability distribution of the population
of X, and comment on the shape of this distribution. [20]
(c) Out of the samples taken over these constant time intervals, in question (2.b), calculate
the following quantities for the population of species X (write your own code to do this,
do not use the intrinsic functions available on Matlab, Python etc.):
i. the estimated mean,
ˆ
ii. the estimated standard deviation,
ˆ
iii. the skewness from the following estimator:
(d) Finally, present your Matlab or Python codes as Appendixes to your report. Make sure
these are copy-pasted as text, not as figures. [10]
References
[1] Gillespie, D. T. (1977). “Exact Stochastic Simulation of Coupled Chemical Reactions”. The Journal
of Physical Chemistry. 81 (25): 2340-2361. doi: 10.1021/j100540a008
[2] I. Prigogine, From Being to Becoming: Time and Complexity in the Physical Sciences, New York:
W.H. Freeman and Company, 1980
4 of 4
Appendix: Pseudocode for Gillespie’s Algorithm
1. Start
2. Perform initialisation procedures
2.1. Input values for kinetic constants cj, j = 1…m
2.2. Input initial populations of all species Xi, i = 1…n
2.3. Specify the maximum time tmax that the simulation will run for, as well as the maximum
allowed number of steps nmax
2.4. Specify the sampling interval tsample
2.5. Initialise the simulated time: t = 0, as well as the reaction counter: cnt = 0
2.6. Initialise the sampling time: tsample = 0, as well as the number of samples: Nsamples = 0
2.7. Initialise uniform pseudo-random number generator
3. Loop while t < tmax and cnt < nmax
3.1. Calculate the reactions’ propensity functions:
αj j j X X c h for j 1,2,...,m
3.2. Calculate the total propensity:
3.3. Generate two uniformly distributed pseudo-random numbers r1 and r2 (0,1)
3.4. Calculate the time for the next reaction:
3.5. Find which reaction will be next by solving:
μ μ
3.6. Realize the reaction and perform sampling:
3.6.1. Update time t = t +
3.6.2. Loop while t > tsample
3.6.2.1. Update number of samples: Nsamples = Nsamples + 1
3.6.2.2. Save the time and the populations of species in a file (or in an array)
3.6.2.3. Update sampling time: tsample = tsample + tsample
3.6.3. Update species populations Xi to reflect the occurrence of reaction R
3.6.4. Update reaction counter: cnt = cnt + 1
4. Stop
Note: on steps 3.6, notice that sampling is done after advancing the time and before updating the
molecular populations. This is so that the correct state is saved to the file (or array), since up until
time t + the system is “in the old” state. The new state “takes effect” immediately after time t + .
如有需要，请加QQ：99515681 或WX：codehelp