Large-scale Virtual Clinical Trials
of Closed-loop Treatments
for People with Type 1 Diabetes

Tobias K. S. Ritschel Asbjørn Thode Reenberg John Bagterp Jørgensen Department of Applied Mathematics and Computer Science,
Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark (e-mail: [email protected]).

Abstract

We propose a virtual clinical trial for assessing the safety and efficacy of closed-loop diabetes treatments prior to an actual clinical trial. Such virtual trials enable rapid and risk-free pretrial testing of algorithms, and they can be used to compare different treatment variations for large and diverse populations. The participants are represented by multiple mathematical models, consisting of stochastic differential equations, and we use Monte Carlo closed-loop simulations to compute detailed statistics of the closed-loop treatments. We implement the virtual clinical trial using high-performance software and hardware, and we present an example trial with two mathematical models of one million participants over 52 weeks (i.e., two million simulations), which can be completed in 2 h 9 min.

keywords:

Virtual clinical trials \sepDiabetes \sepStochastic modeling \sepHigh-performance computing

^†^†thanks: This work was partially funded by the IFD Grand Solution project ADAPT-T2D (9068-00056B).

1 Introduction

Clinical trials ensure the safety and efficacy of medical treatments, but they are also time-consuming and expensive. Furthermore, they might result in a negative outcome, e.g., that the proposed treatment is not safe. Therefore, prior to the actual clinical trial, it is important to 1) evaluate the potential treatment performance, 2) identify shortcomings and risks, and 3) assess the advantages of the treatment over alternative treatments. This is the purpose of virtual clinical trials which involve virtual participants that are represented by a mathematical model (often consisting of differential equations). By using high-performance computing (HPC) software and hardware, such virtual trials can involve large populations and long-term protocols, which allows for extensive and fast testing of different treatment variations.

In this paper, we consider virtual clinical trials of closed-loop diabetes treatment systems. These are also referred to as artificial pancreases (APs). Worldwide, one in ten adults suffer from diabetes, and according to the International Diabetes Federation (2021), it accounted for 9% of the 2021 global health expenditure (USD 966 billion). Specifically, we consider type 1 diabetes (T1D) where, due to autoimmune destruction of the $\beta$ -cells, the pancreas does not produce any insulin. People with T1D require daily treatment with exogenous insulin in order to avoid high blood glucose concentrations (hyperglycemia). Prolonged hyperglycemia can lead to a number of health complications and chronic conditions, e.g., chronic kidney disease, cardiovascular disease, and damage to the eyes and nerves. Conversely, the insulin treatment can lead to low blood glucose concentrations (hypoglycemia), which can result in acute complications such as loss of consciousness and seizures.

Clearly, given the risks associated with hyper- and hypoglycemia, it is not straightforward for people with T1D to manage their insulin treatment. Therefore, over the last few decades, there have been significant developments within AP systems which can decrease this burden (Lal et al., 2019). APs typically consist of 1) a continuous glucose monitor (CGM), 2) an insulin pump, and 3) a control algorithm implemented on a smartphone or a dedicated device. The control algorithm repeatedly computes an appropriate insulin flow rate based on measurements from the CGM device and communicates it to the insulin pump. There exist a variety of control algorithms for computing the insulin flow rate, e.g., based on heuristics, fuzzy logic, proportional-integral-derivative (PID) control (Huyett et al., 2015; Sejersen et al., 2021), and model predictive control (MPC) (Boiroux and Jørgensen, 2018; Chakrabarty et al., 2020). All of these algorithms contain algorithmic parameters which must be tuned based on simulation, i.e., based on a virtual clinical trial. As the human physiology and behavior vary significantly between people and over time, this is a nontrivial task. In spite of this, the tuning is typically based on short-term simulations of one or a few virtual participants who are only represented by a single mathematical model. In contrast, if large-scale long-term virtual clinical trials (involving multiple mathematical models) are used to identify candidate algorithms and algorithmic parameters, the chances of a successful real-world clinical trial increase significantly.

