Computational modelling of the piglet brain to simulate near-infrared spectroscopy and magnetic resonance spectroscopy data collected during oxygen deprivation

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.


INTRODUCTION
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 [1]. 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 [4]. One review found 23 per cent of 292 animal studies of perinatal hypoxic ischaemic encephalopathy had used piglets [5].
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 [8]. Changes in the total haemoglobin concentration in the brain correspond to changes in cerebral blood volume. NIRS has a variety of clinical applications [9], and has frequently been used to study the newborn brain [10].
NIRS can also monitor redox changes at the Cu A centre of cytochrome c oxidase (CCO). The difficulty in interpreting the NIRS Cu A measurement was the main motivation for the development of a previous model by Banaji et al. [7], 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 Cu A signal. CCO is responsible for oxygen consumption in mitochondrial respiration; therefore, the redox state of Cu A 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 [11]. Several different nuclei can be measured including 31 P and 1 H. 31 P MRS measures concentrations of the phosphorus-containing metabolites ATP, inorganic phosphate (P i ) 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 31 P; (ii) test the model by comparing its simulations of anoxia with experimental results from Springett et al. [12] involving simultaneous NIRS and 31 P 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 CMRO 2 and lactate concentrations, which were not measured, and propose a hypothesis regarding brain autoregulation in 1-day old piglets.

METHODS
A schematic diagram of the model is shown in figure 1. It was based on BrainSignals [7], 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 31 P MRS measured variables ATP, PCr and P i , and the lactate concentration measured by 1 H 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 [34].

Circulation
The circulatory part of the model comprises three compartments: arteries and aterioles, capillaries and veins. The arterial and venous compartments have variable volumes V a and V v , normalized such that the sum of their normal values (V a,n and V v,n ) is 1. For simplicity, the volume of the capillaries is ignored, as it is assumed to be small [35]. R O is a parameter controlling the magnitude of the response. The muscular tension in the arterial wall T max depends on h (which is the sum of the O 2 term shown and three similar terms for P a CO 2 , arterial blood pressure and the demand parameter) in the following way: where T max0 is a constant, k aut is a parameter representing the level of autoregulation, and m is a sigmoidal function of h. 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 G v . The pressure at its boundary with the arteries P v is calculated by equating the cerebral blood flow (CBF) through the arteries and veins giving where P vs is the pressure in the venous sinuses. The normalized venous volume V v is then calculated by where C v is the venous compliance, including a normalizing factor. All blood compartments have a fixed haemoglobin concentration [Hb tot ]. 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 DHHb and DHbO 2 are calculated from the concentrations of oxy and deoxyhaemoglobin in the arteries and veins as follows: À ½HHb a n V a;n À ½HHb v n V v;n Þ ð2:6Þ Computational modelling of the piglet brain T. Moroz et al. 1501 and where [HHb] and [HbO 2 ] are the concentrations of deoxy-and oxyhaemoglobin in the blood, and subscripts a and v refer to the arterial and venous compartments.

