## Abstract

We present a mechanistic mathematical model of immune checkpoint inhibitor therapy to address the oncological need for early, broadly applicable readouts (biomarkers) of patient response to immunotherapy. The model is built upon the complex biological and physical interactions between the immune system and cancer, and is informed using only standard-of-care CT. We have retrospectively applied the model to 245 patients from multiple clinical trials treated with anti–CTLA-4 or anti–PD-1/PD-L1 antibodies. We found that model parameters distinctly identified patients with common (*n* = 18) and rare (*n* = 10) malignancy types who benefited and did not benefit from these monotherapies with accuracy as high as 88% at first restaging (median 53 days). Further, the parameters successfully differentiated pseudo-progression from true progression, providing previously unidentified insights into the unique biophysical characteristics of pseudo-progression. Our mathematical model offers a clinically relevant tool for personalized oncology and for engineering immunotherapy regimens.

## INTRODUCTION

Immunotherapy has emerged as a promising therapy for multiple cancers, and its utility is projected to grow as combinatorial effects with existing modalities of cancer treatment become elucidated (*1*, *2*). Notably, some patients with traditionally dire prognosis have been found to respond to immunotherapy treatment with notable efficacy as measured with current canonical metrics of tumor burden response (*3*): the gold standard of overall survival, followed by progression-free survival, and a consensus measure of tumor burden [Response Evaluation Criteria in Solid Tumors (RECIST 1.1)] (*4*, *5*). However, to date, only a minority of patients have been reported to receive clinical benefit from immunotherapy (*6*). Recognizing the need for more refined metrics of treatment response by immunotherapy, the modified iRECIST has been proposed (*7*). While iRECIST provides a guideline that helps maintain standards of interpretation of clinical data from prospective trials, it only provides prognostic value after extended monitoring (*7*), and widespread consensus on standardized guidelines for its implementation remains elusive (*8*). The cost of this delayed response feedback is that by the time immunotherapy is definitively confirmed to be of unsatisfactory benefit to the patient, it is often too late to adapt treatment to affect the patient’s outcome.

To address the need for earlier determinants of treatment response, multiple groups have investigated circulating tumor DNA (*9*), complete “-omics” (*10*) mRNA profiling (*11*), and immune cell population shifts (*12*, *13*). However, most of these metrics require either invasive procedures or technically challenging analyses, which may prove to be a hurdle for some patients and clinical centers, respectively. To achieve the goal of personalized medicine for all patients, a noninvasive and objective metric to determine who will benefit from immunotherapy is needed. To this end, we have developed a mechanistically based biophysical model that describes the immunotherapy response. In this study, we use tumor growth rate obtained from standard-of-care computed tomography (CT) imaging, and model-derived quantifications of tumor cell kill rate (μ) and antitumor immune state (Λ) to predict treatment outcomes for a mixed population of patients undergoing checkpoint inhibitor treatment for metastatic disease of various primary tumor types, and demonstrate methodologies whereby the model may be implemented prospectively as a predictor of patient response and survival at times as early as first restaging (α_{1}). Last, we offer model-derived insights into the mechanistic underpinnings of pseudo-progression, suggesting new directions for further study of this phenomenon.

## MATERIALS AND METHODS

### Assumptions about the tumor microenvironment and immune response

Immune cell presence in the tumor microenvironment is necessary for potential cancer cell kill. Immune cells (killer T cells) respond to diffusion of immunotherapy drug (antibody) molecules and consequently become activated, in the presence of immunotherapy drug binding, through cell signaling pathways and chemotaxis. The underlying model hypothesis is that these drug molecules diffuse within the tumor microenvironment and block the interaction of immune checkpoint inhibitory ligands expressed on immune T cells with their complementary proteins on tumor cells or antigen-presenting cells (APCs). This binding renders effector immune T cells potentially active for cancer cell killing and ultimately implies that the therapeutic site of action occurs within the tumor microenvironment. The equations and functional relationships of the model parameters are thus universal across the spectrum of checkpoint inhibitors, whereas the parameter values are cancer type and immunotherapy drug specific as they reflect specific underlying mechanisms of action.

