Consider a generic chemical reaction such as:

 \(\require{\mhchem} \ce{aA +bB <=>[k_F][k_R] cC + dD}\)

where A, B, C, D are chemical species, while a, b, c, d are integers (the stoichiometry coefficients.)  k_F and k_R are, respectively, the forward and reverse rate constants, as detailed below. 

A and B are generally referred as the reactants, while C and D are the products of the reaction.  However, to be noted that the reaction goes in both directions at the same time – imagine two adjacent rooms, where some people go from one to the other, while others go in the opposite direction.

As customary, we'll use square brackets, as in  \([A]\) , etc, to indicate the concentrations.  Generally, those are functions of time, \([A](t)\) , but we'll often omit the \((t)\)

If the concentrations change with time, the reactants as well as the reaction products must change in lockstep to satisfy the stoichiometry (i.e. maintain the relative balances.)

Upon shifting our attention to large numbers of molecules in a liquid solution, one typically starts using real number and concentrations; for example, at time t₀,  [A](t₀) = 12.425 billion molecules per liter.

Using primes to indicate the time derivatives (i.e. rate of change) of the Concentrations functions, and omitting the (t) parts, the above example generalizes to:

\({1 \over a} [A]' = {1 \over b} [B]' = - \ {1 \over c} [C]' =- \ {1 \over d} [D]' \)

The expressions with the [A] and [B] concentrations indicate the rate of change of the reactants; the expressions with the [C] and [D] concentrations indicate the rate of change of the reaction products.

Solving the Reaction Equations

From the Kinetic Theory of chemical reactions, the foward rate is proportional to \([A]^{O_A} * [B]^{O_B} \) , and the proportionality factor is called forward rate constant.

The exponents O_A and O_B are, respectively, the orders of the reaction relative to A and B.  The orders have to do with the detailed mechanism of the reaction, and they are measured experimentally, or estimated by means of Molecular Dynamics computations.  IMPORTANT: the orders of the reaction are NOT necessarily the same as the stoichiometry coefficients!

\(\require{\mhchem} \ce{aA +bB <=>[k_F][k_R] cC + dD}\)

Intuitively, the higher the concentrations of A and B are, then the higher the odds that they will collide and move the reaction forward.

\(Forward \ Rate = k_F * [A]^{O_A} * [B]^{O_B}\)

similarly:

\(Reverse \ Rate = k_R * [C]^{O_C} * [D]^{O_D}\)

Intuitively, the higher the concentrations of C and D are, then the higher the odds that they will collide and move the reaction in the reverse direction.

Note the somewhat unfortunate nomenclature: the forward or reverse rates should not be confused with the rate constants (k_F and k_R )

We'll define:  \(Delta\ Rate = Forward\ Rate - Reverse\ Rate\).  The "Delta Rate" will measure the imbalance towards the forward direction of the reaction (i.e., the tendency of the reaction to proceed forward.)

Substituting from above:

\(Delta\ Rate = \left( k_F * [A]^{O_A} * [B]^{O_B} \right) - \left( k_R * [C]^{O_C} * [D]^{O_D} \right)\)

Putting it all together, for the reaction products:

\({1 \over c} [C]' = {1 \over d} [D]' = Forward \ Rate -\ Reverse\ Rate = Delta\ Rate \)

and, similarly for the reactants:

\({1 \over a} [A]' = {1 \over b} [B]' = Reverse \ Rate -\ Forward \ Rate = -\ Delta\ Rate \)

For the reactants, for example for A, and changing the notation to \({d \over dt }[A](t)\) to emphasize the ordinary differential equation (ode):

\({d \over dt }[A](t) = -a * \ Delta\ Rate = -a * \left\{ k_F * [A](t)^{O_A} * [B](t)^{O_B} -k_R * [C](t)^{O_C} * [D](t)^{O_D} \right\}\)

For the reaction products, for example for C, the sign is reversed:

\({d \over dt }[C](t) = -c * \ Delta\ Rate = c * \left\{ k_F * [A](t)^{O_A} * [B](t)^{O_B} -k_R * [C](t)^{O_C} * [D](t)^{O_D} \right\}\)

