Perspective Chapter: A New Kind of Chemical Kinetics

Abstract

After reviewing the results for the Michaelis–Menten enzyme mechanism, both from the usual deterministic coupled differential equations of Chemical Kinetics and from the stochastic model of Gillespie, the first conclusion is that both, the smoothness of the concentration changes from the first model and the chaotic concentration fluctuations from the second model, are implied by the kind of mathematics used. I consider that neither the smoothness nor the chaotic fluctuations of the concentrations are real facts. In the new model developed here, the timeline is a sequence of equally spaced time points, at which concentration changes can occur; the time interval, τ, is to be selected by analyzing the results. The coupled algebraic equations resulting from the linear integration of the differential equations of the first model, instead of being solved, are used to extract the constraints of the Objective Function whose minimization renders the collective optimum values of the concentrations along the reaction path. One advantage of this model is that by adding the conservation of mass as an additional constraint in the Objective Function, a self-organized behavior is observed in this prebiotic system along with the chemical dynamics, which I consider real.

1. Introduction

Fortunately, if you need or like to learn about how the science is constructed, there are good introductions to the history of Chemical Kinetics; two of those are A Brief Introduction to the History of Chemical Kinetics, by Ptáček, Opravil, and Šoukal [1] and Three Waves of Chemical Dynamics, by Gorban and Yablonsky [2].

The science of Chemical Kinetics is not connected with the establishment of the equations expressing the transformation and combination of several compounds in order to create new ones, but with the establishment of the equations that express the velocities at which those transformations take place as well as the establishment of the steps involved in such transformations.

In a solution of sugar, sucrose, fructuose, sacarose, and so forth, by measuring the optical property of the transmitted light called polarization, the concentration of the reactant can be deduced; in 1850, the German chemist Ludwig Whilhelmy determined that in a solution of sucrose (Z) and some acid (S), the instantaneous velocity of decrease of the sucrose (dZ(t)/dt) depends on the product of both instantaneous concentration (Z(t)) and the constant concentration of S; the first chemical kinetic equation ever written is:

$$\mathrm{dZ}\left(\mathrm{t}\right)/\mathrm{dt}=‐\mathrm{k}\mathrm{Z}\left(\mathrm{t}\right)\mathrm{S}$$E1

Separating variables and integrating, we get: log Z(t) = − k S t + C, where C is the integration constant and k is the so called reaction rate constant, which are determined experimentally from the horizontal interception and the vertical interception in a plot of log Z(t) versus t.

Later on, the chemist Peter Waage and his brother-in-law, the mathematician Cato Maximilian Guldberg, from Norway, worked together from 1864 to 1879 to establish what is known the law of mass action: the rate of any chemical reaction is proportional to the product of the masses (concentrations) of the reacting substances, with each mass raised to a power equal to the coefficient that occurs in the chemical equation.

According to Gorban and Yablonsky [2], the first wave of chemical kinetics was concluded in 1901 when Jacobus H. van’t Hoff was awarded the Nobel Prize in Chemistry by the (re)discovery of the laws of chemical dynamics—the (re) is mine—and osmotic pressure in solutions; the second wave is the Semenov–Hinshelwood; in 1956, they were awarded the Nobel Prize in Chemistry for their research in the mechanism of chemical reactions, and the third wave started in 1965 when Aris proposed a new synthesis and understanding of the chemical reaction networks, which has been of paramount importance for the latest advances in chemical engineering.

Currently, we have deal only with the first kind of chemical kinetics, the deterministic model based on extracting from the structure of the chemical reactions a system of coupled differential equations that have to be solved and interpreted.

There is also a second kind of chemical kinetics, which is not deterministic, but stochastic, which means that uses probabilistic arguments to interpret what happens during the course of the chemical reactions. It was the American Scientist Gillespie [3] who in 1976 presented this mathematical model, which today is being used, together with the deterministic model, in the attempts to model the chemical evolution of a whole cell [4].

A comparison of those two models and results of chemical kinetics is presented in Section 2 based on the classic Michaelis–Menten enzyme mechanism; in Section 3, the new kind of chemical kinetics will be introduced.

2. The two old kinds of chemical kinetics

I will introduce the deterministic approach by presenting the classic Michaelis–Menten enzyme mechanism (from the book Chemical BioPhysics [5], page 49):

