Virtual clinical trials identify effective combination therapies in ovarian cancer

A major issue in oncology is the high failure rate of translating preclinical results in successful clinical trials. Using a virtual clinical trial simulations approach, we present a mathematical framework to estimate the added value of combinatorial treatments in ovarian cancer. This approach was applied to identify effective targeted therapies that can be combined with the platinum-taxane regimen and overcome platinum resistance in high-grade serous ovarian cancer. We modeled and evaluated the effectiveness of three drugs that target the main platinum resistance mechanisms, which have shown promising efficacy in vitro, in vivo, and early clinical trials. Our results show that drugs resensitizing chemoresistant cells are superior to those aimed at triggering apoptosis or increasing the bioavailability of platinum. Our results further show that the benefit of using biomarker stratification in clinical trials is dependent on the efficacy of the drug and tumor composition. The mathematical framework presented herein is suitable for systematically testing various drug combinations and clinical trial designs in solid cancers.

the same phenotype emerges, ii) unfaithful division where during division new cell type appear with rate u per cell division, and iii) death where cell dies and is removed from the system.
We have included three main platinum resistance mechanisms (pre-target, on-target and post-target). As we consider all possible cell types, the model with three platinum resistance mechanisms contains eight subclones: 000, 100, 010, 001, 101, 110, 011, 111, where 0(1) indicates the presence (absence) of a platinum resistance mechanism. The first, the second and the third digit indicates the presence (absence) of pre-target, on-target, and post-target platinum resistance, respectively. In the absence of treatment intervention, all cells grow exponentially with the same birth rate b. It means that we consider the evolution of platinum resistance as neutral. In addition, all cells die with death rate d leading to net growth rate λ = b − d. In the presence of chemotherapy treatment, however, sensitive cells (000) die with additional rate d chemotherapy proportionally to birth rate leading to new death rate d = d+b·d chemotherapy . Fully resistant cells are not affected by chemotherapy, i.e., d chemotherapy = 0. Cells with one and two platinum resistance mechanisms have d chemotherapy decreased by α 1 and α 2 , respectively. Thus, in the presence of treatment intervention, partially and fully resistant cells have selection advantage over sensitive ones.
In addition to the faithful division, where division leads to the occurrence of the second cell with the same phenotype, a cell can divide where new subclone with an additional resistance mechanism appears. In the model, we do not take into account cross-resistance, and as a result, during one division, only one additional platinum resistance mechanism can be acquired. In addition, we also assume that the acquisition of platinum resistance mechanism is irreversible, leading to transitions between all eight cell types, as shown in Figure 1. As innate resistance is known as the dominant one, and acquired resistance has a minor effect, transition rate u is the same with or without chemotherapy and targeted therapy treatment intervention.
In the model, we also included three types of treatment interventions: debulking surgery, platinumbased chemotherapy, and targeted treatment. Surgery is modeled as the removal of a fraction of β cancer cells, and all cell types can be eliminated by the surgery with an equal probability. Chemotherapy is modeled by increasing the death rate of sensitive and partially resistant cells proportionally to the division rate. In the model, chemotherapy is administered as a finite time-continuous treatment with maximum tolerated dose (MTD), which is the current standard-of-care in high-grade serous ovarian cancer (HGSOC). Targeted treatment is modeled in combination with platinum-based chemotherapy. The details of inclusion of targeted treatment in the model are presented in the next section.

Modeling targeted treatment
Here, we focus on modulators of platinum resistance, which sensitize platinum-resistant cells to platinum-based chemotherapy. As in the model we consider three platinum resistance mechanisms; we included three classes of modulators where each one target different platinum resistance mechanism. The first class modulates platinum transport by increasing platinum concentration in the cell and thus targets pre-target platinum resistance. The second class regulates DNA damage response (DDR) by inhibition of the efficient DNA repair machinery in platinum-resistant cells. Thus, this class of modulators targets on-target platinum resistance. The last class of platinum sensitivity modulators targets post-target platinum resistance which is responsible for enhancing apoptosis rate.
One the most promising method to overcome platinum resistance is a combination of platinum-based chemotherapy with drugs targeting given platinum resistance mechanism [6]. In addition, the majority of the HGSOC patients have subclones that confer resistance to platinum at the time of diagnosis. These lead to the hypothesis that the best strategy to overcome platinum resistance would be a combination of platinum-based chemotherapy with a targeted therapy that hits the dominant resistance mechanisms as the first-line treatment.