In this work, we develop an approach for performing large-scale long-term virtual clinical trials of AP algorithms (i.e., Monte Carlo closed-loop simulations). The approach is an extension of our recent work (Reenberg et al., 2022; Wahlgreen et al., 2021) which allows the virtual participants to be represented by multiple mathematical models. Specifically, we use 1) the model described by Hovorka et al. (2002) combined with a CGM model and 2) a modification of the model described by Dalla Man et al. (2007) and Colmegna et al. (2020). Furthermore, we extend both models to incorporate uncertainty about the physiology and the CGM measurements, i.e., we formulate the models as stochastic differential equations (describing the dynamics) and stochastic algebraic equations (describing the CGM measurements). We implement the approach using high-performance C code, and we use an HPC cluster (DTU Computing Center, 2021) to carry out the computations. Finally, we demonstrate the utility of the virtual clinical trial by comparing the performance of an AP algorithm for the two models. Specifically, we simulate one million participants over 52 weeks with both models (i.e., two million in total) which takes 2 h 9 min.

The remainder of this paper is structured as follows. In Section 2, we describe the virtual clinical trial, and in Section 3, we describe the two models used in this work. Next, we present the results of the virtual clinical trial in Section 4, and finally, we present conclusions in Section 5.

2 Virtual clinical trial

The virtual clinical trial involves a population of people with T1D, a protocol describing, e.g., the size and duration of meals, mathematical models, model parameter values, and AP algorithms. Additionally, meals and other activities can be incorrectly announced to the AP algorithm (as is often the case in reality). Furthermore, the mathematical models can be either deterministic (no uncertainty) or stochastic (uncertain dynamics and measurements). The uncertainty can, for instance, represent physiological phenomena that are not included in the model or unknown model parameters.

2.1 Virtual people

The virtual clinical trial contains one million virtual people. Each person is represented by the same pieces of information as real people, e.g., a unique ID, given name, family name, date of birth, place of birth, sex, height, and body weight. Additionally, each person is associated with a set of mathematical models and model parameter values.

2.2 Protocols

A protocol describes the participants’ activities during the virtual trial. Typical examples include the time and carbohydrate contents of meals as well as the duration and intensity of physical activity. Furthermore, each protocol has its own ID, and it contains IDs, time stamps (start and end time), and type and size for each activity.

In Section 4, we demonstrate the virtual clinical trial using a protocol based on a Northern European lifestyle with respect to 1) meal times, 2) work weeks, 3) vacation weeks, 4) public holidays, and 5) seasons. The year is divided into 4 seasons consisting of 13 weeks. In total, there are 8 weeks of vacation (2 of them represent public holidays). All weeks are either a standard week, an active week, or a vacation week. Similarly, all days are categorized as either a standard day, an active day, a day with a movie night, or a day with a late night. Each type of week contains a different combination of the days, and each season contains a different combination of the weeks. Table 1 gives an overview of the weeks and seasons as well as the carbohydrate contents of the different meals in the trial. During autumn, winter, and vacation weeks, the participants are less active and eat more. Furthermore, in addition to the meals consumed during a standard day, active days have an additional exercise session, days with a movie night have an additional snack in the evening, and days with a late night have two additional snacks in the evening.

Table 1: The seasons, weeks, and meal carbohydrate contents (Reenberg et al., 2022).

Compositions of the seasons
Season	Standard weeks	Active weeks	Vacation weeks
Winter	6	4	3
Spring	6	6	1
Summer	7	3	3
Autumn	9	3	1

Compositions of the weeks
Week type	Standard day	Active day	Movie nights	Late nights
Standard	4	1	1	1
Active	3	3	1	0
Vacation	5	0	0	2

Body weight-dependent meal carbohydrate contents
Meal size	Amount of carbohydrates	For a 70 kg person
Large meal	1.29 g CHO/kg	90 g CHO
Medium meal	0.86 g CHO/kg	60 g CHO
Small meal	0.57 g CHO/kg	40 g CHO
Snack	0.29 g CHO/kg	20 g CHO

Refer to caption — Figure 1: Schematic of the different types of days during autumn and winter (Reenberg et al., 2022). Top: A standard day. Second from the top: An active day. Third from the top: A day with a movie night. Bottom: A day with a late night. In the spring and summer, the snack is consumed between breakfast and lunch, and the dinner is a medium meal.

2.3 Database

