The dynamic interaction of systemic inflammation and the hypothalamic-pituitary-adrenal (HPA) axis during and after major surgery

Major surgery and critical illness produce a potentially life threatening systemic inflammatory response. The hypothalamic-pituitary-adrenal (HPA) axis is one of the key physiological systems that counterbalances this systemic inflammation through changes in adrenocorticotrophic hormone (ACTH) and cortisol. These hormones normally exhibit highly correlated ultradian pulsatility with an amplitude modulated by circadian processes. However, these dynamics are disrupted by major surgery and critical illness. In this work, we characterise the inflammatory, ACTH and cortisol responses of patients undergoing cardiac surgery and show that the HPA axis response can be classified into one of three phenotypes: single-pulse, two-pulses and multiple-pulses dynamics. We develop a mathematical model of cortisol secretion and metabolism that predicts the physiological mechanisms responsible for these different phenotypes. We show that the effects of inflammatory mediators are important only in the single-pulse pattern in which normal pulsatility is lost – suggesting that this phenotype could be indicative of the greatest inflammatory response. Investigating whether and how these phenotypes are correlated with clinical outcomes will be critical to patient prognosis and designing interventions to improve recovery.


INTRODUCTION
Major surgery and critical illness elicit a systemic inflammatory response of both clinical and scientific relevance 1 , which when uncontrolled leads to major morbidity and/or death 2,3 .
One of the key physiological systems that regulates the inflammatory response in humans is the hypothalamic-pituitary-adrenal (HPA) axis, which controls the rhythmic secretion of adrenocorticotrophic hormone (ACTH) and cortisol. The dynamics of these hormones are transiently disrupted by stressors, with HPA axis activity following major surgery 4  The secretion and effects of these hormones are also non-linear, making the overall effect at the level of an organism difficult to elucidate. Understanding how ACTH and cortisol are controlled during acute systemic inflammation (e.g. major surgery and critical illness) is vitally important because both excess and deficiency of these hormones are associated with death and major morbidity for patients 7 .
Under normal physiological conditions, ACTH and cortisol are in a state of continuous dynamic equilibration 8 and exhibit highly-correlated ultradian pulsatility with an amplitude modulated by upstream circadian processes (Fig. 1A). ACTH and cortisol dynamics are disrupted after major surgery: mean concentrations of cortisol are higher, mean concentrations of ACTH are lower, and the correlation between ACTH and cortisol is altered 4 . A combination of animal physiology experiments and mathematical modelling has also shown that acute inflammation can disrupt the high correlation between ACTH and cortisol, a phenomenon called dynamic dissociation 9 . However, the interaction between the HPA axis and the inflammatory cascade in humans is unknown (Fig. 1B), in part because the lack of high-resolution time-course data has precluded accurate modelling.
The most likely mediators of the interactions between inflammation and the HPA axis are cytokines such as IL1α, IL6 and TNFα, and the glucocorticoid hormone cortisol. The isolated effects of cytokines on the HPA axis at the cellular level are well known (e.g. IL6 increases cortisol synthesis 10,11 ). However, this information is not informative at the systems level due to the bi-directional feedback that all components exert on each other, in addition to the dynamic regulation of cortisol in target tissues (Fig. 1C). To address the question of how the inflammatory and HPA axis response interact at a systems level, we performed highfrequency serial blood sampling during and after coronary artery bypass grafting (CABG) surgery to generate profiles of ACTH, cortisol and inflammatory mediators IL6, IL8, IL10 and TNFα (Fig. 1D). These profiles were statistically classified according to dynamic phenotype and used to calibrate a mathematical model to characterise the underlying changes in HPA physiology.