Trientine
Cooper transporter receptor 1 (CTR1) is a major copper influx transporter that mediates transport of platinum to the cells. Clinical studies demonstrated that expression of CTR1 correlates with intratumor platinum concentration leading to better clinical outcome [10]. Copper-lowering agents, such as trientine, increase the expression of human CTR1 (hCTR1), which leads to resensitization of tumor cells to platinum [5].
In the model, trientine is included by increasing simultaneously two parameters α 1 and α 2 in partially resistant cells which are also pre-target resistant to platinum, i.e., 100, 110, and 101. The parameter α 1 and α 2 defines how much chemotherapy killing effect is reduced in partially resistant cells in comparison to sensitive cells. As trientine increases platinum accumulation resulting in enhanced killing effect of platinum-based chemotherapy in pre-target resistant cells.

Wee1 inhibitor
An essential mechanism of action of platinum is introducing interstrand crosslinks to DNA, which need to be recognized and repaired by the DNA repair machinery. One important platinum resistance mechanism is the upregulation of the DNA damage response to cope with the increased replication stress induced by platinum. Preclinical and early clinical studies indicate that cell cycle checkpoint inhibitors, such as Wee1 inhibitors, sensitize platinum-resistant cells to platinum-based chemotherapy. The predominant mechanism-of-action of Wee1 inhibitors is to cause the failure of G2-M checkpoint due to inappropriate CDK1/cyclin B complex activation resulting in mitotic catastrophe.
In the model, we included Wee1 inhibitors by adding a reverse transition from fully-resistant to partially resistant cells and from partially resistant cells to the sensitive cell. There is a probability u reverse per cell division that the cell will reverse back from fully resistant cell to partially resistant cell and from partially resistant cell to sensitive cell (e.g., from x011 to x010). Wee1 inhibitor has no single-agent activity in primary platinum-sensitive tumors [12]. Thus, we assume that Wee1 inhibitor does not kill platinum-resistant cells, but only sensitize cancer cells to platinum chemotherapy by reversion transition.

Birinapant
Cancer cells can resist programmed cell death, for example, by overexpression of proteins blocking proapoptotic pathways. Inhibition of apoptosis is a mechanism of platinum resistance, and thus activation of cancer cell death through apoptosis is an important process to resensitize cells to platinum-based chemotherapy. Upregulation of the gene inhibitor of apoptosis (IAP) is one mechanism leading ti apoptosis evasion. Thus, inhibition of IAP by an IAP inhibitor, such as birinapant, could sensitize cancer cells to platinum by activation of cancer cell death.
We included birinapant in the model by adding additional death with the rate d a called apoptotic cell death. Thus, cells with active post-target platinum resistance have death rate equal to d = d + b · d chemotherapy + d a during treatment phase. As a result, cancer cells with active post-target platinum resistance mechanism have a higher death rate than platinum-sensitive ones.
3. Platinum-sensitive cells can acquire one of three platinum resistance mechanism with equal probability of u per cell division.
4. During each cell division, sensitive or partially-resistant cells can acquire only one platinum resistance mechanism, leading to combinatorial network describing transitions between subclones as shown in Figure 1.
5. In the model, only innate platinum resistance is considered.
6. The platinum-resistant cells grow at the same rate as wild-type cells, i.e., they are neutral in the absence of treatment.
7. Cancer is diagnosed when the colony reaches a certain size (M ).
8. Rate of accumulation of platinum resistance in the absence of treatment is the same as during the treatment, i.e., u is the same in all three phases of simulations.
9. Two types of treatment modalities (debulking surgery and platinum-based chemotherapy) are included in standard-of-care simulations.
10. The effect of chemotherapy is modeled by the drug-induced death rate, d chemotherapy , which leads to negative net growth rate: 11. Chemotherapy toxicity is not considered in the current model implementation.
12. Chemotherapy administration is included in the model as a continuous-time maximum tolerated dose (MTD).
13. Surgery is modeled by the removal of a fraction of β cancer cells. All cell types can be removed with an equal probability.
14. Tumor recurrence is observed when the number of cancer cells reaches M relapse after end of primary treatment.