The initial conditions are the given concentrations at time t=0 :  \([A](0) , [B](0)\) , etc.  We'll refer to those values as A₀, B₀, etc

One straightforward way to solve the above ode is to use the forward Euler method to start with the values at t=0, namely A₀, B₀, C₀, D₀, and then compute all the values at a later time step :

If we introduce subscripts to indicate values at discrete times, and define: 

\(\text{delta_forward}(n) = \left( k_F * [A]_n^{O_A} * [B]_n^{O_B} - k_R * [C]_n^{O_C} * [D]_n^{O_D} \right)\)

then, for the reactants:

\([A]_{n+1} = [A]_n - (\Delta t) * a * \text{delta_forward}(n)\)

\([B]_{n+1} = [B]_n - (\Delta t) * b * \text{delta_forward}(n)\)

and for the reaction products the sign flips:

\([C]_{n+1} = [C]_n + (\Delta t) * c * \text{delta_forward}(n)\)

\([D]_{n+1} = [D]_n + (\Delta t) * d * \text{delta_forward}(n)\)

where \(\Delta t\) is the desired step size, and n = 0, 1, 2, ....

How Many Equations Do We Have?

At the end of last section, it might appear that we're lugging around 4 equations but, in fact, if we compute the change of concentrations of any species, then the other ones are simply determined by the stoichiometry.

For example, let's say we know \(\color{blue}{[A](t)}\) at some time t , as well as all the initial concentrations  A₀, B₀, C₀, D₀.

Then:
\([B](t) = B_0 + \Delta B = B_0 + b { \Delta B \over b} = B_0 + b { \Delta A \over a} = B_0 + { b \over a} (\color{blue}{[A](t)} - A_0) = { b \over a} \color{blue}{[A](t)} + (B_0 - {b \over a} A_0)\)

So, [B](t) is just a linear function of [A](t)

For [C](t) and [D](t) there's a sign flip, going from reactant to product of reaction:

\([C](t) = C_0 + \Delta C = C_0 + c { \Delta C \over c} = C_0 \color{red}- c { \Delta A \over a} = C_0 - { c \over a} (\color{blue}{[A](t)} - A_0) = -{ c \over a} \color{blue}{[A](t)} + (C_0 + {c \over a} A_0)\)

\([D](t) = D_0 + \Delta D = D_0 + d { \Delta D \over d} = D_0 \color{red}- d { \Delta A \over a} = D_0 - { d \over a} (\color{blue}{[A](t)} - A_0) = -{ d \over a} \color{blue}{[A](t)} + (D_0 + {d \over a} A_0)\)
 

Numerical Instability

Numerical instability is a bane of solving ODE's! One advantage of dealing with variables that represent chemical concentrations, is that they must be always be non-negative. This constraint can be utilized to help detect instability: chemical concentrations assuming negative values may be used as an instability warning - not enough to catch all instances of instability, but helpful nonetheless.

3 scenarios may be distinguished.  Illustrative examples are shown; "Delta_1" and "Delta_2" are the changes to the concentration of a particular chemical, arising, respectively, from two separate reactions.

Detection and Automated Correction of Some Instability

At present, Life123 uses negative concentrations, plus multiple heuristics, as a proxy for numerical instability.  This appears to catch "the worst offenders", but not all cases, of instability.

Whenever an individual reaction - or group of reactions - causes any chemical concentration to go negative (under any of the 3 scenarios above), a python exception of custom type ExcessiveTimeStep (defined in class ReactionDynamics) is raised.

That exception gets caught in the method ReactionDynamics.reaction_step_orchestrator() - which results in the following:

1. the (partial) execution of the last step gets discarded
2. the step size gets reduced by a given amount (which the user may modify)
3. the run (simulation of reactions) automatically resumes from the previous time

The detection and automatic correction of negative concentration values is showcased in the experiment negative_concentrations_1

Macromolecules

Version Beta 28 introduced initial support for macromolecules.

Occupancy of binding sites on a macromolecule, as a function of Ligand Concentration, can now be modeled. (For example, Transcription Factors binding to DNA.)

Example:

In upcoming versions, the fractional occupancy values will regulate the rates of reactions catalyzed by the macromolecules...

Examples of usage may be found under Experiments (in the "Single-compartment Reactions" section.)