Hypoglycemia Prevention via Personalized Glucose-Insulin Models Identified in Free-Living Conditions

Abstract

Background:

The objective of this research is to show the effectiveness of individualized hypoglycemia predictive alerts (IHPAs) based on patient-tailored glucose-insulin models (PTMs) for different subjects. Interpatient variability calls for PTMs that have been identified from data collected in free-living conditions during a one-month trial.

Methods:

A new impulse-response (IR) identification technique has been applied to free-living data in order to identify PTMs that are able to predict the future glucose trends and prevent hypoglycemia events. Impulse response has been applied to seven patients with type 1 diabetes (T1D) of the University of Amsterdam Medical Centre. Individualized hypoglycemia predictive alert has been designed for each patient thanks to the good prediction capabilities of PTMs.

Results:

The PTMs performance is evaluated in terms of index of fitting (FIT), coefficient of determination, and Pearson’s correlation coefficient with a population FIT of 63.74%. The IHPAs are evaluated on seven patients with T1D with the aim of predicting in advance (between 45 and 10 minutes) the unavoidable hypoglycemia events; these systems show better performance in terms of sensitivity, precision, and accuracy with respect to previously published results.

Conclusion:

The proposed work shows the successful results obtained applying the IR to an entire set of patients, participants of a one-month trial. Individualized hypoglycemia predictive alerts are evaluated in terms of hypoglycemia prevention: the use of a PTM allows to detect 84.67% of the hypoglycemia events occurred during a one-month trial on average with less than 0.4% of false alarms. The promising prediction capabilities of PTMs can be a key ingredient for new generations of individualized model predictive control for artificial pancreas.

Introduction

Patients with type 1 diabetes (T1D) have to face the problem of the optimal management of diabetes every day. In particular, hypoglycemia avoidance is one of the main barriers to the glucose control, so in the last years several alarm systems were designed to detect hypoglycemia risks and, eventually, prevent them.^1-5

The so-called artificial pancreas (AP)^6-8 accelerated the development of these alarm systems since the long trials in free-living conditions^9-13 needed for AP testing, and requires safety or alarm systems to prevent hypoglycemia events or at least to detect them in advance.

The alarm systems can be divided in two categories on the basis of the method used for the alarm generation: low-threshold detection and prediction. The first notifies the crossing of a critical blood glucose (BG) level, ¹⁴ while the second tries to foresee the hypoglycemia risk to allow the user to act in advance in order to avoid this event.^15-18 The latter typically requires the use of patient models to perform glucose trend predictions.

Different glucose-insulin models can be used in the alarm systems: for example, minimal, maximal, or black-box linear models. The first class is accepted both as clinical tool and as an approach to understand the composite effects of insulin on glucose tolerance^19,20; the second is able to better represent the significant interpatient variability that characterizes T1D population, eg, the one included in the UVA/Padova simulator^21,22 equipped with 100 vectors of model parameters, called “virtual patients.” The third class includes models typically identified from real-life data.

Recently, new identification techniques have been investigated^23-33; for a comprehensive literature review, we refer the reader to Zarkogianni et al ³⁴ and Oviedo et al ³⁵ Promising results have been obtained on this topic by our group.^36-41 The main goal of this paper is to use the algorithm described in Wang et al ³² and Soru et al ³⁶ to identify individualized models from free-living data for seven patients studied at the University of Amsterdam Medical Centre (AMC) and to validate individualized hypoglycemia predictive alerts (IHPAs) on a rather long period (one month) ⁹ for all these patients.

It is worth noting that the identification of reliable models on real data is particularly difficult. In fact, free-living conditions are much more challenging than the highly controlled experimental conditions of in-hospital studies, due to many confounding factors affecting BG in real life. Each patient can experience different problems during the trial and all possible factors of uncertainty have to be taken into account. Some examples are the physical exercise or differences in daily activities, human errors in patient-provided information or technical issues affecting the AP prototype adopted during the trial. All these aspects require adaptations of the identification techniques, originally developed for in silico data to deal with them.

The identification technique adopted in this paper is the impulse-response (IR) identification technique introduced in Soru et al ³⁶ and extended in Toffanin et al ⁴⁰ to cope with free-living data. Given the promising results obtained on a single patient, ⁴⁰ the technique has been applied here to seven patients with T1D studied at AMC to show its effectiveness in front of interpatient variability.

