Mathematical Modeling of Tumor–Tumor Distant Interactions Supports a Systemic Control of Tumor Growth

Concomitant tumor resistance (CR) is a phenomenon by which the presence of a tumor negatively influences the appearance and growth of another implant (see ref. 1 for a review). It has been reported in numerous experimental studies spanning over more than one century, and implementing a large variety of animal models (2–12). These studies investigated either the concomitant or subsequent implantation of a second graft after a primary injection (3, 5, 6, 8, 9, 12), or the inhibition of secondary tumors arising from the primary (metastases; refs. 4, 13–16). They consistently evidenced a systemic growth suppression effect, demonstrated by the occurrence of postremoval growth acceleration (5, 6, 10, 13, 14, 16–19). In the clinic, suppression of the growth of metastases by the presence of the primary tumor has yet to be appreciated in general therapeutic planning, although it has been reported in patients (17, 20–22). Despite the abundance of reports of this phenomenon, the precise determinants of CR remain poorly understood, and only qualitative theories have been advanced.

CR is of direct clinical relevance insofar as it implies that removal of a primary tumor, with the resultant release of its inhibitory pressure on occult secondary sites, could be followed by postsurgery metastatic acceleration (PSMA). PSMA has been demonstrated to occur in numerous animal experiments (5, 6, 13, 16, 19), as well as in clinical case reports (20, 21, 24). Further support for the occurrence of PSMA in a notable fraction of patients was also provided by the observation of two peaks in the hazard relapse rate of a large cohort of breast cancer patients (22–25).

Several hypotheses have historically been proposed for the explanation of the underlying mechanism of CR. The first was due to P. Ehrlich and consisted in athrepsia, i.e., that two (or more) tumors in the same organism would compete for nutrients and that the growth of one tumor would leave less nutrients available for the other (7, 30). However, this was challenged by the observations that CR was decreased when the number of inoculated cells was increased (9). Another popular theory, first introduced by Bashford in 1908, was based on immunologic mechanisms and stipulated that the presence of a first tumor would activate an immune response preventing the second graft to take or grow (3, 7). However, several studies demonstrated the occurrence of CR in tumor models with no or weak immunogenicity, or in immune-deprived mice, thus challenging this explanation (4, 7–9). This implies that, although immunologic factors might contribute to CR, they cannot explain it entirely.

In the 1990s, a team led by J. Folkman discovered endogenous inhibitors of angiogenesis (such as endostatin or angiostatin) by demonstrating that injection of these factors could substitute for the suppressive effect on lung metastases exerted by the primary tumor (14, 31). This led the investigators to link CR to distant systemic inhibition of angiogenesis. Their theory is based on the fact that locally produced angiogenesis inhibitors would spread systemically using the vascular system, reach secondary lesions, and impair angiogenesis at the distant sites (11, 12, 17, 32), eventually overcoming the influence of angiogenesis stimulators also produced by the tumors. This idea is supported by the fact that endogenous angiogenesis inhibitors have long half-lives compared with stimulators, which allows them to diffuse, reach the vascular system, and accumulate (14, 33).

The idea of circulating inhibiting factors due to the presence of a primary tumor had been proposed and confirmed in earlier studies (1, 4, 7), but their precise mode of action has remained elusive. Distinct from the angiogenesis inhibitors previously mentioned, another research group identified other blood-borne factors with direct antiproliferative action, namely meta- and ortho-tyrosine, that would reduce proliferation by driving tumor cells into a G₀-phase state of dormancy or induce an S-phase arrest (27, 28).

So far, arguments and theories about CR have remained qualitative in nature. In the present study, comparing alternative mathematical formalisms, we demonstrate that a simplified model with well-motivated parameters that addresses concomitant resistance specifically was able to capture features of coupled tumor growth and may shed more light on the understudied subject of systemic controls in cancer, a potentially critical step toward eventually understanding metastatic control.