Patient Characteristics
Ten patients underwent isolated CABG surgery via median sternotomy (See Table 1). CABG can take place with or without cardiopulmonary bypass (CPB). Five cases had their surgery performed with the use of CPB and five did not. When CPB was used, the duration was 67 ± 20.7 minutes. Patients and controls were not matched for age, height and weight, but were broadly similar in these parameters. No patient suffered major complications and this was indicated in the median critical care and hospital length of stay, which were shorter than the UK median (3.2 and 7.2 days respectively) 12,13 .  1D).
In the healthy controls, ACTH and cortisol exhibited regular ultradian pulsatility (T u = 2.5-3 hrs) with a circadian modulated amplitude (Fig. S1). S1 and S3), and supports previous work that demonstrates a dynamic dissociation between these hormones in this context 4,9 . The time-dependent Pearson correlation coefficient showed longer-lasting correlated and anti-correlated events between ACTH and cortisol compared to controls, while the angle and IPS showed a loss of phase coherence between ACTH and cortisol rhythms (Figs. S1 and S3). Taken together, these results suggest that the dynamic dissociation between these hormones could be explained by a misalignment of secretion between the pituitary and adrenal glands (due to signals from outside the HPA axis such as inflammatory mediators), an increase in cortisol half-life (due to reduced inactivation into cortisone), or a combination of both.

Identifying three distinct phenotypes of HPA axis activity following CABG
We investigated the dynamic dissociation between ACTH and cortisol levels by means of cross-correlation algorithms that quantify the synchrony between time series data 14,15 .
Specifically, we calculated: (1) the time-lagged cross correlation (TLCC) to estimate the peak association between ACTH and cortisol; and (2) the rolling window time-lagged cross correlation (RWTLCC) for a time-dependent quantification of the lag between ACTH and cortisol. These algorithms differ from traditional cross-correlations in that they consider the non-stationary nature of variables (see Methods), and allow us to characterise the interindividual variability of HPA axis responses in patients.
The TLCC and RWTLCC showed that the healthy control group had a high peak association between ACTH and cortisol, with ACTH leading cortisol by a mean lag of µ lag = 10 min which remained stable across all epochs ( Fig. 2A and Fig. S1). In contrast, we found that the patients experienced different lags and degrees of phase synchrony loss between ACTH and cortisol. Combining the period, TLCC and RWTLCC allowed us to group the dynamic responses of the HPA axis into three categories (Figs. 2B-D) according to the type of dynamic dissociation between ACTH and cortisol elicited by CABG as follows: -Two-pulse group: These patients showed two pulses of ACTH and cortisol with longer than physiological periodicity (T u = 5-6 hrs) but preserving a near-physiological peak association (µ lag = 13 min) and a stable phase synchrony across the entire 12 hr long sampling (Fig. 2B).
-Multiple-pulse group: These patients showed multiple pulses of ACTH and cortisol (T u = 2 hrs) with strong peak dissociation (µ lag = 106 min) and partial phase synchrony during a portion of the 12-hr long sampling (Fig. 2C).
-Single-pulse group: These patients showed a single large excursion of ACTH and cortisol with strong peak dissociation (µ lag = 77 min) and unstable phase synchrony across the entire 12-hr long sampling (Fig. 2D).
These techniques allow us to describe the inter-individual variability of dynamic responses to CABG in quantitative terms. They suggest that patients may fall into different categories depending on the type and degree of the dynamic dissociation between ACTH and cortisol following surgery. Namely, the two-pulse group experiences lesser dynamic dissociation, the

A model of cortisol concentration accounts for the dynamic changes of the HPA axis following CABG
The model assumes a delayed (10 minutes) ACTH input to a hypothetical adrenal gland underpinning cortisol production (implicitly including peripheral conversion from cortisone), cortisol turnover, and the adrenal sensitivity to ACTH stimulation. It also assumes a twocompartment, open-loop architecture where each compartment accounts for fast and slow processes contributing to cortisol dynamics, respectively 16 . In the model ( Fig. 3A and Methods), parameters p f and p s account for fast and slow cortisol production respectively (including adrenal maximum secretory capacity and peripheral conversion from cortisone); 1/l f and 1/l s are the fast and slow cortisol turnover rates, respectively; K A is the adrenal sensitivity to ACTH stimulation; and m is the Hill coefficient denoting the steepness of the sigmoidal function used to represent the adrenal response to ACTH.
We performed computer simulations of the model and calibrated it to physiological data. To do this, we used the ACTH data as an input and assessed the ability of the model to accurately capture the cortisol output trajectories over 1 million sets of parameter values.
The motivation for this was two-fold. First, we aimed to determine whether parameters that fit the control cortisol trajectories could also explain the cortisol trajectories in patients undergoing CABG without the cytokine input. If this were the case, then we could assume that any potential influence of inflammatory cytokines would occur at the level of the pituitary rather than the adrenal glands ( Fig. 1B). Second, if the parameters from the controls were not a good fit to the data from patients having CABG, then we could assess the changes required to achieve a good fit. This would help identify potential targets where cytokines could be influencing cortisol activity, for example by regulating secretion from the adrenal glands or activity in peripheral tissues (Fig. 1C).

