Predicting the Response of Breast Cancer to Neoadjuvant Therapy Using a Mechanically Coupled Reaction–Diffusion Model

Neoadjuvant therapy is often administered to patients with locally advanced breast cancer in order to reduce the tumor burden prior to surgery (10). Patients exhibiting pathologic complete response (pCR; defined as no residual viable tumor on histologic analysis in breast or nodes at the completion of therapy) demonstrate significantly enhanced disease-free progression/survival; thus, this metric has emerged as the gold standard for evaluating response to neoadjuvant therapy (11, 12). With numerous neoadjuvant therapy regimens available, and many more treatments emerging, reliable prediction of patients that will go on to achieve pCR is highly significant. If accurate predictions that the regimen is not effective can be made early in a patient's treatment, the treatment could be changed to another, potentially more effective approach, thereby avoiding unnecessary treatment-related toxicities and providing a greater chance of achieving pCR (13).

Assessing early treatment response with the goal of predicting pCR remains a challenging clinical problem with no current uniform approach. Because of their quantitative and noninvasive nature, imaging-based biomarkers are attractive surrogate metrics for response assessment. Several imaging-based metrics have been investigated as biomarkers for predicting response in the neoadjuvant setting (14–22), with many efforts focusing on monitoring changes in morphologic characteristics. However, morphology-based characteristics are downstream from the underlying treatment-induced biologic response, and are, therefore, indicative only of late stage response. More recent efforts have focused on using functional parameters available from quantitative magnetic resonance imaging (MRI) techniques. For example, in the recent multisite American College of Radiology Imaging Network (ACRIN) 6657/Investigation of Serial Studies to Predict Your Therapeutic Response with Imaging And moLecular Analysis (I-SPY TRIAL) (22), Hylton and colleagues examined traditional clinical assessment combined with imaging morphologic measurements and contrast enhanced signal enhancement ratio (SER). A four-predictor multivariate statistical model (clinical size, longest diameter, volume, and SER) was shown to have receiver operator characteristic (ROC) area under the curve (AUC) of 0.73 for prediction of pCR at an early imaging time point, prior to the second cycle of neoadjuvant therapy. A preliminary extension of this approach to combine quantitative region of interest (ROI) and voxel-level analysis metrics from both dynamic contrast-enhanced (DCE) and diffusion-weighted (DW) MRI, efflux rate transfer constant (k_ep), and apparent diffusion coefficient (ADC), respectively, in a multivariate regression statistical analysis showed promise for achieving a higher predictive ability, with AUC for pCR prediction of 0.87 at an early imaging time point, after the first cycle of neoadjuvant therapy (17). These studies extract a single parameter or sets of parameters from quantitative imaging data analysis and perform a statistical correlation to clinical outcomes. However, the physiologic interpretation of statistical combinations of extracted parameters as they pertain to pCR or reduction in tumor burden is often largely unknown.

One promising way to maximize the impact of quantitative noninvasive imaging data is to use mathematical models of tumor growth driven by early time point imaging data in order to determine and predict therapeutic response. Major advances in mathematical modeling of tumor growth have occurred, where sophisticated simulations can be generated that recapitulate many cellular and bulk-level aspects of tumor growth and treatment response (23–25). Unfortunately, many such models are fundamentally limited in that they cannot readily incorporate clinical data and require knowledge of many different parameters of tumor growth/status that are either impractical or impossible to measure clinically (e.g., immune response, extracellular matrix status, genetic mutations, etc.). As complex, multiparametric models often require an extensive array of model parameters either obtained from literature values or empirically derived, they cannot be practically applied to individual patients. Thus, the field of mathematical oncology has been limited to more fundamental investigations and patient-specific predictive modeling has only seen limited clinical application. There is a clear need to develop clinically relevant oncological models that are amenable to available patient data and yet retain the most salient features of response prediction (13).