Murine Lewis lung carcinoma (LLC) cells, originally derived from a spontaneous tumor in a C57BL/6 mouse (34), were obtained from the ATCC in the period 2005–2006. The LLC cells were cultured under standard conditions (34) in high-glucose DMEM (Gibco Invitrogen Cell Culture) with 10% FBS (Gibco Invitrogen Cell Culture) and 5% CO₂. They were expanded through 4 passages (LLC p4) and then aliquoted and stored in liquid nitrogen. At that time, they underwent Molecular Diagnostics Infectious Disease PCR testing (Mouse Essential Panel) at Charles River Laboratories. An aliquot of the above LLC p4 cell line was thawed and tested for Mycoplasma at Bionique Laboratories for Mycoplasma testing. They were found to be negative for mycoplasma. The cells have not been characterized. An aliquot of LLC p4 was used in these experiments. LLC p4 cells were thawed and passaged through an additional two cycles prior to injection into the mice.

C57BL/6 male mice were used with an average lifespan of 878 days (35). At time of injection, mice were 6 to 8 weeks old (Jackson Laboratory). Subcutaneous injections of 10⁶ LLC cells in 0.2 mL PBS were performed on the caudal half of the back for the control group and on the two lateral sides of the caudal half of the back for the group bearing two tumors, in anesthetized mice. Tumor size was measured regularly with calipers to a maximum of 1.5 cm³ when mice were sacrificed. Animal tumor model studies were performed in strict accordance with the recommendations in the Guide for Care and Use of Laboratory Animals of the NIH and according to guidelines of the Institutional Animal Care and Use Committee at Steward Research and Specialty Projects Corp.

Simulations of the mathematical models were performed using Matlab with statistics and optimization toolboxes (The Mathworks Inc., 2015). Fits of the models to the data were executed using an internally developed software that utilizes the function fminsearch for weighted least squares minimization. The objective function was computed using weights defined by an error variance model previously established on the same animal model and measurement technique (36). Statistical analyses of the fits (computation of the goodness-of-fit metrics and standard errors of the parameter estimates) were implemented in our software as previously described (36, 37).

We studied the phenomenon of CR by combining experiments and mathematical modeling, informed by pre-existing theories in the literature. For the experiments, two groups were considered. The first group (control) consisted of 20 mice in which single implants were performed. In the second group (double tumors, abbreviated as DT) consisting of 10 mice, 2 grafts with identical load (10⁶ cells) were performed on day 0, at the same locations on opposite flanks of the mice.

We first performed a direct (i.e., not model-based) statistical analysis of the data. We compared control tumor growth kinetics in mice bearing single implants with the growth curves of tumors in a double-tumor–bearing host (Fig. 2). Observations of the kinetics of each tumor in the DT group suggested that in each mouse, one tumor was growing faster than the other, possibly inhibiting the second one (Fig. 2A; Supplementary Fig. S1A and S1B). This behavior was observed consistently in all the animals of the DT group except in two of them (animals 2 and 9 in Fig. 2A) and did not seem to result from the lateral location (left or right) of the tumors (Fig. 2A). Intriguingly, the two mice where the phenomenon was not observable were found to have a connecting blood vessel joining the two tumors and were the only ones to exhibit this macroscopic vascular structure. Direct sharing of same vasculature seemed to equilibrate tumor expansions. One possible explanation for the absence of cross growth suppression could be that the production of inhibitors was negligible in these tumor-host systems. This could also explain the formation of the connecting blood vessel due to increased neo-angiogenesis (under the theory of angiogenesis inhibition-driven growth suppression).