The performance achieved on real-life data by patient-tailored glucose-insulin models (PTMs) are reported for each single patient and as mean of the entire cohort. The quality of an IHPA to detect in advance hypoglycemia phenomena based on PTMs is evaluated in terms of sensitivity, accuracy, and false positive rate. The new alarm system shows better results in hypoglycemia detection in terms of both true and false positive with respect to the previously published results.

Given the prediction capabilities of PTMs used in the IHPAs, they have the potential to be used to synthesize a new generation of individualized model predictive control (MPC) for AP.

Methods

In this section, the IR approach ⁴⁰ is summarized together with the technique used to develop the IHPAs based on PTMs.

Impulse Response Identification

Fluctuations in glucose concentration can mainly be attributed to insulin infusions, $i\left(t\right)$ , and carbohydrate intake, $m\left(t\right)$ . Since these signals are continuous in time, the dynamics of glucose concentration in response to insulin injections and ingested meals is captured by a continuous-time linear model with the following structure ⁴⁰ :

$$\mathrm{CGM}\left(s\right)=G_{\mathrm{i}}\left(s\right)I\left(s\right)+G_{\mathrm{m}}\left(s\right)M\left(s\right)$$

where $\mathrm{CGM}\left(s\right)$ , $I\left(s\right)$ , and $M\left(s\right)$ are the Laplace transformations of the deviation of the glucose profile, the injected insulin, and the carbohydrate intake with respect to their basal values, respectively. The transfer functions $G_{\mathrm{i}}\left(s\right)$ and $G_{\mathrm{m}}\left(s\right)$ describe the insulin-glucose and meal-glucose interactions, respectively. The major issue in identifying these transfer functions is the input identifiability. Since the insulin boluses are usually very close in time to the ingested meals, these inputs can be considered as simultaneous, and their effects of glucose concentration are superimposed. Since these two different inputs cannot be temporarily separated without compromising the patient safety, identifying the single transfer function on a single individual can be problematic.

In order to overcome this limitation, the identification procedure is divided into two steps. In the first step, the average in silico patient (Av) of the UVA/Padova simulator ²² is used to identify the transfer functions parameters by performing experiments with a sufficient excitation to the overall system, and therefore, dangerous for a human subject. The first step provides the transfer functions $G_{\mathrm{i}}\left(s\right)$ and $G_{\mathrm{m}}\left(s\right)$ estimated for the average patient. The goal of the second step is to adjust their parameters to patients based on their data.

Step 1: linear average model

The goal of step 1 is the identification of the parameters of the transfer functions $G_{\mathrm{i}}\left(s\right)$ and $G_{\mathrm{m}}\left(s\right)$ with the following structure:

$$G_{\mathrm{i}}=\frac{{\mu }_{\mathrm{i}}}{\left(1+s{T}_{\mathrm{i}1}\right)\left(1+s{T}_{\mathrm{i}2}\right)\left(1+s{T}_{\mathrm{i}3}\right)\left(1+s{T}_{\mathrm{i}4}\right)}$$

$$G_{\mathrm{m}}=\frac{{\mu }_{\mathrm{m}}}{\left(1+s{T}_{\mathrm{m}1}\right)\left(1+s{T}_{\mathrm{m}2}\right)\left(1+s{T}_{\mathrm{m}3}\right)}$$

where the orders of $G_{\mathrm{i}}\left(s\right)$ and $G_{\mathrm{m}}\left(s\right)$ are designed by trial and error. The parameter vector to be identified is defined as follows:

$$\theta =\left[\begin{array}{c}{\theta }_{\mathrm{i}}\\ {\theta }_{\mathrm{m}}\end{array}\right],{\theta }_{\mathrm{i}}=\left[\begin{array}{c}\begin{array}{c}{\mu }_{\mathrm{i}}\\ T_{\mathrm{i}1}\end{array}\\ T_{\mathrm{i}2}\\ \begin{array}{c}{T}_{\mathrm{i}3}\\ T_{\mathrm{i}4}\end{array}\end{array}\right],{\theta }_{\mathrm{m}}=\left[\begin{array}{c}{\mu }_{\mathrm{m}}\\ T_{\mathrm{m}1}\\ T_{\mathrm{m}2}\\ T_{33}\end{array}\right]$$