We and others (4, 26–31) have recently suggested using quantitative in vivo imaging data to drive mathematical models of tumor growth for predicting response. Following parameterization, these models can then be projected forward in space and time to predict eventual tumor distribution outcomes on a patient-specific basis. Rather than taking a classical reductionist approach where each appropriate physiologic constituent tumor-associated process is described mathematically and coupled into a comprehensive, yet increasingly complex representation, we have chosen to begin with a more traditional spatiotemporal growth model that captures first-order effects. In this approach, we use patient-specific quantitative imaging data in two ways: (i) to fit the model-associated parameters and (ii) to guide decisions whether to add complexity to the model based not purely on the existence of more sophisticated tumor-dynamic processes but rather how well incorporation can match patient outcome as determined at the time of surgery. We posit that systematically incorporating (additional) appropriate complexity to maximize clinical prediction, rather than constructing from the myriad of reduced fundamental interactions, has considerable merit for the ultimate realization of practical mathematical oncology in the clinical setting.

Recently, we have developed and shown preliminary results for a biomechanical modeling framework where tumor response is parameterized using data from both before and after the first cycle of neoadjuvant therapy and the model is driven forward in time to predict tumor burden at the time of surgery (9). The model is based on the standard reaction–diffusion equation of tumor growth/response, with terms to describe tumor cell diffusion and proliferation; in our formulation, the model is initialized and constrained using the serial imaging data. We extend this basic macroscopic tumor model to incorporate a mechano-inhibitory effect of surrounding tissue stroma on tumor cell diffusivity/motility. In this study, we build upon the previous efforts and apply patient-specific predictive modeling, constrained by quantitative imaging data, to a group of patients exhibiting a varying degree of response to therapy. We adopt a patient-scale spatiotemporal tumor growth model framework and allow readily available clinical imaging data to guide the determination of model parameters in a best-fit sense. The goal of this study is to determine if a patient-specific mechanically coupled reaction–diffusion model of tumor response to therapy can be used as a predictive indicator of pCR.

An existing database of patients with breast cancer receiving quantitative MRI before (t₁), after one cycle (t₂), and at the conclusion (t₃) of neoadjuvant therapy from a previous study (17) was utilized for this retrospective analysis. Briefly, 33 patients with stage II/III breast cancer who underwent neoadjuvant therapy as a component of their clinical care were selected for this study. Inclusion criteria were as follows: no prior systemic therapies for breast cancer and histologically documented invasive carcinoma of the breast with a sufficient risk of recurrence, based on pretreatment clinical parameters of size, grade, age, and nodal status, to warrant the use of neoadjuvant therapy. Neoadjuvant therapy was administered at the discretion of the treating medical oncologist based on pretreatment clinical parameters.

Tumor type was classified based on expression status for estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor 2 (HER2) and categorized as ER and/or PR positive, HER2 positive, or triple-negative (ER, PR, and HER2 negative); see Table 1. Most patients with HER2 positive tumors received paclitaxel, carboplatin, and trastuzumab every 3 weeks for six cycles. Most patients with ER and/or PR positive tumors received doxorubicin and cyclophosphamide every 2 weeks for four cycles followed by 12 weekly cycles of paclitaxel. Most patients with triple-negative tumors participated in an ongoing clinical trial and received weekly cisplatin and paclitaxel with or without everolimus for 12 weeks. Response status was determined at the time of surgery by a pathologist as either pCR, which is defined as no residual viable tumor on histologic analysis of the breast or nodes at the completion of neoadjuvant therapy, or non-pCR, which designates any residual invasive tumor on histologic analysis in breast or nodes at the completion of neoadjuvant therapy. Participating patients provided informed written consent to an Institutional Review Board approved study. For more information on the patient population, the reader is referred to Li and colleagues (17).

