Computational Models of Thermal Cycling in Chemical Systems

Computational models of chemical systems provide clues to counterintuitive interactions and insights for new applications. We have been investigating models of chemical reaction systems under forced, thermal cycling conditions and have found that some hypothetical processes generate higher yields under thermal cycling than under single, fixed temperature conditions. A simple kinetic model of an actual process, the two-temperature polymerase chain reaction that replicates DNA, is used to simulate the important features of a chemical system operating under thermal cycling. This model provides insights into the design of other chemical systems that may have important applications in chemistry, biochemistry and chemical engineering.


Introduction
Computational models of chemical systems provide insight into new reaction processes and mechanisms.The models described here were inspired by a mathematical concept from game theory, known as Parrondo's paradox, whereby two losing games can produce a winning outcome if appropriately alternated (Abbott D. 2010;Harmer G. P. & Abbott D. 1999;Martin H. & von Baeyer H. C. 2004; The University of Adelaide, School of Electrical Engineering, Parrondo's Paradox Official Website).To mimic alternating game strategies, we have been investigating the effect of thermal cycling on chemical reaction systems, and have indeed discovered that certain chemical processes can generate significantly higher yields of product under thermal cycling than under single, fixed-temperature conditions (Osipovitch D. C., Barratt C. & Schwartz P. M. 2009).While there are few actual chemical processes that are carried out under thermal cycling, there are several examples of reactions carried out in microreactors that have been found empirically to occur more rapidly under fast, forced thermal cycling (Brandner J. J., Emig G., Liauw MA & Schubert K. 2004;Hansen H. A., Olsen J. L., Jensen S., Hansen O. & Quaade U. J. 2006; Jensen S., Thorsteinsson S., Hansen O. & Quaade J. J. 2008; Luther M., Brandner J. J., Schubert K., Renken A. & Kiwi-Minsker L. 2008; Luther M., Brandner J. J., Kiwi-Minsker L., Renken A. & Schubert K. 2008; Silveston P. L. & Hudgins R. R. 2004).To motivate a more general investigation of the utility of thermal cycling, we developed a simple kinetic model of the two-temperature polymerase chain reaction (2T-PCR) as an actual chemical system that replicates DNA.
The polymerase chain reaction is a fundamental tool in molecular biology, which permits amplification of double stranded DNA (dsDNA) in a thermal cycling process (Kramer M. F. & Coen D. M. 2001;Lo Y. M. D. &. Chan K. C. A. 2006).Traditional PCR involves three temperature steps: denaturation of dsDNA at a high temperature to two complementary, single-stranded DNA (ssDNA) templates (a 5'-3'ssDNA and a 3'-5'ssDNA); annealing of specific primers to the ssDNA templates at an intermediate temperature; and extension and replication of DNA at a lower temperature catalyzed by a DNA polymerase in the presence of deoxyribonucleotide substrates (dATP, dGTP, dCTP and TTP).Under appropriate conditions, the traditional three-temperature protocol for PCR can be replaced by a two-temperature cycle (Cha R. S. & Thilly W. G. 1993;Horton J. H., Hagen M. D. & Ko K. S. 1994;Neuzil P., Zhang C., Pipper J., Oh S. & Zhuo L. 2006).With two-temperature PCR, there is a cyclic process between a short period at high temperature (90-95 o C) to denature dsDNA into two ssDNA templates and then a longer period at 60-75 o C in which both annealing and extension of complementary primers occur.Every cycle of the PCR reaction theoretically doubles the amount of dsDNA, i.e. is 100% efficient so that with n temperature cycles, 2 n copies of the original dsDNA are generated, n being determined by how much substrate is present: doubling falls off as substrates (primers and/or deoxyribonucleotides) are depleted.
The goal of creating chemical systems that operate like PCR has eluded investigators.PCR as a chemical system is driven by forced thermal cycling and is unlike spontaneously oscillating reactions (such as the BZ reaction or Turing systems).Like PCR, some very interesting self-replicating, autocatalytic systems have been devised (Beutel K. M. & Peacock-Lopez E. 2007;Kamioka S., Ajami D. & Rebek J. 2010;Patzke V. & von Kiedrowski G. 2007;Rebek J. 1992;Robertson A., Sinclair A. J. & Philp D. 2000;Sievers D. & von Kiedrowski G. 1998;Vidonne A. & Philp D. 2009;Wang B. & Sutherland I. O. 1997;Wintner E. A., Conn M. M. & Rebek J. 1994), but these reactions have not been investigated under thermal cycling conditions.Yoon and Mirkin recently reported a PCR-like cascade reaction using a rhodium coordination complex (Yoon H. J. & Mirkin C. A. 2008); however, as pointed out by D. Blackmon, "the sigmoidal profiles observed in that case may be attributed to an in situ catalyst activation process …. and is not related either to self-replication or to PCR-type amplification" (Blackmond D. G. 2009).We suggest that computational modeling of the kinetic features of two-temperature PCR may provide essential clues to creating actual chemical processes that, when carried out under thermal cycling, would possess the advantages of a self-replicating chemical system.

