MRI-Based Digital Models Forecast Patient-Specific Treatment Responses to Neoadjuvant Chemotherapy in Triple-Negative Breast Cancer

Neoadjuvant systemic therapy (NAST) is widely considered the standard-of-care for treatment of stage II and III, locally advanced triple-negative breast cancer (TNBC). NAST increases the success rate for breast conservation surgery by reducing tumor burden and provides the opportunity to treat micrometastases in a naïve state, thereby improving progression-free survival of patients (1, 2). Importantly, patients with TNBC who achieve a pathologic complete response (pCR) in the neoadjuvant setting have a favorable recurrence-free survival; in contrast, patients who have residual disease after NAST are at increased risk of early recurrence and death (3, 4). Unfortunately, based on recent pooled analysis of 52 studies from 1999 to 2016, only 32.6% of TNBC patients treated with standard taxane/anthracycline-based NAST have a pCR or minimal residual disease at the time of surgical resection (4).

Development of novel neoadjuvant treatment regimens has provided opportunities to tailor treatment for individual patients to improve outcomes for TNBC (5–7). Thus, it is becoming increasingly important to develop techniques that can accurately predict the individual TNBC patient's response to NASTs. If it could be definitively determined that a therapeutic regimen is unlikely to achieve pCR for a patient, then risk-adapted therapy could be adopted, with addition of rationally based treatments or removal of unnecessary components, potentially improving outcomes and lowering side effects. In particular, the importance of being able to remove patients from ineffective therapies, as early as possible, is difficult to overstate given their significant toxicities (8), including increased likelihood of hospitalization, cardiac damage, leukemia, and even death (9). Moreover, accurately predicting response to NAST would enable identification of exceptional responders who might benefit from treatment de-escalation (10), including the possibility of nonsurgical management of their disease (11).