Modelling controls and CABG patients
To calibrate the model, over 1 million sets of candidate parameter values were randomly selected within a latin hypercube and tested to fit the control data. The cost function (see Methods) used during the fitting process was evaluated for each of the controls and each of the CABG patients. Table S1 includes the ranges of tested values for each model parameter and the distribution of optimal parameter values that fit the control group. In order to establish a baseline set of parameters, we first identified the best fit parameter set (single lowest cost function evaluation) as well as the group of parameter sets whose cost function evaluation was within 10% of the lowest evaluation. We found that most of the fits to controls corresponded to scenarios where the Hill coefficient was fixed to m = 2; therefore, we used this value for the rest of the analysis. We also found that the error was low for all controls (median error = 0.07 across the fit parameter values) (Fig. 3B), showing that the model accurately captures the control data. Having many parameter sets that fit the control data with low error implies that there are many parameter sets that can fit all controls simultaneously. The dark red line in Fig. 3B shows the cortisol trajectories predicted by the model when using the median values of the distribution of parameters that best fit all controls ( = 0.068).
The median error for the fits to the CABG patient data was = 0.12, nearly double that of the fits to the controls. Fig. 3C shows the predicted cortisol trajectories using these fits against the CABG patient data. This highlights that while these trajectories were not as good as the fits to the controls, they still follow the cortisol trend observed in CABG patients. The predicted trajectories using the fits to the patient data perform better than when using the median parameter values fitted to the controls (Table S1)  better characterise the differences between the control and patient groups, we used violin plots to represent the distributions of optimal parameter values resulting from fitting the model to data (Fig. 3D). The cortisol production rates p f and p s , as well as the fast cortisol half-life l f , had similar statistical moments and shape of their distributions. In contrast, while the slow cortisol half-life l s was found to have similar statistical moments, its distribution was different between controls and patients, leaning towards higher values (implying slower cortisol turnover) in CABG patients. A more striking difference was observed in the adrenal sensitivity K A , where the estimated distribution had a lower, narrower distribution in the CABG case compared to controls. Taken together, this suggests that the same model explaining the cortisol dynamics observed in controls can also explain that in patients having CABG by assuming a slower cortisol turnover and an increased adrenal sensitivity to ACTH stimulation.

Modelling CABG patients with different phenotypes of cortisol profile
To investigate the origin of the distinct dynamic dissociation between ACTH and cortisol observed across the three groups of patients (two-pulses, multiple-pulses, and single-pulse), we fitted the model to identify optimal parameter values for each group (Table S2). The fits were marginally improved compared to the fits to all CABG patients pooled together, with the error decreasing to = 0.10 for the two-pulses group, = 0.08 for the multiple-pulses group, and = 0.10 for the single-pulse group ( Table 2). The shapes of distributions of parameter values associated with fast processes (p f , l f ) were also very similar across all CABG groups and between these and the control group (Fig. 3E). For parameters associated with slow processes, the cortisol production rate p s and half-life l s were predicted to remain within the same dynamic range (although slightly higher) as controls for the two-pulses and multiple-pulses CABG groups (Table S2). The shape of the l s distribution tending towards higher values in these patient groups. In contrast, the single-pulse CABG group was predicted to have a higher cortisol production rate p s , but a lower half-life l s . The adrenal sensitivity K A was similar between the control and the two-pulses CABG group, but was predicted to have lower values (implying higher adrenal sensitivity) for the single-pulse and multiple-pulses groups (Fig. 3E). Taken together, these results suggest that the degree of cortisol pattern change observed across the different CABG groups could originate from a combination of disrupted physiological processes. Specifically, the estimated distributions of parameter values for the two-pulses group maintain similar statistical moments to the healthy controls, only changing the shape of the l s distribution. The distributions for the multiple-pulses group are similar to the two-pulses group, except that the model also predicts an increased adrenal sensitivity (lower K A ). Finally, the distributions for the singlepulse group show that not only is the adrenal sensitivity increased, but also the cortisol secretion p s and its turnover rate 1/l s . The dynamic range and shape of the distribution of K A predicted for the two-pulses CABG group was very close to controls, suggesting that the adrenal sensitivity to ACTH is the same for these patients.

