We describe a computational model to simulate measurements from near-infrared spectroscopy (NIRS) and magnetic resonance spectroscopy (MRS) in the piglet brain. Piglets are often subjected to anoxic, hypoxic and ischaemic insults, as experimental models for human neonates. The model aims to help interpret measurements and increase understanding of physiological processes occurring during such insults. It is an extension of a previous model of circulation and mitochondrial metabolism. This was developed to predict NIRS measurements in the brains of healthy adults i.e. concentration changes of oxyhaemoglobin and deoxyhaemoglobin and redox state changes of cytochrome c oxidase (CCO). We altered and enhanced the model to apply to the anaesthetized piglet brain. It now includes metabolites measured by 31P-MRS, namely phosphocreatine, inorganic phosphate and adenosine triphosphate (ATP). It also includes simple descriptions of glycolysis, lactate dynamics and the tricarboxylic acid (TCA) cycle. The model is described, and its simulations compared with existing measurements from piglets during anoxia. The NIRS and MRS measurements are predicted well, although this requires a reduction in blood pressure autoregulation. Predictions of the cerebral metabolic rate of oxygen consumption (CMRO2) and lactate concentration, which were not measured, are given. Finally, the model is used to investigate hypotheses regarding changes in CCO redox state during anoxia.
In preclinical neonatology research, piglets are often used as experimental models for human neonates. At birth, their brains have a similar level of maturity to that of human brains . Many studies have used piglets to help understand the effects of oxygen deprivation on cerebral blood flow and metabolism [2,3]. Piglets have also been used to investigate encephalopathy following hypoxia–ischaemia (HI), a major cause of perinatal brain injury . One review found 23 per cent of 292 animal studies of perinatal hypoxic ischaemic encephalopathy had used piglets .
The physiological response of the brain during insults such as HI is complex. To help understand it, we have adopted a systems biology approach, and thus developed a mathematical model of circulation and metabolism in the neonatal piglet brain. The model has been created by extending and adapting previous physiological models of the brain [6,7]. Like these, it includes a biophysical sub-model of blood flow on a macroscopic scale, linked to a biochemical sub-model of cellular metabolism. We aim to use the model to help interpret experimental measurements, in particular, non-invasive measurements from near-infrared spectroscopy (NIRS) and magnetic resonance spectroscopy (MRS), both of which have been used extensively in piglet studies.
NIRS uses light to measure concentration changes of oxyhaemoglobin and deoxyhaemoglobin in the blood . Changes in the total haemoglobin concentration in the brain correspond to changes in cerebral blood volume. NIRS has a variety of clinical applications , and has frequently been used to study the newborn brain .
NIRS can also monitor redox changes at the CuA centre of cytochrome c oxidase (CCO). The difficulty in interpreting the NIRS CuA measurement was the main motivation for the development of a previous model by Banaji et al. , which we will refer to as BrainSignals. We have used this as the basis of our piglet model, since we also aim to investigate interpretation of the CuA signal. CCO is responsible for oxygen consumption in mitochondrial respiration; therefore, the redox state of CuA changes with mitochondrial oxygen tension. However, it is also affected by other factors such as adenosine triphosphate (ATP) and adenosine diphosphate (ADP) concentrations.
MRS uses the principle of nuclear magnetic resonance to measure the concentrations of a range of metabolites . Several different nuclei can be measured including 31P and 1H. 31P MRS measures concentrations of the phosphorus-containing metabolites ATP, inorganic phosphate (Pi) and phosphocreatine (PCr). ADP is not present at high enough concentrations to be measurable. Among the metabolites detectable by proton MRS are lactate and N-acetyl-aspartate.
In this paper, we (i) describe our new development of a piglet brain model, and its addition of simulating the MRS measurable concentrations of lactate, PCr, ATP and 31P; (ii) test the model by comparing its simulations of anoxia with experimental results from Springett et al.  involving simultaneous NIRS and 31P MRS measurements; (iii) present hypotheses regarding the interpretation of the reduction in CCO observed during anoxia in these experiments, and the subsequent hyperoxidation during the resuscitation phase; and (iv) predict CMRO2 and lactate concentrations, which were not measured, and propose a hypothesis regarding brain autoregulation in 1-day old piglets.
A schematic diagram of the model is shown in figure 1. It was based on BrainSignals , a model designed to simulate NIRS signals in the healthy adult brain. It was adapted to the neonatal piglet brain by changing the values of appropriate parameters. Of the 107 explicitly set numerical parameters in our model, 78 are taken from BrainSignals, and 11 of their values have been changed as shown in table 1. As well as these parameter changes, the model was expanded from BrainSignals to simulate variables that are measured by MRS. These include the 31P MRS measured variables ATP, PCr and Pi, and the lactate concentration measured by 1H MRS.
The approach for this model expansion followed the approach taken in BrainSignals i.e. the model was constructed to represent the underlying physiology. However, processes that are not commonly followed by non-invasive measurements necessarily require simplification to prevent the model from becoming overly complicated. Therefore, some processes and parameters are included which are chosen to fit relevant experimental data. The expansion required the addition of 29 numerical parameters shown in table 2. Where possible, these were set from experimental measurements, but some values were taken from other models. In addition, several rate constants were set by flux balance analysis, subject to the constraints of the explicitly set parameter values. A description of the piglet model is given below, concentrating on those aspects that differ from, or are additional to BrainSignals. Full details of the piglet model can be found on the model website .
The circulatory part of the model comprises three compartments: arteries and aterioles, capillaries and veins. The arterial and venous compartments have variable volumes Va and Vv, normalized such that the sum of their normal values (Va,n and Vv,n) is 1. For simplicity, the volume of the capillaries is ignored, as it is assumed to be small . The arterial/arteriolar compartment has a variable conductance G which is sensitive to four inputs: the capillary oxygen concentration [O2,c], the mean arterial blood pressure Pa, the arterial pressure of carbon dioxide PaCO2, and a parameter representing neuronal activation. The oxygen concentration affects conductance via the following equations. 2.1and 2.2
Here, is a time constant and represents a low-pass filtered version of [O2,c] with normal value . RO is a parameter controlling the magnitude of the response. The muscular tension in the arterial wall Tmax depends on η (which is the sum of the O2 term shown and three similar terms for PaCO2, arterial blood pressure and the demand parameter) in the following way: 2.3where Tmax0 is a constant, kaut is a parameter representing the level of autoregulation, and μ is a sigmoidal function of η. An average vessel radius is calculated from the balance of pressures and tensions in the vessel wall. This in turn determines the conductance of the arterial/arteriolar tree via Poiseuille's law. A decrease in oxygen concentration therefore leads to an increase in the vessel radius and consequently an increase in blood volume and blood flow.
The venous compartment has a fixed resistance Gv. The pressure at its boundary with the arteries Pv is calculated by equating the cerebral blood flow (CBF) through the arteries and veins giving 2.4where Pvs is the pressure in the venous sinuses. The normalized venous volume Vv is then calculated by 2.5where Cv is the venous compliance, including a normalizing factor.
All blood compartments have a fixed haemoglobin concentration [Hbtot]. In each compartment, a fraction of this haemoglobin is oxygenated. These fractions are determined from the arterial oxygen saturation (a model input) and from the rate of oxygen transport. Dissolved oxygen from the capillaries diffuses directly to the mitochondrial compartment. The simulated NIRS haemoglobin signals ΔHHb and ΔHbO2 are calculated from the concentrations of oxy and deoxyhaemoglobin in the arteries and veins as follows: 2.6and 2.7where [HHb] and [HbO2] are the concentrations of deoxy- and oxyhaemoglobin in the blood, and subscripts a and v refer to the arterial and venous compartments.
The metabolic part of BrainSignals is limited to oxidative phosphorylation. Particular attention is given to the CCO complex. This part of BrainSignals has been preserved, and extended to include ATP use and production, glycolysis and lactate dynamics. Many of the additional processes take place in the cytoplasm. However, with the exception of protons, no transport processes between cytoplasm and mitochondria are explicitly modelled. All other metabolites are considered in only one of these compartments.
An overview of the processes in the model involving the adenosine phosphates is shown in figure 2a. ATP is produced by ATP synthase in the mitochondria and by glycolysis in the cytoplasm. The rate of flow of protons through ATP synthase depends on the electrochemical potential Δp and the Gibbs free energy of ATP hydrolysis ΔG as follows: 2.8where LCV,max is the maximum rate of proton flow, rCV is a parameter controlling the relative forward and backward rates, and F is the Faraday constant. This is altered from BrainSignals to give a direct dependence on the energy state. The parameter na is the number of protons required to phosphorylate one ATP, and includes the proton used in the exchange of ADP and ATP across the mitochondrial membrane. The Gibbs free energy of ATP hydrolysis ΔG is calculated by 2.9where ΔG° is the standard Gibbs free energy of ATP hydrolysis, Z is a thermodynamic constant, and the ratio gp is known as the phosphorylation potential. All concentrations are taken to be cytoplasmic concentrations. The ATP synthesis rate is equal to LCV/na.
ATP use is modelled as a single Michaelis–Menten process with a half-maximum concentration (km,ATP) much lower than the normal ATP concentration. The reaction rate is therefore close to its maximum at normal ATP concentrations, and does not begin to decline significantly until ATP is low.
The final two processes influencing ATP concentration are buffering by phosphocreatine and AMP which occur with the following rates. 2.10and 2.11
Both are modelled as mass action reactions and are assumed to reach equilibrium quickly. The ratio of forward and backward rate constants kPCr and kPCr− is equal to the effective equilibrium constant for the reaction Keq,PCr*.
The remaining metabolic processes relate to glycolysis and lactate dynamics, and link these to the electron transport chain via the production of nicotinamide adenine dinucleotide (NAD) as shown in figure 2b. The net rate of glucose transport is 2.12where [gluc] is the glucose concentration in the cytoplasm and [glucc] is the glucose concentration in the capillary, which is fixed. The maximum rate vglut is set so that when glucose concentration is at its normal value [gluc]n, the net transport rate is equal to the normal metabolic rate of glucose consumption CMRgluc,n. The same form is used for lactate transport, with the net transport rate given by 2.13where [lac] is the lactate concentration in the cytoplasm, and [lacc] is the capillary lactate concentration which is fixed. Transport of lactate between brain cells has been the subject of much recent modelling, owing to the interest in the neuron–lactate aspartate shuttle [36–38]. Neurons and astrocytes have different expression of the MCT family members, leading to different kMCT . In our model however, there is no distinction between the two cell types, and we use an intermediate value for kMCT. The parameter vMCT is calculated from the normal cerebral metabolic rate of lactate consumption CMRlac,n. If the model is to be at a steady state under normal conditions, the stoichiometry of the reactions gives an expression for CMRlac,n
Satisfying this equation with the chosen values for CMRgluc and CMRO2 results in a net production of lactate at baseline conditions.
Glycolysis is modelled as a single-stage Michaelis–Menten process with a rate dependent on glucose, ADP and Pi concentration. The maximum rate vglyc is given by 2.14where vglyc,n is the normal maximum rate of glycolysis. This expression represents the regulation by ATP and AMP of phosphofructokinase, an enzyme that catalyses an important regulatory step in glycolysis . The parameter I controls how much the rate of glycolysis can change for given changes in ATP and AMP concentrations.
The interconversion between pyruvate and lactate is modelled as a mass action reaction with rate 2.15
This reaction involves a proton, and the rate constants are related by 2.16
The magnitude of the rate constants is set to give a ratio that adapts effectively instantly to pH changes. This reaction, and glycolysis, also involve NAD, but NAD in the cytoplasm is not modelled. However, the reactions cause a conversion between mitochondrial NAD and the reduced form of NAD (NADH) to give the correct stoichiometry.
The transformation of pyruvate to acetyl-CoA and the whole tricarboxylic acid (TCA) cycle are combined into a single reaction that consumes one molecule of pyruvate and produces five reducing equivalents. These are all assumed to be NADH for simplicity. The rate is given by 2.17
The km values were estimated from the detailed model of the TCA cycle by Wu et al. . The km for pyruvate is small, since this model suggests that TCA cycle rate is not sensitive to pyruvate concentration unless it falls very low. The rate of the TCA cycle is also influenced by ATP and ADP concentrations. However, the main controlling factor is thought to be the NAD/NADH ratio .
The electron transport chain is described by three reactions. The first is the transfer of four electrons from NADH to the CuA centre of CCO. The second is the transfer of electrons from here to the cytochrome a3 centre, and the final reaction is the reduction of oxygen to water. The rates are calculated based on thermodynamic principles. The rate of the third reaction, which is proportional to CMRO2, is given by 2.18where k3,0 is a constant, f a function of the proton motive force Δp, [O2] the mitochondrial oxygen concentration and [a3r] the concentration of reduced cytochrome a3.
Model development, simulations and analysis were carried out using the BRAINCIRC modelling environment . Steady-state simulations were used to analyse the model behaviour, and dynamic simulations were used to compare the model with experimental data.
The steady-state predictions by the model of CBF with changes in blood pressure, oxygen saturation and arterial carbon dioxide tension (PaCO2) are shown in figure 3. The changes in PCr and ATP concentration, CMRO2 and ΔCuA with arterial oxygen saturation are also shown. These simulations are compared with the equivalent simulations with the parameter values for adults shown in table 1.
3.1. Simulation of anoxia
The model was used to simulate a previous experiment involving brief anoxias in six newborn piglets . The piglets were anaesthetized with isoflurane and artificially ventilated. Their inspired oxygen was reduced from 40 per cent to 0 for 105 s. After returning to 40 per cent O2 for 10 min, the anoxia was repeated six times for each piglet. Mean arterial blood pressure (MAP) was monitored throughout, and continuous NIRS and 31P MRS measurements were recorded. The results showed a decrease in ΔHbO2 and increase in ΔHHb during anoxia. Immediately after anoxia there was an increase in total haemoglobin concentration indicating an increase in cerebral blood volume. There was a reduction of CuA during anoxia and a small hyperoxidation upon reoxygenation. The MRS results showed a decrease in PCr and inorganic phosphate concentration, but no changes in ATP concentration were seen.
The MAP and arterial oxygen saturation (SaO2) were used as inputs to the model, and its outputs were compared with the averaged NIRS and MRS measurements, as illustrated in figure 4.
Carbon dioxide levels were not reported during the challenge and assumed to remain constant. Oxygen saturation measurements were not available, so SaO2 was estimated from the experimental protocol. This estimate is shown, along with the measured MBP, in figure 5. Comparisons between modelled and measured results showed some significant differences. By changing four parameters that were known to affect relevant parts of the model, the simulations were improved. These changes are listed in table 3 and are explained as follows:
—the autoregulation capacity in the model was decreased by reducing the parameter kaut from 1.0 to 0.5. After this change, the maximum venous volume increase was approximately 45 per cent. The effect of this change on the autoregulation curve is shown in figure 6. The change in CBF as a percentage of baseline between 30 and 70 mmHg increased from 0.8% mmHg−1 with kaut = 1.0 to 2.1% mmHg−1 with kaut = 0.5. The improvement in the simulation of total haemoglobin concentration changes is shown in figure 7;
—the glycolysis rate was made more sensitive to changes in AMP and ATP concentration by increasing the parameter I from 3 to 50. This resulted in an approximately sevenfold increase in glycolysis rate during anoxia. Also, when I = 50, the PCr concentration at the end of anoxia was 55 per cent of its baseline value compared with 25% when I = 3;
—the normal ratio of PCr to Pi was increased from 1.5 to 2.7. In the model, if ATP concentration is constant, changes in PCr and Pi concentration are equal and opposite (ignoring AMP which is present only at very small concentrations). The experimental results therefore require a ratio of [PCr]n to [Pi]n of 2.7. With this ratio, the model will accurately simulate the Pi measurements providing the PCr measurements are accurately simulated;
—the normal NAD/NADH ratio was decreased from 9 to 1.5. This resulted in a decrease in the rate of CuA reduction.
The modelled and measured results after implementing these changes are shown in figure 8. The model's predictions for lactate concentration and CMRO2 are shown in figure 9. CMRO2 was predicted to decrease to 15% of its baseline value during anoxia and increase to 108% during reoxygenation.
In addition to these changes, the model was altered to explore the simulation of the CuA redox state. This was done to investigate the hyperoxidation seen following anoxia, which is seen in the data but is not well understood. Initially, equation (2.18) was changed to 3.1where k = 5 µM. The model was then further changed such that there was no dependence of ATP synthesis rate on phosphorylation potential by replacing ΔG with its normal value in equation (2.8). The results of these changes are shown in figure 10. When equation (3.1) is used, the simulation of ΔCuA is significantly improved as SaO2 begins to drop (although the simulation is slightly worse at lower SaO2). Also, the magnitude of hyperoxidation of CuA during recovery is decreased. When, in addition, the dependence of ATP synthesis rate on phosphorylation potential is removed, no hyperoxidation is seen.
The model has been applied successfully to simulate oxygenation and metabolic changes in the brains of newborn piglets during anoxia. It correctly predicted no significant changes in ATP concentration throughout the experiment. To simulate the other measured signals, changes to the original model were required. These changes offer an insight into the physiological interpretation of the experimental results as discussed below. In addition, the model allows predictions of variables that are more difficult to measure, including CMRO2 and lactate concentration.
During recovery from anoxia, a large increase in MAP accompanied by an increase in [HbO2] was measured. As there was not a corresponding reduction in [HHb], this implies an increase in cerebral blood volume. The model suggests that this hyperaemia cannot be attributed to changes in the arterial volume alone. A dilation of the cerebral arteries to this extent would be predicted to cause a several-fold increase in CBF, which, unless accompanied by a corresponding increase in CMRO2, would lead to a decrease in [HHb] via a washout effect. It is more likely, therefore, that the increased blood volume following anoxia has a significant venous contribution, which is driven by the increase in blood pressure. However, the autoregulation behaviour of the model causes changes in arterial pressure of this magnitude to be damped by changes in arterial resistance before they are felt at the veins. A reduction of the parameter kaut from 1.0 to 0.5 reduced this cancelling effect, and so led to a better replication of the observed hyperaemia. Previous studies in piglets have found that CBF does not change in the range of 50–80 mmHg [43–45]. However, these studies used piglets older than 24 h. In one study involving piglets of different ages, it was found that for piglets less than 4 days old, the average change in CBF as blood pressure was varied was 1.1%/mmHg, whereas for piglets older than 4 days old it was 0%/mmHg . The interpretation of the anoxia measurements using our model suggests that cerebral autoregulation was impaired or not fully developed in these newborn piglets.
The relationship between PCr and ATP concentration in the model as oxygen saturation falls shows the expected buffering of ATP by PCr, with ATP not decreasing until PCr is low. The model predicts that ATP concentration will begin to decrease when SaO2 falls below 40 per cent. However, the model predicted a drop in PCr concentration during anoxia larger than that observed. Increasing the parameter I allowed a greater increase in glycolysis rate, and therefore a greater rate of ATP synthesis during anoxia. This glycolysis rate could not be sustained for longer anoxias, because the model predicts that glucose is used faster than its maximum rate of transport into the cell. The value of I was increased from 3 to 50. Above this value, there is almost no further change in glycolysis rate, because changes in ATP and AMP concentrations become limiting. For large I, the glycolysis rate in the model increased approximately sevenfold. This is greater than the fivefold increase calculated to occur in foetal rats during ischaemia . It should be noted that the extent of PCr concentration decrease in the model is also sensitive to the normal PCr concentration [PCr]n. There are two features of particular interest in the CuA experimental results: the delay between the decrease in HbO2, and the reduction of CuA, and the hyperoxidation of CuA during recovery. The former is not reproduced by the model in its original form, but the latter is. It has been observed previously in piglets that CuA redox state is not affected by mild hypoxia  i.e. CCO has a low apparent km for oxygen. A similar effect has also been observed in adult rats [49,50], but not in adult humans [51,52]. The mechanisms for this effect are unknown, and it does not emerge from our alterations of the BrainSignals model to simulate the neonatal piglet brain. But, when included directly in the model, it improved the simulation of CuA reduction at the onset of anoxia. As expected, it also decreased the magnitude of CuA hyperoxidation following anoxia, since this is partly caused by the increased oxygen tension at the mitochondria as a consequence of the increased oxygen delivery. An oxidation is still seen however, because the increased concentrations of ADP and Pi during and following anoxia increase the rate of ATP synthesis via equations (2.8) and (2.9). The proton motive force across the mitochondrial membrane is decreased which increases CMRO2 and, with the parameter values used, results in an oxidation of CuA. If this dependency on phosphorylation potential is removed, no oxidation is seen. This supports the suggestion of Springett et al.  when reporting the experimental results, that the hyperoxidation was a result of an increased CMRO2. A CuA hyperoxidation following anoxia was also observed in an earlier study . The authors suggest that a pH drop could partly be responsible for the oxidation via a direct affect of pH on CuA. There is no mechanism for this effect in our model, although a small change in pH was observed following anoxia, as measured by MRS.
Another possible explanation for the oxidation is a decrease in the substrate supply to the TCA cycle. Our model suggests the opposite is true for pyruvate, whose concentration increases during anoxia, and takes several minutes to return to normal. However, there is an increase in NADH/NAD ratio which will tend to decrease the rate of the TCA cycle; but this appears to have only a small effect on the oxidation of CuA. The NAD/NADH ratio does however have a significant effect on the rate of reduction of CuA. Although changes in mitochondrial NADH concentration can be measured by fluorescence, the NAD/NADH ratio is difficult to measure in vivo . Estimates of this ratio differ, for example in rat liver cells a ratio of 5–10 was measured , but in myocytes isolated from newborn piglets a ratio of 1.2 was found . We found a decrease in the model's normal NAD/NADH ratio from 9 to 1.5 led to the simulated rate of CuA reduction matching the experimental rate. However, there are other factors that affect this rate, including the form of the TCA cycle rate dependence on NAD/NADH ratio, and the sensitivity of the reactions of the electron transport chain to changes in the proton gradient. This part of our model requires further investigation.
There are some limitations of our model. In general, we wish to keep the model as simple as possible, while still accurately replicating experimental results in a way that is representative of the underlying physiology. We believe that this allows a greater understanding of the model's behaviour. However, some aspects may have been oversimplified. In particular, it may be necessary to model NAD, NADH and pH in the cytoplasm. The ratio of NAD to NADH in the cytoplasm is much larger than in the mitochondria . Changes in these variables affect the lactate to pyruvate ratio, and the rate of glycolysis. Also, we have not distinguished between astrocytes and neurons. Their different functions are important when considering lactate dynamics. Finally, for the anoxia experiment simulated here, the assumption of constant arterial CO2 concentration is unlikely to be correct. CO2 changes have an important effect on blood flow, and may have affected the simulation of the NIRS haemoglobin signals.
We have used a physiological approach to construct our model. The main strength of this type of modelling is that it allows hypotheses to be generated and tested. However, it can lead to more complex models, and this makes it more difficult to analyse. Therefore, some components that have been simplified are treated phenomenologically. In this paper, we relied on knowledge of the model and experimentation to identify parameters with important effects. In the future, we will implement a more rigorous sensitivity analysis to analyse the effects of all parameters on particular model outputs. We are currently using the model to simulate MRS and NIRS measurements during HI in piglets. These measurements have complex relationships with one another, and we will use the model to help understand them in the context of our biological knowledge. We will also use it to simulate new treatments such as xenon-augmented hypothermia which is being investigated in our group .
The authors would like to thank The Wellcome Trust (088429/Z/09/Z) for the financial support of this work.
- Received November 8, 2011.
- Accepted January 4, 2012.
- This journal is © 2012 The Royal Society
This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.