Formal Kinetics with Applications

János Tóth

The usual deterministic model of a complex chemical reaction is a polynomial differential equation with as many variables as the number of species. Mathematical tools called reaction network theory are known for the qualitative investigation of such systems. However, application of these results needs lengthy routine symbolic calculations. A Mathematica package is under development and testing to overcome these difficulties starting from the input, through the construcion of different graphs corresponding to the reaction, up to qualitative, and (less automatic) quantitative and graphical study of complex models. The usual stochastic model of the given reaction is also simulated. The speciality of the present programs is that they are really general, and work independently form the number of chemical species and reaction steps as opposed to the Mathematica programs found in the literature (including the web).

Canonical form of a reaction

The reaction we shall investigate as an example is the Michaelis-Menten model

E + S⟷C⟶E + 
      P

where E is the enzyme, S is the substrate (the material to be transformed by the enzyme), C is an intermediate complex formed by the enzyme and its substrate, and P is the product of the reaction. The model above can also be interpreted as a description of receptor-ligand interaction.

In order to make further investigations we dissect the above reaction into the elementary steps

E + S -> C,  C⟶E + S,  
    C⟶E + P

among the species E,  S,  C,  P. This is done in the following obvious way.

\(ToCanonical[chains_] := 
      Module[{rulebreak = {

              x_ ↔ y_ ↔ 
                  z__ -> {x ↔ y, 
                  y ↔ z}, \((x_ ↔ y_ -> 
                    z__)\) -> {x ↔ y, y -> z}, 

              x_ ↔ y_ ← 
                  z__ -> {x ↔ y, 
                  y ← z}, \((x_ -> 
                    y_ ↔ z__)\) -> {x -> y, 
                  y ↔ z}, 

              x_ ← y_ ↔ 
                  z__ -> {x ← y, 
                  y ↔ 
                    z}, \((x_ -> \(y_ -> 
                      z__\))\) -> {x -> y, 
                  y -> z}, 
\((x_ -> 
                    y_ ← z__)\) -> {x -> y, 
                  y ← z}, \((x_ ← y_ -> 
                    z__)\) -> {x ← y, 
                  y -> z}, 
\((x_ ← y_ ← 
                    z__)\) -> {x ← y, 
                  y ← z}}, 		

          rulerev = {\((x_ ← y_)\) -> \((y -> 
                    x)\), \((x_ ↔ y_)\) -> {x -> y, 
                  y -> x}}}, 

        Flatten[\(chains //. rulebreak\) /. rulerev]]\)

First, longer chains are broken up, then reversible reaction steps are dissected into pairs of irreversible steps, or inappropriately written arrows are reversed. Let us see a small part of the program in working.

intermediateres = Ε + S ↔ 
            C -> Ε + 
            P /. (x_ ↔ y_ -> 
              z__) -> {x ↔ y, y -> z}

intermediateres[[1]] /. (x_ ↔ y_) -> {x -> y, y -> x}

Now it is easy to understand the total effect of the rules. rulebreak must be applied repeatedly. The unnecessary brackets are eliminated by Flatten. Note that ToCanonical 
works in the same way on a single (however complicated) chain, or on a single element list with the same chain.

canMM = 
      ToCanonical[Ε + S ↔ 
            C -> Ε + P]

Another example might be the well known Volterra-Lotka model

A + X -> 2  X,  X + Y -> 2  Y,  
      Y -> B,

where A and B are external species.

canVL = 
      ToCanonical[{A + X -> 2  X, X + Y -> 2  Y, 
          Y -> B}]

Stoichiometric characteristics

Now we name and calculate the number of different chemical species, and of reaction steps, collect the elements of the matrix of stoichiometric coefficients (numbers before the signs of species names, mainly 0 or 1), and we create names to describe the quantities of the chemical species. Before showing the program let us define a few quantities in theory. Let us consider the reaction

(∑_(m = 1)^M) α(m, r) X(
            m) -> (∑_(m = 1)^M) β(m, r) X(
            m)    (r = 1, 2, ⋯, 
            R),

where M is the number of chemical species X(m), R is the number of reaction steps, the nonnegative integers α(m, r) and β(m, r) are the stoichiometric coefficients. In some cases we may also have external species of which the concentration does not change; if this is the case then the species we are interested in are the internal species. The formal linear combinations on both sides of the reaction steps are complexes; reactant and product complexes respectively. All the complexes can be described by vectors in (Subscript[ℕ, 0]^M), \(TraditionalForm\`\(α(\(.\)\(\(,\)\(r\)\))\)\ and\ \(β(\(.\)\(\(,\)\(r\)\))\)\) are the reactant and product (complex) vectors. These quantites are determined from the canonical form of a reaction by the program below in the form as it is used in the package Statistics. Finally, we introduce concentration variables constructed from the names of the species. (Another alternative would be to denote the concentration of the cpecies X by [X] as ti is ususal in chemical kinetics, or in physical chemistry.)

\(ToStoichiometry[r_, externalspecies_List: {}] := 

      Module[{RemoveExternal, 
\ \ ReactionSteps, 
          Species, 
	Complexes, ComplexVectors, 
	Reactants, 
          ReactantVectors, 
	Products, ProductVectors, 
	
          RemoveSigns = {\((x_ -> y_)\) -> {x, 
                  y}, 
			\ \ \ \ \ n_Integer\ x_ -> x, 
				\ \ \ \ \ \ \ x_ + y_ -> {x, 
                  y}}, 
ff}, 

        RemoveExternal = {A_ /; MemberQ[externalspecies, A] -> 0}; 

        ReactionSteps = Sort[r //. RemoveExternal]; 

        Complexes = Union[Cases[ReactionSteps, _, {2}]]; 

        Species = 
          Sort[Complement[
              Flatten[r //. RemoveSigns], 
				externalspecies]]; 

        Reactants = First /@ ReactionSteps; 

        Products = 
          Last /@ ReactionSteps; 
{ComplexVectors, ReactantVectors, 
            ProductVectors}
	 = \((Outer[Coefficient[#1, #2] &, #1, 
                  Species, 1] &)\) /@ 
{Complexes, 
              Reactants, Products}; 
{alpha, beta} = 
          Transpose /@ {ReactantVectors, ProductVectors}; 

        ff[x_Symbol] := ToExpression[ToLowerCase[ToString[x]]]; 

        ff[x_String] := ToLowerCase[x]; 	

        ff[x_] := MapAt[ff, x, {{1}}]; 
	{species -> Species, 
          M -> 
            Length[Species], 
\ \ \ complexes -> 
            Complexes, Ν -> Length[Complexes], 
	
          reactionSteps -> ReactionSteps, 
          R -> Length[ReactionSteps], 
	α -> 
            alpha, β -> beta, 
	
          variables -> ff /@ Species}]\)

More than one examples can be seen for the application of the use of list processing functions and term rewriting.

See, how external species are removed and the reaction steps (as opposed to the original genuine reaction steps) are obtained.

extVL = canVL //. {a_ /; MemberQ[{A, B}, a] -> 0}

Complexes are expressions of any form at level 2.

Level[extVL, {2}]

The result should contain each complexes only once.

cxVL = Union[Cases[extVL, _, {2}]]

cxMM = Union[Cases[canMM, _, {2}]]

In order to collect species, first we discard  numbers and arrows, then we get rid off external species,

nosignsVL = 
      canVL //. {(x_ -> y_) -> {x, y}, 
          n_Integer x_ -> x,  x_ + y_ -> {x, y}}

spVL = Complement[Flatten[nosignsVL], {A, B}]

It seems to be more complicated to first discard external species and then get rid off numbers and arrows.

Obviously, reactants and products are those complexes which stand on the sides of the reaction arrows.

reVL = First /@ extVL

prVL = Last /@ extVL

Stoichiometric coefficients are obtained if one determines the coefficients of the species in the formal linear combinations representing complexes.

Map[Coefficient[#, spVL] &, cxVL]

Map[Coefficient[#, spVL] &, reVL]

Map[Coefficient[#, spVL] &, prVL]

The need to use Outer arises from the fact that we would like to calculate the coefficients in all the complexes, in the reactant complexes and in the product complexes with a single command.

Outer[Coefficient[#1, #2] &, #1, spVL, 1] & /@ {cxVL, reVL, 
        prVL}

The names of the concentrations of the individual species is the corresponding lowercase letter.

ff[x_Symbol] := ToExpression[ToLowerCase[ToString[x]]];

ff[S]

However,this is not so trivial, because species names can be very complicated in chemistry.

\({ff[\(H\_2\) O], ff[\*""\!\(H\_2\)O""]}\)

\(\(\(ff[x_String] := ToLowerCase[x];\)\(	\)\)\)

\({ff[\(H\_2\) O], ff[\*""\!\(H\_2\)O""]}\)

Had we chosen the alternative mentioned at the end for the notation of the concentrations, we would have used the following function instead of ff in the Module.

SetAttributes[h, Listable];
h[x_] := 
        StringJoin["[", 
          ToString[x], "]"];
h[{E, S, C, P}]

Note that the right hand side

ToString["[x]"]

will never result in, say, [E].

stMM = ToStoichiometry[canMM]

stVL = ToStoichiometry[canVL, {A, B}]

The right-hand-side of the induced kinetic differential equation

The right-hand-side of the usual deterministic model of a chemical reaction is constructed in the following way: reaction steps have a certain rate depending on the quantities of the individual species, species change of their quantity depending on these rates of the steps and on the differences of the stoichiometric coefficients. The rate of the individual steps is most often also calculated from the stoichiometric coefficients, this is the case of mass action kinetics, or it can be given (almost) arbitrarily, following only a few qualitative restrictions.

The most common way is to use mass action type kinetics which means that the reaction rate of the reaction step

E + S -> C

being (Subscript[k, 1]) e s, is proportional with the concentrations of the reactants, and the proportionality coefficent is a positive real number. In case of the reaction (4) one has the following reaction rate

\(TraditionalForm\`\(k(
        r)\) \(∏\+\(p = 1\)\%M\( x(p)\)\^\(α(p, r)\)\ \ \ \(\((r = 1, 2, ⋯, 
              R)\)\(.\)\)\)\)

Here are the reaction steps and variables of the Michaelis Menten reaction.

{reactionSteps, variables} /. stMM

List operations and listability of some Mathematica functions make it easy to calculate the reaction rates even in the most general case.

{(Subscript[k, 2]), (Subscript[k, -1]), 
        (Subscript[k, 1])} Times @@ {c, p, s, ϵ}^(α /. 
              stMM)

As it is most common, we consider the mass action case as default, and this will be given as the default value of the option MassAction. If no option value is specified or if it required to be True, then the reaction rates will be calculated by the mass action formula (6) above.

We have used Input, and before asking for the reaction rates we present the reaction steps to see their ordering. Using files for reactions with many steps and species is much more useful. Now let us collect everything into a single Mathematica fucntion.

Options[RHS] = {MassAction -> True};
\(RHS[stres_, opts___?OptionQ] := 
      Module[{w, rate}, 

        Print[reactionSteps /. stres, variables /. stres]; 

        rate = Input["Provide the list of rates in the order of reaction steps."]; 

        w = If[\(MassAction /. {opts}\) /. Options[RHS], 

            rate\ Times @@ \((variables\^α /. stres)\), rate]; 

        If[\(MassAction /. {opts}\) /. 
            Options[RHS], 
\((\((β - α)\) /. stres)\) . 
            w, \((\((β - α)\) /. stres)\) . w]]\)

It may be interesting to consider the above definition as an example of defining an option to our own function, and also that of using external input and repeatedly using term rewriting again.

rhs = RHS[st]

Numerical solution of the induced kinetic differential equation of a reaction

Once we have the right-hand-side, we might wish to solve the induced kinetic differential equation of a reaction. This is almost always only numerically possible.

Timer[lis_List, timvar_:  t] := Map[#[timvar] &, lis]

der = D[Timer[variables /. st], t];

rhst = 
        rhs /. Thread[
            variables /. st -> Timer[variables /. st]];

sol=NDSolve[Join[{c'[t]==-k2 c[t] - km1 c[t] + 
                          k1 s[t] ϵ[t],p'[t]==k2 c[t],s'[t]==km1 c[t] - k1 s[t] ϵ[t],ϵ'[t]==k2 c[t] + km1 c[t] - 
                          k1 s[t] ϵ[t]},{c[0] == 0, p[0] == 0, 
                    s[0] == 0.001, ϵ[0] == 0.01}]/.{k1 -> 10, km1 -> 1, 
                k2 -> 1},{c[t], p[t], s[t], ϵ[t]},{t, 0, 5}];

Plot[Evaluate[{s[t], c[t]} /. sol], {t, 0, 5}, 
        PlotStyle -> {{Thickness[0.015], Red}, {Thickness[0.015], 
              Green}}];

Plots are to be drawn by hand in all the cases because one cannot foresee what is going to happen.

Deficiency

The deficiency of a reaction is defined to be the difference of the number of complexes N, and the sum S of the dimension of the linear space generated by the vectors of stoichiometric coefficients and the number L of connected components of the Feinberg-Horn-Jackson graph of the reaction of which the complexes are the vertices, and the reaction arrows are the edges. This number can easily be computed from the quantities provided by ToStoichiometry.

NN = {Ν /. stMM, Ν /. stVL}

We need a function to calculate the rank of a matrix which can be done like this.

MatrixRang[a_] := 
      Length[Transpose[a]] - Length[NullSpace[a]]

SS = 
      MatrixRang /@ {β - α /. stMM, β - α /. 
            stVL}

We need DiscreteMath to deal with graphs.

<< DiscreteMath`

We form  ordered pairs from the reactant and product complexes of the reaction steps.

opMM = Flatten[Apply[List, reactionSteps /. stMM, 1]]

We replace the complexes by their position.

orMM = 
      Partition[Flatten[Position[cxMM, #] & /@ opMM], 2]

ShowGraph[FromOrderedPairs[orMM]];

We do the same preparations for the Volterra Lotka model, then calculate the number of connected components.

opVL = Flatten[Apply[List, reactionSteps /. stVL, 1]]

orVL = 
      Partition[Flatten[Position[cxVL, #] & /@ opVL], 2]

ShowLabeledGraph[FromOrderedPairs[orVL]];

LL = Length[ConnectedComponents[FromOrderedPairs[#]]] & /@ {orMM, 
          orVL}

δ = NN - LL - SS

We learned that the deficiency of the Michaelis-Menten model is zero, thus no exotic behaviour can be expected from this, but the deficiency of the Volterra-Lotka model is one implying that some kind of exotic behaviour is not excluded.

Characteristics of the Vol'pert graph

The Vol'pert graph (the graph often used in textbooks on biochemistry) is obtained if one forms the set of species and elementary steps and draws an arrow from a species to a reaction step, if the species takes part in the reaction step, or from a reaction step to a species, if the species is produced in the reaction step. The stoichiometric coefficient is also assigned to each arrow. Properties of this graph reflecxt important characteristics of the reaction, see Vol'pert and Hudyaev (1975) and Érdi and Tóth (1989). We construct here the indices of the reaction steps and species, respectively. Initially present species have zero index, reaction steps starting from initial species also have zero index. Larger indices are defined recursively.

\(Volpert[α_, β_, MM_, RR_, MM0_, tiltottak_: {}] := 

      Module[{M = Length[MM], R = Length[RR], 
          M0 = Flatten[\(Position[MM, #] &\) /@ MM0],  (*\(reachable\)\(,\)*) 
            speciesindex, reactionindex, R0, Mold, Rold, Mnew, Rnew, 
          Mtoreach, Rtoreach, a, b}, 
Print["Initial species: ", M0]; 

        reachable[r_, m_List] := 
          And @@ \((\(\((\((α[[#, 
                            r]] > 0)\) \[Implies] 
                      MemberQ[m, #])\) &\) /@ Range[M])\); 

        reachable[r_List, 
            m_] := \ \ Or @@ \((\(\((β[[
                        m, #]] > 0)\) &\) /@ r)\); 	

        speciesindex = Array[\(+∞\) &, M]; 
        reactionindex = Array[\(+∞\) &, R]; 

        speciesindex[[M0]] = 0; \ 

        R0 = Select[Range[R], 
            reachable[#, 
                M0] &]; 
reactionindex[[
            R0]] = 0; 
Mold = M0; Rold = R0; 
        Mtoreach = M0; 
        Rtoreach = R0; 
		While[
          True, 
	\((Mnew = 
              Select[Range[
                  M], \((speciesindex[[#]] == \(+∞\))\) ∧ \((reachable[
                        Rtoreach, #])\) &]; 
		speciesindex[[Mnew]] = 
              Max[speciesindex[[
                    Mtoreach]]] + 1; 
		Rnew = 
              Select[Range[
                  R], \((reactionindex[[#]] == \(+∞\))\) ∧ \((reachable[#, 
                        Mtoreach])\) &]; 
		reactionindex[[Rnew]] = 
              Max[reactionindex[[
                    Rtoreach]]] + 1; 
		Mold = Mnew; 
            Rold = Rnew; Mtoreach = Join[Mtoreach, Mnew]; 
            Rtoreach = Join[Rtoreach, Rnew])\); 
			If[
            Mnew == {} ∧ Rnew == {}, Break[]]]; 

        a = TableForm[
            Sort[Transpose[{MM, 
                  speciesindex}], #2[[2]] > #1[[2]] &], 
            TableHeadings -> {{}, {"species", index}}]; 

        b = 
          TableForm[
            Sort[Transpose[{RR, 
                  reactionindex}], #2[[2]] > #1[[2]] &], 
            TableHeadings -> {{}, {"reaction step", index}}]; 

        Print[a]; Print[b];]\)

Volpert[α /. stMM, β /. stMM, species /. stMM, 
      reactionSteps /. stMM, {Ε, S}]

You may not see that enzyme is denoted by Esc+E+Esc, in order to make distinction between enzyme and the base of natural logarithm.

Volpert[α /. stVL, β /. stVL, species /. stVL, 
      reactionSteps /. stVL, {X}]

Before turning to reactions let us consider the usual Poisson process. Its trajectories behave in the following way: the process starts for state zero, stays in that state for an exponentially distributed random time, then it goes to state one. Again an exponential waiting time follows, and again a jump, now to state two, and so on. One realization of such a process vcan be seen in the figure below produced by the anticipated use of the final stochastic simulation program.

Kurtz[{{0}}, {{1}}, 1000, {0.001}, 1, 100, 1000]

Let us make a few remarks on the Poisson process. First, it is the stochastic analogue of the simplest possible differential equation x' (t) = λ.Second, it is the stochastic model of the reaction A -> X, where A is an external species of constant concentration.

The usual stochastic model of a chemical reaction is a Markovian pure jump process of a similar, but more complicated structure. The state of the process is a vector of nonnegative integers: the numbers of molecules of the individual species. The process stays at a state for an exponentially distributed random time, then one of the reaction steps is selected. We have to define the parameter of the exponential distribution and the distribution according to which the reaction step is selected. Let us start with the second one: each step is selected with a probability proportional to the actual reaction rate. Now the parameter of the waiting time between two jumps caused by a reaction step is just the sum of the probabilities of the reaction rates. There is one minor problem to solve: how to calculate the stochastic reaction rate constants once the deterministic reaction rate constants are known. This is solved at the beginning of the program.

<< Graphics`

Suppose that the initial concentration of the species in the deterministic model is given by the vector (x^0), and it is given that the concentration of X(m) is to be repesented by nmolecules in the stochastic model. Then, obviously, the initial state of the stochastic model (the initial number of molecules) are to be calculated by the following formula:  \(TraditionalForm\`\(x\^0\/x\_m\%o\) n\). At the beginning this will be the state of the stochastic model.

The rate constants are also to be transformed, and their transformation depends on the order of the reaction too, and is to be carried out in the following way.

\(TraditionalForm\`k\^stoch = \(\(k\^det\)(n\/x\_m\%o)\)\^\(∑\+\(m = 1\)\%M α(m, r) - 1\)\)

The simulation will go on until the sum of the waiting times is smaller than the prescribed limit. The reaction rates are calculated a bit differently from the deterministic case. If we have n molecules from the species X, than the rate of the reaction 2  X -> ⋯  will be k n(n - 1), that of the reaction 3  x -> ⋯ will be k n(n - 1) (n - 2).The built-in function Pocchammer  can be used to calculate these numbers.

{Pochhammer[n, 3], Pochhammer[n - 2, 3]}

Exponentially distributed random variables can be obtained from the package Statistics.

<< Statistics`

Values of a discrete random variable can be obtained by creating the partial sums (practically the cumulative distribution function) from the probabilities. This is needed to select from the reaction steps and is carried out using FoldList. The present version of the program has an automatic plotting. A more useful version should provide the simulation results for further elaboration (statistical evaluation) by the user.

\(Kurtz[alpha_, beta_, detrateconstants_, detini_, stochiniindex_, 
        stochini_, timelimit_, ε_:  10\^\(-6\), opts___] := 
      Module[{waitingtimes = {}, 
          jumps = {}, 
reactionvectors = 
            Transpose[beta - alpha], 
initialstate = 
            detini\ stochini\/detini[[stochiniindex]], state, 
          jump, 
w = Array[1 &, Length[detini]], wt, jp, 
          rf, 
selectedstep, 
          stochrateconstants}, 
state = 
          initialstate; 
stochrateconstants = \(\((stochini\/detini[[stochiniindex]])\)\^\(\((Plus @@ 
                      alpha)\) - 1\)\) 
            detrateconstants; 
While[
          Plus @@ waitingtimes ≤ timelimit, 
          w = Apply[Times, Transpose[Pochhammer[state - alpha + 1, alpha]], 
              1]; 
If[Plus @@ w < ε, 
            Break[]]; 
wt = 
            Random[ExponentialDistribution[λ = 
                  stochrateconstants . w]]; 
AppendTo[
            waitingtimes, 
            wt]; 
jp = \(stochrateconstants\ w\)\/λ; 
rf = 
            Rest[FoldList[Plus, 0, 
                jp]]; 
selectedstep = \(Position[
                rf, \(Select[rf, # > Random[] &, 
                    1]\)[[1]]]\)[[1, 1]]; 
jump = 
            reactionvectors[[
              selectedstep]]; 
AppendTo[
            jumps, jump]; state += jump]; 
Plot[
          Evaluate[
            initialstate + 
              Transpose[
                  jumps] . \((\(UnitStep[x - #] &\) /@ 
                    Rest[FoldList[Plus, 0, 
                        waitingtimes]])\)], 
{x, 0, 
            1.01\ timelimit}, opts, 
          PlotStyle -> {{Thickness[0.015], Red}}, 
          AxesLabel -> {time, \ "No. of mol."}];
]\)

Suppose we are given that it is the substrate S (second species, here) that should start with as many as 10 particles.

Kurtz[{{1, 0, 0}, {1, 0, 0}, {0, 1, 1}, {0, 0, 0}}, {{0, 1, 1}, {0, 1, 
          0}, {1, 0, 0}, {0, 0, 1}}, {10.0, 1.0, 1.0}, {0.01, 0.001, 0, 0}, 
      2, 10, (10^(-6))]

The curves (not printed here) show the similar qualitative characteristics as the deterministic solutions.

A short description of Mathematica programs useful in formal reaction kinetics has been given. The interested reader will find extended version of the present paper, including more programs and more references on the homepage of the author: www.math.bme.hu/~jtoth

Acknowledgements

The present work describes tools developed mainly for use in toxicokinetics and supported by the French Ministry of Territorial Planning and Environment (BCRD AP2001) . Partial support by WRI is also acknowledged. The author is grateful for useful discussions on the topics with F. Bois (INERIS, Verneuil-en-Halatte), M. Moreau (PMCU, Paris), J.-P. Rospars (INRA, Versailles), B. Vizvári (ELTE, Budapest), and Ms. Céline Brochot (INERIS, Verneuil-en-Halatte), and also his students, Ms. Andrea Halmschlager and Mr. Z. Rózsa.

Érdi, P., Tóth, J.: Mathematical models of chemical reactions, Theory and applications of deterministic and stochastic models. Princeton University Press, Princeton, 1989.

Feinberg, M.: Mathematical aspects of mass action kinetics. In: Chemical reactor theory: A review, eds. N. Amundson, L. Lapidus, Prentice Hall, Englewodd Cliffs, NJ, pp. 1-78.

Sipos, T., Tóth, J., Érdi, P.: Stochastic simulation of complex chemical reactions by digital computer. I. The model, React. Kinet. Catal. Lett. 1 (1974) 113-117.

Tóth, J., Rospars, J.-P.: Dynamic modeling of biochemical reactions with applications to signal transduction: principles and tools (in preparation)

Vol'pert, A. I., Hudyaev, S. I.: Analysis classes of discontinuous functions and equations of mathematical physics.  Nauka. Moscow. 1975 (In Russian, English version: Dordrecht, Nijhoff, 1985).