Model predictions under fixed adrenal sensitivity
To further investigate the origin of the dynamic dissociation between ACTH and cortisol, we used the assumption that CABG surgery does not markedly change the adrenal sensitivity to ACTH stimuli. This implies that such changes are unlikely to take place at the time scale of hours after surgery but may result from long term critical illness 17,18 . To do this, we fixed the adrenal sensitivity parameter to the median value estimated for controls (K A = 50.28) and identified the distributions of parameter values as before (Table S3). The model predicted low error trajectories for both the controls ( = 0.07) and the CABG group ( = 0.14). This marginally changed with respect to the baseline parameter fitting (compare Figs. 4A and 4B with Figs. 3B and 3C). The optimal parameter sets, when compared between the control and CABG group, showed these to be strikingly similar. The exception was the shape of the distribution for the half-life parameter l s , which tended towards higher values (Fig. 4C). This suggests that, if adrenal sensitivity does not change after CABG, the key physiological process underpinning dynamic dissociation in patients would likely be a slower cortisol turnover rate 19 .
Next, we investigated the model predictions for each CABG group under the assumption of fixed adrenal sensitivity. As before, the shapes of distributions for optimal parameter values associated with fast processes (p f , l f ) were very similar across all CABG groups and the control group (Fig. 4D). For parameters associated with slow processes, the cortisol production rate p s and half-life l s were predicted to remain within the same dynamic range (albeit slightly higher) as controls for the two-pulses and multiple-pulses CABG groups (Table S3), with the shape of the l s distribution leaning towards higher values in these patient groups. Consistent with the variable K A scenario, the single-pulse CABG group was again predicted to have a higher cortisol production rate p s , but a shorter half-life l s .  0.07 0.14 ---

Inflammatory mediators may underpin the dynamic changes in the interaction of ACTH and cortisol following CABG
While the model can predict cortisol dynamics driven by ACTH inputs in both healthy controls and CABG patients, minor discrepancies between the model predictions and data persist. We now explore whether these discrepancies may be related to the effect of inflammatory mediators, which have not been explicitly included in the mathematical model to this point. To do this, we first performed a principal component analysis (PCA) on the trajectories of inflammatory mediators following CABG (Fig. S2). The PCA identified IL6 and TNFα as the cytokines containing the most information about the inflammatory response ( Fig. S4). IL6 and TNFα have previously been shown to be involved in the regulation of 11b-HSD1/2 enzymatic cortisol to cortisone inter-conversion 20-22 (Fig. 1C). We