In order to statistically confirm that this observation was not purely due to intrinsic randomness in experimental conditions (such as the number of initial cells that “take” from the injection) that would by themselves generate different growth curves for the two implants in the same mice, we performed a statistical analysis. It aimed at testing the null hypothesis that both tumors would be identically and independently growing (i.i.g.). We artificially generated 10 couples of i.i.g. tumors by subdividing the control group of 20 animals into two groups of 10, randomly picking tumors from each subgroup and pairing them together. We then picked each small tumor from these pairs (and similarly, each large tumor), by choosing the one with smallest final volume. This yielded two samples of 10 “control small” and “control large” tumors that can be considered as what would have emerged from randomness only in initiation and growth. These two samples could then be statistically compared with the experimental samples of small and large tumors from the DT group. We observed significant differences between the small tumors from the DT group and the “control small” tumors from day 12 until the end (P < 0.05, Student t test), with the exception of days 18 and 19 where differences were not significant due to large variability (Fig. 2B and C). On the other hand, no statistical difference was observed between the large tumors from the DT group and their control counterparts.

These results demonstrate that in a mouse where two tumors are simultaneously growing, the larger tumor is growing at the same speed as would a single tumor, whereas the other has significantly slower kinetics.

To investigate what caused the behavior established in the previous section, we fitted various growth models previously shown able to describe syngeneic LLC single tumor growth with identifiable parameters (36). Models that were considered here were: exponential growth (with free initial volume |{V_0}$|⁠), i.e., |{\frac{{dV}}{{dt}}} = aV,\ V( {t = 0} ) = {V_0}$|⁠, power law, i.e., |{\frac{{dV}}{{dt}}} = \alpha {V^\gamma },\ V( {t = 0} ) = 1$|⁠, and Gompertz, i.e., |{\frac{{dV}}{{dt}}} = ( {{\alpha _0} - \beta \ln ( {\frac{V}{{{{10}^{ - 6}}}}} )} )V,\ V( {t = 0} ) = 1$|⁠. For the two last ones, for identifiability reasons, the initial volume was fixed to the number of injected cells (10⁶ cells ≈ 1 mm³). The objective was to see if any of these models was able, through estimation of its coefficients, to elucidate the observed differences in the growth of the two tumors in the same animal. In this first context, each tumor was fitted independently, resulting in one parameter set per tumor, and we looked for significant differences in the parameters between the group of small (or large) tumors and the group of small (or large) tumors built from pairings of tumors from the single-tumor control group.

As expected from their accurate descriptive power, all the models performed well in terms of goodness of fit and were able to reproduce the tumor growth curves (Supplementary Figs. S2–S4).

Parameter estimates are reported in Fig. 3. The only parameter that exhibited a significant difference between the simultaneous and control groups was the proliferation rate |a$| of the exponential model. This suggests an alteration of proliferative abilities in the altered tumor due to the presence of the other one. On the other hand, the time t_D to reach a given volume V_D (of 100 mm³ for instance) resulting from simulations of the calibrated models was steadily found significantly larger in the small tumors of the two-tumor animals, as compared with the small control tumors (Supplementary Fig. S5).

In summary, single tumor growth models were able to give an appropriate description of the growth curves taken independently and suggest alteration of the proliferative abilities in the small tumors group. However, they do not offer any insight into the dynamics of interactions between the two implants.