MRI was performed using a Philips 3T Achieva MR scanner (Philips Healthcare). A four-channel receive double-breast coil was used for 20 patients (Invivo), and a 16-channel double-breast coil was used for 13 patients. THRIVE (T₁ high-resolution isotropic volume examination) anatomical data were acquired via a 400 × 400 × 129 acquisition matrix over a 20 cm × 20 cm × 12.9 cm transverse field of view (FOV) with one signal acquisition, and TR/TE/α = 6.43 ms/3.4 ms/10°. Dynamic contrast enhanced MRI (DCE-MRI) utilized an acquisition matrix of 192 × 192 × 20 (full breast) over a sagittal square FOV (22 cm × 22 cm) with slice thickness of 5 mm, one signal acquisition, and a SENSE factor of 2. Dynamic scans used a flip angle of 20°. Each dynamic set was collected in 16 s at 25 time points. A catheter placed within an antecubital vein delivered 0.1 mmol/kg (9–15 mL, depending on patient weight) of gadopentetate dimeglumine, Gd-DTPA (Magnevist) at 2 mL/s, followed by a saline flush, via a power injector (Medrad) after the acquisition of three baseline dynamic scans. DW-MRI was acquired with a single-shot spin echo (SE) echo planar imaging sequence in three orthogonal diffusion encoding directions, with b-values of 0 and 500 or 600 s/mm², FOV of 19.2 cm × 19.2 cm, and an acquisition matrix of 96 × 96 reconstructed to 144 × 144. SENSE parallel imaging (acceleration factor = 2) and spectrally selective adiabatic inversion recovery fat saturation were implemented to reduce image artifacts. Subjects were breathing freely with no gating applied. The DW-MRIs consisted of 12 sagittal slices with slice thickness of 5 mm (no slice gap), TR = 3080 ms, TE = “shortest” (41 or 60 ms), Δ = 19.8 or 29 ms, and δ = 10.7 or 21 ms, NSA = 10. The total scan time for THRIVE, DCE-MRI, and DW-MRI data was approximately 2.7, 6.7, and 4.7 minutes, respectively.

A mechanically coupled reaction diffusion model of tumor response to therapy that used the MRI data obtained at time points t₁ and t₂ was used to predict tumor burden at the conclusion of neoadjuvant therapy. Although a brief description follows, the interested reader is referred to Weis and colleagues for details of the modeling methodology (9). Briefly, THRIVE, DCE-MRI, and DW-MRI data were longitudinally coregistered (32) across all time points. Central-slice images through the midpoint of the tumor were extracted from MRI volumes for further analysis. DCE-MRI data were used to define a tumor region-of-interest for voxels satisfying 80% signal enhancement threshold (9) between the precontrast images (i.e., the average signal intensity from DCE-MRI time points 1–3) and the postcontrast images (i.e., the average signal intensity from DCE-MRI time points 4–25). Voxels that showed a signal intensity increase more than the 80% threshold between the pre- and postcontrast data were defined as tumor (1). ADC maps were calculated from DW-MRI datasets using Eq. A and transformed to estimates of tumor cell number using Eq. B (3, 4). A two-dimensional, mechanically coupled reaction diffusion tumor model, Eqs. C–E, was used to describe tumor cell logistic growth/decay and diffusion based on the central slice of the MR data. This model incorporates a tumor growth “mass effect” where accumulated mechanical distortion energy has an inhibitory effect on tumor cell invasiveness (5). Tumor cellularity at baseline (t₁) and after one cycle of neoadjuvant therapy (t₂) was used to estimate the spatially varying proliferation rates and a global tumor cell diffusion parameter by a conjugate-gradient nonlinear optimization. Following parameter fitting, the model was projected forward in time (using the estimated proliferation and tumor cell diffusion parameters) to the final time point at the conclusion of neoadjuvant therapy (t₃), and the tumor cellularity was calculated as a model prediction. A schematic of the model parameter estimation and prediction framework is shown in Fig. 1. In addition, the cellularity maps for each patient at each time point for both the experimentally observed (t₁ and t₂) and model predicted (t₃) cellularity data were summed to yield an estimate of the total tumor cellularity. To eliminate volumetric bias, these numbers are expressed as a percent change relative to a reference time point; for example, a percent change in cellularity from t₁ to t₂ is defined as 100(t₂ − t₁)/t₁.

Statistical analysis between non-pCR and pCR patient groups for percent changes in each analysis metric was performed using the Wilcoxon test. Statistical significance was set at P < 0.05.

Receiver operating characteristic (ROC) curves were generated for prediction of pCR for both the mechanically coupled model and a standard reaction–diffusion model. ROC curves were analyzed and area under the curve (AUC) and 95% confidence intervals (CI) were generated for each model. The optimal cutoff was selected utilizing the Youden index, which maximizes the difference between the true positive rate and the false positive rate (33). Data were analyzed using GraphPad Prism (GraphPad Software).

