Article Text

Original research
Quantitative systems pharmacology model predictions for efficacy of atezolizumab and nab-paclitaxel in triple-negative breast cancer
  1. Hanwen Wang1,
  2. Huilin Ma1,
  3. Richard J Sové1,
  4. Leisha A Emens2 and
  5. Aleksander S Popel1,3
  1. 1Department of Biomedical Engineering, Johns Hopkins School of Medicine, Baltimore, Maryland, USA
  2. 2Department of Medicine, University of Pittsburgh Medical Center, Hillman Cancer Center, Pittsburgh, Pennsylvania, USA
  3. 3Department of Oncology, Johns Hopkins Medicine Sidney Kimmel Comprehensive Cancer Center, Baltimore, Maryland, USA
  1. Correspondence to Mr Hanwen Wang; hwang163{at}


Background Immune checkpoint blockade therapy has clearly shown clinical activity in patients with triple-negative breast cancer, but less than half of the patients benefit from the treatments. While a number of ongoing clinical trials are investigating different combinations of checkpoint inhibitors and chemotherapeutic agents, predictive biomarkers that identify patients most likely to benefit remains one of the major challenges. Here we present a modular quantitative systems pharmacology (QSP) platform for immuno-oncology that incorporates detailed mechanisms of immune–cancer cell interactions to make efficacy predictions and identify predictive biomarkers for treatments using atezolizumab and nab-paclitaxel.

Methods A QSP model was developed based on published data of triple-negative breast cancer. With the model, we generated a virtual patient cohort to conduct in silico virtual clinical trials and make retrospective analyses of the pivotal IMpassion130 trial that led to the accelerated approval of atezolizumab and nab-paclitaxel for patients with programmed death-ligand 1 (PD-L1) positive triple-negative breast cancer. Available data from clinical trials were used for model calibration and validation.

Results With the calibrated virtual patient cohort based on clinical data from the placebo comparator arm of the IMpassion130 trial, we made efficacy predictions and identified potential predictive biomarkers for the experimental arm of the trial using the proposed QSP model. The model predictions are consistent with clinically reported efficacy endpoints and correlated immune biomarkers. We further performed a series of virtual clinical trials to compare different doses and schedules of the two drugs for simulated therapeutic optimization.