The Chemical System under Investigation -2T-PCR
A simplified PCR system that we are investigating is represented by the following reaction equations: in which ssDNA1 and ssDNA2 designates the two complementary ssDNA templates of dsDNA; SUB represents the substrates, both the deoxyribonucleotides and primers; E act represents the active enzyme that functions at the lower temperature and E inact is the inactive form of the enzyme at the higher temperature.

Parameter selection and assumptions
For the model systems, we chose parameters that reflect typical conditions for real-time, two temperature PCR (Cha R. S. & Thilly W. G. 1993;Horton J. H., Hagen M. D. & Ko K. S. 1994;Neuzil P., Zhang C., Pipper J., Oh S. & Zhuo L. 2006).A time course for the thermal cycling was chosen in which the PCR commenced at a high temperature (370K, 97 o C) to denature the initial dsDNA.Then, cycles were chosen such that a period of 75 sec.at a low temperature (330K, 57 o C) was followed by a short time, 15 sec., at the higher temperature.The temperature protocol is typical for 2T-PCR.The simulations were analyzed for 25 cycles.SUB starts at a high concentration and becomes limiting for synthesis of dsDNA.In the initial models, the concentration of dsDNA is taken as 1 x 10 -5 M. The models make several assumptions including no loss in enzyme activity over time and no side reactions (such as primer-dimer interactions).This allowed us to focus on the core reactions of PCR as a chemical system.

Kinetic Description and Theoretical Analysis of the Chemical System
Kinetic parameters, A (the Arrhenius constant) and E a (the activation energy), which yielded a working model, are shown in Table 1.Rate constants, k 1 through k 6 at 330K and 370K, were obtained in the usual way using the Arrhenius equation.The initial concentrations are set to be: The governing differential equations for the simplified PCR reactions shown above are: To simplify the notation, we define variables X,Y,A,B,C , representing the concentrations at time t of the various chemicals: Then the governing equations for X, Y, A, C and B now read: Eqs.
(1) and (2), governing the enzyme concentrations E inact and E act are easy to solve to give For the k 1 and k 2 used in the present investigation, Y fairly quickly (within a few seconds) flatlines after a change in temperature from 330K to 370K or vice versa, to the "constant" value of Y = (X 0 + Y 0 )/(1+k 2 /k 1 ).
Eqs.( 3),( 4) and ( 5) are coupled and nonlinear and thus difficult to solve in general.However, they do yield to analysis in three distinct regions, viz: a) Constant temperature, T=330K or T=370K, for A~constant.
b) The "doubling" region, in which temperature oscillates between 330K and 370K and also A~constant.This consists of the first ~10 cycles, i.e. 0<t<900 secs.The cycle period is taken to be 90s, of which the first 15s is at T 1 =370K, the remaining 75s at T 2 =330K.
c) The "plateau" region which begins after ~15 cycles, by which time the substrate has depleted to A~0.

Constant temperature region (a)
Numerical solution reveals that throughout the integration region, A~1, so that K= k 5 AY is essentially constant.Also, we find that KB >> k 4 B 2 .Thus the equations for C and B are to very good approximation linear and easily solved, yielding, for .Thus, over the period of interest, 0 < t < 2340s (26 cycles), C increases steadily from the initial value C 0 = 10 -5 to 1.12 x 10 -5 .
For constant T =370K, we obtain C ~ (e -2.96t + 2 x 10 -3 e 0.00296t ) C o .In this case, the fast exponential causes C to drop precipitously from the initial C 0 = 10 -5 to a value C ~ 2 x 10 -8 in the initial few seconds.Thereafter, the slow exponential dominates the solution and C increases steadily with time constant (1/3) x 10 3 s to a value C ~ 2 x 10 -5 at t = 2340s.In fact, this is a slight overestimate because, as indicated by numerical solution of the equations, at T=370K, A has depleted from A=0.5 to A=0.487 at t=2340 secs, and C is actually ~ 1.73 x 10 -5 .Either way, an important point, as will be confirmed below, is that under constant temperature conditions the production rate for C is significantly lower than that obtained when the temperatures are cycled.