$$\begin{array}{l}\mathrm{A}+\mathrm{E}\left(\mathrm{k}‐1=\mathrm{km}1\right)<‐‐‐‐>\left(\mathrm{k}+1=\mathrm{kp}1\right)\mathrm{C}\\ \mathrm{C}\left(\mathrm{k}‐2=\mathrm{km}2\right)<‐‐‐‐>\left(\mathrm{k}+2=\mathrm{kp}2\right)\mathrm{A}+\mathrm{B}\end{array}$$E2

E is the enzyme that converts the substrate A into the product B, through the participation of the substrate–enzyme complex C; the rate constants {km1, kp1, km2, kp2} are experimentally determined.

To apply the law of mass action, we should consider a left to right or first reaction and a right to left or second reaction, although both reactions occur simultaneously. Consider the first reaction; if the velocity of the direct (left to right) reaction is defined as (vd = kp1 a e), where a and e are the molar concentrations of A and E, respectively, and the velocity of the reverse (left to right) reaction is defined as (vr = km1 c), where c is the molar concentration of C, regardless of the initial concentrations, the system will evolve to a final, or equilibrium, state in which (vd = vr).

The equilibrium constant of the first reaction, Keq1, comes directly from vd = vr:

$$\text{kp1}\mathrm{a}\left(\mathrm{eq}\right)\mathrm{e}\left(\mathrm{eq}\right)=\text{km1}\mathrm{c}\left(\mathrm{eq}\right),\text{then}\mathrm{kp}1/\text{km1}=\mathrm{c}\left(\mathrm{eq}\right)/\left(\mathrm{a}\left(\mathrm{eq}\right)\mathrm{e}\left(\mathrm{eq}\right)\right)=\text{Keq1}.$$

For the right reaction, kp2 c(eq) = km2 b(eq) e(eq), then Keq2 = kp2/km2 = b(eq)e(eq)/c(eq).

For the whole system:

$$\mathrm{Keq}=\mathrm{Keq}1\mathrm{Keq}2=\mathrm{kp}1\mathrm{kp}2/\left(\mathrm{km}1\mathrm{km}2\right)=\mathrm{b}\left(\mathrm{eq}\right)/\mathrm{a}\left(\mathrm{eq}\right)$$

The time rates of change of the concentrations {A, C, B, E} = {a, c, b, e} are expressed as sums of products of the concentrations according to the stoichiometry by the law of mass action:

$$\begin{array}{l}\mathrm{da}/\mathrm{dt}=‐\mathrm{kp}1\mathrm{a}\mathrm{e}+\mathrm{km}1\mathrm{c}\\ \mathrm{dc}/\mathrm{dt}=‐\left(\mathrm{km}1+\mathrm{kp}2\right)\mathrm{c}+\left(\mathrm{kp}1\mathrm{a}+\mathrm{km}2\mathrm{b}\right)\mathrm{e}\\ \mathrm{db}/\mathrm{dt}=+\mathrm{kp}2\mathrm{c}–\mathrm{km}2\mathrm{b}\mathrm{e}\\ \mathrm{de}/\mathrm{dt}=‐\mathrm{kp}1\mathrm{a}\mathrm{e}+\mathrm{km}1\mathrm{c}+\mathrm{kp}2\mathrm{c}–\mathrm{km}2\mathrm{b}\mathrm{e}=\mathrm{da}/\mathrm{dt}+\mathrm{db}/\mathrm{dt}\end{array}$$

It is assumed that the system is closed; that at all times, all the concentrations are homogenous; and that the sum of free plus substrate bound enzyme is constant, so e + c = Constant = e(0) + c(0) = Eo. Introducing (e = Eo – c) in the first three equations, we have:

$$\begin{array}{l}\mathrm{da}/\mathrm{dt}=‐\mathrm{kp}1\mathrm{Eo}\mathrm{a}+\left(\mathrm{km}1+\mathrm{kp}1\mathrm{a}\right)\mathrm{c}\\ \mathrm{dc}/\mathrm{dt}=‐\left(\mathrm{km}1+\mathrm{kp}2\right)\mathrm{c}+\left(\mathrm{kp}1\mathrm{a}+\mathrm{km}2\mathrm{b}\right)\left(\mathrm{Eo}–\mathrm{c}\right)\\ \mathrm{db}/\mathrm{dt}=‐\mathrm{km}2\mathrm{Eo}\mathrm{b}+\left(\mathrm{kp}2+\mathrm{km}2\mathrm{b}\right)\mathrm{c}\end{array}$$E3