The virtual participants, model parameters, protocols, and results from the virtual clinical trial are stored in a database. We use PostgreSQL which is an open-source database system. The database enables direct comparison of the performance of different AP systems on different populations, e.g., people with high body weights or low insulin sensitivity (provided that this sensitivity is a model parameter). Additionally, the database contains basic components for building protocols, e.g., the 4 types of days used in the protocol described in Section 2.2. The user can construct their own protocols based on these basic components. Finally, the database can be used together with a graphical user interface to select and visualize the results of the virtual clinical trial, show statistics of certain demographics and protocols, and add new elements (for instance, virtual people, mathematical models and parameters, or protocols).

3 Models

In this section, we describe the two models used in the virtual clinical trial. In Section 3.1, we describe the model developed by Hovorka et al. (2002) combined with a CGM model, and in Section 3.2, we describe the model developed by Dalla Man et al. (2007) and Colmegna et al. (2020) where the meal model is replaced with that proposed by Hovorka et al. (2002). In order to simplify the demonstration of the virtual clinical trial in Section 4, we do not include exercise in the models or in the protocol.

3.1 An extension of Hovorka’s model

The insulin subsystem is described by


$\displaystyle\dot{S}_{1}(t)$	$\displaystyle=u_{I}(t)-\frac{S_{1}(t)}{\tau_{S}},$	(1a)
$\displaystyle\dot{S}_{2}(t)$	$\displaystyle=\frac{S_{1}(t)}{\tau_{S}}-\frac{S_{2}(t)}{\tau_{S}},$	(1b)
$\displaystyle\dot{I}(t)$	$\displaystyle=\frac{1}{V_{I}}\frac{S_{2}(t)}{\tau_{S}}-k_{e}I(t),$	(1c)

where $S_{1}$ and $S_{2}$ [mU] constitute a two-compartment chain representing the absorption of insulin, $I$ [mU/L] is the insulin concentration in the plasma, $\tau_{S}$ [min] is the insulin absorption time constant, $V_{I}$ is the volume in which insulin is distributed, $k_{e}$ [1/min] is an elimination rate, and $u_{I}$ [mU/min] is the insulin infusion rate. The insulin action subsystem is described the three compartments


$\displaystyle\dot{x}_{1}(t)$	$\displaystyle=k_{b1}I(t)-k_{a1}x_{1}(t),$	(2a)
$\displaystyle\dot{x}_{2}(t)$	$\displaystyle=k_{b2}I(t)-k_{a2}x_{2}(t),$	(2b)
$\displaystyle\dot{x}_{3}(t)$	$\displaystyle=k_{b3}I(t)-k_{a3}x_{3}(t),$	(2c)

where $x_{i}$ [1/min] are the insulin effect on the glucose distribution ( $i=1$ ), the disposal of glucose ( $i=2$ ), and the endogenous glucose production ( $i=3$ ). Furthermore, $k_{bi}$ [(L/mU)/min²] are activation rates and $k_{ai}$ [1/min] are deactivation rates (for $i=1,2,3$ ). The meal subsystem is described by a two-compartment chain:


$\displaystyle\dot{D}_{1}(t)$	$\displaystyle=A_{G}D(t)-\frac{D_{1}(t)}{\tau_{D}},$	(3a)
$\displaystyle\dot{D}_{2}(t)$	$\displaystyle=\frac{D_{1}(t)}{\tau_{D}}-\frac{D_{2}(t)}{\tau_{D}}.$	(3b)

Here, $D$ [mmol/min] is the meal carbohydrate content (per minute), $D_{1}$ and $D_{2}$ [mmol] represent the meal absorption, $A_{G}$ [–] is the bioavailibility of the carbohydrates, and $\tau_{D}$ [min] is a time constant. The glucose subsystem is also described by two compartments, i.e.,


$\displaystyle\dot{Q}_{1}(t)$	$\displaystyle=\frac{D_{2}(t)}{\tau_{D}}-F_{01,c}(t)-F_{R}(t)-x_{1}(t)Q_{1}(t)$
	$\displaystyle\quad+k_{12}Q_{2}(t)+EGP_{0}(1-x_{3}(t))),$	(4a)
$\displaystyle\dot{Q}_{2}(t)$	$\displaystyle=x_{1}(t)Q_{1}(t)-k_{12}Q_{2}(t)-x_{2}Q_{2}(t),$	(4b)

where $Q_{1}$ and $Q_{2}$ [mmol] are the accessible and non-accessible glucose compartments, $k_{12}$ [1/min] is a transfer rate, $EGP_{0}$ [mmol] is the endogenous glucose production (extrapolated to an insulin concentration of zero). Furthermore,


