The availability of a representative population of virtual patients, i.e., a population large enough to represent all relevant human patient behaviours, is a key enabler for the design of In Silico Clinical Trials (ISCTs), that is trials following simulation-based approaches for the safety and efficacy assessment of pharmacological treatments and biomedical devices. This involves the development of Virtual Physiological Human (VPH) models able to represent the whole phenotypes spectrum of the human physiology of interest. Usually, such models are open-loop models, i.e., their behaviour depends also on exogenous inputs (such as, e.g., pharmacological drugs). In this paper, we propose a methodology to convert an open-loop VPH model into a closed-loop model. As a case study, we apply our methodology to a state-of-the-art VPH model defining the human glucose regulation system of individuals with Type 1 Diabetes Mellitus (T1DM), and show how we generate a representative population of T1DM virtual patients.
The advance of biomedical engineering has promoted the design of new, revolutionary, biomedical devices. The purpose of such devices is to improve the quality of life in different kinds of patients by making the treatments they follow completely, or partially, automated. An In Silico Clinical Trial (ISCT), in opposition to the more traditional in vivo clinical trial, is used in the development of products such as drugs and biomedical devices. More specifically, a clinical trial is performed by means of computer simulations over a population of virtual patients (see, e.g., http://paeon.di.uniroma1.it).
In order to generate virtual patients, a Virtual Physiological Human (VPH) model is required. A VPH model is a mathematical mechanistic model of the (patho-) physiology of interest and reactions to drugs (Pharmaco-Kinetics/Dynamics, PKPD) in the form of a parametric system of, e.g., Ordinary Differential Equations (ODEs). Many of the parameters in the system usually model intersubject variability, which means that different assignments to these parameters determine the behaviour of different patients. Unfortunately, many assignment to the parameters can lead to model evolutions that are not coherent with the laws of biology. As a consequence, the population of virtual patients can not be built by picking assignments arbitrarily, but a smart search in the parameter space is needed.
Given a VPH model M , we can define the corresponding population of
Virtual Patients (VPs) as the set S = f j is an assignment to the
parameters in M and defines a behaviour that is coherent with the laws of biology}
[
Thanks to the use of in silico clinical trials, risky and costly experiments on humans can be reduced and postponed to when a deeper knowledge of the sideeffects of the product has been acquired (http://paeon.di.uniroma1.it). ISCTs basically entail the following steps. First, generate a set of VPs whose predictions are in agreement with (patho-) physiology, PKPD, and in vivo clinical data. Second, use the above in vivo validated VPs to assess, through simulation, safety and efficacy of a candidate drug or biomedical device. 1.2
In this paper, we consider, as a case study, patients with Type 1 Diabetes
Mellitus (T1DM), and generate a population of VPs by exploiting the VP
generator developed by the Model Checking Laboratory (MCLab) (http://mclab.di.
uniroma1.it) and originally presented in [
For our case study, we implemented the Medtronic VPH model of the
human glucose regulation system [
To overcome this obstruction, we set up a methodology which allows us to
use our VP generator tool [
The paper is organised as follows. After recalling the state of the art on VPH models and ISCTs (Section 2), in Section 3 we introduce the Medtronic VPH model of the human glucose regulation system and outline how we implemented it in Modelica. Then, in Section 4 we describe how we extended our VPH model to enable the computation of a representative population of T1DM VPs, by computing and providing proper input patterns for carbohydrates intake and insulin injections. Finally, in Section 5 we analyse the results of our computation and in Section 6 we draw conclusions. 2
A key enabler to carry out an ISCT is the availability of a suitable cohort of
VPs. Such a cohort must be complete, i.e., large enough to represent all relevant
human patient phenotypes. In hormonal regulatory systems (as is our case), the
whole spectrum of possible phenotypes can be quite large, since such systems
typically occur within a complex network of endocrinological, neurological, and
psychological factors (see, e.g., [
It is worth noting that achieving simultaneously completeness and minimality is presently out of reach. For this reason, usually, completeness is privileged. This way, an adequate in silico representation of all human patient phenotypes is guaranteed.
In the following we briefly review the state of the art on VPs generation methods.
VPH models are usually developed using medical knowledge from the
literature and from pathway databases such as KEGG (Kyoto Encyclopedia of
Genes and Genomes) or Reactome. Open formats, such as the Systems
Biology Markup Language (SBML), allow exchange of computational models across
different platforms and their integration within open-standard general-purpose
simulation ecosystems [
Unfortunately, physiology knowledge is often qualitative. As a result, the
model behaviour is hard to define and is not uniquely determined. To overcome
such an obstacle, many qualitative as well as quantitative approaches have been
devised. An overview is in [
Qualitative approaches, often referred to as logic based, discretise the values
of interest. For example, using a Boolean approach, a gene up-regulation
(respectively down-regulation) may be represented with a logical 1 (respectively 0).
Boolean models have been widely studied (see, e.g., [
Quantitative models typically focus on representing dynamics of chemical
reactions through ODEs with stoichiometric parameters and reaction rates
estimated from clinical data using model identification techniques, as, e.g., in [
In such a setting, two approaches are typically used, respectively based on optimisation and statistical techniques.
Both global and local optimisation-based approaches, relying on AI as well
as on numerical methods, have been widely studied. Global optimisation–based
approaches are computationally heavy, however convergence is guaranteed to a
global optimum (see, e.g., [
Typically, both global as well as local optimisation approaches rely on model
identifiability, i.e., different values for the model parameters lead to different
model behaviours (see, e.g., [
When clinical data are scarce, e.g., when clinical assays are costly, risky and
invasive, identification approaches can be applied by averaging data coming from
different patients. However, in this case the resulting parameter value only yields
an inter-patient model behaviour (see, e.g., [
Since the goal of the above mentioned approaches is to generate a model parameter value that fits available experimental data, a huge amount of data per patient is needed in order to generate a virtual population that is representative of all human phenotypes. As a result, generating a complete set of VPs using model identification techniques would require considering a large amount of patients (ideally at least one for each relevant phenotype) along with a huge amount of clinical data for each of such patients. Unfortunately, this is exactly one of the obstacles ISCTs aim at overcoming. Thus, while model identification techniques are essential to validate models against clinical data and to develop patient-specific models, they can be hardly used to generate a complete set of VPs.
Moreover, VPH models are often globally or partially unidentifiable, i.e.,
wide ranges of parameter values lead to very similar model behaviours (see, e.g.,
[
Statistical approaches are also widely used. In such approaches a
parameter probability distribution, rather than a single value, is inferred (see, e.g.,
[
Limited knowledge about VPH model parameters calls for methods that handle models whose parameters are partially unknown. This guarantees completeness, possibly at the price of minimality. In other words, we use an overapproximation of the set of all physiologically admissible (i.e., whose behaviour is plausible with respect to physiology) VPs.
Since model-checking techniques enable exploration of all possible behaviours
of a model, it is not surprising that such approaches have been investigated in
order to generate VPs for both qualitative as well as quantitative VPH models.
In particular, for qualitative models (e.g., Boolean models) the problem of
finding parameter values yielding a biological meaningful behaviour is reduced to the
problem of finding stable states of the model, i.e., attractors. In such a setting,
symbolic model checking is widely used (see, e.g., [
Unfortunately, even using model approximation techniques like those in [
VPH models of human glucose regulation give a quantitative description of
pancreatic hormones-glucose interactions. These models can be used to predict the
evolution of the species of interest (e.g., blood glucose concentration) under
certain conditions (e.g., the amount of ingested carbohydrates), and are usually
validated by comparing the results with data collected during clinical trials. In
this paper, we use the Medtronic VPH model described in [
We implemented the Medtronic VPH model using the Modelica language, by
organising it in 4 compartments: glucose, insulin action, gastrointestinal
absorption, subcutaneous insulin absorption. Also, we validated our Modelica
implementation by comparing our simulation outcomes with the results shown in [
In order to study the behaviour of the patients, an individualisable VPH model
in an executable form is fundamental. Parametric ODE systems are usually
employed to define individualisable VPH models. Changes in the parameter values
entail the definition of patients with different behaviours. One way to generate
VPs could be to instantiate the VPH model by choosing the parameter values at
random. Unfortunately, this strategy usually leads to modelling the behaviour of
a patient which is incompatible with respect to the physiological knowledge, i.e.,
is not Biologically Admissible (BA). Behaviours that are not BA are typically
obtained by choosing parameter values ignoring inter-dependency constraints
among them. This might happen since such constraints are usually not
explicitly known [
In order to compute a population of patients representative of the real
population, we need to define a global search strategy that filters out all the virtual
subjects with a non-BA behaviour. The criterion proposed in [
In the following, we describe the tool we used to generate the population of patients with T1DM (Section 4.1), and present a new methodology that applies in the case of VPH models that require inputs (Section 4.2).
Time (a) Glucose
Time (b) Insulin
Time
(c) Glucose plasma appearance
We compute the population of T1DM patients using the tool presented in [
In order to identify the elements of S, the tool iteratively selects a random parameter (i.e., the model parameter vector) from the search space and, if passes the BA-test, adds to S.
Search effectiveness is achieved by sampling the parameter space using a (AIbased) strategy that conjugates exploration of the search space and exploitation of the knowledge acquired from the VPs already generated.
Search efficiency is achieved by using a HPC infrastructure in a
masterslave architecture, where a single master process (Orchestrator) is responsible
of intelligent sampling, storage of set S, and bookkeeping, and many slaves (BA
Verifiers) are devoted to simulation of candidate VPs and assessment of their
admissibility [
Assessment of biological admissibility (BA) of a parameter vector is
performed by checking whether the predicted time evolution obtained using is
qualitatively similar to the evolution obtained using the reference parameter
vector 0. The concept of similarity is quantified by the use of signal processing
Key Performance Indicators (KPIs) such as cross-correlation, normalised average
differences, and normalised squared norm differences [
The algorithm termination condition is based on statistical hypothesis
rejection techniques. Namely, given two user-defined values (confidence ratio) and
" (error bound ) in (0; 1), the algorithm, at any state of the search (when the
current population is stored in set S), formulates the following null hypothesis
H0(S): “the probability to find, by further sampling, a VP not in S is at least
"”. When the algorithm fails to find such a new VP within N = l lolgo(g1( )") m
consecutive attempts, it rejects the current null hypothesis H0(S), thus inferring
that the probability to find, by further sampling, a new VP is actually < ", and
terminates. Our algorithm implements a one-sided error procedure, in that it
might reject H0(S) when it holds. However, the probability of committing such
an error (which is called a type-II error ) is upper-bounded by [
In this section, we present a new methodology to exploit our VP generator, even in the presence of open VPH models, i.e., models that require external inputs. Then, we apply such a methodology to generate a population of VPs for the Medtronic VPH model.
To do so, we need to take into account the effect of inputs when performing the BA-test. First, we define a portfolio U of representative input patterns suitable for the open VPH model at hand, M . Then, for each input pattern u 2 U , we instantiate the given VPH model by feeding it with u. This leads to jU j instances of the model M , denoted mu, u 2 U . Last, we create a new VPH model, M 0, defined as the composition of such instances. Model M 0 is closed, i.e., no external inputs are required. In this way, by giving M 0 as input to the tool, the BA-test implicitly checks whether the behaviour entailed by a parameter vector is qualitatively similar to the behaviour entailed by 0, for all model instances mu; u 2 U .
Input Generators. The input to the model are the daily dose of insulin and the daily amount of ingested carbohydrates, together with meal times and timeof-maximum appearance rate of glucose. In order to enable the simulation of the behaviour of T1DM VPs, we define the two input generators, describing daily dose of injected insulin and the amount of daily ingested carbohydrates, respectively.
Insulin Dose Generator. In order to get analysable results, we need an element in the model with the capability to supply each patient with the recommended dose of insulin. According to the Diabetes Teaching Center at the University of California (https://dtc.ucsf.edu), the correct dose can be calculated with the following formula: dose = 0:55 body_weight insulin_factor where dose is measured in units of insulin per day (U/day) and body_weight is in kilograms. Therefore, we introduce the body weight in the set of patientspecific parameters. Moreover, to enable the study of behaviour of subjects who received too much or too little insulin, we add a multiplication factor, namely the insulin_factor, having its nominal value equal to 1.
Meals Generator. Our meals generator is based on the assumption that a subject with T1DM has to follow a balanced diet, or at least does not eat an amount of carbohydrates which significantly deviates from the recommended dose. For this reason, we feed the VPs with 5 meals a day, each one at a fixed time.
The time-of-maximum appearance rate of glucose ( m) characterising a meal depends on the complexity of the ingested carbohydrates, and can vary from 10 minutes to more than 2 hours. We set m to a different fixed value for each meal, and we add a multiplication factor in order to model inter-subject variability of time-of-maximum glucose absorption.
Finally, we define the percentage of carbohydrates in each meal as follows.
Our meals generator includes three different functions for calculating the recommended daily calories, and when invoked chooses the best (i.e., more precise) formula according to the available information about the patient. The simplest function in the generator has been proposed by Harvard Medical School (https://hms.harvard.edu) and takes in input only the body weight and the level of physical activity (e.g., sedentary, moderate, active) of the subject.
When further information, such as age and sex, are available, the generator computes the daily calories amount using the World Health Organization (http: //www.who.int/en) equation.
The most precise formula in our meals generator is known as the Mifflin St. Jeor equation and requires 5 inputs: body weight, level of physical activity, age, sex and height. The Mifflin - St. Jeor equation was judged by the American Dietetic Association as the most reliable formula in predicting actual energy expenditure.
Since subjects with diabetes are suggested to consume carbohydrates in the amount of 45–65% of the total calories intake (as estimated by the American Diabetes Association, http://www.diabetes.org) and 1 gram of carbohydrates provides 4 kCal, we define the daily amount of carbohydrates using the following equation:
carbo_factor 55 daily_calories daily_carbos = 400
In the equation above, carbo_factor is used to alter the amount of carbohydrates provided to subjects in order to simulate the behaviour of over- and under-eating patients.
Data. In order to provide a reference parameter vector value, 0, (i.e., a
reference BA VP) to our VP generator, we use the values of the subjects shown
in [
According to our methodology presented above, to perform a closed-loop simulation of the model we combine the VPH model together with our input generator. In particular, we define a portfolio of inputs, U , containing the following three elements: Input 1 The recommended amount of carbohydrates and the recommended insulin dose (with respect to the patient in exam).
Input 2 An excessive amount of carbohydrates and an insufficient insulin dose (with respect to the patient in exam).
Input 3 An insufficient amount of carbohydrates and an excessive insulin dose (with respect to the patient in exam).
In order to compute a population of VPs whose admissibility condition is based on qualitative similarity with respect to the reference patient under all our inputs (i.e., normal, hyperglycemic, and hypoglycemic conditions), we define an outer Modelica model which encapsulates three instances of the Medtronic VPH model, one for each input. The three Medtronic model instances where configured with insulin_factor and patient_carbo_factor set to, respectively: (1; 1) (normal conditions); (0:5; 3) (hyperglycemic conditions); (3; 0:5) (hypoglycemic conditions). Figure 2 shows our extended Medtronic VPH model.
In this section we show the outcome of the execution of our VP generator on the input described in Section 3. Simulations have been executed on 16 nodes, each with 16 cores, of partition A1 of the Marconi CINECA cluster (located in Bologna, Italy). We set both our confidence ratio and error bound " to 10 3. In other words, our algorithm terminates when, with statistical confidence at least 99:9%, the probability to find, via further sampling, a new VP is less than 0:1%. The overall computation took 5 days.
The importance of the obtained computed population resides in its size.
For example, the population computed using Subject 7 from [
Figure 3 shows the time evolution of blood glucose concentration (left) and
exogenous glucose appearance (right) in the three different scenarios (i.e.,
normal, hyperglycemic, and hypoglycemic conditions) in the VPs found by our
generator using Subject 7 as the initial reference. It can be seen that our computed
population covers a wide spectrum of behaviours, with glucose concentrations
assuming values in a much larger range than Subject 7 (physiological thresholds
are part of our tool configuration), although keeping a clear qualitative similarity
with respect to the time evolutions of the reference behaviour. On the other hand,
excursions of the exogenous glucose appearance in the population (right plots
of Figure 3) are way more limited. This is expected, as our input generator was
indeed explicitly designed to generate input patterns properly individualised on
the needs of each of the VPs (e.g., carbohydrates intake computation considers
information such as the body weight of each VP). Populations of VPs generated
starting from the other subjects of [
To make ISCTs effective, the first step to perform is the generation of a representative population of VPs, whose predictions are in agreement with (patho-) physiology, PKPD, and in vivo clinical data.
In this paper, we consider, as a case study, patients with T1DM, and
compute a representative population of VPs exploiting the Medtronic VPH model
of the glucose regulation mechanism [
Since the Medtronic VPH model requires a time profile for carbohydrates intake and insulin injection as inputs, we equipped it with a new module that produces such inputs in a principled-yet-individualised way, by defining both normal conditions (proper carbohydrates and proper insulin intake for the VP
Time (min) Blood glucose concentration
Time (min)
Exogenous glucose appearance
Time (min) Blood glucose concentration Time (min) Exogenous glucose appearance Hyperglycemic conditions
Time (min) Blood glucose concentration Time (min) Exogenous glucose appearance
at hand) and pathological conditions (hyperglycaemic and hypoglycaemic conditions). Our approach is general, as it can be employed to transform any openloop VPH model (i.e., a model requiring inputs to function properly) into a closed-loop model.
The computed population of VPs defines a starting point for ISCTs of new insulin-injection patterns as well as new medical devices (e.g., autonomous insulin pumps or artificial pancreas devices).
Acknowledgements. This work was partially supported by the following research projects/grants: Italian Ministry of University & Research (MIUR) grant “Dipartimenti di Eccellenza 2018–2022” (Dept. Computer Science, Sapienza Univ. of Rome); EC FP7 project PAEON (Model Driven Computation of Treatments for Infertility Related Endocrinological Diseases, 600773); INdAM “GNCS Project 2019”. The experimental part has been run on the Marconi CINECA cluster, thanks to Class C ISCRA Project n. HP10COBFWG.
The authors are grateful to Prof. Marianna Maranghi, MD (Dept. Translational and Precision Medicine, Sapienza Univ. of Rome) for useful discussions.