The value of $\theta$ is identified by using in silico data collected with impulse response experiments performed on the Av patient ( ${\theta }^{\mathrm{Av}}$ ). The two components of ${\theta }^{\mathrm{Av}}$ are identified by solving separately two optimization problems which provide an optimal value of ${\theta }_{\mathrm{i}}$ ( ${\theta }_{\mathrm{i}}^{\mathrm{Av}}$ ) and of ${\theta }_{\mathrm{m}}$ ( ${\theta }_{\mathrm{m}}^{\mathrm{Av}}$ ).

The optimization problem is based on the minimization of the sum of squared residuals (SSR) formulated as follows:

$${\theta }_{*}={\mathrm{argmin}}_{\theta }\left(\mathrm{SSR}\left(\theta \right)\right)={\mathrm{argmin}}_{\theta }\Vert \mathrm{CGM}\left(t\right)-\widehat{CGM}\left(t,\theta \right)\Vert$$

where $\mathrm{CGM}\left(t\right)$ contains the reference glucose measurements and $\widehat{\mathrm{CGM}}\left(t,\theta \right)$ the estimated glucose according to the considered model.

A stability constraint is introduced in the optimization problem to guarantee non-negative values of the time constants. The SSR function is nonlinear in the parameters, so a nonlinear least-squared method is applied. A critical point is the algorithm initialization since the initial conditions can highly affect the final result. In this paper, a reasonable initialization is required for the gains and the slow time constants of the transfer functions, the other unknown parameters are initialized to 1.

The gains are set equal to the area under the curves (AUC) of the impulse response data, while the slow time constant (T₁) is estimated from the final part of these data, $t\in \left[t_{\mathrm{start}},t_{\mathrm{end}}\right].$ In particular, in this time interval, the output of the model can be expressed as an impulse response of a first-order dynamic system as follows:

$$\mathrm{CGM}\left(t\right)=\frac{\mu }{T_1}{\mathrm{e}}^{-t/T_1}$$

Applying a log transformation to the signals, the model can be rewritten as

$$Y=\Phi \theta +V$$

where

$$\begin{array}{l}\theta =\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right]=\left[\begin{array}{c}\ln \left(\frac{\mu }{T_1}\right)\\ -\frac{t}{T_1}\end{array}\right],Y=\left[\begin{array}{c}\ln \left(\mathrm{CGM}\left(t_{\mathrm{start}}\right)\right)\\ ⋮\\ \ln \left(\mathrm{CGM}\left(t_{\mathrm{end}}\right)\right)\end{array}\right],\\ \Phi =\left[\begin{array}{cc}1& t_{\mathrm{start}}\\ ⋮& ⋮\\ 1& t_{\mathrm{end}}\end{array}\right]\end{array}$$

and $V$ is the measurement error vector.

The linear minimal mean squared error (MMSE) estimator is used to identify the time constant $T_1$ :

$${\theta }^{\mathrm{MMSE}}={\mathrm{argmin}}_{\theta }{\left(Y-\Phi \theta \right)}^T\left(Y-\Phi \theta \right)$$

and the solution can be expressed in closed form without compromising stability constraints:

$${\theta }^{\mathrm{MMSE}}=\left({\Phi }^T{\Phi }^{-1}\right){\Phi }^{T}Y.$$

with $T_1=-1/{\theta }_2^{\mathrm{MMSE}}$ .

Step 2: linear individualized model

The first step provides the model parameters that describe the glucose dynamic of the Av patient, which is the starting point for step 2. The second step exploits the patient data to personalize the identified Av model. This step is mainly based on the hypothesis that the acquired knowledge on glucose dynamic is mostly reliable. Specifically, the patient-tailored parameters are identified by minimizing the SSR with respect to the patient data.

A delay term ( ${\tau }_{\mathrm{I}}$ ) is added in the glucose-insulin transfer function $G_{\mathrm{i}}\left(s\right)$ , since it is well-known that the insulin response is affected by absorption delays. Hence, the glucose-insulin transfer function is reformulated as follows:

$$G_{\mathrm{i}}^{\mathrm{d}}\left(s\right)={\mathrm{e}}^{-{\tau }_{\mathrm{I}}s}\cdot G_{\mathrm{i}}\left(s\right).$$

In this paper, the parameter ${\tau }_{\mathrm{I}}$ =15 minutes is set constant on the basis of the results obtained in Toffanin et al. ⁴⁰