$\displaystyle F_{01,c}(t)$	$\displaystyle=\begin{cases}F_{01}&G(t)\geq 4.5~{}\mathrm{mmol/L},\\ F_{01}G(t)/4.5&\text{otherwise},\end{cases}$	(5a)
$\displaystyle F_{R}(t)$	$\displaystyle=\begin{cases}0.003(G(t)-9)V_{G}&G(t)\geq 9~{}\mathrm{mmol/L},\\ 0&\text{otherwise},\end{cases}$	(5b)
$\displaystyle G(t)$	$\displaystyle=\frac{Q_{1}(t)}{V_{G}},$	(5c)

where $F_{01,c}$ and $F_{01}$ [mmol/min] are the corrected and nominal total non-insulin dependent glucose fluxes, $F_{R}$ [mmol/min] is the renal glucose clearance, $G$ [mmol/L] is the glucose plasma concentration, and $V_{G}$ [L] is the volume in which the glucose is distributed. Finally, the CGM subsystem is

\dot{G}_{I}(t)=\frac{G(t)}{\tau_{IG}}-\frac{G_{I}(t)}{\tau_{IG}}.

(6)

Here, $G_{I}$ [mmol/L] is the interstitial glucose concentration measured by the CGM and $\tau_{IG}$ [1/min] is a time constant.

3.2 A modification of the UVA/Padova model

In the UVA/Padova model, the glucose subsystem is described by


$\displaystyle\dot{G}_{p}$	$\displaystyle=EGP(t)+Ra_{m}(t)-U_{ii}(t)-E(t)$
	$\displaystyle\quad-k_{1}G_{p}(t)+k_{2}G_{t}(t),$	(7a)
$\displaystyle\dot{G}_{t}(t)$	$\displaystyle=-U_{id}(t)+k_{1}G_{p}(t)-k_{2}G_{t}(t),$	(7b)
$\displaystyle G(t)$	$\displaystyle=\frac{G_{p}(t)}{V_{g}},$	(7c)

where $G_{p}$ and $G_{t}$ [mg/kg] are the plasma glucose in rapidly and slowly equilibrating tissues, respectively, $EGP$ [(mg/kg)/min] is the endogenous glucose production, $Ra_{m}$ [(mg/kg)/min] is the glucose rate of appearance, $U_{ii}$ and $U_{id}$ [(mg/kg)/min] are the insulin-independent and insulin-dependent glucose utilization, $E$ [(mg/kg)/min] represents renal excretion, $k_{1}$ and $k_{2}$ [1/min] are rate parameters, $G$ [mg/dL] is the glucose concentration, and $V_{g}$ [dL/kg] is the glucose distribution volume. Next, the insulin subsystem is represented by


$\displaystyle\dot{I}_{\ell}(t)$	$\displaystyle=-(m_{1}+m_{3})I_{\ell}(t)+m_{2}I_{p}(t),$	(8a)
$\displaystyle\dot{I}_{p}(t)$	$\displaystyle=-(m_{2}+m_{4})I_{p}(t)+m_{1}I_{\ell}(t)+Ra_{Isc}(t),$	(8b)
$\displaystyle I(t)$	$\displaystyle=\frac{I_{p}(t)}{V_{i}},$	(8c)

where $I_{\ell}$ and $I_{p}$ [pmol/kg] are insulin in the plasma and liver, $m_{i}$ [1/min] for $i=1,\ldots,4$ are rate parameters, $Ra_{Isc}$ [(pmol/kg)/min] is the insulin rate of appearance, $I$ [pmol/L] is the plasma insulin concentration, and $V_{i}$ [L/kg] is the insulin distribution volume. Furthermore,


$\displaystyle m_{2}$	$\displaystyle=\frac{3CL}{5HE_{b}V_{i}BW},$	(9a)
$\displaystyle m_{3}$	$\displaystyle=\frac{HE_{b}m_{1}}{1-HE_{b}},$	(9b)
$\displaystyle m_{4}$	$\displaystyle=\frac{2CL}{5V_{i}BW}.$	(9c)

Here, $CL$ [L/min] is the insulin clearance, $HE_{b}$ [–] is the basal hepatic insulin extraction, and $BW$ [kg] is the body weight. The insulin-independent and insulin-dependent glucose utilization are given by