Conclusions This study provides a QSP platform, which can be used to generate virtual patient cohorts and conduct virtual clinical trials. Our findings demonstrate its potential for making efficacy predictions for immunotherapies and chemotherapies, identifying predictive biomarkers, and guiding future clinical trial designs.

  • computational biology
  • systems biology
  • immunotherapy
  • breast neoplasms
  • tumor microenvironment

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article and the supplementary material. MATLAB scripts for model generation and in silico clinical trials are publicly available at and on GitHub (

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See

Statistics from

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.


Triple-negative breast cancer (TNBC) is an aggressive subtype of breast cancer that lacks expression of estrogen receptor, progesterone receptor, and human epidermal growth factor receptor-2 (HER-2).1 In the past decade, immune checkpoint blockade therapies have shown promising efficacy in patients with melanoma and non-small cell lung cancer, but there have been no approved immunotherapy strategies for patients with TNBC until recently.2 In 2019, the combination of a monoclonal anti-PD-L1 antibody, atezolizumab, and a chemotherapeutic agent, nanoparticle albumin-bound (nab)-paclitaxel, received accelerated approval from the US Food and Drug Administration as first-line treatment for patients with unresectable locally advanced or metastatic TNBC that is PD-L1 immune cell positive, based on results from the IMpassion130 trial (NCT02425891).3 4 Whereas response rates for patients with advanced TNBC ranged from 0% to 24% with single-agent atezolizumab, the IMpassion130 trial reported a response rate of 55% with the combination of atezolizumab and nab-paclitaxel.5 Importantly, the trial showed clinically meaningful improvements in median progression-free and overall survival in PD-L1 immune cell (IC)-positive patients with the addition of atezolizumab to nab-paclitaxel, identifying PD-L1 expression on tumor-infiltrating ICs as a clinically relevant biomarker associated with benefit from this combination regimen.3 4 While the results of IMpassion130 established immunotherapy as a treatment option for some patients with advanced TNBC, it is critical to investigate additional cytotoxic agents to be combined with atezolizumab and other immunotherapies, and further evaluate biomarkers in both PD-L1-positive and PD-L1-negative patient cohorts.6

Nab-paclitaxel is a taxane-based chemotherapeutic agent, which is commonly used as a first-line treatment for patients with metastatic breast cancer.7 The addition of albumin facilitates the transport of nab-paclitaxel into the tumor interstitium through albondin (gp60), an albumin receptor, on endothelial cells. Furthermore, the overexpression of osteonectin (SPARC) on TNBC cells, which has a high affinity with albumin, enhances the uptake of nab-paclitaxel and the release of paclitaxel near the cancer cells.8 9 These mechanisms are thought to contribute to the higher intratumoral concentration of nab-paclitaxel compared with paclitaxel and its other derivatives in TNBC, which was confirmed by population pharmacokinetic (PK) models and preclinical studies.10 11 Various pharmacodynamic (PD) effects of nab-paclitaxel on TNBC have been observed in preclinical studies, including its cytotoxicity towards cancer and endothelial cells, and upregulation of vascular endothelial growth factor A (VEGF-A) expression, which provides a rationale for combination therapies with anti-VEGF agents such as bevacizumab.12–15 In IMpassion130, a dose of 100 mg/m2 nab-paclitaxel was administered days 1, 8 and 15 every 28 days along with 840 mg atezolizumab or placebo days 1 and 15 to 902 patients with advanced TNBC.3 4 This combination was supported by the following rationale: the cytotoxic activity of nab-paclitaxel enhances tumor-specific antigen release, with and uptake by and subsequent maturation of antigen-presenting cells (APCs) in the tumor and tumor-draining lymph nodes (TDLNs), and atezolizumab reverses T cell suppression by selectively blocking interactions between PD-L1 and programmed death 1 (PD-1) on ICs and cancer cells, as PD-L1 can be expressed in TNBC.16

To accompany clinical efforts to develop therapies and informative biomarkers that improve clinical outcomes, multiple quantitative systems pharmacology (QSP) models have been developed using various approaches.17–19 Hardiansyah et al19 proposed a model to simulate cellular kinetics and cytokine profiles in patients with chronic lymphocytic leukemia receiving chimeric antigen receptor T cell therapy, which was trained by patient-level data. Betts et al18 presented a model to predict the efficacy of CD3-bispecific antibody therapy using in silico, in vivo and in vitro data across species. Gong et al17 built a spatial model to simulate interactions between cancer cells and stromal cells and predict the efficacy and biomarkers of immunotherapy based on digital pathology data. In our previous studies, we developed an ordinary differential equation-based QSP platform to predict results of a clinical trial in HER-2 negative breast cancer using anti-PD-1 and anti-CTLA-4 antibodies with an epigenetic modulator, using in vitro and in vivo data from preclinical studies.20 By expanding our previous model,20 we conducted a virtual clinical trial of nab-paclitaxel and atezolizumab and compared our simulated results with published population-level data from IMpassion130. The virtual clinical trial aims to generate a virtual patient cohort that can be mapped to a clinical cohort in a specific trial. While there exist various methodologies in the virtual patient generation, the optimization of these algorithms is under active investigation.21–25 In this study, we aim to determine the relationship between our QSP platform and the virtual cohort with the patient cohort and results of the clinical trial. We discuss the limitations related to our choice of methodology, which need to be taken into account while interpreting the present numerical results and comparisons.


Model overview

The present model is adapted from our previously published QSP platform20 using the SimBiology toolbox in MATLAB (MathWorks, Natick, Massachusetts, USA), which comprises four compartments: central, peripheral, tumor, and TDLNs. As a modular model, it was introduced with eight modules that describe the kinetics and dynamics of effector T cells (Teff), regulatory T cells (Treg), cancer cells, APCs, tumor-specific neoantigens and tumor-associated self-antigens, immune checkpoints, myeloid-derived suppressor cells (MDSCs), and therapeutic agents.20 In this study, we added a chemotherapeutic module for nab-paclitaxel, modified tumor growth and PD-L1 dynamics, and further expanded the T cell module for naïve T cell kinetics in this study. The dynamics of the major species in the model are illustrated in figure 1A, including: immune activation of naïve CD8+ and CD4+ T cells in the TDLNs by tumor-specific neoantigens and tumor-associated self-antigens; T cell trafficking throughout the body; immune evasion mediated by immune checkpoints, Treg, MDSCs, and anti-inflammatory cytokines; and PK/PD of drugs of interest. The modules were written using MATLAB scripts and incorporated 271 parameters, 140 ODEs, and 54 algebraic equations in total. Details of all the modules are elaborated in the online supplemental methods, and full lists of model parameters, reactions, algebraic equations, and cellular and molecular species are presented in the online supplemental tables S1-S6).

Supplemental material