To next study and quantify the possible underlying mechanisms of tumor–tumor interactions, we investigated models of simultaneous growth of two tumors in the same organism. A virtually infinite number of models can be conceived for description of CR, both in terms of structural shape (equations) and links between the parameters. We report here only on the results from analysis of five informative models. Supplementary Table S1 details the equations of these models. Interestingly, several models were found unable to fit the data, suggesting rejection of (at least one of) the hypotheses that they rely on. These include a model for the athrepsis theory (competition for nutrients), as shown in Supplementary Fig. S6. The other models were based on either systemic inhibition of angiogenesis (SIA)—formalized using a model of interactions between tumor growth and vascular support (38)—or induction of quiescence due to a cytostatic seric factor, as experimentally evidenced by Ruggiero and colleagues (9, 28). We termed the models based on this last theory proliferation inhibition (PI) models. The SIA model was able to give a reasonably accurate description of our data (Supplementary Fig. S7). For PI models, which have similar structures as equation A, three hypotheses were investigated (see Supplementary Table S1 for the detailed expressions of the models): (i) direct effect (a given quantity of IFs induces a given number of cells going to quiescence; equation A), (ii) log-kill effect (a given quantity of IFs induces a given fraction of cells going to quiescence), and (iii) total number of cells (P_i + Q_i) as a source of IFs. Notably, models 2 and 3 were unable to fit the data and had to be rejected (Supplementary Figs. S8 and S9). Residuals analysis supporting these results is shown in Supplementary Fig. S10A–S10D. On the other hand, the model 1 generated a particularly good fit to the data (Fig. 4A and B). Table 1 summarizes statistical quantitative metrics of goodness-of-fit models that allow comparison of their descriptive power, whereas Fig. 4C shows the distribution of the residuals. Table 2 reports the parameter values of all the models estimated from the best fits, together with their interanimal variability (parameters were individually fitted for each mouse) and standard errors for the estimates.

These results demonstrate that a mathematical modeling approach was able, by confrontation of the best-fits of the model, to discriminate among qualitatively equally likely theories of CR and suggest a PI model as having attributes that may explain (or at least describe) this particular phenomenon.

The PI model 1, formalized by equation A, consists in assuming a direct and mutual growth rate decrease between the two tumors, due to passage to quiescence (Fig. 1). Goodness of fit was found excellent (Fig. 4), as well as identifiability of the parameters (see standard errors in Table 2). Notably, while being fitted directly on the two tumors growths, the predicted behavior when simulating no interactions was in full agreement with the control growth curves. Indeed, the dashed lines in Fig. 4A are close to the growth curves of the large tumors, which were found to be not significantly different from “naturally happening” large tumors. Hence, the model was able to learn and identify the unaltered growth part from altered growth curves, highlighting its reliability.

In mouse number two (second plot in the top row of Fig. 4A), consistently with the observation of identical growth kinetics between the two tumors, the model identified a value of parameter |\gamma $| not significantly different from zero. On the other hand, the model did identify interactions between the two tumors in the other animals, as emphasized by 95% confidence intervals of parameter |\gamma $| inferred from the parameter estimation process that did not contain 0 [(0.0374, 0.0436); 95% confidence interval for the median value of γ in our estimation]. In turn, this translated into substantial differences in the kinetics of the two implanted tumors (see Fig. 4A where growth curves are plotted with or without interactions). Of important note, this difference in the kinetics was mostly due to the interaction between the two tumors, rather than the initial difference in cell loss. This is demonstrated in Fig. 4A, where the growth curves of the small tumors with only a different initial volume (dashed curves) were considerably higher than the curves, where interaction was taken into account. Moreover, these curves were both close to the large tumor growth curves, indicating that the difference in |{V_0}$| had only a negligible impact on the difference between the two growth curves, the major determinant being the tumor–tumor cross-inhibition effect. Critically, the differences for the large tumor curves were much smaller than for the small tumor curves (although the interaction parameter was the same for both tumors), indicating that the model gives a valid quantitative theory of why only one tumor was affected by CR.

Together, our results provide a biologically based and minimally parameterized mathematical model for tumor growth kinetics interactions in a two-tumor–bearing host. The model confirmed a significative nonzero value for the interaction parameter in 9 of 10 mice, which gave a quantitative measure of the phenomenon. Asymmetry between the two tumors was explained by an initial difference in the take between the two implants.

In addition to being able to describe double tumor growth and CR, the elementary model that we propose also offers a simple formalism for single tumor volume growth. The model consists in the division of the cancer cells into two subpopulations: proliferative and quiescent cells (which could also account for necrotic tissue still present in the total volume measurement). Indeed, in the absence of a secondary tumor, the model 1 becomes:

This provides a valid and simple mathematical construct able to describe the growth of single tumors (Fig. 5A and B; Supplementary Fig. S11). It sheds new lights on general tumor growth laws as it demonstrates that classical Gompertzian growth—which is able to describe accurately Lewis Lung tumor growth curves (36)—can be reproduced by these equations, with no significant differences (i.e., a difference in Akaike Information Criterion less than 2, see Supplementary Table S2). Indeed, it had remained elusive why the Gompertzian curve, which was originally designed not even for growth processes (39), describes tumor growth curves and their consistent relative growth rate decrease with such important accuracy, while being only phenomenological and not biologically grounded. Interestingly, we obtained that a model where growth deceleration was assumed to result (only) from passage to quiescence due to the production of factors by the proliferative tumor cells themselves was able to explain single tumor growth curves as accurately as the Gompertz model, or other models such as the power law (Fig. 5A and B). Parameters identifiability of this new model was also very good (Supplementary Table S3).

To quantitatively inform on CR, we performed an integrative study that combined experimental and theoretical investigations. At the experimental level, we found that when two identical tumors were implanted in the same organism, distant interactions occurred. Specifically, in one (and only one) of the two tumors, growth was suppressed, whereas the other remained unaltered. Several theories formalized by mathematical models were assessed against the data by means of appropriate mathematical constructs. Our findings revealed that direct inhibition of proliferation was able to appropriately describe the data. Under this hypothesis, both tumors release systemic factors that mutually induce quiescence in the other tumor, which is in line with experimental results identifying these factors as meta- and ortho-tyrosine (27, 28). In our study, the origin of the asymmetry between the two tumors was hypothesized to emerge from stochastic fluctuations in the initial number of cells that take between the two tumors, with variations of about 25% between the two. Critically, although simulating this difference in initial condition resulted in only small changes when no interaction was taken into account, it generated substantial discrepancies when it was. Our model thus gives a dynamical explanation of the observation of one tumor “winning on the other,” because one and only one of the two tumors gets substantially suppressed despite equal growth and interaction parameters. Notably, within the context of PI, other models such as a cytotoxic (log-kill) effect or production of inhibitory factors by the entire tumor burden had to be rejected. This suggests experimentally testable predictions in the mode of production (only by proliferative cells) and action of the inhibitory factors.

Arguments disregarding the competition theory had already been put forward by others (3, 7, 40). For example, Gorelik had argued that under this theory, the intensity of CR should be an increasing function of the amount of cells implanted in a subsequent graft, in contradiction with experimental findings, thus disqualifying the theory (7). However, these arguments had remained qualitative. Our formal study adds a quantitative basis to these considerations by showing that, under the modeling assumptions we operated, this theory was unable to accurately describe our data (Supplementary Figs. S6 and S10A).

More elusive in the literature had remained the question to discriminate between angiogenesis inhibition, as evidenced by the work of Folkman and colleagues (14, 32), and direct induction of quiescence by seric factors, as proposed by Ruggiero and colleagues (9, 27, 28). Our results suggest that the latter theory, when considered alone, could be sufficient to drive CR, insofar as it exhibited good match to the data. However, this does not preclude SIA to occur, because the two theories are not mutually exclusive. Mutual nonexclusivity also applies to the competition theory: it cannot be completely disregarded that a combination of the three phenomena happens in the occurrence of CR. However, it is beyond the scope of the current study to be able to disentangle between a combination of the phenomena and one phenomenon alone.

In our analysis, we did not address the implication of the immune system in CR, despite the reasonable belief that immune players could be involved. In fact, CR was originally termed "concomitant immunity," because it was believed that CR occurred principally by triggering an immune response from the presence of the first tumor (3). However, several subsequent studies demonstrated that the phenomenon of CR could occur in a broad range of situations where the immune system could not be the main driver, including studies using the animal model employed here (4, 7–9). For example, in ref. 11, the authors demonstrated that, in T-cell–depleted mice bearing LLC, the CR effect was equally effective as in normal mice, and it was tumor nonspecific. Similar results were obtained when studying CR in mice with functions of macrophages and natural killer cells suppressed (8). Therefore, given the amount of evidence that nonimmunological phenomena are involved in the generation of concomitant resistance, and for the sake of simplicity, we decided to focus here on theories that did not require the intervention of immune players and postpone this to future work.