The biophysical processes briefly outlined below and in Fig. 1, and in particular molecule and cell transport physics, effectively set the values of several parameters and thus effectively determine immunotherapy outcome. For example, the physics of transport of immune cells via chemotaxis (χ) and of drug molecule diffusion (*D*_{A}) affect the total immune cell counts (ψ_{k}) in the tumor microenvironment, and the total drug (σ) delivered to and bound (λ) to their respective ligands on T cells (e.g., Fig. 1); the molecular biology and tumor mutational burden (among other factors) are likely to play a major role in the binding, as well as in the values of the pharmacodynamic coefficient λ_{p} representing the specific kill rate effected by immune cells, and in the coupling between immune and cancer cells. This rate is a property of the cancer type, the immune cell types (T cells and dendritic cells), and the immunotherapy drug, and, in principle, can be measured experimentally. It phenomenologically describes the drug-mediated activation of T cells by the dendritic cells, after which the T cells can recognize and kill cancer. Together, these biological properties (and others, see below and the Supplementary Materials) are described mathematically according to the physics of their movement and behavior within the tumor microenvironment, leading to the development of the model presented here.

### Mathematical model

Building on our established models of patient response to chemotherapy (*14*–*21*), we have developed a mathematical model that describes the total tumor burden (ρ) over time and is built upon the key mechanistic biological factors or processes [e.g., concentration of therapeutic T cells (ψ_{k}), intratumoral immunotherapy antibody concentration (σ), cytokine secretion (Λ_{C}), and ratio of immune to tumor cells over time (Λ_{ψ})] and physical factors or processes [e.g., rates of tumor cell proliferation (α) and death (μ), antibody-target binding rates (λ), diffusion of antibodies and cytokines (*D*_{A}, *D*_{C}), specific death rate of cancer cells (λ_{p}), chemotaxis (χ), and mass conservation] that underlie immunotherapy intervention in many or all solid tumors (Fig. 1).

By explicitly combining and mathematically linking the relationships between these biological and physical processes (Fig. 1 and Eqs. 1 to 4), our model quantifies their combined effects (and the feedback processes between them) on the time-dependent change in tumor burden under immunotherapy intervention in a way that is advantageous over descriptions based exclusively on only biological or physical parameters. In words, the four key underlying equations shown in Fig. 1 state that (i) the concentration of viable cells within a tumor over time is a function of the tumor’s intrinsic growth rate reduced by tumor cell kill rate, which occurs due to antibody binding (calculated or estimated; in this case, checkpoint inhibitors binding to their respective ligands) and the time history of antibody uptake and binding within the tumor, and (ii) the time-dependent intratumoral concentration of therapeutic immune cells is a function of chemotaxis-mediated migration into and within the tumor via cytokine signaling and immune cell coupling with tumor cells and (iii and iv) describe the steady-state diffusion balances of antigen-antibody interaction and cytokine concentrations within the tumor microenvironment. Solving these four equations (detailed in the Supplementary Materials) leads to our master equation (Fig. 1 and Eq. 5) that mechanically describes tumor burden over time under immunotherapy intervention based on only three parameters, where α_{0} is the intrinsic tumor growth rate without treatment intervention, and superparameters (e.g., each represents a combination of several biophysical factors) μ and Λ represent key biological and physical processes involved in the immunotherapy mechanism (Fig. 1). μ is the tumor kill rate by immunotherapy (mathematically represented as the product of immune cell fitness, penetration into tumor, and antibody binding), and Λ is the patient’s antitumor immune state, here defined as the coupling of immune cell activity and the tumor cell kill [immunogenicity of an individual tumor, mathematically described as the number of cancer cells that an immune cell can kill (dimensionless) scaled by the ratio of tumor cells to intratumoral immune cells at the time of treatment]. Definitions of these superparameters and other key parameters are provided in Table 1. Together, these superparameters serve as mathematical biomarkers of immunotherapy response. Additional details on model development, derivations, and underlying assumptions are provided in the Supplementary Materials.

At early times subsequent to treatment initiation, i.e., from initiation of treatment up to the first restaging, the short-term solution of Eq. 5 is

Thus, the early stages of treatment depend primarily on the tumor cell kill rate (μ) representing the decrease in the tumor growth rate due to immunotherapy, which may be calculated as the difference between baseline (pretreatment) and long-term (over the course of treatment; calculated retrospectively after time of last follow-up) growth rates (Fig. 2 and Eq. 7). These parameters are, in principle, directly measurable from histopathology and imaging (e.g., CT). All parameters used in the results presented here are listed and defined in Table 1.