Thirty-three patients (median age, 46; range, 28–67) completed scanning at t₁ and at t₂. The median time between t₁ and t₂ was 14 days (range, 7–29 days). Table 1 summarizes the salient features of the patient data set, including their tumor type and pathologic response. After neoadjuvant therapy, 12 of 33 patients (36%) were defined as having achieved pCR.

Imaging data and model predictions for tumors from two representative patients are shown in Figs. 2 and 3. Fig. 2 displays data from a patient that achieved pCR at the conclusion of neoadjuvant therapy that was scanned before, after the first cycle, and at the conclusion of neoadjuvant therapy in the form of anatomical images (Fig. 2A–C), ADC (from DW-MRI) values superimposed on the averaged postcontrast images (Fig. 2D–F), and the estimated cellularity (using Eq. B) as a fraction of the carrying capacity (Fig. 2G–I). The mechanically coupled reaction–diffusion model is used to reconstruct parameter estimates of global diffusion and spatial proliferation rate (Fig. 2J) based on observed cellularity data from the baseline and post-one cycle time points. The model is then evaluated forward in time based on the estimated parameters and observed cellularity data in order to predict cellularity at the conclusion of neoadjuvant therapy (Fig. 2K). Model prediction strongly agrees with data observations from imaging and pathology of complete response to therapy. Qualitatively, the estimated cellularity map is shown to exhibit a significant trend of decrease between t₁ and t₂. The prediction of the mechanically coupled reaction–diffusion model is seen to continue the trend of decreasing cellularity at t₃, with a prediction of near complete response. Figure 3 displays similar data for a representative patient with residual disease at the time of surgery (i.e., a non-pCR). The estimated cellularity map for this patient exhibits an opposite trend from that of the pCR patient, with increasing estimated cellularity observations from t₁ to t₂. The mechanically coupled reaction–diffusion model prediction is then shown to predict a lack of response to therapy with an increase in cellularity at t₃. Compared to observed cellularity data, the model overpredicts cellularity at the final time point. Combined, representative images and associated estimated parameter maps of proliferation rates depict overall prediction of regression of disease for a responding patient and significant progression of disease for a nonresponding patient.

Table 2 reports the overall performance for predicting pCR from a group of 33 patients exhibiting a varying degree of responsiveness to neoadjuvant therapy for analysis metrics based on RECIST (i.e., measurement of tumor longest dimension), observed cellularity (based on ADC values from DW-MRI), reaction–diffusion model prediction, and mechanics coupled reaction–diffusion model prediction. Area under the receiver operating characteristic curve, sensitivity, and specificity are reported for each metric, along with 95% confidence intervals. Figure 4 depicts the average value, 95% confidence interval as indicated by error bars, and P value for comparisons between non-pCR and pCR tumors for each analysis metric, along with ROC curves.

Simple analysis of differences in tumor morphologic and cellularity observations before and after the first cycle of neoadjuvant therapy exhibit a modest ability to predict pCR status at the conclusion of therapy. The change in tumor longest dimension between imaging time points t₁ and t₂, based on RECIST criteria, reflect an AUC and sensitivity of 0.71% and 92%, respectively; however, RECIST exhibits a poor specificity of 53% in this study. Similarly, the estimated tumor cellularity based on ADC values from DW-MRI is shown to have modest predictive performance, with ROC AUC and sensitivity of 0.73% and 100%, respectively, and poor specificity of 47%.

Using image-based mathematical modeling in order to predict tumor response to therapy is seen to improve pCR predictive performance. Reaction–diffusion model predictions of tumor cellularity at the conclusion of neoadjuvant therapy, expressed as a percent change from cellularity after one cycle of neoadjuvant therapy, exhibit an AUC, sensitivity, and specificity of 0.75%, 100%, and 55%, respectively. Here, model predictions are expressed as a percent change to eliminate volumetric bias from prediction results. Incorporating mechanical coupling into the reaction–diffusion model significantly improves pCR predictive performance, achieving an ROC AUC, sensitivity, and specificity of 0.87%, 92%, and 84%, respectively.