Supplemental material

Figure 1

QSP model diagram (A) and workflow (B). The model is composed of four compartments: central, peripheral, tumor, and tumor-draining lymph node, which describe cycles of immune activation in lymph nodes, T cell trafficking to the tumor, killing of cancer cells, immune evasion, and antigen release and lymphatic transport. ARG-I, arginase I; AT, activated T cell; MAPC, mature antigen presenting cell; NO, nitric oxide; QSP, quantitative systems pharmacology; NT, naïve T cell; TEFF, effector T cell; TH, T helper cell; Treg, regulatory T cell. Modified from ref 47.

Model initiation and simulation settings

The model workflow is illustrated in figure 1B with five steps: (1) load baseline parameters into the MATLAB workspace and initialize the SimBiology model with initial compartments and simulation settings; (2) add modules that include model species, rules, and events by calling the corresponding MATLAB functions; (3) simulate to find the initial condition with baseline parameters and a desired initial tumor diameter; (4) create SimBiology dose objects for drugs of interest and simulate therapy; and (5) generate a number of parameter sets to perform virtual clinical trials and subsequent analyses. The baseline parameters are estimated or fitted based on in vitro and in vivo experimental data generated using TNBC cell lines from literature (online supplemental methods). After model initialization with baseline parameters and required modules, a simulation is performed with a small number of cancer cells and a desired initial tumor diameter. The desired initial tumor diameter corresponds to the pretreatment tumor size of a patient in the clinical trial. Once the initial tumor diameter is reached, the amount of model species and parameter values are saved as the initial condition, which corresponds to the pretreatment condition of the patient in the clinical trial. The model is then reinitialized with the previously saved initial condition and simulated with drugs of interest. A sample of model outputs is presented in online supplemental figure S1.

Virtual patient generation and virtual clinical trial

The virtual clinical trial, or in silico clinical trial, aims to predict the efficacy of drugs of interest in patients with particular conditions (eg, melanoma and TNBC) using the virtual patient cohort.22 26 To generate a virtual patient cohort that resembles the clinical population, a subset of model parameters is selected to be varied, while others are kept at the baseline level. The selected parameters, such as antigen binding affinity, the number of tumor-specific T cell clones, initial tumor diameter, and PK/PD parameters of nab-paclitaxel aim to capture the interindividual variabilities in real patients. The distributions of these parameters are estimated based on available literature data on TNBC. With the majority of the model parameter values estimated based on in vitro and in vivo data, we use data from clinical trials to further calibrate the distribution of varied model parameters. Specifically, the published clinical results from NCT01375842 and the placebo comparator arm of the IMpassion130 trial are used for model calibration, and the results from the experimental arm of the IMpassion130 trial are used for model validation.

To perform a virtual clinical trial, the model is first initialized as described by steps 1–2. The values of selected parameters are randomly generated based on the calibrated parameter distributions using Latin Hypercube Sampling (LHS), with each parameter set representing a potential virtual patient. Each randomly assigned parameter set is then plugged into the model to simulate tumor response to the therapy following steps 3–4. For postprocessing after the simulation, the potential virtual patients are filtered by the acceptance criteria based on their pretreatment T cell densities so that their values fall into the physiologically plausible range in patients with TNBC,21 which are reported by clinical measurements.27 28 For comparison purpose, the two-dimensional (2-D) T cell density reported by the clinical measurements are converted to three-dimensional (3-D) density using the equation from Mi et al.29 In this study, the total number of 900 virtual patients are generated on calibration for efficacy predictions and statistical analyses. The study serves as the first retrospective analysis of the IMpassion130 trial using virtual patient cohorts.4 To match the clinical settings, the simulation time is set to be 400 days, which corresponds approximately to the median follow-up time of 12.9 months in the IMpassion130 trial.4 Additionally, the simulated tumor diameters are recorded every 8 weeks, which corresponds to the frequency of tumor measurements in the clinical trial.

Statistical analysis

Global uncertainty and sensitivity analysis are performed with LHS and Partial Rank Correlation Coefficient methods to examine the impact of varied parameters on model observations.30 The objective response rate (ORR) and duration of response (DOR) are predicted based on RECIST V.1.1,31 and the 95 percentile bootstrap CIs are calculated for comparison between model predictions and clinical results. In subgroup analysis, 95% Agresti-Coull CIs are estimated for the ORR predictions based on the normal approximation for the binomial distribution. For comparison of model observations in subgroups of different response status and treatment regimens, the Wilcoxon test is performed using ggpubr package in RStudio V.1.2.32


