Physically-Informed AI Closed-Loop Sensor Glucose Forecasting Modeling Methodology with Application to Type 1 Diabetes

1. Introductionh

Type 1 diabetes (T1D) is a condition that renders the pancreas unable to produce beta cells that are needed to produce insulin to keep blood glucose concentrations (BGC) within a certain range to maintain glycemic homeostasis. High BGC (hyperglycemia) can lead to permanent damage to organs and body functions such as vision, limbs, etc., whereas low BGC (hypoglycemia) can cause seizures, brain damage, loss of consciousness, death, etc. [1,2]. Thus, maintaining healthy BGC is critical and quite challenging for people with type 1 diabetes.

The concept of an artificial pancreas (AP) is a closed-loop glucose control system with the following three main components: 1. a device and/or methodology to infer sensor glucose concentration (SGC), the controlled variable (CV), sufficiently, frequently, and accurately for 24-h online automatic feedback forecast control (AFFC); 2. a suitably precise, effective, and robust manipulated variable (MV) and delivery system (e.g., insulin pump); and 3. an intelligent system for fully 24-h feedback forecast control, diagnosis, and other critical AP tasks. Recently, progress in online knowledge-based (i.e., expert system) algorithms appears to have contributed to modest improvements in classification and aiding in reducing the frequency and magnitude of SGC spikes [3,4,5]. The focus of this work is the first component above. More specifically, the goal of this research is the development of an effective online SGC forecast modeling methodology for effective 24-h AFFC.

The typical behavior of a physical (i.e., nonbiological) automatic control system is that a change in MV will immediately impact its CV. However, when insulin (the MV) is injected into the body, it takes time for it to impact SGC. This time is the true deadtime for this insulin injection (II). The time it takes for SGC to decrease from an II we defined as the true observable deadtime (θ_MV). One of the major and unique AP control challenges is the existence of θ_MV. (Bi-hormonal modeling or control is not in the scope of this work. Our assumption is that advancements in SGC forecast modeling can contribute to advancements in bi-hormonal forecast modeling.) Successful (i.e., sufficiently tight) SGC requires accurate II. This objective requires sufficiently accurate cause-and-effect forecast modeling of SGC (the explanatory variable) a θ_MV distance in the future.

“Forecast” and “Prediction” are often treated as synonyms in casual conversations as well as in the mathematical modeling literature. However, in the context of this work, it is important to distinguish and define them to avoid subtle and critical misunderstandings. Prediction is a word that is commonly used in the statistical literature [6]. However, this use is applied in a static (i.e., nondynamic) context. More specifically, in the statistics literature, a prediction is the value a static (i.e., nondynamic) stochastic variable may become when sampled in the future regardless of how far into the future. This type of modeling (i.e., static) is not in the scope of this work.

In contrast, the value of a forecast (i.e., dynamic) stochastic variable is dependent on the forecast (i.e., time) distance. There are at least two types of forecasting scenarios. The first type, e.g., weather forecasting, seeks to predict future weather conditions but not change them. We call this “forecast monitoring” (FM). AFFC, our objective, seeks to accurately predict the future value of CV, and then effectively and safely manipulate conditions in the present, by changing MV, to drive CV to its target (i.e., set point, SP) after this future distance is reached. Thus, AFFC models SGC (the CV) to determine II (the MV) sufficiently accurately for acceptable SGC automatic control. We strongly believe that this model must be adequately “cause and effect” and developed and evaluated in outpatient studies. Our premise is that an FM methodology that models a correlation structure will fail in an outpatient, free living, AFFC application if it is not sufficiently “cause and effect”. All the SGC models that we have found in this literature model the correlation structure (e.g., see [30,31,32,33]). In contrast, the motivation of our work is to advance AFFC modeling of SGC that we feel will also be helpful in hormone modeling for AFFC for increasing SGC.

There are two ways to develop cause and effect models -- statistical design of experiments (SDOE) and theoretically-based modeling. We classify any experimental design that uses orthogonal (or near orthogonal) input changes for model identification as SDOE. In this context, SDOE does not seem like a practical alternative in general because a significant number of input changes will be highly undesirable and/or impractical (e.g., 100% fat meals). In addition, we feel it is critically important to develop SGC models from free living data to limit subject stress and model under the most comfortable conditions of the subject since studies can run for weeks to get adequate data size to build sufficiently acceptable models.

