New model of glucose-insulin regulation characterizes effects of physical activity and facilitates personalized treatment evaluation in children and adults with type 1 diabetes

Accurate treatment adjustment to physical activity (PA) remains a challenging problem in type 1 diabetes (T1D) management. Exercise-driven effects on glucose metabolism depend strongly on duration and intensity of the activity, and are highly variable between patients. In-silico evaluation can support the development of improved treatment strategies, and can facilitate personalized treatment optimization. This requires models of the glucose-insulin system that capture relevant exercise-related processes. We developed a model of glucose-insulin regulation that describes changes in glucose metabolism for aerobic moderate- to high-intensity PA of short and prolonged duration. In particular, we incorporated the insulin-independent increase in glucose uptake and production, including glycogen depletion, and the prolonged rise in insulin sensitivity. The model further includes meal absorption and insulin kinetics, allowing simulation of everyday scenarios. The model accurately predicts glucose dynamics for varying PA scenarios in a range of independent validation data sets, and full-day simulations with PA of different timing, duration and intensity agree with clinical observations. We personalized the model on data from a multi-day free-living study of children with T1D by adjusting a small number of model parameters to each child. To assess the use of the personalized models for individual treatment evaluation, we compared subject-specific treatment options for PA management in replay simulations of the recorded data with altered meal, insulin and PA inputs.

. Glucose levels and GU and GP rates are shown. The PA session is marked by solid lines and onset of depletion by dashed lines. PA-driven changes in GU and GP are separated into insulin-dependent (id) and insulin-independent (ii) contributions. Table A. Model parameters. Healthy: healthy subjects. V1-V3: T1D subjects under euglycemia low-insulin (V1), euglycemia high-insulin (V2), and hyperglycemia low-insulin (V3) clamp conditions [3,4]. We used data by Romeres et al. [3,4] to re-calibrate the exercise model to people with T1D. GU and GP were studied during 60 min of exercise at 65% VO max 2 under three glucose and insulin conditions (V1: euglycemia -low insulin, V2: euglycemia -high insulin, V3: hyperglycemia -low insulin). When studying glucose uptake, Romeres et al. found that parameter p 2 and the glucose distribution volume V g were similar across all study conditions, while the participants' insulin sensitivity (SI, corresponding to p 3 ) and glucose effectiveness p 1 varied. Hence, we fixed V g and p 2 (p 2 determined from GIR required after insulin aspart injection [5]) as well as the exchange rates p 4 and p 5 (since glucose fluctuations are small), and only estimated parameters p 1 and p 3 of the glucose-insulin model. In addition, we estimated exercise parameters b, q 1 , q 2 , q 3l and q 4l as well as basal levels of glucose, G b , and insulin, I b . We fixed τ Z = 600 min, which defines the return of insulin sensitivity to its baseline, since it is difficult to estimate such a time scale correctly based on the comparatively short recovery period included in the data.
Our parameter estimation strategy entails two steps. First, we identified parameters for condition V1 from GU and GP data and used glucose and insulin measurements as known inputs. This condition represents the physiologically 'normal' scenario and hence, we used the resulting parameters as our baseline description of T1D populations for the remainder of this work. Second, we tested whether our model is applicable also to other glucose and insulin scenarios, and estimated parameters for conditions V2 and V3. We used the V1 parameter estimates as prior information to constrain the parameter estimation for these conditions, since they share the same study population.
Condition V1: euglycemia -low insulin We determined parameters for condition V1 using maximum likelihood estimation, minimizing the negative log-likelihood which corresponds to a Gaussian distribution of data points y i around model predictions f (t i , θ) at time t i with parameter set θ and residual variance σ 2 . For the minimization, we used the Nelder-Mead algorithm in Python's scipy.optimize package. Parameters p 1 and p 3 are not identifiable in the unconstrained optimization problem (see below). However, insulin-independent and -dependent contributions to GU are similar at rest in the original study and we exploited this fact by enforcing the relation p 1 = X b = p 3 /p 2 · I b during optimization, which makes the parameters identifiable. We also detected an outlier in GU at t = 190 min, which we removed to improve estimation.
Conditions V2: euglycemia -high insulin, and V3: hyperglycemia -low insulin We used the parameter estimates θ V 1 j from condition V1 to constrain the estimation problem for conditions V2 and V3. For this, we added a penalty term to the log-likelihood function that corresponds to a log-normal distribution of parameter values around the V 1 estimates with about 60% relative standard deviation. The penalized log-likelihood for the optimization is then We used the residual standard deviation σ V 1 from the previous estimate. We did not employ the previous constraint on p 1 and p 3 , which are freely estimated during this optimization.
The resulting glucose production and uptake rates are shown in Figure B. For all conditions, the model captures the exercise-driven changes in GU and GP well. Parameter values are reported in Table A and our results agree qualitatively with the assessment of exercise effects in the original study, also in terms of the differences observed between conditions.

