Systems Biology Model of Cerebral Oxygen Delivery and Metabolism During Therapeutic Hypothermia : Application to the Piglet Model

Hypoxic ischaemic encephalopathy (HIE) is a significant cause of death and disability. Therapeutic hypothermia (TH) is the only available standard of treatment, but 45–55% of cases still result in death or neurodevelopmental disability following TH. This work has focussed on developing a new brain tissue physiology and biochemistry systems biology model that includes temperature effects, as well as a Bayesian framework for analysis of model parameter estimation. Through this, we can simulate the effects of temperature on brain tissue oxygen delivery and metabolism, as well as analyse clinical and experimental data to identify mechanisms to explain differing behaviour and outcome. Presented here is an application of the model to data from two piglets treated with TH following hypoxicischaemic injury showing different responses and outcome following treatment. We identify the main mechanism for this difference as the Q10 temperature coefficient for metabolic reactions, with the severely injured piglet having a median posterior value of 0.133 as opposed to the mild injury value of 5.48. This work demonstrates the use of systems biology models to investigate underlying mechanisms behind the varying response to hypothermic treatment.


Introduction
Hypoxic-ischaemic encephalopathy (HIE) is a significant cause of death and morbidity amongst neonates with around 700,000 deaths attributed to HIE alone annually [1]. Following neonatal HI injury, infants are treated with therapeutic hypothermia (TH) at a body temperature of 33.5 °C [2]. However, 45-55% of cases treated with hypothermia end with death or moderate to severe neurodevelopmental disability [1,3]. During hypothermia, a close neuromonitoring is in place combining clinical electroencephalography (EEG) [4] and broadband NIRS [2,5] as a research tool in the neonatal unit in University College London Hospital (UCLH). After completion of TH, infants undergo magnetic resonance imaging (MRI) and spectroscopy (MRS) [6]. The collected multimodal data has the potential to provide not only diagnostic and prognostic information but also insights on the mechanisms of the injury.
Our approach to analysis of this multimodal data has been multifaceted, with one key facet being the development and application of a physiology-informed "mathematical model" of the cerebral circulation under a systems biology approach, which is specially designed for the interpretation of bNIRS signals [7][8][9][10][11][12].
The first model BRAINCIRC was developed in 2005 [7]. This was later simplified and extended to include metabolism in 2008 via the BrainSignals model [8]. This was used to develop the BrainPiglet model [9], with the piglet being a common clinical model of human neonates. BrainPiglet was extended and used to investigate the effect of hypoxic ischaemia in the piglet model in the BrainPigletHI model [10]. Separate to the piglet models, BrainSignals was simplified further in the BrainSignals Revisited (BSRV) model to improve run time [11]. BSRV was extended to include extracerebral blood flow in BSX and looked at the confounding effects of the scalp on measurements [12]. Recently, significant work has been undertaken in developing the "BP Hypothermia" model (BPH1) which extended the BrainPiglet HI model to include temperature as an input [13]. This is needed to properly interpret data collected during therapeutic hypothermia. For example, it has been observed that both cerebral metabolic rate of oxygen (CMRO 2 ) and cerebral blood flow (CBF) in piglets decrease with reduced body temperature [14]. For reliable inferences to be made using systems biology models, they must be able to simulate this behaviour.
We present here an extended version of this model, BP Hypothermia 2 (BPH2), able to incorporate the effect of temperature separately for both metabolic and haemodynamic processes and reactions. This model is then validated against data collected from piglets undergoing TH following HIE.