### Patients

All studies were conducted in accordance with U.S. Common Rule and with Institutional Review Board Approval at MD Anderson (PA14-0852), including waiver of informed consent. We reviewed all patients that had completed participation in the clinical trials and were eligible for inclusion (calibration cohort, *n* = 28; institutional validation cohort NCT02239900, *n* = 93) by the time of our study; some patients were excluded from analysis (calibration, *n* = 2; institutional validation, *n* = 3) due to unavailable pretreatment CT imaging. Regarding the calibration cohort, we note that data for a total of 58 patients were obtained for the calibration study; however, only 28 were useable, as 17 patients had received nonimmune checkpoint inhibitor immunotherapy, 11 had been concurrently treated with complimentary standard (i.e., nonimmunotherapy) or noncheckpoint inhibitor immunotherapy, and 2 were missing pretreatment measurements needed to quantify α_{0}. Out of the total 121 patients included, for the calibration cohort, 14.3% (4 of 28) were responders (tumor burden reduced at last restaging, i.e., ρ ≤ 1), 2 of whom were pseudo-progressors (showed initial tumor burden increase followed by subsequent reduction in tumor burden; also responders), and 85.7% (24 of 28) were nonresponders (tumor burden increased at last restaging, i.e., ρ > 1), while in the institutional validation cohort, 22.6% (21 of 93) were responders (of these, 6 were pseudo-progressors) and 77.4% (72 of 93) were nonresponders. Patient characteristics are described in tables S1 and S2 for the calibration and institutional validation cohorts, respectively.

### Determining normalized total tumor burden by CT analysis

All patients underwent triple-phase (precontrast, arterial, and portal venous phases) CT scans at baseline. For postcontrast phases, 2.5-mm-thick slices were obtained. Arterial and portal venous phase scanning were initiated with 20- to 25-s and 50- to 60-s delay, respectively. At each restaging, routine abdomen, pelvis, and lung CT scans were done. Lesion measurements were taken on postcontrast CT scans at baseline and at each restaging (restagings ranged from 1 to 12; median, 2). Selection of indexed lesions and follow-up guidelines adhered to standard RECIST 1.1 procedures, and the long and short axes of each indexed lesion (total indexed lesions ranged from 1 to 9) were determined at each follow-up time point (*4*). Axes were averaged geometrically to obtain a reasonable approximation of the mean lesion diameter (*d*), which was then converted to the associated lesion volume (mm^{3}) by approximating lesions as three-dimensional (3D) spheres of radius *d*/2. Any nonindexed lesions were omitted from our calculations. This method of estimation of tumor volume calculated from orthogonal diameter measurements is based on published methodologies, which have been shown to have high correlation with region of interests (ROIs) determined from 3D scans in patients (*22*), and are also standard and accurate in murine xenograft models (*23*, *24*). If lesions were obscured or not seen in images at a restaging time point, long and short axes were estimated via linear interpolation or extrapolation from measured axes in the immediately preceding and following restaging. In all work presented here, the time of immunotherapy initiation is set as *t* = 0, with pretreatment events being *t* < 0 and all events after treatment initiation as *t* ≥ 0. At each time point, we calculated a representative total tumor burden for each patient by summing the volumes of all indexed lesions at each time point divided by the total burden at beginning of treatment. We refer to this normalized quantity as total patient tumor burden (ρ) in this article. Representative time-course data are shown in fig. S1.

### Measuring baseline tumor growth rate (α_{0}), long-term tumor-cell killing rate (μ), and antitumor immune state (Λ) from imaging

Equation S2 was fit numerically to these time-course data using the built-in Mathematica function NonLinearModelFit (*25*) (numerically optimizes a defined function to data) to retrospectively obtain best-fit, “exact” model parameter values Λ and μ for each patient.

This was a two-step process for each patient (Fig. 2); in step 1, lesion data from before treatment initiation (*t* < 0) and treatment initiation (*t* = 0) were interpolated to determine the pretreatment growth kinetic rate α_{0} for each patient assuming exponential growth kinetics before initiation of therapy according to Eq. 8 (see also Eq. 6 and its related considerations). Then, α_{0} was inputted into eq. S2, leaving only two unknowns: Λ and μ, whose values were then obtained in step 2 from the nonlinear fitting of eq. S2 to the patient tumor burden data (ρ) measured from imaging at *t* ≥ 0 (Table 1 and fig. S1, D to F).

