It remains unclear how localized radiotherapy for cancer metastases can occasionally elicit a systemic antitumor effect, known as the abscopal effect, but historically, it has been speculated to reflect the generation of a host immunotherapeutic response. The ability to purposefully and reliably induce abscopal effects in metastatic tumors could meet many unmet clinical needs. Here, we describe a mathematical model that incorporates physiologic information about T-cell trafficking to estimate the distribution of focal therapy–activated T cells between metastatic lesions. We integrated a dynamic model of tumor–immune interactions with systemic T-cell trafficking patterns to simulate the development of metastases. In virtual case studies, we found that the dissemination of activated T cells among multiple metastatic sites is complex and not intuitively predictable. Furthermore, we show that not all metastatic sites participate in systemic immune surveillance equally, and therefore the success in triggering the abscopal effect depends, at least in part, on which metastatic site is selected for localized therapy. Moreover, simulations revealed that seeding new metastatic sites may accelerate the growth of the primary tumor, because T-cell responses are partially diverted to the developing metastases, but the removal of the primary tumor can also favor the rapid growth of preexisting metastatic lesions. Collectively, our work provides the framework to prospectively identify anatomically defined focal therapy targets that are most likely to trigger an immune-mediated abscopal response and therefore may inform personalized treatment strategies in patients with metastatic disease. Cancer Res; 76(5); 1009–18. ©2016 AACR.

Major Findings

Using a mathematical model, we show that the distribution of activated T cells varies significantly between different metastatic sites and depends on the immune system activation site. Therefore, which metastases are chosen for localized therapy contributes to the potential for inducing systemic metastatic regression. We present a framework that could inform about patient-specific treatment targets, single or multiple in metastatic patients, thereby supporting personalized clinical decision making. The developed model offers additional insights into the dynamics of systemic disease, including metastases-mediated acceleration of primary tumor growth and rapid metastatic progression after surgical intervention.

Quick Guide to Equations and Assumptions

A general metastatic cancer setting with N distinct tumors located in different organs, such as lung, liver, breast, or kidney is considered. We advance a minimally parameterized model of immune system interactions within a single tumor site that was successfully compared with experimental data (1). The cancer population at each metastatic site, Ci(t) for i = 1,…,N, is assumed to follow logistic growth with site-dependent carrying capacity, Ki, and growth rate, ri. These populations are modulated by the predation of immunocompetent effector cells, Ei(t). The equation governing growth of each metastatic site is

formula

A detailed description of all parameters with specific values can be found in Table 1. To limit model complexity and in the absence of evidence to the contrary, we assume the same parameter values for each metastatic site, with the exception of possible differences in growth rates and carrying capacities.

Effector cells are recruited endogenously as well as in response to tumor burden and may undergo spontaneous decay and exhaust from interactions with cancer cells. Cells that are recruited in response to the tumor are mainly cytotoxic T cells, a key component of the adaptive immune system. We extend Kuznetsov's local model (1) to metastatic disease such that the equation governing the population of tumor-infiltrating effector cells at each of the metastatic sites is

formula

where functions ωij(C1(t),…,CN(t)) describe the probabilities that a cytotoxic T-cell activated at site j will infiltrate the ith tumor site. The initial delay in establishing a robust immune response is incorporated in the time-dependent parameter f (compare Table 1) for each metastatic site. Probabilities ωij depend on the anatomic distribution of metastatic sites and current tumor sizes.

After being activated in the tumor-draining lymph node, T cells travel through the lymphatic system, enter the blood circulation system, and travel in cycles through the network of arteries, capillaries, and veins. To quantify T-cell movement through the blood system, we distinguish four major coarse-grain flow compartments (COMP): pulmonary circulation (LU); gastrointestinal tract and spleen (GIS); liver (LI); and all other organs in the systemic circulation (SO), such as breast or kidney. The distinct consideration of LU, GIS, and LI is required, as venous blood from GIS flows through LI via the hepatic portal vein, before being reoxygenated in LU (Fig. 1A). T-cell trafficking between these compartments is based on the anatomy, with rates being assigned as physiologic values in the form of blood flow fractions (BFF; percentage of cardiac output reaching the compartment). We denote HCOMP as the probability that the T cell will extravasate at one of the metastatic sites after entering a given compartment, i.e., |$H_{COMP} = \sum\nolimits_{i = 1}^{N_{COMP}} {P_{ij}}$|⁠, where Pij is the probability that a T cell activated at site j will infiltrate the ith tumor site after entering a given compartment (Fig. 1B). We calculate the probability of T-cell absorption by a given compartment (WCOMP) when starting from the LU compartment using a Markov chain approach

formula

where Δ is a normalizing constant such that the sum of absorption probabilities over all four compartments is equal to one (see Supplementary Material for details). Without evidence of T-cell decay in the circulating compartment of the blood system, each T cell will eventually be absorbed in one of the compartments, with the average number of systemic circulation cycles before absorption

formula

Probability Pij is equal to the probability that the T cell will flow through the tumor site multiplied by hij, the probability of extravasation from blood into the tissue. We approximate the former using the relative blood flow through a specific tumor–bearing organ multiplied by the fraction of organ populated by the tumor, and, thus, the equation describing Pij is

formula

Experimental studies show that T cells are programmed during the activation process in the lymph node to express homing molecules that favor trafficking back to the site of the activation (2, 3). Thus, we assume hij = ha if the T cell was activated in the given organ (organi = organj), and hij = hn otherwise (1 ≥ ha > hn). The probabilistic model defined above allows the calculation of |$\omega _{ij} = W_{COMP_i} P_{ij}/H_{COMP_i}$| in the effector cell dynamics equation (Eq. B) for different anatomic distributions of metastatic disease (see Supplementary Material for details).