We previously introduced a novel mathematical framework for predictive modeling that is centered on integrating quantitative in vivo imaging data with biomechanical models of tumor growth (9). This approach represents a significant departure from existing paradigms in mathematical oncology modeling by using patient-specific data to parameterize a model that can then be used to make a prediction at a future time point that can then be checked directly against patient response. In this work, we retrospectively apply this framework to a database of patients exhibiting varying degrees of responsiveness to neoadjuvant therapy. We compared model predictions of responsiveness to pathologic response assessment at the time of surgery. Our results demonstrate a significant predictive ability of our mechanics coupled modeling approach in identifying patients that would go on to achieve pCR, as indicated by a strong AUC of 0.87. Our mechanics coupled reaction–diffusion modeling approach significantly outperformed other, commonly used metrics, including analysis of tumor longest dimension (i.e., RECIST), early-cycle changes in DW-MRI (used as empirical input data to the models), and standard reaction–diffusion predictive modeling. Compared with the standard reaction–diffusion model, the mechanically coupled model exhibits significantly improved prediction specificity, reflected by improved prediction of true non-pCR patients (55% vs. 84%, respectively). Further analysis of the patients for which mechanically coupled modeling improved predictability over standard modeling (data not shown), failed to reveal correlations with any other known factor (including breast density, patient age, and breast cancer subtype, among others).

Our results also confirm previous reports of the use of early cycle DW-MRI ADC data alone in predicting neoadjuvant therapy response (34–39). Trends of decreasing/increasing cellularity are clearly visible between the observed cellularity before and after the first cycle of therapy, as shown in Figs. 2 and 3, with predictive AUC of 0.73. However, analysis of these data is seen to exhibit limited specificity, as has been previously demonstrated (35, 40).

The results from this study are promising, but we recognize that there are several limitations in our current approach. First, we estimate the key model parameters (i.e., proliferation and tumor cell diffusion) from data obtained before and after the first cycle of neoadjuvant therapy; this is an admittedly (temporally) sparse set of measurements for attempting to capture the dynamics through a model fitting procedure. In addition, we have assumed that reconstructed model parameters do not vary in time; this is essentially equivalent to assuming a constant antitumor effect during neoadjuvant therapy treatment and this is clearly a simplification of the response of tumor cells to neoadjuvant therapy. True response undoubtedly exhibits strong nonlinear temporal response to therapy; indeed, the evolution of a resistant phenotype in the presence of continued chemotherapeutic administration speaks directly to a nonlinear effect that is not currently taken into account in our approach. However, it is clinically infeasible to acquire sufficient imaging time points throughout the course of therapy to accurately predict the temporal changes in proliferation rates and tumor cell diffusion with the current model. We are currently exploring acquisition of one additional early imaging time point, collected after two cycles of therapy, in order to better understand the temporal evolution of parameters and more accurate prediction of response.

A second limitation, as also noted in Li and colleagues (17), is the variation in the breast cancer subtype, therapeutic regimen (including the specific therapeutic agent and dosing frequency), and number of days elapsed between imaging time points, t₁ and t₂. Thus, it is possible that the modeling approach could be influenced by the underlying dynamics and biology of specific breast cancer subtypes as well as the particular therapy regimen utilized. This source of uncertainty due to patient-specific biologic variability is particularly difficult to control, due to the nature of clinical research data acquisition constraints. It will be important, as a greater number of patients are analyzed with this methodology, to statistically analyze the model relative to cancer subtype, therapeutic regimen, and imaging time point delay to determine these effects. In addition, future work toward model validation in preclinical and in vitro systems will facilitate greater experimental control and provide further understanding toward model robustness in the presence of such biological variability. Of course, the encouraging preliminary results achieved with our modeling approach in spite of this variability points to the importance, and potential generalizability, of using models based on patient-specific information.