Computer simulation
We performed Monte Carlo computer simulation using an approximation of an exact stochastic simulation algorithm developed by Gillespie [7]. We used an approximation algorithm due to its computational efficiency. First, slow events (unfaithful division) are separated from fast events (faithful division, cell death). The propensities of slow reactions are: a i (x), for i = 1, 2 . . . n, where n is the number of slow events. The propensities of fast reactions are: where m is the number of fast events. Next, time of the next slow reaction is calculated: where R 1 is a random number from a unit-interval uniform distribution. Next, the index of next slow event is given as the smallest k satisfying: where R 2 is a random number from a unit-interval uniform distribution. The system is then updated and fast events that occur until time of the first slow event τ . Fast events are updated using tauleaping, i.e., for every fast event the number of times each event occurs during the time interval [t, t + τ ) is: and the system is updated: where v j is the state change vector.

Text S2 Standard-of-care simulations
To investigate resistance to platinum-based chemotherapy and suggest novel methods to resensitize HGSOC patients to standard-of-care treatment, we performed computer simulations of the model. The simulation consists of three phases: 1) pre-treatment phase in which the tumor grows until diagnosis, 2) treatment phase that includes platinum-based chemotherapy and debulking surgery, and 3) post-treatment phase where the tumor grows until the first relapse.
The first phase starts with a single platinum-sensitive cell and proceeds with exponential tumor growth in the absence of treatment interventions. In addition to the growth of the platinum-sensitive cells, accumulation of platinum resistance mechanisms is also simulated by the appearance and the growth of partially and fully resistant clones during unfaithful cell division. The phase ends when the total number of tumor cells reaches M cells. The tumor is then diagnosed, and the treatment phase starts.
Treatment phase starts after the patient is diagnosed with HGSOC, i.e., we assume no delay in treatment initiation. Standard-of-care includes two types of interventions: platinum-based chemotherapy and debulking surgery. Chemotherapy reduces the net growth rate of sensitive and partially-resistant cells. As in the model, we do not consider pharmacokinetics and pharmacodynamics of platinumbased chemotherapy; we simulate platinum-based chemotherapy as a continuous injection. Debulking surgery is modeled by removing a fraction of tumor cells at a one-time point. Every cell can be removed during debulking surgery with the same probability. The treatment phase is simulated in a framework that follows the first-line standard-of-care therapy guidelines for HGSOC. Thus, in our simulations, the HGSOC patient is treated with three cycles of neoadjuvant chemotherapy (NACT). Next, a debulking surgery is performed, aiming for optimal cytoreduction, followed by another three cycles of adjuvant chemotherapy.
The last phase of simulation starts after the treatment phase and ends at the time of the first recurrence. In this phase of simulation, the tumor grows without treatment intervention until the total number of cells reaches M relapse cells.