DISCUSSION
A maladaptive response to acute systemic inflammation underlies most of the morbidity and mortality from major surgery 23 and critical illness (most notably sepsis 24 and more recently COVID-19 25 ). Therefore, understanding the mechanisms that modulate this is key to designing clinical diagnostics and interventions that reduce the risk of patient death and morbidity. The HPA axis is one such system that is thought to overall counterbalance inflammation -although this has been difficult to elucidate due to the multiple levels of interaction and sites of feedback between the two systems. High-frequency blood sampling techniques and mathematical models have advanced our understanding of the interaction between the two systems -but previously not in a real-world environment 9,26 . Using a combination of high-frequency blood sampling and statistical analysis, this study has shown that at least three distinct phenotypes of cortisol dynamic profiles emerge after cardiac surgery; single-pulse, two-pulses and multiple-pulses. The mechanisms underpinning these three different phenotypes were inferred using mathematical modelling, with parameter estimations used to predict the contribution of distinct cortisol control mechanisms in patients. Lastly, we identified that inflammatory mediators may only be important in specific patterns of HPA axis activity, paving the way to identify those most at risk of serious harm from inflammation.
We postulated a two-compartment mathematical model of cortisol activity that discriminates between fast and slow processes, with model parameters corresponding to the adrenal sensitivity to ACTH, cortisol production (a contribution of adrenal secretory capacity and peripheral conversion from cortisone), and cortisol turnover and metabolic rates (arising from degradation and inactivation into cortisone). We determined the models with the best fits by assuming a similar adrenal sensitivity between patients undergoing CABG and controls. During fitting, only the parameters associated with fast processes remained invariant between the control and patient groups. This suggests that CABG may not particularly affect cortisol regulatory processes associated with short time scales (e.g. rapid cortisol synthesis and distribution out of the central (plasma) compartment), but only affects slower processes such as enzymatic degradation. When examining the slow parameters, the model predicted that patients in the single-pulse group had a greater slow cortisol production rate and higher slow turnover rate compared to controls and other patient groups. The single-pulse group also had the highest dynamic dissociation between ACTH and cortisol (strong peak dissociation and unstable phase synchrony). The model predicted that patients in the two-pulses group had a similar cortisol production rate and an increased slow turnover rate compared to controls. This group also had the lowest dynamic dissociation between ACTH and cortisol following surgery. Lastly, the model predicted that patients in the multiple-pulses group had an increased slow cortisol production rate and an increased slow turnover rate when compared to controls. This group also had a dynamic dissociation between ACTH and cortisol falling in the middle of the other two phenotypes.
Taken together, these results suggest which mechanisms may underpin the three different phenotypes of cortisol activity observed in patients undergoing CABG. Our model also suggests that inflammatory mediators do not play a significant role in regulating circulating levels of cortisol in all patients having CABG within the 12 hrs following surgery, as previously postulated 19 . However, in some cases the effect of inflammatory mediators may cause significant disruption to HPA axis dynamics, leading to a large sustained pulse of cortisol.
There are several limitations to consider when interpreting the results of the mathematical model. While the fast and slow time scale separation in the model has been described previously 16,26 , the optimisation is agnostic with respect to the physical mechanisms underlying this separation and only attempts to minimise the cost function. Therefore, some of the fast mechanisms are captured by the slow variable of the model and vice-versa, but this is likely to be minimal. On the other hand, we consider ranges of parameter sets that fit the data well, rather than only considering the single best fit parameter set. This is to minimise concerns of overfitting the data. This means that we can only assess the difference between groups, but not assume that a single group in isolation provides specific estimates of parameter values. Our sample size of patients and healthy participants was small, so we cannot be sure that we have captured the full range of possible HPA axis response patterns that occur after cardiac surgery. It also means we were unable to fully investigate the patient and operative factors that may cause or are associated with the different patterns of physiological processes. We envisage that novel blood-free biosampling technologies 27 will facilitate sampling a larger number of patients. This will allow a better characterisation of CABG phenotypes and investigate the effects of cytokines over a longer timeframe.
Our model does not consider inflammatory cytokines either as state variables or model inputs (as was the case with ACTH). Nor does it consider any explicit modulation by inflammatory mediators on ACTH -our open loop model architecture assumes such modulation is already accounted for in the resulting ACTH dynamics. Instead, we try to capture potential time-varying modulations of cortisol using parameter changes that remain fixed over time. This could explain the increased error in the CABG patient fits relative to the control fits. This is to be expected given the inflammatory response is only present in CABG patients. However, including inflammatory cytokines as model inputs would not only require additional parameters and therefore increase the risk of overfitting, but would also require detailed knowledge of the interactions between the HPA axis and inflammatory pathways. Although some advances have been made in modelling these interactions 9,28 , uncovering the network architecture between these pathways after major surgery and in critical illness has so far not been achieved. To do so will require a larger cohort of patients and a combination of mathematical modelling and machine learning techniques 29 .
We have been able to show that that there is not a simple graded HPA response to cardiac surgery. There are major dynamic changes involving not only the concentration of ACTH and cortisol but also their pattern of secretion. Using novel mathematical techniques, we show that CABG patients exhibit three different patterns of HPA axis response, which reflect different underlying physiological changes in adrenal sensitivity, cortisol production and turnover. Inflammatory mediators appear to be driving changes in only one of these patterns -the pattern of a single-pulse implying loss of ultradian pulsatility -despite being postulated as the likely cause for the elevated cortisol seen in all types of systemic inflammation 19,28 . We suggest that the different patterns we have observed may represent different "severities" of response to the surgery. An adequately powered study is now required to investigate whether and how these patterns are correlated with clinical outcomes. This will be critical to establish whether the patterns could be used for risk stratification after surgery. This study also shows that the existing model of corticosteroid physiology used for diagnosis and prognosis after major surgery and in critical illness 30 may only represent mean values within a population rather than the responses of individuals or groups of individuals. Improved diagnostics based on individual or subgroup responses is likely to lead to greater precision in diagnosis and more targeted interventions.