Metabolism
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 Dp and the Gibbs free energy of ATP hydrolysis DG as follows: where L CV,max is the maximum rate of proton flow, r CV 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 n a 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 DG is calculated by where DG8 is the standard Gibbs free energy of ATP hydrolysis, Z is a thermodynamic constant, and the ratio g p is known as the phosphorylation potential. All concentrations are taken to be cytoplasmic concentrations. The ATP synthesis rate is equal to L CV /n a . ATP use is modelled as a single Michaelis -Menten process with a half-maximum concentration (k m,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.
k PCr ½PCr½ADP À k À PCr ½ATP½Cr ð 2:10Þ and k AK ½ADP 2 Àk À AK ½ATP½AMP: ð2:11Þ Both are modelled as mass action reactions and are assumed to reach equilibrium quickly. The ratio of forward and backward rate constants k PCr and k PCr 2 is equal to the effective equilibrium constant for the reaction K eq,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 v glut ½gluc c ½gluc c þ k glut À ½gluc ½gluc þ k glut where v glyc,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 [40]. 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 k pl ½Py À k À pl ½lac: ð2:15Þ This reaction involves a proton, and the rate constants are related by 2CMR lac;n þ k À pl ½lac n ½Py n : ð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 v TCA ½Py½NAD ðk m;tcaP þ ½PyÞðk m;tcaN þ ½NADÞ : ð2:17Þ The k m values were estimated from the detailed model of the TCA cycle by Wu et al. [32]. The k m 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 [41]. The electron transport chain is described by three reactions. The first is the transfer of four electrons from NADH to the Cu A centre of CCO. The second is the transfer of electrons from here to the cytochrome a 3 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 CMRO 2 , is given by

Implementation
Model development, simulations and analysis were carried out using the BRAINCIRC modelling environment [42]. Steady-state simulations were used to analyse the model behaviour, and dynamic simulations were used to compare the model with experimental data.

RESULTS
The steady-state predictions by the model of CBF with changes in blood pressure, oxygen saturation and arterial carbon dioxide tension (P a CO 2 ) are shown in figure 3. The changes in PCr and ATP concentration, CMRO 2 and DCu A 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.

Simulation of anoxia
The model was used to simulate a previous experiment involving brief anoxias in six newborn piglets [12]. 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 O 2 for 10 min, the anoxia was repeated six times for each piglet. Mean arterial blood pressure (MAP) was monitored throughout, and continuous NIRS and 31 P MRS measurements were recorded. The results showed a decrease in DHbO 2 and increase in DHHb 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 Cu A during anoxia and a small Computational modelling of the piglet brain T. Moroz et al. 1503 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 (SaO 2 ) 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 SaO 2 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 k aut 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 21 with k aut ¼ 1.0 to 2.1% mmHg 21 with k aut ¼ 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 P i was increased from 1.5 to 2.7. In the model, if ATP concentration is constant, changes in PCr and P i 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 [P i ] n of 2.7. With this ratio, the model will accurately simulate the P i 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 Cu A reduction.
The modelled and measured results after implementing these changes are shown in figure 8. The model's   predictions for lactate concentration and CMRO 2 are shown in figure 9. CMRO 2 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 Cu A 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 ðk þ ½O 2 n Þ½a3rf ðDpÞ; ð3:1Þ where k ¼ 5 mM. The model was then further changed such that there was no dependence of ATP synthesis rate on phosphorylation potential by replacing DG 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 DCu A is significantly improved as SaO 2 begins to drop (although the simulation is slightly worse at lower SaO 2 ). Also, the magnitude of hyperoxidation of Cu A during recovery is decreased. When, in addition, the dependence of ATP synthesis rate on phosphorylation potential is removed, no hyperoxidation is seen.

DISCUSSION
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 CMRO 2 and lactate concentration. During recovery from anoxia, a large increase in MAP accompanied by an increase in [HbO 2 ] 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 CMRO 2 , would lead to a decrease in [HHb] via a washout effect. It is more likely, therefore, that  [12]. Blood pressure was taken from the measurements from one piglet. SaO 2 was estimated from the experimental protocol.   Computational modelling of the piglet brain T. Moroz et al. 1505 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 k aut 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][44][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 [46]. 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 SaO 2 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 [47]. 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 Cu A experimental results: the delay between the decrease in HbO 2 , and the reduction of Cu A , and the hyperoxidation of Cu A 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 Cu A redox state is not affected by mild hypoxia [48] i.e. CCO has a low apparent k m 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 Cu A reduction at the onset of anoxia. As expected, it also decreased the magnitude of Cu A 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 P i 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 CMRO 2 and, with the parameter values used, results in an oxidation of Cu A . If this dependency on phosphorylation potential is removed, no oxidation is seen. This supports the suggestion of Springett et al. [12] when reporting the experimental results, that the hyperoxidation was a result of an increased CMRO 2 . A Cu A hyperoxidation following anoxia was also observed in an earlier study [53]. The authors suggest that a pH drop could partly be responsible for the oxidation via a direct affect of pH on Cu A . 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 Cu A . The NAD/NADH ratio does however have a significant effect on the rate of reduction of Cu A . Although changes in mitochondrial NADH concentration can be measured by fluorescence, the NAD/ NADH ratio is difficult to measure in vivo [54]. Estimates of this ratio differ, for example in rat liver cells a ratio of 5 -10 was measured [55], but in myocytes isolated from newborn piglets a ratio of 1.2 was found [56]. We found a decrease in the model's normal NAD/ NADH ratio from 9 to 1.5 led to the simulated rate of Cu A 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 [55]. 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 Computational modelling of the piglet brain T. Moroz et al. 1507 lactate dynamics. Finally, for the anoxia experiment simulated here, the assumption of constant arterial CO 2 concentration is unlikely to be correct. CO 2 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 [57].