In this short paper we present preliminary results on computing, in silico, personalised therapies for Colorectal Cancer (CRC), one of the deadliest form of tumour for adult humans. We exploit a recent SBML (Systems Biology Markup Language) model of the tumour growth, which also models the Pharma-coKinetics/Pharmacodynamics (PK/PD) of two immunotherapic drugs that may be used in combination treatments.
The recent availability of quantitative models of the human patho-physiology (Virtual Physiological Human, VPH) has inspired and facilitated new approaches to the design of pharmacological treatments and the safety and efficacy assessment of therapies and biomedical devices. These kinds of approaches, collectively referred to as In Silico Clinical Trials (ISCTs), hold the promise to enable precision medicine, in which Artificial Intelligence (AI) methods support the design of personalised therapies and the assessment of their efficacy by means of simulation (i.e., in silico). ISCTs require computational VPH models that take into consideration the physiological differences between different human individuals and that capture the kinetics and dynamics of pharmacological drugs. Quantitative VPH models are typically defined as hybrid systems, where the dynamics is described by systems of Ordinary Differential Equations (ODEs) and the inter-patient variability is encoded by parameters. These models are often designed and distributed in open-standard modelling languages, often the Systems Biology Markup Language (SBML).
The automatic synthesis of therapies that may involve several drugs is a computationally complex task. In fact, it is necessary to determine which drugs to use and their associated dosing regimen, i.e., the drug amounts and the frequency of administration. Also, as it is the case for most pathologies and drugs, not all possible therapies are actually feasible; most notably, there exist known constraints that limit the quantity of drugs that may be administered in a certain period, in order to deal with drug toxicity. In the general case, the number of possible therapies is infinite. We study an approach to the problem that is based on the parametrisation of therapies. In order for the approach to be both feasible and effective, the number of parameters must be small enough to enable an efficient search over the possible therapies but still large enough to allow the modelling of all (or most of) the therapies of interest. Our goal, then, is the synthesis of therapies (assignment to the parameters) that meet a given set of constraints and optimise a user-defined objective function, which usually combines a set of performance metrics that include efficacy (i.e., how well the therapies cures the pathological condition) and the total quantity of administered drugs (generally to be minimised, so to reduce costs and health risks). We propose a method that uses VPH models described in the SBML language and that defines the whole optimisation problem inside the model itself. In particular, the parameters of the treatment will be modelled as additional parameters of the VPH model and the objective function and the therapy constraints will be modelled as observable model outputs. The fact that the whole problem is defined in pure SBML, a standard language supported by a large suite of software and with a large community, has the goal of showing that computational methods have the potential to be used in an easy and immediate way in clinical contexts.
This section describes the steps we followed to set up an ISCT to evaluate the effectiveness of our approach in the synthesis of personalised therapies for the CRC model proposed in [1]. Generation of a virtual population. The starting point to perform an ISCT is the availability of a complete and representative population of virtual patients (VPs), i.e., assignments to the VPH model parameters. The CRC model presents 23 real-valued parameters that encode the inter-patient variability. We used the approach described in [2,3,4] to generate a population of VPs that are of interest (i.e., do actually develop a tumour) and show evolution of the model observables that are physiologically admissible (i.e., do not violate the laws of biology). Therapy modelling. We modelled 60-week-long therapies as assignments to a set of 32 parameters. For each of the two drugs, atezolizumab and cibisatamab, each one of 15 parameters governs the amount of drug that can be administered during a 28-day period and one additional parameter defines the number of weeks between two consecutive administrations of the drug. Therapy constraints and objective function. We modelled two constraints, each of which limits the total amount of each drug cumulatively administered in each 28-day period according to known toxicity levels and past clinical trials (https://clinicaltrials.gov/ct2/show/NCT03866239). The objective function combines the total amount of administered drug with a measure of inefficacy of the treatment. This last measure is computed based on the volume of the tumour at the end of the treatment and its initial volume. An high value for the inefficacy measure means that the therapy is not able to keep the tumour growth under control.
We chose to carry out the experimental evaluation of our approach using COPASI [5], one of the most well-known software tools supporting both plain simulation and simulation-based parameter optimisation of SBML models via iterative improving algorithms. The starting point of our ISCT was the random sampling of 35 VPs from the previously computed complete population. We compared the optimisation performance of 3 algorithms implemented in COPASI: the standard Genetic Algorithm (GA), Particle Swarm Optimisation (PSO), and the Levenberg-Marquardt (LM) algorithm, a gradient-descent based method that combines Steepest Descent and the Newton Method. The hyper-parameters of the algorithms were chosen as to perform the optimisation in a reasonable time for every patient (within 2 hours on an Intel Xeon E5430 @ 2.66GHz (8 cores) 32GB RAM machine). In order to compare the quality of the personalised therapies synthesised by each algorithm with a common baseline, we considered the dosing regimen of [1] as a a reference. Given a therapy, we define the Drug Amount Percentage as 𝑡amount 𝑟amount × 100%, where 𝑡 amount and 𝑟 amount are the total amounts of drugs administered by the given therapy and the reference one, respectively. For each VP, the Inefficacy Percentage for a given therapy is defined as 𝑡 ineff 𝑟 ineff × 100%, where 𝑡 ineff and 𝑟 ineff are the inefficacy values for the given and the reference therapy, respectively. The objective function is a linear combination of these two metrics, where the coefficients are chosen so to balance the search for effective treatments and the minimisation of administered drug amounts. In our experiments, only 11 VPs out of 35 showed a response to the drugs. This is in agreement with the results from [1].
For such 11 patients, LM was not able to find a therapy that improves the reference one, while GA and PSO show, on average, good results. The average reduction of the amount of administered drug is as high as 96.9% for GA and 98.62% for PSO, while the average reduction of the inefficacy is around 35% for both algorithms. Nonetheless, such personalised therapies still manage to reduce the tumour growth significantly.
Many attempts have been proposed in the literature to solve large optimisation problems defined via logic-, automata-or constraint-based formalisms (e.g., [6,7,8,9,10,11,12] among others). However, such approaches cannot be applied when the problem model (a complex ODE-based VPH model, as in our case) cannot be accurately defined within such formalisms and is available only as a simulator. Indeed, although such VPH models are hybrid systems whose inputs are discrete event sequences [13,14], to find an optimal treatment means to find an optimal plan in hybrid domains. Although symbolic approaches exist to model and solve planning problems in hybrid domains [15,16], the typical complexity of the ODEs of clinically-relevant VPH models makes such models out of reach for them, and appoints numerical integration as the only viable means to compute (black-box) the model evolutions under a given input function.
The synthesis of personalised therapies exploiting black-box VPH models is addressed in, e.g., [17,18]. In [19,20] the authors propose a intelligent backtracking simulation-based algorithm guided by multiple heuristics to seek, on a patient digital twin, an optimal robust personalised treatment for assisted reproduction, an area with many factors hard to control [21,22,23] and for which treatment personalisation is crucial for success.
One of the main problems in system biology is the estimation of unknown model parameters that fit a series of experimental data. Various optimisation algorithms are studied in [24,25,26] and applied in real-world case studies in, e.g., [27,28,29]. Many of the available tools rely on SBML simulators (see www.sbml.org). Among such simulators is SBML2Modelica [30], which focuses on the interoperability between system biology and (hybrid) CPS domain by translating SBML models to Modelica and the FMI/FMU open standard. This enables the seamless exploitation of tools and methodologies already established for CPS optimisation and verification, in particular backtracking-based search and optimisation via the efficient storing and retrieval of intermediate simulator states [31], verification of closed-loop systems also in presence of uncontrollable events (e.g., [32,33,34,35,36,37,38]).
The good results of the GA and PSO algorithms show the potential of the approach in the synthesis of personalised therapies. We interpret the failure of the LM algorithm as a clue of the fact that purely gradient-based optimisation is not suited for this problem, due to the strong constraints enforced on the therapies. Conversely, population-based algorithms show good results thanks to their ability to widen their focus throughout the search space.