Introducing the known parameters {km1, kp1, km2 kp2}, the system of three coupled ordinary, nonlinear, and differential equations (3) in the three unknown functions {a(t), c(t), b(t)} will be numerically solved starting from a given initial condition {Eo, ao, co, bo}, using the software Mathematica© by writing just two lines of code:

{kp1, km1, kp2, km2, Eo} = {1000, 1.0, 0.1, 10.0, 0.0001}; s = NDSolve[{{a’[t] == − kp1 Eo a[t] + (km1 + kp1 a[t]) c[t], a[0] == 0.001}, {c’[t] == (km1 + kp2) c[t] + (kp1 a[t] + km2 b[t])(Eo – c[t]), c[0] == 0}, {b’[t] == − km2 Eo b[t] + (kp2 + km2 b[t]) c[t], b[0] == 0}}, {a, c, b}, {t, 0, 500}];

The command NDSolve is a proprietary differential equation solution machine developed by the staff of mathematicians of Wolfram Scientific Corp, which after analyzing the problem to solve, automatically decides the methods and time steps to use, unless the user provides the necessary instructions; the solution is given as (Figure 1).

Figure 1.

Plots of the Mathematica© solution of the system of differential equations (3).

s = {{a - > InterpolatingFunction[], c - > InterpolatingFunction[], b - > InterpolatingFunction[]}, from which the following graphs are obtained using Plot[{Evaluate[a[t]/.s, Evaluate[c[t]/.s, Evaluate[b[t]/.s}, {t, 0, 500}}] (x axis units are seconds and y axis units are M (moles per liter)).

If we look carefully the graphs above, it seems that at t = 0, the a and c graphs have a sort of discontinuity. The plots of those concentrations for time {0 to 5} are shown below (Figures 2 and 3).

Figure 2.

Plot of the concentration of a for the first 5 s; y axis units are M.

Figure 3.

Plot of the concentration of c for the first 5 s; y axis units are M.

One can see that there are no discontinuities; however the substrate a reduces its initial concentration by only 5% in the first 2 seconds, while the substrate–enzyme complex c goes from its initial value of 0 at t = 0 to its maximum value of approximately 0.000045.

Then, what is wrong with this kind of Chemical Kinetics? According to Gillespie [3], the physical basis for this approach leaves something to be desired because it assumes that the time evolution of a chemically reacting system is continuous and deterministic. Time cannot be considered continuous, he says, as it is assumed by the use of differential equations, because the molecular populations change only when the reactions occur, and those changes are in the number of molecules of the species, which are whole numbers.

He claims that the chemical reactions cannot be deterministic because it is not possible to predict the exact molecular population levels at some future time unless we take account of the precise positions and velocities of all the molecules in the system. In some cases, as in microscopic biological systems, the inability of the reaction-rate equations to describe the fluctuations in molecular population levels can be a serious shortcoming, he concludes.

In the stochastic chemical kinetics, instead of trying to obtain the evolution in time of the concentrations of a number of reacting molecular species {c1, c2, c3, …}, one looks for the evolution of the probabilities of those concentrations; mathematically, the so-called Chemical Master Equation, CME, whose solution gives the probability that at time t, the values of the concentrations {c1, c2, ….} are exactly {x1, x2, …}, is expressed as:

$$\mathrm{p}\left(\left\{\mathrm{x}1,\mathrm{x}2,\dots ,\mathrm{t}\right\}\right)=Pr\left[\mathrm{c}1\left(\mathrm{t}\right)=\mathrm{x}1,\mathrm{c}2\left(\mathrm{t}\right)=\mathrm{x}2,\dots .\right\}]$$E4

Many times, particularly when the reactive systems have small numbers of molecules, the CME uses the numbers of molecules instead of concentrations, expressing the probability that at time t, the numbers of the reactive molecules {N1, N2, …} have certain given values {n1, n2, …}, which is written as:

$$\mathrm{p}\left(\left\{\mathrm{N}1,\mathrm{N}2,\dots ,\mathrm{t}\right\}\right)=Pr\left[\mathrm{N}1\left(\mathrm{t}\right)=\mathrm{n}1,\mathrm{N}2\left(\mathrm{t}\right)=\mathrm{n}2,\dots \right\}]$$

Given a system of coupled chemical reactions, Gillespie developed an algorithm by which it is possible to draw a diagram for the reaction jumps of the system starting from a given initial state. Using two random numbers between 0 and 1, by statistical arguments, one of them is used to determine the reaction time, τ, it takes to jump, and the other is used to determine, again using statistical arguments, the two states involved in the jump; by iteration of this procedure, the stochastic trajectories of the reactants will be obtained.

Using some indications and the parameters given in ([5] page 271), I will develop the Mathematica© code to model and simulate the stochastic kinetics of the previously considered Michaelis–Menten system:

$$\mathrm{E}+\mathrm{A}\left(\mathrm{\kappa m}1\right)<‐‐‐‐‐>\left(\mathrm{\kappa p}1\right)\mathrm{C}\left(\mathrm{\kappa m}2\right)<‐‐‐‐>\left(\mathrm{\kappa p}2\right)\mathrm{E}+\mathrm{B}$$

Be aware that the stochastic jump constants {κp1, κm1, κp2, κm2}, all of them having units inverse of time and are not equal to the previously used rates constants {kp1, km1, kp2, km2}.

Suppose that at time (t), we have m molecules of substrate A and n molecules of enzyme–substrate complex C, so the state of the system is given by [A(t), C(t)] = [m, n]. If at the beginning of the reactions, (t = 0), we had m = 100 molecules of A, 0 molecules of C, 0 molecules of B, and 10 molecules of free enzyme, Eo, then, at time (t), there will be (10 – n) molecules of free enzyme, E; remember that Eo = e + c = constant, and because Mo = a + b + c = constant = 100, there will be Mo – m – n = 100 – m – n molecules of B; then [A(t), C(t), B(t), E(t)] = [m, n, 100 - m – n, 10 – n].

At time (t), the stochastic state [m, n] of the reacting system is:

$$\mathrm{E}\left(10–\mathrm{n}\right)+\mathrm{A}\left(\mathrm{m}\right)<‐‐‐‐‐‐>\mathrm{C}\left(\mathrm{n}\right)<‐‐‐‐‐>\mathrm{E}\left(10–\mathrm{n}\right)+\mathrm{B}\left(100–\mathrm{m}–\mathrm{n}\right)$$

The set of 5 reaction states involved in the 8 jumps to and from the state [m, n] is shown below:

$$\begin{array}{l}\left[\mathrm{m}‐1,\mathrm{n}+1\right]\left[\mathrm{m},\mathrm{n}+1\right]\\ \left[\mathrm{m},\mathrm{n}\right]\\ \left[\mathrm{m},\mathrm{n}‐1\right]\left[\mathrm{m}+1,\mathrm{n}‐1\right]\end{array}$$

We observe that the values of m increase to the right—creation of a molecule of A-, and those of n increase to the top—creation of a molecule of C-.

I will first examine all the jumps that end at [m, n]:

The jump [m - 1, n + 1]➔[m, n]: A increases, C decreases, comes from the reaction C ---- > E + A, and has a jump strength (instead of reaction velocity km1 c) κm1 (n + 1), κm1 instead of km1 and the number of molecules of C, n + 1, instead of the concentration c.

The jump [m, n + 1]➔[m, n]: - A constant, C decreases, comes from C ----- > E + B, (instead of a reaction velocity km2 c), and has a jump strength κp2 (n + 1).

The jump [m, n-1]➔ [m, n]: -A constant, C increases-, comes from E + B ----- > C, and has a jump strength κm2 (10 – (n − 1))(100 – m – (n - 1)) (instead of km2 e b).

Finally, the jump [m + 1, n - 1]➔[m, n]: -A decreases, C increases-, comes from E + A ----- > C, and has a jump strength κp1 (10 – (n - 1))(m + 1) (instead of kp1 e a).

On the other hand, the possible jumps that start from the state [m, n] are:

[m, n]➔[m – 1, n + 1]: -A decreases, C increases- from E + A ----- > C has a jump strength κp1 (10 – n)m.

[m, n]➔ [m, n + 1]: -A constant, C increases- from E + B ➔ C has a jump strength κm2(10 – n)(100 – m – n),

[m, n]➔[m, n – 1]: -A constant, C decreases from C ----- > E + B has a jump strength κp2 n.

And [m, n]➔[m + 1, n - 1]; A increases, C decreases from C ----- > E + A has a jump strength κm1 n.

Using the previous jump strengths, Gillespie would arrive to the following equation for the Chemical Master Equation, CME, expressing the evolution of the probability of the general state [m, n] at time t, p(m, n, t):

dp(m, n, t)/dt = − (κp1 (10 – n)m + κm2(10 – n)(100 – m – n) + κp2 n + κm1 n)p(m, n, t) + κp1 (10 – (n - 1))(m + 1)p(m + 1, n-1, t) + κm1 (n + 1)p(m-1, n + 1, t) + κp2 (n + 1)p(m, n + 1, t) + κm2(10 –(1))(100 – m – (n − 1))p(m, n − 1, t).

The first term, with four contributions, is negative because it represents the jumps out of the state [m, n] that diminish the probability of finding the system at [m, n] at a later time, and the last four terms are positive because they increase that probability.

Now comes the simulation of the mathematical model given by dp(m, n, t)/dt. As mentioned before, each step will require two random numbers, {r1, r2}, obtained from a uniform distribution in the interval [0, 1].

Considering the jumping out of the state [m, n] as a first order decay process whose time distribution is given by f(t) = q e^{(−q t)} where q is equal to the negative term in the right side of the CME, q = κp1 (10 – n)m + κm2(10 – n)(100 – m – n) + κp2 n + κm1 n, and as 1/q represents the mean lifetime of the state, which obviously has units of time, the value of the random jump time, τ, would be given by: τ = − log (r1)/q.

The second random number, r2, is used to determine to which one of the four possible states the system would jump. The probability of each state is given by its corresponding jump strength divided by the sum of the jump strengths of the three states; for example:

p1 = κp1(10 – n)m/q = κp1(10 – n)m/(κp1(10 – n)m + κm2(10 – n)(100 – m – n) + κm1 n + κp2 n).

p2 = κm2(10 – n)(100 – m – n)/q, p3 = κm1 n/q, p4 = κp2 n/q. Evidently, p1 + p2 + p3 + p4 = 1.

The Mathematica© program to simulate Gillespie’s Algorithm for this system is the following—considering the case of an irreversible system; κm2 = 0.

{κp1, κm1, κp2, κm2} = {1.0, 0.5, 2.0, 0};

{mo, no, Z, t} = {100, 10, 1000, 0};

{a, c, b, e} = {Table[{0, 0}, {i, Z}],Table[{0, 0},{i, Z}],Table[{0, 0},{i,Z}], Table[{0, 0},{i, Z}]};

{m, n} = {mo, 0};

{a[[1]], c[[1]], b[[1]], e[[1]]} = {{t, m}, {t, n}, {t, mo - m - n}, {t, no}};

j = 2;

Do[

q = κp1 (10 - n) m + κp2 n + κm1 n;

{r1 = RandomReal[], r2 = RandomReal[]};

t = t - Log[r1]/q;

If[r2 < p1, m = m – 1; n = n + 1, If[r2 < p1 + p2, m = m + 1; n = n - 1; n = n − 1]];

{a[[j]], c[[j]], b[[j]], e[[j]]} = {{t, m}, {t, n}, {t, mo – m - n}, {t, no -n}};

j = j + 1,

{i, Z-1}].

Using the commands ListLinePlot[{a, b, c} and {ListLinePlot{Take[a, 20], Take[c, 20], Take[c, 20]}, we get the graphs below.

Figure 4 graphs are the stochastic {a, c, b} Michaelis–Menten set but now for a 6 seconds stochastic simulation.

Figure 4.

Plots from the Mathematica© code of the Gillespie’s Algorithm, x seconds, y # of molecules.

The difference in the time values from the deterministic set comes from the values of the jump constants, which are different from the rate constants, as well as the values of the numbers of each of the reactants, which are different from their corresponding concentrations. Figure 5 shows the main changes in the first 0.45 seconds.

Figure 5.

Plots from the Mathematica© code of the Gillespie’s Algorithm, first 0.45 s.

From the previous analysis, it is clear that neither the smoothness of the concentrations nor their chaotic fluctuations are real phenomenological predictions; instead, they are the result of the mathematics used in the models.

3. A new kind of chemical kinetic

This new kind of Chemical Kinetics will bring to the understanding of this phenomenon some new mathematical concepts and structures, hopping to capture predictions with more of its real phenomenological essence.

To begin, I will not consider that the timeline in which all the reactions occur is a continuum in which there may exist intervals from an infinitesimal to an infinite number or a sequence of intervals each one determined in a probabilistic ad hoc manner. Rather, I will consider that the timeline is a sequence of time points where each consecutive pair is separated by an adjustable value τ.

Within this timeline, I will pay attention to a pair of consecutive points numbered (n − 1), which I will rename m, because of the minus and n; {m,n}. Returning to the Michaelis–Menten system of differential equations (3), I will integrate the first one of them from time tm to time tn, in Mathematica© language:

Integrate [da/dt, {t, tm, tn}] = Integrate[− km1 Eo a + km1 c + kp1 a c, {t, tm, tn}]

The left-hand term is equal to Integrate[da, {t, tm,tn}] = a(tn) – a(tm) = an – am. In order to integrate the terms of the right hand side, I will make the concentration of any reactant to change linearly in time, from the beginning to the end of any time interval; thus, a(t), from tm to tn, should obey the equation: a(t) = a(tm) + t(a(tn) – a(tm))/(tn – tm) = am +t (an - am)/τ, so:

Integrate[am + t (an–am)/τ, {t, tm, tn}] = am τ + (an–am)/τ Integrate[t, {t, tm, tn}] = τ(an – am)/2, which is exactly the area under the line from am to an and the horizontal time axis.

Similarly, the integration of the second term will give + km1 τ(cn – cm)/2. For the last term kp1 Integrate[a(t) c(t), {t, tm, tn}] = kp1 Integrate[(am+t(an-am)/τ)(cm + t(cn - cm)/τ), {t,tm,tn}] gives: kp1 Integrate[am cm + t (am(cn – cm) + cm(an – am))/τ + t²(an – am)(cn – cm)/τ².

Doing the three integrations, applying the limits, and simplifying and arranging all the remaining terms, we get: kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6.

Using the previous results, the linear integration of the system of the three differential equations (3) gives:

an – am = − kp1 Eo (am + an) τ/2 + km1 (cm + cn) τ/2 + kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6

cn – cm = − (km1 + kp2) (cm + cn) τ/2 + Eo kp1 (am +an) τ/2 + Eo km2 (bm + bn) τ/2–kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6 – km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6

bn – bm = − km2 Eo (bm + bn) τ/2 + kp2 (cm + cn) τ/2 + km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6

Suppose that the initial conditions for the reactants are {A, C, B} = {ao, co, bo} = {100, 0, 0} and suppose that eo = 10, so Eo = eo + co = 10. If we introduce those values in the previous three equations, along with the other known values, we can solve the algebraic problem of three equations with the three unknowns {a1, c1, b1} and in the same way, from {a1, c1, b1}, determine {a2, c2, b2} and so on and so forth.

However, instead of the previous mathematical solution procedure, typical of the scientist who aims for 100% accuracy, I will make use of a mathematical optimization procedure, typical of the engineer who aims for the best possible accuracy to handle a given problem.

In other words, I will work under the hypothesis that the kinetics of any system of coupled chemical reactions, with many chemical components, can be visualized not as a set of trajectories of points in the concentration space, but rather as a collection of micro-tubes, which allows a coordinated minimal transversal motion apart from the usual coordinated forward motion coming from the kinetic equations.

For the optimization procedure, we need to construct the Objective Function, OF, which contains a collection of variables, together with many other constants and known parameters, and instead of finding the values of all the variables that render OF = 0, we try to find the values of the variables that render OF equal to the minimum possible value that it can admit. That value can be very close to 0, but that small difference has an important physical meaning; it is connected with the diameter of the previously mentioned micro-tubes.

I will construct the Objective Function with the following mathematical structure:

$$\text{OF}=\mathrm{f}{1}^2+\mathrm{f}{2}^2\dots ..+{\mathrm{fn}}^2,$$E5

where each of the functions fj—also called constraints—is supposed to be equal to 0, and thus, OF should be equal to 0. In our case at hand, the fj comes from the integrated differential equations, so:

OF = (− an + am - kp1 Eo (am + an) τ/2 + km1 (cm + cn) τ/2 + kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6)²+

(− cn + cm - (km1 + kp2) (cm + cn) τ/2 + Eo kp1 (am +an) τ/2 + Eo km2 (bm + bn) τ/2 – kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6 – km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6)² +

(− bn + bm -km2 Eo (bm + bn) τ/2 + kp2 (cm + cn) τ/2 + km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6)²

The Mathematica© program to obtain the optimum micro-tube trajectories of the concentrations {a, c, b} from the initial condition {100, 0, 0} in our Michaelis–Menten system is the following:

{kp1, km1, kp2, km2, Eo} = {1000, 1.0, 0.1, 10.0, 0.0001};

{am, cm, bm} = {0.001, 0, 0};

{τ, Z} = {1, 500};

{a, c, b, d} = {Table[0, {i,Z}], Table[0, {i,Z}], Table[0, {i,Z}], Table[0, {i,Z}]};

{a[[1]], c[[1]], b[[1]], d[[1]]} = {am, cm, bm, 0};

Do[

OF = (− an + am - kp1 Eo (am + an) τ/2 + km1 (cm + cn) τ/2 + kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6)² +.

(− cn + cm - (km1 + kp2) (cm + cn) τ/2 + Eo kp1 (am +an) τ/2 + Eo km2 (bm + bn) τ/2 – kp1 (2 am cm + 2 an cn + an cm + am cn) τ/6 – km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6)² +.

(− bn + bm - km2 Eo (bm + bn) τ/2 + kp2 (cm + cn) τ/2 + km2 (2 bm cm + 2 bn cn + bn cm + bm cn) τ/6)²;

Opt = FindMinimum[OF, {an, am}, {cn, cm}, {bn, bm}}];

d[[j]] = Opt[[1]];

{a[[j]], c[[j]], b[[j]]} = {an, cn, bn}/.Opt[[2]];

{am, cm, bm} = {a[[j]], c[[j]], b[[j]]};

j = j + 1,

{i, 2, Z}].

After running this program, we get four Tables, or Lists, with 500 values each; three of them, {a, c, b}, for the trajectories of the optimal collective concentrations as points separated by 1 sec., the fourth List is for d, the overall kinetic offset, during the reaction time, which is the sum of the squared kinetic offsets of the three constrain functions, fk, in the Objective Function (Figures 6 and 7).

Figure 6.

The set of Michealis–Menten collective optimal concentrations {a, b, c} from the new kind of chemical kinetics, which later on will be compared with the results obtained from the differential equations (3).

Figure 7.

The plot of the overall kinetic offset d, which is the collection of minimum values of the objective function for all the iterations of the optimization.

The graph above, for the overall kinetic offset d, shows that at the beginning, when the systems is the most far from equilibrium, after a spike, the values first rise abruptly and then decrease asymptotically to zero.

One first discovery is that the value of d, at any given time while the reactions are taking place, reflects how far from equilibrium is the system. The graphs of the optimal values of individual constraints, {f1, f2, f3} top to bottom, are shown below:

Remember that d = f1² + f2² + f3².

A second discovery is that the three f’s have the same general shape (Figures 8–10), although one with a positive offset and the other two with a negative offset, with values different several orders of magnitude. However, the most important kinetic offsets, or dynamical concentration offsets, are the ones occurring for each of the reagents. The graph below shows the difference between the values of the old deterministic kind of chemical kinetics and the new kind that I will call quasi-deterministic (Figure 11).

Figure 8.

Plot of the optimum f1 constraint values, much larger than the values for d above.

Figure 9.

Plot of the optimum f2 constraint values, also much larger than the values for d above.

Figure 10.

Plot of the optimum f3 constraint values, smaller than both f1 and f2.

Figure 11.

Plot of the difference between the old kind and the new kind of chemical kinetics; the positive peak is for a, and the negative peak is for b; the very small peak near the x axis is for c.

In order to determine which results are more accurate, from the two models, the old and the new Michaelis – Menten sets of plots, the error in the conservation of mass, a + b + c – 0.001, for t = {0, 5000} was calculated, and it is shown below.

In the old model, we see that as time increases the error in the conservation of mass grows quasi-monotonically in the negative half plane, while in the new formulation, that error is smaller and after one cycle, finishing at about t = 1000 sec.; the error settles at about - 4. 10⁽⁻¹⁹⁾.

One extremely valuable operation that we can do with the new formulation, and cannot be done with the old one, is to add an f4 term to the Objective Function of the form (an + cn + bn – 0.001)²; obviously this additional constraint represents the conservation of mass—remember that the previous three constraints {f1, f2, f3} come from the three reactions.

The new concentration sets of graphs look very similar including the third one for the difference with the old kind, but the conservation of mass versus time graph shown below is absolutely different (Figures 12–14).

Figure 12.

Plot of the error in the conservation of mass in the old deterministic chemical kinetics.

Figure 13.

Plot of the error in the conservation of mass in the new quasi- deterministic model.

Figure 14.

Plot of the error in the conservation of mass in the new quasi-deterministic model including a fourth constraint, f4, expressing the conservation of mass.

From the last graph, we discover that all the spikes have the same height of 2 × 10⁻¹⁹ (as those in the previous Figure); if at a certain time, the result of the mass conservation is not zero, that situation is corrected in the next iteration. This means that the new model, with the mass conservation constraint, imposes to the collection of optimum concentrations of the reactive system some kind of self-organization.

Karsenti [6], in his brief history about self-organization in cell biology, concludes that the dynamic order of the cell results from a combination of complex stereospecific interactions which result in self-assembly and other extremely varied dynamical interactions between molecules that result in self-organization. Perhaps some of those dynamical interactions occur in a real process to select the sets of collective optimum concentrations, even in small closed reactive systems.

If the initial concentrations {ao, eo} are increased 100 fold, that is, {ao, eo} = {0.1, 0.01}, and all other parameters are kept equal, we will see fluctuations and self-organization occurring together. From top to bottom, the first graph comes from the new model, the second one comes from the old model, and the third one is the overlapping of the old and the new concentration of c at the beginning of the reaction time (Figure 15).

Figure 15.

If the initial concentrations of a and e increase 100 times, the set of the old chemical kinetics does not change, while fluctuations appear in the a and the c plots of the new model. The amazing difference between the c plots in the old and the new formulation at the beginning of reactions is shown above.

The false smoothness imposed by the old model contrasts with the real decaying fluctuations of the enzyme–substrate complex, c, emerging from the new model from the synchronization of its chemical reactions with a as well as from the self-organizing behavior from the mass conservation (Figure 16).

Figure 16.

Mass conservation errors of the old and new models, the last showing self-organization.

4. Conclusions

A New Kind of Chemical Kinetics, and beyond, in 12 steps.

1. Divide the time line in intervals of a conveniently chosen value τ.

2. Consider the interval from (n - 1) = m and n: {m- > n}.

3. Linearize each one of the concentrations, for example, a, in the reactive system inside the interval {m- > n}: a = am + t (an – am)/τ.

4. Introduce the linearizations and integrate in the interval {tm - > tm} each of the kinetic differential equations from the law of mass action.

5. In the previous equations, the left hand side, for example, da/dt, will become (an – am) and similarly for the rest of the reactants.

6. If the reactive system has N reactants, after the N integrations, we will get N algebraic equations, with N—unknown—concentrations {an, bn, cn,…}, while the previous concentrations {am, bm, cm,…) are known.

7. Instead of solving the system of N algebraic equations, we will find the set of concentrations that collectively minimizes the value of the Objective Function, OF, constructed with the N algebraic equations as constraints, fi, raised to power of 2, so OF = fa² + fb² + ….. fN².

8. Starting with the concentrations at t = 0, in a chain of optimizations, we can get the optimal values of the concentrations after each interval.

9. Apart from that, we can also get the values of OF and {fa, fb, fc,..} at the collective optimizations, whose analysis could provide new information.

10. It is possible to add a term to OF expressing the conservation of mass of the system; doing this, it is possible to interpret some changes in the resulting graphs as the appearance of self-organization.

11. For concentrated systems, the appearance of synchronization is a fact.

12. I consider that in the future, the main contribution of this new kind of Chemical Kinetics could be its extension to the case of living systems.