## Using Matrix Notation in Simulating Metabolic Pathways

For large metabolic networks, it is often simpler to express the set of simultaneous differential equations that describe the reaction system in matrix form. The use of matrix notation is also useful for analysing the existence and stability of steady states and for performing MCA (see Chapter 5).

Thus the set of differential equations describing the rate of production and utilization of all metabolites, S,, in a metabolic network of reactions is given by where vj (j = 1,...,r) is the rate of reaction j, and nij is the stoichiometric coefficient with i and j referring to the metabolite and reaction, respectively.

In matrix notation, Eqn [4.2] is given by where v and S denote the vectors of reaction rate expressions and concentrations, respectively, while N is termed the stoichiometry matrix with each element being «j. The stoichiometry matrix has one row for each reactant and one column for each reaction in the scheme. The number at the intersection of a row and a column defines the number of molecules of that reactant to engage in the specified elementary reaction, namely, its stoichiometry. Although this representation may seem overly complex at present, expressions like Eqn [4.2] and [4.3] are almost invariably much simpler than they look at first sight because in a large model most of the «j are zero and nearly all the rest are either - 1 or 1.

To analyse systems of differential equations that are represented in matrix notation requires a moderate amount of computer programming. However, once the general procedure has been set up it can be used for any reaction scheme. Because of its general usefulness we have written an add-on package called MetabolicControl-Analysis that contains the matrix representation of reaction schemes and several other useful functions that are described below and in Chapter 5. The details of the programming are not given in the text but the programs in Appendix 2 contain many comment statements to help you, the reader, understand the algorithms. Before this package is used for the first time, the reader must evaluate all the Cells in Appendix 2

so that the appropriate .m file is created and the functions can be called in a

Mathematica session.

 << read in the add - MetabolicControlAnalysis on package MetabolicControlAnalysis NDSolveMatrix[S, N, v, Uses the functionNDSolve to find a numerical initial conditions, {t, tmin, tmax}] solution for the metabolite concentrations, S, with time in the range tmin to tmax, for a system of ordinary differential equations defined by the matrices S, N, and v, and subject to the initial conditions

Reading in the package MetabolicControlAnalysis.

### Reading in the package MetabolicControlAnalysis.

Q: By using the matrix representation of a reaction scheme, plot the 10 h time course of the concentrations of the reactants in the following linear sequence,

where the vi are simple irreversible Michaelis-Menten rate equations (Section 2.3.1) with all Fmax values set to 1 mmol L-1 h-1; all Km values are 1 mmol L-1 and S1 [0] = 1 mmol L-1, S2 [0] = S3 [0] = S4 [0] = 0 mmol L-1.

A: The first step in answering the question is to load the MetabolicControlAnalysis package.

<< MetabolicControlAnalysis*

Then the matrices and vectors of Eqn [4.3] must be defined. A simple way to define the concentration vector is as follows:

Similarly, the stoichiometry and rate equation vectors are defined as follows:

Vmax,1,f = 1; Vmax,2,f = 1; Vmax,3,f = 1; Km,1,Si = 1; Km,2,s2 = 1; Km,3,s3 = 1;

With the S, N, and v matrices defined and the parameter values set, NDSolveMatrix is used to obtain the numerical solution.

{S1 [0] == 1, S2[0] == 0, S3 [ 0 ] == 0, S4 [ 0 ] == 0}, {t, 0, 10}]

{{si 0InterpolatingFunction[{{ 0., 10.}}, <>] , s2 0InterpolatingFunction[{{0., 10.}}, <>] , s3 0InterpolatingFunction[{{0., 10.}}, <>] , s4 0InterpolatingFunction[{{ 0., 10.}}, <>]}}

The following is a plot of all the concentrations of all the Si that are listed in the vector S:

AxesLabel -> {"Time (h)", "Concentration (mM)"}, PlotRange -> {0, 1}] ;

Concentration (mM)

Concentration (mM)

Figure 4.2. Time course of the reaction scheme shown in Eqn [4.4]. The figure shows Si declining from an initial concentration of 1 mM and S4 accumulating to ~1 mM by the end of the time course. The metabolites S2 and S3 are seen to build up and then decline during the time course.

And here is the value of each of the reactants listed in S at 2 h: s /. t ->2 /. sol

The list of four numbers gives the concentrations in mM for the reactants in the order in which they are defined in S. Thus, Si is 0.28 mM, S2 is 0.33 mM, etc.

In the worked example above the parameter values were individually specified to be i, with units that depend on the nature of the parameter. By setting the parameter values with the = command, whenever Mathematica sees the parameter it automatically replaces it with the value. For some analyses, these parameter value assignments will need to be unset or cleared and we shall see examples of this later in Chapter 5. This can be tedious when simulating complex networks of reactions. To obviate this problem, the function NDSolveMatrix can accept parameter values as a separate parameter vector.

 NDSolveMatrixftf, N, v, Uses the functionNDSolve to find a numerical initial conditions, {t, tmin, tmax}, p\ solution for the metabolite concentrations, S, with time in the range tmin to tmax, for a system of ordinary differential equations defined by the matrices S, N, v, and the parameter matrix p

Finding numerical solutions to a system of differential equations with parameter values specified in a parameter matrix.

Q: Set all the Km values in the previous worked example to 0.5 mmol L 1 and re-simulate the time course using a matrix representation of the Km and other parameter values.

A: First the previous values should be cleared and then the matrix denoted p is specified. Clear[Subscript];

P = {{Vmax,1,f, 1}, {Vmax,2,f, 1}, {Vmax,3,f, 1}, {Km,i,s1, 0.5}, {Km,2,s2, 0.5}, {Km,3,s3, 0.5}} ; p // MatrixForm

Vmax,2,f 1

Vmax,3,f 1

The matrix has the parameter name in the first column and its numerical value in the second. Often, a more convenient way of inputting the parameter matrix is to do as follows:

p = {Vmax,1,f , Vmax,2,f , Vmax,3,f , Km,1,s1 , Km,2,s2 , Km,3,s3 } ;

pv = {1, 1, 1, 0.5, 0.5, 0.5}; p = Transpose[{p, pv}]

The numerical solution and graphical output of the system are then obtained by using sol = NDSolveMatrix[S , N, V,

{s1 [0] == 1, s2 [0] == 0, s3[0] == 0, s4[0] == 0}, {t, 0, 10} , p] ;

AxesLabel -> {"Time (h)", " Concentration (mM)"}, PlotRange -> {0, 1}] ;

Concentration (mM)

Concentration (mM)

Figure 4.3. Time course of the reaction scheme shown in Eqn [4.1]. S1 declines from an initial concentration of 1 mmol L-1 and S4 accumulates to ~1 mmol L-1 by the end of the time course. The metabolites S2 and S3 are seen to build up and then decline during the time course.

Note the decline in [S1 ] with time, the rise to maximum concentrations and then the decline of the two intermediate species, S2 and S3, and the sigmoidal monotonically strictly increasing concentration of S4.