Our findings not only shed light on the dynamics of CR but also proposed a simple and biologically based model of (double and) single tumor growth able to describe the ubiquitously observed growth retardation with larger volumes (41–44), usually modeled by means of the Gompertz equation (41, 45). Although several attempts of deriving the Gompertz equation from basic principles have been performed in the literature (46, 47), the model we propose here benefits from its simplicity. For single tumor growth, it has only two (aggregated) parameters, which quantify two phenomena: (i) proliferation of the active cancer cells and (ii) production by these active cells of factors that drive them to quiescence. Interestingly, this model brings new light on the so-called paradox of CR (40), which can be expressed as follows: if distant inhibition occurs, potentially driving other tumors to dormancy, then why does the primary tumor continue growing? Our general model gives a way to quantitatively formalize this. Indeed, the same factors act as local and distant inhibitors, and we showed that, under appropriate values of the parameters (but identical for the two tumors), one could obtain at the same time almost unaltered growth of the large (primary) tumor and significant suppression of the growth of the small (secondary) tumor. The presence of endogenous molecules with inhibitory activity thus challenges a naïve view where growth retardation would only be due to interactions dictated by competition (for space or nutrients). Consistently, such a model (logistic growth) had already been shown unable to adequately fit experimental tumor growth curves (36). Considering the implications, the mere fact a tumor would produce both angiogenesis stimulators and inhibitors at the same time, with near and far ranges, does not readily reconcile with a purely localized purpose, but instead speaks to tumor control being manifestly a systemic phenomenon, quite distinct from the naïve concept of an entity governed by local conditions alone, independent of other tumor sites. Implicit in this, and as proposed by others (40, 48), a vision of tumor growth as an integrated, organ-like development could bring sense to this seeming paradox.

Concomitant resistance and consequent PSMA have important clinical implications (21–23). In some instances, surgery might not be the optimal therapeutic option, and it might be more beneficial to try to control the patient total tumor burden (49). In addition, it has been experimentally demonstrated that preoperative administration of chemotherapy or radiotherapy could prevent this acceleration, thus giving a strong rationale for neoadjuvant therapy (50). Although PSMA has been reported experimentally since more than one century (2), consideration of this phenomenon in the clinic has remained limited. After testing of more than 40 distinct mathematical models of CR, we dispose now of a validated model for the interaction of two tumors. In parallel, we have developed since several years mathematical models for the systemic dynamics of metastases (51–54). Our next step is thus to integrate in these models the nontrivial interactions between the primary tumor and the metastases, or among the metastases themselves, as captured by the mathematical model defined here. When adequately validated, this will provide a computational tool of valuable clinical interest, as it could lead to personalized quantification of the impact of CR in patients. In turn, this will bring the opportunity to simulate various therapeutic scenarios and predict the potential occurrence of PSMA, paving the way to patient-specific adaptation of neoadjuvant and adjuvant therapy.

No potential conflicts of interest were disclosed.

Conception and design: S. Benzekry, D. Barbolosi, L. Hlatky, P. Hahnfeldt

Development of methodology: S. Benzekry, D. Barbolosi, L. Hlatky

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Benzekry, C. Lamont

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Benzekry

Writing, review, and/or revision of the manuscript: S. Benzekry, D. Barbolosi, L. Hlatky, P. Hahnfeldt

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Benzekry

Study supervision: L. Hlatky, P. Hahnfeldt

We thank Dr. Sauveur Merlenghi, Dr. Marie Dominique Battesti, and N. Spinosi from the Ligue Contre le Cancer de Corse du Sud for their support.

This project was supported, in part, by the NCI, under award Number U54CA149233 to L. Hlatky.

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.