Efficacy predictions of the PD-L1 inhibitor

In the previous study, we made prospective predictions of response rates for PD-1 and CTLA-4 inhibitors, and an epigenetic modulator, using the same dose regimens of an ongoing clinical trial in patients with HER-2 negative breast cancer at the time.20 Our predicted response rate of anti-PD-1 monotherapy fell into the range reported by multiple clinical trials of PD-1 inhibitors, such as nivolumab and pembrolizumab, in patients with breast cancer.5 To calibrate our modifications made in this study, we first conduct a virtual clinical trial using a PD-L1 inhibitor, atezolizumab, in virtual patients with TNBC and compare with the results from a phase I clinical trial reported by Emens et al.33 In this phase I trial, 116 patients with metastatic TNBC were enrolled and received atezolizumab treatment. Among the 116 patients, 93 patients received 15 mg/kg atezolizumab, 1 patient received 20 mg/kg atezolizumab and 22 patients received 1200 mg atezolizumab, every 3 weeks. One hundred fifteen of the 116 patients had evaluable disease and were included to calculate the ORR by RECIST V.1.1.31 In summary, 11 patients had partial or complete response, 15 patients had stable disease, 73 patients had progressive disease and 16 patients died within 6 weeks after treatment started. The results correspond to an ORR of 10%, a percentage of stable disease of 13%, and a percentage of progressive disease of 77%.

To perform a virtual clinical trial that is comparable with the phase I trial, we aim to reproduce the clinical settings in our model simulations. Since most of the patients received 15 mg/kg or 1200 mg atezolizumab, we use 1200 mg atezolizumab every 3 weeks in the virtual clinical trial due to a lack of information of patients’ body weights. On calibration of the parameter distributions in the virtual patient generation, we generate 900 virtual patients based on our filter criteria specified in the Methods. To match the number of patients enrolled in the clinical trial, we perform 1000 bootstrap resampling of 99 virtual patients from the virtual patient cohort and calculate the bootstrap medians and the 95 percentile bootstrap CIs for efficacy predictions. The resulting median ORR is 13.1% with a CI from 6.1% to 20.2%; the median percentage of stable disease is 64.7% with a CI from 52.5% to 73.7%; and the median percentage of progressive disease is 23.2% with a CI from 14.1% to 32.3%. By comparison, the clinically observed ORR from the phase I trial falls into the bootstrap CI reported by the model prediction. However, the model overestimates the percentage of stable disease and underestimates the percentage of progressive disease when compared with clinical observations. This shift of model prediction from progressive disease to stable disease is, at least partially, due to the lack of model prediction for metastasis. As reported by Emens et al,33 a number of patients are considered to have progressive disease due to newly detected lesions (ie, newly detected metastases), even with a reduction of overall tumor burden.

To better virtualize the difference between the model predictions and the clinical results, we plot the rate of response and the best overall response as the spider plot and waterfall plot, respectively, in figure 2A. The figure demonstrates the capability of virtual clinical trials to capture the interindividual variabilities, as the tumor dynamics of the virtual patient cohort resemble those reported by the clinical trial. To further explore the results, we plot the distributions of potential predictive biomarkers in figure 3 and compare their differences between responders and non-responders. As shown in the figure, the number of tumor-specific T cell clones, binding affinity (KD) of neoantigens, CD8+ and CD4+ T cell densities, CD8+/Treg and CD4+/Treg ratios are significantly higher in responders, which shows consistency with our previous model predictions and clinical evidence.6 20 34–36

Figure 2

Rate of response (left) and the best overall response (right) in model-predicted tumor diameter of 100 randomly selected virtual patients. Response is assessed by RECIST V.1.131 in atezolizumab monotherapy (A), nab-paclitaxel group (B), and atezolizumab+nab-paclitaxel group (C). Median (thick lines) and individual (thin line) rate of response are shown in PD (red), SD (purple), and PR/CR (blue) subgroups. CR, complete response; nab, nanoparticle albumin-bound; PD, progressive disease; PR, partial response; SD, stable disease.

Figure 3

Pretreatment distributions of potential predictive biomarkers in responders and non-responders. Statistical significance is calculated by Wilcoxon test. Atezo, atezolizumab monotherapy 1200 mg every 3 weeks; Combo, atezolizumab 840 mg every 2 weeks+nab-paclitaxel 100 mg/m2 Q3/4W; MDSC, myeloid-derived suppressor cell; Nab-P, nab-paclitaxel 100 mg/m2 Q3/4W; NR, non-responders; R, responder.