In order to provide a reliable initial condition for the optimization, the gains µ_i and µ_m are initialized using some available clinical parameters, ⁴⁰ such as correction factor (CF) and the carbohydrate-to-insulin ratio (CR), in particular:

$$\begin{array}{c}{\mu }_{\mathrm{i}}^{\mathrm{Cp}}=\mathrm{CF}\\ {\mu }_{\mathrm{m}}^{\mathrm{Cp}}=\frac{\mathrm{CF}}{\mathrm{CR}}\end{array},$$

while the time constants are initialized to the estimated time constants of the Av model:

$${\theta }_{\mathrm{i}}^{\mathrm{Cp}}=\left[\begin{array}{c}\begin{array}{c}{\mu }_{\mathrm{i}}^{\mathrm{Cp}}\\ T_{\mathrm{i}1}^{\mathrm{Av}}\\ \begin{array}{c}{T}_{\mathrm{i}2}^{\mathrm{Av}}\\ T_{\mathrm{i}3}^{\mathrm{Av}}\end{array}\end{array}\\ T_{\mathrm{i}4}^{\mathrm{Av}}\end{array}\right],{\theta }_{\mathrm{m}}^{\mathrm{Cp}}=\left[\begin{array}{c}{\mu }_{\mathrm{m}}^{\mathrm{Cp}}\\ T_{\mathrm{m}1}^{\mathrm{Av}}\\ T_{\mathrm{m}2}^{\mathrm{Av}}\\ T_{\mathrm{m}3}^{\mathrm{Av}}\end{array}\right].$$

Then, the optimization problem is split in two parts to identify separately the effect of insulin and meals on glycemia. This choice is due to the higher variability of the insulin parameters with respect to meal ones. So, the parameter ${\theta }_{\mathrm{i}}$ is estimated without involving ${\theta }_{\mathrm{m}}$ in the optimization problem, specifically by solving

$${\theta }_{\mathrm{i}1}^{*}={\mathrm{argmin}}_{{\theta }_{\mathrm{i}}}\left(\mathrm{SSR}\left({\theta }_{\mathrm{i}}\right)\right)$$

with ${\theta }_{\mathrm{i}}^{\mathrm{init}1}={\theta }_{\mathrm{i}}^{\mathrm{Cp}}$ and ${\theta }_{\mathrm{m}}^{\mathrm{init}1}={\theta }_{\mathrm{m}}^{\mathrm{Cp}}$ .

Then, the estimated ${\theta }_{\mathrm{i}1}^{*}$ is exploited to initialize ${\theta }_{\mathrm{i}},$ ie, ${\theta }_{\mathrm{i}}^{\mathrm{init}2}={\theta }_{\mathrm{i}1}^{*},$ and a second optimization problem is solved to obtain ${\theta }_{\mathrm{i}}$ and ${\theta }_{\mathrm{m}}:$

$${\theta }_2^{*}={\mathrm{argmin}}_{\theta }\left(\mathrm{SSR}\left(\theta \right)\right)$$

with ${\theta }_{\mathrm{i}}^{\mathrm{init}2}={\theta }_{\mathrm{i}1}^{*}\mathrm{and}{\theta }_{\mathrm{m}}^{\mathrm{init}2}={\theta }_{\mathrm{m}}^{\mathrm{Cp}}$ .

The optimal parameters of PTM will be ${\theta }_2^{*}=[{\theta }_{\mathrm{i}},{\theta }_{\mathrm{m}}].$

Discrete-time model

The continuous-time model identified in the previous section is discretized via zero-order hold method obtaining $G_{\mathrm{i}}\left(z\right)$ and $G_{\mathrm{m}}\left(z\right)$ with sample time $T_{\mathrm{s}}$ = 1 minute.

So, the final model has the following structure:

$$\mathrm{CGM}\left(z\right)=G_{\mathrm{i}}\left(z\right)I\left(z\right)+G_{\mathrm{m}}\left(z\right)M\left(z\right)+E\left(z\right)$$

and involves a stochastic part that describes the residual error, $E\left(z\right)$ . This error is modeled as an AR process of order $n$ :

$$e\left(k\right)=a_{1}e\left(k-1\right)+\dots +a_{n}e\left(k-n\right)+\epsilon \left(k\right)$$

