Modeling Boundary Conditions for Balanced Proliferation in Metastatic Latency

The dormancy dilemma impacts millions of cancer patients worldwide. Certain epithelial carcinomas such as from breast or prostate tissues can recur years or decades after primary tumor ablation (1). This late recurrence afflicts nearly 50% of breast cancer patients that develop metastases; other carcinomas exhibit similar late recurrence rates (2). This recurrence is most problematic within metastatic organs such as the liver, lung, bone, and brain where tumor resection may not be possible and chemotherapeutic agents seem to be ineffective in most cases (3, 4). For instance, aggressive treatment (e.g., adjuvant chemotherapy) of breast cancer after surgical removal of the evident primary tumor may reduce 10-year recurrence by less than one third. Consequently, ostensibly cured cancer patients harbor unseen metastases for decades that eventually emerge into untreatable, lethal tumors. This long period between primary tumor treatment and metastatic appearance is referred to as metastatic latency though the tumor cell behavior during this period is unknown.

Metastatic latency allows disseminated tumor cells to avoid clinical detection and withstand aggressive and toxic neoadjuvant treatment (5). How the disseminated carcinoma cells maintained a long-term, clinically undetectable state in an ectopic site such as the liver is unknown (6). Although extravasation of carcinoma cells from the primary tumor is common, but most of these cells seem not to establish metastases (7, 8), it is likely that the premetastatic microenvironment contributes to whether carcinoma cells will seed, survive, and either proliferate to prompt tumor formation or exhibit a dormant phenotype.

Because in vitro and in vivo experimental systems are not sufficiently malleable to consistently and reliably recapitulate the metastatic latency phenotype, in silico models have evolved to explore the mechanisms of metastatic competency, dormancy, and emergence (9–11). Many of these in silico models predominately simulate metastatic latency according to the theory of balanced proliferation among micrometastases of greater than one million cells. Therefore, they introduce various biological factors such as angiogenesis, immune response, hypoxia, and growth factor availability and deduce through partial differential equations what impact these have on the ability for the micrometastases to emerge or maintain a dormant state. A major limitation to this approach is that little empirical biological data are available to make these models translatable to the in vivo response. Moreover, such studies often assume that latent metastases experience balanced proliferation, meaning additional biological data, even if available, would not resolve latency as being due to balanced proliferation, cellular quiescence, or a combination of both.

Notwithstanding these limitations, in silico models are becoming increasingly important to offering insights in prioritizing the landscape of contributing metastatic factors. By using an Approximate Bayesian Computation model applied to empirical patient data Willis and colleagues inferred that metastatic breast cancer late relapse evolved from just 1 to 6 micrometastases that escape from dormancy (12). Other approaches used the Gompertzian growth function to calculate periods of metastatic growth and growth arrest to fit clinical data as presented in the Munich Cancer Registry (13). These results predicted that metastatic seeding occurred before the clinical detectability of the primary tumors. Computational models have also been used to predict the metastatic response according to primary tumor removal and front line chemotherapy (14).

Biologically, the issue of metastatic dormancy has been approached as balanced proliferation and death, quiescence, or a combination of both. In arguing for balanced proliferation and death, investigators invoke various constraints, such as the “angiogenic switch,” before which the metastatic nodule is constrained to a certain mass in the absence of a vascular supply (2, 15). In other considerations, multiple feedback mechanisms are invoked to maintain a small, but cycling cell population (16). The cellular and/or molecular mechanisms behind tumor dormancy present with an absence of solid clinical data, although emerging evidence of phenotypic plasticity suggests some other possible routes to dormancy and may suggest quiescence of the disseminated tumor cells.

Thus, the situation remains unsettled as to which behavior predominates during these latent periods. To take an unbiased approach, we modeled metastatic growth dynamics by invoking a 2 state, Markov chain Monte Carlo (MCMC) simulation where only cell-survival probability and starting metastasis populations were varied. As depicted in Fig. 1, the MCMC method assigns a survival probability for each cell in each metastasis during each replication cycle. Cells will then either proliferate or die as statistically influenced by the survival probabilities allowing for metastases to exhibit a range of outcomes not predetermined. If metastatic dormancy can emerge across a wide spectrum of survival probabilities then balanced proliferation may be the dominant contributor. However, if the latent phenotype via balanced proliferation is only manifest within a narrow probability window we propose that cellular quiescence would likely participate nontrivially in protracted metastatic survival.

