BENG0091编程代写、代做R编程代写、R程序设计调试

日期：2021-02-26 10:41

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]

(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 + ?.