Transformed cancer cells are confronted by both innate and adaptive immune surveillance, and it is believed that tumors have evolved mechanisms to evade the immune system even before they become clinically apparent (4). The notion of increasing immune system efficacy by therapeutic intervention to systemically eradicate cancer cells has long been a vision of oncologists and cancer researchers. A particularly exciting development, hailed by the editors of Science as the scientific breakthrough of 2013 (5), is that novel immunotherapeutic strategies show remarkable responses in some patients, especially if combined with common cytotoxic agents. Radiotherapy and chemotherapeutic agents have been shown to substantially enhance tumor-specific immune responses in well-established tumors (4, 6–11). The synergy between radiotherapy and immunotherapy stems from radiation-induced effects including (i) immunogenic cell death that locally exposes a wealth of tumor antigens and (ii) release of danger signals, which act as endogenous immune adjuvants to stimulate the activation of dendritic cells (DC; Fig. 1C; ref. 12).

Figure 1.

Activated T-cell trafficking. A, activated cytotoxic T cells (CTL) enter the blood system via the great veins, flow through the pulmonary circulation, and then continue into systemic circulation. Venous blood from the gastrointestinal tract and spleen goes to the liver through the hepatic portal vein. B, cytotoxic T cells can flow through each compartment without reaching any of the tumor sites (Mi, circles) or flow through one of them. At the metastatic site, a T cell has a certain probability to?extravasate, ha if it has been activated by that site and hn otherwise. We denote by Hcompartment the overall probability that a T cell will extravasate in the given compartment. C, radiotherapy (RT) triggers immunogenic cell death that activates DCs that subsequently travel to the lymph node. DCs transform naïve T cells into cytotoxic T cells, which, in case of nonmetastatic disease, traffic through the blood system back to the tumor site for cancer cell eradication.

Figure 1.

Activated T-cell trafficking. A, activated cytotoxic T cells (CTL) enter the blood system via the great veins, flow through the pulmonary circulation, and then continue into systemic circulation. Venous blood from the gastrointestinal tract and spleen goes to the liver through the hepatic portal vein. B, cytotoxic T cells can flow through each compartment without reaching any of the tumor sites (Mi, circles) or flow through one of them. At the metastatic site, a T cell has a certain probability to?extravasate, ha if it has been activated by that site and hn otherwise. We denote by Hcompartment the overall probability that a T cell will extravasate in the given compartment. C, radiotherapy (RT) triggers immunogenic cell death that activates DCs that subsequently travel to the lymph node. DCs transform naïve T cells into cytotoxic T cells, which, in case of nonmetastatic disease, traffic through the blood system back to the tumor site for cancer cell eradication.

Close modal

Most fascinating is the observation that the stimulation of the immune system by localized radiotherapy may modulate systemic regression of metastatic nodules, known as the radiation-induced abscopal effect (13). Such abscopal responses triggered by focal radiotherapy have been reported for multiple cancer types including renal cell carcinomas (14) and papillary adenocarcinomas (15), yet due to the high volume of patients with metastatic disease being treated with radiotherapy, many still consider these reports to be nothing more than anecdotal (16). However, with a combination of both irradiation and immunotherapy, abscopal responses can be triggered more reliably (12). A strong systemic response against squamous cell carcinoma in a murine model was observed when DCs were administered intratumorally following local irradiation (17).

The possibility of rationally inducing abscopal effects using radiotherapy and immunotherapy has the flavor of the long-sought “magic bullet” of cancer therapy, yet generating the synergy required to provide local control and also induce the abscopal effect is difficult. A myriad of factors contribute to this challenge, including (i) heterogeneity of mechanisms a tumor can use to resist immune attack; (ii) heterogeneity of the capacity for an anticancer response by different components of the immune system; and (iii) anatomic- and time-dependent treatment effects. Numerous mathematical models have been developed to describe tumor–immune interactions at different phases of tumor progression (18–20) or to look at different pathways that can be exploited for immunotherapy (21–25). However, to our knowledge, models of metastatic cancer and both local and systemic immune system interactions currently do not exist. Thus, no prominent inroads have been made to decipher the contribution of local therapy to systemic disease modulation in metastatic patients. Herein, we propose a tumor–immune system interaction modeling framework that incorporates a mathematical model of activated T cells trafficking between metastatic sites. In addition to numerous biologic bottlenecks, such as heterogeneities in immune infiltration, immune repertoire, and antigen presentation, it is conceivable that an abscopal response can only be achieved if locally activated T cells are systemically distributed among the metastatic sites, and in numbers sufficient to tip immune surveillance back in favor of tumor eradiation at each of these sites (26). We show that trafficking and distribution of activated T cells strongly depends on the anatomic distribution of metastatic sites, physiologic blood flow fractions to tumor bearing organs, tumor burden in each metastatic tissue, and the strength of cues that impact immune-cell homing (2, 3). It follows that different metastatic sites may have varying potential to trigger a systemic response following focal immune-activating therapy. On the basis of the T-cell homing distribution, we determine optimum treatment targets in a virtual patient cohort. We integrate an extended tumor–immune interaction model (1) into the framework to simulate growth of each metastasis. This provides nonintuitive insights into how the seeding of a new metastatic site can promote the growth of a primary tumor and offers an elegant explanation of the observation of rapid progression of metastatic disease after surgical removal of a primary tumor (27–30).

The different components of the mathematical framework discussed in the Quick Guide to Equations and Assumptions are calculated or solved numerically in Matlab (www.mathworks.com). A detailed description of the computational methods can be found in Supplementary Material.

Model parameter estimation

The physiologic blood flow fraction (BFF) to the spleen is estimated to be 3% (31). The gastrointestinal tract consists mainly of stomach and esophagus (BFF = 1%), intestines (BFF = 16%), and pancreas (BFF = 1%), which together with the BFF to the spleen, yields BFFGIS = 21%. It is estimated that the internal mammary arteries with an average blood flow of 59.9 mL/minute (32) provide approximately 60% of breast blood supply (33). With an average cardiac output of 300 L/hour (34), we estimate BFFbreast = 2%. BFF to the kidney is estimated to be 8.5% (31).