Text S3 Simulation of virtual HGSOC patients cohort
We performed in silico treatment response analysis using the simulation approach described in Supplementary Text 2. Each virtual HGSOC patient is described with two parameters: tumor burden at diagnosis M and chemotherapy effect d chemotherapy .
For each virtual HGSOC patient, we sampled tumor burden at diagnosis (M ) and chemotherapy effect (d chemotherapy ) from a log-normal distribution. The log-normal distribution was selected because it describes the best distributions for M and d chemotherapy among six tested statistical distributions (details in Supplementary Text S4). The goal of sampling is to reproduce in silico real-life patient's heterogeneity of response to standard-of-care treatment.
Next, we simulated each virtual HGSOC starting from a single platinum-sensitive cell until the first relapse. The goal of the simulations was to compute the platinum-free interval from virtual HGSOC patients as well as the proportion of all cancer cells subclones at the time of diagnosis and after first-line treatment.
In our analysis, we perform simulation of 1,000 HGSOC patients to estimate response to standard-ofcare treatment. The computed PFI estimates from virtual HGSOC patients were applied to perform Kaplan-Meier analysis. As a result, we created platinum-free interval survival plots. Kaplan-Meier analysis was carried out in MATLAB using ecdf function with no censoring as all virtual HGSOC patients relapse after first-line treatment.

Text S4 Parameter calibration
Tumor growth The average cell cycle duration in an ovarian cancer cell is 36 hours [17], and is defined as the time between cell divisions in the absence of cell death. This implies that the division rate is: b = To achieve the observed λ, we set the cell death rate to d = λ − b. In our simulations, the parameters b and d are the same for all cell types and constant in absence of treatment.

Size of tumors at clinical diagnosis and chemotherapy effect
We modeled the probability of diagnosis previously as a function of tumor size as follows [11]. First, we estimated tumor burden using metabolic tumor volume (MTV) values extracted from 18 F-FDG-PET/CT images. Next, the MTV values were converted to a number of cells assuming 10 9 cells in 1cm 3 tumor bulk according to [3]. The number of cells was fitted to normal, exponential, Weibull, log-normal, logistic, and log-logistic probability distributions. The Bayesian Information Criterion (BIC) was calculated to measure the goodness of fit (GoF). The best agreement with the data was obtained for log-normal probability distribution with parameters µ = 26.59 and σ = 0.47. Converting µ of the log-normal distribution to µ of the normal distribution, the total number of cancer cell at diagnosis equals M = 3.9548 · 10 11 cells.
Patients in the calibration cohort received a median of three cycles of NACT where each cycle duration is 21 days. This means τ N ACT = 63 [days] of NACT. Without loss of generality, we can assume that τ N ACT is also the time interval between the acquisition of PET/CT images before and after NACT. Parameter d chemotherapy was thus estimated for each patient separately as follows. We set all parameters except d chemotherapy and M to default values listed in Table 2. Next, for each patient, we extracted a number of cells before (M ) and after (M N ACT ) chemotherapy. In addition, we know that time between two PET/CT scans is τ N ACT = 3 [weeks]. Thus, we estimated d chemotherapy numerically by with the bisection method to find a minimum of the following objective function:

Debulking surgery effect
The optimal debulking surgery is defined as a residual disease of 10 [mm] or fewer [18]. Tumor with 10 mm diameter corresponds to tumor volume of 10 [cm 3 ]. HGSOC is diagnosed when the tumor burden is of order 10 11 cells. Thus, optimal surgery removes at the minimum fraction of β 10 9 10 11 and corresponds to two cell-log kill fraction of cells removed by surgery.

Rate of resistance accumulation
We utilized clinical data from Koz lowska et al. 2018 to estimate the value of u [11]. Using a cohort of 33 HGSOC patients (training cohort), we performed Kaplan-Meier analysis to create platinum-free interval survival plot from the patient data.
Next, for a wide range of transition rate (u), the weight of chemotherapy effect on cells with one platinum resistance mechanism accumulated (α 1 ), and weight of chemotherapy effect on cells with two platinum resistance mechanisms (α 2 ), we performed a grid search. We took 100 values of u between 10 −8 − 10 −4 and 101 values of α 1 and α 2 between 0 − 1 leading to 1,020,100 parameter combinations and performed simulations of 1,000 HGSOC patients for each parameter combination.
For each simulated cohort, Kaplan-Meier analysis was performed, and platinum-free interval survival plot was created. Next, the obtained survival plot from simulations was compared with one from the training cohort. Mean squared error (MSE) was calculated to measure the deviation between the two survival estimates. To take into account stochasticity of model simulations, for each value of u, 100 virtual HGSOC patient cohorts were simulated, and the mean value of MSE was calculated as a GoF estimate. The smallest MSE was obtained for: u = 10 −5 , α 1 = 0.02 and α 2 = 0.01, and was set as a nominal value in our simulations.

