Beyond RECIST: mathematical modeling and Bayesian inference reveal the importance of immune parameters in metastatic breast cancer

Abstract

Immunotherapies that target the host immune system to mount effective responses hold great promise. Yet, overcoming patient- and organ-specific tumor heterogeneities remains a significant challenge. In order to quantify individual patient responses, we fit a tumor-immune mathematical model to patient and site-specific dynamics during combination therapy (nivolumab + ipilimumab + entinostat) informed by RECIST measurements of the tumor dynamics and immune markers measured by spatial proteomics. Bayesian parameter inference of site-specific patient responses revealed that only the immunosuppression parameters were predictive of response; parameters controlling cytotoxicity were uninformative. Via comparison of a large cohort of fitted tumors, we quantified the variability in tumor-immune dynamics to reveal controllable parameter regimes. We developed methods that employed posterior parameter sampling and simulation to create virtual tumor populations, enabling extrapolation beyond the data to predict probabilities of response in metastatic lesions, even when no data exist at a site. We also showed that scans in the week immediately following treatment are particularly valuable to identify the tumor dynamics. Our modeling and inference framework can thus be used to overcome sample size limitations to create virtual patient cohorts that give new insights into mechanisms of disease progression.

Figure 1:. Target lesion tumor measurements and IMC data from tumor tissue shows response or non-response to combination treatment.

A. Relative change in tumor volume for 55 individual tumors (from 20 patients) plotted by patient RECIST response classification. 5 patients were classified as responders (partial/complete response) and 15 patients were classified as non-responders (progressive/stable disease). B. Relative change in tumor volume for 55 individual tumors plotted by site of tumor. See the SI for further details. C. Proportions of high-level immune cell populations identified via IMC in tumor tissues from 10 patients at baseline, C1D1 (post entinostat), and week eight (post combination therapy).

Figure 2:. Bayesian tumor fitting to mathematical model provides algorithmic pipeline to infer direct tumor-specific immune response to treatment.

A. Diagram of analysis pipeline. B. Posterior parameter distributions and model fits for a tumor responsive to combination treatment. C. Posterior parameter distributions and model fits for a tumor not responsive to combination treatment.

Figure 3:. MDSC inhibition of anti-tumor populations and CTL stimulation control tumor growth.

A. Parameter table and notation. B. Posterior distributions (green) and prior distributions (gray, dashed) are shown for five tumors. ${\beta }_3$ and ${\beta }_4$ decrease as tumor response improves, whereas ${\alpha }_6$ increases for responsive tumors. ${\beta }_1$ and ${\beta }_2$ show no change (posterior distribution compared to prior distribution) for both responsive and non-responsive tumors.

Figure 4:. Site-specific patient-specific tumor prediction enables the analysis of a large virtual patient population undergoing combination treatment.

A. Individual tumor data from liver (purple) and patient R-18 (light green). B-C. Posterior parameter distributions from liver and patient R-18. D. 10² tumor simulations for liver tumors, patient R-18 tumors, and liver and patient R-18 tumors (blue). The fraction of tumors that decrease (green), are stable (orange), and increase (red) are also shown.

Figure 5:. Immune population changes not captured by RECIST control tumor response to treatment.

A. 10² tumor simulations for patient R-18, patient R-2, and patient R-6. The fraction of tumors that decrease (green), are stable (orange), and increase (red) are also shown. B. Same as A but tumor simulations include an ${\alpha }_6$ (CTL stimulation by tumor-NK cell interaction) modulation as quantified by IMC data.

Figure 6:. Parameter identifiability is constrained by the temporal density of tumor data.

Left panels show tumor data (red dots) along with the mathematical model simulation with ${\beta }_3$, ${\beta }_4$, ${\alpha }_6$ at MLE values. Red lines are tumor cells, yellow lines are MDSCs, green lines are NK cells, and blue lines are CTL cells. Right panels show univariate profile likelihood functions for the three parameters. The likelihood function (blue line) is superimposed with a vertical green line at the MLE and a horizontal red line at the asymptotic 95% threshold. A. Actual tumor data from clinical trial (tumor 1). B. Simulated tumor data every day. All three parameters are identifiable and well-constrained. C. Simulated tumor data at 11 equally spaced points (every 30 days) over the timespan. D. Simulated tumor data every day for the first week after treatment initiation and then three equally spaced points after that, same number of data points as panel C.