with $\epsilon \left(k\right)$ being a zero-mean white noise with variance $\lambda$ . The parameters $a_1,\dots ,a_n$ and $\lambda$ are estimated from the data by minimizing the one-step ahead prediction. The complexity $n$ of the AR model is fixed to five by trial and error. If the glucose measurements sample time ( $T_{\mathrm{CGM}}$ ) is greater than $T_{\mathrm{s}},$ an interpolated version of the glucose data can be used in the procedure. Finally, the complete model can be resampled with $T_{\mathrm{s}}$ = $T_{\mathrm{CGM}}=$ 5 minutes.

Performance metrics for identification

The prediction performance is assessed comparing the prediction of the future glucose ( $\widehat{\mathrm{CGM}}$ ) with the measurements collected during the trial. Various prediction horizons (PHs) are considered, and with $\widehat{\mathrm{CGM}}(k|k-\mathrm{PH})$ denoting the PH-steps ahead prediction of a model, three performance indices have been considered: FIT, coefficient of determination (COD), and Pearson’s correlation coefficient ( $\rho$ ).

The FIT is defined as

$$\mathrm{FIT}\left(\mathrm{PH}\right)=100\left(1-\frac{\left|\right|\widehat{\mathrm{CGM}}(k|k-\mathrm{PH})-\mathrm{CGM}\left(k\right)||}{\left|\left|\mathrm{CGM}\left(k\right)-\overline{\mathrm{CGM}}\right|\right|}\right)$$

with the mean of the glucose signal referred as $\overline{\mathrm{CGM}}$ .

The COD is defined as

$$\mathrm{COD}\left(\mathrm{PH}\right)=100\left(1-\frac{\left|\right|\widehat{\mathrm{CGM}}(k|k-\mathrm{PH})-\mathrm{CGM}\left(k\right)|{|}^2}{{\Vert \mathrm{CGM}\left(k\right)-\overline{\mathrm{CGM}}\Vert }^2}\right).$$

Both FIT and COD are equal to 100% if and only if $\widehat{\mathrm{CGM}}\left(k\right)=\mathrm{CGM}\left(k\right)\forall k=1,\dots ,N$ , ie, the prediction is perfect, and smaller than 100% and possibly negative otherwise.

$\rho$ is defined as

$$\rho \left(PH\right)=\frac{\begin{array}{l}{\sum }_{k=\mathrm{PH}}^N\left(\mathrm{CGM}\left(k\right)-\overline{\mathrm{CGM}}\right)\\ \left(\widehat{\mathrm{CGM}}(k|k-\mathrm{PH})-\overline{\widehat{\mathrm{CGM}}}\right)\end{array}}{\left|\left|\mathrm{CGM}\left(k\right)-\overline{\mathrm{CGM}}\right|\left|\right|\right|\widehat{\mathrm{CGM}}(k|k-\mathrm{PH})-\overline{\widehat{\mathrm{CGM}}}||}$$

with $\overline{\widehat{\mathrm{CGM}}}$ being the sample mean of $\widehat{\mathrm{CGM}}$ .

All these metrics are functions of PH, so the primary outcome to evaluate the quality of the models has been defined as the average value of each metric ( $\overline{\mathrm{FIT}}$ , $\overline{\mathrm{COD}}$ , and $\overline{\rho }$ ) for the considered PHs. A secondary outcome is the average value of each metric ( $\overline{\overline{FIT}}$ , $\overline{\overline{COD}}$ , and $\overline{\overline{\rho }}$ ) for every PH over the patient cohort.

For each index, the mean (±standard deviation, SD) for normally distributed data or the median [25th-75th percentiles] for the non-Gaussian case has been computed over the patient cohort as cumulative index of the entire population.

Individualized hypoglycemia predictive alert

In this section, we introduce the algorithm for predictive hypoglycemia alarms: the main goal is to forewarn the patient of these potentially dangerous events ahead in time in order to prevent them. The IHPA algorithm has been introduced in Toffanin et al ⁴⁰ and its core is a PTM used for the glucose prediction.