In a study performed by Mikhak and colleagues (2), ovalbumin (OVA)-specific T cells isolated from transgenic mice were activated in vitro with DCs isolated from different mouse tissues including lung, thoracic lymph nodes, and skin. Populations of activated T cells were then injected into naïve mice challenged with aerosolized OVA. At 24 hours after the last challenge, mice were sacrificed and lung tissue harvested to measure OVA-specific T-cell counts. T cells that were activated with lung DCs were present in the lung in numbers three times greater than T cells activated with DCs from different organs. Thus, we estimate the extravasation probability of T cell into the tissue of activation (ha) to be approximately three times larger than extravasation into nonactivation tissue sites (hn; hn/ha = 1/3). Values of other parameters associated with T-cell trafficking, that is, those that are necessary to calculate values of ωij(C1,…,CN), were taken from the literature and are summarized in Table 2.

Table 1.

Parameter values for the differential part of the model, equations A and B, taken from ref.?1 

ParameterValueUnitDescription
r 0.19 (for tumor located in breast) 0.38 (for tumor located in lung) 1/day Maximal tumor growth rate 
K 531.9 × 106 cells Carrying capacity 
a 0.14 × 10−6 1/day/cells T cell–cancer cell interactions constant 
p 0.998 Nondimensional Probability that during the interaction between T cell and cancer cell, the latter will be killed 
λ 0.59 1/day Effector cells decay rate 
E* 0.3 × 106 Cells Physiologic level of effector cells 
f 0.53 if time since metastatic site creation > 28 days and 0 otherwise 1/day/cells Magnitude of immune system stimulation by the presence of cancer cells 
g 0.16 × 106 Cells Immune stimulation damping coefficient 
ParameterValueUnitDescription
r 0.19 (for tumor located in breast) 0.38 (for tumor located in lung) 1/day Maximal tumor growth rate 
K 531.9 × 106 cells Carrying capacity 
a 0.14 × 10−6 1/day/cells T cell–cancer cell interactions constant 
p 0.998 Nondimensional Probability that during the interaction between T cell and cancer cell, the latter will be killed 
λ 0.59 1/day Effector cells decay rate 
E* 0.3 × 106 Cells Physiologic level of effector cells 
f 0.53 if time since metastatic site creation > 28 days and 0 otherwise 1/day/cells Magnitude of immune system stimulation by the presence of cancer cells 
g 0.16 × 106 Cells Immune stimulation damping coefficient 

NOTE: C(t) measures cell number, but in the T-cell trafficking part of the model, its value is used in terms of tumor volume in milliliters according to the formula V(C(t)) = C(t)×(4/3 π r3)×10−12 mL, where r is the diameter of the cell assumed to be 10 μm.

Table 2.

Parameter values for the T-cell trafficking part of the model

ParameterValueDescriptionReference
BFFLI 6.5% Blood flow fraction to liver (31) 
BFFGIS 21% Blood flow fraction to gastrointestinal tract Estimated 
BFFSO 72.5% Blood flow fraction to SO compartment By definition = 100% − BFFLIBFFGIS 
BFFbreast 2% Blood flow fraction to the breast Estimated 
BFFkidney 8.5% Blood flow fraction to the kidney (31) 
Vliver 1,493 mL Average liver volume (48) 
Vbreast 500 mL Average breast volume (34) 
Vlungs 3,679 mL Average lung volume (49) 
Vkidney 249 mL Average kidney volume (50) 
ParameterValueDescriptionReference
BFFLI 6.5% Blood flow fraction to liver (31) 
BFFGIS 21% Blood flow fraction to gastrointestinal tract Estimated 
BFFSO 72.5% Blood flow fraction to SO compartment By definition = 100% − BFFLIBFFGIS 
BFFbreast 2% Blood flow fraction to the breast Estimated 
BFFkidney 8.5% Blood flow fraction to the kidney (31) 
Vliver 1,493 mL Average liver volume (48) 
Vbreast 500 mL Average breast volume (34) 
Vlungs 3,679 mL Average lung volume (49) 
Vkidney 249 mL Average kidney volume (50) 

The basic structure of Eqs. A and B are adapted from the Kuznetsov model (1), from which we adopt the parameter values that were estimated therein (Table 1).

Virtual case studies and immunogenicity index quantification

We create a cohort of 40 virtual metastatic patients with arbitrary combinations of breast, liver, kidney, and lung metastases of random sizes between 50 and 300 mL. Given static information about patient-specific metastatic anatomy prior to therapy, we investigate the homing distribution of T cells activated in response to localized therapy, such as radiotherapy applied to one of the tumors. To this extent, we calculate values of ωij for i = 1,…,N when activation occurred at the jth site (see Eq. B) for all sites of activation (without simulating the full temporal evolution model). We determine T-cell dissemination quality for activation at a given site by calculating the ratio of the entropy of the established homing distribution to the entropy of uniform distribution

formula

where j is the site of activation, and ωij is the probability that a T cell activated at site j will infiltrate the tumor at site i. Higher values of Sj indicate T-cell homing distributions closer to uniform. The number of activated T cells depends on the number of tumor cells undergoing immunogenic cell death after local therapy, which in turn is dependent on tumor volume. Thus, we define the immunogenicity index of the ith tumor, Ii, by a combination of its homing distribution entropy with its size relative to the largest metastatic site.

formula

The immunogenicity index can be considered to represent the systemic reach of activated T cells, which contributes to the probability of triggering a systemic abscopal effect. The maximum immunogenicity index of Ix = 1 can only be achieved theoretically, that is, if the largest tumor Vx is treated and the corresponding distribution of activated T cells is uniform between the metastatic sites (Sx = 1). For expected nonuniform distributions of activated T cells (Si < 1), however, tumor volumes alone are insufficient to predict immunogenicity indices.

Simulation of metastatic disease progression

In simulations of tumor progression (solving the full model described by Eqs. A and B), we arbitrarily initiate a breast tumor with 0.5 × 106 cancer cells. After 200 days of simulated growth of this primary tumor, we simulate the onset of a metastasis by initiating a population of 0.5 × 106 cancer cells in the lung, following which, we simulate the growth of both tumor sites for an additional 400 days. For the purpose of illustration, we assume the growth rate of the metastatic site in the lung to be twice as fast as the growth rate of the breast tumor. We assume extravasation probabilities ha = 0.6 and hn = 0.2.