The "doubling" region, 0<t<900 secs (b)
For the first 10 cycles or so, the concentration, A, of substrate, remains almost fixed at A=0.5.Also, as in a) we find numerically that KB >> k 4 B 2 and because of the quick flatlining of Y at each temperature, K=k 5 AY is also constant throughout each partial cycle.The general solution of the kinetic equations is: where, again, We here note that C 0 , B 0 represent the values of the concentrations C and B at the beginning of each T=constant partial cycle; they are the values of C and B acquired in the previous partial cycle.

The "plateau" region (>20 cycles) (c)
The Eqs 3-5 yields: so that when the system enters the "plateau" region defined by A~0, which occurs (for the parameters used in this paper) after ~20 cycles, we have thereafter C + B = σ, where the constant σ = C o + B o + A o /2 = 10 -5 + 0 + 0.25.
Using this constraint between C and B, the dC/dt equation now simplifies (recall A~0) to It should be noted that when using the above solution, t is again reset to t = 0 at the beginning of each "partial cycle", i.It is found that the above solutions a) b) c) reproduce accurately the behavior of the system, as revealed by a Runge-Kutta solution of the governing equations and/or by use of computational methods (see below).By way of summarizing the above analysis, we have been able to solve the nonlinear governing kinetic equations in the special cases a) T = constant = 330K or 370K, b) the "doubling" period, 0<t<900 secs, during which period the concentration of dsDNA approximately doubles in each 370-330 cycle, and c) the "plateau" region the onset of which occurs when the substrate A has depleted (A=0) and in which region the concentration of dsDNA oscillates between 0.00736 M and 0.208 M. What is clear though, whether from the numerical or analytic solutions, is that cycling the temperatures between two values leads to substantially higher production of dsDNA than maintaining the temperature fixed at either temperature.

Deterministic methods
Kintecus (version 3.953) was developed by James Ianni as a powerful simulation program for chemical dynamics (Ianni J. C. 2009); it is free for academic use.As a deterministic, Arrhenius-based program the inputs include: the reaction steps, energies of activation (E a ) and Arrhenius constant, A. The program assumes elementary reaction steps and numerically solves the differential equations for the related rate laws.The concentrations of participating species are calculated and displayed over time at either a fixed temperature or under varying temperature conditions.

Deterministic model of 2T-PCR
A model of two-temperature PCR was developed using Kintecus, a computational package based on deterministic kinetic equations.A model was designed with kinetic parameters, initial concentrations of reactants, a temperature-sensitive catalyst and an appropriate protocol for the temperature cycles (Table 1 and Fig. 1a.).E act represents the active enzyme that functions at 330K; E inact is the inactive form of the enzyme at the higher temperature (370K).SUB represents the substrates, both the deoxyribonucleotides and primers; SUB starts at a high concentration and becomes limiting for synthesis of dsDNA.The two ssDNA templates are designated as ssDNA1 and ssDNA2.The output from the model describes a pattern of replication of dsDNA that is typical for experimental results obtained using real-time, two-temperature PCR (Fig. 1b.); there is replication of DNA during the low temperature part of the cycle, denaturation of dsDNA during the high temperature phase, an exponential phase in which there is doubling of dsDNA during each cycle and a plateau phase in which replication is approached as availability of substrates becomes limiting.
We used the deterministic model to simulate the generation of standard curves.The initial concentration of dsDNA was varied from 1 x 10 -5 M to 1 x 10 -8 M (Fig. 2.).The x-axis is time but can be re-calculated as cycle number.The plot shows the peak DNA concentration at the each cycle.The amplification efficiency is related to the slope of the curve in the log-linear exponential phase; theoretically, the efficiency of PCR is 100%, i.e. the concentration of dsDNA is doubled at each thermal cycle.In the models (Fig. 1b and Fig. 2) the calculated efficiency is 100%.

CKS (Chemical Kinetics Simulation
) is a program provided by IBM for developing stochastically based models (Hinsberg W. & Houle F. Chemical kinetics simulator 1.01).Inputs include the reaction steps, kinetic constants for the temperature dependent reactions, the initial concentrations of reactants and time course for the thermal cycling; stochastic parameters including the random number seed, the time course for the simulation and the number of "molecules" in the model.

