1 of 4 BENG0091 Coursework 1 To be submitted on Moodle by 5th 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: X dX c2 Y dY c 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: α dX dX X X 1 X c 2 α dY dYY c Y 2 of 4 (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: 1 dX 0 1 X c t X 0 dYY Y exp c t [5] 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.): X 1c X Y2c 2X+Y 3X3 c X 4c Again, symbol 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: c1 = 4800 S c2 = 20 c3 = 3 10– 5 / S2 c4 = 6 3 of 4 The propensities are given as: α 1 1X,Y c α 2 2X,Y c X α 3 3 X X 1 X,Y c Y 2 α 4 4X,Y c X (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: 3 1 ˆ ˆ ˆ1 2 N i i XN N N Comment on your results. For instance, is the mean representative of the population and why (or why not)? Do the values of ˆ and ˆ make sense? [15] (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 jc h for j 1,2,...,mX X 3.2. Calculate the total propensity: α α m 0 j j 1 X X 3.3. Generate two uniformly distributed pseudo-random numbers r1 and r2 (0,1) 3.4. Calculate the time for the next reaction: τ α0 1 1 1 ln r 3.5. Find which reaction will be next by solving: μ μ α α α 1 j 2 0 j j 1 j 1 r 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 + .
欢迎咨询51作业君