MCMC simulation is a computational method that samples from a probability distribution to assign system outcomes commensurate with the underlying distribution (17). This model is described as a branching process whereby each cell in cell cycle k will randomly give rise to 0 or 2 cells in the k + 1 cell cycle. We applied MCMC simulation as a means to model the fate of extravasated circulating tumor cells into the metastatic site. We developed 2 MCMC approaches to assign cell fates (Fig. 1). The first, a fixed probability model, assigned static survival probabilities to every cell in each metastatic nodule; the probability of survival for each individual cell was not dependent on the other cells in that metastatic nodule (Fig. 1).

The second approach incorporated a stochastic element where the survival probability for any given cell was sampled randomly between ±10 percentage points around the specified mean survival probability for the group. We refer to this variable probability as P∼ where appropriate to distinguish it from the assigned probability P. All results herein were developed from the second, stochastic model, diagrammed in Fig. 1A. Both formulations, however, are examples of branching processes, which allow for convenient analytic expressions for some aspects of tumor behavior. In our 2-state MCMC system, every cell in each metastasis could either divide or die per cell cycle, thus mimicking the balanced proliferation phenotype. Notably, cells were individually assigned a new survival probability at each cycle with the stochastically assigned probability adjustment within 10 percentage points above or below the initially assigned survival probability. For example, simulations at a fixed 60% survival probability will allow for random fluctuations between 50% and 70% survival. Although these survival probabilities coalesced to the expected means, this additional percentage point adjustment further imitated the dynamic survival conditions of the metastatic niche.

These simulated micrometastases could devolve to 0 cells (died out), exceed 1 million cells (clinical metastasis) or exceed 1,218 cycles with fewer than 1 million cells (dormant). Our data represent 100 compiled patient trials (100,000 micrometastases) for each combination of cell survival rate and starting cell number. For the progressive clinical metastases, we capture the range of cycles at which this occurred. For the dormant metastases we document the mean cell numbers at the time of dormancy. Although this model does not achieve pseudo-equilibrium given all simulated metastases terminate under one of 3 conditions, dormancy is nonetheless modeled by means of achieving neither the 0 or 1 million cell boundary before reaching 1,218 cycles; although it is likely that should there be no cycle limit, the cell boundaries will be reached.

We simulated patients in groups of 100 person cohorts across a matrix of cell survival probabilities (0.3 to 0.7, given as percentages) and number of initially extravasated cells. For each patient we modeled 1,000 carcinoma deposits to the metastatic organ at the time of primary tumor resection. The probability of an outgrowth is based on individual micrometastasis and not on patient load of total micrometastases. Each deposit started with 1, 2, 4, or 8 cells (later at 1,000 and 4,000 cells for boundary validation) to separate initial extravasation and seeding from subsequent metastatic behavior. We also assumed a one-and-one-half to 3-day cell cycle so that 1,218 cell cycles represented 5 to ten person years (dormancy) at a cycle time of 1.5 to 3 days per mitosis. Micrometastases that grew beyond 1 million cells were classified as approaching clinically evident and were not simulated further. Finally, we assumed that the initial 1,000 micrometastatic deposits did not give rise to additional deposits and were free to either die out, grow out to greater than 1 million cells, or become dormant. Representative traces of micrometastatic-tumor growth for 4 survival probabilities (0.500, 0.502, 0.504, and 0.506) are shown in Fig. 1B.

The metastatic fate across all simulated patient trials for 1, 2, 4, and 8 cell starting metastases revealed that the majority of metastases died out until the survival probability exceeded 60% (Fig. 2). Sampled time courses of metastatic tumor growth illustrate this phenomenon, even for probabilities above 50% (Fig. 1B). Eight-cell starting metastases conferred a grow-out advantage particularly above the 60% survival probability whereas 1-cell starting metastases still exhibited metastatic abatement in that range. Of those metastases that die out, they do so quickly in fewer than 60 cycles (Fig. 3). Dormant metastases represented a small fraction and were tightly bounded by a small survival probability.

Metastatic dormancy was a rare outcome that manifested only between 49.7% and 50.8% survival probability (Fig. 4, top). Starting with 1 versus 8 cells in the micrometastases did not impact this survival probability range. However, the mean number of dormant metastases increased as the starting cell number increased, up to a maximum of 80 mean dormant metastases per patient (still less than 10%) for 8-cell simulations (Fig. 4, bottom). In a small series of 1000- or 4000-cell micrometastases, the boundary percentages were within the same narrow survival probability range.

Although dormancy was tightly bounded proximal to the 50% survival probability, a dramatic variance of the mean dormant metastatic size was highly sensitive to small changes in survival probability (Fig. 5). Half of the survival range (between 49.7% and 50.2%) yielded dormant metastases composed of fewer than a few thousand cells. Interestingly, between 50.3% and 50.7% survival probabilities the populations peaked at around 500,000 cells. Not until the 50.8% survival probability did the dormant metastases (a single metastasis for 2-cell initialization, 3 metastases for the 8-cell case) approach the outgrowth condition—reaching nearly 1 million cells.

