Initial condition values in moles of carbon for the 10 compartment model shown in Figure 1. The flux is the sum of the flows into and out of each pool during the first time step.
Northern elephant seals (
Mathematical models of biochemical systems are a prerequisite for a true understanding of the complexity of metabolic and physiologic systems. A model can be defined in both a physical and mathematical sense as a set of equations that describe the behavior of a dynamic system, and the response of the system to a given stimulus (Jeffers 1982). Many types of models exist. The models that most closely approximate reality are often the most complex, and it is often difficult to derive unbiased or valid estimates of model parameters. Matrix models offer a way to sacrifice some of the “reality” to gain the advantages of mathematical deduction and prediction (Jeffers 1978). Matrix models are ideally suited to simulate the results of isotope tracer experiments and linear compartment analysis models (Shipley and Clark 1972).
2. Materials and methods
The amount of carbon residing in each pool in Figure 1, and the fluxes between the pools, were determined using single injection radiotracer methods as described (Pernia,
where k is the fractional turnover rate, t is time, and e is the base of the natural logarithms. These parameters may be estimated by injecting a known amount of tracer (q) labeled in some way (such as 14C or 3H) but which is metabolically indistinguishable from the metabolite of interest (tracee).
Blood samples are taken over time, and the specific activity of the tracer determined and plotted on semi-log paper. If the assumption of first order kinetics holds, the plot will be linear, and the slope of the line is k. The assumption of instantaneous mixing allows extrapolation of the line back to time = 0, and the estimation of the specific activity at time = 0 (SA0). The size of the pool (Q) can then be estimated (Katz, et. al. 1974):
The magnitude of the entry rate (R0 ) can then be estimated:
If the specific activity curve is plotted on regular paper and integrated, the Stewart-Hamilton equation provides a stochastic estimate of the irreversible loss rate (L) (Shipley and Clark 1972):
Recycling rate (R) is the part of the entry rate (R0) which leaves the pool of interest and returns to it during the experiment. It is the difference between the entry rate (R0) and the irreversible loss rate (L) (Nolan and Leng 1974):
Once the sizes of the pools and the flux rates between them are known, differential equations can be written to describe the rate of change of each compartment (Shipley and Clark 1972):
Where: Qj = Size of pool j (source of flux to pool i).
kij = Rate constant for flow to i from j.
Qi = Size of pool i (pool of interest).
kii = Sum of all rate constants for flows from compartment i:
If the coefficients of the differential equations (kij) are entered into a source/destination matrix (A), the matrix may be used to predict the magnitude of the flux between each compartment:
Where F is the vector of the sum of the fluxes to and from all pools, Q0 is the vector of initial pool sizes, and A is the coefficient matrix derived above. The future state of the system then becomes:
Iteration of the above two operations allows computation, in difference equation format, of the future state of the system at any time desired. However, use of the matrix exponential permits analytical determination of the state of the system at any future time in one operation:
Initially the system dynamics were simulated for an 8-week period, approximately as long as the duration of the fasting period of the northern elephant seal pups. In most cases, the simulation revealed monotonic (linear) declines in the sizes of the pools, with the obvious exception of the sink pool. However, in the case of the carbohydrate pools, i.e. glucose, glycerol and lactate, the pool sizes increased rapidly during the first part of the simulation, and then declined monotonically. This suggested that the estimates for the initial condition sizes of the pools were too low. In order to correct for this, a second, shorter (24 hr) simulation was conducted, after correcting the initial condition sizes of these pools by extrapolating the linear part of the pool size decline later in the simulation back to time zero.
A matrix may be described in terms of its characteristic equation, which will have the same order as the number of rows (= columns) in a square matrix (Jeffers 1978, Swartzman 1987). The roots of this characteristic equation are the eigenvalues (λ), which can be used to assess the stability of the system described by the matrix (Heinrich et al. 1977, Edelstein-Keshet 1988). Simply put, if all of the eigenvalues, or their real parts, are negative, the system is stable. If one or more eigenvalues are positive, the system is unstable. Zero value eigenvalues indicate a closed system (Edelstein-Keshet 1988, Halfon 1976, and Swartzman 1987). Other matrix parameters relevant to stability analysis are the trace (τ) and determinant (Δ) of the matrix. The trace is the sum of the eigenvalues, and the determinant in the product of the eigenvalues (Heinrich et al. 1977). A dimensionless τ-Δ parameter plane may be envisioned which relates the magnitude of these two parameters to model stability, or the type of instability (Edelstein-Keshet 1988). All matrix calculations described herein were conducted using MATLAB®.
Table 1 shows the initial conditions of the vector Q0 which contains the initial conditions of the system, and the first vector of fluxes calculated by multiplying Q0 by the source/destination matrix A, which is contained in Table 2. The characteristic equation of this matrix is: 0 = 1.0x10 – 0.186x9 + (9.34x10-3)x8 + (9.50x10-5)x7 + (3.78x10-7)x6 + (6.45x10-10 )x5 + (4.18x10-13)x4 + (7.95x10-17)x3 + (7.88x10-22)x2 + (2.15x10-28)x + 0.
|2. Glycerol||2.3100e-0 01||1.3800e-005|
The roots of this equation are the eigenvalues (λ) of the matrix: 0.00; -8.10x10-4; -4.89x10-3 ; -2.70x10-3; -3.25x10-3; -1.06x10-1; -1.02x10-5; -3.00x10-4; -6.80x10-2; -2.81x10-7. Notice that there is one zero value, indicating that this is a closed system. The remainder of the eigenvalues are all negative, indicating a stable system. The trace (τ), or sum of the eigenvalues, of matrix A is -0.186 and the determinant (Δ), or product of the eigenvalues of matrix A is zero. Figure 2 shows the dimensionless τ-Δ parameter plane. The point for matrix A would lie on the y-axis, just below the x-y intercept, indicating that the system approaches a saddle point of stability. Notice that τ2 is 0.0346, which is greater than 4Δ, which is zero, indicating again that the matrix lies in the portion of the phase-plane corresponding to a saddle point condition (Heinrich et al. 1977).
Figure 3 shows the changes in the size of the adipose tissue and sink pools over an 8-week simulation. These pools are plotted together because they were of similar size. As expected, the adipose tissue pool declined at a constant rate, while the sink pool accumulated carbon asymptotically.
Figure 4 shows the changes in the size of the various lipid and lipid-derived pools over an 8-week simulation. The ketone pool rose slightly early on, and then declined slowly, as did the palmitate pool. At this resolution, it appeared that the glycerol pool declined at a constant rate throughout the simulation but when compared to the glucose and lactate simulations (Figure 5) it became clear that this pool was poorly initialized.
Figure 5 shows the changes in carbohydrate and carbohydrate-derived pools over an 8-week simulation. Apparently the glucose and lactate pools were initialized incorrectly, as reflected in their rapid increase in the first week of the simulation. Once they became stable, they both declined at the same rate, which was equivalent to the rate of decline of the glycerol pool. This similarity in rates of decline suggests that these pools are closely linked through recycling pathways, and may in fact represent a subsystem of the model.
Figure 6 shows the changes in nitrogen containing pools over an 8-week simulation. The tissue protein pool declined at a very slow rate, suggesting almost no protein catabolism. This is supported by the decline and continued low level of the urea and alanine pools. The continued low urea level indicates again that protein catabolism is occurring at a very low rate.
Figure 7 shows the time course of the contents of the palmitate, glycerol, and ketone pools over a 24-hour simulation. Figure 8 shows the time course of the contents of the alanine, urea, and tissue protein pools over an 8-week simulation. Notice that the tissue protein pool declined almost imperceptibly, while the urea and alanine pools declined initially and then became stable.
Figure 9 shows the time course of the contents of the glucose, glycerol, and lactate pools over a 24-hour simulation. The glucose and lactate pools were apparently poorly parameterized, and grew rapidly to a zenith, and then declined steadily. The glycerol pool declined continuously.
Figure 10 shows the time course of the contents of the glucose, glycerol, lactate, alanine and ketones pools after re-estimation of the initial conditions by extrapolating the linear parts of Figures 7, 8 and 9 to time zero.
As initially formulated, the value of the rate constant for efflux from the sink pool (k10,10) was set at zero, because there should be no efflux from the sink. Eigenanalysis of this matrix yielded 9 negative eigenvalues and one zero eigenvalue, indicating a stable, closed system (Keith 1999). An eigenvalue of zero is to be expected from a matrix with a zero on the main diagonal (Edelstein-Keshet 1988, Swartzman 1987). For this reason, the zero eigenvalue can be considered an artifact of matrix construction, and not truly representative of the system being simulated. Therefore, a second analysis was conducted in which the value of the rate constant for efflux from the sink pool was set to one, which allowed all of the contents of the sink pool to exit at every time step. Eigenanalysis of this matrix yielded 10 negative eigenvalues, indicative of a now open and still stable system. The trace (τ) of this matrix was -0.186, and the determinant (Δ) was 2.174 x 10-32. In this case τ2 was 0.0346 and 4Δ was 8.696 x 10-32. Taken together, these values indicate that the system without a sink lies near a stable node condition in the phase-plane. However, this is difficult to reconcile with the biological reality that a fasting elephant seal with no food or water inputs is not at equilibrium, and cannot survive forever (Ortiz et. al. 1978).
The genesis of this apparent contradiction may lie in differences in time scale or time constants. Differential equations with widely different time scales are “stiff” (Heinrich et al. 1977) and will have eigenvalues of different orders of magnitude. This is apparent here where the eigenvalues range from -1.06 x 10-1 to -2.81 x 10-7. The reciprocals of the eigenvalues are the relaxation times (Heinrich et al. 1977) and these likewise vary over six orders of magnitude, indicating that there are fast-reacting variables and slow-reacting variables in the simulation. Such hierarchical time structure may obscure predictions of model stability because the eigenvalues only characterize the system in the close time-neighborhood of the steady state where linear approximation is appropriate (Heinrich et al. 1977). Thus, predictions of model stability based on the signs of the eigenvalues may contradict a prediction of model instability based on relaxation times and slow-moving versus fast-moving variables in the system (Heinrich et al. 1977) over the duration of the actual fast of the animal.
Elevated palmitate levels are consistent with field data indicating that the major energy substrate during fasting (Castellini, et. al. 1987, Keith 1984). The decline in ketone levels through the simulation is consistent with field data indicating low levels of ketone bodies in the plasma of fasting northern elephant seal pups (Costa and Ortiz 1982). Declines in alanine, tissue protein, and urea levels are also consistent with field data (Pernia,