All the metabolic pathways that we have considered so far have involved the assumption that the reactions take place in free solution and no consideration has been given to partitioning of reactants and enzymes into various compartments. To actually simulate such a system introduces another order of complexity into formulating the model. The principal 'trick' in such models is to express the rate equations in terms of the rate of change of amounts (moles) rather than concentrations. If we do this Eqn [4.3] becomes d (V.S)

dt where N is the stoichiometry matrix, v is the vector of reaction velocities containing rate equations expressed in terms of the time rate of change of amounts, S is the vector of substrate concentrations, and V is a diagonal matrix where the ith element contains the compartment volume of the ith substrate in S. Note that in this formulation, although the rate equations are expressed in terms of the rate of change of amounts, the rate equations are also expressed as functions of substrate concentrations. In other words, all we do for this analysis is to take the usual rate expressions and multiply them by the appropriate compartment volumes.

If we expand Eqn [4.14] we obtain dS dV

dt dt

Now for many systems of interest we can assume that the compartment volumes remain constant with time. In this case Eqn [4.15] simplifies to dS

dt or dS 1

This may seem a little abstract at the moment but it should become clearer with the following question/answer.

Q: Consider the following simple two-compartment model (Eqn [4.18]) consisting of an intracellular and an extracellular compartment with volumes Voli and Vole, respectively.

The model contains only two metabolites A and B; the subscript e denotes that the metabolite is in the extracellular compartment, while the subscript i refers to the intracellular compartment. Assume that all reactions are described by first-order rate constants. Simulate a time course for this system and determine the steady-state concentrations of A and B in each compartment.

A: The reaction scheme represented in Eqn [4.18] can be described by the following equation list:

{Ai[t] ^ Ae[t]}, {Ae[t] ^ Be[t]}, {Be [t] ^ Bi[t]}, { Bi[t] ^ Ai[t]}};

We model each of these transport and metabolic reactions with simple first-order rate equations as indicated in the question. Note that because we want to express the reaction/transport rates in terms of the number of molecules per unit time, it is necessary to multiply each rate equation by the appropriate compartment volume. Suppose that the compartment volumes are

Hence, the requisite rate equations are v[1] := ki Ai[t] Voli; ki = 1;

Now use the MetabolicControlAnalysis functions to derive the matrices and vectors of this system, in the manner used above.

S = SMatrix[eqns] ; N = NMatrix[eqns] ; ■v = VMatrix[eqns] ; S // MatrixForm N // MatrixForm ■v // MatrixForm

1 |
-1 |
0 |
0 |

-1 |
0 |
0 |
1 |

0 |
1 |
-1 |
0 |

0 |
0 |
1 |
-1 |

It remains to define the V matrix of Eqn [4.14]; this is simply the diagonal matrix where the ith element contains the compartment volume of the ith substrate in S. Hence for this system, V is given by

V // MatrixForm

When constructing this matrix, make sure that the order of the volumes in the diagonal matrix matches the order of the substrates in the substrate list.

Having constructed all the appropriate matrices we can solve this system of matrices and vectors by using the NDSolveMatrix function in the manner described above, except that we use V-1 N instead of simply N.

{Ai [0] + 1, Ae[0] + 0, Bi[0] + 0, Be[0] + 0}, {t, 0, 10}]