The cumulative density function in Fig. 6 compares 1-cell versus 8-cell starting metastases to the outgrowth condition. At the 70% survival probability, all of the metastases that grew out did so in less than 50 cycles. In contrast, half the metastases grew out in approximately 400 cycles for the 51.5% survival probability. Interestingly, at just above the 50% survival probability, a one-percentage point increase in survival probability resulted in the required number of cycles for the surviving metastases to achieve outgrowth being reduced by two thirds.

Research into metastatic latency has rapidly evolved over the past decade given the high prevalence of this condition and the lack of effective clinical interventions. Unfortunately, this phenomenon is difficult to recapitulate across in vitro and in vivo experimental systems. Moreover, it is not certain whether metastatic latency results from small deposits of quiescent carcinoma cells that had recently extravasated into the metastatic site or whether small- to mid-sized micrometastases maintain clinically undetectable sizes based on balanced proliferation and cell death (or a combination of explanations). This distinction is critical as therapeutic approaches would be different for either condition: that of targeting proliferating or quiescent cells.

Because of the metastatic latency prevalence, patients are confronted with the dilemma of undergoing extensive and toxic systemic treatments for putative metastases based on the characteristics of the primary nodule and statistical considerations of similar patients that experienced metastatic latency. Although the ability to predict which patients will suffer metastatic recurrence has improved (18), the current approach to these dreaded sequelae only reduces the rate by a third, while evoking untoward toxicity including hair loss, neurological impairments, opportunistic infections, bleeding, and death. As most of the therapy is ineffective, we need to characterize the status of the disseminated tumor cells before their frank metastatic emergence so as to more effectively target them. What is not known is whether these few surviving cells exist in cellular quiescence or undergo balanced proliferation to engage the dormant phenomenon. This distinction has immense implications, as predominantly only cycling cells are susceptible to chemotherapy, and is therefore the focus of this article.

Given these limitations, many in silico models have been developed that may help prioritize or frame the metastatic etiology (19, 20). For example, Michor and colleagues computed the probability of metastatic dissemination from primary tumors based on oncogenic mutations in RAS, ERBB2/NEU, or MYC (21). Implementing a stochastic Moran process the authors simulated a single genetic mutation that gave rise to neutral, advantageous, or disadvantageous metastatic potential. These simulations computed the temporal accumulation of metastases as a function of primary tumor size. Further model development eventually may help clinicians validate whether to treat more aggressively with adjuvant therapy. Another study implemented a Voronoi tessellation approach to model the impact heterogeneous or homogenous extracellular matrix (ECM) stiffness had on metastatic and primary tumor growth kinetics (9). Whereas nonrigid, homogenous ECM gave rise to isotropic tumors and rigid homogenous ECM invoked anisotropic morphologies, a heterogeneous ECM stiffness exhibited both properties, but in a predicable manner. Consequently this modeling if matured may assist clinicians in determining which metastatic tumor-adjacent tissue (in critical metastatic organs such as the brain) may likely harbor disseminated tumor cells—leading to either more or less aggressive tissue resection. Although these models have promise to improve clinical treatments there is a dearth of modeling within the very early stages of metastatic extravasation.

The model presented here invokes explicit representation of each cell and leverages repeated computed trials to sample the distributions of die out, grow out, and dormancy phenotypes exactly. However, analytic formulations called branching processes can also be applied to derive convenient predictions for this type of growth (22). These models exploit the fact that daughter cells are independent from one another and thus behave as if individual metastases (branching processes) themselves, lending to a compact, recursive formulation: |$h_n (t) = [f^n (t)]^N$|⁠. Here, h(t) is the “probability generating function” (pgf) at cell cycle n, and is computed as the product of N identical pgf's (for N starting cells in the micrometastatic tumor; ref. 23). The pgf allows us to quickly compute the expected cell population after n cell cycles |$\left({N\left[ {(2p)^n} \right]} \right)$|⁠, the expected variance |$\left({\frac{{\left({2p} \right)^{(n - 1)} \left({\left({2p} \right)^n - 1} \right)}}{{\left({2p} \right) - 1}}}{\left({4p - 4p^2} \right)} \right)$|⁠, and the probability of eventual extinction |$\left({\frac{p}{{\left({1 - p} \right)}}} \right)^N$|⁠. In all expressions, the inputs are the number of starting cells, N, the survival probability, p, and the number of cell cycles, n, just as in our stochastic simulation model. The recursive formula above cannot yield the exact population distribution for any large n (23), such as our cell cycle limit of 1,218. That is, although we can predict dying out and survival, we cannot exactly predict dormancy with the model, for which simulations are required.