In AP terms, θ_MV is the forecast distance, i.e., the amount of future time it takes an insulin injection to start decreasing SGC. However, the forecast,$S{G}_{t+{\theta }_{MV}},$is the predicted value of SGC a θ_MV distance into the future. Thus, the AP feedback error (FBE) at the current time (t) is e_t, where e_t = SGC Setpoint -$S{G}_{t+{\theta }_{MV}}.$ It is highly critical that e_t is accurate for all t since it determines the amount of insulin that is injected into the user when the automatic control system is online, i.e., controlling insulin infusion. If the forecast of$S{G}_{t+{\theta }_{MV}}$is too high, this can lead to hypoglycemia (dangerously low SGC) or if it is too this can lead to hyperglycemia (dangerously high SGC).

Therefore, AP effectiveness is strongly dependent on the estimation accuracy of SG_t+θMV, the overall objective of this work. To achieve this objective, this work sets the following conditions (i.e., scope). The first one is that data collection is out-patient (non-hospital) and free-living (determined completely by the subject). Secondly, the training and validation periods are approximately one week each or more. The idea here is that there will likely be some subjects where both data sets include days that have similar and different activities (e.g., working days and non-working days). Thirdly, input data sets must not have missing data. A protocol must be defined and used to estimate missing data such as using the average of the value before the gap and after the gap to fill in the missing data in the gap. This is the method we used to fill in missing input data gaps. No missing input data means no missing forecast (i.e.,$S{G}_{t+{\theta }_{MV}},$ values. Missing input data results in large gaps in $S{G}_{t+{\theta }_{MV}},$values, and thus, in the ability of the controller to obtain the FBE and, thus, being online. Fourthly, the structure of $S{G}_{t+{\theta }_{MV}},$ must be sufficiently parsimonious (low parametrization) and physically interpretable to guard against “curve fitting” to maximize cause and effect behavior, since this is a control application and not a monitoring one. Fifthly, SGC measurements are used to “correct” $S{G}_{t+{\theta }_{MV}},$model bias must be at least a θ_MV distance in its past (i.e., at t, the current time, or less than t). Lastly, dynamic modeling must be physically-informed. Empirical modeling has a strong tendency to fit to correlation structures in data. Free-living data inherently has correlation structures, as humans are creatures of habit. Experimentally designed data sets with orthogonal input changes can produce cause-and-effect, empirical models. However, in this application, experimentally designed data is not realistic as mentioned above. Physically-informed modeling can provide cause-and-effect modeling from free-living data when the variables with physical interpretation are estimated accurately, and physical constraints are met. More specifically, now, the objective of this work is the development of a sufficiently accurate control application SG_t+θMV approach meeting these six criteria.

As mentioned above, in-patient (hospital clinical) or simulation studies are not in the scope of this work. However, they can provide valuable insights into the development of outpatient studies. In-patient and simulation studies commonly have durations of a few hours to a couple of days (see [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]).

There are several interesting outpatient modeling studies in the literature as well. For example, see [34,35,36,37,38,39,40,41,42,43,44]. However, we have not found any study in this literature meeting all six of our criteria. While there are a few studies meeting several of these criteria, no methodology, but our proposed approach, meets the sixth one, “dynamic modeling must be physically-informed.” All the ones that we have found are empirical modeling approaches.

We strongly advocate for the following three performance statistics for quantitative comparative methodological evaluation (the formulas for each will be given later). The measure of performance that is the least found in this literature is r_fit. This statistic quantifies how well the fitted model increases and decreases with SGC data, with an upper limit of r_fit = 1. Results that only show fit visually are considerably less informative than ones that give both the visual results and the quantitative results, which makes it possible to compare methods quantitatively and to set a numerical goal. For example, our r_fit validation goal is 0.9. We have found only a few studies in this literature reporting r_fit results [45,46,47].

The second type of statistic that we use measures spread, i.e., variability of the difference between the measured and fitted values. Statistics are typically given for at least one measure of spread in the SGC modeling literature. A popular one has been called the root mean squared error (RMSE) [48,49]. This statistic is a measure of spread, and more specifically, the standard deviation of the difference between measured SGC and fitted SGC. The drawbacks of this statistic include not being easily interpretable quantitatively since it is distributionally dependent, which can vary considerably across studies and the methodology used. Thus, similar values, quantitatively, can be very different inferentially based on the true probability distribution of the statistic. Notwithstanding, the measure of spread that we prefer to use is the average of the absolute differences (AAD) of measured and fitted values. AAD is easily interpretable across studies since it is an average of the magnitude of the differences, i.e., a measure of its centrality of absolute deviation. For example, for a data set with differences from its mean value of -30, 30, 30, -30, AAD is 30; in other words, the average distance from its mean value of zero is 30.

The final statistic that we advocate to be used is a measurement of model bias and is the average difference (AD) of measured and fitted SGC values. We have not seen this statistic reported much in the SGC modeling literature. However, its reporting is crucial to an overall assessment of a methodology. For example, if r_fit = 0.99 and AD is large positively or negatively, the model was excellent in following the increasing and decreasing behavior of the observed SGC but is quite biased. For the AAD example above, AD = 0. However, for the data set of -30, -30, -30, -30, AAD is still 30, but AD is -30, and, therefore, has a higher estimated modeling bias.

Thus, for each subject, our validation modeling goal is high r_fit,val (r_fit for validation), low AAD, and low AD. In practice, these values will be subject-specific. However, until we test our methodology in closed-loop studies, we have set a quantitative modeling goal for fit only. It is r_fit,val $\ge$ 0.90, which is quite ambitious for large, highly complex, outpatient, free-living SGC datasets.

As mentioned above, methodologies should not model an input with a modeled dead time less than θ_MV unless its value is known at t + θ_i where θ_i < θ_MV, e.g., for meal-announced carbohydrates. This is because a real time dynamic θ_MV forecast prediction model for a control application can only use current, or earlier than current, input data (excluding announcement inputs) as mentioned above.

Classical ANN in this context is a one-block, completely empirical modeling approach as illustrated in Figure 1a. All the inputs enter this block, and a mathematical structure with adjustable parameters (i.e., coefficients or “weights) is changed to obtain an acceptable agreement with measured response data. When the functions use lagged variables, it is an empirical dynamic model structure.

The work by Pun et al. is some of the earliest work on “Physically-Informed-Neural-Network (PINN)” that was found [50]. However, the work in [51] presented PINN with a two-block structure shown in Figure 2, aligning more clearly with the model representations in this work. As shown, the inputs enter a static ANN block (i.e., one with no lag variables), and its output variable enters a theoretically (i.e., physically) based block (i.e., model structure). We recognized Figure 2 as a Hammerstein Block-Oriented Structure [52] when its second block is a physically informed (i.e., interpretable) linear dynamic structure (G), that we have named H-PINN (as shown in Figure 1b). When physically interpretable linear dynamic blocks (G_i) are followed by a static ANN block, as in Figure 1c, this is a Wiener Block-Oriented Structure, which we have given the name W-PINN.

The aims of this study are to develop and evaluate two, 2-stage W-PINN methodologies, using the v_i’s of eleven (11) Type 1 diabetes models in a manuscript currently in review and posted at https://drollins9.wixsite.com/derrickrollins. Our arbitrary, but ambitious, goal is to achieve an r_fit,val $\ge$ 0.90 on at least one of the data sets using input data only (defined as Model 1). The first methodology estimates f(V) (see Figure 1c) using the ANN package in JMP^®. The second methodology estimates f(V) using a Python^® coded methodology that we developed.

2. Materials and Methods

The eleven free-living, Type 1, glucose data sets used in this work uses were first modeled in [53] and later by [54]. However, neither dynamically modeled these data sets for a closed-loop control application since θ_MV equaled zero in both. Thus, these studies are not applicable to AFFC, the application for this work.

[53] generated these data sets in 2011. Their SGC sampling rate, Δt, is 5 min, and all twelve (12) inputs have this sampling rate as required for discrete-time modeling. The inputs include lifestyle, activity, food consumption, stress, physiological changes, and time of day (i.e., the 24-h clock). When these data sets were obtained, the Medtronic glucose sensor required replacement every three to four days and had a much longer warm-up period than current devices. Due to variations in SGC missing data, the data set sample size (n) varied from 3,733 to 4,837. Subject-reported food logs for carbohydrates, fats, and proteins were used. We note that while the activity tracker was innovative for this study at the time, it did not have the most critical sensor, heart rate, and its technology is considerably less advanced than today’s wearable activity devices and, thus, obsolete by today’s standards. We combined bolus and basal insulin and reduced the number of inputs by one, to twelve. These data sets are available to the public on the website of the first author at https://drollins9.wixsite.com/derrickrollins. For additional details of the inputs and the data collection, see [53].

2.1. W-PINN

Our W-PINN approach uses backward difference derivatives (BDD) to discretize second-order-plus-dead-time-plus-lead (SOPDTPL) (see [55]) theoretical dynamic systems, the only type used in this work, as given in Eq. 1, below. The dynamic system does not have to be initially at a steady state for our W-PINN modeling methodology since this initial condition is also estimated. Note that, Eq. 1, below, is the expression for each of the p-inputs.

$$\begin{array}{l}{\tau }_i^2{\displaystyle \frac{d^2{v}_i(t)}{d{t}^2}}+2{\tau }_i{\zeta }_i{\displaystyle \frac{d{v}_i(t)}{dt}}+v_i(t)\\ \hspace{1em}\hspace{1em}={\tau }_{ai}{\displaystyle \frac{d{x}_i(t-{\theta }_i)}{dt}}+x_i(t-{\theta }_i)\end{array}$$

(1)

with

$$E[y(t)]=f(\mathrm{V}(t))$$

(2)

where $t\ge {\theta }_i\ge 0$and${\tau }_i\ge 0$for i = 1, …, p, x_i(t) is the value of the ith input variable at t, and v_i(t) is the value of the ith output variable at t, in the units of x_i, y(t) is the output variable in its units at t, $E[y(t)]$means the expected value (i.e., true mean) of y(t),$f(\mathrm{V}(t))$is the true output (gain) function of V(t), the vector of the v_i(t)’s. When$f(\mathrm{V}(t))$is a nonlinear function of V(t), as in the ANN (i.e., W-PINN) case, Eqs. 1 and 2, taken together, have a Wiener block-oriented structure [47], as shown in Figure 1c.

The lead term is the first term on the right side of the equal sign in Eq.1. This term tends to “speed up” the response and provides what the process modeling and control community has termed “numerator dynamics” [51,52]. [47] developed a second-order, multiple-input, single-output, discrete-time, nonlinear Wiener dynamic approach using BDD based on Eq. 1. More specifically, using BDD approximation applied to a sampling interval of Δt, an approximate discrete-time form of Eq. 1 is

$${\widehat{v}}_{i,t}=\left\{\begin{array}{l}{\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t-\Delta t}+{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t-2\Delta t}+\hspace{0.33em}{\widehat{\omega }}_{1,i}{x}_{i,t-{\widehat{\theta }}_i-\Delta t}\\ \hspace{1em}\hspace{0.33em}\hspace{0.33em}+{\widehat{\omega }}_{2,i}{x}_{i,t-{\widehat{\theta }}_i-2\Delta t},\hspace{1em}\hspace{1em}t>{\widehat{\theta }}_i=\widehat{m}{}_i\Delta t\\ {\widehat{v}}_{{\widehat{\theta }}_i}={\widehat{v}}_{{\widehat{m}}_i\Delta t}\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{1em}\hspace{0.33em}t={\widehat{\theta }}_i={\widehat{m}}_i\Delta t\\ \hspace{1em}\hspace{0.33em}is\hspace{0.33em}undefined,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}t<{\widehat{\theta }}_i={\widehat{m}}_i\Delta t\end{array}\right.$$

(3)

with

$${\widehat{\delta }}_{1,i}={\displaystyle \frac{2{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(4)

$${\widehat{\delta }}_{2,i}={\displaystyle \frac{-{\widehat{\tau }}_i^2}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(5)

$${\widehat{\omega }}_{1,i}={\displaystyle \frac{({\widehat{\tau }}_{ai}+\Delta t)\Delta t}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(6)

such that ${\widehat{\omega }}_{2,i}=1-{\widehat{\delta }}_{1,i}-{\widehat{\delta }}_{2,i}-{\widehat{\omega }}_{1,i}$to satisfy the unity gain constraint. From Eq. 3 with $t>{\widehat{\theta }}_i+2\Delta t,$

$$\begin{array}{l}{\widehat{v}}_{i,t}-{\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t-\Delta t}-{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t-2\Delta t}\\ \hspace{1em}\hspace{1em}={\widehat{\omega }}_{1,i}{x}_{i,t-\left({\widehat{\theta }}_i+\Delta t\right)}+{\widehat{\omega }}_{2,i}{x}_{i,t-\left({\widehat{\theta }}_i+2\Delta t\right)}\end{array}$$

(7)

$$\begin{array}{l}\Rightarrow \hspace{1em}\left(1-{\widehat{\delta }}_{1,i}B-{\widehat{\delta }}_{2,i}{B}^2\right){\widehat{v}}_{i,t}\\ \hspace{1em}\hspace{1em}=\left({\widehat{\omega }}_{1,i}{B}^{{\widehat{\theta }}_i+1}+{\widehat{\omega }}_{2,i}{B}^{{\widehat{\theta }}_i+2}\right)x_{i,t}\end{array}$$

(8)

$$\Rightarrow \hspace{1em}{G}_{i,t}={\displaystyle \frac{{\widehat{v}}_{i,t}}{x_{i,t}}}={\displaystyle \frac{{\widehat{\omega }}_{1,i}{B}^{{\widehat{\theta }}_i+1}+{\widehat{\omega }}_{2,i}{B}^{{\widehat{\theta }}_i+2}}{1-{\widehat{\delta }}_{1,i}B-{\widehat{\delta }}_{2,i}{B}^2}}$$

(9)

After obtaining${\widehat{v}}_{i,t}$for each input i, the modeled output value, at time t, is determined by entering these results into$f({ \mathrm{V} ^}_t)$, a static ANN in this application, i.e.,

$${\widehat{y}}_t=f\left({\widehat{V}}_t\right)$$

(10)

2.2. Stage 1 Modeling Method

As mentioned above, Stage 1 of this two-stage methodology is completed. This work is in review and posted at https://drollins9.wixsite.com/derrickrollins. We now give important Stage 1 details and results to aid in the understanding of the Stage 2 methodology. While missing output (i.e., SGC) measurements are acceptable, missing input values are not for discrete-time modeling. Activity tracker data were the only missing input data. These missing values were estimated by averaging the two values on both sides of a gap and filling in the gap with this value. Some gaps were several hours long. Cross-validation [57] was used to guard against overfitting, with the first week as the training (Tr) data set, and the second week as the validation (Val) data set.

In Stage 1, all inputs were first modeled separately on their own Excel^® worksheet with a first-order linear regression function. For each case, insulin was modeled first. The estimated deadtime, i.e., ${\widehat{\theta }}_{MV},$was set at 60 min and was varied one Dt forwards and backwards to find the value that gave the best fit. For all the inputs, the estimate of θ_MV, ${\widehat{\theta }}_{MV}$, was determined to be 60 min, i.e., 12 Dt. The food variables were the only ones with announcements, and the carbohydrate input was the only one found to have a deadtime less than ${\widehat{\theta }}_{MV}$. The “time of day” input has no deadtime, and the dead time for all other inputs was${\widehat{\theta }}_{MV}$, except for fats that had deadtimes that were much larger than${\widehat{\theta }}_{MV}$, as determined by model estimation. After determining the dynamic estimates for each input (i.e., Eqs. 4 to 6), these values were copied to an Excel^® worksheet as the dynamic structure starting values for fitting the SOPDTPL dynamic, and first-order static, multiple-input model (i.e., Eq. 10). With f(V) as a first-order multiple linear regression static function, the fixed estimate of V, $\widehat{V}$, for Stage 2 was determined.

2.3. Stage 2 Model Development Modeling Methods

The objective of this work is the development, evaluation, and comparison of two, Stage 2, W-PINN modeling approaches using ${ \mathrm{V} ^}_t$from Stage 1 to obtain$f({ \mathrm{V} ^}_t)$for each of the eleven cases for the two approaches. The first approach uses the JMP^® ANN toolbox to approximately find the smallest SSE (i.e., SSR) by fitting many cases and selecting the best one. The second approach is a newly developed, confidential, ANN structure that is coded using Python^®. Both approaches significantly improved fit over the first-order regression function$f({ \mathrm{V} ^}_t)$. A diagram illustrating the general ANN Stage 2 approach for obtaining${\widehat{y}}_t$ is given in Figure 1c.

2.4. Three Input Models

There are three types of input model structures we have developed for this application. The first one we call the “input only model” or “Model 1.” All the inputs in this structure have a deadtime$\ge \hspace{0.33em}{\widehat{\theta }}_{MV}$except for announcement inputs that can have deadtimes less than$\hspace{0.33em}{\widehat{\theta }}_{MV}$like carbohydrates, equal to$\hspace{0.33em}{\widehat{\theta }}_{MV}$like proteins, greater than$\hspace{0.33em}{\widehat{\theta }}_{MV}$like fats, and zero like time of day.

The second one we call the “input-output model” or “Model 2.” It combines the input-only structure of Model 1 (i.e., Eq. 10) with a model of weighted residuals, a minimum of $\hspace{0.33em}{\widehat{\theta }}_{MV}$distance in the past, as shown in Eq. 11 below (see [45] for the derivation).

$$\begin{array}{l}{\widehat{y}}_t=f\left({ \mathrm{V} ^}_{\mathrm{t}}\right)+{\widehat{\varphi }}_1\left(y_{t-{\widehat{\theta }}_{MV}}-{\widehat{y}}_{t-{\widehat{\theta }}_{MV}}\right)\\ \hspace{1em}\hspace{1em}+{\widehat{\varphi }}_2\left(y_{t-{\widehat{\theta }}_{MV}-\Delta t}-{\widehat{y}}_{t-{\widehat{\theta }}_{MV}-\Delta t}\right)+\dots \\ \hspace{1em}\hspace{1em}=f\left({ \mathrm{V} ^}_{\mathrm{t}}\right)+{\widehat{\varphi }}_1{e}_{t-{\widehat{\theta }}_{MV}}+{\widehat{\varphi }}_2{e}_{t-{\widehat{\theta }}_{MV}-\Delta t}+\dots \end{array}$$

(11)

where Eq. 11 has no value if any residual is not determinable due to missing output measurements. Thus, unlike Model 1, which has estimates for all t since it uses only input data, Model 2 will not have estimates when outputs are missing.

The final model, Model 1-2, is a combination of the strengths of Models 1 and 2. More specifically, for Model 1-2,

$${\widehat{y}}_t=\left\{\begin{array}{l}\mathrm{Model} 2 {\widehat{y}}_t, \mathrm{if} \mathrm{its} \mathrm{value} \mathrm{exists}\\ \mathrm{Model} 1 {\widehat{y}}_t, \mathrm{if} \mathrm{its} \mathrm{Model} 2 \mathrm{value} \mathrm{does} \mathrm{not} \mathrm{exists}\end{array}\right.$$

(12)

2.5. Forecast Structures

Equation 10,${\widehat{y}}_t=f\left({\widehat{V}}_t\right)$, is the model development version of Model 1 as it is used to estimate the output, y, at the current time, t. There are no missing input values in Eq. 10. This is why missing armband data had to be estimated. In addition, non-announcement input values to obtain Eq. 10 must be at least a distance of ${\widehat{\theta }}_{MV}$in the past. This requirement is because the model development input lag must be the same as the forecast input lag, i.e., $\widehat{V}$at t, for forecasting a${\widehat{\theta }}_{MV}$distance into the future.

After obtaining ${\widehat{V}}_t$, its transformation into the kΔt forecast, i.e., the online version, is given by Eq. 13 below:

$${\widehat{y}}_{t+k\Delta t}=f\left({\widehat{V}}_{t+k\Delta t}\right)$$

(13)

where, from Eq. 7, with t = t + kΔt,

$$\begin{array}{l}{\widehat{v}}_{i,t+k\Delta t}={\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t+(k-1)\Delta t}+{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t+(k-2)\Delta t}\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\widehat{\omega }}_{1,i}{x}_{i,t-\left({\widehat{\theta }}_i-(k-1)\Delta t\right)}+{\widehat{\omega }}_{2,i}{x}_{i,t-\left({\widehat{\theta }}_i-(k-2)\Delta t\right)}\end{array}$$

(14)

Note that, if k = 12 and Δt = 5 min, the Eq. 14 model is forecasting 60 min into the future. Thus, all the non-announcement inputs must have a model building and forecast prediction dead time of at least 60 min.

2.6. Statistical Analyses

The first, and most important, modeling statistic is r_fit (which is bounded between -1 and 1), the fitted correlation of the measured SGC,$y_i$, and the fitted SGC,${\widehat{y}}_i$, and given in Eq. 15 below.

$$r_{fit}=r_{y_t,{\widehat{y}}_t}={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left(y_i-\overline{y}\right)\left({\widehat{y}}_i-\overline{\widehat{y}}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}\cdot \sqrt{{\displaystyle \sum _{i=1}^n{\left({\widehat{y}}_i-\overline{\widehat{y}}\right)}^2}}}}$$

(15)

where n is the number of samples in the set and the bar above a statistic means that it is its sample mean value. The equations to determine AAD and AD are, respectively,

$$AAD={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left|y_i-{\widehat{y}}_i\right|}}{n}}$$

(16)

$$AD={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left(y_i-{\widehat{y}}_i\right)}}{n}}$$

(17)

The equation for SSE (i.e., SSR, the sum of squared residuals), the more common name and used by JMP^®, is

$$SSE=SSR={\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(18)

3. Results

Training and Validation Stage 1 and Stage 2 summary statistics for the eleven subjects are given in Table 1. All the results are r_fit unless indicated otherwise. Recalling that, each subject has a fixed${\widehat{V}}_t$that was determined in Stage 1 using the first-order static structure as given in Eq. 19, below:

$${\widehat{y}}_t={\widehat{a}}_0+{\widehat{a}}_1{\widehat{v}}_{1,t}+\dots +{\widehat{a}}_p{\widehat{v}}_{p,t}$$

(19)

As shown, Stage 1, Model 1, r_fit,val results varied from 0.59 to 0.77, with a mean of 0.68. Moreover, Stage 2 Model 1 r_fit,val results improved significantly over the Stage 1 results for both ANN approaches. As shown, JMP^® Stage 2 Model 1 r_fit,val results varied from 0.60 to 0.85, with a mean of 0.74. However, Python^® Stage 2 Model 1 r_fit,val results are significantly better than JMP^®, varying from 0.72 to 0.93, with a mean of 0.82. As a result, Model 2 training and validation results, and Models 1-2 validation results are given in Table 1 for Python^® only. From Model 1 to Model 2, the Python^® mean r_fit,val increased from 0.82 to 0.87, the minimum from 0.72 to 0.80, and the maximum of 0.93 did not change. In summary, Python^® Stage 2 results improved considerably over Stage 1 results and are significantly better than JMP^® Stage 2 results.

Graphical Python^® Stage 2 fitted and measured SGC results for Subject 2 (the best case) are given (i.e., plotted) in Figure 3, Figure 4 and Figure 5. Figure 3 is Model 1 training and validation. Figure 4 is Model 2 training and validation. Figure 5 plots are the combined validation results, with Model 1 plotted when there is no output data, and Model 2 is plotted when there is output data, i.e., the Model 1-2 validation plot. The Model 1-2 plots are associated with the results in the last three columns in Table 1. Figure 4 shows an excellent fit of Model 2 and the highly realistic behavior of Model 1 when Model 2 results are not possible because of missing SGC data (see Eq. 11).

4. Discussion

In this study, two physically-based virtual forecasting sensor approaches were developed for obtaining the value of the controlled variable (CV), SGC, a θ_MV time distance in the future, for a two-stage closed-loop process control application. The first stage, a physically (i.e., theoretical) based dynamic modeling approach [47] estimates the physically interpretable dynamic parameters from the measured inputs (x_i’s) with multiple physical constraints to obtain dynamic outputs (v_i’s). The v_i’s are the inputs to the second stage, a static ANN structure. For the first method, this structure was determined by using the ANN toolbox in JMP^®. For the second method, this structure was determined by using a confidential method that this work developed and coded using Python^®. Both methods resulted in large average improvements over the Stage 1 results using a first-order linear regression static structure (see Table 1). In addition, a critical advantage of these two approaches is that the modeling is much easier and much less time-consuming than the 2nd order multiple linear regression (MLR) approach. Thus, we strongly recommend ANN over MLR for the static model structure. Note that ANN modeling is just a particular class of multiple nonlinear regression. The MLR model was applied to Subject 11 to compare the performance with the ANN. The r_fit,val had a modest improvement from Stage 1 alone, going from 0.74 to 0.79, but much less than the 0.85 obtained for P-ANN.

We were pleasantly surprised by the P-ANN achievements of Models 1, 2, and 1-2. Model 1 has two subjects over the r_fit,val goal of 0.90 and a mean of 0.82. Model 2 has four subjects meeting or exceeding the r_fit,val goal and a significant increase in the mean r_fit,val of 0.87. The combined Model 1 and Model 2 approach, i.e., Models 1-2, had essentially the same summary statistics results as Model 2, as shown in Table 1. Thus, combining Models 1 and 2, to have continuous forecasting without missing fits did not adversely affect r_fit,val relative to Model 2, which had missing fits due to missing SGC measurements.

Insulin is the process variable that is changed to keep SGC close to its set point, i.e., it is the manipulated variable (MV). For the control system to do this well in a forecast feedback control scheme, the controlled variable, SG_t+θMV, must be accurately estimated. An empirical method could possibly control SG_t+θMV online accurately if the correlation structure remains the same as it was when the model was developed. However, it is not possible for the correlation structure of an empirical forecast modeling approach in this context to remain intact, i.e., fixed, in online forecast feedback control because the correlation structure changes each time the controller signal to the manipulated variable is transmitted. Thus, it is prudent to restrict free-living empirical modeling to monitoring open-loop processes but not to make decisions on how much to change a manipulated variable to make changes in the control variable.

During his time as a professor, the second author gained valuable insight into the limitations of empirical modeling through a real-world industrial application. A BS Chemical Engineering student, also pursuing an MS in Statistics, undertook a summer project at a leading Midwest chemical company, which was approved as the basis for her MS thesis. The project focused on developing a multivariate Statistical Process Control (SPC) monitoring methodology for a process line. Data were collected, and an SPC control chart was developed, resulting in an excellent model fit. However, when the process exceeded control limits, adjustments to the manipulated variable based on this model failed to restore control. A subsequent attempt with new data and a revised control chart, despite another excellent fit, similarly failed to correct deviations when applied in a feedback control scenario. This experience highlighted that the control chart, designed for monitoring, was unsuitable for feedback control due to its reliance on empirical correlation rather than cause-and-effect relationships. Empirical SGC modeling, which uses free-living data and non-physiological structures, faces similar limitations, as it cannot adequately capture cause-and-effect dynamics critical for model-based control applications like automatic forecast control. In contrast, physically-informed modeling, which integrates physiological information and structure with free-living data, offers inherent intelligence and robust structure for developing effective models for control applications. The W-PINN methodology proposed in this manuscript exemplifies such an approach, enabling cause-and-effect modeling suitable for closed-loop SGC control.

A limitation of the proposed two-stage W-PINN approach is the sequential estimation of parameters, with dynamic modeling parameters determined in Stage 1 and static modeling parameters in Stage 2. A more robust, yet considerably more complex, alternative is a one-stage approach that estimates all parameters simultaneously. Our research group is currently developing this method, with implementation in Python, aiming to complete the work and draft a manuscript within approximately one month. We anticipate that this one-stage approach will yield significant but moderate improvements in model performance compared to the two-stage W-PINN methodology.

A drawback of the proposed two-stage W-PINN approach is that the dynamic modeling parameters are estimated in Stage 1 and the static modeling parameters are estimated in Stage 2. A better but significantly more challenging approach is to estimate all modeling parameters in a One Stage approach. Our research group is working on this development and hopes to complete this work and write this manuscript in a month or so. It will be coded using Python and our expectation is a significant but modest improvement over our Two-Stage approach.

Two very popular empirical dynamic modeling approaches are Nonlinear Autoregressive Moving Average with eXogenous variables (NARMAX) [43,58,59] and Long and Short Term Memory (LSTM) [60]. Our research group evaluated NARMAX in [61] using real, freely-existing distillation column data. The ten (10) data sets, covering a period of three years, were generated by undergraduate chemical engineering students for their unit operation lab course. For this data set, [54] was not able to find adequate starting values to obtain a NARMAX lagged-based fit using Matlab^®. They transformed NARMAX into a physically-informed (N-PINN) structure and compared the results of the eight (8) test cases. The mean testing r_fit (i.e., r_fit,ts) for W-PINN and N-PINN, were 0.84 and 0.28, respectively. We plan to evaluate LSTM against our one-stage W-PINN approach using the diabetes data sets in this work and these distillation column data sets.

5. Conclusions

Type 1 diabetes SGC modeling for monitoring can be effective (i.e., informative) using empirical or physically-informed dynamic modeling approaches. Closed-loop Type 1 SGC automatic control is inherently forecast automatic control because a change in MV, injected insulin, will take a time of θ_MV to start lowering SGC. For automatic closed-loop control, empirical dynamic modeling approaches are not likely to succeed because they lack a cause-and-effect relationship, unlike Physical Informed Neural Network (PINN) approaches. The W-PINN approach, developed in this work, is particularly powerful because each input x_i is dynamically transformed to its vi and is the input to a static ANN. The proposed two-stage W-PINN approach greatly improved the SGC model fit for eleven historical diabetes data sets. A one-stage W-PINN approach, in its evaluation stage, is the next step in this research.

There are several drawbacks to the data sets used in this work. First, they are nearly a decade and a half old, particularly in terms of the amount of missing or lost data. Glucose sensor technology and activity tracker technology have improved considerably and particularly in terms of the amount of missing or lost data, also. In addition, there are advancements in ways to get accurate consumption of food nutrients. Thus, one future goal is to evaluate the methodology using modern technology.

Since this is a forecast control application modeling free-living data, to get a more informative evaluation of a methodology, we feel it is critical to collect a minimum of 3 weeks of free-living outpatient data (one week each for training, validation, and testing). However, we feel that four would be more ideal with two weeks of testing data.

The most important input is SGC. We feel that it is an open question as to whether it is the only needed input because it is the MV. Since several inputs, like the three (3) nutrients, have a big influence on the level of SGC but no role as contributors to the MV, we may discover that the number of critical inputs for effective automatic feedback forecast control (AFFC) is a very short list and fit performance will be relative to the maximum possible abilities of the critical inputs.