Numerous efforts have been devoted to investigating approaches that can accurately distinguish pCR from non-pCR patients early in NAST. Imaging biomarkers derived from MRI, PET, and ultrasound imaging have been shown to be strongly correlated with the response of breast tumors to NAST (12–14). In particular, MRI measurements before and during NAST have been valuable predictors of pCR, especially the functional tumor volume (FTV) derived from dynamic contrast-enhanced (DCE) MRI and the apparent diffusion coefficient (ADC) derived from diffusion-weighted (DW) MRI (15–18). More recently, the methods of artificial intelligence have been used to extract features from high-dimensional data to build predictive models for differentiating pCR from non-pCR in breast cancer (19, 20). However, the majority of these approaches for predicting or assessing response have the intrinsic limitation of being population-based. Importantly, population-based approaches, which rely exclusively on statistical inference from properties of large populations, inevitably obscure conditions specific to the individual patient over time (21), especially for a disease as heterogeneous as cancer (22). Conversely, biologically based models employing patient-specific data have the potential to shift the paradigm from population- to individual-based approaches (7). Furthermore, biologically based mathematical modeling of tumor response can not only predict the changes in global metrics summarizing tumor burden (e.g., total tumor volume and cellularity), but also reveal biologically specific information (e.g., spatially resolved maps of proliferation, pharmacokinetics, and each patient's sensitivity to the administered therapies). It promises unique opportunities to characterize tumor pathophysiology, to rigorously forecast long-term outcomes, and even optimize treatment plans on a patient-specific basis (23, 24).

In this contribution, we develop a clinical-computational approach to establish patient-specific models to make early predictions on the spatiotemporal development and response of individual patients with TNBC to standard-of-care NAST. As this model can represent a physical object (i.e., tumor), predict the behavior of the object given influences (i.e., treatments), and enable decision-making to optimize the future behavior of the object (i.e., improve the treatment outcomes; refs. 25, 26), we posit that our methodology represents a practical manifestation of digital twins in oncology. The approach requires no population-based training of models; instead, it integrates an individual patient's multiparametric MRI data obtained at multiple time points during their treatment (Fig. 1A) with biologically based mathematical modeling. In particular, two frameworks were constructed to determine the predictive utility of each patient's digital twin (Fig. 1B). In Framework 1, we seek to employ digital twins to predict the outcome of a single NAST regimen (i.e., A/C); specifically, we evaluate the accuracy of digital twins to predict global metrics related to change in tumor burden, as well as spatiotemporally resolved tumor dynamics at the end of A/C (Fig. 1C). In Framework 2, we seek to employ digital twins to predict the outcome of the entire NAST (i.e., both A/C and paclitaxel); specifically, we evaluate the accuracy of digital twins to predict individual pCR or non-pCR status at the completion of NAST (Fig. 1D).

Treatment-naïve patients with biopsy-confirmed TNBC were enrolled in the prospective, Institutional Review Board–approved clinical trial, “A Robust TNBC Evaluation FraMework to Improve Survival” (ARTEMIS, NCT02276443). Among patients who signed informed consent between June 2018 and January 2020, those who had clinical stage I to III disease, had completed NAST with a complete series of MRI scans, and had known postsurgical pathology were included in this study (n = 56).

Each patient underwent multiparametric MRI acquisitions before treatment (baseline/V1), after two cycles (V2), and after four cycles (V3) of the standard-of-care (neoadjuvant) doxorubicin and cyclophosphamide (A/C). (Each “cycle” of A/C is 2 weeks; see Fig. 1A.) Patients with disease progression, or <70% reduction in tumor volume at the end of A/C (27), were offered the opportunity to enroll in a biomarker-guided clinical trial using targeted bio/chemotherapy to complete therapy (n = 9). Patients not meeting the criteria for suboptimal response to A/C were recommended to continue standard-of-care paclitaxel weekly for 12 cycles (n = 43), or double-dosed every 3 weeks for four cycles (n = 4). (The exact paclitaxel regimen for each patient's was determined by their physician.) All patients received surgery after NAST. The postsurgical pathology was used to classify patients as pCR or non-pCR; pCR was defined as the absence of residual invasive and in situ cancer on hematoxylin and eosin evaluation of the complete resected breast specimen and all sampled regional lymph nodes following completion of NAST (28).

MRI was performed on a GE Discovery MR750 or MR750w whole-body scanner (GE Healthcare) with an eight-channel bilateral breast coil. In particular, the DCE-MRI data were acquired using a 3D DISCO sequence with bipolar readouts (Fig. 2A; ref. 29) and the following scan parameters: field-of-view = 30 × 30 cm², matrix size = 320 × 320, slice thickness/spacing = 3.2/−1.6 mm, number of slices = 140, flip angle = 12°, repetition/echo time = 8/2 milliseconds. After one precontrast phase was obtained, a single bolus of contrast agent (Gadovist, Bayer HealthCare) was injected (0.1 mL/kg at ∼ 2 mL/second followed by saline flush) at the start of the postcontrast acquisition. The temporal resolution of the DISCO series ranged from 8 to 15.5 seconds (median = 11 seconds), depending on the slice coverage, which resulted in the number of postcontrast phases varying from 32 to 64 (18). The DW-MRI was acquired using a 2D spin-echo sequence of the diseased breast with the following scan parameters: field-of-view = 16 × 16 cm², matrix size = 80 × 80, slice thickness/spacing = 4/0 mm, number of slices = 16, flip angle = 90°, repetition/echo time = 4,000/70 milliseconds. The b values used were 100 and 800 second/mm². The ADC map was calculated using a GE AW server (v3.2; GE Healthcare). The tumors were manually segmented by two board-certified breast radiologists with 4 to 12 years of experience (authors MB, RMM). Tumor segmentations were reviewed by two breast fellowship trained radiologists with 19 and 20 years of experience (authors GMR, BEA). The entire tumor volume was segmented on two early phases (60 and 150 seconds) of the DCE-MRI using an in-home software package developed in MATLAB (R2021b; Mathworks). All segmentations were further refined using the thresholding tool of the package to exclude nontumor voxels that were determined by radiologists. Necrotic regions and artifacts from the biopsy clip were manually segmented and excluded by two radiologists (authors MB, RMM).

All MRI data from each patient were processed through a pipeline that consists of three components: (i) preprocessing, (ii) inter-visit registration, and (iii) post-processing (Fig. 2B–D). This highly automated pipeline allows efficient processing of multi-visit, multi-parametric MRI with minimal user input.

First, the multiparametric images are colocalized to the same imaging grid to align slices and voxel locations. Specifically, DCE-MRI collected bilaterally were trimmed to the DW-MRI field-of-view covering the diseased breast, and the slices of DW-MRI and ADC maps were linearly interpolated to match the slice locations of DCE-MRI. A rigid registration was applied on the DCE-MRI to align all phases in one scan to the precontrast phase (MATLAB function, imregtform). A rigid registration was applied between the colocalized DCE-MRI and DW-MRI data to remove the small mismatches between the image volumes.

Second, an inter-visit image registration was performed to account for the change of breast tissue shape and patient position across MRI visits. Specifically, the registration was performed to align the images from V1 and V3 to the images from V2. The algorithm consisted of a rigid registration of the tumor ROIs for initial alignment, followed by a deformable b-spline, nonrigid registration on the whole breast with a rigidity penalty on the tumor regions (30–32). This rigidity penalty was imposed to preserve the tumor volume and shape across all visits. This registration was developed based on MATLAB and an open-source, command-line software, Elastix (33).

Third, the post-processing was performed as preparation for the subsequent predictive modeling. Specifically, a semiautomatic segmentation of the breast contour was performed on the precontrast frame of DCE-MRI based on a manually chosen intensity threshold, followed by a smoothing of the segmented mask edge (MATLAB function, imgaussfilt). A two-class k-means clustering (MATLAB function, kmeans) was used to segment fibroglandular and adipose tissues in each precontrast DCE-MRI (34). The enhancement of each DCE-MRI was calculated by subtracting the precontrast phase from the average of the postcontrast phases. A tumor cellularity map, N(x, t), was estimated based on the measured ADC map of each MRI visit (35):

where ADC_w is the ADC of free water (3 × 10⁻³ mm²/s), ADC(x, t) is the ADC value for the voxel at position x and time t, and ADC_min is the minimal (positive) ADC value in the tumor for the patient across all visits. The carrying capacity, θ, describes the maximum number of tumor cells that can physically fit within a voxel, which was determined by assuming a spherical packing density of 0.7405 and a nominal tumor cell radius of 10 μmol/L (30).

We have previously developed a biologically based mathematical model to represent the spatiotemporally resolved dynamics of tumor growth and response to NAST (30, 36). Specifically, a reaction-diffusion partial differential equation is used to describe the evolution of tumor cells, N(x, t), in response to the therapies, as shown in Eq. B:

where a detailed list of variables, parameters, their definitions, and assignments are given in Supplementary Table S1. The first term on the right-hand side of Eq. B describes the movement of tumor cells, as well as the compression of surrounding tissue, by a diffusion process coupled to tissue mechanical properties (37) via Eq. C:

where D₀ is the tumor cell diffusion coefficient in the absence of external forces. The exponential term reduces tumor cells’ mobility due to the surrounding tissue stiffness via the von Mises stress, σ(x, t), and an empirical coupling constant, γ. The von Mises stress was calculated for the fibroglandular and adipose tissues within the breast, with the fibroglandular tissue assigned a greater stiffness than the adipose tissue. Technical details on the mechanical-coupled diffusion can be found in ref. 37. The proliferation of the tumor cells (second term on right-hand side of Eq. B) is described by logistic growth with a spatially varying proliferating rate, k(x), and a global carrying capacity, θ. The effects of administered therapies (third term on the right-hand side of Eq. B) is modeled as the treatment-induced death rates of tumor cells, λ_i(x, t). The death rates are determined by the concentration of the administered drugs and the exponential decay of their efficacies:

where α_i is the efficacy of the ith administered drug, and i = 1, 2, and 3 refers to doxorubicin, cyclophosphamide, and paclitaxel, respectively. As each drug was administered multiple times during the therapeutic regimen, the total effectiveness of drugs at time t was an accumulation of all administered cycles and their decays; that is, the jth administration of the ith drug was conducted at time τ_i,_j at a total of J_i times. The decay of drug efficacy from each administration is represented by the decay rate of β_i. The spatial distribution of ith drugs caused by the jth administration of this therapy, C(x, τ_i,j), is determined by the enhancement of DCE-MRI (30). Specifically, from the last DCE-MRI data collected before the injection at τ_i,j, the area under the DCE-MRI time course is calculated at each voxel and normalized by the maximum value within the tumor. The voxel-wise normalized AUC represents the map of drug concentration induced by this injection (see Supplemental Section 1.1 for details).

Equations B to D are personalized by the imaging and clinical data of each patient. Specifically, the geometry of the computational domain is determined by the segmented breast contour, tumor, and fibroglandular/adipose tissues. The tumor cellularity map at each imaging time point is determined from the ADC map obtained at the corresponding time point via Eq. A. The sequential cellularity maps are then used for patient-specific calibration of the model parameters [i.e., k(x), D₀, α_i, and β_i]. In addition, the DCE-MRI acquired during NAST were used for updating the spatial distribution of the administered drugs and the mechanical properties of breast tissues. The model constrained by patient-specific data was implemented in MATLAB and solved with the finite difference method. Details of the numerical implementation can be found in refs. 30 and 36.

As shown in Fig. 1C, Framework 1 focuses on employing digital twins to predict the outcome of a single NAST regimen (i.e., A/C). Specifically, for each patient, the processed images from V1 and V2 are imported with the treatment regimen and calibrated to Eqs. B–D. The calibrated model is the digital twin as it represents patient-specific pathophysiological properties of tumor growth and response, including pretreatment tumor shape and cellularity, the proliferation rate and mobility of tumor cells, and the efficacy and decay of administered drugs (A/C). As the efficacy and decay rates were strongly coupled and challenging to simultaneously calibrate with only two-time-points, the decay rates of A/C were randomly sampled five times from ranges found in the published literature (38): β₁ ∈ [0.01, 0.6] day^–1, β₂ ∈ [1.0, 5.4] day^–1. With each set of sampled β₁ and β₂, the efficacy of A/C was calibrated, resulting in five parameter sets (see Supplementary Section S2.1). By applying the remainder of the A/C (i.e., all A/C after V2; Fig. 1C) to the digital twin, the patient-specific spatiotemporally resolved tumor cell distributions are predicted up to the end of A/C; one prediction is given by each parameter set, resulting in a median and a range of predictions. Furthermore, for each patient, global metrics—total tumor cellularity (TTC) and total tumor volume (TTV)—are derived from the spatiotemporally resolved predictions.

The predictive accuracy of Framework 1 is evaluated both temporally and spatially. The temporal accuracy is assessed by the agreement between the predicted and measured global metrics (i.e., TTC and TTV). In particular, the concordance correlation coefficient (CCC; see Supplementary Section S1.2) is computed between the predicted and measured changes of TTC at the end of A/C; similarly, the CCC is computed between the predicted and measured changes of TTV. The spatial accuracy is assessed by the difference between the predicted and measured spatially resolved tumor cell distributions. In particular, for each patient, the percent change of tumor cell counts from baseline (V1) to the end of A/C (V3) can be calculated at each location, x. We calculate the change from predicted and measured tumor cell distributions [i.e., ∆TCD_p(x) and ∆TCD_m(x), respectively], then the spatially resolved difference is given by ∆TCD_p(x) – ∆TCD_m(x). In each patient, this difference is reported as the median and interquartile range within the original tumor region; across the cohort, the differences are evaluated by the mean and 95% confident interval (CI) of the median difference from each patient. (Considering five predictions are given for each patient, all the evaluations above are based on the median of the five predictions.)

As show in Fig. 1D, Framework 2 focuses on employing digital twins to predict the outcome of the entire NAST (i.e., both A/C and paclitaxel). Specifically, for each patient, the processed images from V1, V2, and V3 are imported into the mechanism-based model for initialization and calibration. The calibrated model is the digital twin. (Note that the three-time-point data enables calibration of the efficacy and decay rates simultaneously within the same ranges assumed in Framework 1, so the sampling scheme in Framework 1 is not needed for Framework 2.) As no imaging data were available during the paclitaxel regimen, the efficacy of paclitaxel was set to literature value, α₃ = 0.3 day^–1 (38), and the decay rate of paclitaxel was assumed as the average of calibrated A/C decay rate; β₃ = (β₁ + β₂)/2. (See Supplementary Section S2.1 for details.) By applying paclitaxel to the digital twin, it predicts the patient-specific spatiotemporally resolved tumor cell distributions, TTC, and TTV at the end of NAST.

The output of Framework 2 is evaluated by the ability of the digital twins to differentiate pCR and non-pCR. Specifically, ROC analysis is performed on the predicted TTC and TTV. We report the area under ROC curve (AUC), sensitivity (i.e., the ability to correctly identify residual tumor at final surgical pathology), and specificity (i.e., the ability to correctly identify pCR at final surgical pathology) based on the optimal cut-off. In addition, the AUC/sensitivity/specificity from the predicted TTC and TTV were compared with those obtained by the measured TTC and TTV. The 95% CIs of the AUCs were computed and compared via DeLong method (39), with P < 0.05 considered statistically significant.

Raw data for this study were generated at the University of Texas MD Anderson Cancer Center. Raw data are not publicly available due to IRB restrictions of data containing information that could compromise research participant privacy and/or consent. The derived data that support the findings of this study are available from the corresponding author upon reasonable request.

A cohort of 50 patients was used in Framework 1. Six patients were excluded from the entire patient cohort (n = 56) due to image acquisition error or artifacts (n = 1), unsuccessful inter-visit registration caused by a large change of breast shape between visits (n = 1), and a complete response of the tumor at V2, which led to no data for model calibration (n = 4).

Each patient's digital twin provided a range of estimated treatment efficacies of A/C (e.g., Fig. 3A and B), simulating a range of treatment outcomes (e.g., Fig. 3C and D). Specifically, Fig. 3C presents a patient who had a suboptimal response to A/C (i.e., V3 imaging showed a <70% reduction in tumor volume). The digital twin predicted a mean (range) TTC and TTV of 3.17 × 10⁸ (3.05 × 10⁸ – 3.29 × 10⁸) cells and 3.19 × 10³ (3.15 × 10³ – 3.23 × 10³) mm³ at V3, respectively. This corresponds to predicted decreases of 49.55% (47.60%–51.50%) and 36.14% (35.44%–36.84%) for TTC and TTV at V3, respectively. In comparison, the measured decreases are 38.99% and 32.58% for TTC and TTV at V3, respectively. In contrast, Fig. 3D presents a patient who had a good response to A/C (i.e., V3 imaging showed >70% reduction in tumor volume). The digital twin predicted a TTC and TTV of 3.82 × 10⁶ (1.69 × 10⁶ – 5.94 × 10⁶) cells and 12.63 (0.00–25.27) mm³ at V3, respectively. This corresponds to predicted decreases of 98.00% (96.88%–99.11%) and 99.24% (98.48%–100.00%) for TTC and TTV at V3, respectively. In comparison, the measured decrease is 100% for both TTC and TTV at V3. Across the cohort, the CCC between the predicted and measured changes in TTC at V3 was 0.95 (Fig. 3E); the CCC between predicted and measured changes in TTV was 0.94 (Fig. 3F). These results indicate a high predictive accuracy and precision (i.e., uncertainty in the model's prediction; see Supplementary Section S2.2 for interpretation) of the temporal dynamics of the response of TNBC to neoadjuvant A/C (Table 1).

Importantly, the personalized digital twins provide not only the global metrics summarized in the previous paragraph, but also the spatiotemporal evolvement of each patient's tumor. Figure 4A and B show the measured and predicted tumor cell distributions, respectively, from the central slice of the same two example patients; and Fig. 4C and D present the 3D rendering of measured and predicted tumor volumes, respectively. In both cases, the digital twins successfully capture the lack (Fig. 4A) or presence (Fig. 4B) of response. Quantitatively, for the first patient, the difference between the predicted and measured change of tumor cell distributions at V3 has a median (interquartile range) of −3.30% (−22.07%–0.00%); for the second patient, the difference has a median (interquartile range) of 0.00% (0.00%–0.00%). The median differences across all patients had a mean (95% CI) of 0.20% (−20.35%–20.75%) at V3 (Fig. 4E). These results indicate high predictive accuracy and precision of spatially resolved predictions of tumor cell distributions in patients with TNBC in response to neoadjuvant A/C (Table 1).

A cohort of 37 patients (18 pCR, 19 non-pCR) was used in Framework 2. After excluding the 6 patients from the entire cohort (n = 56) as conducted in Framework 1, 13 more patients were excluded due to presumed nonresponse to further standard-of-care chemotherapy and enrollment into clinical trials (n = 9; part of ARTEMIS schema), missing the schedule of paclitaxel (n = 2), and errors in image acquisition (n = 2; ADC maps did not cover the entire tumor at V3).

For each patient, the personalized digital twin estimated treatment efficacies based on the V1 to V3 images (Fig. 5A), and represented tumor dynamics in response to A/C. These results were then used to predict the response to paclitaxel and, therefore, final treatment outcome after all NAST (Fig. 5C). For example, Fig. 5C shows the digital twin was able to be calibrated during the A/C regimen and used to predict no regrowth during paclitaxel for a patient that did, in fact, achieve pCR after NAST. In contrast, Fig. 5D shows the digital twin captured an initial response and then subsequent regrowth during A/C, and predicted strong regrowth before and during paclitaxel, resulting in a predicted TTC and TTV values of 9.03 × 10⁸ and TTV of 5.21 × 10³ mm³, respectively, after NAST.

Across the cohort, as shown in Table 2, the predicted TTC and TTV (based on the model calibrated with the V1–V3 data) from the digital twins both achieved an AUC of 0.89 (0.78–0.99) for distinguishing pCR from non-pCR (with postsurgical pathology as the ground truth). In comparison, the measured TTV or TTC (based on the V3 data) achieved an AUC of 0.78 (0.62–0.94) for distinguishing pCR from non-pCR. Furthermore, using the predicted TTC and TTV, the specificity was 0.95 and 0.89, respectively. In comparison, if using the measured TTV or TTC, the specificity was only 0.79. These results demonstrate that compared with the directly measured data, the digital twins improved the prediction of final response. Specifically, the AUC was improved by 14.28% for TTC (P = 0.04; significant) and 13.83% for TTV (P = 0.07). The specificity was improved by 20.25% and 12.66% for TTC and TTV, respectively; the sensitivity was unchanged.

We have developed a digital twin approach to achieve early, patient-specific, spatiotemporally resolved predictions of the response of patients with TNBC to neoadjuvant doxorubicin, cyclophosphamide, and paclitaxel. This approach was based on a biologically based mathematical model calibrated with multi-visit, multiparametric MRI acquired for the individual patient. Thus, the methodology represents a significant step away from population-based predictions, and towards individual-based predictions.

Framework 1 shows how the digital twin employs the imaging data from individual patients acquired early in the course of neoadjuvant A/C to make very accurate predictions of tumor status at the conclusion of A/C. The CCC between the predicted and measured values of total tumor burden and total tumor volume were 0.95 and 0.94, respectively. This strongly indicates that early changes during the A/C regimen contain sufficient information to calibrate the digital twin and confidently predict tumor response at the end of A/C. This observation aligns with previous reports demonstrating that metrics from early treatment MRI are strong predictors of NAST response in breast cancer (15–18, 40). This provides strong support that our approach could be used to adjust treatment regimens on a patient-specific basis. For example, our approach could be applied after the first two cycles of A/C to predict if further dosing with A/C should be continued, or if an alternative intervention should be considered.

Framework 2 shows how the digital twin employs the imaging data from individual patients acquired during the A/C portion of NAST to make accurate predictions of their final pathologic status (i.e., pCR or non-pCR) at the completion of all NAST, with an AUC = 0.89. Importantly, using the tumor volume measured at the end of A/C yielded only an AUC = 0.78; thus, the digital twin provided a substantial improvement on the AUC (14%). Comparing with previous MRI-based predictions of breast cancer response to NAST (17, 18, 40, 41), the digital twin also shows an improved accuracy. Both our previous study (18) and an I-SPY study (41) reported the best discrimination between pCR/non-pCR in TNBC using the optimized mid-treatment FTV, with AUC = 0.85. Moreover, the addition of post-NAST ADC to FTV showed an improvement for predicting response in TNBC, increasing AUC from 0.71 to 0.81 (17). Combining a pharmacokinetic parameter (i.e., k_ep) with ADC measured after one-cycle NAST yielded an AUC = 0.88 in breast cancer (40). Moreover, the predictive accuracy of the digital twin is comparable with state-of-the-art, machine learning (ML)-based predictions. For example, Ravichandran and colleagues applied a convolutional neural network (CNN) to predict pCR from pretreatment DCE-MRI and achieved an AUC of 0.77 in a total of 166 patients with breast cancer (42). A more recent CNN-based study using both pre- and posttreatment DCE-MRI achieved an AUC of 0.91 in a cohort of 42 patients with breast cancer (43).

Importantly, our digital twin approach has several inherent advantages comparing to ML algorithms. First, ML methods rely on access to a large patient population to train the algorithm and this training dataset must include all relevant pathophysiologic features of the disease under investigation, have high-quality annotations, and be generalizable from one population to the next. In contrast, our approach calibrates a biologically based model using patient-specific data to make patient-specific predictions, which does not require population-based training or annotation labels. Second, ML can be difficult to interpret biologically due to complex modeling features. In contrast, the digital twins provide an accurate prediction not only of pCR status at the conclusion of NAST, but also of the mechanistic interpretation of the tumor development during NAST. For example, our modeling framework can capture the initial and subsequent responses to A/C as depicted in Fig. 3A to D for patients with very different response dynamics. Third, because the digital twin parameters quantify and elucidate the observed tumor response dynamics, it provides another potential application: predicting response to multiple candidate therapeutic regimens. Thus, digital twins built on quantitative imaging data can provide—early in the course of NAST—a practical way to optimize individual treatment plans and hasten truly personalized cancer care (44).

There are, however, a few areas of our digital twin that require further investigation. First, instead of utilizing DCE-MRI to inform the spatial distribution of delivered drug, coupling the tumor-response model (i.e., Eq. B) with drug delivery and/or tumor angiogenesis models (45, 46) could estimate drug distribution more accurately. Second, instead of assuming cell death is merely proportional to drug concentration (i.e., the last term on the right-hand side of Eq. B), a model accounting for detailed therapeutic mechanisms and pharmacodynamics could improve predictive accuracy (47, 48). Third, our methodology is not currently capable of predicting invasion to the lymph nodes or axilla, which limits the accuracy for predicting pathological response or residual cancer burden (see Supplementary Section S2.3 for additional analysis). Of course, building more comprehensive models necessitates either more assumptions or more measurements, so a careful balance of model complexity and predictive accuracy must be sought.

Although the dataset employed in this study is much larger and more homogeneous than our previous preliminary study (36), there are still a few points about the cohort that could affect the analysis. Framework 1 excluded four patients with no visible tumor at V2. This decision was made so we did not overestimate the accuracy of our predictions as such patients also show no visible tumor at V3, thereby resulting in a 100% accurate prediction without really testing the modeling. In addition, the absence of pathological evaluation at the intermediate time point led to the lack of “ground-truth” for Framework 1. Framework 2 excluded nine patients due to enrollment in other trials, causing the cohort to be enriched with pCR patients. This could lead to an overestimation of accuracy for differentiating pCR/non-pCR (see Supplementary Section S2.4 for additional analysis). Future efforts will seek to apply our digital twin formalism to larger patient population, with a more typical partition of pCR/non-pCR, for further validation.

Although we (and others) have worked hard to characterize errors in the measurements themselves (34, 49, 50), the processing (e.g., segmentation, registration) steps are also potential sources of error that can be propagated through the modeling pipeline and lead to bias in the prediction. A detailed investigation suggested that an un-anticipated changes in ADC values and distribution, as well as the appearance of necrotic regions, are potential sources of error (see Supplementary Section S2.5). In addition, Framework 1 involved sampling the drug decay rates, which introduces an uncertainty in the predicted tumor response and limits the accuracy for predicting final pathology (see Supplementary Section S2.6). However, compared with previous attempts of simultaneously calibrating the drug efficacy and decay (36), this procedure not only ensures a more robust model calibration, but also allows the uncertainty in tumor dynamics to be quantified and interpreted. Another source of uncertainty is the setting of drug efficacy and decay rate of paclitaxel in Framework 2, due to the lack of imaging during the paclitaxel portion of NAST. One solution is to incorporate more measurements during NAST, especially after alternating the therapies, so that the digital twin can be updated to preserve an accurate prediction. Of course, it is important to note the goal of the digital twin is not to provide a perfect reproduction of the patient's situation. Rather, a realistic goal is to provide an accurate and practical formalism that provides clinically actionable insights. In the present contribution, we have achieved this goal in the context of predicting the response of early-stage TNBC to neoadjuvant systemic chemotherapy.

This study also supports the value of longitudinal MRI in cancer care. Currently, only pretreatment MRI is standard for assessing patients with breast cancer. While increasing to multiple MRIs will increase the cost for imaging, the benefit of early detection of, e.g., chemoresistance would help to avoid unnecessary toxicity and costs. Both this and prior studies showed that follow-up imaging after the first one to three cycles (15–18, 40, 41), of NAST is helpful for the early prediction of response. Including even one follow-up MRI (enough to enable model calibration) would be of great benefit.

We have developed patient-specific digital twins via integrating longitudinal, multiparametric MRI data with biologically based mathematical modeling. This technique accurately captures the spatiotemporal response of TNBC to NAST and achieves high accuracy and specificity for predicting the final pathologic status for each individual patient. The success of this approach demonstrates the potential of digital twins to shift the paradigm from assessing response to predicting and, eventually, optimizing response.

C. Wu reports grants from NIH and grants from CPRIT during the conduct of the study; also has a patent for U.S. Provisional Patent Application No. 63/247233 filed September 22, 2021, pending. A.M. Jarrett reports grants from NIH and CPRIT during the conduct of the study; also has a patent for U.S. Provisional Patent Application No. 63/247233 filed September 22, 2021, pending. L. Huo reports other support from MD Anderson Cancer Center during the conduct of the study. J.K. Litton reports grants from Genentech, AstraZeneca, and Novartis during the conduct of the study; grants from Zenith, GSK, EMD-Serono, and Pfizer/Medivation outside the submitted work. C. Yam reports grants from Amgen, Merck, GSK, Conquer Cancer Foundation, and Genentech outside the submitted work. J. Ma reports grants from GE Healthcare and Siemens Healthineers outside the submitted work; also has a patent for UT MD Anderson Cancer Center issued, licensed, and with royalties paid from GE Healthcare and a patent for UT MD Anderson Cancer Center issued, licensed, and with royalties paid from Siemens Healthineers; and is a consultant for C4 Imaging, LLC. T.E. Yankeelov reports grants from NIH and CPRIT during the conduct of the study; also has a patent for U.S. Provisional Patent Application No. 63/247233 filed September 22, 2021, pending. No disclosures were reported by the other authors.

C. Wu: Conceptualization, resources, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. A.M. Jarrett: Conceptualization, resources, software, funding acquisition, investigation, methodology, writing–review and editing. Z. Zhou: Resources, data curation, formal analysis, investigation, methodology, writing–review and editing. N. Elshafeey: Resources, data curation, formal analysis, investigation, methodology, writing–review and editing. B.E. Adrada: Data curation, investigation, writing–review and editing. R.P. Candelaria: Data curation, investigation, writing–review and editing. R.M. Mohamed: Data curation, investigation, writing–review and editing. M. Boge: Data curation, investigation, writing–review and editing. L. Huo: Data curation, investigation, writing–review and editing. J.B. White: Data curation, investigation, writing–review and editing. D. Tripathy: Data curation, investigation, writing–review and editing. V. Valero: Data curation, investigation, writing–review and editing. J.K. Litton: Data curation, investigation, writing–review and editing. C. Yam: Data curation, investigation, writing–review and editing. J. Son: Data curation, investigation, writing–review and editing. J. Ma: Investigation, methodology, writing–review and editing. G.M. Rauch: Funding acquisition, investigation, methodology, writing–review and editing. T.E. Yankeelov: Conceptualization, supervision, funding acquisition, investigation, methodology, writing–original draft, writing–review and editing.

The authors thank the NIH for funding through NCI U01CA142565, U01CA174706, and U24CA226110. They thank the Cancer Prevention and Research Institute of Texas for support through CPRIT RR160005. T.E. Yankeelov is a CPRIT Scholar in Cancer Research. The authors thank the Oncological Data and Computational Sciences collaboration between the Oden Institute for Computational Engineering and Sciences at The University of Texas at Austin, the MD Anderson Cancer Center, and the Texas Advanced Computing Center for providing seed funding on this project. The authors thank the generous philanthropic contributions of the University of Texas MD Anderson Cancer Center Moon Shots Program. J.K. Litton discloses the following financial relationships: Grant or research support from Novartis, Medivation/Pfizer, Genentech, GSK, EMD-Serono, AstraZeneca, Medimmune, Zenith, Merch; Speaker's Bureau for MedLearning, Physician's Education Resource, Prime Oncology, Medscape, Clinical Care Options, Medpage; Royalty from UpToDate.

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.