Virtual clinical trial of atezolizumab and nab-paclitaxel

To specifically calibrate the parameters related to chemotherapy, we performed another virtual clinical trial using the treatment regimens from the IMpassion130 trial. The trial contains two treatment arms: placebo plus nab-paclitaxel (placebo comparator) arm and the atezolizumab plus nab-paclitaxel (experimental) arm. In the placebo comparator arm, 100 mg/m2 dose is administered days 1, 8, and 15 every 28 days (Q3/4W), and in the experimental arm, a placebo or 840 mg atezolizumab is administered days 1 and 15 every 28 days (Q2W). As reported by Schmid et al,4 the ORRs of nab-paclitaxel with placebo and nab-paclitaxel in combination with atezolizumab are 45.9% and 56.0%, and the DOR for the two regimens are 5.6 and 7.4 months, respectively. To simulate the efficacy of nab-paclitaxel, we introduce a new module to the QSP platform, which incorporates the current knowledge of mechanisms of action of nab-paclitaxel. In summary, the nab-paclitaxel PK is adapted from the published PK model by Chen et al,10 and its PD is calibrated to experimental and clinical data (online supplemental figure S2). By varying the PK/PD parameters in the virtual patient generation, we are able to account for the interindividual variabilities in the drug delivery and resistance.

Similar to the atezolizumab monotherapy, the virtual patient cohort is generated on calibration to perform the virtual clinical trial. As shown in figure 2B and C, the rate of response and the best overall response of the randomly selected virtual patients are plotted. The tumor dynamics and the shape of the waterfall plots also resemble those reported by clinical trials of atezolizumab and nab-paclitaxel.37–39 To better compare the model predictions of ORR and DOR with clinical results, we generate an efficacy prediction table using a similar format to that reported by Schmid et al.4 In table 1, we perform 1000 bootstrap resampling of 450 virtual patients from the whole virtual cohort (900 virtual patients in total) and calculate the bootstrap medians and the 95 percentile bootstrap CIs for the endpoint predictions in each treatment regimen. While most of the efficacy predictions overlap with the ranges reported by the clinical trial (table 1; online supplemental figure S2), the complete response rate is underestimated. Although we predict a tumor smaller than 2 mm as a complete response by assuming limited detectability by imaging modalities, a tumor that is smaller than 2 mm can be detected in some cases.40 Additionally, the inhibitory effect of MDSC, Treg, immune checkpoints, and cytokines on Teff may be overestimated by the model, which together lead to the underestimation of complete tumor eradication.

Table 1

Efficacy prediction for the virtual patient cohort generated based on calibrated parameter distribution

As shown in figure 3, the distributions of potential predictive biomarkers show different trends in the three treatment regimens. Unlike atezolizumab monotherapy, the pretreatment MDSC density, the number of tumor-specific T cell clones, neoantigen KD, and CD4+/Treg ratio are not significantly affected by the response status. This prediction suggests that the efficacy of chemotherapy may not depend on the reversal of the inhibition on pre-existing Teffs, but on the therapy-induced activation of immune response by newly released tumor antigens. Interestingly, figure 3 shows that the pretreatment CD8+/Treg ratio is significantly lower in the responders of nab-paclitaxel monotherapy, which does not match our expectation. In figure 4, we divide the virtual patient cohort into subgroups by their pretreatment values of selected biomarkers and calculate the corresponding ORR with the 95% Agresti-Coull CI. The CIs for subgroups based on the CD8+ and CD4+ T cell levels show significant differences in the corresponding response rates, while those for other subgroups overlap. These results suggest that the response to combination therapy of PD-L1 inhibitors and nab-paclitaxel correlates with high T cell levels, although this was not observed for CD8+ T cells in IMpassion 130.41 42 In addition, subgroups with high PD-L1 expression and CD4+/Treg ratio in the tumor also correspond to notably higher ORRs than those with low PD-L1 expression and CD4+/Treg ratio.

Figure 4

Subgroup analysis of the combination therapy in virtual patient cohort. The total 900 virtual patients are divided into eight subgroups based on the pretreatment values of selected biomarkers, and the objective response rates in each subgroup are calculated with 95% Agresti-Coull CIs. MDSC, myeloid-derived suppressor cell.

Performance of potential predictive biomarkers

