Abstract
The usefulness of mechanistic models to disentangle complex multiscale cancer processes, such as treatment response, has been widely acknowledged. However, a major barrier for multiscale models to predict treatment outcomes in individual patients lies in their initialization and parametrization, which needs to reflect individual cancer characteristics accurately. In this study, we use multitype measurements acquired routinely on a single breast tumor, including histopathology, MRI, and molecular profiling, to personalize parts of a complex multiscale model of breast cancer treated with chemotherapeutic and antiangiogenic agents. The model accounts for drug pharmacokinetics and pharmacodynamics. We developed an open-source computer program that simulates cross-sections of tumors under 12-week therapy regimens and used it to individually reproduce and elucidate treatment outcomes of 4 patients. Two of the tumors did not respond to therapy, and model simulations were used to suggest alternative regimens with improved outcomes dependent on the tumor's individual characteristics. It was determined that more frequent and lower doses of chemotherapy reduce tumor burden in a low proliferative tumor while lower doses of antiangiogenic agents improve drug penetration in a poorly perfused tumor. Furthermore, using this model, we were able to correctly predict the outcome in another patient after 12 weeks of treatment. In summary, our model bridges multitype clinical data to shed light on individual treatment outcomes.
Mathematical modeling is used to validate possible mechanisms of tumor growth, resistance, and treatment outcome.
Introduction
Current personalized cancer treatment is based on a few biomarkers that allow assigning each patient to a subtype of the disease for which treatment has been established (1). In breast cancer, multigene tests like Mammaprint or PAM50 give prognostic information to guide clinical decisions (2, 3). Such stratified patient treatments represent a first important step away from one-size-fits-all treatment. However, the accuracy of disease classification comes short in the granularity of the personalization: it assigns patients to one of the few classes within which heterogeneity in response to therapy is still large (4). In each of these classes, randomized clinical trials (RCT) can be run to compare a few treatment regimens and to identify the best on-average one. As there is a combinatorial explosive quantity of combinations of cancer drugs, doses, and regimens, only very few can be explored by RCTs. The concept of “one disease in each patient” challenges diagnosis and therapy.
Mechanistic mathematical modeling and simulation have emerged as a powerful approach to investigate the influence of biological factors on tumor progression and therapy response (5–7). Current models are able to account for complex interactions at the cellular and molecular level and are capable of bridging multiple spatial and temporal scales in ways that would be impossible using experimentation (8, 9). Successful multiscale models can describe, with acceptable approximation, the dynamics of tumors under the effect of a specific therapy. The current computing capacity allows exploration of a large number of treatment regimens by running multiple model simulations in parallel. Yet, considerable challenges hinder the use of computational modeling to guide patient's treatment in clinical practice (7, 10). The choice of the level of biological detail represented in the mathematical model can always be questioned and a judicious balance between biological and model complexity has to be found. An essential question is whether mathematical models can be personalized to predict the effect of therapies in each specific patient (7, 11). This requires individualized initialization and parameterization of such a model that are typically difficult to perform.
For this article, our team of oncologists, pathologists, molecular biologists, medical imaging physicists, statisticians, and mathematicians examined whether a computational modeling approach can effectively simulate individual outcomes of a group of patients with invasive breast cancer receiving different treatment regimens (described later). Is it possible to design a mathematical model that integrates individual patient's data collected in routine clinical practice to reproduce patients’ outcome? Can such a model be flexible enough to simulate a wide range of possible responses by capturing fundamental biological processes at an acceptable level of approximation? Are the available data sufficient to personalize such a mathematical model? This study answers positively the first two questions. We developed a model that allows in silico simulation of the treatment outcome. Regarding the third question, we show the need to measure more precisely certain individual drug and tumor characteristics, to estimate certain parameters that modulate the dynamics.
Barbolosi and colleagues (7) underlined the importance of capturing pharmacokinetics (the fate of drugs in biological tissues) and pharmacodynamics (their mechanisms of action and effect) in a mathematical model. Both are present in our model. Furthermore, ref. 7 distinguishes between phenomenological (descriptive) models and mechanistic (explicative) ones. Phenomenological models are parametric models of the time evolution of the tumor that do not rest on explicit biological mechanisms. Mechanistic models instead account for specific biological details and dynamics. We present here a mechanistic model that captures nonlinear, multiscale dynamics in space and time. We consider discrete individual cells and blood vessels together with continuous quantities like oxygen and drug concentrations; such models are termed “hybrid.” Gallo and Birtwistle (12) used the term “enhanced pharmacodynamics” for models that merge multiscale dynamics with pharmacokinetics, as is in our case. An important distinction between the models is whether they are informed by individual patient data; those that are not can therefore not be used for personalized treatment optimization. While there are many studies in the literature that propose models in some of these dimensions (e.g., refs 13–25), our study is possibly the first one toward personalized computer simulation of breast cancer treatment incorporating relevant biologically specific mechanisms and multitype individual patient data in a mechanistic and multiscale manner.
We inform our model with data from five breast tumors collected in a recent neoadjuvant clinical phase II trial (26). Patients included in this study were randomized in two arms to receive chemotherapy with or without bevacizumab. Histologic, MRI, and molecular data were collected before, during, and at the end of neoadjuvant treatment. We developed precise pipelines for clinical data preprocessing, model initialization, and personalization. Besides individual histologic and MRI data, our model also makes use of some genomic patient data. By means of extensive numerical simulations, we show, as a proof of concept, that patient-specific and multiscale modeling allows us to reproduce treatment outcomes of the 4 patients. We also correctly predicted the outcome after 12 weeks treatments of one further patient. In addition, we investigate if and how alternative treatment protocols would have produced different outcomes. This is a first step toward virtual treatment comparison. Finally, we suggest which additional data and experimentation with tumor material is needed to improve the accuracy of the simulation toward precise outcome prediction. Our study shows that simulation-based personal treatment optimization is feasible and powerful and should be developed further as a promising avenue of personalized treatment.
Patients and Methods
Patients and treatment
We selected 5 patients with HER2-negative mammary carcinomas from the NeoAva cohort (26), a randomized, phase II clinical trial that evaluated the effect of bevacizumab in combination with neoadjuvant treatment regimens for 24 weeks. Written informed consent was obtained from all patients prior to inclusion. The study was approved by the Institutional Protocol Review Board, the regional ethics committee, the Norwegian Medicines Agency, and carried out in accordance with the Declaration of Helsinki, International Conference on Harmony/Good Clinical practice. The study is registered in the http://www.ClinicalTrials.gov/ database with the identifier NCT00773695. See section S2.1 in Supplementary Material for more details about the trial inclusion criteria.
For simplicity, we analyzed only the first 12 weeks, where patients were treated with chemotherapy (FEC100, 4 courses) and randomized 1:1 to receive bevacizumab (15 mg/kg every third week, given concurrently with chemotherapy) or not. Five patients were selected to belong to both arms of the trial and to have either a complete or no response by clinical examination and caliper measurements at 12 weeks of treatment. In addition, MRI was used to facilitate precise validation of simulated outcomes in small tumor portions. An overview of baseline characteristics, treatment, and response of the 5 patients is shown in Table 1. Further details on clinical data can be found in Supplementary Table S1 and Supplementary Fig. S1–S10.
Patient overview
Patient ID . | 1 . | 2 . | 3 . | 4 . | 5 . | |
---|---|---|---|---|---|---|
Age | 48 | 31 | 40 | 41 | 53 | |
Histologic tumor type | IDCa | IDC | IDC | IDC | IDC | |
Histologic tumor grade | 2 | 3 | 3 | 2–3 | 3 | |
Hormone receptor status | ER | Positive (100%) | Positive (1%) | Negative | Negative | Negative |
PR | Positive | Negative | Negative | Positive | Negative | |
Lymph node (LN) status | Week 0 | Positive | Negative | Positive | Negativec | Positive |
PAM50 subtype | Luminal B | Basal | HER2 | Basal | Basal | |
Clinical trial arm | FEC100 + bevacizumab | FEC100 + bevacizumab | FEC100 only | FEC100 only | FEC100 + bevacizumab | |
Week 0 | 3,987 | 14,940 | 11,449.5b | 1,934 | 3,327b | |
Tumor volume (mm3) | Week 1 | 3,392 | — | 11,460 | 1,279 | 2,115b |
Week 12 | 1,662b | 89 | 11,780 | 0 | 0 | |
Response at week 12 | NR | CR | NR | CR | CR |
Patient ID . | 1 . | 2 . | 3 . | 4 . | 5 . | |
---|---|---|---|---|---|---|
Age | 48 | 31 | 40 | 41 | 53 | |
Histologic tumor type | IDCa | IDC | IDC | IDC | IDC | |
Histologic tumor grade | 2 | 3 | 3 | 2–3 | 3 | |
Hormone receptor status | ER | Positive (100%) | Positive (1%) | Negative | Negative | Negative |
PR | Positive | Negative | Negative | Positive | Negative | |
Lymph node (LN) status | Week 0 | Positive | Negative | Positive | Negativec | Positive |
PAM50 subtype | Luminal B | Basal | HER2 | Basal | Basal | |
Clinical trial arm | FEC100 + bevacizumab | FEC100 + bevacizumab | FEC100 only | FEC100 only | FEC100 + bevacizumab | |
Week 0 | 3,987 | 14,940 | 11,449.5b | 1,934 | 3,327b | |
Tumor volume (mm3) | Week 1 | 3,392 | — | 11,460 | 1,279 | 2,115b |
Week 12 | 1,662b | 89 | 11,780 | 0 | 0 | |
Response at week 12 | NR | CR | NR | CR | CR |
NOTE: Baseline characteristics are shown together with assigned clinical trial arms and tumor volumes measured from DW-MRI segmentation at weeks 0, 1, and 12 after initiation of the treatment. By clinical examination and caliper measurement, at week 12, patients 1 and 3 were classified as nonresponders (NR) while patients 2 and 4 were classified as complete responders (CR).
Abbreviations: ER, estrogen receptor; PR, progesterone receptor.
aInvasive ductal carcinoma
bVolumes were unavailable due to problems such as fat suppression and instead computed from DCE segmentation. As patient 1 was treated with bevacizumab, which typically reduces DCE-MRI signal, maximum tumor diameters were used (30 mm at week 0 and 31 mm at week 12).
cA few possibly malignant cells not further classified were seen in lymph node aspirate before treatment start.
Histopathology
The histopathologic analysis was performed on needle biopsies taken from the breast tumors at week 0. See details and used images in section S2.2 of the Supplementary Material. Fiji and its cell counter plugin (27) were used for manual identification of cancer and stroma cells in each image.
MRI
Patients were examined on a 1.5 T MRI Scanner (ESPREE, Siemens) at weeks 0, 1, and 12 of the neoadjuvant treatment. See section S2.3 in Supplementary Material for the description of the imaging protocol and data, including dynamic contrast–enhanced (DCE-MRI) and diffusion-weighted (DW-MRI) imaging. For DCE-MRI analysis, an extended Tofts model was used, which included the determination of the contrast enhancement curve of the contrast agent in each individual voxel (volume 1 mm × 1 mm × 1.5 mm). DW-MRI data were analyzed using a simplified IVIM model–based analysis as described in ref. 28.
Molecular data
To estimate parameters representing subcellular processes, we assessed several molecular features for each tumor. mRNA levels of VEGFA and TP53 in the tumor samples were determined using one-color SurePrint G3 HumanGE 8 60 k Microarrays (Agilent Technologies) as described in ref. 26, available in the ArrayExpress database, accession number E-MTAB-4439. The PAM50 subtyping (29) was used to assign a subtype to each sample in the NeoAva cohort and proliferation score was derived for each tumor by computing mean expression values of the 11 proliferation-related PAM50 genes (26, 29). TP53 mutation status was determined by sequencing the entire coding region (exons 2–11), including splice junctions as described in ref. 26. Furthermore, pathway deregulation score (PDS; ref. 30) of the hypoxia-inducible factor 1-alpha (HIF1α) pathway was calculated for all 5 patients at screening and at week 12 and normalized against 50 tumor-free samples.
Mathematical model
To model the response of a cross-section of tumor tissue to a combination of chemotherapeutic and antiangiogenic drugs, we used a multiscale hybrid cellular automaton model. We combined cellular, extracellular, and intracellular dynamics and informed them by multitype, individual patient data. The model is detailed in section S1 and is illustrated in Supplementary Fig. S11 of the Supplementary Material. Briefly, we describe individual cells and cross-sections of functional blood vessels in a 2D section of tumor tissue as discrete agents on a regular grid. Cell division and death as well as blood vessel formation and removal are controlled in the cellular automata by intracellular and environmental factors, described by ordinary and partial differential equations. For instance, the concentration of each drug in the blood and in the tumor tissue is described by ordinary (pharmacokinetics) and partial (reaction-diffusion) differential equations, respectively. As FEC chemotherapeutic agents can only kill cycling cells, we considered a simple cell-cycle model for each cell and modeled the effect of low oxygen tension as delaying cell division. We took into account in each cell the effect of hypoxia on TP53 and VEGF expression using a system of ordinary differential equations. This allowed us to model the amount of VEGF produced by each cell and the inhibition of VEGF molecules by the antiangiogenic drug. To model the effect of VEGF on the tumor vasculature, we used a stochastic model, where the probability of formation or removal of vessels was influenced by the local VEGF concentration.
Tumor simulation and validation workflow
The workflow to run personalized simulations of drug response of a breast tumor portion is outlined in Fig. 1. Box A sketches the timeline of treatment application and tumor screening during the first phase of the NeoAva protocol. Tumor screening data at week 0, including MRI, histopathology, and molecular data (box B) were used to initialize and parametrize the mathematical model (box E). Model initialization and parametrization are explained in detail in section S3 and are illustrated in Supplementary Fig. S12 in the Supplementary Material: (i) patient-derived parameters were obtained directly from the clinical data (see box B in Fig. 1; Supplementary Table S2); (ii) common parameters to all studied patients were obtained from public data (see box C in Fig. 1; Supplementary Table S3); (iii) finally, there were model parameters that were calibrated because they could not be inferred from the clinical data directly, nor relevant quantities were found in the literature. Calibration was done by finding ranges of parameter values where the simulated and patient outcome qualitatively agree (see Supplementary Table S4 and Supplementary Fig. S13–S18 for parameter exploration of effects on treatment outcomes). Specifically, the exact drug schedule used in the clinical trial for each patient (box D) was simulated by increasing the amount of the corresponding drugs (FEC100 ± bevacizumab) in the blood at the time points of administration according to the dosage. Each computer simulation then ran a complete cycle of 12 weeks of the spatiotemporal dynamics of the considered tumor portion under the effect of the applied drugs. Computational details and a link to the open-source computer code can be found in section S4 of the Supplementary Material. Simulated outcomes include the spatial distribution of cells, blood vessels, and all considered molecules at any time between the start of the therapy and week 12 and in particular at the end of the study period (box F). To calibrate the three model parameters, we compared the simulated outcomes with the actual ones (box H). For that, we used clinical tumor volumes and apparent diffusion coefficients (ADC; ref. 26) calculated from DW-MR imaging at weeks 1 and 12 (box G). As we were simulating only very small portions of the tumor bulk, we only selected patients who were complete responders (number of tumor cells in all simulated tumor portions at week 12 should be around zero) or complete nonresponders at week 12 (approximately the same number of tumor cells in week 12 as in week 0). Moreover, changes in ADC histograms of the segmented tumor volume at weeks 0, 1, and 12, which relate to changes in the tumor cellularity (20), should qualitatively agree with the number of tumor cells in our simulations.
Results
Comparison of patient characteristics and treatment responses
We considered 4 patients (1–4; see Table 1 and section S2 in Supplementary Material). Patients 1 and 2 received FEC100 plus bevacizumab with the same dose and schedule. Patient 1 did not respond to the therapy after 12 weeks, while patient 2 responded well. At baseline, patient 2 had a higher tumor cell density than patient 1, as shown by histologic images and DW-MRI. This coincided with higher PDS of HIF1α pathway. Their tumors were also classified differently as basal-like and luminal B, respectively, and the estimated proliferative capacity of patient 2 was much higher (see parameter Tmin in section S3.5.1 of the Supplementary Material), in agreement with the estimated PAM50 proliferation scores (31) of both patients. Other molecular and MRI-deduced parameters used in this study did not differ significantly. Specifically, expression levels of VEGF and TP53 were almost identical, and both tumors were TP53 wild-type. The typical values of perfusion parameters ktrans and vp estimated from DCE-MRI were similar too.
Patients 3 and 4 received only FEC100 with the same dose and schedule. Patient 3 did not respond to the therapy after 12 weeks, while patient 4 responded completely. The tumors were classified as HER2-enriched and basal-like, respectively, but they had similar proliferative capacity as shown by their estimated PAM50 proliferation scores. They differed in tumor morphology, vessel perfusion, and TP53 status and expression. Comparing histologic slices and DW-MRI, the tissue of patient 3 was densely packed with cells, while patient 4 showed more heterogeneity. DW-MRI of patient 4 also showed heterogeneity. At week 0, by DCE-MRI, the tumor in patient 3 had a very poorly perfused core while it was highly perfused on the outer edge, resulting in cross-sections with ring-like patterns. Interestingly, DW-MRI analysis suggested that the tumor core was not necrotic. The tumor from patient 3 was TP53 wild-type while the tumor from patient 4 was TP53-mutated. Expression levels of TP53 were higher in patient 4 while VEGF expression was comparable and very high in both patients. Patient 3 had much higher HIF1α PDS at week 0, exhibiting signs of a denser and more hypoxic tumor.
Model simulations reproduce treatment outcomes
For each of the 4 patients, we ran personalized simulations of tumor portions under the treatment received in the clinical trial as outlined in Fig. 1 and described in the Supplementary Material.
Low versus highly proliferative tumor in bevacizumab arm
Model simulations of the outcomes of patients 1 and 2 are shown in Fig. 2. We used two biopsy portions for each patient (biopsy A and B). For each portion, we ran ten independent stochastic simulations of the 12-week treatment using the same parameterization, but different random events such as births and deaths of vessels. We plotted the time evolution of the cancer cell density, defined as the proportion of grid points occupied by cancer cells at any given time, for the ten simulations. Simulated cell densities for patient 1, shown in Fig. 2A, decreased moderately after drug administration and grew in the period between consecutive administrations. The four grid snapshots show, at different time points, the spatial distribution of cancer and stroma cells together with the oxygen level in one representative simulation. All 20 simulated experiments presented moderate degrees of hypoxia, as seen by the white and light blue background displayed in the snapshots of Fig. 2A. This is in agreement with the observed hypoxic pathway activity and VEGF expression, shown in Supplementary Fig. S10B and S10D. In fact, despite the production of VEGF by hypoxic cells, the applied dose of bevacizumab reduced VEGF concentration in the tumor tissue to very low levels, as recorded above each snapshot in Fig. 2A. In our model, blood vessels disappeared with a certain probability under low VEGF concentration. For patient 1, if this probability is high enough, the tumor can persist, because bevacizumab reduces tumor perfusion and the chemotherapeutic agents cannot reach the tumor tissue. We could reproduce the outcome at week 12 using any value for the probability for vessels to become dysfunctional under low VEGF and an intermediate value of the chemosensitivity parameter, as shown in Supplementary Figs. S13 S16. In Fig. 2A, even when using a very low probability of vessel death (Pdeath = 0.0001), the killing of cancer cells after each dose administration did not prevent cancer cells from proliferating again, resulting in a net balance between killing and proliferation in 20 of 20 simulations. This is in agreement with the available ADC data of patient 1 at week 0 and 1, showing little difference in tumor density (Supplementary Fig. S9A). By the end of the simulations, the cancer cell density was approximately the same as it was at the beginning.
A, Simulations of NeoAva drug schedule in patient 1. B, Simulations of NeoAva drug schedule in patient 2. Simulated time evolution of the cancer cell density starting from two different cell configurations at week 0 corresponding to biopsy portions A and B. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Bottom, the spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for each patient. Background color represents oxygen pressure in mm Hg.
A, Simulations of NeoAva drug schedule in patient 1. B, Simulations of NeoAva drug schedule in patient 2. Simulated time evolution of the cancer cell density starting from two different cell configurations at week 0 corresponding to biopsy portions A and B. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Bottom, the spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for each patient. Background color represents oxygen pressure in mm Hg.
Patient 2 was a responder. Many blood vessels remained in the tissue to allow chemotherapeutic drugs to be distributed and thus killed most of the tumor cells. To reproduce the available ADC data, where virtually no tumor remained at week 12 (Supplementary Fig. S9B), the following was required: first, a low probability of vessel disappearing (Pdeath < 0.01), so that part of the tumor vasculature was resistant to the VEGF inhibitor; second, a chemosensitivity high enough (β > 6000), so that few cells escaped the therapy (see Supplementary Figs. S13 and S16). In the simulations of patient 2 shown in Fig. 2B, all tumor cells were killed after 12 weeks in 19 of 20 experiments. Because of higher cell densities at week 0, hypoxia was more severe than for patient 1, in agreement with the HIF pathway (Supplementary Fig. S10D) and ADC data (Supplementary Fig. S9B).
Heterogeneous perfusion condition without bevacizumab
Patient 3 was a nonresponder. As revealed by MRI data (Supplementary Fig. S7), the core of the tumor was poorly perfused, with much higher perfusion at the tumor edge. Because the location of our biopsy was not identifiable, we show simulations for both perfusion profiles in Fig. 3. ktrans = 0.0067min−1 and vp = 1.73% were estimated from the core to reflect its poorly perfused condition (Fig. 3A). Because the number of vessels was small and their permeability was low, the simulated drug concentration arriving in the tissue was very low. Therefore, free spaces in the simulation grid were occupied by the cancer cells (see Fig. 3A). This is in agreement with DW-MRI data (Supplementary Fig. S9C) and HIF1α PDS that explain patient 3 being a nonresponder. At the highly perfused tumor edge, ktrans and vp were estimated to be 0.2145/minute and 15.35%, respectively, and drugs arrived in the tissue more efficiently. Provided cell chemosensitivity is high enough, simulated cancer cell densities can be reduced (as shown in Fig. 3B). Although cells in the edge could be killed, this area could have been repopulated by the highly proliferative cells from the core.
A, Simulations of NeoAva drug schedule in patient 3 under poorly perfused condition. B, Simulations of NeoAva drug schedule in patient 3 under well-perfused condition. Time evolution of cancer cell density for two different cell configurations, biopsy portions A and B, and drug schedules under two different perfusion profiles. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Panel A represents the core of the tumor of patient 3, while the bottom panel of B shows the peripheral spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for each patient. Background color represents oxygen pressure in mm Hg.
A, Simulations of NeoAva drug schedule in patient 3 under poorly perfused condition. B, Simulations of NeoAva drug schedule in patient 3 under well-perfused condition. Time evolution of cancer cell density for two different cell configurations, biopsy portions A and B, and drug schedules under two different perfusion profiles. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Panel A represents the core of the tumor of patient 3, while the bottom panel of B shows the peripheral spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for each patient. Background color represents oxygen pressure in mm Hg.
Patient 4 was a responder with also a heterogeneous perfusion profile as seen by DCE-MRI. Contrary to patient 3, a scattered pattern with two representative perfusion profiles was observed, one highly perfused with low permeability, the second with lower perfusion but high permeability. We show simulations of the first perfusion profile in Fig. 4A, with the estimated parameters vp = 6.0% and ktrans = 0.0067/minute. Simulations of the second profile are in Fig. 4B, where we estimated vp = 3.99% and ktrans = 0.1107min−1. To reproduce the DW-MRI data in both cases, where we observed that the tumor was drastically reduced at week 1 and almost disappeared by week 12 (Supplementary Fig. S9D in the SM), a higher rate of vessel creation compared with the other 3 patients was required (Supplementary Figs. S14 and S15B). This was achieved by increasing pbirth and HighV (see Supplementary Table S4). This allowed for transport of the drug efficiently to the tissue, even when the permeability was very small. The simulated reduction of severe hypoxia to normoxia shown in the snapshots of Fig. 4 is also in agreement with the observed significant decrease in HIF1α PDS, comparing week 0 with week 12 (Supplementary Fig. S10D).
A, Simulations of NeoAva drug schedule in patient 4 under perfusion profile 1. B, Simulations of NeoAva drug schedule in patient 4 under perfusion profile 2. Time evolution of cancer cell density for two different cell configurations, biopsy portions A and B, and drug schedules under two different perfusion profiles. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Bottom panels of A and B show the spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for the patient. Background color represents oxygen pressure in mm Hg.
A, Simulations of NeoAva drug schedule in patient 4 under perfusion profile 1. B, Simulations of NeoAva drug schedule in patient 4 under perfusion profile 2. Time evolution of cancer cell density for two different cell configurations, biopsy portions A and B, and drug schedules under two different perfusion profiles. Each line represents the average cancer cell density of 10 independent stochastic simulations. The corresponding color band indicates the 95% bootstrap confidence interval. Bottom panels of A and B show the spatial distribution of cancer and stroma cells for a representative simulation of the biopsy portion B for the patient. Background color represents oxygen pressure in mm Hg.
Simulating alternative drug regimens
In addition to the regimens used in the clinical study, we simulated alternative drug regimens (schedules and doses) for the two nonresponders.
Low dosage frequent chemotherapy dosing for slowly proliferative tumor
The simulations of patient 1 in Fig. 2A suggest that a main reason behind the survival of cancer cells after each drug administration was their low proliferation rate. As the interplay of treatment frequency and cell proliferation rate can contribute to the outcome, we hypothesized that administering a smaller dose more frequently while keeping the overall amount would be beneficial. We therefore simulated four chemotherapy regimens: every week with one third of the original dose, every one and a half week with half of the original dose, every 2 weeks with two thirds of the original dose, and the original 3-week schedule used in the clinical trial. In Fig. 5A, we compare the outcomes after 12 weeks of treatment in terms of change in cell density. We used the same initialization and parameters of the personalized simulations of patient 1. To see the effect of cell proliferation rate, we repeated the same simulations while changing the cell-cycle length parameter over a wide range. Simulation results showed that the interplay between cell-cycle length and schedule is complex and can exhibit nonmonotonic behaviors. Interestingly, we observed that drug administrations every week or every week and a half can improve the treatment outcome of patient 1 (marked with a asterisk in Fig. 5A). In Fig. 5B, we show simulations of the most successful schedule for patient 1, where we reduced FEC100 and bevacizumab to a third of their original dose, but administered every week instead of every 3 weeks. The tissue was free of cancer cells after approximately 6 weeks.
A, Effects of different chemotherapy schedules on patient 1 and effect of cell proliferation. B, Simulations of an alternative schedule in patient 1. In A, each colored dot represents the difference in cancer cell density between start and end of therapy under a specific drug schedule and minimal cell-cycle length. Tmin of daughter cells was randomized to introduce variability in cell-cycle duration. Colored lines represent the locally weighted smoothed curve smoothing fitting of the simulations of different cell-cycle length under the same drug schedule. The simulation of patient 1 corresponding to the value Tmin = 14.69 is labeled with a star. B, The temporal dynamics of the alternative drug schedule providing better outcome. Solid lines represent the average cell density (n = 10) of the alternative therapy obtained by reducing FEC100 and bevacizumab to a third of its original dose and administrating them every week instead of every third week.
A, Effects of different chemotherapy schedules on patient 1 and effect of cell proliferation. B, Simulations of an alternative schedule in patient 1. In A, each colored dot represents the difference in cancer cell density between start and end of therapy under a specific drug schedule and minimal cell-cycle length. Tmin of daughter cells was randomized to introduce variability in cell-cycle duration. Colored lines represent the locally weighted smoothed curve smoothing fitting of the simulations of different cell-cycle length under the same drug schedule. The simulation of patient 1 corresponding to the value Tmin = 14.69 is labeled with a star. B, The temporal dynamics of the alternative drug schedule providing better outcome. Solid lines represent the average cell density (n = 10) of the alternative therapy obtained by reducing FEC100 and bevacizumab to a third of its original dose and administrating them every week instead of every third week.
Low dosage bevacizumab for poorly perfusive tumor
Patient 3 was a nonresponder in the trial arm with chemotherapy only. Our analysis suggests that the main reason for the negative outcome was that the drugs did not penetrate enough in the tumor core. We hypothesized that by administering an appropriate amount of bevacizumab, VEGF levels could be reduced appropriately, and tumor perfusion could be improved with the creation of new functional vessels in the core (32, 33). Thereafter, chemotherapy could be delivered efficiently. We added different bevacizumab regimens to the personalized simulations of patient 3, using the same initialization and parameters. Moreover, to see the effect of VEGF expression on the outcome, we repeated the same simulations while changing the parameters modulating VEGF expression levels from low to high (Fig. 6A). Because bevacizumab has a long half-life, the same schedule was used for the chemotherapy (every 3 weeks) while scaling the amount with respect to what was originally administered in the other arm of the trial. As seen in Fig. 6A, applying full dose of bevacizumab led to no improvement compared with chemotherapy only and the simulated outcomes strongly depended on the applied bevacizumab regime and the VEGF expression. For patient 3 (VEGF level = 2.32), tumor burden was reduced with doses lower than a fourth of the full bevacizumab dose. In Fig. 6B, we show simulations using the optimal bevacizumab regime for patient 3. Comparing this simulation with its actual clinical outcome, the new therapy appeared to improve the outcome by inhibiting the growth of the tumor, reducing the density of the tumor by 50% on average.
A, Effects of different bevacizumab schedule on patient 3 with different estimated VEGF expression. B, Simulation of an alternative drug schedule for patient 3. In A, the y-axis represents the simulated difference in cancer cell density between week 0 and week 12 under a 3-week administration interval. The x-axis indicates the hypothetical quantity of bevacizumab administered to the patient as fractions with respect to the full dose, corresponding to 15 mg/kg. Each colored band summarizes 10 independent simulations for patient 3, with VEGF expression from low to high given in the legend. These four expression levels correspond to those observed in our 4 patients. However, all other parameters are fixed as for patient 3. In B, we showed temporal dynamics of an alternative drug schedule providing better outcome. Solid lines represent the average cell density (n = 10) of the alternative therapy obtained by administrating a reduced bevacizumab dosage of 0.5 mg per kg body weight, equivalent of 3.33% of the dosage in the experimental arm alongside of chemotherapy following a 3-week interval. The corresponding simulation is labeled with a star in A.
A, Effects of different bevacizumab schedule on patient 3 with different estimated VEGF expression. B, Simulation of an alternative drug schedule for patient 3. In A, the y-axis represents the simulated difference in cancer cell density between week 0 and week 12 under a 3-week administration interval. The x-axis indicates the hypothetical quantity of bevacizumab administered to the patient as fractions with respect to the full dose, corresponding to 15 mg/kg. Each colored band summarizes 10 independent simulations for patient 3, with VEGF expression from low to high given in the legend. These four expression levels correspond to those observed in our 4 patients. However, all other parameters are fixed as for patient 3. In B, we showed temporal dynamics of an alternative drug schedule providing better outcome. Solid lines represent the average cell density (n = 10) of the alternative therapy obtained by administrating a reduced bevacizumab dosage of 0.5 mg per kg body weight, equivalent of 3.33% of the dosage in the experimental arm alongside of chemotherapy following a 3-week interval. The corresponding simulation is labeled with a star in A.
Other alternative treatments
We tested several regimens removing fluorouracil from the treatment. We found no difference in patient's outcome compared with the FEC100 regime administered every 3 and every 2 weeks. Interestingly, this was in agreement with the findings in a clinical study (see Supplementary Fig. S19; ref. 34).
We investigated the effect of prolonging the 12-week FEC regimens by continuing to administer the same chemotherapy and bevacizumab for another 12 weeks for patient 1. We saw no benefit in terms of patient's treatment outcome. This confirmed the appropriateness of the actual clinical decision of switching to taxane-based chemotherapy after week 12 (Supplementary Fig. S20).
Predicting treatment outcome
Here we predict the response to treatment of a new patient using our model. We chose a patient, patient 5, with the same subtype and in the same trial arm as one of the other 4 patients, namely patient 2. Screening data were then collected as we did for the other patients and as described in section S3 of the Supplementary Material. The simulations were initialized with parameters estimated from the same patient 5 or from the literature, but we did not calibrate any parameter. Instead we used exactly the same chemosensitivity, probability of vessel birth and death, and the VEGF thresholds as for patient 2. At week 0, patient 5 had a moderately sized tumor with necrotic region shown on DCE and DW-MR images. Both biopsies showed a high cancer cell density. As a result of this, HIF pathway score was high, comparable with patient 2, indicating hypoxia.
For each of the two digitally captured biopsy images, we ran 20 stochastic simulations and obtained a bootstrap confidence interval for the cell density at the end of week 12. At the end of week 12, the simulated tumor disappeared, and HIF PDS was close to 0. The two-sided 95% bootstrapped confidence interval for biopsies A and B were 0–0.0154 and 0–0.0049, respectively. This was true, as patient 5 experienced a clinical complete response.
Discussion
This study is a new step toward the simulation of individualized tumor response to therapy. We have shown that it is possible to design a multiscale mathematical model, which integrates different types of individual patient's data collected in clinical practice and simulates the observed patient outcome. Our model is flexible enough to simulate different outcomes based on individual data by capturing fundamental biological processes at an acceptable level of approximation. Importantly, our model suggests possible mechanistic explanations of individual treatment outcomes and allows virtual testing of alternative treatment plans. For example, the administered amount of bevacizumab did not benefit patients. Simulations show that moderate amounts of bevacizumab can improve the outcome of patients with extreme hypoxia. As for FEC, we found that for highly proliferative tumors, cell-cycle–specific drugs are very effective, while for low proliferative tumors, a more frequent but lower dosing of FEC can be advantageous.
We are well equipped to recognize limitations of our approach that need to be addressed in the future. First, we have only simulated a few tumor sections of 200 × 300 μm, which might not be representative of the whole tumor bulk. Although the current computational power allows multiscale model simulations of larger tumor portions, the approach is still limited by the availability of clinical data to inform them. For instance, clinical information at the cellular level, such as the number, position, and type of cells, is currently possible only for biopsies. We are extending our algorithms to run cross-sections of full biopsies. A second limitation of our model relates to tumor heterogeneity. Some cell clones can be more proliferative than others, produce more VEGF, or be more resistant to therapy. Vessels can be different in size, have unequal functionality and permeability, and contain different oxygen and drug concentrations. Because of the current impossibility of characterizing cell and vessel heterogeneity from the available clinical data, our simulations assume all cells and vessels to be of the same type. Extending our model to multiple clones, competing for resources, would rely on a deeper picture of the patients’ tumor than routinely available today.
Third, we have assumed that model parameters are constant during the simulated treatment. For instance, we used constant drug chemosensitivity, which is a simplification and neglects the possible evolution of resistant phenotypes. To incorporate such details in the current approach, monitoring of tumor evolution is needed.
We incorporated some degree of heterogeneity in our studies by running simulations with different cell configurations (observed in the biopsies), and different perfusion and permeabilities (observed in MR images of the whole tumor). We account for spatiotemporal heterogeneity in the simulated tissue section where cell and vessel numbers, as well as oxygen, VEGF, and drug concentrations change in space and time.
Fourth, VEGF is the only angiogenic factor present in our model and does not explicitly reflect the redundancy of angiogenic stimulants and inhibitors known to be present in a tumor. The mode of action of the VEGF depletion treatment was assumed to be a reduction in the number of functional capillaries, resulting in local reduced supply of oxygen and treatment. However, the effect may be more complex, related to the perfusion and leakage of the tumor influenced by the VEGF, as seen in refs. 26 and 35. See also ref. 22 for a different susceptibility of VEGF withdrawal dependent on the immune activation in the tumor. Future version of the model should incorporate first principles of such mechanisms.
An important contribution of our study is the identification of three parameters that cannot be estimated precisely enough from the present clinical data: the chemosensitivity of tumor cells and the two parameters encoding the sensitivity of vessels to the local VEGF concentration. In our feasibility study, we calibrated these parameters, within a range of realistic values, by choosing the ones that allow simulation of a trajectory compatible with the true patient endpoint. A systematic sensitivity analysis is beyond the scope of this study. We intend to address it in the future. We were able to predict the treatment outcome for 1 patient without calibration. More tests are needed.
There are several ideas on how the three parameters could be estimated in the future. Drug chemosensitivities can be measured using ex vivo patient material; sensitivity of vessels can be perhaps estimated using in vivo experiments on mice xenografts (35, 36). Another possibility is to use approximate Bayesian computation (37) to estimate individual parameters from observing the first cycle of treatment. Preliminary results in this direction are promising.
We highlight the importance of rich longitudinal data to improve the accuracy of the simulated therapies. Repeated MR images at different time points matched to each other (11) would allow to relate tumor features such as volume, perfusion, and cell density along therapy. In addition, matching biopsies and MRI data would be crucial to locate the exact location of the extracted biopsy in the MR images, allowing more accurate estimation of parameters such as the vessel permeability. Furthermore, instead of using average drug pharmacokinetic models, individual longitudinal measurements of drug concentrations in the blood could be easily performed, allowing further personalization (7).
In this study, we used some genomic data, but much more should be done. It is known that mutations, gene expression, and copy number variations portray properties of tumor cells, some of which may imply treatment success. For example, the PAM50 gene expression signature (29) allows classification of breast cancer tumors into five distinct classes. As these genomic features modulate the efficacy of therapy, our study suggests that the key parameters in our model, proliferative capacity and chemosensitivity of tumor cells, should depend on such genomic features of the patient.
This study opens the possibility to build models that allow the simulation of other solid tumors. This would require extensive work on the mathematical modeling and inferential side, and the identification of informative clinical data. In conclusion, our work shows realistic possibilities for simulation-guided personalized therapy and indicates where further research should focus to make this possible.
Disclosure of Potential Conflicts of Interest
A.-L. Børresen-Dale has ownership interest (including stocks and patents) in Arctic Pharma AS and is a consultant/advisory board member for Arctic Pharma AS, Saga Diagnostic AS, and PubGene AS. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: X. Lai, T. Seierstad, A.-L. Børresen-Dale, V.N. Kristensen, O. Engebråten, A. Köhn-Luque, A. Frigessi
Development of methodology: X. Lai, O.M. Geier, S.W. Funke, M.E. Rognes, T. Seierstad, A. Köhn-Luque, A. Frigessi
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): O.M. Geier, Ø. Garred, E. Borgen, T. Seierstad, A.-L. Børresen-Dale, V.N. Kristensen, O. Engebråten
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): X. Lai, O.M. Geier, T. Fleischer, S. Kumar, T. Seierstad, V.N. Kristensen, O. Engebråten, A. Köhn-Luque, A. Frigessi
Writing, review, and/or revision of the manuscript: X. Lai, O.M. Geier, T. Fleischer, Ø. Garred, E. Borgen, S.W. Funke, M.E. Rognes, T. Seierstad, A.-L. Børresen-Dale, V.N. Kristensen, O. Engebråten, A. Köhn-Luque, A. Frigessi
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): O.M. Geier, V.N. Kristensen, O. Engebråten, A. Frigessi
Study supervision: S.W. Funke, V.N. Kristensen, O. Engebråten, A. Köhn-Luque, A. Frigessi
Acknowledgments
A. Köhn-Luque would like to acknowledge the funding from the European Union Seventh Framework Programme (FP7-PEOPLE-2013-COFUND) under grant agreement number 609020 - Scientia Fellows, and M.E. Rognes would like to acknowledge the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement 714892 for its support. We also would like to acknowledge the help received from the Department for Research Computing at USIT, the University of Oslo IT-department.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.