A Matrix Model of Fasting Metabolism in Northern Elephant Seal Pups

1. Introduction

Northern elephant seals (Mirounga angustirostris) undergo regular periods of aphagia during their annual life cycle, as do many other phocids (Le Boeuf and Laws 1994). After nursing for about 30 days, the weaned pup fasts for 6-8 weeks, maintaining a fasting hyperglycemia, hyperlipidemia, hypoketonemia, and hypoinsulinemia (Champagne, et al. 2005). In most mammals, fasting is accompanied by hypoglycemia, and thus the fasting hyperglycemia in these animals is paradoxical. Previous studies of glucose metabolism in these animals (Keith and Ortiz 1989) indicate that the hyperglycemia results from both low rates of glucose utilization, due to very low insulin levels (Kirby, et. al.1987), as well as high rates of glucose carbon recycling through both lactate and glycerol. Other studies indicate that fatty acids are the major energy substrate during this time (Castellini, et. al. 1987), and that these animals conserve nitrogen by having very low urea turnover and excretion rates (Houser and Costa 2001). Figure 1 shows a 10 compartment conceptual flow diagram of metabolite flux in fasting northern elephant seal pups as simulated in this study.

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, et. al.1980, Keith and Ortiz 1989, and Castellini, et. al. 1987). If the animal is in steady state, the change in pool size (Q) for each compartment will be:

$${\text{Q}}_{(\text{t})}={\text{ Q}}_{(0)}{\text{e}}^{-\text{kt}}$$E1

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 ¹⁴C or ³H) 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 (SA₀). The size of the pool (Q) can then be estimated (Katz, et. al. 1974):

$$\text{Q }=\frac{\text{q}}{{\text{SA}}_0}$$E2

The magnitude of the entry rate (R₀ ) can then be estimated:

$${\text{R}}_0=\text{ k }*\text{ Q}$$E3

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):

$$\text{L }=\frac{\text{q}}{{\displaystyle \underset{0}{\overset{\inf }{\int }}SA(t)dt}}$$E4

Recycling rate (R) is the part of the entry rate (R₀) which leaves the pool of interest and returns to it during the experiment. It is the difference between the entry rate (R₀) and the irreversible loss rate (L) (Nolan and Leng 1974):

$$\text{R }={\text{ R}}_0-\text{ L}$$E5

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):

${\text{dQ}}_{\text{i}}/\text{dt }=({\displaystyle \sum _{j=1}^{10}Qj}*{\text{ k}}_{\text{ij}})+\left({\text{Q}}_{\text{i}}*{\text{ k}}_{\text{ii}}\right)\text{ for i }=\text{ 1}\to \text{1}0$E6

Where: Q_j = Size of pool j (source of flux to pool i).

• k_ij = Rate constant for flow to i from j.

• Q_i = Size of pool i (pool of interest).

• k_ii = Sum of all rate constants for flows from compartment i:

$$\text{i}.\text{e}{.}_{\text{ji}}\text{for j }=\text{ 1}\to \text{1}0$$E7

If the coefficients of the differential equations (k_ij) are entered into a source/destination matrix (A), the matrix may be used to predict the magnitude of the flux between each compartment:

$$\text{<munder accentunder="true">F<mo>\_</mo></munder>}=\text{<munder accentunder="true">A<mo>\_</mo></munder>}*{\text{<munder accentunder="true">Q<mo>\_</mo></munder>}}_0$$E8

Where F is the vector of the sum of the fluxes to and from all pools, Q₀ is the vector of initial pool sizes, and A is the coefficient matrix derived above. The future state of the system then becomes:

$${\text{<munder accentunder="true">Q<mo>\_</mo></munder>}}_{\text{t}+\text{1}}=\text{<munder accentunder="true">F<mo>\_</mo></munder>}*{\text{<munder accentunder="true">Q<mo>\_</mo></munder>}}_0$$E9

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:

$${\text{<munder accentunder="true">Q<mo>\_</mo></munder>}}_{\text{t}+\text{1}}={\text{ e}}^{\text{<munder accentunder="true">A<mo>\_</mo></munder>}}{}^{\text{t}}*{\text{<munder accentunder="true">Q<mo>\_</mo></munder>}}_0$$E10

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^®.

Figure 1.

A 10-compartment flow diagram of metabolite flux in fasting northern elephant seals as simulated in this study. The boxes (or pools) represent the moles of carbon present in the animal as each metabolite, and the arrows represent the interchange of carbon between metabolites. There are no inputs to the model because the animal is fasting.

3. Results