Simulation of primary tumor removal

After simulating growth of the primary breast tumor and a lung metastasis, we mimic complete surgical removal of the primary tumor by instantaneously removing both cancer cell and T-cell populations present in the breast. Parameters describing intrinsic lung metastasis dynamics remain unchanged.

Variation in immunogenicity between metastatic sites

Intuitively, if T cells are unable to extravasate at any site other than the tissue of activation, i.e., hn/ha = 0, then no systemic response is possible. On the other hand, if extravasation occurs at all sites with equal strength (i.e., hn/ha = 1), the systemic T-cell distribution is independent of the tissue of activation. Thus, it is important to investigate the intermediate values of the hn/ha ratio. We begin with a detailed analysis of a virtual case study of breast (113 mL), liver (220 mL), and lung (270 mL) tumors. Figure 2 shows the homing distribution of T cells for the different sites of immune system activation with ha = 0.6 and hn/ha = 1/3. It can be seen that the distribution of activated T cells depends on the site of activation. The probability that an activated T cell will home to the breast metastasis is only larger than 10% if the immune system was activated at that site (or possibly another site located also in the breast). The highest value of normalized entropy of T-cell homing distribution is obtained for activation in breast (Sbreast = 0.85). Activation in lung or liver attracts about 98% of all activated T cells back to these two organs (Slung = 0.43, Sliver = 0.66). Figure 3A–C shows the distribution of T cells for immune system activation in breast, liver, and lung, respectively, for different values of parameters ha∈[0,1] and hn/ha∈[0,1]. Simulations reveal that the actual extravasation probability at the tissue of activation, ha, has a negligible influence on the final systemic T-cell distribution but has a large impact on the temporal component of T-cell trafficking, i.e., the average number of circulatory cycles performed before extravasation. Experimental derivation of the extravasation probability becomes imperative for the design of optimum local therapy fractionation protocols for immune response induction. Final T-cell distribution is mostly determined by the ratio of extravasation probabilities at nonactivation sites versus activation sites, hn/ha (Fig. 3A–C). For intermediate values of hn/ha, including estimated hn/ha = 1/3, the T-cell dissemination patterns vary nonintuitively across activation sites.

Figure 2.

Model-predicted homing distribution for different activation sites. Results for a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases with ha = 0.6, hn/ha = 1/3. Activation in breast yields the most uniform distribution (Sbreast = 0.85 vs. Sliver = 0.66 vs. Slung = 0.43). Values of other parameters used in calculations are reported in Table 2.

Figure 2.

Model-predicted homing distribution for different activation sites. Results for a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases with ha = 0.6, hn/ha = 1/3. Activation in breast yields the most uniform distribution (Sbreast = 0.85 vs. Sliver = 0.66 vs. Slung = 0.43). Values of other parameters used in calculations are reported in Table 2.

Close modal
Figure 3.

T-cell homing distribution for different activation sites and extravasation probabilities. Model-predicted homing distributions between metastatic sites in a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases for different sites of activation: breast (A), liver (B), lung (C), and for different values of extravasation probabilities hn and ha. Rectangles correspond to the narrow ranges around the value of hn/ha estimated from the literature. D, average number of transitions between model compartments before extravasation at one of the metastatic sites for different sites of activation and different extravasation probabilities hn and ha. Calculations were performed using parameters reported in Table 2.

Figure 3.

T-cell homing distribution for different activation sites and extravasation probabilities. Model-predicted homing distributions between metastatic sites in a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases for different sites of activation: breast (A), liver (B), lung (C), and for different values of extravasation probabilities hn and ha. Rectangles correspond to the narrow ranges around the value of hn/ha estimated from the literature. D, average number of transitions between model compartments before extravasation at one of the metastatic sites for different sites of activation and different extravasation probabilities hn and ha. Calculations were performed using parameters reported in Table 2.

Close modal

Given the complexity of the arterial tree, trafficking through a specific site is a relatively rare event (evident from the BFFs reported in Table 2). In addition, ha and ratio hn/ha strongly impact the number of cycles that the T cells will remain in transit in the circulatory system before extravasating at a tumor-harboring site (Fig. 3D). It follows from Eq. C that the average number of cycles increases with smaller values of ha. In the case of immune system activation in the breast, if ha = 1 and hn/ha = 1, T cells will perform an average of approximately 8 circulatory transit cycles before extravasation. With a blood recirculation time of 10 to 25 seconds (35), this translates to around 2 minutes. For lower values of ha = 0.05 and hn/ha = 0, our model suggests that a T cell could complete approximately 4,500 circulatory cycles (∼32 hours) before extravasation.

Comparing the entropy between T-cell homing distributions (Eq. D) shows that, in this virtual patient example, activated T cells distribute most uniformly after treating the breast tumor metastasis with focal radiotherapy as the immune system activating agent (Fig. 4A), regardless of extravasation probabilities ha and hn. Targeting the lung tumor, however, yields a heavily skewed distribution of T cells toward the lung, given the large BFF to this site. Integration of tumor volumes to calculate the immunogenicity index Eq. E, the metric we propose to support clinical decision making, reveals that the liver metastasis has the highest immunogenicity index for hn/ha values between 0.15 and 0.6, including the estimated value of hn/ha = 0.33 (Iliver = 0.49, Ibreast = 0.38, Ilung = 0.34; Fig. 4B).

Figure 4.

T-cell trafficking and extravasation for different activation sites in a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases. A, normalized entropy values (Eq. D). B, immunogenicity indices (Eq. E) for different extravasation probabilities hn and ha. C, analysis of 40 virtual case studies of possible metastatic tumors in lung, liver, kidney, and the breast. Circles denote existence of the metastatic site, and radii correspond to tumor sizes. Black background identifies the tumor with the highest immunogenicity index for hn/ha = 1/3 and ha = 0.6. Calculations were performed using parameters reported in Table 2.

