Immune interconnectivity of anatomically distant tumors as a potential mediator of systemic responses to local therapy

Complex interactions occur between tumor and host immune system at each site in the metastatic setting, the outcome of which can determine behavior ranging from dormancy to rapid growth. An additional layer of complexity arises from the understanding that cytotoxic T cells can traffic through the host circulatory system. Coupling mathematical models of local tumor-immune dynamics and systemic T cell trafficking allows us to simulate the evolution of tumor and immune cell populations in anatomically distant sites following local therapy and thus computationally evaluate immune interconnectivity. Results suggest that the presence of a secondary site may either inhibit or promote growth of the primary, depending on the capacity for immune recruitment of each tumor and the resulting systemic redistribution of T cells. Treatment such as surgical resection and radiotherapy can be simulated to estimate both the decrease in tumor volume at the local treatment-targeted site, and the change in overall tumor burden and tumor growth trajectories across all sites. Qualitatively similar responses of distant tumors to local therapy (positive and negative abscopal effects) to those reported in the clinical setting were observed. Such findings may facilitate an improved understanding of general disease kinetics in the metastatic setting: if metastatic sites are interconnected through the immune system, truly local therapy does not exist.


S1. Reduction of model complexity
The original model framework of Kuznetsov et al. introduced a term describing the rate (1-p) at which immune cells die in interactions with tumor cells, the optimal value of which was found to be 0.002. 1 Repeating the model fitting to experimental data presented in the original manuscript 1,2 in the presence and absence of this term (i.e. with (1-p) equal to 0.002 and 0, respectively), demonstrated an almost identical ability of the model to qualitatively represent the data when this term is removed (Figures 1, 2). To quantify the benefit of omitting this term the Akaike information criterion 3 (AIC) could be used to measure the relative "quality" of each version of the model (a = original 7 parameter model, b = new 6 parameter model). The AIC evaluates the trade-off between model complexity in terms of number of parameters, and the accuracy of model fit for a specified data set, and can be calculated as follows: where n is the number of data points (observations), M is the number of parameters in the model, and the latter term provides a bias-adjustment for small sample sizes. The first term represents the log likelihood of each model, approximated by the residual sum of squares. AIC analysis suggests that the simpler model with one less parameter is the "better" model with a relative likelihood of , R the number of models being compared, and Δ ! = ! − min ( ). 3 This suggests that the 6 parameter model (model b) is comparably representative of the data while reducing risk of overfitting. Parameters obtained from fitting the original data set to this new model are provided in Table S1, along with the comparable values from the original model of Kuznetsov et al. Based on this comparison, the simpler model was used for the purposes of this manuscript.

S2. Identification of growth parameters
Sensitivity analysis was conducted to identify the parameters for which a small perturbation induces the largest change in the dynamics of the system. 4 Relative local sensitivities S i were calculated for each parameter i using the standard definition where C is the cancer cell population and p i the ith parameter. 5 Approximations of the derivative were calculated using the second-order central finite differencing method. These logarithmic sensitivities were calculated locally (based on perturbations to each parameter of ±1%) for the default parameter values and for 10,000 random parameter sets generated using Latin hypercube sampling. 4,6 All sensitivities were calculated at 10 evenly spaced time points during a 1000 day simulation period and averaged prior to ranking. Parameters δ, γ, and σ and were found to lie within the top 3 most sensitive parameters in 92%, 90% and 58% of the 10,000 simulations, respectively ( Figure S3), with δ most commonly ranking 1 st, γ most commonly ranking 2 nd , and α most commonly ranking 3 rd .
A parameter sweep was then conducted to identify the values of these parameters that could induce constant growth behavior over our time period of interest (1000 days). This involved allowing each parameter to vary from 0 to 1 in increments of 0.05. Simulations were run with all combinations of the three parameters (N=200^3), and combinations for which + 1 > ∀ ∈ 0,1000 were identified. These parameter combinations and the distribution of each parameter within these combinations are shown in Figure S4. We arbitrarily select δ = 0.25, γ = 0.75 and σ ∈ 0.02,0.2 for presentation of results in the present manuscript, although other combinations within these distributions could have been equally demonstrative. Figure S3. Proportion of all simulations for which each respective parameter ranked 1 st , 2 nd , or 3 rd . Figure S4. Distribution of all combinations of parameters (α, γ, δ) for which continuous growth behavior was observed over our time period of interest.