Another limitation is that our current approach utilizes central-slice MR data to initialize and constrain two-dimensional mathematical models, rather than using full image volumes and three-dimensional models. We do note, however, that we are utilizing somewhat thick slices, 5 mm, due to the spatio-temporal acquisition requirements of the trial for which the MRI data were originally acquired. Extension of our approach to three-dimensional analysis is straightforward, and we would expect an enhanced predictive capability through the use of additional volumetric data. However, the current computational demands for iterative reconstruction of volumetric parameter maps significantly limit the feasibility of extending our approach to three-dimensional analysis for an entire patient database at this time. Currently, efforts are underway to refine our approach by exploring significant computational speed enhancements though the use of high-performance computing cluster and/or graphics processing unit processing in addition to algorithmic enhancements.

The mechanically coupled reaction diffusion model adds two additional degrees of freedom to each node in the finite element mesh (in addition to the existing degree of freedom for tumor cell number) due to the need for calculation of a displacement solution. This added complexity increases the computation time for the forward model by approximately 40%, without parallelization of matrix assembly and/or matrix inversion. In addition, the impact on total wall-clock time is relatively limited as the entire mechanically coupled reaction diffusion model prediction approach typically requires only 1 to 2 hours for each patient.

At present, our mechanics coupled approach relies on the use of literature values for the elastic properties of the breast (9). In addition, our approach currently ignores other important biologic factors beyond simple growth and mechanoinhibitory diffusion. Other ongoing efforts are designed to address these limitations through incorporation of additional quantitative patient-specific imaging methodologies to inform more complex models (13).

Given these limitations, we note that our approach compares very favorably, with similar AUC, to previous multimodal statistical regression analysis of this same patient population in previous work that combines ROI and voxel-level DW-MRI and DCE-MRI analysis (17). Importantly, our approach currently includes data from only one quantitative imaging source, DW-MRI. Contrast-enhanced MRI data are used in this work only to identify tumor ROI, with no sensitivity to pharmacokinetic parameters fit from dynamic DCE-MRI data. It is therefore very encouraging that our model-based analysis framework, based on unimodal and two-dimensional central-slice input data, achieves similar predictive results to multimodal full-volume analysis.

In summary, we have demonstrated a model-based analysis method that provides predictions of therapeutic response using patient-specific early time point imaging data within the course of therapy. The framework provides good agreement with pathologic response observations as demonstrated by a strong predictive ability. The results suggest that an imaging-driven biophysically constrained modeling approach that combines reaction–diffusion with mechanical effects is a predictive indicator of response to therapy. Future work will seek to further enhance the predictive capacity of our approach by expanding the mathematical model of tumor response to include model terms that can be populated with other sources of quantitative patient-specific noninvasive imaging, including DCE-MRI, fludeoxyglucose positron emission tomography (FDG-PET), and magnetic resonance elastography (MRE).

No potential conflicts of interest were disclosed.

Conception and design: J.A. Weis, M.I. Miga, T.E. Yankeelov

Development of methodology: J.A. Weis, M.I. Miga, T.E. Yankeelov

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.A. Weis, L.R. Arlinghaus, V. Abramson, A.B. Chakravarthy, T.E. Yankeelov

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.A. Weis, X. Li, T.E. Yankeelov

Writing, review, and/or revision of the manuscript: J.A. Weis, M.I. Miga, L.R. Arlinghaus, X. Li, V. Abramson, A.B. Chakravarthy, T.E. Yankeelov

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.A. Weis, A.B. Chakravarthy, P. Pendyala

Study supervision: J.A. Weis, A.B. Chakravarthy, T.E. Yankeelov

The authors offer their sincerest appreciation to the women who volunteered to participate in this study. The authors thank the National Institutes of Health for funding through NCI 1U01CA142565, NCI U01CA174706, NCI R25CA092043, NCI 1P50 098131, NCI P30CA68485, NCI R01CA138599, and NINDS 2R01NS049251. The authors also thank the Kleberg Foundation for generous support of the imaging program at Vanderbilt and the support provided by the Vanderbilt Initiative in Surgery and Engineering Pilot Award Program.

This work was supported by the National Institutes of Health: NCI 1U01CA142565, NCI U01CA174706, NCI R25CA092043, NCI 1P50 098131, NCI P30CA68485, NCI R01CA138599, and NINDS 2R01NS049251

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.