Towards personalized computer simulations of breast cancer treatment
VVPH2020 Conference – Paris 26-28 August 2020
Towards personalized computer simulation of breastcancer treatment
Alvaro K ¨ohn-Luque a ∗ , Xiaoran Lai a and Arnoldo Frigessi a a Oslo Centre for Biostatistics and Epidemiology, University of Oslo, Norway,
Keywords multiscale modeling, hybrid cellular automaton, pharmacokinetics, pharmacodynamics, parallelcomputing, approximate Bayesian computation
1. Introduction
Cancer pathology is unique to a given individual, anddeveloping personalized diagnostic and treatment pro-tocols are a primary concern. Mathematical model-ing and simulation is a promising approach to person-alized cancer medicine. Yet, the complexity, hetero-geneity and multiscale nature of cancer present se-vere challenges. One of the major barriers to usemathematical models to predict the outcome of ther-apeutic regimens in a particular patient lies in their ini-tialization and parameterization in order to reflect indi-vidual cancer characteristics accurately.Here we present a study where we used multitypemeasurements acquired routinely on a single breasttumor, including histopathology, magnetic resonanceimaging (MRI), and molecular profiling, to personalizea multiscale hybrid cellular automaton model of breastcancer treated with chemotherapeutic and antiangio-genic agents. We model drug pharmacokinetics andpharmacodynamics at the tumor tissue level but withcellular and subcellular resolution. We simulate thosespatio-temporal dynamics in 2D cross-sections of tu-mor portions over 12-week therapy regimes, resultingin complex and computationally intensive simulations.For such computationally demanding systems, for-mal statistical inference methods to estimate individ-ual parameters from data have not been feasible inpractice to until most recently, after the emergenceof machine learning techniques applied to likelihood-free inference methods. Here we use the inferenceadvances provided by Bayesian optimization to fit ourmodel to simulated data of individual patients. In thisway, we investigate if some key parameters can be es-timated from a series of measurements of cell densityin the tumor tissue, as well as how often the measure-ments need to be taken to allow reliable predictions.
2. Methods
To develop and validate our mathematical model weused data from HER2-negative mammary carcinomasfrom the NeoAva cohort [3], a randomized, phase IIclinical trial that evaluated the effect of bevacizumab ∗ Corresponding author. Email: [email protected] in combination with neoadjuvant treatment regimes.Four patients were selected to belong to both arms ofthe trial and to have either a complete or no responseat 12 weeks of treatment. Details of baseline char-acteristics, treatment, response and used clinical datafrom all patients can be found in reference [2].
We model the response of a cross-section of tumortissue to a combination of chemotherapeutic and an-tiangiogenic drugs using a multiscale hybrid cellularautomaton. We combined tissue, cellular, extracel-lular, and intracellular dynamics and informed themby multitype, individual patient data. Figure 1 showsan overview of the model structure and the differentmodel formalisms used. Full details are available inreference [2].
Avastin FEC
Model Formalism
Stochastic Cellular Automaton Partial Differential Equations Stochastic Cellular Automaton Ordinary Differential Equations
ODEs
Oxygen Avastin FECCancer CellStroma Cell
Cell cycle p53VEGFVEGF p53 Blood VesselsVEGF
Model Modules
Vascular Module Extravacular-Extracellular Module Cellular Module Subcellular Module Intravascular Module
Ordinary Differential Equations
Model Components
Figure 1: Multiscale hybrid cellular automaton struc-ture
For parameter inference we use an approximateBayesian computation method called Bayesian op-timization for likelihood-free inference (BOLFI) [1].BOLFI uses Gaussian processes to model the un-known discrepancy function that compares the simu-lations with the data via Bayesian optimization. Oncethe probabilistic model has been inferred, it can beused to approximate the true posterior distribution.
The open-source code to simulate the multiscalemathematical model is written in Python 3.6.2 and a r X i v : . [ q - b i o . T O ] J u l PH2020 Conference – Paris 26-28 August 2020available at bitbucket.org/xlai/codeavastin-py.git
3. Results
We first use our model to reproduce and elucidate thetreatment outcome of four patients. For each of them,we run personalized simulations of different tumor por-tions under the treatment received in the trial. Tumorscreening data at baseline (histopathology, MRI andmolecular data) are used to initialize and parameterizethe model together with other data available in the lit-erature. We simulate 12 weeks of the spatio-temporaldynamics of the considered tumor portion under theeffect of the applied drugs. We then use MRI examina-tions at weeks 1 and 12 to compare with the simulatedoutcomes. As an example, we illustrate in Figure 2-I asimulation of a non-responder patient with low prolif-erative cells able to scape the standard chemotherapyregime administered every three weeks.
Mathematical Oncology Cancer Research
Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (a)
Simulations of NeoAva drug schedule in patient 1 . Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (b)
Simulations of NeoAva drug schedule in patient 2.
Figure 2: Simulated time evolution of cancer cell density starting from di↵erent cell configura-tions at week 0. Each line represents the average cancer cell density of 10 simulations, di↵eringonly in variables randomly drawn with certain probabilities, such as initialization of cell cycleand killing in the model. The corresponding color band indicates the 95% bootstrap confidenceinterval. Panel B of each figure shows the spatial distribution of cancer and stroma cells for arepresentative simulation of the second cell configuration for each patient. Background colorrepresents oxygen concentration in mmHg.Patient 3 is a non-responder. As we saw from DCE-MRI data, fig. 15 in the SupplementaryMaterial, the core of the tumor of this patient is poorly perfused, while the border is wellperfused, as seen from DCE-MRI data. Therefore we show, in fig. 3, model simulations thatreproduce the outcome of patient 3, representing a heterogenous perfusion profile.We used the estimated values k trans = 0 . / min and v p = 1 .
73% to reflect poorly-11
Mathematical Oncology Cancer Research
Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (a)
Simulations of NeoAva drug schedule in patient 1 . Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (b)
Simulations of NeoAva drug schedule in patient 2.
Figure 2: Simulated time evolution of cancer cell density starting from di↵erent cell configura-tions at week 0. Each line represents the average cancer cell density of 10 simulations, di↵eringonly in variables randomly drawn with certain probabilities, such as initialization of cell cycleand killing in the model. The corresponding color band indicates the 95% bootstrap confidenceinterval. Panel B of each figure shows the spatial distribution of cancer and stroma cells for arepresentative simulation of the second cell configuration for each patient. Background colorrepresents oxygen concentration in mmHg.Patient 3 is a non-responder. As we saw from DCE-MRI data, fig. 15 in the SupplementaryMaterial, the core of the tumor of this patient is poorly perfused, while the border is wellperfused, as seen from DCE-MRI data. Therefore we show, in fig. 3, model simulations thatreproduce the outcome of patient 3, representing a heterogenous perfusion profile.We used the estimated values k trans = 0 . / min and v p = 1 .
73% to reflect poorly-11
Mathematical Oncology Cancer Research
Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (a)
Simulations of NeoAva drug schedule in patient 1 . Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (b)
Simulations of NeoAva drug schedule in patient 2.
Figure 2: Simulated time evolution of cancer cell density starting from di↵erent cell configura-tions at week 0. Each line represents the average cancer cell density of 10 simulations, di↵eringonly in variables randomly drawn with certain probabilities, such as initialization of cell cycleand killing in the model. The corresponding color band indicates the 95% bootstrap confidenceinterval. Panel B of each figure shows the spatial distribution of cancer and stroma cells for arepresentative simulation of the second cell configuration for each patient. Background colorrepresents oxygen concentration in mmHg.Patient 3 is a non-responder. As we saw from DCE-MRI data, fig. 15 in the SupplementaryMaterial, the core of the tumor of this patient is poorly perfused, while the border is wellperfused, as seen from DCE-MRI data. Therefore we show, in fig. 3, model simulations thatreproduce the outcome of patient 3, representing a heterogenous perfusion profile.We used the estimated values k trans = 0 . / min and v p = 1 .
73% to reflect poorly-11
Mathematical Oncology Cancer Research
Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (a)
Simulations of NeoAva drug schedule in patient 1 . Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (b)
Simulations of NeoAva drug schedule in patient 2.
Figure 2: Simulated time evolution of cancer cell density starting from di↵erent cell configura-tions at week 0. Each line represents the average cancer cell density of 10 simulations, di↵eringonly in variables randomly drawn with certain probabilities, such as initialization of cell cycleand killing in the model. The corresponding color band indicates the 95% bootstrap confidenceinterval. Panel B of each figure shows the spatial distribution of cancer and stroma cells for arepresentative simulation of the second cell configuration for each patient. Background colorrepresents oxygen concentration in mmHg.Patient 3 is a non-responder. As we saw from DCE-MRI data, fig. 15 in the SupplementaryMaterial, the core of the tumor of this patient is poorly perfused, while the border is wellperfused, as seen from DCE-MRI data. Therefore we show, in fig. 3, model simulations thatreproduce the outcome of patient 3, representing a heterogenous perfusion profile.We used the estimated values k trans = 0 . / min and v p = 1 .
73% to reflect poorly-11
VEGF max = − µ gmL − VEGF max = − µ gmL − VEGF max = − µ gmL − VEGF max = − µ gmL − Cells
CancerStroma05101520
O2/mmHg
Week C an c e r C e ll D en s i t y Administration
FEC+Avastin
Biopsy Cell densityCase ACase B (I)(II)
Mathematical Oncology Cancer Research
Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (a)
Simulations of NeoAva drug schedule in patient 1 . Week C an c e r C e ll D en s i t y Biopsy AB Administration
FEC+Avastin A VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − VEGF max = x − µ g mL − B Cells
CancerNormal05101520
O2/mmHg (b)
Simulations of NeoAva drug schedule in patient 2.
Figure 2: Simulated time evolution of cancer cell density starting from di↵erent cell configura-tions at week 0. Each line represents the average cancer cell density of 10 simulations, di↵eringonly in variables randomly drawn with certain probabilities, such as initialization of cell cycleand killing in the model. The corresponding color band indicates the 95% bootstrap confidenceinterval. Panel B of each figure shows the spatial distribution of cancer and stroma cells for arepresentative simulation of the second cell configuration for each patient. Background colorrepresents oxygen concentration in mmHg.Patient 3 is a non-responder. As we saw from DCE-MRI data, fig. 15 in the SupplementaryMaterial, the core of the tumor of this patient is poorly perfused, while the border is wellperfused, as seen from DCE-MRI data. Therefore we show, in fig. 3, model simulations thatreproduce the outcome of patient 3, representing a heterogenous perfusion profile.We used the estimated values k trans = 0 . / min and v p = 1 .
73% to reflect poorly-11
AdministrationSingle cellsCancerStromaFEC + AvastinO2 mmHg
Figure 2: Personalized treatment simulations of onepatient: (I) Failed regime received in the trial. (II) Al-ternative regime with improved outcome. Simulatedevolution of the cancer cell density starting from twoobserved cell configurations at week 0, case A andcase B. Each line represents the average cancer celldensity of 10 independent stochastic simulations withtheir 95% bootstrap confidence interval. Lower panelof each figure shows the spatial distribution of can-cer and stroma cells for a representative simulationof Case B. Background color represents oxygen pres-sure in mmHg.
In addition to the regimes used in the clinicalstudy, we also tested multiple hypothetical alterna-tive drug regimes (schedules and doses) for the non-responders. Figure 2-II shows a simulation of themost successful schedule obtained for the same non- responder patient shown in Figure 2-I. In that casethe drugs doses are reduced to a third of its orig-inal ones, but administrated every week instead ofevery 3 weeks. By means of the higher drug fre-quency, low proliferative cells can no longer escapethe chemotherapy and the tissue becomes free of can-cer cells after approximately 6 weeks of therapy.
Using the simulated number of cells from the patientsas measurements, BOLFI is able to accurately inferup to three of the most relevant model parameters. Infact, parameters were estimated well enough by usingonly a few measurements made in the very beginningof the therapy. Importantly, this was sufficient to pre-dict the effect of the therapy regime at the end of the12 weeks with good accuracy.
4. Conclusion
The use of mathematical models and simulations todecide the best therapy for a cancer patient starts tobecome a possibility. For breast cancer, we have beenable to develop and implement a first multiscale phar-macokinetic and pharmacodynamics model, whichshows encouraging pilot runs on a few patients. Anysuch model needs to be informed by data from theunique patient for whom we wish to determine the besttherapeutic regime. We have shown that it is possibleto estimate personalized parameters with data com-patible with clinical practice and use them to make re-liable predictions of the effect of therapy for a uniqueand single patient.
5. References [1] Michael U Gutmann, Jukka Corander, et al.Bayesian optimization for likelihood-free inferenceof simulator-based statistical models.
Journal ofMachine Learning Research , 2016.[2] Xiaoran Lai, Oliver M Geier, Thomas Fleischer,Øystein Garred, Elin Borgen, Simon W. Funke,Surendra Kumar, Marie E. Rognes, Therese Seier-stad, Anne-Lise Børressen-Dale, Vessela N. Kris-tensen, Olav Engebr ˚aten, Alvaro K ¨ohn-Luque,and Arnoldo Frigessi. Towards personalizedcomputer simulation of breast cancer treatment:a multi-scale pharmacokinetic and pharmacody-namic model informed by multi-type patient data.
Cancer Research , 2019.[3] Laxmi Silwal-Pandit, Silje Nord, Hedda von derLippe Gythfeldt, Elen K Møller, Thomas Fleischer,Einar Rødland, Marit Krohn, Elin Borgen, ØysteinGarred, Tone Olsen, et al. The longitudinal tran-scriptional response to neoadjuvant chemother-apy with and without bevacizumab in breast can-cer.