### Measurements of model parameters from imaging at first restaging

A patient-specific, accurate estimation of the tumor growth rate after immunotherapy α_{1} (and thus of parameter μ_{1} from Eq. 10) at time of first restaging during the course of treatment was calculated for each patient by fitting the short-term model solution between the measured tumor burden at time of treatment initiation and at the time of first restaging. The exponential tumor growth rate was measured via Eq. 9 (Fig. 2); note that this definition is consistent with Eqs. 6 and 7.

### Categorizing patients into response groups

For each patient, we analyzed the total normalized tumor burden (ρ) at each restaging time point, including from the time of first restaging to the end of treatment. We define response based on the total tumor burden measured at the time of last patient follow-up relative to baseline tumor burden and thus classify responders (ρ ≤ 1) versus nonresponders (ρ > 1).

### Statistics

All statistical analyses were conducted in Excel, GraphPad Prism version 8, and RSWE (*26*). Significance between groups was calculated using Wilcoxon rank sum (two-tailed). Survival curves were constructed using the Kaplan-Meier method; *P* values were reported as log-rank test (Mantel-Cox). Multivariate survival analyses were performed using Cox proportional hazards model. Overall survival under immunotherapy treatment was defined from time of treatment initiation to time of death or last follow-up.

## RESULTS

### Rate of tumor cell killing (μ) and antitumor immune state (Λ) associated with immunotherapy response

Retrospective model analysis was performed on a patient set obtained from various immunotherapy clinical trials treated with checkpoint inhibitor therapy (*n* = 28, see table S1) to verify that model parameters were able to reliably identify patient response based on quantification of the underlying biophysical and mechanistic forces and processes, to determine how model parameters vary with patient response, and to establish quantified thresholds for predicting patient outcome. Because of our concerns about the small size of this sample, and the misbalance between responders (*n* = 4) and nonresponders (*n* = 24), this analysis was then repeated on, and the results were validated in, two additional patient cohorts: one from a clinical trial performed at MD Anderson (*n* = 93; table S2), and the second obtained from the literature (*n* = 124; table S3).

Our mathematical parameters μ (rate of tumor cell killing) and tumor growth rate after treatment (α_{1}) were found to have statistically significant separation (*P* < 0.01) when we compared them between populations that either had reduced tumor volumes (ρ ≤ 1.0; responders) or increased tumor volumes (ρ > 1.0; nonresponders) at last restaging following immunotherapy treatment (Fig. 3, A and B). We note that Λ was not found to separate patient response significantly (*P* = 0.1991) in the calibration set, and it was found to be only minimally significant in both validation cohorts (Figs. 4 and 5).

### Model parameters associated with individual patient response

We investigated whether model parameters α_{∞} (tumor growth rate after treatment), μ (tumor cell killing rate by immune cells), and Λ (antitumor immune state) would correlate with patient response on a per-patient basis. Receiver operator characteristic (ROC) analysis of patient separation between responder and nonresponder groups in the calibration cohort revealed high accuracy between parameter values μ and α_{∞} and patient response (accuracy = [(true positives + true negatives)/total sample size] = 0.893 and 0.964, respectively), while accuracy was reduced (accuracy, 0.750) when parameter Λ was used. Cutoff values that maximize accuracy of predicting response were found to be μ = 0.01 day^{−1}, Λ = 0.144, and α_{∞} = 0.0 (Fig. 3 and fig. S6). The cutoff value for μ is consistent with the average growth kinetics parameter in the calibration set α_{0} = 0.0120 day^{−1} (i.e., it is roughly equal and opposite the intrinsic growth rate, where μ > α_{0} is required for tumor burden reduction). We note that, for both clinical patient cohorts examined, the intrinsic tumor growth rate (α_{0}) was not significantly associated with response (fig. S4).

Measurement of any of these parameters for a given patient, and comparison to the corresponding cutoff value, leads to determination of the probability of that patient being classified as responders or nonresponders, with an accuracy as high as 88% (for α_{1}) and greater than 89% (if μ is use; table S5). The probability of a patient being in a response group may be plotted as percentage of responder or nonresponder patients with a parameter value greater than the value measured in that patient (fig. S6), allowing the probability of response for new patients to be read off the curve; for example, there is ca. 10% probability of values of μ greater than the cutoff in nonresponder patients and ca. 90% in responder patients (calibration cohort; fig. S6A). Combined, these results provide valuable information on the likelihood of response to immunotherapy for the patient in question.

