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 has been developed 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.

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]]\)

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

It would be nicer to formulate these rules together rather than separately.

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). We also create names to describe the quantities of the chemical species.

\(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], 
	reactionSteps -> ReactionSteps, 
          R -> Length[ReactionSteps], 
	α -> 
            alpha, β -> beta, 
	
          variables -> ff /@ Species}]\)

st = ToStoichiometry[can]

One may observe how useful are the list processing functions and term rewriting.

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.

Options[RHS] = {MassAction -> True};
\(RHS[stres_, opts___?OptionQ] := 
      Module[{w, rate}, 
Print[reactionSteps /. 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]]\)

rhs = RHS[st]

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 repeated term rewriting.

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 possible to do numerically.

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}}];

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[α /. st, β /. st, species /. st, 
      reactionSteps /. st, {Ε, S}]

The usual stochastic model of a chemical reaction is a Markovian pure jump process of the following 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 wherein 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. A minor issue is 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.

\(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).