Stochastic model of 2T-PCR
A kinetic model of 2T-PCR using a stochastic chemical kinetics simulation program, CKS (Hinsberg W. & Houle F. Chemical Kinetics Simulator 1.01), is shown in Figure 3.The input parameters for this model were those used from the deterministic model described above (Table 1 and Fig. 1).The stochastic model of 2T-PCR is equivalent to the above deterministic model for 2T-PCR.

Discussion
Systems chemistry is a new and exciting area that explores the interaction of molecules in multi-reaction systems often displaying behavior not expected from individual components.The models of 2T-PCR recapitulate the real behavior of this system; from the numerical or computational methods, cycling the temperatures between two values leads to substantially higher production of dsDNA (target product) than maintaining the temperature fixed at either temperature.Viewing 2T-PCR as a chemical system provides insight into features of more general chemical processes that might use thermal cycling advantageously.The model highlights two catalytic events, an enzyme catalyzed reaction and an autocatalytic process in which DNA acts as a template for its replication.The interaction of these components is an essential feature of the non-linear behavior of this chemical system under thermal cycling.To extend this to a more general pattern, we suggest the following analogy: where X and Y are active and inactive catalytic forms, C is a complementary polymer of two templates, B1 and B2 and A is substrates.
Our simple kinetic model of thermal cycling is likely to suggest many possible applications where thermal cycling would be advantageous for product formation in a chemical system.One possible application of our model is the replication of RNA in a prebiotic environment (  Figure 2. dsDNA Standard Curves.The standard curves were generated from the model described in Table 1 and Fig. 1.The curves correspond with initial concentration of dsDNA: ♦ 1 x 10 -5 M, ■ 1 x 10 -6 M, ▲ 1 x 10 -7 M, • 1 x 10 -8 M Figure 3. Output from the Stochastic 2T PCR Model.Simulation conditions includes: total number of events = 9 x 10 8 ; total number of particles = 2 x 10 6 ; equilibrium enabled with recording intervals of 200 events.Note: early points are unreliable due to low initial concentrations of dsDNA, which is typical of a stochastic model e. when k 3 and k 4 change values, and C 1 is the value of C at the beginning of the next partial cycle.For the parameters used in our model, we obtain, for T = 330K, β ~ 0.25, α ~ 0.00675 and for T = 370K, β ~ 4.25, α ~ 4.2426.Notice that, for T=370K, the large value of αk 4 means that within just a few seconds, C flatlines to C = βα ~ 0.00736, the value C 1 at the beginning of the T=330K partial cycle, at the end of which (75s later) C has the value of 0.208, the maximum value that it achieves.
Ferris J. P., Joshi P. C., Wang K-J., Miyakawa S. & Huang W. 2004; Ferris J. P. 2006; Ferris J. P. & Delano J. W. 2008; Joshi P. C., Aldersley M. F., Delano J. W. & Ferris J. P. 2009; Lincoln T. A. & Joyce G. F. 2009).Our models suggest that catalytic surfaces and substrates might interact with complementary RNA templates under thermal cycling conditions to replicate RNA.deRoss has proposed the hypothesis that an PCR-like process might have replicated dsRNA in a non-enzymatic fashion as an early chapter in the RNA World hypothesis (de Roos D. G. 2007).Kinetic models of chemical systems operating under thermal cycling may provide clues to creating innovative applications in chemistry and chemical engineering.Table 1.Description of kinetic parameters for model of 2T-PCR E act and E inact represent the concentrations of the catalyst, i.e.DNA polymerase, dsDNA the target molecule, and ssDNA1 and ssDNA2 the templates.SUB represents the four deoxynucleotides and primers for PCR.The kinetic parameters (A and E a ) are used by the Arrhenius equation (k = A e (-Ea/RT) using R = 0.008315 kJ/mmol) to calculate rates and determine concentrations of components over time.The initial concentrations are: [dsDNA] o = 1.0x10 -5 M, [ssDNA] o = 0.0M, [SUB] o = 0.5M, [E inact ] o = 0.0M, [E act ] o = 1.0x10 -3 M.

Figure 1 .
Figure 1.Deterministic, Two-temperature PCR Model.(a) The temperature versus time profile.The temperature profile is typical of PCR, starting at the higher temperature, 370K, and then cycling between 330K (75sec) and 370K (15sec) for 25 cycles.(b) Formation of dsDNA (log scale) versus time at ( ) for 330K and ( ) for 370K and ( ) for thermal cycling as in Fig 1a.Initial dsDNA is 1 x 10 -5 M