### Model parameters identify unique characteristics of pseudo-progressors

It has been observed following immunotherapy administration that some patients have an initial increase in tumor size, followed by a decrease in tumor size [e.g., ρ(*t _{i}*) > 1.0, then ρ(

*t*) ≤ 1.0;

_{j}*t*

_{0}<

*t*

_{i}<

*t*]. This clinical phenomenon is called pseudo-progression and represents a challenge when deciding whether to continue immunotherapy or not, especially when patients report feeling better after starting immunotherapy but demonstrate enlarging tumors (

_{j}*27*). Now, there is no reliable method to delineate patients with pseudo-progression from patients who are true nonresponders. Using our model, we are able to estimate tumor kill rate at early time (first restaging) after treatment (μ

_{1}; fig. S5) from α

_{1}as

Our analysis found that, for the clinical validation set, the tumor cell kill rate at the first restaging (μ_{1}) (Eq. 10) for pseudo-progressors is consistent with nonresponders, but at a later time, the tumor cell kill rate (μ) is consistent with responders. Thus, the early measured value of the tumor cell kill rate parameter μ_{1} is lower than the long-term value μ that is obtained from the fits (fig. S5D); we refer to this ramping-up of immune response as Δμ = μ – μ_{1}. This behavior was not observed in normal responder or nonresponder patients, which were observed to have unsubstantial Δμ between first restaging and the later retrospectively determined value (fig. S5E). Moreover, our analysis found that Λ (antitumor immune state under immunotherapy) was statistically similar to responder and nonresponder patients in the calibration cohort but was statistically different from nonresponders in the validation cohort (fig. S5I), suggesting in our mechanistic analysis that pseudo-progression is caused by a delayed ramp-up of tumor cell killing by the immune system [this is explored further in the Supplementary Materials; we note that this was not observed in the calibration cohort, likely due to high patient heterogeneity, small sample size, and low number of pseudo-progressors (*n* = 2)].

### Model parameters are quantifiable at early times after start of treatment

We investigated whether patient response could be predicted at the first restaging, using the growth rate α_{1} directly measured from imaging (see Materials and Methods and Eq. 9) as a surrogate for the exact model parameter α_{∞}. We observed ca. 82% accuracy, as compared to the 96% accuracy of the exact parameter (notably, α_{1} and α_{∞} were found to be strongly correlated; Fig. 3E), using the same cutoff for the exact model parameter α_{∞} (fig. S6 and table S5).

### Model validation

To confirm the results obtained from analysis of the calibration dataset, validation was performed on two additional patient cohorts, first using one obtained from an in-house clinical trial using the CTLA-4 inhibitor ipilimumab (NCT02239900, *n* = 93) and then on a second set obtained from the literature, consisting of data from four clinical trials of anti–PD-1 (nivolumab) or anti–PD-L1 (atezolizumab) applied to four different cancer types (renal cell carcinoma, non–small cell carcinoma, melanoma, and urothelial cell carcinoma; *n* = 124; table S3) (*28*–*31*).

Analysis of the first (institutional) validation cohort included long-term estimation (see Materials and Methods; Fig. 2) of α_{∞}, μ, and Λ, and the measurement of tumor growth after treatment at first restaging (α_{1}). In all cases, model parameter values showed significant separation between responder and nonresponder patients (Fig. 4, A to D and F), and the estimate of tumor growth rate at first restaging α_{1} was found to be strongly correlated with the retrospectively determined true long-term growth rate α_{∞} (Fig. 4E). Further, the optimal values to sort patient response identified in the calibration set were found to also perform well in this institutional validation cohort, with high accuracies of 0.774 for μ, 0.731 for Λ, 1.0 for α_{∞}, and 0.849 for α_{1}.