Section 1.2.2 Parameter Identifiability
Next, we considered practical identifiability of parameters and computed profile likelihoods (PL) to obtain appropriate confidence intervals (CI) [6]. The profile likelihood is the log-likelihood function evaluated for a specific parameter θ j with values p and maximized over all remaining parameters: The (1 − α) confidence interval for θ j is then found from the quantiles of the χ 2 -distribution with one degree of freedom.
Based on the available data, some parameters are not identifiable for condition V1 (Fig. C (a)). Glucose and insulin levels are constant during the study period, impeding for example the separation of insulin-dependent and -independent processes already at rest. The parameter correlations ( Fig. C (b)) reveal more insights: p 1 and p 3 correlate negatively under condition V1 and can compensate for each other, confirming the observed identifiability issues. The exercise parameters q 1 and q 2 , and similarly q 3 and q 4 , also show strong correlations. This is expected, since the steady state values of the GU and GP rates can be maintained when the ratio of the two corresponding parameters is held constant. The profile likelihoods for conditions V2 and V3 (Fig. C (c),(d)) show more confined parameter ranges. This is due to changes in glucose or insulin concentration during the study, allowing to distinguish between insulin-dependent and -independent contributions, and was further improved by the penalty term used during optimization.

Section 4 Definition of Standard Patient
For our full-day simulation results, we defined a standard patient with parameter values taken over from the model calibration for T1D (condition V1). However, data did not include glycogen depletion and high intensity for this condition, and we therefore used depletion parameters β and q 6 from healthy subjects, assuming that the contribution of glycogenolysis to overall GP during PA is similar for healthy and T1D subjects. In condition V3, GP exceeds GU during high-intensity PA by a factor of 1.6 at steady state. We assumed the same ratio here to determine high-intensity parameters q 3h and q 4h .

Section 5 Replay Simulations Section 5.1 In-Silico Performance Assessment of Model Personalization
We tested the performance of the personalized model in replay simulations in-silico, re-creating the experiments of the original study [14]. We used the UVa/Padova simulator [15] (Python implementation [16]) to generate individual data of 10 virtual subjects. We then personalized our model on these data and performed replay simulations with altered meal and insulin bolus inputs. We followed Hughes et al. [14] and evaluated the replay performance using the mean absolute relative difference (MARD) between 'true' and replayed glucose trajectories. Specifically, each individual received a meal containing 0.7gCHO/kg at 12:00 in a 24h simulation using the UVa/Padova simulator. A meal bolus was injected at the same time, where the bolus dose was computed according to a personal insulin-to-carbohydrate ratio. We personalized the model on these data to each subject. We fixed the parameters at their T1D population values and only estimated p 1 , p 3 , G b , f and τ m , with τ m constrained between 10 and 150 min.
In the first experiment, we changed the meal size from 0 to 200% of the original meal size in steps of 20%. For each subject, we simulated 'true' data with the UVa/Padova simulator, and we simulated replay predictions according to our personalized models. The resulting MARD between 'true' and replayed data starting from meal time at 12:00 is shown in Fig. G (a).
To evaluate the impact of insulin inputs, we varied the insulin bolus from 50 to 150% of the original bolus in steps of 10% in a second experiment (Fig. G (b)).
In both simulation experiments, we found an acceptable MARD of less than 10% in most cases, but exceeded this threshold slightly for more extreme alterations. Note that in contrast to the replay simulations shown in Hughes et al. [14], we did not consider further deconvolution to feed unaccounted dynamics as additional inputs into the simulation. Rather, we relied solely on adjusting model parameters to individual subjects, and accepted slightly higher MARD values in return.

Characteristics of Study Participants
We selected five pediatric T1D participants for model personalization. Baseline characteristics are presented in Table C.

.2 Data Preparation
The recorded study data consist of glucose measurements, accelerometer counts, timing and dosing of insulin injections and timing and carbohydrate content of meals. Glucose data were recorded continuously by a CGM device. Continuous exercise data were obtained using an Actigraph model GT3X+ accelerometer (ActiGraph, LLC; Pensacola, Florida, USA) that is worn on the right hip. Acceleration is measured along three axes, and data of the vertical axis were used to quantify intensity of movement. In contrast, participants or their caregivers manually reported insulin injections and meals in logbooks. Discrepancies between the provided information and the measured glucose levels indicate partly inaccurate or incomplete logbook data. While incorrectly estimated meal sizes were compensated for by allowing arbitrary values of the bioavailability factor f in the meal model, we had to manually account for incorrect meal times and for missing meals. To this end, we shifted meal times and the corresponding insulin bolus times to adjacent glucose minima if required. On few occasions, we found small glucose peaks during the day without reported meals or other disturbances, and introduced additional snacks without insulin bolus in these cases. All meal adjustments are given in Table E. To extract PA intervals, we smoothed the AC data using a median filter with a window size of 15 min. We identified intervals exceeding 1500 counts/min for at least 10 min as PA sessions, according to the threshold for minimum PA intensity set in the model. We merged sessions that were separated by less than 10 min into a single PA period.
Finally, we needed to define the basal insulin infusion rate u b for each individual. For study participants on multiple daily injection (MDI) therapy, we assumed a constant basal rate, which we determined from the person's total daily insulin requirements and the contribution from basal injections. For participants on continuous subcutaneous insulin infusion (CSII) therapy, we either used the reported values for u b , or applied the same strategy as for MDI therapy if u b was not reported (Table C).

Section 5.2.3 Personalized Model
Parameter values of the personalized models for each of the five children are given in Tables D and E. The resulting model fits are shown in Figure 5 of the main text for two participants, and in Figure H for the remaining participants. Due to missing information, we were unable to fit a small number of glucose excursions either in the early morning or late evening; these are listed in Table F and were ignored for model personalization.