In-Silico Clinical Trials (ISCT), i.e., clinical experimental campaigns carried out by means of computer simulations, hold the promise to decrease time and cost for the safety and efficacy assessment of pharmacological treatments, reduce the need for animal and human testing, and enable precision medicine. In this paper we present a case study aiming at quantifying, by means of a multi-arm ISCT supervised by intelligent search, the potential impact of precision medicine approaches on a real pharmacological treatment, namely the downregulation phase of a complex clinical protocol for assisted reproduction.
Model-based approaches to safety and efficacy assessment of pharmacological treatments (In-Silico Clinical Trials, ISCT) hold the promise to decrease time and cost for the needed experimentations, reduce the need for animal and human testing, and enable precision medicine, where personalised treatments optimised for each patient can be designed before being actually administered. This is to be achieved by developing computational models for human physiology, pathophysiology, PharmacoKinetics (PK) and PharmacoDynamics (PD), which also define the possible physiological differences between different individuals (i.e., the possible phenotypes ).
Research in Virtual Physiological Human (VPH) provides mechanistic
quantitative models of the human physiology at levels of scale ranging from body
compartments (see, e.g., [
Independently of the level of scale used, most VPH models are defined by means of complex, highly non-linear differential equations. By means of numerical integration, such models can be executed within simulators in order to provide quantitative information about the time evolution of the modelled biological quantities (e.g., hormone concentrations).
Also, VPH models of different body compartments can be integrated in larger
models, eventually defining whole human models. Important attempts in this
direction are given by HumMod [
Current pharmacological treatments are often designed with the average patient in mind. A key topic in precision medicine is to develop pharmacological treatments optimised for any given individual (personalised treatments). Several optimisation criteria can be defined. A typical criterion is minimisation of the overall amount of administered drug, which often also reduces the probability of adverse side-effects and contributes to decreasing the treatment cost.
With their amenability to define different phenotypes, VPH models of proved accuracy are a key enabler for precision medicine. 1.2
In this paper we present a case study aiming at quantifying, by means of
extensive computer simulation–based experimental campaigns (In-Silico Clinical
Trials, ISCT) guided by intelligent search, the impact of precision medicine
approaches on a real pharmacological treatment, namely the downregulation phase
of a complex clinical protocol for assisted reproduction in humans. To do this,
we exploit a large VPH model of the Hypothalamic–Pituitary–Gonadal (HPG)
axis [
Here we select, from the population of VPs computed in [
Our multi-arm ISCT has been conducted on a large HPC infrastructure using a backtracking-based search algorithm to seek for the lightest but still effective treatment for each VP. Our search algorithm drives a (black-box) simulator of the VPH model at hand to intelligently explore the space of possible time-series of drug administrations.
We note that a clinical trial with a so high number of arms is not even conceivable in the classical in vivo setting (it would require a distinct human patient for each phenotype and for each possible alternative time-series of drug administrations). This undoubtedly shows the revolutionary potential of artificial intelligence for model-based (in silico) precision medicine. Paper outline. The paper is organised as follows. After describing some background knowledge on the HPG axis model used in our case study and its phenotypes in Section 2, we outline our reference downregulation protocol in Section 3. In Section 4 we state our problem and describe our algorithm to solve it. Section 5 describes and analyses the outcome of our multi-arm ISCT. Finally, Section 6 discusses related work and Section 7 draws conclusions. 2
In this section we give the necessary background knowledge. 2.1
We focus on the GynCycle Virtual Physiological Human (VPH) model, initially
presented in [
The model defines the time evolution of overall 33 biological quantities (mostly blood concentration of hormones) and the Pharmacokinetics/Pharmacodynamics (PKPD) of several pharmaceutical drugs used in assisted reproduction (in particular, GnRH analogues like Triptorelin, or FSH/LH–based drugs) by means of highly non-linear differential equations.
As such, GynCycle is a hybrid system whose inputs represent drug administrations. As it happens with VPH models of practical interest as GynCycle, the complexity of the differential equations makes their symbolic analysis (by, e.g., a model checker for hybrid systems) very challenging. Indeed, the model can be simulated by means of numerical integration in order to compute the time evolution of the hormones of interest upon administration of a sequence of doses of one or more drugs. Thus, GynCycle is used as an executable black-box model. 2.2
VPH models like GynCycle typically take into account inter-subject variability (i.e., the physiological differences among different individuals) by including suitable parameters in their equations. Different value assignments to model parameters yield different model time evolutions and/or different reactions to drug administrations, thus defining different Virtual Phenotypes (VPs). Intuitively, each VP represents a class of patients behaving similarly to each other.
Computing a complete set of VPs for the model at hand is the starting point to obtain a representative population of virtual patients (hence, ideally showing all possible phenotypes). In turn, this is the key enabler to perform In-Silico Clinical Trials (ISCT) to assess, e.g., safety and/or efficacy of a pharmacological treatment.
Unfortunately, computing the complete set of VPs defined by a complex
VPH model is all but easy from the combinatorial point of view, and may
require weeks of massive simulation on large High Performance Computing (HPC)
infrastructures. In [
The availability of such a representative population of VPs is the key enabler for our ISCT. 3
In this section we define the treatment we used as our case study: the downregulation phase of an assisted reproduction protocol currently administered at the Department of Reproductive Endocrinology of University Hospital Zurich.
At each menstrual cycle of a human female, among the several follicles initially present in the ovaries, only one, unless in exceptional cases, grows and reaches maturity. A complex competition among follicles takes place, driven by a hormone-based feedback loop, in order to inhibit the maturation of multiple follicles.
In a so-called “long protocol”, the treatment aims at neutralising (through drugs) such a feedback loop (downregulation phase), in order to achieve a controlled and an as-simultaneous-as-possible growth of 5–15 ovarian follicles (stimulation phase). When the first three follicles reach a size where a mature oocyte can be expected (about 18 mm in diameter), ovulation is induced, the oocytes are retrieved, those which are mature are tried to be fertilised in vitro, and implanted either immediately within the treatment cycle or later after cryopreservation back into the uterus.
Assisted reproduction treatments are complex and challenging, with low
average success rates (around 30%) even in the top clinics, and with many factors
that, to date, can be hardly kept under full control. Indeed, as hormonal
regulatory systems occur within a complex network of endocrinological, neurological
and psychological factors [
Our case study is one of the worldwide classically used downregulation protocols aiming at suppressing the usual hormonal oscillations of the menstrual cycle (as in Figure 1) and preparing the patient to the following stimulation. At the Deptartment of Reproductive Endocrinology of University Hospital Zurich, this protocol currently consists of a sequence of daily administrations of 0.1 mg of Triptorelin (a GnRH analogue). For physiological reasons, downregulation is started within a precise time window of the menstrual cycle, namely between 450 400 350 300 /pgLm220500 150 100 50 0 20 18 16 14 /L12 I10 U 8 6 4 2
TriptorelinNo (a) E2
TriptorelinNo (c) FSH day 21 and day 25. The treatment might have different duration in order to address different patient reactions. In particular, a downregulation treatment is considered successful (effective) for a patient in case the blood concentrations of a given set of hormones and other physiological quantities go below certain thresholds within 9 days from the first drug administration, and stay always below such thresholds for the following 21 days. As a consequence, a downregulation treatment might last up to 30 days. 4
In the following, we denote with R, R 0, R+ and N+ the sets of, respectively, all real, non-negative real, strictly positive real, and strictly positive natural numbers. 4.1
VPH models (as GynCycle) are hybrid systems. In the typical case (as ours) in
which model inputs are drug administrations, VPH models can be abstracted
into Discrete Event Systems (DESs) (see, e.g., [
Events for a DES defining a VPH model as GynCycle represent clinical actions, and take values from a finite alphabet A (in our case, defining all the possible doses for the treatment drug). We assume that a distinguished element nop exists in A, representing the “null” event (“no action”).
A discrete event sequence (defining an input function for our VPH model) is a function u(t) associating to each (continuous) time instant t 2 R 0 an event (i.e., a clinical action) in A, and such that u(t) differs from nop only in time points multiple of a given time quantum 2 R+. Intuitively, in our setting, a discrete event sequence u(t) defines a series of drug administration actions (including “no action”) occurring at time points multiple of .
Given an assignment to the VPH model parameters (i.e., a Virtual
Phenotype, VP) and a discrete event sequence u, the associated VPH model trajectory
is a continuous-time function x( ; t) representing the time evolution of the
biological quantities defined by the model when fed with u starting from its initial
state and with its parameters set to . Of course, as the VPH model is causal,
its trajectories up to any time point t 2 R 0 only depend on the restriction of
the input discrete sequence up to time point t [
Treatment invariants and goals define conditions that must be, respectively,
always and eventually satisfied by a successful treatment. Invariants and goals
for the treatment being sought can be modelled as continuous-time monitors
embedded within the VPH (DES) model, along the lines of [
In particular, given a model trajectory x( ; t) under a given input function u(t), our monitor output is Undet as long as x( ; t) satisfies the invariants (in other words, the monitor decision is “undetermined” as there is hope to extend the current treatment into a successful treatment) and goes to and stays at value Fail as soon as invariants are violated. When a Undet model trajectory satisfies the goal conditions, the monitor output turns to Success, informing the caller that the input function u(t) defines a successful treatment (i.e., a effective treatment which always satisfies invariants).
The use of continuous-time monitors embedded in the VPH model gives us a
flexible way to model both bounded safety and bounded liveness properties (see,
e.g., [
In our setting, the properties of interest are the conditions of successful downregulation treatments (Section 3). In particular, our invariant requires that the value of all the biological quantities under observation go and stay beyond their thresholds from the 9-th day after the first drug administration. Our goal condition instead requires that values for those quantities stay below their thresholds for 21 consecutive days.
invariants violated
goal reached
As a consequence, our monitor output is Undet in the initial state, and eventually turns either to Fail or to Success (see Figure 2). Model output turns to value Fail as soon as, from the 9-th day after the first drug administration, the value of any of the physiological quantities under observation goes beyond its threshold. Model output turns instead from value Undet to value Success as soon as all the thresholds are satisfied for 21 consecutive days. 4.3
As anticipated in Section 2 the complexity of the (highly non-linear) differential equations typically occurring in actual VPH models hinders the possibility for their symbolic analysis. Numerical simulation is often the only means to compute the time evolution of the model quantities of interests (mainly blood hormone concentrations in GynCycle) upon a sequence of clinical actions (administration of one or more drugs with their associated doses). Also, as VPH model equations are often stiff, simulation can be an expensive process from a computational point of view.
Our algorithm for optimal personalised treatment computation regards the
input VPH model (GynCycle in our case study) as an executable black-box
model. In order to drive the VPH model simulator, our algorithm assumes that
it accepts a set of basic commands (which are available or can be readily
implemented within most modern simulators [
In this section we describe the search-based algorithm that we used at the core of our In-Silico Clinical Trial (ISCT).
Given a VP of our GynCycle VPH model and a set of possible alternative doses for the drug at hand defining the set A of possible clinical actions (including action “no action”, i.e., dose = 0), our backtracking-based algorithm performs a depth-first search in the space of possible treatments (i.e., discrete sequences of clinical actions in A) lasting at most h = 25 + 9 + 21 = 55 days from day 1 of the patient menstrual cycle (i.e., the latest cycle day, 25, when the treatment can start, plus the latest day, 9, within which the safety conditions must be met, plus the number of consecutive days, 21, in which the safety conditions must be always satisfied in order to declare success).
The search algorithm is implemented in C and drives a GynCyle model
simulator (implemented within the Limex Solver [
The initial state of the GynCycle model represents a patient of VP on day 1 of her menstrual cycle. Furthermore, as described in Section 4.2, the model is equipped with a continuous-time monitor which checks at any time whether the current model trajectory satisfies the treatment safety and goal conditions or not.
The algorithm initialises the current treatment to the empty treatment and advances the simulator to cycle day 21 (the earliest day of the first clinical action) with no drug administrations. From this point, at each step it extends the current treatment with a new clinical action a 2 A (possibly nop, i.e., “no drug”).
Checking safety and goal conditions and backtracking After having appended to the current partial treatment a new action a 2 A, the action is injected into the simulator, the simulator is advanced by time , and the monitor output at the end of the last simulation step is retrieved, thus evaluating the safety and goals treatment conditions (as defined in Section 4.2) after one day from the last clinical action. The simulated (partial) treatment is said safe if the monitor output is Undet, and unsafe if the monitor output if Fail. Whenever the monitor output turns from Undet to Success, the current (safe) treatment represents a successful treatment. As soon as the monitor returns Fail, a backtracking is performed, with the previous simulator state being restored (together with the state of the monitor and the monitor output value in the restored state). If no further actions can be executed in the current node of the search tree, a backtrack is triggered which restores the simulator (plus monitor) state one level up in the stack.
Optimality constraint Our algorithm does not stop at the first found
successful treatment. Indeed, along the lines of [
At each node of the search tree, any action extending the current partial treatment with an administration of a drug dose which would make the overall administered drug amount reaching or exceeding Dmin is regarded as not applicable. Preference ordering among actions In a black-box setting as ours, no inference can be made on the effects of each candidate action without performing a simulator run to actually advance the model and then querying the model monitor output. Hence, approaches to compute, through inference, a dynamic preference order among the candidate actions to be tried during search (as those exploited in, e.g., classical planning, planning for white-box hybrid systems –see Section 6– or CSP, SAT or local search solvers) cannot be applied.
The only information available to the algorithm without running the simulator is value Dmin and the overall cost of the current partial treatment (overall amount of drug administered). Hence, given that our algorithm searches for a successful treatment of minimum cost and that any action has a non-negative contribution to the overall treatment cost, it not surprisingly tries the clinical actions on each node in the search tree in ascending order of their associated dose (i.e., cost), hence performing a greedy, optimistic choice.
The optimality constraint (whose threshold is updated each time a new optimal treatment is found), the presence of a bounded horizon h and the fact that our algorithm does not stop at the first successful treatment clearly guarantee that a global optimum is always returned (if any successful treatment exists). 5
In this section we outline how we set up and conducted our In-Silico Clinical Trial (ISCT). 5.1
The GynCycle model defines around 104 VPs (Section 2.1) whose time behaviour
(under no treatment) is shown in Figure 3 (light curves). In order to carry out
our multi-arm ISCT, we clustered such VPs by merging in the same cluster
phenotypes behaving similarly, where similarity has been defined (along the lines
of [
Not surprisingly, the reference downregulation treatment of Section 3 does not succeed on all possible phenotypes. Hence, we excluded from our ISCT those VPs for which the reference treatment fails. Our exclusion criterion directly stems from domain knowledge on downregulation treatments (namely: for these protocols, if the reference treatment, which envisions one full drug dose per day, fails for a patient, a lighter treatment will fail as well).
As a result of the application of clustering and of our exclusion criteria, we obtained a population of 98 VPs representative of the spectrum of all model VPs for which the reference downregulation treatment works (dark curves in Figure 3). Any two different VPs in our selection differ significantly for the time behaviour of at least one of the 33 biological quantities defined by the model (although they might appear similar on others, e.g., those in Figure 3).
Each of the 98 VPs in our selected population defines a distinct arm of our multi-arm ISCT. (c) FSH (d) LH We ran our 98-arm ISCT using a large High Performance Computing (HPC) infrastructure (the Marconi cluster) kindly provided by the Cineca consortium.
For each VP in the selected population, we ran our algorithm of Section 4 on an independent node of the cluster searching for an optimal (i.e., lightest) downregulation treatment for that VP, thus in an embarrassing parallel fashion. 5.3
Figure 4a shows the distribution of the computation times of each parallel process implementing a single arm of our multi-arm ISCT.
In particular, the average completion time of the arms of our ISCT is 19.09 h. Although standard deviation is very large (21.31 h), from Figure 4a it can be seen that the vast majority (around 90%) of ISCT arms terminate in less than 48 hours.
Figure 4b shows the distribution of the number of computation nodes among our ISCT arms. On average, 918 770 nodes have been expanded by each parallel search process in order to find an optimal treatment for the VP at hand. Also here, standard deviation is very large (1 144 123 nodes).
The high standard deviations of the distributions of computation time and number of expanded nodes among VPs can be explained by observing that our algorithm deterministically orders clinical actions to be selected in each search tree node in ascending order of the administered drug dose (i.e., on their contribution to the overall cost function). In Section 4.4 we have argued that any choice of the actions based on inference of the prospective effect of each action is a not viable option, because of the black-box nature of the Virtual Physiological Human (VPH) model. That said, the left skew of the distributions in Figures 4a and 4b suggests that our chosen action preference order does behave very well in the vast majority of cases.
Overall, Figures 4a and 4b show the revolutionary potential of ISCT: an ISCT over a population of 98 VPs has been conducted in a few days of parallel computation, whilst an equivalent in vivo clinical trial would have been not even conceivable. In fact, we would need to recruit in the clinical trial a distinct human patient for each phenotype and for each of the tested treatment variations. The overall number of human patients to be recruited would have been several millions, thus clearly making the approach infeasible.
Individualised treatments have the potential value to reduce costs and improve
outcomes of standard clinical treatments. In recent years data-driven techniques
have been investigated thanks to the availability of big data [
Computation Cumulative
Nodes < for lung cancer. Unfortunately in presence of scarce experimental data, the above approaches cannot be applied. For example, in our case study hormones blood concentrations are not measured every day, since those measurements are costly and invasive.
Model-based approaches, exploiting PharmacoKinetics (PK), as e.g., [
In our model-based setting, we have to face with complex Virtual
Physiological Human (VPH) models, e.g., HumMod [
As outlined in Section 4.1, such VPH models are hybrid systems that can be
abstracted into Discrete Event Systems (DESs) (see, e.g., [
However, as outlined in Section 4.3, the complexity of the differential equations of VPH models relevant for in silico clinical practice makes such models out of reach for symbolic approaches like those mentioned above, and appoints numerical integration as the only viable means to compute (black-box) the model evolutions under a given input function.
In particular, even considering that clinical actions have constant and equal
duration (as, e.g., in [
We presented a case study showing how model-based approaches coupled with intelligent search techniques can be used to support precision medicine by computing, for each given patient, a personalised treatment maximising effectiveness while minimising cost as well as likelihood and severity of adverse effects.
By exploiting a complex Virtual Physiological Human (VPH) model of the human female Hypothalamic–Pituitary–Gonadal (HPG) axis, its set of Virtual Phenotypes (VPs), and a backtracking-search algorithm that intelligently drives a simulator for the VPH model, we set up and conducted a multi-arm In-Silico Clinical Trial (ISCT) which computes, for each VP satisfying our inclusion criteria, a personalised variation of a reference pharmaceutical treatment which, preserving effectiveness, minimises the overall amount of drug used.
Our ISCT can be carried out with a few days of parallel computation on a regular High Performance Computing (HPC) infrastructure, while an equivalent (classical) in vivo trial would require the recruiting of several millions of human patients in order to test all treatment variations assessed in silico, thus resulting just infeasible.
Although our trial shows extremely favourable results (with computed personalised treatments using only a tiny fraction of the drug quantity employed by the reference treatment), there are sill open issues in the area of ISCT. In particular, mapping each human patient to her corresponding VP calls for novel and highly inter-disciplinary approaches at the intersection of Artificial Intelligence, Medicine, Mathematics and Biology.
Acknowledgements. This work was partially supported by the Italian Ministry of University and Research (MIUR) under grant “Dipartimenti di eccellenza 2018–2022” of the Department of Computer Science of Sapienza University of Rome and by the EC FP7 project PAEON (Model Driven Computation of Treatments for Infertility Related Endocrinological Diseases, 600773).