Figure 4.

T-cell trafficking and extravasation for different activation sites in a virtual patient with breast (113 mL), liver (220 mL), and lung (270 mL) metastases. A, normalized entropy values (Eq. D). B, immunogenicity indices (Eq. E) for different extravasation probabilities hn and ha. C, analysis of 40 virtual case studies of possible metastatic tumors in lung, liver, kidney, and the breast. Circles denote existence of the metastatic site, and radii correspond to tumor sizes. Black background identifies the tumor with the highest immunogenicity index for hn/ha = 1/3 and ha = 0.6. Calculations were performed using parameters reported in Table 2.

Close modal

We additionally calculated which site has the highest immunogenicity index in a virtual patient cohort for the estimated values of hn/ha = 1/3 and ha = 0.6 (Fig. 4C). Calculations suggest that largest tumor volume or a specific combination of metastases is not necessarily always an indicator of maximum system response or optimal treatment target for an individual patient. A much smaller breast tumor (53 mL) has a higher immunogenicity index than the much larger kidney tumor (254 mL) in virtual patient 5 (Fig 4C). Virtual patients 20 and 21 have the same combination of virtual sites of metastatic tumors, but the tumor with the highest immunogenicity index is in different organs. Therefore, the identification of treatment targets with the highest likelihood of inducing a systemic abscopal response is highly patient-specific and not predictable from the average response of historical cases.

Systemic interdependence of metastatic sites through activated T-cell trafficking

Simulation of primary tumor growth for 200 days shows that a primary tumor can be modulated by a cytotoxic immune surveillance, keeping the tumor in a dormant state in equilibrium with the immune system (26) and below the imposed carrying capacity (Fig. 5A). Seeding of a metastatic site disrupts the dynamic equilibrium between T cells and cancer cells in the primary site. Before the immune response against the metastasis is established locally, some of the T cells generated in the primary breast tumor travel and extravasate at the metastatic site, allowing for transient immune escape and progression of the primary tumor (Fig. 5B). After T cells activated by the lung metastasis begin trafficking in the circulatory system, the total number of T cells increases, and an increase in T cells surveilling the primary breast cancer can be observed. However, because of the attained equilibrium of homing distributions (Fig. 5C and D), the now lower proportion of T cells penetrating the breast allows the tumor to attain dormancy with cell numbers larger than during the time period prior to metastasis formation. Interestingly, the lung metastasis is also kept at a dormant state despite its larger growth rate, mostly due to T cells that are generated in the breast and trafficking to the lung. It is worth pointing out that due to its assigned large proliferation rate, a single tumor site in the lung without a preexisting breast tumor would escape immune surveillance and grow close to its imposed carrying capacity (Supplementary Fig. S2). Thus, surgical removal of the breast tumor results in removal of locally infiltrated T cells as well as prevention of future T cells activated in the breast. This leads to the escape of the lung metastasis from immune surveillance due to the rapid decrease of the effector-to-cancer cell ratio in the lung (Fig. 5A).

Figure 5.

Partial removal of T-cell surveillance and metastatic progression after surgery. A, model-predicted growth of a primary breast tumor before and after the onset of a lung metastasis at t = 200 days, and surgical removal of the primary breast tumor and its local effector cells at day 600. B, corresponding T-cell/tumor cell ratios for the tumors shown in A. C and D, homing distributions of T cells that are activated in the breast tumor and lung metastasis, respectively. Simulations were performed with ha = 0.6, hn = 0.2 and parameters reported in Tables 1 and 2. rbreast = 0.19/day, rlung = 0.38/day.

Figure 5.

Partial removal of T-cell surveillance and metastatic progression after surgery. A, model-predicted growth of a primary breast tumor before and after the onset of a lung metastasis at t = 200 days, and surgical removal of the primary breast tumor and its local effector cells at day 600. B, corresponding T-cell/tumor cell ratios for the tumors shown in A. C and D, homing distributions of T cells that are activated in the breast tumor and lung metastasis, respectively. Simulations were performed with ha = 0.6, hn = 0.2 and parameters reported in Tables 1 and 2. rbreast = 0.19/day, rlung = 0.38/day.

Close modal

Progression from localized to metastatic disease sharply diminishes patient prognosis. For example, in 10% to 15% of breast cancer patients, the second most common cancer diagnosed worldwide, metastasis will develop within 3 years of diagnosis (36), decreasing the 5-year survival from >80% to around 26% (37). When targeting both local and metastatic disease, immunotherapy has been shown to synergize with local radiation (6–12), which is currently being explored in more than 10 active clinical trials (12). The biologic and physical principles underlying this synergy are yet to be completely understood. The complexity arises from the dynamic interplay between the immune system and the tumor, which is further confounded when examining local as well as systemic dynamics. The increasing number of case reports of abscopal effects in metastatic tumors distant to the areas targeted by radiation (14, 15) encourages further investigation into the possibility of intentionally triggering such responses. We hypothesized that different metastatic sites in individual patients may have varying potential to induce a systemic response, and that mathematical modeling could help to identify the most promising treatment targets prior to intervention on a per patient basis.

A myriad of complex biologic cascades are without doubt involved in the development of abscopal effects. Although T-cell trafficking may only be one contributing mechanism to an abscopal response, it follows physiologic blood flow and hence has an element of predictability. Identifying local treatment targets with the highest systemic reach of activated T cells promises to increase the likelihood of purposely inducing abscopal effects. We propose a systemic tumor–immune interaction framework, which accounts for activated T-cell trafficking through the host circulatory system. Patient-specific distribution of metastatic sites can be acquired from routine radiology scans including CT, MRI, and PET/CT. When combined with physiologic information about T-cell trafficking, the model estimates the distribution of focal therapy–activated T cells for each metastatic site. Using a virtual patient cohort, we showed that the activated T-cell distribution is dependent on (i) the anatomic distribution of metastatic sites, (ii) the tumor volume of each metastasis, and (iii) the site of activation. In the clinic, once patient-specific anatomic metastatic distribution and individual tumor volumes are obtained from radiology, the presented framework can calculate proposed immunogenicity indices for each of the metastatic sites. Those values may support clinical decision making as to which metastatic sites serve as the most promising local treatment targets to induce systemic abscopal responses, which can be further enhanced through various immunotherapeutic approaches (12).

