On the Web
A Convenient Tool for the Stochastic Simulation of Reaction Mechanisms Nikolaos Bentenitis Department of Chemistry and Biochemistry, Southwestern University, Georgetown, TX 78626;
[email protected] Given the chemical reaction
A
k
B
(1)
most physical chemistry students will attempt to find the amount, or concentration of A as a function of time and at a certain temperature, by integrating the following rate equation
d k dt
(2)
where [A] is the molar concentration of A. A careful student however will note that eq 2 is valid only when several requirements are met: First, the reaction must be elementary, otherwise there will be an exponent on [A], on the right-hand side of eq 2, different from 1. Second, the system must be closed, otherwise some A might escape. And third, the system volume must be constant, otherwise the concentration of A might change because of expansion or compression. Even a careful student however, might not realize that eq 2 is deterministic, so called because its solution leads to a formula, [A]/[A]0 = e‒kt, that determines the concentration of A at any time during the reaction from the concentration of A at the beginning of the reaction. And for a deterministic equation to be valid, a large number of molecules, usually greater than 1015, has to participate in the reaction. For most reactions studied in undergraduate physical chemistry courses, a large number of participating molecules is assumed. It is not surprising, therefore, that physical chemistry textbooks and many articles in this Journal emphasize the integration of deterministic equations, such as eq 2. Physical chemistry textbooks (1–8) focus on integrating relatively simple rate equations and simple coupled rate equations using analytical methods. Additionally, Andraos (9) has used Laplace transforms to analytically integrate the rate equations of a few more coupled rate equations. However, because the analytical integration of coupled rate equations can be challenging or even impossible, many articles in this Journal have used numerical integration. The approaches proposed in such articles include writing computer code in a programming language, such as C++ (10) or Basic (11), entering special commands in a computer algebra system, such as Mathematica (12–14) or Mathcad (15), using a dynamic system simulator, such as Stella (16) or Powersim (17), and using a stand-alone numerical integrator, such as Acuchem (18) or React (19). In the future, similar contributions might include using the dynamic system simulator Simulink (20) or the stand-alone numerical integrator Kintecus (21).1 Instructors sometimes choose to integrate simple rate equations in class analytically and coupled rate equations in the laboratory numerically.2 Although deterministic equations are appropriate for reactions studied in undergraduate physical chemistry, they have several limitations. First, they are difficult to solve for reactions that never attain steady-state and for reactions for which the 1146
temperature, pressure, or volume are not constant. Second, they are inappropriate when chemical instabilities are important, and when the range of rate constants and concentrations is large. And third, they contradict our knowledge that chemical reactions occur between individual species whose interactions follow probabilistic or stochastic laws. But there is an approach that overcomes the limitations of deterministic equations. It is stochastic simulation. The most important difference between the stochastic and the deterministic approaches is that deterministic approaches model the concentrations of the reacting species whereas stochastic approaches model the probabilities that species will react (the concentrations of the species are calculated after the number of molecules that have reacted are known.) The stochastic simulation of chemical reactions was introduced in this Journal in 1969 by Rabinovitch, with a series of paper-and-pencil exercises (22). Since 1969, several articles have emphasized the importance of stochastic modeling (23–25). That most of the articles on kinetic modeling in this Journal have focused on the numerical integration of rate equations (10–19) may be explained by the fact that such an integration is a convenient way to obtain concentration profiles of reactant or product species. Surprisingly though, there has been a similarly convenient way to perform stochastic simulations for chemical reactions since 1996: the Chemical Kinetics Simulator (CKS) program (26, 27). This program performs stochastic simulations of chemical reactions, and it is based, on one hand, on algorithms developed by Bunker et al. (28), and Gillespie (29, 30) in the 1970s, and on the other hand, on MSIM, a program written by Houle and Bunker (The program was first presented in this Journal in 1981; ref 31). CKS is free, easy-to-use, available on the Internet, and runs even on older computers. The goal of this article is twofold: to show how CKS can be applied to a few chemical and biochemical mechanisms that can be easily studied in the physical chemistry laboratory and to demonstrate that CKS can be a simpler and faster alternative to the analytical and numerical integration of rate equations. To achieve these goals, the article will briefly introduce the stochastic simulation of chemical reactions, and then it will describe the CKS-assisted modeling of four reaction mechanisms.3 These mechanisms are (i) two competing first-order reactions that serve as a model for the kinetic versus thermodynamic control of reaction products; (ii) an enzyme-catalyzed reaction that demonstrates the limitations of the commonly used steady-state approximation; (iii) the Lotka–Volterra mechanism, an oscillatory reaction whose study with Mathematica has been presented before in this Journal (12, 13); and (iv) the biochemical cycle that controls circadian rhythms in cyanobacteria, a mechanism that was studied numerically by Mehra et al. (32). With this approach, physical chemistry students will be introduced to the stochastic modeling of chemical reactions, the limitations of the steady-state approximation, oscillatory reactions, and to CKS, a versatile program that they can use in their future career.
Journal of Chemical Education • Vol. 85 No. 8 August 2008 • www.JCE.DivCHED.org • © Division of Chemical Education
On the Web 1.0
1.0
0.8
0.8
Concentration / M
[A]/[A]0
C
0.6
0.4
0.2
0.0
0.6
0.4
0.2
0.0 0
1
2
3
4
5
B
0
25
50
75
100
Time / s
kt Figure 1. Concentration profile of reactant A for elementary reaction 1. The filled circles correspond to the results of a 100-molecule Monte Carlo simulation and the solid line is plot of the integrated rate equation, [A]/[A]0=e∙kt. The results of a 1000-molecule Monte Carlo simulation coincide with the solid line.
Figure 2. Concentration profiles of products B and C according to a 104-molecule Monte Carlo simulation of mechanism 3, when [A]0 = 1 M, k1 = 1 s∙1, k∙1 = 0.5 s∙1, k2 = 0.1 s∙1, and k∙2 = 0.005 s∙1. Species B is formed faster (kinetically controlled product), but species C is more stable (thermodynamically controlled product.)
Stochastic Modeling of Chemical Reactions
Since both reactions in mechanism 3 are elementary and reversible, their equilibrium constants are defined as the ratio of their forward to their reverse rate constants, namely, K1 = k1/k‒1 and K2 = k2/k‒2. If molecule C is thermodynamically more stable than molecule B (K2 > K1), the final reaction mixture will contain more C than B. However, if molecule B is produced faster than molecule C (k1 > k2), the initial reaction mixture will contain more B than C. In other words, the composition of the final reaction mixture is thermodynamically controlled, but that of the mixture at the initial stages of the reaction, is kinetically controlled. This qualitative argument is not hard to make, but a quantitative argument based on the integration of the coupled rate equations is more challenging. However, a quantitative argument based on a Monte Carlo simulation is very simple with CKS. In CKS one can vary the values of k1, k‒1, k2, and k‒2 and observe which reaction product is the most abundant at various stages of the reaction. Figure 2, for example, shows the results of a 104-molecule Monte Carlo simulation for [A]0 = 1 M, k1 = 1 s‒1, k‒1 = 0.5 s‒1, k2 = 0.1 s‒1, and k‒2 = 0.005 s‒1, and K2 = k2/k‒2 = 20 > K1 = k1/k‒1 = 2. Figure 2 shows that if the reaction is terminated early (within the first 10 s), B is the main product (kinetic control), whereas if the reaction is allowed to reach completion, C is the main product (thermodynamic control). One can vary the values of the rate constants to see that products B and C can be present in equal amounts in the final mixture (when k‒2 = 0.05 s‒1) and that product B is both the thermodynamically and the kinetically controlled product (when k‒2 = k‒1 = 0.01 s‒1.)
The stochastic modeling of chemical kinetics has been described in detail in the past (22, 33, 34). As a brief illustration, consider the elementary reaction (eq 1), and assume that at t = 0 there are 100 molecules of species A. The Monte Carlo method, a common algorithm for stochastic simulation, starts by “labeling” 100 A-molecules with numbers from 0 to 99. After a time step, Δ t, the method generates a two-digit random number, and the A molecule that corresponds to this number is allowed to “react” and become a B molecule. The generated B molecule is labeled with the same two-digit number as the A molecule that produced it. After another time step, a new two-digit random number is generated and if the new number does not correspond to an existing A molecule, nothing happens. If the new number corresponds to an existing A molecule however, the reaction proceeds. After a sufficient number of time steps, the system has a new composition. Figure 1 shows the disappearance of species A due to eq 1, as predicted by both a Monte Carlo simulation (filled circles) and from the integrated rate equation, [A]/[A]0 = e‒kt (solid line). The Monte Carlo simulated concentration profile is more “noisy” than that predicted by the analytical expression when only 100 molecules participate in the simulation. However, when the number of molecules in the simulation is increased to 1000, the two concentration profiles coincide.4 In general, the noise of a stochastic simulation is inversely proportional to the square root of the number of molecules in the simulation. Example Mechanisms to Study with CKS Competing First-Order Reactions and Kinetic Versus Thermodynamic Control A reaction mechanism that demonstrates the kinetic versus thermodynamic control of reaction products is
A
k1 kļ1
B
and
A
k2 kļ2
C
(3)
Enzyme-Catalyzed Reactions and the Limitations of the Steady-State Approximation A common mechanism that describes enzyme-catalyzed reactions consists of the following two elementary steps
A E
k1 kļ1
B
k2
C E
(4)
where A is the substrate, E the enzyme, B the enzyme–substrate complex, and C the product of the reaction. Using the steady-
© Division of Chemical Education • www.JCE.DivCHED.org • Vol. 85 No. 8 August 2008 • Journal of Chemical Education
1147
On the Web C
0.5
0.3
0.2
0.1
0.8 0.7 0.6 0.5 0.4
B 0.0
X
0.9
0.4
Concentration / M
Concentration / M
1.0
0
2
4
6
Time / s
8
Y
0.3 10
0
20
40
60
80
100
Time / s
Figure 3. Concentration profiles of intermediate B and product C, according to a 104-molecule Monte Carlo simulation of mechanism 4, when k1 = 1 M∙1 s∙1, k∙1 = 0.1 s∙1, and k2 = 1.1 s∙1. The dashed line corresponds to the concentration of C, as predicted by eq 5, which is based on the steady-state approximation.
Figure 4. Concentration profiles of the intermediates X and Y according to a 106-molecule Monte Carlo simulation of the Lotka– Volterra mechanism, eq 6. The concentrations of both species oscillate with a period of approximately 20 s.
state approximation, Goodman (35) has shown that the concentration profile of the product, [C], is given by
Autocatalytic Reactions and Chemical Oscillations Reactions that show oscillations in the concentration of one or more species as a function of time are responsible for the formation of patterns on animal skin, the beating of mammalian hearts, and circadian rhythms in mammals (36, 37). Oscillating reactions were studied theoretically in the 1960s, and recently physical chemists have introduced the kinetics of these complicated reactions into the curriculum. A model oscillating reaction is described the Lotka–Volterra mechanism5 (34):
k2 0 t
0 k2 k1 ln k1 0
(5)
Goodman (35) also showed that when k1 = 1 M‒1 s‒1, k‒1 = 0.1 s‒1, k2 = 1.1 s‒1, and [A]0 = [E]0 = 0.5 M, the predictions of eq 5 differ quantitatively from those of a numerical integration of the rate equations. The dashed line in Figure 3 shows the concentration profile of [C], according to eq 5, and the solid line the predictions of a 104-molecule Monte Carlo simulation. The Monte Carlo predictions also differ quantitatively from those of eq 5. The difference between the predictions of eq 5 and those from numerical integration and stochastic simulation can be explained by the observation that the initial concentration of the intermediate B is comparable to that of the product C. Therefore, the assumption of the steady-state approximation, that the concentration of intermediate B is much lower than that of the product C, is not valid in this case. As expected, when the concentration of B diminishes, the predictions of the steady-state approximation agree more with those of the Monte Carlo simulation. However, the initial concentrations of the enzyme and the substrate used by Goodman, and reproduced above, are unrealistic because they are both too high and also because the concentrations of enzymes under physiological conditions are several orders of magnitude lower than that of substrates. When [E]0 is set to a more reasonable value, for instance such that [E]0/[A]0 = 10‒3, a stochastic simulation shows that the concentration of the intermediate is vanishingly small, and the steady-state approximation is valid. Using CKS students can therefore understand the physical interpretation of the steady-state approximation, appreciate its limitations, but also verify its appropriateness for enzymatic reactions.
1148
A X
X Y Y
k1 k2 k3
2X 2Y
(6)
P
Because species X is both a reactant and a product in the first reaction, as is species Y in the second reaction, the two reactions are autocatalytic in X and Y respectively. It can be shown that if the concentration of reactant A is kept constant during the reaction, the concentrations of the intermediates X and Y will exhibit steady oscillations with time. Such behavior has been demonstrated in this Journal with Mathematica in two articles (12, 13); however, the same task can be accomplished more easily with CKS. Figure 4 shows a 106-molecule Monte Carlo simulation of the Lotka–Volterra mechanism with k1 = 0.3 M‒1 s‒1, k2 = 0.6 M‒1 s‒1, k3 = 0.4 s‒1, [A]0 = 1.0 M, [X]0 = 0.8 M, [Y]0 = 0.8 M, and [P]0 = 0.0 M (the concentrations are the same as those used by Ferreira et al.; ref 12) Figure 4 shows that the concentrations of X and Y oscillate with a period of 20 s. Although a relatively large number of molecules (106) is necessary to obtain smooth concentration profiles in this case, the CKS simulation finishes in approximately five minutes.
Journal of Chemical Education • Vol. 85 No. 8 August 2008 • www.JCE.DivCHED.org • © Division of Chemical Education
On the Web 1.0
1.0
KaiB
Concentration / M
Concentration / NM
KaiB
0.5
0.5
KaiAC* 0.0
20
30
40
50
60
KaiAC* 0.0
70
Time / h
20
30
Time / h
40
50
Figure 5. Concentration profiles of protein KaiB and of the activated complex KaiAC* according to a 106-molecule Monte Carlo simulation of mechanism 7, at 25 oC. (The values of the rate constants are given in the text.) The concentrations of both species oscillate with a period of approximately 24 h.
Figure 6. Concentration of protein KaiB and of the activated complex KaiAC* according to a 106-molecule Monte Carlo simulation of mechanism 7, at 35 oC. (The values of the rate constants at 35 oC were calculated from those at 25 oC assuming Arrhenius behavior and an activation energy of 34 kJ/mol.) The concentrations of both species oscillate with a period of approximately 14 h.
Oscillatory Biochemical Reactions and the Circadian Rhythms of Cyanobacteria CKS can be used to model the system of oscillatory biochemical reactions that was studied in a 2006 article by Mehra et al. (32). Mehra et al. modeled the circadian rhythms in cyanobacteria as a sequence of reactions that involve three Kai proteins, KaiA, KaiB, and KaiC. The mechanism given in Figure 1 of the Mehra et al. article, is reproduced in the following scheme:
× 106 M‒1 h‒1, k5 = 4.0 h‒1, k6 = 0.90 h‒1, and k7 = 0.18 h‒1, and the initial concentrations of [KaiA] = 3.0 × 10‒6 M, [KaiB] = 1.0 × 10‒6 M, and [KaiC] = 3.5 × 10‒6 M. The result of a 106–molecule Monte Carlo simulation is given in Figure 5, which shows that the concentration KaiAC* increases when the concentration of KaiB decreases, a result that is consistent with the proposed mechanism. Figure 5 also shows that the concentrations of KaiB and KaiAC* oscillate with a period of approximately 24 hours. If one assumes, as Mehra et al. did in their article, that each step of the mechanism has an Arrhenius activation energy of 34 kJ/mol and calculates the pre-exponential factor, the circadian period can be found at a different temperature. Figure 6 shows the results of a 106–molecule Monte Carlo simulation at 35 °C, and indicates that the circadian period is reduced to 14 hours, a result that agrees with the predictions by Mehra et al.
KaiA KaiC KaiAC KaiAC* KaiA KaiC
KaiAC* KaiB
KaiABC* KaiBC* KaiC*
k1 k2 k3 k4 k5 k6 k7
KaiAC KaiAC* 2KaiAC* KaiABC*
(7)
KaiA KaiBC* KaiB KaiC* KaiC
In the above mechanism, KaiAC is a complex of KaiA and KaiC, and KaiAC* is the corresponding active complex. (The nomenclature is similar for KaiABC*, KaiBC*, and KaiC*.) According to the mechanism, proteins KaiA and KaiC bind to form the activated complex KaiAC*, which is released by protein KaiB. Mehra et al. used the 25 °C rate constant values of k1 = 100 M‒1 h‒1, k2 = 0.40 h‒1, k3 = 4.5 × 1011 M‒2 h‒1, k4 = 3.65
Summary The Monte Carlo simulation of chemical reactions, as implemented by CKS, introduces students to the stochastic nature of chemical reactions. It can also help students appreciate the physical meaning of the steady-state approximation, observe the kinetic versus thermodynamic control of chemical reaction, model oscillatory reactions, and apply kinetics to a complicated biochemical mechanism. Additionally, CKS provides a simpler, and sometimes faster, alternative to the analytical and numerical integration of rate equations. Notes 1. There are several programs on specific reaction mechanisms, like KinSimXP for mechanism A + B C + D (38, 39), Kinetica for reaction nA → B (40), a Java simulator for mechanism A B → C (41), and an Excel worksheet for enzyme kinetics (42).
© Division of Chemical Education • www.JCE.DivCHED.org • Vol. 85 No. 8 August 2008 • Journal of Chemical Education
1149
On the Web 2. To choose from the various numerical integration approaches, an instructor has to consider two factors: price (not all approaches are free) and ease-of-use (C++ and Basic require good programming skills; Mathematica and Mathcad require some programming skills; and dynamic system simulators and numerical integrators require some knowledge of the programs). 3. Toby and Toby (18) used a similar approach in their article on Acuchem, which is a stand-alone numerical integrator of rate equations. The Toby and Toby article can be used as a source of additional example mechanisms to be simulated with CKS. 4. Modeling the reaction with 100 or 1000 molecules in CKS takes less than a second, in both cases. 5. Although the mechanism is commonly used to describe predator–prey population dynamics, it does not describe any known chemical reaction (34).
Literature Cited 1. Levine, I. N. Physical Chemistry; McGraw-Hill: New York, 2002. 2. Atkins, P. W.; De Paula, J. Physical Chemistry, 8th ed.; Oxford University Press: New York, 2006. 3. Ball, D. W. Physical Chemistry; Thomson-Brooks/Cole: Pacific Grove, CA, 2003. 4. Silbey, R. J.; Alberty, R. A.; Bawendi, M. G. Physical Chemistry, 4th ed.; Wiley: Hoboken, NJ, 2005. 5. Laidler, K. J.; Meiser, J. H. Physical Chemistry, 2nd ed.; Houghton Mifflin: Boston, 1995. 6. Metiu, H. Physical Chemistry. Kinetics; Taylor and Francis Group: New York, 2006. 7. Raff, L. M. Principles of Physical Chemistry; Prentice Hall: Upper Saddle River, NJ, 2001. 8. McQuarrie, D. A.; Simon, J. D. Physical Chemistry: A Molecular Approach; University Science Books: Sausalito, CA, 1997. 9. Andraos, J. J. Chem. Educ. 1999, 76, 1578–1583. 10. Pavlis, R. R. J. Chem. Educ. 1997, 74, 1139–1140. 11. Macomber, R. S.; Constantinides, I. J. Chem. Educ. 1991, 68, 985–988. 12. Ferreira, M. M. C.; Ferreira, W. C. J.; Lino, A. C. S.; Porto, M. E. G. J. Chem. Educ. 1999, 76, 861–866. 13. Francl, M. M. J. Chem. Educ. 2004, 81, 1535. 14. Mulquiney, P.; Kuchel, P. W. Modelling Metabolism with Mathematica; CRC Press: Boca Raton, FL, 2003. 15. Harvey, E.; Sweeney, R. J. Chem. Educ. 1999, 76, 1309. 16. Steffen, L. K.; Holt, P. L. J. Chem. Educ. 1993, 72, 991–993. 17. Ricci, R. W.; Van Doren, J. M. J. Chem. Educ. 1997, 74, 1372– 1374. 18. Toby, S.; Toby, F. S. J. Chem. Educ. 1999, 76, 1584–1590. 19. Manka, M. J. Chem. Educ. 2000, 77, 165. 20. Simulink, version 6.6.1; The MathWorks, Inc.: Natick, MA, 2007.
1150
21. Kintecus, version 3.90; Vast Technologies Development, Inc: Lansdowne, PA, 2007. 22. Rabinovitch, B. J. Chem. Educ. 1969, 46, 262–268. 23. Dixon, D. A.; Shafer, R. H. J. Chem. Educ. 1973, 50, 648–650. 24. Freeman, G. R. J. Chem. Educ. 1984, 61, 944. 25. Mira, J.; Fenandez, C. G.; Urreaga, J. M. J. Chem. Educ. 2003, 80, 1488–1493. 26. Hinsberg, W. D.; Houle, F. A. Chemical Kinetics Simulator, version 1.01. http://www.almaden.ibm.com/st/computational_science/ ck/ (accessed Mar 2008). 27. Houle, F. A.; Hinsberg, W. D. Real-World Kinetics via Simulations. In Annual Reports in Computational Chemistry, Vol. 2; Spellmeyer, D. C., Ed.; Elsevier: Amsterdam, 2006. 28. Bunker, D. L.; Garrett, B.; Kleindienst, T.; Long, G. S., III. Combust. Flame 1974, 23, 373–379. 29. Gillespie, D. T. J. Comput. Phys. 1976, 22, 403–434. 30. Gillespie, D. T. J. Phys. Chem. 1977, 81, 2340–2361. 31. Houle, F. A.; Bunker, D. L. J. Chem. Educ. 1981, 58, 405– 407. 32. Mehra, A.; Hong, C. I.; Shi, M.; Loros, J. J.; Dunlap, J. C.; Ruoff, P. PLoS Comput. Biol. 2006, 2, e96. 33. Connors, K. A. Chemical Kinetics: The Study of Reaction Rates in Solution; VCH: New York, 1990. 34. Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics; Prentice Hall: Englewood Cliffs, NJ, 1989. 35. Goodman, J. M. J. Chem. Educ. 1999, 76, 275–277. 36. Atkins, P. W. Atoms, Electrons, and Change; W. H. Freeman: New York, 1991. 37. Prigogine, I.; Kondepudi, D.; Dewel, G. Chemistry Far from Equilibrium: Thermodynamics, Order and Chaos. In The New Chemistry; Hall, N., Ed.; Cambridge University Press: Cambridge, 2001. 38. Allendoerfer, R. D. J. Chem. Educ. 2002, 79, 638. 39. Allendoerfer, R. D. J. Chem. Educ. 2003, 80, 110. 40. Vera, L.; Ortega, P.; Guzman M. J. Chem. Educ. 2004, 81, 159. 41. Haustedt, L. O.; Goodman, J. M. J. Chem. Educ. 2003, 80, 839. 42. Bruist, M. F. J. Chem. Educ. 1998, 75, 372–375.
Supporting JCE Online Material
http://www.jce.divched.org/Journal/Issues/2008/Aug/abs1146.html Abstract and keywords Full text (PDF) Links to cited URLs and JCE articles Supplement
Student handouts
Instructor notes
The CKS reaction files of the mechanisms
Journal of Chemical Education • Vol. 85 No. 8 August 2008 • www.JCE.DivCHED.org • © Division of Chemical Education