S3. Recruitment term
The following description is taken from Poleszczuk et al. and the reader is directed to [7] for additional detail. Four compartments (COMPs) of the blood system are described: pulmonary circulation (LU); gastro-intestinal tract and spleen (GIS); liver (LI); and all other organs in the systemic circulation such as breast or kidney (SO). Matrix !" describes the probabilities of T cells activated at site j infiltrating site i, and is given by where !"#$ = /Δ with and Δ is a normalizing constant such that the sum of absorption probabilities over all four compartments is equal to one. BFF COMP is the blood flow fraction to the respective compartment, and V i is the volume of the i th organ. Probability !" is the probability that a T cell activated at site j will infiltrate the i th tumor site after entering a given compartment. This is calculated by the probability that a T cell will flow through a tumor site (approximated by relative blood flow through the tumor bearing organ multiplied by the fraction of the organ populated b the tumor), multiplied by ℎ !" , the probability of extravasation from the blood into the tissue: Here, ℎ !" = ℎ ! if the T cell was activated in the given organ (organ i = organ j ) and ℎ !" = ℎ ! otherwise (1 ≥ ℎ ! ≥ ℎ ! ). Finally, !"#$ is the probability that a T cell will extravasate at one of the metastatic sites after entering a given compartment: Relevant parameters for our tumor sites of interest are as follows: (V lung = 3679 mL, V kidney = 249 mL, V breast = 500 mL) [8][9][10][11] and oxygenated arterial blood flow fraction (BFF lung = 100%, BFF kidney = 8.5%, BFF breast = 2%) 7,12 .

S4. Radiation-associated parameters
Based on observations in the breast, α = 0.3 and α/β = 10 were arbitrarily chosen for demonstrative purposes. 13 Figure S5 demonstrates the influence of varying these radiationassociated parameters on model output. While the model appears to be relatively insensitive to the value of the ratio α/β, the specific α value can have a more significant effect. For future testing, site-specific LQ parameters should be determined for each respective metastatic site of interest for maximum accuracy.  Figure S5. Influence of variation of each respective radiation-associated parameter on model outcomes.  (c,d). Resection of a breast primary results in an increased dormant volume of a lung metastasis (e), attributable to a decrease in immune surveillance at the lung site. Resection of the lung metastasis results in significantly decreased dormant volume of a breast primary (f), attributable to an increase in immune surveillance at the breast site. Irradiation of the breast site (g) results in not only shrinkage of this site, but also an initial, transient abscopal effect prior to the lung established an equilibrume volume greater than that in the presence of the breast tumor. Irradiation of the lung site (h) allows the previously redirected immune response to return to its site of activation in the breast, inducing a prolonged abscopal response here. In all dormancy modeling, parameters were as defined in Table 1 for the 6 parameter model. Figure S7. In the case of a primary breast tumor having reached a dormant equilibrium with the host immune system, the seeding of an additional site in the lung induces both an increase in the number of activated effector cells circulating systemically, and the redistribution of effector cells previously surveilling the breast site to the new lung metastasis. This has the potential to reduce the immune surveillance at the breast site, inducing a transient period of growth until dormancy is re-established at a higher equilibrium volume. This pattern can be repeated until such point that the immune surveillance of the breast site reaches a minimum, or the physiologic base level of T cell surveillance in healthy tissue. This demonstrates the feasibility of an metastasis-seeding induced escape from dormancy mechanism, distinct from previously described processes of tumor growth.  Metastatic seeding increases the total systemic number of effector cells, and increases the number of immune cells extravasating at the primary lung tumor due to its higher blood flow fraction. These results with initial tumor volume of C i (t=0)=5x10 9 are qualitatively similar to those observed with an initial tumor volume of C i (t=0)=5x10 5 (Figure 2), suggesting the results presented herein do not depend on choice of initial tumor volume.