Radiotherapy may be the most intuitive local treatment choice, as the abscopal effect has been observed clinically after focal radiation both with and without immunotherapy (14, 15). The framework may be equally applicable to other potentially immune-activating therapies such as radiofrequency ablation (high frequency current to heat and coagulate tissues), cryoablation (extreme cold to cause tumor destruction), or photodynamic therapy (drug activated by light of suitable wavelength; ref. 38). Before translation into a prospective clinical trial to validate clinical decision support, however, the model needs to be validated with retrospective clinical data, and estimated parameters need to be calibrated with the results of experimental studies. Of utmost importance is the tracking of labeled T cells and measurement of their distribution among different tissues and metastatic sites, as well as quantification of the homing rate. This information would allow the determination of optimum time intervals among radiotherapy fractions to facilitate a strong immune response. In future work, the systemic organs compartment may be further divided into anatomically distinct and important areas such as bone, skin, muscle, and central nervous system. Bone in particular is often treated with focal radiation for definitive or palliative control (39), and thus the explicit simulation of T-cell trafficking to the bone is of clinical importance.

The presented framework reproduced several established key features of local tumor–immune interactions including long-term tumor dormancy (19, 26) but additionally revealed a previously unappreciated interdependence of metastatic sites through activated T-cell trafficking. First, the onset of a new metastatic site can boost the growth of the primary tumor by diverting T cells to the distant metastasis and therefore facilitating a transient escape from immune-induced dormancy. This might offer yet another angle to the Norton hypothesis that tumors may not metastasize because they are large, but may be large because they are self-metastatic (40). Second, immune system activation and activated T-cell trafficking from the primary tumor may impede the growth of a metastatic site by mounting an immune response in the form of movement of T cells from the primary to the metastasis, a process that is degraded at the time of resection of the primary. This augments an understanding of the modulation of distant tumors from the primary, a phenomena recognized more than 100 years ago and termed "concomitant immunity" (41). Third, surgical removal of the primary tumor may trigger progression of existing metastatic sites, which has been seen in the clinic (27–30) as well as in animal experiments (42, 43) and mathematical models (44, 45). With 93% of breast cancer (early stage I and II), 98% of colon cancer (stage I, II, and III), and 71% of non–small cell lung cancer (early stage I and II) patients undergoing surgery as a part of their treatment (46), systemic perturbation by the presumed local treatment needs to be explored further.

Herein, we present the first attempt to quantify how local tumor–immune interactions may propagate systemically, which may lead to the prediction of patient-specific treatment targets to trigger abscopal effects. For illustration purposes, we utilized the established Kuznetsov model of the interaction of local tumor with the immune system, starting with parameter values developed using data obtained from a mouse lymphoma model (1). This module of the framework, however, can be interchanged for other models of tumor growth and immune modulation. In future work, either the Kuznetsov model or its replacement will need to be calibrated with organ-specific tumor growth kinetics and immune surveillance to fully integrate local and systemic dynamics and confidently support personalized decision making in the clinic. Furthermore, the model of tumor–immune interactions for each local site may be expanded to include other cell types, as well as tissue specific parameterizations of carrying capacity, growth rates, and immune interaction parameters. The blood flow associated trafficking model can be also utilized to explore the metastatic dissemination patterns, which is currently being investigated by others (47).