Metastatic lesions are visible through imaging modalities such as MRI when they approach 0.1 cm³. A tumor of this small size actually consists of approximately 100 million cells (24). Therefore, disseminated tumor cell deposits of less than this cannot be detected without invasive biopsy and subsequent analyses such as immunohistochemistry. Because metastatic latency can be explained by 2 very different mechanisms—that of cellular quiescence or proliferation offset by cell death (or a combination of both), we sought to model the second hypothesis and determine what survival probability range would give rise to the dormant phenotype. We simulated 100 patients each with 1,000 metastatic deposits (1-, 2-, 4-, and 8-cell starting sizes) across stochastic survival probabilities discontinuous between 30% and 70%. The modeling revealed a very narrow survival probability range between 49.7% and 50.8% that gave rise to dormant metastases. The reason why this dormancy probability window is so small is based upon how sensitive simulated metastases are to small changes in survival probability. If the probability of survival is less than 50%, the metastases all die out with the exception of the few dormant clusters because the probability favors the boundary condition of dying out. Only a rare escape from this tendency will therefore evoke a dormant metastasis. Once the probability of survival exceeds 50.8%, the grow-out condition of reaching 1 million cells begins to dominate. Thus, we conclude that it is unlikely that metastatic dormancy is governed predominantly by internally balanced proliferation and cell death. As shown in Fig. 3, it may seem unexpected that micrometastases will die out even at survival probabilities exceeding 70%. The reason for this phenomenon is closely related to the number of starting cells in the metastases. It is difficult for a simulated metastasis to live through its first few generations especially when it has 2 starting cells; because when a cell dies the metastasis evolves closer to the boundary of dying out. As Fig. 2 shows, the chance of a metastasis dying out is reduced when there are 8 starting cells versus 2 starting cells. Of the metastases that die out, they do so in a short time (Fig. 3) as this relates to the fact that the metastases start near the null boundary and thus are easy to involute early before attaining a greater cell number. An often suggested model involves micrometastases growing and then being constrained at a subclinical size by the need for external supports; this is the basis of the postulated “angiogenic switch” (25). However, even here the constraints must be tightly controlled; in an initial series of examining metastases of 1,000 or 4,000 starting cells, the range for staying dormant remained constrained similarly to the smaller micrometastases at between 49.5% and 51% survival probabilities; if the survival percentage drops below this, the metastatic nodules quickly die out (data not shown). Thus, even this tumor extrinsic constraint would need to be tightly controlled or frequently vacillate between outgrowth and involution. In the spirit of Occam's razor, we propose that extended periods of cellular quiescence most likely accounts for this dormant period.

Although there are numerous models for dormancy these suffer from significant differences in the animal hosts, such as life-span and time-scale to development, and a lack of information about the micrometastatic state in actual cancer patients. This approach does not account for the richness of interactions and signals that may constrain or support micrometastatic nodules. For instance, it is possible that an “angiogenic switch” is required for outgrowth once these nodules reach a certain size. However, the experimental data are lacking currently to build models incorporating such inputs. Rather, such models would presuppose these constraints and be built to account for the limitations and interventions (26). Herein, we chose an a priori unbiased calculation to simply define the boundary conditions of metastatic outgrowth regardless of the actual biological and physiological networks.

The value of mathematical models lies in their ability to suggest avenues for investigation. Recent findings of tumor cells in circulation or bone marrow despite no evidence of metastases may provide some insights into the behavior of the early metastatic cells. Interestingly, ongoing proliferation is evident only in few of these cells (27) consistent with quiescence of these cells. If these cells were truly quiescent then therapies targeted to killing growing cells would not be effective and may even cause outgrowth if the dormant microenvironment is perturbed. Still such findings are only suggestive, and highlight the need to examine human micrometastases with innovative technologies to determine the actual proliferative rates in these hidden metastases. As our currently available chemotherapy mainly attacks cycling cells, the mode of quiescence would make this treatment less effective and necessitate new approaches to the problem of metastatic dormancy.

No potential conflicts of interest were disclosed.

Conception and design: D.P. Taylor, J.Z. Wells, A. Soval, A. Wells

Development of methodology: D.P. Taylor, A. Soval

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D.P. Taylor, J.Z. Wells, A. Soval

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D.P. Taylor, J.Z. Wells, A. Soval, A. Wells

Writing, review, and/or revision of the manuscript: D.P. Taylor, J.Z. Wells, A. Soval, A. Wells

Study supervision: A. Wells

The authors thank members of the Wells and Bahar laboratories for helpful discussions.

This work was supported by grants from the VA Merit Award program and NIH (UH2TR000496) to A. Wells.

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.