Table 1 shows the initial conditions of the vector Q₀ which contains the initial conditions of the system, and the first vector of fluxes calculated by multiplying Q₀ by the source/destination matrix A, which is contained in Table 2. The characteristic equation of this matrix is: 0 = 1.0x¹⁰ – 0.186x⁹ + (9.34x10^-3)x⁸ + (9.50x10^-5)x⁷ + (3.78x10^-7)x⁶ + (6.45x10^-10 )x⁵ + (4.18x10^-13)x⁴ + (7.95x10^-17)x³ + (7.88x10^-22)x² + (2.15x10^-28)x + 0.

	Qi	FLUX
1. Adipose	1.2500e+002	-1.2750e-003
2. Glycerol	2.3100e-0 01	1.3800e-005
3. Palmitate	6.0000e-003	1.5000e-006
4. Lactate	7.5000e-003	3.9000e-005
5. Glucose	4.0000e-002	4.6722e-004
6. Alanine	4.4000e-003	-2.4304e-004
7. Urea	4.1700e-001	-1.2070e-004
8. Protein	2.0000e+002	-5.6160e-005
9. Ketones	3.5000e-003	8.1250e-006
10. Sink	1.0000e-000	9.3508e-004

Table 1.

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.

SOURCE POOL
DESTINATION
POOL	ADIPOSE	GLYCEROL	PALMITATE	LACTATE	GLUCOSE	ALANINE	UREA	PROTEIN	KETONES	SINK
ADIPOSE	-1.02e-5	0	0	0	0	0	0	0	0	0
GLYCEROL	5.10e-6	-2.70e-3	0	0	0	0	0	0	0	0
PALMITATE	5.10e-6	0	-1.06e-1	0	0	0	0	0	0	0
LACTATE	0	0	0	-1.20e-3	1.20e-3	0	0	0	0	0
GLUCOSE	0	2.70e-3	0	1.20e-3	-4.50e-3	3.30e-3	0	0	0	0
ALANINE	0	0	0	0	0	-6.80e-2	0	2.808e-7	0	0
UREA	0	0	0	0	0	1.00e-3	-3.00e-4	0	0	0
PROTEIN	0	0	0	0	0	0	0	-2.808e-7,	0	0
KETONES	0	0	3.25e-3	0	0	0	0	0	-3.25e-3	0
SINK	0	0	1.03e-1	0	3.30e-3	3.40e-2	3.00e-4	0	3.25e-3	-1

Table 2.

Source/destination matrix used to simulate the kinetics of metabolite flux in fasting northern elephant seal pups. Non-diagonal elements represent the rate constants for flow from the column pool to the row pool, in terms of time^-1. The diagonal elements are the sum all rate constants for flows from the pool, i.e. the sum of the column.

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 τ² 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 2.

Dimensionless Cartesian phase-plane representing the range of trace (τ) and determinant (Δ) values possible, and the regions of qualitative behavior of linear systems. The point corresponding to the system simulated here lies on the y-axis just to the left of the x-y intercept.

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 3.

Linear plot of 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 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.

Figure 4.

Log-normal plot of 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. The glycerol pool appeared to decline at a constant rate throughout the simulation.

Figure 5.

Linear plot of the changes in carbohydrate and carbohydrate-derived pool 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.

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.

Linear plot of the changes in the size of the various lipid and lipid-derived pools over a 24-hour simulation. The ketone pool rose slightly early on, and then declined slowly, as did the palmitate pool. The glycerol pool declined at a constant rate throughout the simulation.

Figure 8.

Changes in nitrogen containing pools over time. 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 9.

Changes in carbohydrate and carbohydrate-derived pool over 24 hours. Apparently the glucose and lactate pools were initialized incorrectly, as reflected in their rapid increase in the first hours 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 10.

Changes in the glucose, glycerol, lactate, alanine and ketones pools over 24 hours after re-estimation of the initial conditions by extrapolating the linear parts ofFigures 6, 7, and 8 to time zero.

4. Discussion

As initially formulated, the value of the rate constant for efflux from the sink pool (k_10,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 τ² 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, et. al.1980), and are validated by data which show that the lean body mass of the animal doesn’t change during significantly during the fasting period (Ortiz, et. al. 1978). Lack of significant protein catabolism, with concomitant low urea levels, indicates that the animals do not maintain their elevated blood glucose levels at the expense of gluconeogenesis from amino acids (Keith 1984). Close similarities in the rates of decline of the glucose, glycerol, and lactate pools in the later parts of the simulation may be indicative of the extensive glucose carbon recycling which occurs in these animals during fasting. There is extensive interchange of carbon between these three pools, as indicated by high levels of Cori cycle and glucose-glycerol cycle activity (Keith 1984, Keith and Ortiz 1989). It is postulated that this high degree of recycling may be one reason for the ability of these animals to avoid ketoacidosis, a major deleterious consequence of fasting in many other mammals, and thus allow them to undergo a prolonged fast during this vulnerable period in their life history.