A second validation cohort was obtained from literature-reported data (consisting only of normalized tumor burden measurements at and after time of treatment initiation), with the intention of testing the model’s ability to perform under conditions of minimal data availability. Analysis was performed as described above. Because of the lack of availability of pretreatment tumor volume data for this set, we parametrized the baseline tumor growth rate α_{0} based on the fastest growing 10% of tumors (Supplementary Materials; fig. S3) to estimate the tumor kill rate (μ) for this patient set. This parameterization and underlying assumptions were validated through extensive sensitivity analysis, which demonstrated stability in μ (maximum observed variation in μ was 13.5% across all iterations). Significant separation (*P* < 0.0001) was observed between μ values for responders and nonresponders for the parameter perturbation range examined (e.g., Eq. 7; detailed in the Supplementary Materials). In addition to μ, all the other model parameters (Λ, α_{∞}, and α_{1}) were also found to show significant separation between response types (Fig. 5, A to D), and α_{1} was significantly correlated with α_{∞} (Fig. 5E).

### Model parameter α_{1} predicts patient survival

Using the available retrospective data about patient outcome in the calibration cohort (death or censor), we were able to correlate α_{1} (tumor growth rate at first restaging) with patient outcome through a ROC analysis. The binary classifier for ROC analysis was defined as patient survival (or death) to 50% patient survival time, i.e., 447 days in the calibration cohort, with the individual α_{1} value for each patient used as the discrimination threshold across the entire patient set. The optimal threshold value for prediction of patient survival to 50% survival time was found to be (Fig. 6A) α_{1} = 0.002 day^{−1} (Youden’s J statistic, 0.25; F score, 0.72). Comparison of Kaplan-Meier survival curves between favorable and unfavorable groups was statistically significant to 95% (*P* = 0.0347). Confirmation of the optimal threshold value α_{1} = 0.002 day^{−1} using the institutional validation cohort also demonstrated significant patient separation to 50% survival time (416 days in this validation cohort, Youden’s J statistic = 0.185, F score = 0.693; Fig. 6B). This result was further validated via Cox proportional hazards multivariate analysis; this is detailed in the Supplementary Materials (table S4).

### Univariate and multivariate analyses

Univariate analysis of the calibration cohort revealed a significant association between neutrophil-to-lymphocyte ratio [NLR, which is reported to be associated with increased overall and progression-free survival in melanoma (*12*) and other cancers (*32*)] and overall survival [hazard ratio (HR) = 1.1 (1.02 to 1.24), *P* = 0.016], and multivariate analysis confirmed this correlation [HR = 1.17 (1.03 to 1.3), *P* = 0.01]; this was confirmed in the institutional validation cohort analysis (Table 2). Examination of our early estimate of patient response (α_{1}) in this limited dataset did not reveal a significant correlation with patient survival (univariate analysis: α_{1} HR = 1.95 × 10^{11}, *P* = 0.08; multivariate analysis: α_{1} HR = 5.22 × 10^{13}, *P* = 0.054). However, α_{1} was found to be significantly correlated with patient survival in the validation set (*P* < 0.01), likely due to this set being a larger, more uniform patient cohort [e.g., only two cancer types all treated with the same checkpoint inhibitor, as compared to the high cancer type and drug heterogeneity in calibration cohort (table S1)].

## DISCUSSION

In this paper, we apply a mathematical model (Eq. 5) to three patient cohorts: two clinical cohorts and one literature-derived cohort to predict response to treatment with checkpoint inhibitor immunotherapy. Our model uses standard-of-care CT imaging data as inputs, and with these inputs, the model is able to provide reliable quantification for three key mechanistic values: (i) the tumor cell killing rate (μ), (ii) the antitumor immune state (Λ), and (iii) the posttreatment growth rate α_{∞} on a per-patient basis. The parameters μ and α_{∞} indicate immunotherapy response (Figs. 3, 4, A and C, and 5B). However, it is important to note that a positive tumor kill rate (μ > 0) is not necessarily sufficient for favorable tumor response; instead, tumor kill rate must exceed the intrinsic growth rate (e.g., μ > α_{0;} Eq. 7). This suggests valuable insights for personalized medicine, where an individual patient’s α_{0} may be used to select treatment aggressiveness on a per-patient basis. We also found that tumor growth rate at the time of first restaging (α_{1}) is significantly associated with tumor burden at end of treatment (Figs. 3 to 5E) and predicts patient survival (Fig. 6). When applied to our institutional validation set, which contains a larger and more uniform patient cohort relative to our calibration cohort (tables S1 and S2), we observed significant correlation between patient survival and our measure of tumor growth at first restaging (α_{1}). Moreover, our model parameters were found to have significant correlation with patient response, while baseline blood markers examined (neutrophils, lymphocytes, and NLR) and growth rates did not (figs. S2 and S4). Last, we have observed in some patients that tumor growth rate is faster after treatment; this manifests as negative values of μ and Λ. It is known that immunotherapy causes increased growth rates in some tumors; however, the biological underpinnings of this phenomenon remain elusive (*33*). Unfortunately, this prevents our inclusion of these biophysical processes in our model; however, we are currently conducting new experiments in our laboratory to isolate the processes that may be responsible for this behavior; this is discussed further in the Supplementary Materials. This work not only provides an easily implemented clinical tool but also a framework to understand and study the mechanistic biophysical underpinnings of immunotherapy response by quantitatively describing the constitutive biophysical processes involved.