$\displaystyle U_{ii}(t)$	$\displaystyle=F_{cns},$	(10a)
$\displaystyle U_{id}(t)$	$\displaystyle=\frac{(V_{m0}+V_{mx}X(t))G_{t}(t)}{K_{m0}+G_{t}(t)},$	(10b)
$\displaystyle\dot{X}(t)$	$\displaystyle=-p_{2U}X(t)+p_{2U}(I(t)-I_{b}),$	(10c)

where $F_{cns}$ [(mg/kg)/min] is the glucose uptake of the erythrocytes and the brain, $V_{mx}$ [mg L/(kg pmol min)] and $K_{m0}$ [mg/kg] are parameters, $X$ [pmol/L] is the insulin concentration in the interstitial fluid, $p_{2U}$ [1/min] is the rate of the insulin action on the peripheral glucose utilization, $I_{b}$ [pmol/L] is the basal insulin plasma concentration, and $V_{m0}$ [(mg/kg)/min] is


$\displaystyle V_{m0}$	$\displaystyle=\frac{(EGP_{b}-F_{cns})(K_{m0}+G_{tb})}{G_{tb}},$	(11a)
$\displaystyle G_{tb}$	$\displaystyle=\frac{F_{cns}-EGP_{b}+k_{1}G_{pb}}{k_{2}},$	(11b)

where $EGP_{b}$ [(mg/kg)/min], $G_{tb}$ [mg/kg], and $G_{pb}$ [mg/kg] are the basal endogenous glucose production and the basal plasma glucose masses. The model of the oral glucose absorption is similar to that in (3):


$\displaystyle\dot{D}_{1}(t)$	$\displaystyle=A_{G}D(t)-\frac{D_{1}(t)}{\tau_{D}},$	(12a)
$\displaystyle\dot{D}_{2}(t)$	$\displaystyle=\frac{D_{1}(t)}{\tau_{D}}-\frac{D_{2}(t)}{\tau_{D}},$	(12b)
$\displaystyle Ra_{m}(t)$	$\displaystyle=\frac{D_{2}(t)}{BW\tau_{D}}.$	(12c)

Again, $D$ [mg/min] is the meal carbohydrate rate, $D_{1}$ and $D_{2}$ [mg] describe the meal absorption, $A_{G}$ [–] is the carbohydrate bioavailibility, and $\tau_{D}$ [min] is a time constant. The endogenous glucose production is


$\displaystyle EGP(t)$	$\displaystyle=\max\{0,EGP_{b}-k_{p2}(G_{p}(t)-G_{pb})$
	$\displaystyle\quad-k_{p3}(I_{d}(t)-I_{b})\},$	(13a)
$\displaystyle\dot{I}_{d}(t)$	$\displaystyle=-k_{i}(I_{d}(t)-I_{1}(t)),$	(13b)
$\displaystyle\dot{I}_{1}(t)$	$\displaystyle=-k_{i}(I_{1}(t)-I(t)),$	(13c)