{{Ae ^InterpolatingFunction[{{0., 10.}}, <>] ,

Ai ^InterpolatingFunction[{{ 0., 10.}}, <>] ,

Plot [Evaluate [S /. sol] , {t, 0, 10}, PlotRange-> {0, 1.5}, AxesLabel -> {"Time", "Concentration"}];

Concentration

Concentration

Figure 4.6. Time course of the reaction scheme shown in Eqn [4.18]. Ai delines from an initial concentration of 1 to a steady state of 0.25 mmol L-1, while Bi approaches a steady state of 0.25 mmol L-1 from an initial concentration of 0 mmol L-1. Both [Ae ] and [ Be ] approach a steady state of 0.75 mmol L-1 with [Ae ] initially overshooting this steady-state concentration by a modest amount.

Figure 4.6. Time course of the reaction scheme shown in Eqn [4.18]. Ai delines from an initial concentration of 1 to a steady state of 0.25 mmol L-1, while Bi approaches a steady state of 0.25 mmol L-1 from an initial concentration of 0 mmol L-1. Both [Ae ] and [ Be ] approach a steady state of 0.75 mmol L-1 with [Ae ] initially overshooting this steady-state concentration by a modest amount.

From this simulated time course it is seen that the system approaches a steady state. To determine the values of the concentrations we use the function SteadyState; however, we must first determine the conservation of mass relations; these are given by the following function from the MetabolicControlAnalysis package.

ConservationRelations[S, Inverse[V ] .N]

{Ae[t] + 3.Ai[t] + 1.Be[t] + 3.Bi[t]} == {Const [ 1 ]}

From the initial conditions used in the time course simulation, we must set Const[1] = 3.

Hence, the steady state of the system is given by the following function evaluation: SteadyState[S, Inverse[V ] .N, V]

{{Ae[t] ^ 0.75, Ai [t] ^ 0.25, Be[t] ^ 0.75, Bi[t] ^ 0.25}}

In the above example we assumed that cell volumes were constant. In many situations of interest this may not be the case. For example, most cells are unable to sustain a high transmembrane difference in osmotic pressure and, if this occurs, the cell volume will change.

To simulate metabolic systems in which compartment volumes change we could resort to using Eqn [4.15]. However, another method, which often makes model formulation easier, is to represent the changes as changes in amounts rather than concentrations. When doing so the familiar equations are still employed:

However, in this case, S, is a list of substrate amounts (not concentrations), while the rate equations in v are written in terms of the rate of change of amounts. The stoichiometric matrix, N, remains unchanged. This overall approach to modelling these systems is illustrated in the following question/answer.

Q: Simulate the time dependence of the concentrations of Na, Ca, and K in a suspension of red blood cells in which a calcium ionophore is added to the cells. The ionophore stimulates K efflux by delivering Ca into the cells, thus 'switching on' the calcium-activated (Gardos) K channel.

A: In this problem it is assumed that a cell is not able to sustain a substantial difference of osmolality across its membrane. Hence, as ions flow across the membrane, water distribution is altered and this is reflected in a change in the cell volume. Furthermore, since the metal cations permeate the cell membrane, but charged intracellular proteins do not, a potential difference is generated across the membrane. This is called the Donnan potential. Its value is calculated by using the Nernst equation (see below). A graphical illustration of the reaction scheme is as follows:

Figure 4.7. Model of a suspension of cells that take up calcium that, in turn, activates the loss of potassium from the cells. The calcium influx is also accompanied by exchange of sodium between the inside and the outside of the cells. The scheme is a model of human erythrocytes in which calcium activates the Gardos channel. The model consists of a Michaelis-Menten membrane transport protein for each of the cationic species, while it is assumed, as is the case with human erythrocytes, that the passive exchange of anions (chloride and bicarbonate) is very rapid (on the sub-second timescale) via the band 3 (capnophorin) transport protein.

Figure 4.7. Model of a suspension of cells that take up calcium that, in turn, activates the loss of potassium from the cells. The calcium influx is also accompanied by exchange of sodium between the inside and the outside of the cells. The scheme is a model of human erythrocytes in which calcium activates the Gardos channel. The model consists of a Michaelis-Menten membrane transport protein for each of the cationic species, while it is assumed, as is the case with human erythrocytes, that the passive exchange of anions (chloride and bicarbonate) is very rapid (on the sub-second timescale) via the band 3 (capnophorin) transport protein.

The first step in modelling this system is to define carefully all the initial volumes in the system. Suppose that the total sample volume is 3 x 10-3L. Hence we define the values as follows:

We now specify the fraction of the sample volume that is occupied by the cells. This fraction is called the hematocrit. Note that the hematocrit potentially changes in time as the ions exchange across the cell membrane and an input value is the haematocrit at zero time. Thus, ht[0] = 0.65;

The volume fraction of an erythrocyte that is occupied by haemoglobin and the cytoskeleton is usually taken to be 0.29 of the volume of the cell of normal volume. When the cell changes volume, the volume of these proteins stays the same but their volume relative to the cell volume changes. Thus the fraction of the erythrocyte volume that is occupied by haemoglobin is given by volhbn = ht [0] x 0.29 x voltotal;

Now use these values to calculate the initial water volumes in the intracellular (voliw[0]) and extracellular (vole[0])space of the sample.

vole[0] = (1 - ht[0]) X voltotal; voli [0] = ht [0] X voltotal; voli,w[0] = voli [0] - volhbn;

Having defined the initial water volumes in the two compartments of the sample, we must determine how these two volumes change with time. Before we do this it is necessary to define the ion movements that will lead to these volume changes.

The system of transport processes that is described in Fig. 4.7 can be represented with the following equation list:

Each of the cationic transport processes is modelled by a Michaelis-Menten rate equation. For simplicity we assume that intracellular calcium activates the influx and efflux of both sodium and potassium ions. The rate equations and their associated parameters are as follows:

ki [ t ] + K /1 + Ka c vo 1 i , „ ] t ] "li-Wjti- + Km,ki i1 + ca i [t ]-

Vmax'nai voliw[t]