Methods
BPH2 follows the same form as the original BP Hypothermia model [14], incorporating the effect of temperature by adapting work by Orlowski et al. [15]. Reaction rates, Michaelis-Menten rate constants and diffusion rates are modified by the quantities where Q 10, haemo is the Q 10 coefficient for haemodynamic reactions, Q 10, met is the Q 10 coefficient for metabolic reactions, k i, * is the reaction rate for the ith haemodynamic reaction and k j, * is the reaction rate for the jth metabolic reaction at temperatures T new and T previous . Q 10 is the temperature coefficient, defined as the ratio of reaction rates measured for the same reaction at two temperatures 10 °C apart. 0 < Q 10 < 1 indicates that decreasing temperature increases the reaction rate, whilst Q 10 > 1 indicates that decreasing temperature decreases the reaction rate. Figure 5.1 outlines the structure of this model. This model was then used to analyse data from two piglets, as shown in Fig. 5.2. Data was collected as per [16], with the piglets' common carotid arteries occluded for a period of around 25 minutes. Following HI, the piglets were treated with TH at 33.5 °C. HbO 2 , HHb and CCO data was collected via bNIRS, whilst the thalamic lactate/N-acetyl-aspartate (Lac/NAA) ratio was measured at 24 hours using proton MRS. Piglet LWP475 suffered a mild injury with a 10-minute CCO recovery fraction of 122% and a Lac/NAA ratio of 0.21 (Lac/ NAA ≥ 0.39 [6] is associated with a poor outcome). Piglet LWP479 suffered a severe injury with a 10-minute CCO recovery fraction of 69% and a 24-hour Lac/NAA ratio of 1.03. Both piglets received TH with piglet LWP475 responding typically, as per [13], with increased HbO 2 , decreased HHb and decreased CCO, whilst piglet LWP479 responded atypically with an "inverted" response in the bNIRS signals. Model analysis was performed using the BayesCMD framework [17].
Before performing model fitting, it is necessary to reduce the model down to manageable number of parameters. This is done using sensitivity analysis (SA) in order to identify parameters that control the majority of the behaviour we are attempting to simulate. We used the Morris method, summarising model output by the normalised root-mean-square error (NRMSE) between simulation and measured data. This was performed twice per model, once substituting the NRMSE value of failed runs with zero and once replacing it by 10,000,000. This was done to try and find a balance between including sensitive parameters that may produce failed runs and only selecting parameters that overwhelmingly cause failed runs. Failed runs should not be ignored completely as the reason a run fails is likely to be due to an invalid value for a highly sensitive parameter. Parameter count was reduced from an original count of 234 to a fitted count of 11, with this parameters shown in Table 5.1. Temperature is added as an input to the model. Parameter Q 10, met and temporary variable Q temp, met are added to the metabolic compartment, and parameter Q 10, haemo and temporary variable Q temp, haemo are added to the blood flow compartment

Result
Following this, posterior parameter distribution fitting was performed using approximate Bayesian computation via a rejection algorithm with a total of 50,000,000 samples. Each parameter was given an uninformative uniform prior distribution. The top 0.01% of parameter samples, based on NRMSE, were used to generate posterior distributions for each piglet giving a posterior of 5000 samples. These are shown in Fig. 5.3 with piglet LWP475 shown in blue and piglet LWP479 in orange. A distinct difference between the parameter spaces of the two piglets is clearly visible for Q 10, met , Q 10, haemo , normal oxidised fraction of Cu A (a frac, n ) and the normal total concentration of HbO 2 binding sites in blood (X tot, n ), with Q 10, met showing the most obvious and potentially important disparity. For further information about a frac, n and X tot, n , see Table 5.1 and [10]. Figure 5.4 shows the posterior predictive distributions for each piglet generated by sampling repeatedly from the posterior distributions shown in Fig. 5.3. The posterior predictive distributions are shown in blue and include a 95% confidence interval, whilst data is shown in green. The model simulations show good agreement with the measured data with only a small time lag between the model and data for the CCO signal in piglet LWP479. Importantly, the model can correctly reproduce the overall behaviour of the system in both piglets suggesting that the posterior distri-

Discussion
We have successfully expanded our model of hypothermia in the piglet brain to incorporate different temperature effects on metabolic and haemodynamic reactions. Additionally, this has been achieved without introducing excessive complexity into the model, adding only one new parameter as compared to the previous BPH1 model -Q 10 is split into Q 10,met and Q 10,haemo . This model has then been validated against data collected from two piglets suffering the same injury to differing levels of severity. Not only was this able to validate the new model but it was also able to provide some insight into the different haemodynamic and metabolic behaviour seen in both piglets. between the two piglets is clearly visible, with piglet LWP479 having a distribution of values well below 1. For more information about specific parameters, please see Table 5.1 Piglet LWP475 suffered a mild injury, as indicated by both its typical response to hypothermia, as per [13], and its Lac/NAA ratio of 0.21. The Q 10,met and Q 10,haemo values obtained by fitting are all above 1 indicating that decreasing temperature does decrease reaction rates for both haemodynamic and metabolic reactions. In contrast, piglet LWP479 suffered a severe injury and exhibited an atypical response to hypothermia. The HbO2, HHb and CCO signals all produce behaviour opposite to that seen in piglet LWP475. Model analysis then showed that this is likely to be due to a breakdown in how the metabolic reactions responded to hypothermia, with Q 10,met values all within a narrow distribution well below 1, with a median value of 0.133. Thus, for this piglet, hypothermia will increase metabolic reaction rate rather than decrease it, but haemodynamic reaction rate will react "typically", reducing with temperature. A small time lag is seen between measured data and predictive posterior distribution, which may be due to the model reduction removing a parameter nec- and data (green) in both piglets, with the model able to reproduce overall large-scale behaviours Whilst it is not possible to draw clinical conclusions from these results yet, the different parameter spaces for these two piglets does potentially highlight issues with treating all HIE injuries in the same way. Findings here suggest that in certain circumstances the initial injury may impact on how the various reactions and processes respond to cooling and that this change in response may be the opposite of that desired. Further work will look at applying these methods and models to understanding data collected from human neonates that have suffered HIE.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.