The main steps performed by IHPA algorithm are reported here in order to understand this procedure in detail. At each sample time k*, IHPA uses the PTM to predict the trend of the future glycemia $\widehat{\mathrm{CGM}}\left(k^{*}+1\right)$ , . . ., $\widehat{\mathrm{CGM}}\left(k^{*}+\mathrm{PW}\right)$ (ie, with $\mathrm{PH}=\mathrm{PW}$ ) by assuming that no insulin is administrated, ie, setting $\mathrm{i}(k^{*})=\mathrm{i}(k^{*}+1)=\cdots =\mathrm{i}(k^{*}+\mathrm{PW})=0$ . The hypoglycemia alarm is issued if $\widehat{\mathrm{CGM}}$ values are below the threshold G_hypo in the prediction window (PW) for at least $N_{\mathrm{samp}}$ . The main idea is to issue an alarm only if the hypoglycemia event is unavoidable. In fact, it cannot be avoided only by suspending insulin, and the subject must take some rescue carbohydrates. This alarm system is designed to be used in conjunction with an AP, where the insulin suspension is automatically performed by the system to optimize the glucose profile. So, the main purpose of this alarm is to notify when the AP cannot avoid the hypoglycemia by itself by only modifying the insulin delivery. If data are missing in the previous hour for more than 20 minutes, insufficient data are available for a reliable prediction and the IHPA is preventively shut off.

The parameters that characterize the algorithm have been set in this work as follows: PW = 8 (ie, 40 minutes), G_hypo = 70 mg/dL, and $N_{\mathrm{samp}}$ = 2 (ie, 10 minutes).

Performance metrics for alarm system evaluation

In order to introduce the performance metrics, some indices have to be clearly defined:

• - Hypoglycemia event (HE): It starts when the glycemia of the patient falls below the threshold G_hypo and it ends when the glycemia remains above the threshold for more than t_rec minutes

• - True positive (TP): It indicates that an HE really occurred at time k_h and an alarm is activated by the algorithm in the detection window (DW) defined as [k_h − DW_s, k_h − DW_e] (see Figure 1(a)). Note that HE has to be notified at least DW_e minutes before its occurrence and not too far in the past (too far defined by DW_s), so that the alarm is clearly related to that event.

• - True negative (TN): It occurs at $k^{*}$ if no alarm is issued and no HE occurred in the window [ ${\mathrm{k}}^{*}$ , $k^{*}$ + (DW_s − DW_e)] (see Figure 1(b)).

• - False positive (FP): It occurs when an alarm is activated at k* and no HE occurred in the window [k*, k* + DW_s] (see Figure 1(c)). An episode is not counted as a FP if a meal or a hypo treatment is administrated in the window [k*, k* + DW_s]. In fact, the alarm has to inform the patient about an unavoidable hypoglycemia, where unavoidable means that without the administration of carbohydrate by the patient this HE cannot be avoided. If a carbohydrate treatment is administrated, the alarm is not easily assessable. Note also that “late” alarms, ie, alarms occurring after DW_e, are not associated with a FP event even if they do not count as TP.

• - False negative (FN): It indicates that a HE occurred without an alarm issued in the DW interval. An example is reported in Figure 1(d) where DW does not contain any alarm. In fact, the alarm system turns on the alarm too late, just after DW_e, resulting in a FN.

Figure 1.

Example of true positive (a), true negative (b), false positive (c), and false negative (d) of prediction by individualized hypoglycemia predictive alert. In each panel, the continuous glucose monitoring (CGM) measured before the considered time instant $k^{*}$ (magenta), the future real CGM measurements (dashed pink) and the glucose prediction obtained via the patient-tailored glucose-insulin models used in the individualized hypoglycemia predictive alert (dashed blue) are shown. At each sample time, the individualized hypoglycemia predictive alert state is reported (gray/red square); if present, hypoglycemia events with their detection windows are shown.

The parameters related to these definitions have been set in this work as follows ³⁶ : t_rec = 20 minutes, DW_s = 9 (ie, 45 minutes), and DW_e = 2 (ie, 10 minutes).

The performance metrics used to evaluate IHPA approach are as follows:

• - True positive rate (TPR), or sensitivity, and true negative rate (TNR), or specificity, to measure the proportion of positives and negatives that are correctly identified, respectively; these indices can be computed as follows:

$$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}},\mathrm{TNR}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}$$

• - Predictive value (PPV), or precision, and negative predictive value (NPV) to measure the proportion of positives and negatives that are correctly identified over all the positive or negative predictions:

$$\mathrm{PPV}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}},\mathrm{NPV}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FN}}$$

• - False positive rate (FPR) and false negative rate (FNR) to measure the proportion of positives and negatives that are wrongly identified, respectively:

$$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{TN}+\mathrm{FP}},\mathrm{FNR}=\frac{\mathrm{FN}}{\mathrm{TP}+\mathrm{FN}}$$

• - Accuracy (ACC) to measure the proportion of classifications, both positive and negative, that are correctly identified:

$$\mathrm{ACC}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}}$$

• - False omission rate (FOR) to measure the proportion of negatives wrongly predicted over all the negatives:

$$\mathrm{FOR}=\frac{\mathrm{FN}}{\mathrm{TN}+\mathrm{FN}}$$

• - F1 score (F1) to measure the harmonic average of the precision and sensitivity, F1 score reaches its best value at 1 (perfect precision and sensitivity) and worst at 0:

$$F1=\frac{2*\mathrm{TPR}*\mathrm{PPV}}{\mathrm{TPR}+\mathrm{PPV}}$$

For each index, the mean (±SD) for normally distributed data or the median [25th-75th percentiles] for the non-normally distributed data has been computed over the patient cohort as cumulative index of the entire population.

It should be noted that in our population, the time spent in hypoglycemia was 1.9% ± 1.1% (mean ± SD) ⁹ and this poses well-known problems on the interpretation of the metrics reported above. In fact, in this condition, TN is very large with any reasonable alarm system, and it saturates the metrics influenced by this quantity, eg, TNR will be close to 100%, thus making these metrics less informative than in the case of balanced datasets. For instance, in an unbalanced dataset, like the one considered here, any degenerate alarm/classification algorithm that detects always the most common class (not hypoglycemia in our case) scores a very high accuracy without a real good performance.

In view of this, in this work, we will focus on sensitivity (TPR) and precision (PPV), not depending on TN, even if all the performance metrics reported before are computed for completeness.

Results and Discussion

We identified a PTM for each AMC patient. The performance metrics of these models are evaluated on the entire trial excluding the training set intervals used for the identification. The results reported in Table 1 involves the average values of the different metrics ( $\overline{\mathrm{FIT}}$ , $\overline{\mathrm{COD}}$ , and $\overline{\rho }$ ) computed by considering the prediction horizons (PH) varying from 1 to 12 (ie, from 5 to 60 minutes).

Table 1.

Predictive Performance of Patient-Tailored Models Identified for Seven Patients.

Patient		$\overline{\mathrm{COD}}$	$\overline{\rho }$
Patient #1	58.72	76.95	0.89
Patient #2	60.56	79.55	0.9
Patient #3	74.42	91.35	0.96
Patient #4	67.92	86.23	0.93
Patient #5	70.77	88.18	0.94
Patient #6	45.68	58.49	0.83
Patient #7	68.12	86.37	0.93
Population mean (±SD)	63.74 (±9.67)	81.02 (±11.11)	0.91 (±0.04)

Abbreviation: SD, standard deviation.

All indices are normally distributed, so mean ± standard deviation is reported for the entire population.