The performance of the predictive biomarkers identified previously is investigated using a binary classification model. As shown in figure 5, the sensitivity and 1-specificity values from each cut-off were plotted as receiver operating characteristic (ROC) curves. CD8+ and CD4+ T cell densities have higher areas under curves (AUCs) (0.631 and 0.659, respectively) than PD-L1 expression and CD8+/Treg and CD4+/Treg ratios (0.570, 0.507, and 0.554, respectively), further implicating their potential to be predictive biomarkers for this double combination regimen. However, the model-predicted AUCs are not as high as those reported in our previous simulations and clinical analyses for anti-PD-1/PD-L1 monotherapy.20 43 Furthermore, the dynamics of selected biomarkers are investigated under the three treatment regimens in online supplemental figure S3. The post-treatment (week 8) to pretreatment ratios suggest that either atezolizumab or nab-paclitaxel alone is able to significantly reduce tumor volume and increase both CD8+ and CD4+ densities in tumor within 8 weeks on drug administration. The addition of atezolizumab to nab-paclitaxel further significantly increases tumor volume reduction and T cell levels, which shows an additive effect of the double combination therapy. Interestingly, the CD8+/Treg ratio in the tumor is significantly lowered by nab-paclitaxel at week 8, which agrees with the in vivo observations and needs to be validated by future experimental or clinical results.44

Figure 5

ROC analysis of potential predictive biomarkers in combination therapy. Cut-off values are selected based on the range of PD-L1 molecules on APCS, pretreatment effector T cell density, tumor mutational burden, and Teff to regulatory T cell ratio. For each cut-off value, response status (R vs NR) is predicted for each virtual patient by comparing the pretreatment amount of the potential predictive biomarker to the cut-off value. Sensitivity (true positive rate) is plotted against 1 − specificity (true negative rate) for each biomarker. APCs, antigen-presenting cells; AUC, areas under curve; NR, non-responders; R, responders; ROC, receiver operating characteristic.

To further explore the effect of the biomarkers on model predictions, we perform the global uncertainty and sensitivity analysis using the calibrated virtual patient distribution. As shown in online supplemental figure S4, 26 parameters are varied in the virtual patient generation, including parameters with fitted range to experimental data and those with estimated range that are then calibrated by the clinical results. Additionally, for a subset of parameters varied in the virtual patient generation, we sort the virtual patients by each parameter in ascending orders and evenly divide them into multiple subgroups. The ORR of each subgroup in the double combination therapy is plotted against the median parameter values in figure 6. As a result, predicted ORRs show a trend of increase as PD-L1 expression, CD8+ and CD4+ T cell densities increase, while ORRs decrease as tumor growth rate and pretreatment tumor size increase.

Figure 6

Effects of parameters on objective response. For each parameter of interest, 900 virtual patients are sorted by the parameter values in ascending order and evenly divided into nine subgroups. The response status of each subgroup in the combination therapy is plotted against the corresponding median parameter values. MDSC, myeloid-derived suppressor cells.

Sequential therapy simulation

Now that the present QSP platform is validated by its prediction of the additive effect of atezolizumab to nab-paclitaxel on ORR and DOR (table 1; online supplemental figure S2), we aim to explore its potential for dose optimization. To this end, we perform a series of virtual clinical trials using various doses and schedules of nab-paclitaxel and atezolizumab. Specifically, 840 mg atezolizumab is administered every 2 weeks starting on day 1, week 2, or week 4 on reaching the initial tumor diameter, in combination with nab-paclitaxel, and 100 mg/m2 (on days 1, 8, and 15 of a 28-day cycle), 125 mg/m2 (on days 1 and 8 of a 21-day cycle), or 260 mg/m2 (every 3 weeks) of nab-paclitaxel is administered staring on day 1, week 2, or week 4 (figure 7A–C, online supplemental S5A–C). These selected dose schedules have been used in clinical trials of combination therapies of nab-paclitaxel and other therapeutic agents, which are adopted here to avoid potential toxicity concerns.7 The median tumor volume, CD8+ T cell level, and Treg density, CD8+/Treg, and CD4+/Treg ratios in the tumor at week 8 are reported with the ORR for each combination of the treatment regimens in figure 7 and online supplemental figure S5. Although 260 mg/m2 nab-paclitaxel results in the highest predicted ORR in the concurrent therapy (online supplemental figure S5), more frequent dosing with 100 mg/m2 nab-paclitaxel leads to the highest median tumor volume reduction and CD8+ T cell level at week 8 (figure 7). Notably, simulating concurrent therapy with atezolizumab and 260 mg/m2 nab-paclitaxel shows high Treg density in the tumor (figure 7), which corresponds to the lowest CD8+/Treg and CD4+/Treg ratios among the combination regimens (online supplemental figure S5). This upregulation of Treg density was also observed in a preclinical in vivo study of murine TNBC treated with paclitaxel.44 Overall, the virtual clinical trials suggest that concurrent therapies of atezolizumab and nab-paclitaxel lead to better response than sequential therapies, and the model predictions agree with the clinical observation that weekly 100 mg/m2 nab-paclitaxel shows similar efficacy to 300 mg/m2 every 3 weeks dosing but with reduced toxicity.45