For patients exhibiting disease progression (nonresponders), we compared this progression to the growth kinetics before initiation of therapy (α_{0}), projecting that growth would have continued exponentially at the same baseline rate in the absence of treatment (Eq. 8). It was observed that, even for patients who had tumors that continued to grow during treatment (nonresponders), there is a meaningful effect of treatment resulting in a reduction of the growth rate with respect to the rate before treatment: α_{∞} < α_{0}. Reduced tumor growth relative to the projected growth reflects the effect of treatment even for progressive disease (fig. S1, D to F). We have shown that the migration of μ values from those that correspond to nonresponders to those of responders over time is a critical feature of pseudo-progressors (fig. S5). This suggests that their immune system is actually experiencing a gradual but definitive “ramp-up” in patients who are clinically described as pseudo-progressors. This counters the current dogma that pseudo-progression is caused by inflammatory cell infiltrates into the tumor (*27*). This influx of immune cells reaches a steady state (hours/days) in a far shorter time scale than what would likely be observed in restaging (weeks to months later). This has biological implications on the very essence of immunology and warrants further investigation.

With the increasing utility of immunotherapy in the treatment of various forms of cancer, there have been fervent efforts to predict treatment outcome for individual patients. While these studies have yielded positive predictors, they often require input data/samples that are not readily available in all clinical settings. Our contribution to this ongoing effort is the development of a predictor that uses clinically available standard-of-care, noninvasive CT scans to predict patient benefit from immunotherapy as early as the first restaging, making it easily integrated into clinical practice. Furthermore, our predictor can be complementary to other biomarkers to increase accuracy. Conceptually, our model describes the change in tumor mass over time after immunotherapy treatment initiation as a balance of tumor proliferation and tumor-killing effects of immunotherapy, thereby yielding valuable insights into the mechanistic underpinnings of immunotherapy response. While predicting treatment outcome from tumor volume change is not novel in and of itself (*34*, *35*), to our knowledge, this is the first report that mathematically encapsulates the multiple biological mechanisms of immunotherapy and retrospectively applies it to a prospective immunotherapy clinical trial. Furthermore, a key feature of our model is that, by representing tumor volume change as a rate (e.g., per time) instead of by absolute volume change, it allows for variability in terms of time of CT scan acquisition that is often the case in the clinical setting and does not require additional laboratory-based testing. We note that we have intentionally used RECIST v1.1 axial measurements to estimate tumor volumes in the work presented here (in our case, long and short orthogonal axes), instead of other available estimates of tumor volume (e.g., contouring). RECIST gives larger volume weights to tumors that develop extreme aspect ratios relative to tumors that remain relatively uniform in shape, a quantification that is representative of the idea that irregularly shaped tumors often exhibit more aggressive growth; we have intentionally included this clever idea into our methodology as well (mathematically, this is a “virtual volume”). Our model, and the parameters it contains, are built on the fundamental physical law–based mechanistic properties of tumor response, while RECIST 1.1 is empirically derived on the basis of meta-analysis of large patient sets and therefore lacks the mechanistic underpinning our model is able to provide. Hence, our model can provide a platform to study the underlying biological and physical causes of the response outcome, and for prediction of response based on measurements of these quantities (offering a valuable improvement over RECIST or its derivatives).

