Computer Simulation of Metabolism
Introduction
Sims and Folkes (1964) [and references cited therein] describe kinetic equations for calculating the isotope abundance of intermediates in a linear metabolic pathway (A > B > C > D), when the precursor A is
supplied at a known isotope abundance (A') (atom % excess) at time (t) = 0 (Sims, A.P., Folkes, B.F. 1964. A kinetic study of the assimilation of [^{15}N]ammonia and the synthesis of amino acids in an exponentially growing culture of Candida utilis. Proc. Roy. Soc. Lond. B. Biol. Sci. 159: 479502). These equations are as follows:
B' = A'[1e^{bt}]
C' = A'[(c(1e^{bt})b(1e^{ct}))/(cb)]
D' = A'[(cd(cd)(1e^{bt})bd(bd)(1e^{ct})+bc(bc)(1e^{dt}))/((cd)(bd)(bc))]
where :
A' = isotope abundance of precursor (atom %)
B' = isotope abundance of pool B (atom %)
C' = isotope abundance of pool C (atom %)
D' = isotope abundance of pool D (atom %)
b = turnover constant of pool B (h^{1}) [rate of synthesis.pool size^{1}])
c = turnover constant of pool C (h^{1}) [rate of synthesis.pool size^{1}])
d = turnover constant of pool D (h^{1}) [rate of synthesis.pool size^{1}])
e = exponential function
t = time (h)
To illustrate the application of these equations, consider the hypothetical pathway shown in Fig. 1, where we assume a precursor isotope abundance of A' = 90%, pool sizes of B, C and D (B2, C2 and D2, respectively) each of 500
nmol.gfw^{1}, and rates (B3, B4, C3, C4, D3 and D4) as indicated in Fig. 1, in units of nmol.h^{1}.gfw^{1}. The turnover constants of pools B, C, and D are b =
10, c = 4, and d = 2 (h^{1}), respectively (Fig. 1) [gfw = gram fresh weight].
Fig. 1. Model assumptions used to generate simulations shown in Figs. 2, 3 and 4.
In this specific example, A', B', C', and D' are plotted as a function of time (t) over a 1.0 h time interval in Fig. 2.
Fig. 2. Simulations of A', B', C' and D' using kinetic equations
It is important to note that these equations are applicable only to steadystate situations where pools and rates remain constant with time, and are not applicable to the "chase" phase of classical pulsechase labeling experiments.
An alternative approach to simulating this labeling behavior of intermediates in a pathway is to use a simple iterative computer model of the type shown below:
z = Mx / 2000
For t = 0 To Mx Step z
B1 = (B1 * B2 + A1 * z * B3) / (B2 + z * B3)
B2 = B2 + z * B3
B1 = (B1 * B2  B1 * z * B4) / (B2  z * B4)
B2 = B2  z * B4
C1 = (C1 * C2 + B1 * z * C3) / (C2 + z * C3)
C2 = C2 + z * C3
C1 = (C1 * C2  C1 * z * C4)/(C2  z * C4)
C2 = C2  z * C4
D1 = (D1 * D2 + C1 * z * D3)/(D2 + z * D3)
D2 = D2 + z * D3
D1 = (D1 * D2  D1 * z * D4)/(D2  z * D4)
D2 = D2  z * D4
[goto subroutine to plot A1, B1, C1 and/or D1 as a function of t]
Next t
where :
A1 = isotope abundance of precursor (atom %) [supplied]
B1 = isotope abundance of pool B (atom %)
C1 = isotope abundance of pool C (atom %)
D1 = isotope abundance of pool D (atom %)
B2 = pool size of pool B (nmol.gfw^{1})
C2 = pool size of pool C (nmol.gfw^{1})
D2 = pool size of pool D(nmol.gfw^{1})
B3 = rate of synthesis of pool B (nmol.h^{1}.gfw^{1})
C3 = rate of synthesis of pool C (nmol.h^{1}.gfw^{1})
D3 = rate of synthesis of pool D (nmol.h^{1}.gfw^{1})
B4 = rate of utilization of pool B (nmol.h^{1}.gfw^{1})
C4 = rate of utilization of pool C (nmol.h^{1}.gfw^{1})
D4 = rate of utilization of pool D (nmol.h^{1}.gfw^{1})
t = time (h)
z = iteration interval (in this case 1/2000^{th} of total time, Mx = scale of Xaxis (h))
A1, B1, C1 and D1 are plotted as a function of time over a 1.0 h time interval in Fig. 3 using the same rates and pool sizes as assumed in Fig. 2 (see Fig. 1 for details). Note that the output from the iterative model (Fig. 3) is virtually identical to that of the kinetic equations (cf. Fig. 2) (i.e. A' = A1, B' = B1, C' = C1, D' = D1).
Fig. 3. Simulations of A1, B1, C1 and D1 using iterative computer model
As will be described in subsequent pages, important advantages of the iterative model over the kinetic equations described above, are that the iterative model can be used in nonsteadystate situations (e.g. where rates are not constant, and where pools expand or deplete with time), or can be readily applied to consider multiple pools of intermediates with different turnover rates. The iterative model can also be easily adapted to handle the "chase" phase of classical pulsechase labeling experiments. For example, the data of Fig. 4 simulates a "chase" initiated at t = 0.5 h, where the precursor isotope abundance A' = A1 is reduced from a value of 90% to 0% at t = 0.5 h. This is achieved simply by adding the statement "If t > 0.5 Then A1 = 0"
within the "For ... Next" loop.
Fig. 4. Simulations of A1, B1, C1 and D1 assuming a "chase" at t = 0.5 h
The examples above are given in the context of their original use with stableisotopes (e.g. ^{15}NH_{3} supplied to cells, monitoring the ^{15}Nabundance (atom %) of amino acids
derived therefrom) (Sims and Folkes (1964)). However, as will be described in subsequent pages, the iterative computer model is easily applied to radioisotope labeling experiments, where A1, B1, C1, D1 would represent the specific radioactivities (e.g. nCi.nmol^{1}) of the intermediates. The total radioactivity (nCi) in any one intermediate at any one time would be equal to the pool size multiplied by its specific activity (e.g. total radioactivity in pool B would equal B1 * B2 (nCi)).