Figure 7

Model simulation of sequential therapies using various nab-paclitaxel doses and schedules. Top row (A–C) represents the median tumor volume after 8 weeks of each dose regimen; middle row (D–F) represents the median CD8+ T cell density at week 8; and bottom row (G–I) represents the median Treg density in the tumor at week 8. Administration of nab-paclitaxel starts on day 1 (A,D,G), week 2 (B,E,H), and week 4 (C,F,I) on reaching initial tumor diameter.


In this study, we introduce multiple modifications to our previously developed QSP platform, including a new nab-paclitaxel module and modified cancer, T cell, and checkpoint modules.20 46 47 The addition of nab-paclitaxel begins our attempt to make efficacy prediction for combination therapies involving checkpoint inhibitors and chemotherapeutic agents. While nab-paclitaxel has been approved for the treatment of metastatic breast cancer, its efficacy when combined with other types of therapeutic agents is still under investigation. Here, we implement all the known effects of nab-paclitaxel on TNBC cell lines, including its cytotoxic, angiogenic, and antiangiogenic activities, as reported by in vivo preclinical and clinical observations.12–15 As we incorporate all the mechanisms of action into the model, we are able to investigate the overall effect of nab-paclitaxel on tumor dynamics. In addition to the new module, the available clinical data of T cell levels in patients with TNBC allow us to expand the T cell module to better describe the dynamics of naïve T cells, CD4+ helper T cells, and the interactions within the T cell subsets, such as Th to Treg transdifferentiation. As shown in figure 3, CD8+, CD4+, and Treg densities in the tumor all fall within the physiologically realistic ranges reported by clinical measurements.27–29 Furthermore, the upregulation of PD-L1 expression by interferon gamma (IFN-γ) in the tumor is implemented into the checkpoint module to account for the negative feedback mechanism that follows the tumor infiltration of lymphocytes.48 This mechanism allows us to better predict the dynamics of PD-L1 expression during the therapy and how it is associated with the response to the combination therapy.

In this retrospective clinical trial analysis using our proposed QSP platform, we demonstrate its potential for making efficacy predictions of checkpoint inhibitors and chemotherapeutic agents by conducting virtual clinical trials. When the distributions of model parameters are calibrated by data from the two monotherapies, the model-predicted efficacy of atezolizumab in combination with nab-paclitaxel overlaps with the experimental arm of the IMpassion130 trial.4 However, the comparisons between our model predictions and the clinical results also reflect the limitations of the present model. Particularly, the efficacy prediction for atezolizumab monotherapy shows a shift from progressive disease to stable disease. This discrepancy is due to the complexity of the immune system and the clinical settings that cannot be entirely captured or reproduced by model simulations. Using model simulations, we identify the objective response status of the virtual patients by their best overall response assuming a spherical tumor. In other words, virtual patients with a reduction of tumor size more than 30% at any time point of the simulation are considered as partial or complete responders; those with an increase of tumor size more than 20% within 8 weeks on the dose administration are considered to have progressive disease; and the rest are considered to have stable disease.31 In reality, the characterization of response status is also impacted by other conditions, such as patients’ survival, new metastatic lesions, the definition of tumor burden as the sum of largest diameters, and resolution of the imaging modalities, all of which are likely to cause a deviation from model predictions.

In the biomarker analysis, the model first confirms that the number of tumor-specific T cell clones (as a measure of tumor mutational burden) and T cell densities in the tumor are associated with response status in single-agent PD-L1 blockade therapies.6 34–36 For the combination of atezolizumab and nab-paclitaxel, the model confirms that subgroups with high PD-L1 expression have higher ORRs, even though the increase is not as significant as reported by the trial. Besides, the model identifies that CD8+ and CD4+ T cell levels are the best two predictive biomarkers due to their significantly higher medians in responders and the significantly higher ORR in subgroups with high levels of T cells (figures 3 and 4). Although similar correlations are observed in neoadjuvant chemotherapies using other cytotoxic agents,41 42 the performance of these biomarkers needs to be explored in combination treatments. Nonetheless, there exist discrepancies between model-predicted biomarkers and clinical observations in chemotherapy. As shown in figure 3, the pretreatment CD8+/Treg ratio is significantly lower in the responders of nab-paclitaxel monotherapy, but this correlation has not been confirmed by clinical studies. In fact, pathological complete response rate is reported to be significantly higher in patients with early-stage TNBC in high-CD8+/FoxP3+ ratio group during FEC100 treatment.49 Although this discrepancy is possibly due to the differences in treatment regimens and immunobiology between early-stage and metastatic TNBC, the relationship between CD8+/Treg ratio and response status in nab-paclitaxel treatment remains to be confirmed by clinical studies.