Some limitations of this approach are that (i) our study is a retrospective review of a prospective clinical trial; (ii) the population is a heterogeneous cohort of patients with liver and/or lung metastasis; (iii) variations in CT imaging protocols among different institutions and in tumor boundary identification among radiologists have not been investigated at this stage; (iv) some patients also underwent varying combinations of radiation therapy (for example, at least one tumor in each patient was irradiated in the patient cohort trial NCT02239900; these were removed from our calculations of tumor burden), and the overall “out-of-field” or abscopal effects remains unclear; and (v) a small number of patients were not able to remain on study until first restaging, and we were unable to apply our model to these patients due to lack of response data, or to associate model parameter α_{1} to survival in this subset (this was 6% of patients in NCT02239900; detailed values for other cohorts examined here are unavailable). However, we emphasize that our results demonstrate that the model performs satisfactorily in all patient cohorts to which it is able to be applied, and this is not diminished by this small set of patients. We understand that treatment response is affected by additional factors not included in our model in its current form, for example, drug resistance, hypoxia, angiogenesis, and stroma. These may also be supplemented by other imaging-derived bulk tumor quantities that could be informative to our model, such as tumor texture or morphology; however, we are not completely sure how these parameters may be of consistent value in such a diverse cohort, and hence, these are not explored in this study. Biophysical factors selected for inclusion in the model as presented here were chosen because they represent the most robust factors that are widely supported by the literature, in terms of both magnitude of effects on immunotherapy intervention and broadness of applicability to a large range of solid tumor types. Furthermore, because of the limited biological knowledge of the basic mechanisms of immunotherapy, we cannot account mathematically for parameters that are still unknown to the field at large. Our use of long- and short-axis measurements to estimate 3D tumor volumes is advantageous as it offers the possibility for rapid transition to widespread clinical application of our model; however, detailed contouring is expected to provide more exact tumor volumes and model parameter quantification; the ability of contouring tumor volumes to increase accuracy of our model predictions is currently under investigation in our laboratory, and one focus of an upcoming publication. Nevertheless, our application provides a useful clinical tool in the emerging landscape of immunotherapy. This successful initial model validation will provide a foundation for continued investigation and further iterations of the model to improve its precision and accuracy.

In conclusion, this paper applied a refined mechanistic mathematical model of immunotherapy response to a prospective clinical trial of checkpoint inhibitor therapy, validated with two separate cohorts representing various immunotherapy drugs, drug mechanisms, and cancer types, to predict the outcomes of patients. In doing so, we demonstrated that our model can accurately describe the variable response patterns of patients and discern pseudo-progressors from traditional responders and nonresponders. Notably, the model’s unique signature of (μ) and (Λ) provides testable hypotheses of a distinct biology in this patient population and warrants further investigation. This early predictor of immunotherapy response may allow for critical adjustments to individual patient treatments and is one step closer to the promise of personalized medicine. Future investigations will focus on expanding the scope of this model to other clinical datasets and further refinement of the complex interactions between the tumor and its milieu.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/18/eaay6298/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**We gratefully acknowledge support from the Andrew Sabin Family Fellowship, Center for Radiation Oncology Research, the Sheikh Ahmed Center for Pancreatic Cancer Research, institutional funds from The University of Texas MD Anderson Cancer Center, GE Healthcare, Philips Healthcare, and Cancer Center Support (Core) Grant CA016672 from the National Cancer Institute (NIH) to MD Anderson. E.J.K. was also supported by Project Purple, NIH (U54CA210181-01, U01CA200468, and U01CA196403), and the Pancreatic Cancer Action Network (16-65-SING). Z.W. and V.C. were also supported by NSF grant DMS-1716737 and NIH grants 1U01CA196403, 1U01CA213759, 1R01CA226537, 1R01CA222007, and U54CA210181. V.C. was also supported by the University of Texas System STARS Award. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Author contributions:**V.C. and E.J.K. supervised the research. V.C., D.S.H., and J.W. conceived and designed the approach. V.C. developed the model. J.D.B., Z.W., and V.C. implemented the model and performed model analysis. J.D.B., D.E., C.X.W., Z.W., S.-H.C., N.F.E., R.P., W.A., D.S.H., J.W., E.J.K., and V.C. analyzed and interpreted the data. D.S.H., J.W., and E.J.K. supervised the clinical data collection. D.E. and J.W. collected the clinical data. J.D.B., D.E., C.X.W., Z.W., S.-H.C., N.F.E., R.P., W.A., D.S.H., J.W., E.J.K., and V.C. wrote or edited the paper.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).