v[naetransD :

v[caetrans] "

Note that each rate equation is specified in terms of the rate of change of amount (units: mol s-1) rather than the rate of change of concentration. Because of this the Vmax parameters in the above equations have units of mol s-1 and are made up as follows: number of molecules of transporter per cell (N) x turnover number per active site of the transporter x number of cells in the sample. The total number of cells in the sample is given by voltOtal x ht[0] / volume of a single cell .

Hence, when determining Vmax values from the literature, we must pay special attention to its physical units. The Vmax values we needed for the above equations are reported in terms of mol (L internal volume)-1 s-1so it is necessary to convert them to the appropriate units by multiplying by the internal volume of the cell.

Having defined the systems of reactions and the associated rate equations, we are now able to define the S, N, and v matrices, and the initial conditions.

S = SMatrix[eqns] ; N = NMatrix[eqns]; v = VMatrix[eqns];

IC = {nag [0] == 0.14 x vole [0] , nai[0] == 0.005 xvoli/W [0] , ke [0] == 0.010 x vole [0] , ki [0] == 0.14 x voli,w [0] , cae [0] == 0.002 xvole [0] , cai [0] == 0.0001 xvoli/W [0]};

The final task before we can simulate a timecourse is to define how cell volume will change with changes in ionic concentrations. Our model of volume change requires the following major assumption: that the passive exchange of anions is rapid compared to the rates of transport of other ions. With this assumption we know that chloride ions (and other anions such as bicarbonate) will distribute between the intracellular and extracellular spaces so that electroneutrality is maintained in each compartment. Thus the following equation must hold:

cle[t_] : = nag [t] + ke[t] + 2 ca<, [t] ; cli [t_] : = nai [t] + ki [t] + 2 cai [t] - 4hb;

Note that we have assumed that the charge on each haemoglobin molecule is - 4. In addition we know that the erythrocyte membrane is not able to withstand an osmotic pressure difference of any great magnitude. So it is appropriate to assume that the osmotic pressure inside the cell is equal to the osmotic pressure outside the cell. This assumption is equivalent to assuming that the sum of the concentrations of all solutes inside is equal to the sum of all those outside, under conditions of thermodynamic 'ideality'. Hence, the following equation will yield the intracellular volume as a function of the concentrations of the ionic species inside and outside the cell:

(voltotai,w[0] (k±[t] + nai [t] + cai[t] + cl±[t] + hb) )/ (k±[t] + nai [t] + cai [t] + cli [t] + hb + ke [t] + nae [t] + ca« [t] + cle [t]) ;

where the mass (mol) of haemoglobin, hb, is calculated by using its partial specific volume (0.73 L kg-1 ) and its molecular weight (64.5 kg).

0.73x 64.5

The volume of aqueous medium outside the cells is simply the total solute-accessible volume in the sample minus the aqueous volume inside the cells.

We are now in a position to solve the time-dependent behaviour of this system using NDSolveMatrix and then to view the time courses of concentrations or amounts graphically.

res = NDSolveMatrix [S, N, v, IC, {t, 0, 5000} , MaxSteps -> 5000] ;

{t, 0.0, 5000} , AxesLabel ^ {"Time (s) ", "K+ (mol)"}, PlotRange^ {0.0, 0.0002}];

0.0002r

0.000125

0.000175

0.00015

0.0001

0.000075

0.00005

0.000025

1000 2000 3000 4000 5000

Figure 4.8. Simulated time course of the amount of intracellular potassium ion in a cell suspension for the reaction scheme shown in Fig. 4.7.

Plot [Evaluate [{cai [t] /voli/W[t] , ca [t] /vole[t]} /. res] , {t, 160, 320}, PlotRange 0 All,

AxesLabel 0 {"Time (s) ", "[Ca+]in and [Ca+]out(mol L-1)"}];

Figure 4.9. Simulated time course of the concentration of intra- (lower curve) and extracellular calcium ions for the reaction scheme shown in Fig. 4.7.

Plot [Evaluate [{voli w [t] + vol^n , vole [t]} /. res] , {t, 0, 5000}, PlotRange0 All,

AxesLabel 0 {"Time (s) ", "Volin and Volout (L) " }] ;

Figure 4.9. Simulated time course of the concentration of intra- (lower curve) and extracellular calcium ions for the reaction scheme shown in Fig. 4.7.

Volin and Volr

0.0018 0.0017 0.0016 0.0015 0.0014 0.0013

!000 3000 4000 5000

Figure 4.10. Simulated time course of the intra- (upper curve) and extracellular volumes for the reaction scheme shown in Fig. 4.7.

Note that the volumes do not reach the same value, primarily because the hemoglobin and other proteins inside the cells occupy volume.

From these results we can also calculate the membrane, based on the anion distribution ratio, by using the Nernst equation as follows:

rR= 8.314; tT= 310; fF = 96500; parNernSt = rR tT / fF;

gph4 = Plot [1000 parNernSt * Log [Evaluate [ cli [t] /cle[t] /. reS]] , {t, 0, 5000}, PlotRange ^ All,

AxeSLabel ^ {"Time (s)", "Membrane Potential(mV)"}];

Membrane Potential(mV)

Figure 4.11. Simulated time course of the membrane potential for the reaction scheme shown in Fig. 4.7.

Figure 4.11. Simulated time course of the membrane potential for the reaction scheme shown in Fig. 4.7.

Was this article helpful?

## Post a comment