Furthermore, due to the promising synergistic effect of sequential therapies using checkpoint inhibitors and other therapeutic agents in multiple cancer types, we performed a series of in silico clinical trials to investigate the optimal dose schedule of atezolizumab and nab-paclitaxel for TNBC. The results suggest that concurrent therapies result in higher ORRs and T cell densities than sequential therapies. In addition, although 260 mg/m2 every 3 weeks nab-paclitaxel corresponds to the highest ORR, more frequent doses with 100 mg/m2 Q3/4W (on days 1, 8, and 15 of a 28-day cycle) have similar ORR and higher CD8+ T cell levels with potentially less toxicity.50 Clinical data from other trials involving checkpoint inhibitors and chemotherapy in TNBC would be valuable for model calibration and validation, such as IMpassion132 and KEYNOTE-355. Notably, the simulated dose regimens are adopted from the existing clinical trial protocols whose safety has been tested clinically to avoid potential toxicity concerns. Future incorporation of toxicity prediction into the QSP platform will allow optimization of dose intensity and frequency for atezolizumab. In a recently published study, Ma et al51 suggest cytokine release syndrome, which is a common dose-dependent adverse event caused by cancer immunotherapy, to be predicted by mechanistic modeling of interleukin-6 expression. However, this method requires mechanistic understandings of immune-related adverse effects and patient-level data from clinical trials for model implementation, which is currently unavailable for atezolizumab.

Overall, the predictive power of virtual clinical trials depends on three factors: translation of molecular and cellular mechanisms into model reactions, parameter estimation using experimental data, and calibration of the virtual patient distribution. First, the QSP model aims to incorporate the current knowledge of molecular and cellular mechanisms to make accurate predictions based on preclinical observations. However, the phenotypic observations may not reflect the true mechanisms behind them. In many cases, hypotheses have to be made to incorporate the observed dynamics into the model. In addition, the degrees of mechanistic detail are limited by the scope of the model. For example, the development of chemoresistance is known to be related to cytokines, membrane proteins, gene expression changes, and many other pathways in cancer cells.52 Since it is not pragmatic to incorporate all known mechanisms into the model in every detail, a constant effort to optimize the level of model complexity is essential. Second, the parameter estimation remains one of the major challenges when building QSP models, given their high levels of complexity. In this study, most of the parameters are estimated using experimental data of breast cancer, as our goal is to simulate tumor dynamics in patients with advanced TNBC. However, a number of issues remain to be explored: how to adjust parameters for spatial heterogeneity or for the tumor that metastasized from its primary tumor site? How are parameters derived from animal data translated to humans? How to account for potential covariance between the model parameters in the virtual patient generation (eg, between body surface area and PK parameters)? Third, the distributions of parameters varied in the virtual patient generation are calibrated by clinical results, which can be biased toward the inclusion criteria of the trial. Notably, there are a number of methods for the virtual patient generation with different algorithms to address potential biases.21–25 Here, we use the methods that are similar to the recently published studies.53 54 The optimal techniques for virtual patient generation based on the availability of clinical data is an active area of research that is undergoing rapid development.

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article and the supplementary material. MATLAB scripts for model generation and in silico clinical trials are publicly available at and on GitHub (

Ethics statements

Patient consent for publication


The authors would like to thank Qingqing Yin for her helpful comments on the statistical analyses that improved the manuscript.


Supplementary materials

  • Supplementary Data

    This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.


  • Twitter @HanwenWang95, @EmensLeisha, @HopkinsPopelLab

  • Correction notice This article has been corrected since it first published. The provenance and peer review statement has been included.

  • Contributors ASP, LAE, and HW designed and planned the project. HW, HM, and RJS built the model. HW modified the model, performed all simulations, analyzed the simulation data, and prepared a draft of the manuscript. ASP and LAE critically revised the manuscript. All authors have read and approved the final manuscript.

  • Funding Supported by NIH grant R01CA138264.

  • Competing interests None declared.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.