where $k_{p2}$ [1/min] and $k_{p3}$ [mg L/(kg pmol min] are the liver glucose effectiveness and the amplitude of the insulin action of the liver, $I_{d}$ and $I_{1}$ [pmol/L] constitute a two-compartment delayed insulin signal chain, and $k_{i}$ [1/min] is a rate parameter. Next, the renal excretion is given by

E(t)=\max\{0,k_{e1}(G_{p}(t)-k_{e2})\}.

(14)

Here, $k_{e1}$ [1/min] and $k_{e2}$ [mg/kg] are the glomerular filtration rate and the renal glucose threshold. The subcutaneous insulin delivery is described by


$\displaystyle\dot{I}_{sc1}(t)$	$\displaystyle=-k_{d}I_{sc1}(t)+k_{a1}I_{sc1}(t)+\frac{u_{I}(t)}{BW},$	(15a)
$\displaystyle\dot{I}_{sc2}(t)$	$\displaystyle=k_{d}I_{sc1}(t)-k_{a2}I_{sc2}(t),$	(15b)
$\displaystyle Ra_{Isc}(t)$	$\displaystyle=k_{a1}I_{sc1}(t)+k_{a2}I_{sc2}(t),$	(15c)

where $I_{sc1}$ and $I_{sc2}$ [pmol/kg] are insulin in a non-monomeric and monomeric state, $k_{d}$ , $k_{a1}$ , and $k_{a2}$ [1/min] are rate parameters accounting for subcutaneous insulin kinetics, and $u_{I}$ [pmol/min] is the insulin infusion rate. Finally, the subcutaneous glucose concentration, $G_{sc}$ [mg/dL], is

\displaystyle\dot{G}_{sc}(t)=-k_{sc}G_{sc}(t)+k_{sc}G(t),

(16)

and $k_{sc}$ [1/min] is the inverse of a time constant.

3.3 General mathematical form and simulation

The two models described in this section are in the form


$\displaystyle x(t_{0})$	$\displaystyle=x_{0},$	(17a)
$\displaystyle dx(t)$	$\displaystyle=f(t,x(t),u(t),d(t),p)dt$
	$\displaystyle\quad+\sigma(t,x(t),u(t),d(t),p)dw(t),$	(17b)
$\displaystyle z(t)$	$\displaystyle=h(t,x(t),p),$	(17c)
$\displaystyle y(t_{k})$	$\displaystyle=g(t_{k},x(t_{k}),p)+v(t_{k}).$	(17d)

Here, $t$ is time, $t_{0}$ is the initial time, $x$ are the states (i.e., the physiological state), $x_{0}$ are the initial states, $u$ are manipulated inputs computed by the AP algorithm (e.g., the insulin infusion rate), $d$ are disturbance variables (e.g., the meal carbohydrate content), and $p$ are model parameters (specific to each person). The first term in (17) is the deterministic (drift) term, and the second term is the stochastic (diffusion) term. At time $t_{k}$ (e.g., every 5 min), the AP receives a CGM measurement, $y$ , which is corrupted by noise, $v$ . Furthermore, the outputs, $z$ , are used to evaluate the AP algorithm. The measurement noise is normally distributed, and $w$ is a standard Wiener process, i.e., $v(t_{k})\sim N(0,R)$ and $dw(t)\sim N(0,Idt)$ where $I$ is an identity matrix. In between measurements, the manipulated inputs and the disturbance variables are modeled as constant:


$\displaystyle u(t)$	$\displaystyle=u_{k},$	$\displaystyle t$	$\displaystyle\in[t_{k},t_{k+1}[,$	(18a)
$\displaystyle d(t)$	$\displaystyle=d_{k},$	$\displaystyle t$	$\displaystyle\in[t_{k},t_{k+1}[.$	(18b)

Finally, when the AP algorithm (represented by the functions $\kappa_{k}$ and $\lambda_{k}$ ) receives a measurement, it updates its internal state, $x_{k}^{c}$ , and computes the manipulated inputs based on the measurement, $y_{k}=y(t_{k})$ , the targets $\bar{u}_{k}$ and $\bar{y}_{k}$ , an estimate of the disturbance variables, $\hat{d}_{k}$ , and the hyperparameters $p_{\kappa}$ and $p_{\mu}$ :


$\displaystyle x_{k+1}^{c}$	$\displaystyle=\kappa_{k}(x_{k}^{c},y_{k},\bar{u}_{k},\bar{y}_{k},\hat{d}_{k},p_{\kappa}),$	(19a)
$\displaystyle u_{k}$	$\displaystyle=\lambda_{k}(x_{k}^{c},y_{k},\bar{u}_{k},\bar{y}_{k},\hat{d}_{k},p_{\mu}).$	(19b)

In Section 3.1 and 3.2, we described the models without uncertainty, i.e., with $\sigma=0$ . In Section 4, we use uncertain models where we 1) add a stochastic term to the plasma glucose compartments (4) and (7) with (equivalent) diffusion coefficients of 1.5 mmol/min^3/2 and 270.24/ $BW$ mg/(kg min^3/2), and 2) add measurement noise with a variance of 0.1 mmol/L $\approx$ 1.8 mg/dL.

We use the Euler-Maruyama method with a step size of 0.5 min to simulate each participant, and we only store statistics (and not individual simulations). We implement the virtual clinical trial using parallel high-performance C code for shared-memory architectures and two 2.9 GHz AMD EPYC 7542 32-core processors.

4 Virtual clinical trial results

In this section, we demonstrate how the virtual clinical trial can be used to assess the safety and efficacy of an AP algorithm. We use the protocol described in Section 2, and for each model, we include one million participants. In the analysis, we consider the first 4 weeks a titration period and disregard them. In the following, we refer to the model presented in Section 3.1 as model A and the one presented in Section 3.2 as model B. We choose the algorithmic hyperparameters based on deterministic simulations with model A. The total computation time is 2 h 9 min.

We demonstrate the virtual clinical trial using a simple algorithm based on physiological insight and concepts from proportional-integral-derivative (PID) control. It uses two integrators (I-controllers) to estimate the basal and bolus insulin requirements, and a PD-controller to mitigate smaller variations in the blood glucose. It also uses deadbands, error truncation, a hypoglycemia amplification factor, switching logic, as well as a few heuristics. Fig. 2 shows how the controller doses insulin during 4 days of the virtual clinical trial for a single participant.

We evaluate the algorithm using the performance measures and targets described by Holt et al. (2021). The measures are 1) average glucose, 2) glucose management index (GMI), 3) glucose variation (GV) computed as the coefficient of variation, 3) time above range (TAR), 4) time in range (TIR), and 5) time below range (TBR). All of these are based on CGM measurements. The 5 ranges used to compute TAR, TIR, and TBR are as follows (all values are in mmol/L, and the colors are used throughout this section). Level 2 hypoglycemia: $[0,3[$ (red). Level 1 hypoglycemia: $[3,3.9[$ (light red). Normoglycemia: $[3.9,10]$ (green). Level 1 hyperglycemia, $]10,13.9]$ (yellow). Level 2 hyperglycemia: $]13.9,\infty[$ (orange).

For model A, $\tau_{IG}$ has the same value for all participants. All other parameters are sampled from normal distributions based on the sample means and variances of the parameter values presented by Hovorka et al. (2002). For model B, we generate normal distributions based on the means and variances reported by Colmegna et al. (2020). For the remaining parameters, we use the means and variances of the values reported by Kovatchev et al. (2010) and Dalla Man et al. (2007). In both cases, we discard parameter sets that deviate more than one standard deviation from the mean or lead to a basal rate lower than 0.4 U/h (corresponding to a glucose concentration of 6 mmol/L). For model A, we also require that all time constants (including the inverses of rate parameters) are within one order of magnitude of the mean.

Fig. 3 shows the mean and worst-case (i.e., the participant with the lowest CGM measurement) time in ranges (TIRs) as well as a box plot of the distributions of the TIRs. The mean TIRs are remarkably similar for the two models. The TIRs of the worst-case participants can be used to identify shortcomings. Here, the two worst-case participants appear to be problematic for different reasons. The box plot gives a more thorough overview and illustrates that, for model B, the distributions of TIR and the times in level 1 and 2 hyperglycemia are more narrow. Fig. 4 shows the total daily doses (TDDs) of basal and bolus insulin. For model B, the participants get less basal insulin but almost the same amount of bolus insulin. Fig. 5 shows the cumulative distributions of the CGM measurements. It also indicates that the worst-case participants are qualitatively different. Finally, Table 2 shows that the AP performs better for model B. Presumably, the fact that the AP was tuned for model A is outweighed by the (seemingly) smaller variations between the participants for model B, which could be due to the parameter distributions.

Table 2: Comparison based on the glycemic targets described by Holt et al. (2021).

Quantity	Target	A	B
Average glucose	$<154$ mg/dL	76.87%	91.09%
GMI	$<7\%$	77.48%	91.46%
GV	$\leq 36\%$	91.37%	92.38%
TAR (level 2)	$<5\%$	85.80%	92.55%
TAR (level 1 and 2)	$<25\%$	68.92%	70.76%
TIR (normoglycemia)	$>70\%$	78.01%	82.79%
TBR (level 1 and 2)	$<4\%$	100.00%	100.00%
TBR (level 2)	$<1\%$	100.00%	99.98%
All targets		68.10%	70.66%

5 Conclusions

We present an approach for conducting large-scale long-term virtual clinical trials for evaluating APs prior to the actual clinical trials. The participants are represented by multiple mathematical models consisting of stochastic differential equations. We use Monte Carlo closed-loop simulations, implemented with HPC software and hardware, to compute detailed treatment statistics. Finally, we present results from a virtual clinical trial with one million participants (represented by two mathematical models) over 52 weeks, which is conducted in 2 h 9 min.

References

Boiroux and Jørgensen (2018) Boiroux, D. and Jørgensen, J.B. (2018). Nonlinear model predictive control and artificial pancreas technologies. In Proceedings of the 2018 IEEE Conference on Decision and Control (CDC), 284–290.
Chakrabarty et al. (2020) Chakrabarty, A., Healey, E., Shi, D., Zavitsanou, S., Doyle III, F.J., and Dassau, E. (2020). Embedded model predictive control for a wearable artificial pancreas. IEEE Transactions on Control Systems Technology, 28(6), 2600–2607.
Colmegna et al. (2020) Colmegna, P., Wang, K., Garcia-Tirado, J., and Breton, M.D. (2020). Mapping data to virtual patients in type 1 diabetes. Control Engineering Practice, 103, 104605.
Dalla Man et al. (2007) Dalla Man, C., Rizza, R.A., and Cobelli, C. (2007). Meal simulation model of the glucose-insulin system. IEEE Transactions on Biomedical Engineering, 54(10), 1740–1749.
DTU Computing Center (2021) DTU Computing Center (2021). DTU Computing Center resources. 10.48714/DTU.HPC.0001. URL https://doi.org/10.48714/DTU.HPC.0001. Technical University of Denmark.
Holt et al. (2021) Holt, R.I.G., DeVries, J.H., Hess-Fischl, A., Hirsch, I.B., Kirkman, M.S., Klupa, T., Ludwig, B., Nørgaard, K., Pettus, J., Renard, E., Skyler, J.S., Snoek, F.J., Weinstock, R.S., and Peters, A.L. (2021). The management of type 1 diabetes in adults. A consensus report by the American Diabetes Association (ADA) and the European Association for the Study of Diabetes (EASD). Diabetologia, 64, 2609–2652.
Hovorka et al. (2002) Hovorka, R., Shojaee-Moradie, F., Carroll, P.V., Chassin, L.J., Gowrie, I.J., Jackson, N.C., Tudor, R.S., Umpleby, A.M., and Jones, R.H. (2002). Partitioning glucose distribution/transport, disposal, and endogenous production during IVGTT. American Journal of Physiology-Endocrinology and Metabolism, 282, E992–E1007.
Huyett et al. (2015) Huyett, L.M., Dassau, E., Zisser, H.C., and Doyle III, F.J. (2015). Design and evaluation of a robust PID controller for a fully implantable artificial pancreas. Industrial & Engineering Chemistry Research, 54(42), 10311–10321.
International Diabetes Federation (2021) International Diabetes Federation (2021). IDF diabetes atlas, 10th edition. ISBN: 978-2-930229-98-0.
Kovatchev et al. (2010) Kovatchev, B.P., Breton, M.D., Cobelli, C., and Dalla Man, C. (2010). Method, system and computer simulation environment for testing of monitoring and control strategies in diabetes. The U.S. Patent and Trademark Office. U.S. Patent Application Publication No. US 2010/0179768 A1, July 15, 2010.
Lal et al. (2019) Lal, R.A., Ekhlaspour, L., Hood, K., and Buckingham, B. (2019). Realizing a closed-loop (artificial pancreas) system for the treatment of type 1 diabetes. Endocrine Reviews, 40(6), 1521–1546.
Reenberg et al. (2022) Reenberg, A.T., Ritschel, T.K.S., Dammann, B., and Jørgensen, J.B. (2022). High-performance uncertainty quantification in large-scale virtual clinical trials of closed-loop diabetes treatment. In Proceedings of the 2022 American Control Conference (ACC). Accepted.
Sejersen et al. (2021) Sejersen, M., Boiroux, D., Engell, S.E., Ritschel, T.K.S., Reenberg, A.T., and Jørgensen, J.B. (2021). Initial titration for people with type 1 diabetes using an artificial pancreas. IFAC PapersOnLine, 54(15), 484–489.
Wahlgreen et al. (2021) Wahlgreen, M.R., Reenberg, A.T., Nielsen, M.K., Rydahl, A., Ritschel, T.K.S., Dammann, B., and Jørgensen, J.B. (2021). A high-performance Monte Carlo simulation toolbox for uncertainty quantification of closed-loop systems. In Proceedings of the 60th IEEE Conference on Decision and Control (CDC), 6755–6761.

Large-scale Virtual Clinical Trials of Closed-loop Treatments for People with Type 1 Diabetes