Abstract: Hypothermia is often used to treat out of hospital cardiac arrest (OHCA) patients who often simultaneously receive insulin for stress induced hyperglycaemia. Variations in response to insulin reflect dynamic changes in insulin sensitivity (SI), defined by the overall metabolic response to stress and therapy. Thus, tracking and forecasting this parameter is important to provide safe glycaemic control in highly dynamic patients. This study examines stochastic forecasting models of model-based SI variability in OHCA patients to assess the resulting potential impact of this therapy on glycaemic control quality and safety. A retrospective analysis of clinically validated model-based SI profiles identified using data from 240 post-cardiac arrest patients (9988 h) treated with hypothermia, shortly after admission in the Intensive Care Unit (ICU). Data were divided into three periods: (1) Cool (T≥35°C), (2) Idle period of 2 h as hypothermia was removed and (3) Warm (T≥37°C). The stochastic model captured 60.7 and 90.2% of SI predictions within the (25-75th) and (5-95th) probability forecast intervals during cool period. Equally, it is also recorded 62.8 and 92.1% of SI predictions respectively during the warm period. Maintaining the kernel density variance estimator to c = 1.0 yielded 60.7 and 90.2% for the cool period. Similarly, adjusting a variance estimator of c = 2.0 yields 60.4 and 90.1% for the warm period. A cohort-specific stochastic model of SI provided a conservative forecast for the inter-quartile range and was relatively exact for the 90% range. Adjusting the variance estimator provides a more accurate, cohort-speciWc stochastic model of SI dynamics for the 90% range. These latter results show clearly different levels and distribution of forecasted SI variability between the cold and warm periods.
INTRODUCTION
Hyperglycaemia is prevalent in critical care1-4 and increases the risks of further complications and mortality1,4,5. Glycaemic control has shown benefits in reducing mortality4,6,7. However, due in parts to excessive metabolic variability8, other studies have found it difficult to reproduce these results9-11. Out-of-hospital cardiac arrest (OHCA) patients have low survival rates and often experience hyperglycaemia12,13. However, cardiac patients are one group who have more consistently shown benefit from Glycaemic Control (GC), but can be highly insulin resistant and variable, particularly on the first day of stay14.
Hypothermia is often used to treat OHCA patients13. In general, it leads to a lowering of metabolic rate that induces changes in energy metabolism. However, its impact on metabolism and insulin resistance in critical illness is unknown, although one of the adverse events associated with hypothermic therapy is a decrease in insulin sensitivity and insulin secretion15. However, this decrease may not be notable in a cohort that is already highly resistant and variable14. Hence, understanding metabolic evolution and variability would enable safer, more accurate GC in this cohort.
Model-based glycaemic control methods using both insulin and/or nutrition modulation have been employed successfully in the control of hyperglycaemia. These methods allow the derivation of patient metabolic state, SI in this case, by using serial blood glucose (BG) measurements and records of nutrition and insulin administration16. Once the current SI has been identified, prediction of future SI would allow predictions of outcome BG concentration for an intended clinical intervention17.
Variations in the SI parameter reflect the metabolic response to stress2 and drug therapy18. Thus, tracking and forecasting this parameter is important to provide safe glycaemic control in the highly dynamic out-of-hospital cardiac arrest (OHCA) patients, treated with hypothermia. Since stochastic modeling has shown its ability to quantify the probability of a future SI17, the resulting distribution of BG concentrations that would result from a given intervention can be determined17,19,20,21. This information can be used to guide both insulin and/or nutrition interventions, which is the key to avoid unintended hypoglycaemia, improve overall glycaemic control and identify periods of potential high glucose variability that may be indicative of unusual clinical events or cohorts22.
This study presents the adaptation of a stochastic model for SI prediction from adult critical care to the unique clinical and physiological case of OHCA patients, treated with hypothermia. The stochastic model presented is constructed by the distribution of insulin sensitivity variation using a 2-D kernel density method. Clinically validated, model-based insulin sensitivity (SI)20,23 is used to provide more accurate measure of patient metabolic state and its stochastic model during cool and warm periods. Modifications to the initial kernel density estimation model are made to explore and optimize the relationship between the model and the underlying dataset, particularly for the clinically relevant 90% range that can be used to quantify and minimize the risks of hyper or hypo-glycaemia as an integral part of a control protocol20,21.
MATERIALS AND METHODS
Patients and data: A retrospective analysis of glycaemic control data from 240 OHCA patients (9988 h) treated with hypothermia, shortly after admission to the Intensive Care Unit (ICU) at Christchurch Hospital, New Zealand, Erasme Hospital, Belgium and Lausanne Hospital, Switzerland. Patients from Christchurch Hospital (20) were on the SPRINT glycaemic control protocol7, whereas the remaining 160 patients from Erasme (82) and Lausanne (78) Hospitals were on unpublished local glycaemic protocols. Blood Glucose (BG) and temperature readings were taken 1-2 h along with all insulin and nutrition data. Data were divided into three periods: (1) Cool (T≤35°C) (2) Transition period of 2 h as hypothermia was removed and (3) Warm after (T≥37°C). A maximum of 24 h and a minimum of 15 h for each period (1 and 3) were considered, ensuring a balance of contiguous data between periods. Overall demographics are shown in Table 1.
Therapeutic hypothermia: Therapeutic Hypothermia (TH) was applied following a standardized written protocol. All patients were treated with mild TH to 33±1°C for up to 24 h, irrespective of age, initial arrest rhythm and other physiological conditions. The TH was started immediately after admission and was induced with ice-cold packs and/or intravenous cold fluids.
| Table 1: | Demographic data for all patients |
| Data are presented as median (interquartile range) where appropriate | |
Body temperature was maintained at hypothermia using a surface cooling device with a computerized adjustment of patient temperature target. During this time, some short-acting drugs, such as midazolam (0.1 mg kg–1 h–1), fentanyl (1.5 μg kg–1 h–1) and vecuronium (0.1 mg kg–1 boluses) were used to administer sedation, analgesia and control shivering. Rewarming was achieved passively and sedation-analgesia was typically stopped when patient temperature was greater than 35°C.
Metabolic system model: Model-based insulin sensitivity (SI) in this study is a patient-specific parameter describing the whole body effect of insulin. The analysis of patient-specific insulin sensitivity employs the ICING model24 as a glucose-insulin system model developed and clinically validated in critical care glycaemic control and insulin sensitivity studies20,21,23,24. It is defined as:
| (1) |
| (2) |
| (3) |
| (4) |
| (5) |
| (6) |
| (7) |
A summary of parameter values and descriptions and exogenous input variables for the model are listed in Table 2 and 3, respectively.
Insulin sensitivity SI is identified hourly from patient data, producing a step-wise hourly varying profile25. This profile effectively describes patient-specific metabolic behavior under time-varying physiologic conditions. The validity and independence of this patient-specific parameter have been validated using data from independent, clinically matched cohorts23 and in gold-standard insulin sensitivity tests26.
Stochastic model: A 2-D kernel density estimation method is used to construct the stochastic model that describes the hourly transition of SI. The kernel density method combines probability distribution functions for each point of data to generate an overall density function for the dataset. This method has the advantage of producing a smooth, physiologically likely, continuous function across the parameter range to provide continuity when interpolating SI forecasts to account for each particular patient state.
| Table 2: | Parameter values and descriptions for the ICING model |
| Table 3: | Exogenous model input variables |
| Table 4: | Descriptions of 6 h blocks for data analysis |
It also automatically accounts for any possible multimodality where the density of data may show several distinct peaks corresponding to patterns of changes in SI. The overall result is a bivariate probability density function for the potential parameter values. The goal of this statistical model is to quantify the range of SI 1 h ahead in time (SI, n+1) based on available data (SI, n, SI ,n-1, SI, n-2, ..., SI, 0) to guide real-time clinical control22,21. Thus, it is potentially important that the model is also as cohort-specific as possible for greatest accuracy and to minimize over-conservative forecasts22.
A 2-D kernel density method is chosen because the distribution of SI,n+1 varies with SI,n and cannot be simply described with a single standard statistical distribution. Thus, the variations in SI can be treated as a Markov process. A Markov process has the property that the conditional probability density function of future states of the process, given the current state, depends only upon the current state. Therefore, using the Markov property of the stochastic behaviour of SI, the conditional probability density of SI,n+1 taking on a value y can be calculated by knowing SI,n = x. Model equations and derivation were defined17.
In this study, the model was cross validated by splitting the 240 of cool and warm patient cohorts into five groups, each containing 48 patients per period following the method of Lin et al.17. For each group, the model created using the remaining 192 patients of the cohort representing approximately 3990 h of data per period. Out-of-sample SI predictions were generated for the 48 patients of unused group and compared to the actual fitted SI from these 48 patients to assess model accuracy and cohort generality. Based on results from in-sample tests, where the stochastic model is generated from the entire retrospective dataset and tested on the same data and out-of-sample tests, where different subsets of data are used for model generation and testing, the kernel density estimator was modified by multiplying the variance estimators by a constant c (i.e., cσx and cσy) to explore the model bias-variance trade-off for this data between cool and warm periods as well as 6 and 12 h blocks. This adjustment to the variance estimator effectively adjusts the kernel bandwidth and the degree of smoothing over the data, which is used here to assess variability differences between periods or cohorts.
Analyses and metrics: Current SI (n) during the cool (T≤35°C) and warm (T≥37°C) periods were identified hourly using the ICING model24 for each patient. Using current SI (n) data, the predicted SI (n+1) data can be generalized by shifting 1 h forward from the current SI (n) data. These SI vectors will be paired and used to create the stochastic SI model of OHCA cohort. Stochastic SI models were analyzed during both cool and warm periods as follows:
| • | Overall cohort patient |
| • | Analysis of patients in 6 h block |
Overall cohort analysis assessed the stochastic model behavior of insulin sensitivity during both cool and warm periods, which includes percentage of SI within prediction interval and analysis of modifying kernel density estimation. The SI is also analyzed using 6 h blocks as described in Table 4, to capture SI variability and forecasting over time with different resolution. The analysis includes percentage of SI within predicted interval and analysis of modifying the kernel density estimation.
RESULTS AND DISCUSSION
Distribution of hourly variation in SI: Figure 1 presents the distribution of hourly variation in SI and the stochastic model percentile probability bands at c = 1.0 for the 240 OHCA patients during cool (4987 h) and warm (5001 h) periods. Approximately 85% of the values during cool period and 70% during warm period are below 1.0×10-3 L mU–1 min–1. The results show that the hourly variation of SI is wider during the cool period and SI rises from cool to warm.
Table 5 shows the in-sample results of stochastic model prediction widths for both cool (ncool = 4622 predictions) and warm (nwarm = 4832 predictions) periods. The number of predictions is less than the total hours of SI as the patient data records are not always perfectly divisible by 1 h and because predictions can only be computed after the second hour of patient data.
For the cool period, SI predictions (60.7%) were within the (25-75th) probability intervals and within the (5-95th) probability interval. Thus, the proportion of SI values that fell within the forecast (25-75th) probability intervals were measurably higher than the expected 50%, but very much closer for the 90% intervals. Similar results could be observed for the warm period, but were measured higher than cool period for both 50 and 90% intervals.
The percentage SI within prediction interval results of Lin et al.24 of 54.0% is within the (25-75th) probability bound shows that normal adult ICU patients with normal body temperature produce far closer to the ideal 50%. However, the prediction interval percentage results of 60.7 and 62.8% at cool and warm respectively on the same probability bound for the OHCA patients, treated with hypothermia as shown in Table 5 appear to be unique and significant. These results are similar to neonate’s results of Le Compte et al.19, which record 62.6% at the same interval.
Cross-validation comparison studies: Table 6 shows the results of the cross validation comparison study for 240 patient’s cohort during the cool and warm periods, respectively. Generally, these results are consistent between groups, suggesting that the overall model contains sufficient data to account for the range of dynamics observed in this cohort.
Probability-bound determination using local variance estimator: Table 7 shows the effect of modifying the kernel density estimation for several values of c, ranging from 0.1-3.0. For this cohort, the increase of c>1.0, yield better coverage widths for the expected proportions. Thus, for cool period, the value c = 1.0 enables the best 90% interval coverage, but is conservative for the inter-quartile range. Similarly, c = 2.0 is best for the warm period.
| Table 5: | Sample results for in-sample stochastic model prediction widths at c = 1 |
| Data are presented as cohort median (ncool = 4622 predictions and nwarm = 4832 predictions) | |
| Table 6: | Cross-validation comparison study for 240 patient cohort |
| Fig. 1(a-b): | Probability interval and distribution of hourly variation in SI for OHCA patients, treated with hypothermia (a) During cool and (b) Warm (right) periods |
| Fig. 2(a-b): | Probability-bound determination for raw SI data and corresponding BG forecasted values based on an equal-tailed 0.90 probability interval of SI between 5th -95th at local variance estimator, both cool (c = 1.0) and warm (c = 2.0) period. The solid lines represent the 5, 25, 50, 75 and 95% probability bounds |
| Table 7: | Comparison of probability bounds for modifications of kernel density estimator (σ`x = cσx and σ`y = cσy) during both cool and warm periods |
The difference of optimal c values between cool and warm suggests that the variation and stochastic modeling for both periods are different, leading to potentially different control requirements to ensure safe glycaemic control in these highly dynamic patients and states. In particular, c = 1.0 is currently used in the STAR protocol27,28 in medical ICU patients. Hence, no change is required for the cool period. However, during the warm period, the stochastic model requires a wider and smoother probability distribution coverage at c = 2.0, which would necessitate a control change of this known change of state. These results are further reflected in Fig. 2.
The kernel density estimator method employed in this stochastic model provides a layer of safety as wider probability bounds would be more likely to capture dynamics and any changes not observed in the cohort. As the (5-95th) band is what has been used for control previously, these cohorts show that they are closer to ideal 90%. However, wider coverage bands may also have impact on glycaemic control performance. As the wider probability band might be useful to avoid potential hypoglycaemia, it may also force a controller to maintain a mildly hyperglycaemic state. Importantly, different level of smoothing are required for the cool and warm periods, indicating that cool patients are less multi-modal (more smoothed) and warm OHCA patients are more multi-modal (less smoothed), in comparison to broad ICU cohorts. Perhaps this scenario is caused by the suppression of metabolic activities during cool period.
Stochastic analysis by 6 h block: The analysis of stochastic model based on 6 h time block have used the same methodology as the overall cohort patient, except that the patient cohort data are separated in block hours. Table 8 presents the in-sample results of stochastic model prediction widths based on 6 hour block analysis at c = 1.0.
The results show that SI prediction interval coverage meets the prediction expectation despite exceeding its ideal values of 50 and 90% for both (25-75th) and (5-95th) interval width respectively.
| Table 8: | Stochastic model prediction widths based on 6 h block analysis |
| Table 9: | Comparison of probability bounds for modifications of kernel density estimator (σ`x = cσx and σ`y = cσy) based on 6 h block analysis |
It is also observed that the proportion of fitted SI values that fell within the (25-75th) and (5-95th) probability intervals for 6 h block analysis were measured higher than the expected 50 and 90%.
Table 9 shows the effect of modifying the kernel density estimation for several values of c, ranging from 0.1-3.0 based on 6 h block analysis. For this cohort, the results has shown that for block 1, c = 1.0 and followed by block 2-3 (c = 1.5), block 4-5 (c = 2.0), block 6-7 (c = 2.5) and block 8 (c = 3.0). The trend shows that the value of estimator, c is increased as SI increases from cool to warm and match with overall cohort stochastic model analysis as shown in the Table 7. This will lead to another idea of separating stochastic model for each time block to implement stochastic control. The difference of optimal c values between 6 h blocks suggests that the variation and stochastic modeling for each time block is different, leading to different control requirements to ensure safe glycaemic control in the highly dynamic conditions.
Overall, this stochastic method and analysis in this study provides predictions based on a cohort dataset. The prediction bounds for more dynamic patients are difficult to decide since the SI level and variability distribution for this cohort is unique29 and tracked its evolution over time, particularly during cool period. This observation is far differing than for the less dynamic patients who are typically more conservative.
CONCLUSION
Thus, in this study the probability bounds are optimized in a cohort sense, but not necessarily applicable on a per-patient basis.
ACKNOWLEDGMENTS
I would like to thank my supervisory team, headed by Prof. Dr Geoff Chase, Dr Thomas Desaive, Dr. Chris Pretty and my colleagues at the Center for Bio-Engineering, Mechanical Engineering, University of Canterbury, Christchurch, New Zealand who provided insight and expertise that greatly assisted the research.