Targeted treatment parameters
The model includes three targeted treatment drugs: trientine, Wee1 inhibitor, and birinapant. The drugs are incorporated in the model of platinum drug resistance in HGSOC by changing corresponding model parameter value (see Supplementary Text 1), i.e., α 1 and α 2 (trientine), u ddr (Wee1 inhibitor) and d a (birinapant).
To estimate values of the treatment parameters, first, we estimated total relative drug efficacy using efficacy values from: in vitro, in vivo and early clinical trials, which were extracted from the literature. The Table S1 lists relative drug efficacy at three levels of evidence: in vitro (I), in vivo (V ) and early clinical (C), which were applied to calculate total relative drug efficacy (E) using the following formula:  The values of E were applied to fit treatment parameters for selected drugs as follows. First, for each drug, ten parameter values (p i ) were taken uniformly between [p min , p max ]. Next, for each value of p, the total relative efficacy was estimated via simulations of virtual HGSOC patients. The value of p is selected in such a way that E from the model simulations agree with the E extracted from literature (see Table S1).
Estimation of E from the model simulation requires further comment. We assume that the total relative drug efficacy could be measured by the tumor burden (the total number of cancer cells) after the first-line treatment, and estimated with the following formula: where X i is a sum of all cancer cells from the model simulations after first-line treatment (i.e., SOC alone), and X j is the sum of all cancer cells from the model simulations after first-line treatment (i.e., SOC and targeted therapy). Thus, the expected tumor burden after the treatment with SOC and targeted therapy equals:X i =Ê ·X j , and the goal of targeted treatment parameter calibration is to find the parameter value whereÊ from model simulation equals the one estimated and listed in Table S1 (E). For a given drug, the following procedure was performed to estimate tumor burden after primary treatment as a function of p. For each of ten values of p, a cohort of 1,000 HGSOC virtual patients was simulated. As the model stochastic and thus, it is computationally intensive to simulate a large cohort of HGSOC virtual patients, we performed bootstrapping to estimate the value of p for each drug together with confidence intervals. Thus, from a cohort of 1,000 HGSOC virtual patients, 1,000 patients were sampled with replacement, and average tumor burden after primary treatment (X j ) was extracted. The procedure of sampling and extractingX j is performed for all ten values of p, and plot ofX j as a function of p was created. Next, linear regression for performed to estimate coefficients of a linear function. As the last step, the value of p for a given drug was estimated using the formula: where α 1 and α 2 are coefficients of linear function.
The procedure from creation of bootstrap sample to estimation of p using above equation was performed 1,000 times. The value of p for a given drug is then calculated as a mean value over 1,000 estimates and 95% confidence interval around estimated p drug was calculated. The pseudocode describing the method to estimate targeted treatment parameters is presented in Algorithm 1 and the Table S2 presents estimated targeted treatment parameters together with confidence intervals.
Algorithm 1 Estimation of targeted treatment parameters 1: procedure Fit targeted treatment parameter 2: Calculate: X i =Ê ·X j

3:
Set: p min , p max and N patients

4:
Simulate HGSOC virtual patient cohort for each p j :

5:
for p j = linspace(p min , p max , 20) do 6: SimulateVirtualPatientsCohort(p j ) Plot X i as a function of p 16: Perform linear regression to estimate α 1 and α 2 of: X i = α 1 · P j + α 2

17:
Calculate:p i = Xi−α1  Table S2: Parameter values of targeted treatment parameters. Figure S3: Result from virtual clinical trial simulations of two and three drugs combinations. The waffle plots illustrate tumor composition after primary treatment. Blue, yellow, orange and red color represent sensitive, cells with one platinum resistance mechanism active, cells with two platinum resistance mechanisms active and fully resistant cells.