 Research
 Open Access
 Published:
3D kerneldensity stochastic model for more personalized glycaemic control: development and insilico validation
BioMedical Engineering OnLine volumeÂ 18, ArticleÂ number:Â 102 (2019)
Abstract
Background
The challenges of glycaemic control in critically ill patients have been debated for 20Â years. While glycaemic control shows benefits inter and intrapatient metabolic variability results in increased hypoglycaemia and glycaemic variability, both increasing morbidity and mortality. Hence, current recommendations for glycaemic control target higher glycaemic ranges, guided by the fear of harm. Lately, studies have proven the ability to provide safe, effective control for lower, normoglycaemic, ranges, using modelbased computerised methods. Such methods usually identify patientspecific physiological parameters to personalize titration of insulin and/or nutrition. The StochasticTargeted (STAR) glycaemic control framework uses patientspecific insulin sensitivity and a stochastic model of its future variability to directly account for both inter and intrapatient variability in a riskbased insulindosing approach.
Results
In this study, a more personalized and specific 3D version of the stochastic model used in STAR is compared to the current 2D stochastic model, both built using kerneldensity estimation methods. Fivefold cross validation on 681 retrospective patient glycaemic control episodes, totalling over 65,000Â h of control, is used to determine whether the 3D model better captures metabolic variability, and the potential gain in glycaemic outcome is assessed using validated virtual trials. Results show that the 3D stochastic model has similar forward predictive power, but provides significantly tighter, more patientspecific, prediction ranges, showing the 2D model overconservativeâ€‰>â€‰70% of the time. Virtual trial results show that overall glycaemic safety and performance are similar, but the 3D stochastic model reduced median blood glucose levels (6.3 [5.7, 7.0] vs. 6.2 [5.6, 6.9]) with a higher 61% vs. 56% of blood glucose within the 4.4â€“6.5Â mmol/L range.
Conclusions
This improved performance is achieved with higher insulin rates and higher carbohydrate intake, but no loss in safety from hypoglycaemia. Thus, the 3D stochastic model developed better characterises patientspecific future insulin sensitivity dynamics, resulting in improved simulated glycaemic outcomes and a greater level of personalization in control. The results justify inclusion into ongoing clinical use of STAR.
Background
Stress hyperglycaemia, or abnormal elevated blood glucose (BG) concentrations resulting from metabolic stress response to injury, is a common complication in critically ill patients, associated with increased morbidity and mortality [1,2,3,4]. Glycaemic control (GC) using insulin therapy to reduce BG to safer concentrations has shown improved outcomes, reducing organ failure, clinical burden, and costs [5,6,7,8,9,10]. However, other studies failed to replicate these results [11,12,13,14,15,16], showing increased glycaemic variability and higher risk of hypoglycaemia, independently associated with severe complications and death [17,18,19,20,21]. GC has been hard to achieve safely and effectively, often lacking patient specificity and failing to account for inter and intrapatient variability [22], showing the critical need for modelbased patientspecific GC solutions.
To date, the optimal target band for GC is still being debated [23]. Most intensive care units (ICUs) use a higher target band than the normoglycaemic range as a â€˜first do not harmâ€™ approach, hypoglycaemia being more harmful [17] for the patient than the potential benefits from GC. However, these standards are based on studies failing to provide safe, effective control for all patients when targeting a lower glycaemic band [24]. In fact, the association between mortality and glycaemic levels, safety, and variability has been shown a function of the control provided and not patient condition [25].
Hence, GC design is the key factor in patient GC outcomes. Failing to achieve safe, effective control for all, regardless of the target band, could bias study results and conclusions [26, 27]. More recently, new studies showed the possibility to achieve safe GC to lower target bands for reduced workload, without increasing hypoglycaemia [28,29,30]. These recent analyses show intensive GC to lower target bands is possible in critically ill patients, and emphasise, once again, the importance of safe, effective modelbased and patientspecific GC protocol design.
StochasticTARgeted (STAR) is a validated modelbased GC framework, which has shown effective control in three different countries [28]. STAR uses a physiological model to describe the glucoseâ€“insulin dynamics, and a populationbased stochastic model to directly account for patientspecific future metabolic variability [31]. Patientspecific insulin sensitivity (SI) key parameter [22, 32], describing patient metabolic response to insulin, is identified from clinical data and its future potential variability assessed to adjust treatment according to potential risks. This unique patientspecific riskbased dosing approach minimizes the risk of hypoglycaemia while providing effective GC [28, 31].
The quality of control resides in the ability to adapt treatment to patientspecific needs, which is a function of the level of difficulty to control [25]. The difficulty to control is mainly captured by SI variability [33], where variability extremes can lead to hyper or hypoglycaemia for a given insulin intervention. In STAR, SI is considered constant over 1h periods, and the variability is assessed by the hourtohour change in SI levels. The prediction of future SI evolution is thus a key element for the quality of GC, which needs to be well addressed, especially since SI variability is equivalent among patients with different clinical outcomes [25].
STAR predicts future SI evolution using stochastic models [34]. The stochastic model in STAR forecasts future SI (SI_{n+1}) distributions based on the current identified SI value (SI_{n}). However, the current stochastic models [34,35,36] are potentially overconservative due to large prediction bands. Wide prediction bands can limit insulin dosing, resulting in lower insulin doses to avoid stochastically forecasted hypoglycaemic risk. Better control might be obtained from a more detailed, and thus more personalized stochastic models. This study aims to improve SI variability forecasting using additional data information.
Using prior temporal information of SI evolution has shown better prediction accuracy [37]. Using such model gives generally tighter prediction bands at a given SI level, allowing potentially higher insulin rates for the patient. While encouraging, the method in [37] lacks model resolution and definition, making comparison hard with the current 2D stochastic model used in STAR. This study thus more specifically aims to develop a new 3D stochastic model accounting for prior knowledge of SI evolution, using 2 inputs (SI_{n} and SI_{nâˆ’1}) to determine likely future SI_{n+1}. The added input can provide higher patient specificity, allowing more accurate insulin dosage for the patient. More specifically, a wider future prediction range for SI would suggest higher potential variability; thus, lower insulin rates will likely be recommended. In contrast, tighter prediction bands would suggest lower variability and thus, potentially higher insulin recommendation. In this study, the new 3D stochastic model is generated using multivariate kerneldensity estimation method similar to the one used for the current 2D. This will improve model resolution. In addition, this study assesses the impact of this new presented 3D stochastic model on GC performances, using virtual trial simulations.
Results
2D vs. 3D models forward predictive power
A representation of the kerneldensity estimation is shown in Fig.Â 1. The left panel shows the kerneldensity surface using the normal data, whereas the right panel shows the kerneldensity surface when data are transformed into the lognormal space to meet the normal distribution assumption under Silvermanâ€™sÂ rule of thumb (ROT) [38]. Clearly, lognormal data provide increased data density for higher SI ranges, where the raw data are sparser. Hence, this approach, taken for the first time here, potentially improves safety by better characterising SI potential variability for higher SI ranges, where the risk of experiencing hypoglycaemia due to insulin dosing is greater.
Crossvalidation resultsâ€™ summary of the forward predictive power for both models is presented in TableÂ 1. In addition, the resulting 5th and 95th percentile predictions for each model are shown in Fig.Â 2. Both the 2D and 3D models have close to 50% (~â€‰53% vs.â€‰~â€‰51%) and 90% (~â€‰91% vs.â€‰~â€‰90%) predictions in the 25thâ€“75th and 5thâ€“95th percentile ranges, respectively. However, the prediction ranges are generally narrower (~â€‰70% of hours) in the case of the 3D model. An example of the evolution of SI for a patient and the 2D and 3D predictions ranges for a specific virtual patient is shown in Fig.Â 3. In addition, the median [IQR] percentage predictions in the 25thâ€“75th and 5thâ€“95th percentile prediction ranges are closer to the expected 50% and 90% for the 3D model, suggesting that the 2D model is too conservative for most patients. To characterise the difference in prediction ranges from both models, the percentage change in the 5thâ€“95th percentile range widths is computed for every prediction and the median [IQR] of percentage change is reported in TableÂ 1. The high prediction performances are achieved with significantly 15.5â€“24.4% tighter 5thâ€“95th percentile prediction range 69.9â€“73.8% of the time and 14.8â€“22% wider otherwise. The median [IQR] 3D/2D prediction width ratios as a function of the hourtohour percentage change in of SI (%Î”SI) are shown in Fig.Â 4, where clearly, prediction bands are typically tighter when Â %Î”SI is withinâ€‰Â±â€‰20%. Overall, the new model thus better captures patientspecific differences from this more optimal model.
Virtual trialsâ€™ results
Virtual trial results of STAR using the two different stochastic models are summarised in TableÂ 2. Overall, both versions of STAR provided similar performance in terms of median BG [IQR] (6.3 [5.7, 7.0] vs. 6.2 [5.6, 6.9] mmol/L) and percentage time in the 4.4â€“8.0Â mmol/L target band (88%). However, the overall Â %BG measurements shifted toward lower BG ranges using STAR3D, with significantly higher Â %BG within 4.4â€“6.5Â mmol/L and 4.4â€“7.0Â mmol/L (61% vs. 56% and 75% vs. 72%, pâ€‰<â€‰0.01 using Ï‡^{2} statistical test on proportions of measurements). In terms of safety, both models excel similar to only 2% BGâ€‰<â€‰4.4Â mmol/L, 1% BGâ€‰<â€‰4.0Â mmol/L, and 0.03% BGâ€‰<â€‰2.2Â mmol/L, despite STAR3D administering higher median insulin (3.0 [1.5, 5.0] vs. 2.5 [1.5, 4.5] U/h). Slightly lower, but similar, Â %BG in 8â€“10Â mmol/L (mild hyperglycaemia) for STAR3D is also observed (7% vs. 8%). Finally, STAR3D provided higher goal feed (97 [36, 100] vs. 95 [40, 100] Â %GF]).
Discussion
One of the key factors making GC difficult is patient variability. The riskbased dosing approach used in STAR relies on stochastic forecasting of likely future variation of SI, which is used to assess the risks of hyper and, more importantly, hypoglycaemia, and select appropriate treatment based on this risk. If the 5thâ€“95th percentile in SI prediction range is narrowed, control can be improved using more aggressive insulin dosing and/or greater dextrose intake to safely reach the same or lower glycaemic range. On the other hand, if the 5thâ€“95th percentile prediction range for the 3D stochastic model is widened, it suggests greater potential variability in future SI, leading to less aggressive insulin dosing to overcome the higher risk of hypoglycaemia.
The comparison between the 2D and 3D models clearly shows the new model accuracy to predict future SI, with 15.5â€“24.4% tighter prediction range for more than 69.9â€“73.8% of the hours (TableÂ 1). Typically, the prediction range is tighter when Â %Î”SI is withinâ€‰Â±â€‰20% (Fig.Â 4). On the contrary, the prediction range is wider when the variation is larger thanâ€‰Â±â€‰20%. This key outcome thus suggests that the previous patientspecific metabolic variability has a direct impact on future SI forecasting. More specifically, this 3D model shows stable patients, with low previous variation in SI, and tends to remain stable, whereas more variable patients are more likely to have bigger future metabolic variations, as clearly shown in Figs.Â 2 and 3. Hence, the 2D stochastic model is overconservative in terms of insulin intervention for most patients. The 3D approach allows STAR to select more aggressive insulin dosing more than 69.9% of the time, while ensuring safety, using the proven riskbased dosing approach. Therefore, the resulting greater patient specificity implies better GC with lower glycaemic variability, and improved glycaemic outcomes.
Virtual trial results comparing STAR using the 2D and the 3D stochastic models confirmed these observations showing higher percentage time in normoglycaemic ranges, with 5% more time spent in the 4.4â€“6.5Â mmol/L range, for similar incidence of mild hypoglycaemia (BGâ€‰<â€‰4.4Â mmol/L). In addition, the 3D model resulted in more aggressive insulin dosing and higher feed rates for similar intervention workload. Higher caloric intake is associated with improved outcomes [39,40,41,42]. These outcomes confirm the 3D stochastic model, using prior information in SI variability, and achieve effective control for all patients using more aggressive insulin dosing without compromising safety. Hence, STAR3D offers a more patientspecific control, better accounting for either stable or very variable patients, potentially resulting in improved patient outcomes.
More importantly, the slightly lower median BG using STAR3D (6.3 [5.7, 7.0] vs. 6.2 [5.6, 6.9]) was achieved with significantly higher time (61% vs. 56%, pâ€‰<â€‰0.01 using Ï‡^{2} statistical test on proportions of measurements) in the 4.4â€“6.5Â mmol/L band and in the 4.4â€“7.0Â mmol/L band (75% vs. 72%, pâ€‰<â€‰0.01 using Ï‡^{2} statistical test on proportions of measurements). While the low values for these p values could be influenced by the large data set size [43, 44], this difference is also clinically significant, since larger values in these ranges have been associated with improved outcomes and higher odds of living [45,46,47]. In addition, there was a consistent, high, 88% BG in target band (4.4â€“8.0Â mmol/L). High percentage time in these ranges has all been associated with improved clinical outcomes in multiple independent studies [21, 45,46,47]. These results, together with the minimal risk of hypoglycaemia (<â€‰2%) and severe hyperglycaemia (<â€‰3%), prove the STAR framework design to be adapted for GC in critical care, to provide safe, effective control for all patients, and show GC to lower target ranges to be possible without compromising safety.
It is also important to note specific safety benefits of this new model are hard to highlight. First of all, because hypoglycemia, in STAR, is extremely infrequent, unlike many other protocols failing to achieve safe control [12, 15, 48,49,50,51]. Hence, the few hours where the 3D model enables a gain in potentially very harmful hypoglycemia due to highly variable SI is hard to see and overwhelmed in the overall high effectiveness of STAR. Thus, we examine improved performance, which is beneficial to patients, with equivalent safety.
To further show this, the following cumulative distribution functions of the ratio between the 5thâ€“95th percentile range widths of each model when the subsequent SI value is within the predicted range and when the prediction is outside this range is shown below (Fig.Â 5). When SI is within the predicted range (~â€‰90% of hours, TableÂ 1), the 3D model prediction band is tighterâ€‰>â€‰75% of the time. However, when the subsequent SI value is outside the predicted range (~â€‰10% of the time), the 3D model is alreadyâ€‰>â€‰55% of the time wider than the 2D model. This result suggests that when the subsequent SI value is outside the range, the 3D model is generally more conservative (with a wider interval predicted) despite SI being outside predicted range. However, when the subsequent SI is within the predicted range, it is far narrower. Thus, the 3D model is overall safer.
While the difference in the two models shown in TableÂ 1 is quite important, and the virtual trials showed higher performance (TableÂ 2), a greater difference in glycaemic outcomes might have been expected. First, this difference shows how the STAR framework is consistent and manages to control patients in a safe manner. Second, the difference in SI prediction ranges between these models may not be big enough to change the discretized insulin interventions in STAR, as the controller is limited to 0.5 U/h increments. More specifically, in [25], an analysis suggests that a change below 12â€“15% in SI levels can be considered clinically equivalent, limiting some impact on GC recommendations.
STAR treatment selection relies in putting the 5th percentile of predicted BG outcome on the lower target band limit. Hence, it mainly uses the 95th percentile of predicted SI. Looking deeper at the 95th percentile difference between the 2D and 3D model, there is median reduction ofâ€‰~â€‰6%, which may not be enough to significantly change the administration rate of insulin.
As reflected in these results, using more information to better predict how likely patientspecific metabolic conditions will change seems a good approach to improve control in the STAR framework. More specifically, using prior state of identified SI values also allows to suffer less from direct measurement errors or identification errors for the future prediction [52]. While one could think to extend this method to more dimensions, the danger would be to over fit the data and/or suffer from low data density. Therefore, it could result in nondesired behaviour for higher computational costs. However, other parameters could be useful to improve both predictions and GC outcomes. In [53], BG data are used as an entry with current SI level to forecast metabolic variability. In doing so, not only it potentially can improve control safety and efficacy, but it also allows to identify specific behaviour in the data, reflected by the resulting estimated distributions. In particular, [53] observed typical underestimation of SI changes at lower BG values and vice versa. Hence, more work could be done to identify possible critical factors or parameters allowing to further improve prediction of important changes in metabolic variability and SI.
The bivariate kerneldensity estimation method requires much fewer total data to create an effective model for use in clinical practice compared to the trivariate model presented. However, the 3D stochastic model demonstrated better performance and equivalent safety in this study due to the much higher number data triplets (~â€‰60,000 vs.â€‰~â€‰20,000) available from the larger population data set used in this study than in creating the 2D stochastic model [34]. In addition, the equivalence across the virtual trial fivefold cross validation results suggests that the stochastic models were created on enough data to be robust and that the data used were representative of a general ICU population.
The interpretation of these results has some limitations. Virtual trials represent realistic glycaemic outcomes in perfect conditions, fully compliant to the protocol [54]. Glycaemic outcomes will likely differ at least somewhat in a real clinical environment. However, these virtual trials have been validated and shown to well capture the overall potential glycaemic outcomes [55, 56]. In addition, compliance to STAR is very high in regular clinical use [28, 57].
Conclusions
Trivariate kerneldensity estimation methods are used to build a new 3D stochastic model forecasting likely future changes in insulin sensitivity based on its prior 2 states. This 3D stochastic model shows similar, high, forward predictive power compared to the previous 2D version, but achieved with 15â€“25% tighter prediction ranges more than 70% of the time. This suggests that the 3D stochastic model better predicts future SI dynamics and thus offers greater personalisation of care than the prior 2D model.
Virtual trials using this model showed similar glycaemic control safety and better performance based on higher time in the normoglycaemic intermediate ranges (4.4â€“6.5Â mmol/L and 4.4â€“7.0Â mmol/L), resulting in slightly lower median BG levels for similar workload. These improvements are due to greater personalisation of care, and were achieved using higher insulin rates and slightly higher nutrition rates in cases where possible and as enabled by the tighter prediction ranges offered in over 70% of interventions. These results suggest that the implementation of this new 3D stochastic model within the STAR framework could potentially improve patient clinical outcomes resulting from improved glycaemic control.
Methods
Modelbased insulin sensitivity
The validated Intensive Control Insulinâ€“Nutritionâ€“Glucose physiological (ICING) [58] model used in STAR describes the glucoseâ€“insulin pharmacokinetics and can be graphically represented as the 3 compartment model in Fig.Â 6. It is defined:
where G(t) is the blood glucose concentration (mmol/L), I(t) and Q(t) are the plasma and interstitial insulin concentrations (mU/L), P(t) is the glucose appearance in plasma from enteral and parenteral dextrose intakes (mmol/min), and SI is the insulin sensitivity (L/mU/min). Other parameters are listed in TableÂ 3. Clearance rates and constants can be found elsewhere [25, 52, 58].
The modelbased timevarying SI parameter describes the patientspecific glycaemic metabolic response to insulin administration. SI is identified hourly from clinical BG, insulin, and nutritional data using integralbased fitting methods [59]. The validity of this SI metric is demonstrated in [25] for ICU patients, as well as in several clinical studies [60,61,62].
Stochastic modelling
To account for patientspecific metabolic variability, and thus assess unexpected potential changes in metabolic response to insulin, [34] introduced a probabilistic model predicting likely future 1â€“3 hourly change in SI level (SI_{n+1}, SI_{n+2}, SI_{n+3}). These predictions are only based on current identified patient metabolic condition (SI_{n}). This stochastic model was built using a twodimensional kerneldensity estimation method on population data, and led to the emergence of the first successful riskbased dosing approach for GC [31, 63]. The kerneldensity estimation method enables highresolution behaviour estimation of a specific parameter based upon its prior evolution or state, even where specific data points may be scarce.
Using the identified SI_{n}, the 2D stochastic model forecasts likely future distribution of SI_{n+1}, as graphically represented in Fig.Â 7. This likely future SI distribution allows prediction of the corresponding likely future BG distribution for a given insulin and nutrition intervention using Eqs.Â 1â€“3 (Fig.Â 7). Specifically, the 5thâ€“95th percentile range of likely future SI is used to compute the corresponding 5thâ€“95th percentile range of predicted future BG outcomes. STAR then adjusts treatments by ensuring the 5th and 95th percentiles of future BG lie within the clinically specified target range (4.4â€“8.0Â mmol/L in STAR), minimizing the risk of BGâ€‰<â€‰4.4Â mmol/L to 5% [31].
This study extends the bivariate kernel estimation method to trivariate. The predictions of future SI_{n+1}, SI_{n+2}, and SI_{n+3} are thus determined using two inputs to potentially increase patientspecific variability forecasting, which could also result in better overall safety and performance for STAR GC decision making. In particular, this choice of data triplets (SI_{nâˆ’1}, SI_{n})â€‰â†’â€‰SI_{n+1,2,3}, add patient specificity to the SI_{n}â€‰â†’â€‰SI_{n+1,2,3} 2D model by making these distributions a function of prior states. This difference thus includes a greater part of the patientspecific evolution, and thus will further characterise patients, creating greater personalisation in the GC predictions based on thus enhanced stochastic model. It thus assumes that there will be measurable differences in the predicted SI_{n+1,2,3} distributions by this added data, compared to those from the 2D model. Importantly, the 3D approach significantly increases the data requirements for model generation, resulting in the use of a much larger data set size (~â€‰60,000Â h) compared to the previous studies [34].
From data density to conditional probability
SI in this study can be considered a secondorder finite Markov chain, where the current state depends only on its two prior states. Therefore, the conditional probability distribution of the future SI_{n+1} is a function of SI_{n} and SI_{n1} states which can be expressed:
where the righthand expression is derived from the general product rule. Kerneldensity methods are used to estimate the joint probability \(P({\text{SI}}_{n + 1} ,{\text{SI}}_{n} ,{\text{SI}}_{n  1} )\) and \(P({\text{SI}}_{n} ,{\text{SI}}_{n  1} )\) using tri and bivariate Gaussian kerneldensity estimator functions [64]. Therefore, the conditional probability of SI_{n+1} taking a specific value can be calculated using the identified SI_{n} and SI_{n1} values, such that
where \(K_{h} \left( u \right)\) denotes the Gaussian kerneldensity function \(K_{h} \left( u \right) = \frac{1}{{\sqrt {2\pi } h}}e^{{  \frac{1}{2}\left( {\frac{u}{h}} \right)^{2} }}\), centered on u with variance h, constructed using the available N data points [38, 65]. To optimize the approximation of data behaviour, the variance h, or scale factor, is determined using the general Silvermanâ€™s ROT [38, 64], weighted according to local data density:
where m is the number of data point within a radius \(N^{  1/7}\) after orthonormalisation of the data [34], and R is the radius from the origin encompassing Z*N data points (0â€‰â‰¤â€‰Zâ€‰â‰¤â€‰1). This rule assumes that data have an underlying normal distribution [38]. Nonnegativity is ensured by normalizing each Gaussian function to the positive defined domain, such that for each \(({\text{SI}}_{n} = y,{\text{SI}}_{n  1} = x)\) pair, there exists an estimated conditional probability function \(P\left( {{\text{SI}}_{n + 1} = z  {\text{SI}}_{n} = y,{\text{SI}}_{n  1} = x} \right)\), where \(\smallint P\left( {{\text{SI}}_{n + 1} = z  {\text{SI}}_{n} = y,{\text{SI}}_{n  1} = x} \right){\text{d}}z = 1\) is satisfied [64]. Normalization is achieved by dividing each kerneldensity function \(K_{hx_i} \left( u_{x_i} \right), K_{hy_i} \left( u_{y_i} \right), K_{hz_i} \left( u_{z_i} \right)\) in Eq. 5Â by the area under each gaussian curve between zero and infinity:
This forces x, y, and z to beâ€‰â‰¥â€‰0, ensuring thus physiological validity of SI valuesâ€‰â‰¥â€‰0. An example of the resulting uni, bi, and trivariate Gaussian kerneldensity estimation for 10 data triplets is shown in Fig.Â 8.
Patients and cohorts
This study uses clinical data from 606 patients across 3 different clinical trials and ICU settings (STAR Christchurch New Zealand 2011â€“2015, SPRINT Christchurch New Zealand 2005â€“2007, and STAR Gyula Hungary 2011â€“2015) [9, 28]. These data include 819 GC episodes and a total of 68,629Â h of treatment. A patient can have multiple GC episodes, generally because:

1.
Patientsâ€™ glycaemia is stabilized, but then several hours later GC is started again due to dysglycaemia arising from any potential clinical reason, or

2.
Patients are sent out of the ICU for clinical procedures (most commonly imaging or surgery), where GC is stopped and started again as they return (if necessary).
Overall cohort demographics are shown in TableÂ 4.
From the original 819 episodes, only 681 episodesâ€‰â‰¥â€‰10Â h and with initial BGâ€‰â‰¥â€‰7Â mmol/L are considered (Fig.Â 9), corresponding to 59,439Â h of control. These criteria ensure the exclusion of patient data with very short GC episodes, and thus low BG measurement numbers, or uncommonly low starting BG values likely less reflective of general metabolism dynamics. SI is identified hourly for each patient using integralbased fitting method and a total of 58,539, 57,840, and 57,141 data triplets (SI_{n1}, SI_{n}, SI_{n+i}) for iâ€‰=â€‰1, 2, and 3Â h forward, respectively, are created.
Validation and comparison analysis
The 2D and 3D stochastic models are built and compared using fivefold cross validation, where the resulting training (80%) and testing (20%) sets are believed to be statistically representative of the general data set, minimizing bias and variance in the validation [66]. Patients are thus randomly divided into 5 equally sized groups, models are built using 80% of patient episodes (4/5 groups), and the other 20% of patients (1/5 groups) are used for validation. As the Silvermanâ€™s ROT for multivariate kerneldensity estimation assumes data has a Gaussian distribution [38], and SI has a lognormal distribution, the logarithmic domain is chosen here to build the model.
The 25thâ€“75th and 5thâ€“95th percentile ranges are computed for both models. Tighter prediction ranges for future SI_{n+i} would suggest likely lower future variability. In this case, the future potential variation for a given insulin dose being smaller, STAR can provide insulin with greater certainty, and thus potentially more aggressively (higher insulin rates) with equal safety. On the other hand, wider prediction bands would suggest higher future variability and, thus, more conservative dosing of insulin use to avoid hypoglycaemia. Forward predictive power and model accuracy are compared using the percentage of accurate predictions within these two ranges. The expected accuracies are 50% and 90%, respectively, where conformation of an independent cohort to these expected outcomes would indicate the 3D methods accurately capture SI dynamics to predict future SI.
Finally, to assess clinical impact, validated virtual trials on virtual patients are simulated to assess the new models ability to control patients. Such virtual trials enable comparison of glycaemic outcomes from different GC designs, on the same underlying patients. In summary, virtual patients are characterised by their identified patientspecific SI traces generated from clinical data, and can be used to test a range of new protocols or technologies [67, 68]. They are wellvalidated in their independence from the data used to create them and their accuracy [56, 69], their ability to predict trial outcomes [63, 70] and in clinical use to guide care in STAR [28, 31]. The underlying model is also wellvalidated in insulin sensitivity testing and similar clinical studies [60, 71, 72]. These virtual trials have been validated in the previous studies [55, 56], and are used here to simulate STAR using either the 2D (STAR2D) or 3D (STAR3D) stochastic model.
Unlike most GC protocols, STAR has the ability to modulate both insulin and nutrition inputs. Enteral nutrition can be lowered if the maximum allowed insulin is not sufficient to decrease BG levels, often occurring for very resistant patients with low SI and saturation of insulindosing effects. In STAR, insulin is administered as boluses up to a maximum of 6Â U/h, with an additional 3Â U/h continuous infusion for highly resistant patients. Enteral nutrition administration can be modulated between 30 and 100% of the total calorific goal feed (GF) if necessary. The original 100% GF for a patient is computed according to the standard 25Â kcal/kg/day target adapted based on age and sex. Further details are in [28, 42].
Safety and performance, administered insulin, and nutrition delivery are compared from these simulations. BG is resampled hourly, to allow fair comparison across the different measurement intervals. Safety is assessed by the %BG in mild (%BGâ€‰â‰¤â€‰4.4Â mmol/L) and severe (%BGâ€‰â‰¤â€‰2.2Â mmol/L) hypoglycaemia and in hyperglycaemia (%BGâ€‰>â€‰8.0Â mmol/L and %BGâ€‰>â€‰10.0Â mmol/L). Performance is assessed by the %BG in target band (4.4â€“8.0Â mmol/L) and the median [IQR] BG levels achieved. Nutrition is reported as the percentage goal feed (%GF) achieved. In addition, workload is also compared as the number of BG measurement per day, where a higher value indicates increased workload.
Availability of data and materials
The data sets used and/or analysed during the current study are available from the corresponding author on reasonable request. However, a subset of the data is publicly available in another journal: Chase J.G. et al. A benchmark data set for modelbased glycemic control in critical care, J Diabetes Sci Technol 2008, 2(4): 584â€“94.
Abbreviations
 BG:

blood glucose
 GC:

glycaemic control
 ICING:

Intensive Control Insulinâ€“Nutritionâ€“Glucose
 ICU:

intensive care unit
 IQR:

interquartile range
 ROT:

rule of thumb
 STAR:

stochastic targeted
 SI:

insulin sensitivity
References
McCowen KC, Malhotra A, Bistrian BR. Stressinduced hyperglycemia. Crit Care Clin. 2001;17:107â€“24.
Ali NA, Oâ€™Brien JM, Dungan K, Phillips G, Marsh CB, Lemeshow S, Connors AF, Preiser JC. Glucose variability and mortality in patients with sepsis. Crit Care Med. 2008;36:2316â€“21.
Krinsley JS. Glycemic variability and mortality in critically ill patients: the impact of diabetes. J Diabetes Sci Technol. 2009;3:1292â€“301.
Capes SE, Hunt D, Malmberg K, Gerstein HC. Stress hyperglycaemia and increased risk of death after myocardial infarction in patients with and without diabetes: a systematic overview. Lancet. 2000;355:773â€“8.
Van den Berghe G, Wouters P, Weekers F, Verwaest C, Bruyninckx F, Schetz M, Vlasselaers D, Ferdinande P, Lauwers P, Bouillon R. Intensive insulin therapy in critically ill patients. N Engl J Med. 2001;345:1359â€“67.
Krinsley JS. Effect of an intensive glucose management protocol on the mortality of critically ill adult patients. Mayo Clin Proc. 2004;79:992â€“1000.
Reed CC, Stewart RM, Sherman M, Myers JG, Corneille MG, Larson N, Gerhardt S, Beadle R, Gamboa C, Dent D, et al. Intensive insulin protocol improves glucose control and is associated with a reduction in intensive care unit mortality. J Am Coll Surg. 2007;204:1048â€“54 (discussion 1054â€“1045).
Chase JG, Pretty CG, Pfeifer L, Shaw GM, Preiser JC, Le Compte AJ, Lin J, Hewett D, Moorhead KT, Desaive T. Organ failure and tight glycemic control in the SPRINT study. Crit Care. 2010;14:R154.
Chase JG, Shaw G, Le Compte A, Lonergan T, Willacy M, Wong XW, Lin J, Lotz T, Lee D, Hann C. Implementation and evaluation of the SPRINT protocol for tight glycaemic control in critically ill patients: a clinical practice change. Crit Care. 2008;12:R49.
Krinsley JS. Is glycemic control of the critically ill costeffective? Hosp Pract. 1995;2014(42):53â€“8.
Finfer S, Chittock D, Li Y, Foster D, Dhingra V, Bellomo R, Cook D, Dodek P, Hebert P, Henderson W, et al. Intensive versus conventional glucose control in critically ill patients with traumatic brain injury: longterm followup of a subgroup of patients from the NICESUGAR study. Intensive Care Med. 2015;41:1037â€“47.
Brunkhorst FM, Engel C, Bloos F, MeierHellmann A, Ragaller M, Weiler N, Moerer O, Gruendling M, Oppert M, Grond S, et al. Intensive insulin therapy and pentastarch resuscitation in severe sepsis. N Engl J Med. 2008;358:125â€“39.
Signal M, Fisk L, Shaw GM, Chase JG. Concurrent continuous glucose monitoring in critically Ill patients: interim results and observations. J Diabetes Sci Technol. 2013;7:1652â€“3.
Signal M, Pretty CG, Chase JG, Le Compte A, Shaw GM. Continuous glucose monitors and the burden of tight glycemic control in critical care: can they cure the time cost? J Diabetes Sci Technol. 2010;4:625â€“35.
Preiser JC, Devos P, RuizSantana S, Melot C, Annane D, Groeneveld J, Iapichino G, Leverve X, Nitenberg G, Singer P, et al. A prospective randomised multicentre controlled trial on tight glucose control by intensive insulin therapy in adult intensive care units: the Glucontrol study. Intensive Care Med. 2009;35:1738â€“48.
Van Herpe T, De Moor B, Van den Berghe G. Ingredients for adequate evaluation of blood glucose algorithms as applied to the critically ill. Crit Care. 2009;13:102.
Bagshaw SM, Bellomo R, Jacka MJ, Egi M, Hart GK, George C, Committee ACM. The impact of early hypoglycemia and blood glucose variability on outcome in critical illness. Crit Care. 2009;13:R91.
Egi M, Bellomo R, Stachowski E, French CJ, Hart G. Variability of blood glucose concentration and shortterm mortality in critically ill patients. Anesthesiology. 2006;105:244â€“52.
Egi M, Bellomo R, Stachowski E, French CJ, Hart GK, Taori G, Hegarty C, Bailey M. Hypoglycemia and outcome in critically ill patients. Mayo Clin Proc. 2010;85:217â€“24.
Krinsley JS. Association between hyperglycemia and increased hospital mortality in a heterogeneous population of critically ill patients. Mayo Clin Proc. 2003;78:1471â€“8.
Penning S, Pretty C, Preiser JC, Shaw GM, Desaive T, Chase JG. Glucose control positively influences patient outcome: a retrospective study. J Crit Care. 2015;30:455â€“9.
Chase JG, Le Compte AJ, Suhaimi F, Shaw GM, Lynn A, Lin J, Pretty CG, Razak N, Parente JD, Hann CE, et al. Tight glycemic control in critical careâ€“the leading role of insulin sensitivity and patient variability: a review and modelbased analysis. Comput Methods Programs Biomed. 2011;102:156â€“71.
Krinsley JS. Is it time to rethink blood glucose targets in critically ill patients? Chest. 2018;154:1004â€“5.
Roubicek T, Kremen J, Blaha J, Matias M, Kopecky P, Rulisek J, Anderlova K, Bosanska L, Mraz M, Chassin LJ, et al. Pilot study to evaluate blood glucose control by a model predictive control algorithm with variable sampling rate vs. routine glucose management protocol in peri and postoperative period in cardiac surgery patients. Cas Lek Cesk. 2007;146:868â€“73.
KuureKinsey M, Palerm CC, Bequette BW. A dualrate Kalman filter for continuous glucose monitoring. Conf Proc IEEE Eng Med Biol Soc. 2006;1:63â€“6.
Reifman J, Rajaraman S, Gribok A, Ward WK. Predictive monitoring for improved management of glucose levels. J Diabetes Sci Technol. 2007;1:478â€“86.
Uyttendaele V, Knopp JL, Shaw GM, Desaive T, Chase JG. Is intensive insulin therapy the scapegoat for or cause of hypoglycaemia and poor outcome? IFAC J Syst Control. 2019;9:100063.
Stewart KW, Pretty CG, Tomlinson H, Thomas FL, Homlok J, Noemi SN, Illyes A, Shaw GM, Benyo B, Chase JG. Safety, efficacy and clinical generalization of the STAR protocol: a retrospective analysis. Ann Intensive Care. 2016;6:24.
Lunn DJ, Wei C, Hovorka R. Fitting dynamic models with forcing functions: application to continuous glucose monitoring in insulin therapy. Stat Med. 2011;30:2234â€“50.
Schultz MJ, Harmsen RE, Korevaar JC, AbuHanna A, Van Braam Houckgeest F, Van Der Sluijs JP, Spronk PE. Adoption and implementation of the original strict glycemic control guideline is feasible and safe in adult critically ill patients. Minerva Anestesiol. 2012;78:982â€“95.
Evans A, Le Compte A, Tan CS, Ward L, Steel J, Pretty CG, Penning S, Suhaimi F, Shaw GM, Desaive T, Chase JG. Stochastic targeted (STAR) glycemic control: design, safety, and performance. J Diabetes Sci Technol. 2012;6:102â€“15.
Suhaimi F, Le Compte A, Preiser JC, Shaw GM, Massion P, Radermecker R, Pretty CG, Lin J, Desaive T, Chase JG. What makes tight glycemic control tight? The impact of variability and nutrition in two clinical studies. J Diabetes Sci Technol. 2010;4:284â€“98.
Uyttendaele V, Dickson JL, Shaw GM, Desaive T, Chase JG. Untangling glycaemia and mortality in critical care. Crit Care. 2017;21:152.
Lin J, Lee D, Chase JG, Shaw GM, Le Compte A, Lotz T, Wong J, Lonergan T, Hann CE. Stochastic modelling of insulin sensitivity and adaptive glycemic control for critical care. Comput Methods Programs Biomed. 2008;89:141â€“52.
Lin J, Lee D, Chase JG, Shaw GM, Hann CE, Lotz T, Wong J. Stochastic modelling of insulin sensitivity variability in critical care. Biomed Signal Process Control. 2006;1:229â€“42.
Vanhorebeek I, Langouche L, Van den Berghe G. Glycemic and nonglycemic effects of insulin: how do they contribute to a better outcome of critical illness? Curr Opin Crit Care. 2005;11:304â€“11.
Uyttendaele V, Knopp JL, Stewart K, Desaive T, Benyo B, SzaboNemedy N, Illyes A, Shaw G, Chase JG. A 3D insulin sensitivity prediction model enables more patientspecific prediction and modelbased glycaemic control. Biomed Signal Process Control. 2018;46:192â€“200.
Sheather SJ: Density estimation. In: Statistical science. Volume 19: Institute of mathematical statistics; 2004.[JSTOR (Series Editor).
Villet S, Chiolero RL, Bollmann MD, Revelly JP, Cayeux RNM, Delarue J, Berger MM. Negative impact of hypocaloric feeding and energy balance on clinical outcome in ICU patients. Clin Nutr. 2005;24:502â€“9.
Krishnan JA, Parce PB, Martinez A, Diette GB, Brower RG. Caloric intake in medical ICU patients: consistency of care with guidelines and relationship to clinical outcomes. Chest. 2003;124:297â€“305.
Heyland DK, Cahill N, Day AG. Optimal amount of calories for critically ill patients: depends on how you slice the cake! Crit Care Med. 2011;39:2619â€“26.
Stewart KW, Chase JG, Pretty CG, Shaw GM. Nutrition delivery of a modelbased ICU glycaemic control system. Ann Intensive Care. 2018;8:4.
Goodman SN. Toward evidencebased medical statistics. 1: The P value fallacy. Ann Intern Med. 1999;130:995â€“1004.
Motulsky H. Common misconceptions about data analysis and statistics. Br J Pharmacol. 2015;172:1017â€“23.
Signal M, Le Compte A, Shaw GM, Chase JG. Glycemic levels in critically ill patients: are normoglycemia and low variability associated with improved outcomes? J Diabetes Sci Technol. 2012;6:1030â€“7.
Penning S, Chase JG, Preiser JC, Pretty CG, Signal M, Melot C, Desaive T. Does the achievement of an intermediate glycemic target reduce organ failure and mortality? A post hoc analysis of the Glucontrol trial. J Crit Care. 2014;29:374â€“9.
Krinsley JS, Preiser JC. Time in blood glucose range 70 to 140Â mg/dlâ€‰>â€‰80% is strongly associated with increased survival in nondiabetic critically ill adults. Crit Care. 2015;19:179.
Finfer S, Chittock DR, Su SY, Blair D, Foster D, Dhingra V, Bellomo R, Cook D, Dodek P, Henderson WR, et al. Intensive versus conventional glucose control in critically ill patients. N Engl J Med. 2009;360:1283â€“97.
Arabi YM, Dabbagh OC, Tamim HM, AlShimemeri AA, Memish ZA, Haddad SH. Intensive versus conventional insulin therapy: a randomized controlled trial in medical and surgical critically ill patients. Crit Care Med. 2008;36:3190â€“7.
Rosa C, Donado JH, Restrepo AH, Quintero AM, Gonzalez LG, Saldarriaga NE. Strict glycaemic control in patients hospitalised in a mixed medical and surgical intensive care unit: a randomised clinical trial. Crit Care. 2008;12:120.
Treggiari MM, Karir V, Yanez ND, Weiss NS, Daniel S, Deem SA. Intensive insulin therapy and mortality in critically ill patients. Crit Care. 2008;12:R59.
Pretty CG, Signal M, Fisk L, Penning S, Le Compte A, Shaw GM, Desaive T, Chase JG. Impact of sensor and measurement timing errors on modelbased insulin sensitivity. Comput Methods Programs Biomed. 2014;114:e79â€“86.
Marik PE, Preiser JC. Toward understanding tight glycemic control in the ICU: a systematic review and metaanalysis. Chest. 2010;137:544â€“51.
Chase JG, Andreassen S, Jensen K, Shaw GM. Impact of human factors on clinical protocol performance: a proposed assessment framework and case examples. J Diabetes Sci Technol. 2008;2:409â€“16.
Dickson JL, Stewart KW, Pretty CG, Flechet M, Desaive T, Penning S, Lambermont BC, Benyo B, Shaw GM, Chase JG, et al. Generalisability of a virtual trials method for glycaemic control in intensive care. IEEE Trans Biomed Eng. 2018;65:1543â€“53.
Chase JG, Suhaimi F, Penning S, Preiser JC, Le Compte AJ, Lin J, Pretty CG, Shaw GM, Moorhead KT, Desaive T. Validation of a modelbased virtual trials method for tight glycemic control in intensive care. Biomed Eng Online. 2010;9:84.
Dickson J, Chase JG. Clinical Compliance in Personalised Modelbased Medical Decision Support: do computers and interfaces yield better compliance? IFACPapersOnLine. 2019;51:341â€“6.
Lin J, Razak NN, Pretty CG, Le Compte A, Docherty P, Parente JD, Shaw GM, Hann CE, Geoffrey Chase J. A physiological Intensive Control Insulinâ€“Nutritionâ€“Glucose (ICING) model validated in critically ill patients. Comput Methods Programs Biomed. 2011;102:192â€“205.
Dandona P, Mohanty P, Chaudhuri A, Garg R, Aljada A. Insulin infusion in acute illness. J Clin Invest. 2005;115:2069â€“72.
McAuley KA, Berkeley JE, Docherty PD, Lotz TF, Te Morenga LA, Shaw GM, Williams SM, Chase JG, Mann JI. The dynamic insulin sensitivity and secretion testâ€”a novel measure of insulin sensitivity. Metabolism. 2011;60:1748â€“56.
Docherty PD, Chase JG, Lotz T, Hann CE, Shaw GM, Berkeley JE, Mann JI, McAuley K. DISTq: an iterative analysis of glucose data for lowcost, realtime and accurate estimation of insulin sensitivity. Open Med Inform J. 2009;3:65â€“76.
Docherty PD, Chase JG, Te Morenga L, Fisk LM. A novel hierarchalbased approach to measure insulin sensitivity and secretion in atrisk populations. J Diabetes Sci Technol. 2014;8:807â€“14.
Fisk LM, Le Compte AJ, Shaw GM, Penning S, Desaive T, Chase JG. STAR development and protocol comparison. IEEE Trans Biomed Eng. 2012;59:3357â€“64.
Silverman BW. Density estimation for statistics and data analysis. London: Chapman and Hall; 1986.
Scott DW. Multivariate density estimation and visualization. In: Gentle JE, HÃ¤rdle WK, Mori Y, editors. Handbook of computational statistics: concepts and methods. Berlin: Springer; 2012. p. 549â€“69.
James G, Witten D, Hastie T, Tibshirani R. Resampling methods. An introduction to statistical learning. New York: Springer; 2013. p. 175â€“201.
Chase JG, Preiser JC, Dickson JL, Pironet A, Chiew YS, Pretty CG, Shaw GM, Benyo B, Moeller K, Safaei S, et al. Nextgeneration, personalised, modelbased critical care medicine: a stateofthe art review of in silico virtual patient models, methods, and cohorts, and how to validation them. Biomed Eng Online. 2018;17:24.
Chase JG, Benyo B, Desaive T: Glycemic control in the intensive care unit: a control systems perspective. Annual Reviews in Control 2019.
Dickson JL, Stewart KW, Pretty CG, Flechet M, Desaive T, Penning S, Lambermont BC, Benyo B, Shaw GM, Chase G. Generalisability of a virtual trials method for glycaemic control in intensive care. In: IEEE transactions on biomedical engineering 2017, p. 1.
Lonergan T, Le Compte A, Willacy M, Chase JG, Shaw GM, Wong XW, Lotz T, Lin J, Hann CE. A simple insulinnutrition protocol for tight glycemic control in critical illness: development and protocol comparison. Diabetes Technol Ther. 2006;8:191â€“206.
Docherty PD, Chase JG, David T. Characterisation of the iterative integral parameter identification method. Med Biol Eng Comput. 2012;50:127â€“34.
Docherty PD, Chase JG, Lotz TF, Hann CE, Shaw GM, Berkeley JE. Independent cohort crossvalidation of the realtime DISTq estimation of insulin sensitivity. Comput Methods Programs Biomed. 2011;102:94â€“104.
Acknowledgements
Not applicable.
Funding
The study was supported by the FRIAFund for Research and Training in Industry and Agriculture (Belgium), the EUFP7 program, the NZ National Science Challenge 7, Science for Technology and Innovation (#CRSS32016), the MedTech Centre for Research Expertise (CoRE) funded by the New Zealand Tertiary Education Committee (#3705718), the Hungarian National Scientific Research Foundation, Grant No. K116574, and by the BMEBiotechnology FIKP grant of EMMI (BME FIKPBIO).
Author information
Authors and Affiliations
Contributions
All authors contributed to conception, design, and results interpretation of the presented study. VU carried out the main analysis and drafted the manuscript. SD had significant input in methodological design. JLD, JGC, TD, and BB had input in interpretation, manuscript formulation and into redaction process. GMS provided the initial clinical data and ongoing interpretation of the results. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
SPRINT was implemented as standard practice, and data audit and analysis were approved by the New Zealand Health and Disability Ethics Committee Upper South Regional Ethics Committee B (Ref: URB/07/15/EXP).
Consent for publication
Data are deidentified, and consent for its use and publication is approved by the New Zealand Health and Disability Ethics Committee Upper South Regional Ethics Committee B (Ref: URB/07/15/EXP).
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided 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 Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Uyttendaele, V., Knopp, J.L., Davidson, S. et al. 3D kerneldensity stochastic model for more personalized glycaemic control: development and insilico validation. BioMed Eng OnLine 18, 102 (2019). https://doi.org/10.1186/s1293801907208
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1293801907208
Keywords
 Glycaemic control
 Hyperglycaemia
 Blood glucose
 Insulin
 Insulin sensitivity
 Kernel density