From the results, it is possible to note that $\overline{\mathrm{FIT}}$ computed on the entire population remains above 60% (63.74%) with a $\overline{\mathrm{COD}}$ of 81.02%. The results are in line with those in Toffanin et al. ⁴⁰ One patient (Patient #6) performed worse than the others. However, considering the large changes of the patients’ habits and the time-varying nature of the system under study, all the results seem acceptable.

The results reported in Table 2 involve the population values of different metrics ( $\overline{\overline{FIT}}$ , $\overline{\overline{COD}}$ , and $\overline{\overline{\rho }}$ ) computed by considering each single PH varying from 1 to 12 over the patient cohort.

Table 2.

Predictive Performance of Patient-Tailored Models Identified for the Entire Population Across the Prediction Horizon.

PH	$\overline{\overline{FIT}}$	$\overline{\overline{COD}}$	$\overline{\overline{\rho }}$
1	96.83 (±0.84)	99.91 [99.86, 99.94]	1 [1,1]
2	92.47 (±1.69)	99.41 (±0.27)	1 [1,1]
3	86.84 (±2.93)	98.19 (±0.8)	0.99 [0.99, 0.99]
4	80.34 (±4.45)	95.97 (±1.81)	0.98 (±0.01)
5	73.45 (±6.21)	92.62 (±3.45)	0.96 (±0.02)
6	66.53 (±8.14)	88.23 (±5.79)	0.94 (±0.03)
7	59.79 (±10.21)	82.94 (±8.89)	0.92 (±0.04)
8	53.35 (±12.36)	76.93 (±12.71)	0.89 (±0.06)
9	47.24 (±14.52)	70.36 (±17.18)	0.86 (±0.07)
10	41.43 (±16.63)	63.32 (±22.19)	0.83 (±0.08)
11	35.92 (±18.65)	55.96 (±27.6)	0.8 (±0.1)
12	30.72 (±20.56)	48.38 (±33.29)	0.76 (±0.11)

Abbreviation: PH, prediction horizon.

The indices normally distributed are reported in terms of mean ± standard deviation and the others as median [25th-75th percentiles].

The performance decreases with the increase of the PH as expected. Note that these models have been identified to design IHPA systems that uses predictions till 45 minutes (PH = 9), so it is important that the performance indices related to this interval remain above 45% for the FIT and 70% in terms of COD. However, seen the application considered in this work, the performance indices of more interest are the ones related to the hypoglycemia prediction reported in Table 3.

Table 3.

Individualized Hypoglycemia Predictive Alerts Performance Indices Computed on the Selected Cohort of Seven Patients.

Metric	Pat. #1	Pat. #2	Pat. #3	Pat. #4	Pat. #5	Pat. #6	Pat. #7	Population
HE	34	13	29	14	20	16	11	19.57 (±8.73)
TPR	97.06	84.62	72.41	78.57	85.00	75.00	100.00	84.67 (±10.56)
PPV	76.74	31.43	48.84	52.38	44.74	9.52	26.19	41.41 (±21.54)
FPR	0.143	0.328	0.328	0.130	0.264	1.291	0.404	0.328 [0.173-0.385]
FNR	2.94	15.38	27.59	21.43	15.00	25.00	0.00	15.33 (±10.56)
F1	85.71	45.83	58.33	62.86	58.62	16.90	41.51	52.82 (±21.27)
TNR	93.12	97.71	94.76	98.42	96.82	96.15	97.70	96.38 (±1.88)
NPV	99.99	99.97	99.88	99.96	99.96	99.95	100.00	99.96 (±0.04)
ACC	99.84	99.65	99.56	99.83	99.70	98.67	99.60	99.65 [99.57-99.80]
FOR	0.014	0.027	0.119	0.039	0.038	0.046	0.000	0.041 (±0.038)

Abbreviations: ACC, accuracy; F1, F1 score; FNR, false negative rate; FOR, false omission rate; FPR, false positive rate; HE, hypoglycemia event; NPV, negative predictive value; PPV, predictive value; TNR, true negative rate; TPR, true positive rate.

The performance computed on the entire population is reported in terms of mean ± standard deviation if the data are normally distributed and in terms of median [25th-75th percentiles] otherwise.

The patients experienced about 20 HE on average. The performance indices of IHPAs on each patient and on the entire the population are reported in Table 3. The results include the sensitivity (TPR) of the predictor of 84.67% and precision (PPV) of 41.41% with FPR of 0.328% and FNR of 15.33%. F1 reaches 52.82% on average for all seven patients. Regarding the performance metrics affected by the unbalanced nature of the dataset, the specificity (TNR), accuracy (ACC), NPV (all >96%), and FOR (<0.042%) indicate alarms activated mostly when needed.

Conclusion

In this study, the technique described in Soru et al ³⁶ and extended in Toffanin et al ⁴⁰ to deal with free-living data collected without ad hoc clinical protocols has been used to identify patient-tailored models for seven AMC patients. The data used in this work were collected during a long (1 month) AP trial in free-living conditions ⁹ involving three clinical centers. Even if the number of the considered subjects is limited, the importance of the results obtained in this work is enforced by the large changes of the patient habits and the time-varying nature of the system under study, together with the long duration of the observed period. Patient-tailored glucose-insulin models show good predictive performance on the considered test set in terms of glucose prediction in free-living conditions with an average FIT of 63.74% on the entire considered population.

The PTMs obtained with this technique have been used to develop an individualized alarm system for each patient. The IHPAs based on PTMs showed a sensitivity (TPR) of 84.67% with FPR of 0.328% on the entire cohort of patients: on average, about 85% of the hypoglycemia events occurred during the trial have been detected in time to allow a rescue action with a negligible (<0.4%) number of false alarms.

Future works will be focused on the synthesis of individualized MPC algorithms based on these PTMs, but the testing of these controllers requires suitable clinical trials.