M. Fishman reports receiving other commercial research support from BMS. H. Enderling has ownership interest (including patents) in a patent application that has been filed and is in the international PCT stage (serial #: PCT/US2015/024278). No potential conflicts of interest were disclosed by the other authors.

Conception and design: J. Poleszczuk, K.A. Luddy, S. Prokopiou, M. Robertson-Tessi, E.G. Moros, M. Fishman, J.Y. Djeu, S.E. Finkelstein, H. Enderling

Development of methodology: J. Poleszczuk, K.A. Luddy, E.G. Moros, M. Fishman

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.E. Finkelstein

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Poleszczuk, J.Y. Djeu, S.E. Finkelstein

Writing, review, and/or revision of the manuscript: J. Poleszczuk, K.A. Luddy, S. Prokopiou, M. Robertson-Tessi, E.G. Moros, M. Fishman, J.Y. Djeu, S.E. Finkelstein, H. Enderling

Study supervision: H. Enderling

The authors thank Drs. Alexander R. A. Anderson and Tom Sellers for the organization and support of the third Moffitt IMO workshop on Personalized Medicine, where the idea of this project was conceived.

This article is supported by the Personalized Medicine Award 09-33504-14-05 and 09-33000-15-03 from the DeBartolo Family Personalized Medicine Institute Pilot Research Awards in Personalized Medicine (PRAPM).

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.

1.
Kuznetsov
VA
,
Knott
GD
. 
Modeling tumor regrowth and immunotherapy
.
Math Comput Model
2001
;
33
:
1275
87
.
2.
Mikhak
Z
,
Strassner
JP
,
Luster
AD
. 
Lung dendritic cells imprint T cell lung homing and promote lung immunity through the chemokine receptor CCR4
.
J Exp Med
2013
;
210
:
1855
69
.
3.
Calzascia
T
,
Masson
F
,
Di Berardino-Besson
W
,
Contassot
E
,
Wilmotte
R
,
Aurrand-Lions
M
, et al
Homing phenotypes of tumor-specific CD8 T cells are predetermined at the tumor site by crosspresenting APCs
.
Immunity
2005
;
22
:
175
84
.
4.
Zitvogel
L
,
Apetoh
L
,
Ghiringhelli
F
,
Kroemer
G
. 
Immunological aspects of cancer chemotherapy
.
Nat Rev Immunol
2008
;
8
:
59
73
.
5.
Couzin-Frankel
J
. 
Breakthrough of the year 2013. Cancer immunotherapy
.
Science
2013
;
342
:
1432
3
.
6.
Reits
EA
,
Hodge
JW
,
Herberts
CA
,
Groothuis
TA
,
Chakraborty
M
,
Wansley
EK
, et al
Radiation modulates the peptide repertoire, enhances MHC class I expression, and induces successful antitumor immunotherapy
.
J Exp Med
2006
;
203
:
1259
71
.
7.
Lugade
AA
,
Moran
JP
,
Gerber
SA
,
Rose
RC
,
Frelinger
JG
,
Lord
EM
. 
Local radiation therapy of B16 melanoma tumors increases the generation of tumor antigen-specific effector cells that traffic to the tumor
.
J Immunol
2005
;
174
:
7516
23
.
8.
Finkelstein
SE
,
Fishman
M
. 
Clinical opportunities in combining immunotherapy with radiation therapy
.
Front Oncol
2012
;
2
:
169
.
9.
Finkelstein
SE
,
Iclozan
C
,
Bui
MM
,
Cotter
MJ
,
Ramakrishnan
R
,
Ahmed
J
, et al
Combination of external beam radiotherapy (EBRT) with intratumoral injection of dendritic cells as neo-adjuvant treatment of high-risk soft tissue sarcoma patients
.
Int J Radiat Oncol Biol Phys
2012
;
82
:
924
32
.
10.
Finkelstein
SE
,
Rodriguez
F
,
Dunn
M
,
Farmello
MJ
,
Smilee
R
,
Janssen
W
, et al
Serial assessment of lymphocytes and apoptosis in the prostate during coordinated intraprostatic dendritic cell injection and radiotherapy
.
Immunotherapy
2012
;
4
:
373
82
.
11.
Finkelstein
SE
,
Timmerman
R
,
McBride
WH
,
Schaue
D
,
Hoffe
SE
,
Mantz
CA
, et al
The confluence of stereotactic ablative radiotherapy and tumor immunology
.
Clin Dev Immunol
2011
;
2011
:
439752
.
12.
Vatner
RE
,
Cooper
BT
,
Vanpouille-Box
C
,
Demaria
S
,
Formenti
SC
. 
Combinations of immunotherapy and radiation in cancer therapy
.
Front Oncol
2014
;
4
:
325
.
13.
Demaria
S
,
Ng
B
,
Devitt
ML
,
Babb
JS
,
Kawashima
N
,
Liebes
L
, et al
Ionizing radiation inhibition of distant untreated tumors (abscopal effect) is immune mediated
.
Int J Radiat Oncol Biol Phys
2004
;
58
:
862
70
.
14.
Wersall
PJ
,
Blomgren
H
,
Pisa
P
,
Lax
I
,
Kalkner
KM
,
Svedman
C
. 
Regression of non-irradiated metastases after extracranial stereotactic radiotherapy in metastatic renal cell carcinoma
.
Acta Oncol
2006
;
45
:
493
7
.
15.
Ehlers
G
,
Fridman
M
. 
Abscopal effect of radiation in papillary adenocarcinoma
.
Br J Radiol
1973
;
46
:
220
2
.
16.
Kaminski
JM
,
Shinohara
E
,
Summers
JB
,
Niermann
KJ
,
Morimoto
A
,
Brousal
J
. 
The controversial abscopal effect
.
Cancer Treat Rev
2005
;
31
:
159
72
.
17.
Akutsu
Y
,
Matsubara
H
,
Urashima
T
,
Komatsu
A
,
Sakata
H
,
Nishimori
T
, et al
Combination of direct intratumoral administration of dendritic cells and irradiation induces strong systemic antitumor effect mediated by GRP94/gp96 against squamous cell carcinoma in mice
.
Int J Oncol
2007
;
31
:
509
15
.
18.
Kuznetsov
VA
,
Makalkin
IA
,
Taylor
MA
,
Perelson
AS
. 
Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis
.
Bull Math Biol
1994
;
56
:
295
321
.
19.
Wilkie
KP
. 
A review of mathematical models of cancer-immune interactions in the context of tumor dormancy
.
Adv Exp Med Biol
2013
;
734
:
201
34
.
20.
d'Onofrio
A
,
Ciancio
A
. 
Simple biophysical model of tumor evasion from immune system control
.
Phys Rev E Stat Nonlin Soft Matter Phys
2011
;
84
:
031910
.
21.
Arciero
JC
,
Jackson
TL
,
Kirchner
DE
. 
A mathematical model of tumor-immune evasion ans siRNA treatment
.
Discret Contin Dyn S
2004
;
4
:
39
58
.
22.
Kirchner
D
,
Panetta
JC
. 
Modeling immunotherapy of the tumor-immune interaction
.
J Math Biol
1998
;
37
:
235
52
.
23.
Kronik
N
,
Kogan
Y
,
Elishmereni
M
,
Halevi-Tobias
K
,
Vuk-Pavlovic
S
,
Agur
Z
. 
Predicting outcomes of prostate cancer immunotherapy by personalized mathematical models
.
PLoS One
2010
;
5
:
e15482
.
24.
Kronik
N
,
Kogan
Y
,
Vainstein
V
,
Agur
Z
. 
Improving alloreactive CTL immunotherapy for malignant gliomas using a simulation model of their interactive dynamics
.
Cancer Immunol Immunother
2008
;
57
:
425
39
.
25.
Cappuccio
A
,
Elishmereni
M
,
Agur
Z
. 
Cancer immunotherapy by interleukin-21: potential treatment strategies evaluated in a mathematical model
.
Cancer Res
2006
;
66
:
7293
300
.
26.
Dunn
GP
,
Old
LJ
,
Schreiber
RD
. 
The three Es of cancer immunoediting
.
Annu Rev Immunol
2004
;
22
:
329
60
.
27.
Lange
PH
,
Hekmat
K
,
Bosl
G
,
Kennedy
BJ
,
Fraley
EE
. 
Accelerated growth of testicular cancer after cytoreductive surgery
.
Cancer
1980
;
45
:
1498
506
.
28.
De Giorgi
V
,
Massi
D
,
Gerlini
G
,
Mannone
F
,
Quercioli
E
,
Carli
P
. 
Immediate local and regional recurrence after the excision of a polypoid melanoma: tumor dormancy or tumor activation
?
Dermatol Surg
2003
;
29
:
664
7
.
29.
Deylgat
B
,
Van Rooy
F
,
Vansteenkiste
F
,
Devriendt
D
,
George
C
. 
Postsurgery activation of dormant liver micrometastasis: a case report and review of literature
.
J Gastrointest Cancer
2011
;
42
:
1
4
.
30.
Okada
K
,
Makihara
K
,
Okada
M
,
Hayashi
EI
,
Fukino
S
,
Fukata
T
, et al
A case of primary squamous cell carcinoma of the breast with rapid progression
.
Breast Cancer
2000
;
7
:
160
4
.
31.
Valentin
J
. 
Basic anatomical and physiological data for use in radiological protection: reference values: ICRP Publication 89
.
Ann ICRP
2002
;
32
:
5
265
.
32.
Hausmann
H
,
Photiadis
J
,
Hetzer
R
. 
Blood flow in the internal mammary artery after the administration of papaverine during coronary artery bypass grafting
.
Tex Heart Inst J
1996
;
23
:
279
83
.
33.
Vorherr
H
. 
The breast: morphology, physiology, and lactation
.
New York, NY
:
Academic Press
; 
1974
.
34.
Abduljalil
K
,
Furness
P
,
Johnson
TN
,
Rostami-Hodjegan
A
,
Soltani
H
. 
Anatomical, physiological and metabolic changes with gestational age during normal pregnancy: a database for parameters required in physiologically based pharmacokinetic modelling
.
Clin Pharmacokinet
2012
;
51
:
365
96
.
35.
Puskas
Z
,
Schuierer
G
. 
[Determination of blood circulation time for optimizing contrast medium administration in CT angiography]
.
Radiologe
1996
;
36
:
750
7
.
36.
McGuire
A
,
Brown
JA
,
Kerin
MJ
. 
Metastatic breast cancer: the potential of miRNA for diagnosis and treatment monitoring
.
Cancer Metastasis Rev
2015
;
34
:
145
55
.
37.
Howlader
N
,
Noone
AM
,
Krapcho
M
,
Garshell
J
,
Miller
D
,
Altekruse
SF
, et al
SEER Cancer Statistics Review, 1975-2012, National Cancer Institute
.
Bethesda, MD.
Available from: http://seer.cancer.gov/csr/1975_2012/.
38.
O'Brien
MA
,
Power
DG
,
Clover
AJ
,
Bird
B
,
Soden
DM
,
Forde
PF
. 
Local tumour ablative therapies: opportunities for maximising immune engagement and activation
.
Biochim Biophys Acta
2014
;
1846
:
510
23
.
39.
Simone
CB
,
Jones
JA
. 
Multidisciplinary approaches to palliative oncology care
.
Ann Palliat Med
2014
;
3
:
126
28
.
40.
Norton
L
. 
Conceptual and practical implications of breast tissue geometry: toward a more effective, less toxic therapy
.
Oncologist
2005
;
10
:
370
81
.
41.
Demicheli
R
,
Retsky
MW
,
Hrushesky
WJ
,
Baum
M
,
Gukas
ID
. 
The effects of surgery on tumor growth: a century of investigations
.
Ann Oncol
2008
;
19
:
1821
8
.
42.
Qadri
SSA
,
Wang
J-H
,
Coffey
JC
,
Alam
M
,
O'Donnell
A
,
Aherne
T
, et al
Can surgery for cancer accelerate the progression of secondary tumors within residual minimal disease at both local and systemic levels
?
Ann Thorac Surg
2005
;
80
:
1046
51
.
43.
Shafirstein
G
,
Kaufmann
Y
,
Hennings
L
,
Siegel
E
,
Griffin
RJ
,
Novak
P
, et al
Conductive interstitial thermal therapy (CITT) inhibits recurrence and metastasis in rabbit VX2 carcinoma model
.
Int J Hyperthermia
2009
;
25
:
446
54
.
44.
Kim
Y
,
Boushaba
K
. 
Regulation of tumor dormancy and role of microenvironment: a mathematical model
.
Adv Exp Med Biol
2013
;
734
:
237
59
.
45.
Hanin
L
. 
Seeing the invisible: how mathematical models uncover tumor dormancy, reconstruct the natural history of cancer, and assess the effects of treatment
.
Adv Exp Med Biol
2013
;
734
:
261
82
.
46.
Siegel
R
,
DeSantis
C
,
Virgo
K
,
Stein
K
,
Mariotto
A
,
Smith
T
, et al
Cancer treatment and survivorship statistics, 2012
.
CA Cancer J Clin
2012
;
62
:
220
41
.
47.
Scott
J
,
Kuhn
P
,
Anderson
AR
. 
Unifying metastasis–integrating intravasation, circulation and end-organ colonization
.
Nat Rev Cancer
2012
;
12
:
445
6
.
48.
Henderson
JM
,
Heymsfield
SB
,
Horowitz
J
,
Kutner
MH
. 
Measurement of liver and spleen volume by computed tomography. Assessment of reproducibility and changes found following a selective distal splenorenal shunt
.
Radiology
1981
;
141
:
525
7
.
49.
Puybasset
L
,
Cluzel
P
,
Chao
N
,
Slutsky
AS
,
Coriat
P
,
Rouby
JJ
. 
A computed tomography scan assessment of regional lung volume in acute lung injury. The CT Scan ARDS Study Group
.
Am J Respir Crit Care Med
1998
;
158
:
1644
55
.
50.
Poggio
ED
,
Hila
S
,
Stephany
B
,
Fatica
R
,
Krishnamurthi
V
,
del Bosque
C
, et al
Donor kidney volume and outcomes following live donor kidney transplantation
.
Am J Transplant
2006
;
6
:
616
24
.