METHODS
The study was reviewed and approved by the UK National Research Ethics Service (NRES,  Fig. 2A).
Simpler, stationary cross-correlation methods were used to estimate the potential effects of inflammatory mediators not accounted for in the model. These involved a Z-score normalisation of the residual error of the model predictions vs cortisol data, and calculating its cross-correlation with the Z-score normalisation of the inflammatory mediators (Fig. 5A).
The results are summarised in a cross-correlation matrix for each patient, each cytokine, and some key combinations of these cytokines (Fig. 5B). We also performed a pairwise principal component analysis of the hormones ACTH and cortisol, as well as the inflammatory cytokines, to determine how similar their dynamic profiles were to one another. We did this on Z-scored data and used the MATLAB function pca to identify the variance explained by the two principal components, which were identified as IL6 and TNFα.
For each pair, we found the principal components by concatenating all patient data together. Therefore, the percentages of variance were calculated at the patient group level (Fig. S4). The model splits the cortisol dynamics C(t) as the contribution of two compartments, C f and C s , respectively accounting for fast and slow time scales for cortisol dynamics. Parameters p f and p s respectively account for the fast and slow cortisol production rates (including adrenal maximum secretory capacity and conversion from cortisone); 1/l f and 1/l s are the fast and slow cortisol turnover rates, respectively; 1/K A is the adrenal sensitivity to ACTH stimulation; and m is an integer denoting the power of the Hill function. A is the ACTH data, for which we assume a 10-minute delay before its effects are reflected onto cortisol 9 ( Fig.   2A). Overall, the model has six parameters, some of which we fix at different stages in the study. The initial condition for the model is chosen by setting the fast compartment C f to quasi-equilibrium, then 4 0 = max ( ?@A@ 0 − 2 0 , 0).
Optimisation. We use the cost function: where ?@A@ is the cortisol data for a given control or patient and ;B?CD is the model cortisol output for a given parameter set (using the control or patient ACTH data as input).
We derive the sample parameter sets from a latin hypercube of 1024 J points, capture the cost function for each patient or control and for each parameter set. This results in two matrices: = × = × Where = 1024 J the number of parameter sets on the hypercube, = 6 the number of parameters, and = 13 the number of patients (10) plus controls (3). Parcube is the list of parameters and errcube is the list of errors. We can then sort errcube and parcube by an individual control or patient to find the best-fit parameter choices. In some cases, we also sort by the patient group as a whole (or a subgroup of patients) or the control group. In this case, we use the mean-square of the error (MSE) across the given group to sort the parameters. This allows us to identify parameter sets that are best-fit for groups of multiple individuals. In some cases, we also fix some of the parameters and run a new latin hypercube. For example, we identified the best choice of the exponent parameter using an initial hypercube. We found that most good fits were for = 2, so we fixed that parameter and identified a new hypercube of 5 parameters for subsequent simulations. We identify good fits as any solution with an error value less than 10% greater than the best-fit error value. Therefore, we identify a distribution of good-fit solutions to minimize issues related to overfitting.

DATA AVAILABILITY
Raw data on hormone profiles, inflammatory mediators and surgery times, together with the time series analysis implemented in Python 3.6 for CABG phenotyping is available at https://gitlab.com/ezavala1/CABG_phenotypes The full model and optimisation pipeline was implemented in MATLAB and is available at https://github.com/dgalvis/heart_surgery_model