## Multiplecoupled reactions

What if the reaction sequence is more extensive? Is it possible to solve the simultaneous differential equations analytically? The answer is "Yes!" but only in general if the reactions are simple first-order ones. Indeed Mathematica solves these with ease. Consider, for example, the sequence of five reactions with six reactants:

Integrating the six differential equations 'by hand' is very tedious (e.g., reference®), but Mathematica readily yields the exact analytical solutions.

Q: Derive and solve the system of differential equations resulting from the reaction scheme shown in Eqn [1.19].

A: Use DSolve in the same manner as for the previous example.

solution = DSolve[{ a' [t] & -ki a[t] , b' [t] & k1 a[t] - k2 b[t] , c' [t] & k2 b [t] - k3 c [t] , d' [t] & k3 c [t] - k4 d[t] , e' [t] & k4d[t] - k5e[t] , f' [t] & k5e[t] , a & a0, b  & 0, c & 0, d & 0, e & 0, f  & 0} , {a[t] , b[t] , c[t] , d[t] , e[t] , f[t]}, t];

And, with the following initial conditions, we can plot the results.

k1 = 0.5; k2 = 0.55; k3 = 0.6; k4 = 0.65; k5 = 0.7;

Plot[Evaluate[{a[t], b[t], c[t], d[t] , e[t] , f[t]} /. solution], {t, 0, 10}, PlotRange ^ All, AxesLabel -> {"Time ", "Concentration" }] ;

Concentration Figure 1.6. Simulation of the time course of the reaction scheme described by Eqn [1.19] using DSolve with the corresponding differential rate equations. The ordinate denotes concentration in mmol L-1 and the abscissa, time in seconds.

Figure 1.6. Simulation of the time course of the reaction scheme described by Eqn [1.19] using DSolve with the corresponding differential rate equations. The ordinate denotes concentration in mmol L-1 and the abscissa, time in seconds.

Notice that in the solution to this problem we use a different method for evaluating and plotting the results generated by DSolve. Instead of using the replacement rule generated by DSolve to define new functions, we use a combination of the Evaluate function and the /. operator directly within the Plot function itself. This is often a more efficient way of using results which are expressed as replacement rules. In order to evaluate the results of a replacement rule at a specific time (say, t = 4) without defining a new function, we can use the following commands:

{a[t] , b [t] , c[t] , d [t] , e[t] , f [t]} /. solution/. t -> 4

{{1.35335, 2.45321, 2.44581, 1.7734, 1.04475, 0.929452}}

PlotfEvaluate/[x] /. rule], Plotting solutions given as replacement rules {x, xmin, xmax<]

/[x] = /[x] /. rule /. x -> value Evaluating solutions given as replacement rules

### Using solutions expressed as replacement rules.

From the above example we see that the analytical expression for each of the reactants in the scheme in Eqn [1.19] is a function of time which is mathematically described as a sum of exponentials. These complicated expressions can be viewed by removing the at particular value of the independent variable end-of-line semicolons and re-executing the first Cell in the example. In a linear reaction system like that above, the exact number of exponentials in the expression for the concentration of a species is equal to 1 + the number of species that precedes it in the sequence. The exponents are a complicated series of sums of various products of the rate constants, and they are preceded by similarly complicated pre-exponential coefficients. The complexity evident in such analytical solutions does not bode well for a direct analytical assault on the kinetic representation of a metabolic system with hundreds of reactants!