Effects of Surgery and Chemotherapy on Metastatic Progression of Prostate Cancer: Evidence from the Natural History of the Disease Reconstructed through Mathematical Modeling

This article brings mathematical modeling to bear on the reconstruction of the natural history of prostate cancer and assessment of the effects of treatment on metastatic progression. We present a comprehensive, entirely mechanistic mathematical model of cancer progression accounting for primary tumor latency, shedding of metastases, their dormancy and growth at secondary sites. Parameters of the model were estimated from the following data collected from 12 prostate cancer patients: (1) age and volume of the primary tumor at presentation; and (2) volumes of detectable bone metastases surveyed at a later time. This allowed us to estimate, for each patient, the age at cancer onset and inception of the first metastasis, the expected metastasis latency time and the rates of growth of the primary tumor and metastases before and after the start of treatment. We found that for all patients: (1) inception of the first metastasis occurred when the primary tumor was undetectable; (2) inception of all or most of the surveyed metastases occurred before the start of treatment; (3) the rate of metastasis shedding is essentially constant in time regardless of the size of the primary tumor and so it is only marginally affected by treatment; and most importantly, (4) surgery, chemotherapy and possibly radiation bring about a dramatic increase (by dozens or hundred times for most patients) in the average rate of growth of metastases. Our analysis supports the notion of metastasis dormancy and the existence of prostate cancer stem cells. The model is applicable to all metastatic solid cancers, and our conclusions agree well with the results of a similar analysis based on a simpler model applied to a case of metastatic breast cancer.


Two Paradigms of Cancer
According to the conventional paradigm, cancer emerges when one or a few adjacent cells acquire a number of irreversible oncogenic mutations which eventually perturb cell cycle controls and apoptotic regulation. Subsequent proliferation of the initial transformed cell(s) results in a malignant tumor that develops a capillary network and acquires the ability to invade surrounding tissues and metastasize. Within this framework, cancer is viewed as an alien entity that progresses sequentially through stages characterized by the extent of its anatomic spread -local, regional and distant. Metastases are considered independently growing tumors that arise from malignant cells shed by the primary tumor and seeded at various secondary sites. An upshot of this view is that cancer treatment should consist of eliminating all cancer cells and the earlier and more aggressive the treatment of the primary tumor the better the prognosis.
An alternative paradigm of cancer has recently started to crystallize on the basis of more than 100 years of extensive clinical observations, epidemiological studies and animal experiments (see [1][2][3][4] for a comprehensive review). According to this new paradigm, a malignant tumor is an organ-like entity that exists in a dynamic state of homeostasis with its microenvironment which can act either as a promoter or suppressor of cancer progression. An important aspect of such a homeostasis is the state of dormancy of small avascular primary and secondary tumors and solitary cancer cells. In particular, large numbers of cancer cells may circulate in the blood and lymph vessels. The state of dormancy is maintained through the balance between the factors of growth and angiogenesis, on the one hand, and inhibitors of these processes, on the other. Seeding of metastases may occur long before primary cancer manifests clinically, which makes cancer a systemic disease already at its very early subclinical stages. Primary and secondary tumors at various sites engage in a complex biochemical interaction that typically results in suppression of the growth of small tumors by larger tumors. Accordingly, resection of the primary tumor will weaken this inhibitory effect. Additionally, wound healing processes following resection of the primary may promote growth of secondary tumors and angiogenesis. This leads to disruption of the dormant state of metastases and causes their accelerated growth and vascularization. Therefore, surgery may not always be a good treatment option. Maintaining the state of homeostasis or dormancy may prove a better strategy of cancer control.
The goal of the present work is to confirm or confute, in the case of metastatic prostate cancer, the principal biomedical hypotheses that lie behind this alternative paradigm of cancer and to extend them to chemotherapy (see below). Three of these hypotheses are related to the natural history of the disease (group A) and the other three to the effects of surgery and chemotherapy on metastatic progression (group B). Our findings will show that, contrary to the commonly accepted views, these two modes of treatment of the primary tumor have only minor effect on the rate of metastasis shedding but have a dramatic accelerating effect on the rate of metastatic growth. We believe the same is true for radiotherapy. As yet another outcome, our analysis lends further support, if only indirect, to the notions of tumor dormancy and cancer stem cells.
Testing of the hypotheses will be accomplished through reconstruction of the individual natural history of cancer and the effects of its treatment on the basis of a comprehensive mathematical model of cancer progression developed in [5][6][7] and extended in this work. Parameters of the model are estimated from the data on volumes of primary tumor and metastases for a cohort of 12 prostate cancer patients diagnosed and treated at the Memorial Sloan-Kettering Cancer Center (MSKCC).

Working Hypotheses
A1. Metastatic dissemination off the primary tumor is an early event in the natural history of the disease that may occur long before primary tumor becomes clinically detectable.
A2. Prior to the start of irreversible proliferation in a secondary site, micrometastases or solitary cancer cells spend an extended period of time in a state of dormancy or free circulation.
A3. The primary tumor has a small subpopulation of "cancer stem cells" of relatively constant size characterized by self-renewal, capacity for fast proliferation and high metastatic potential. B1. Treatment of the primary tumor has only limited effect on the process of metastasis shedding. B2. Extirpation of the primary tumor by surgery may boost the proliferation of dormant or slowly growing metastases, trigger their vascularization and accelerate growth of vascular secondary tumors.
B3. Chemotherapy may also accelerate the growth of metastases, although not necessarily by the same mechanism.
These hypotheses are mainly focused on metastatic progression because metastases are responsible for about 90% of cancer-related deaths. Below we discuss the status of these hypotheses, briefly review supporting biomedical evidence and discuss underlying biological mechanisms. For a more extensive discussion, the reader is referred to the reviews [1-4].

Hypotheses The Whys and Wherefores
A1. This hypothesis has been discussed in the medical literature for several decades [8][9][10]. It was estimated in [10] that more than 70% of cancer patients have occult metastases at presentation. Because hypothesis A1 deals with unobservable events, its direct confirmation may prove elusive. However, a wealth of indirect evidence points to its validity. For example, how else could one explain the fact that a significant fraction of patients whose primary tumor was diagnosed and removed at the earliest stages of cancer progression still develop distant metastases?
A2. The earliest report on circulating cancer cells goes back to 1869 [11] (see also [12]). Numerous modern studies of various types of cancer [13][14][15][16] showed that significant quantities of circulating tumor cells can be present in blood and lymph channels without clinically manifest metastases. The next critical step in the multi-stage process of metastasis formation is invading a host site and establishing there. High sensitivity of tumor cells to the conditions of host microenvironment and the resulting selective affinity of tumor cells to specific organs and tissues has been known since 1889 under the name of "seed and soil hypothesis" [17]. However, it was uncovered by subsequent pioneering studies [18,19] that even if seeded into a fertile "soil," a tumor cell may remain dormant in a secondary site yet retain its clonogenic capacity. Direct evidence of prevalent breast tumor cell dormancy was presented in [20]. The technique of in vivo video microscopy enabled direct observation and quantitative study of dormant cancer cells [21,22]. Recently, tumor cell dormancy was recognized as a major new direction in cancer research [23]. It was found that about 1/3 of breast cancer patients 7-22 years after mastectomy and without any evidence of the disease had circulating tumor cells [24]. Their presence in peripheral blood was also confirmed in 24% of prostate cancer patients before prostatectomy [25] as well as in patients with undetectable PSA levels more than five years postprostatectomy [25,26].
A balance between proliferation, apoptosis and dormancy of cancer cells brings about the possibility that primary or secondary cancer remains subclinical for an extended period of time; as an example, breast cancer recurrence was reported to occur after 20 to 25 years of disease-free period [27].
The last critical step on the pathway leading to a detectable metastasis is induction of angiogenesis [28]. It was estimated that only about 4-10% of actively growing metastases eventually develop the capillary network which enables further growth [29].
The time period between shedding of a metastatic cell by the primary tumor and the beginning of its irreversible proliferation in a host site resulting in a clinically detectable secondary tumor will be referred to as metastasis latency. The latency period comprises the stages of free circulation, dormancy or slow avascular growth in a secondary site and induction of angiogenesis. In what follows, we will estimate the expected metastasis latency times for 12 prostate cancer patients.
A3. "Cancer stem cells" were first discovered in the case of acute myeloid leukemia [30]. Later their existence was confirmed for other types of hematologic cancers as well as for several solid cancers, including prostate cancer [31]. The hallmarks of cancer stem cells are self-renewal, pluripotency (i.e., the ability to produce various types of cancer cells through differentiation), fast proliferation and high metastatic potential. The existence of cancer stem cells is highly consequential for cancer treatment, for it implies that targeting and eliminating, or at least controlling, a very small subpopulation of tumor cells is critical for the success of cancer therapy.
B1. Arguments in support of the hypothesis that local control may have only a limited effect on the probability and timing of distant failure have been so far mostly of three kinds: (a) assessing whether local and distant failure are statistically correlated events and whether the same clinical variables and risk factors predict for both of them; (b) observing various relationships between the age at local recurrence and distant failure (for example, patients who fail locally may display an increase in the hazard rate of the time from treatment to detection of metastases as well as larger number and volumes of metastases, as compared to locally-controlled patients); and (c) in the case of radiotherapy, examining the effects of dose escalation on cancer-specific survival. What these phenomenological considerations tend to neglect is heterogeneity of biological mechanisms underlying the effects of various modes of treatment of the primary tumor (such as surgery, external beam radiation, brachytherapy, and chemo-and hormonal therapy) on metastases. The probability and age at distant failure depend on four important characteristics of metastasis: (1) the rate of metastasis shedding by the primary tumor; (2) the fraction of metastases shed by the primary tumor that may potentially give rise to detectable secondary tumors in a given site; (3) the duration of metastatic latency; and (4) the sitespecific rates of growth of metastases. The structure of the model employed in this work and the scarcity of data available for estimation of model parameters allowed us to study the effects of treatment of the primary on only two of these four characteristics -the rates of metastasis shedding (hypothesis B1) and growth (hypotheses B2 and B3). The model can be extended to include the effects of treatment of the primary on the duration of metastatic latency at the cost of introducing additional model parameters. However, parameter estimation for such an extended model would require much larger sample sizes or longitudinal data.
B2. Cancer patients who present even at late stages of the disease rarely have clinically manifest metastases; typically, they surface after the start of treatment. That this phenomenon is deeply rooted in basic cancer biology is postulated in hypothesis B2. This hypothesis addresses one of the multiple effects that a primary tumor exerts on other primary or secondary tumors. Experimental studies of these effects on animal models were conducted as early as the beginning of the 20 th century [32][33][34][35]. Numerous later works confirmed these early findings although changed their interpretation, see [2,3] and references therein. The most important discovery within this realm of research is that larger tumors inhibit growth of smaller ones and, as a result, resection of large primary or secondary tumors accelerates the growth of smaller tumors. In particular, it was found that extirpation of the primary tumor triggers aggressive proliferation of dormant or slowly growing metastases and their vascularization [36][37][38]. This important finding was further supported by a number of clinical case studies including eight cases of non-seminomatous germ-cell testicular cancer [39] and three cases of melanoma [40,41]. Additionally, and most importantly for the present study, substantial evidence of post-resection progression of metastatic disease in prostate cancer patients was presented in [42]. The validity of hypothesis B2 is supported by epidemiological analyses of the time course of post-surgery recurrence for various categories of cancer patients [43][44][45], propped by a simple mathematical model of breast cancer progression [46] and corroborated by an experimental study of accelerated growth of metastases after resection of the primary Lewis lung carcinoma in mice [47]. The hypothesis was also employed to explain the "mammography paradox", viz. an unexpected increase in the mortality of premenopausal node-positive women aged 40-49 diagnosed with breast cancer as a result of screeningbased early detection in a number of large-scale clinical trials [48][49]. The extent of accelerated growth of metastases was found to be proportional to the extent of surgery. As one example, tumor recurrence in non-metastatic colon cancer patients was markedly lower in a group that had laparoscopic surgery than for the open colectomy group [50]. Finally, even biopsy was reported to result in a measurable increase in the incidence of lung metastases in mice [51].
What is the mechanism of accelerated growth of metastases following resection (and possibly radiation treatment) of the primary tumor? Briefly, as hypothesized in [52] and confirmed in a host of other studies, many of which are reviewed in [1][2][3][4], the primary tumor and its microenvironment produce tumor growth factors as well as growth inhibitors. The growth factors are more easily degradable than growth inhibitors and propagate mostly by diffusion thus acting locally and promoting primary tumor growth. By contrast, growth inhibitors are more stable; as a result, they may reach remote secondary sites and impede the growth of metastases. The latter is also limited by various circulating angiogenesis inhibitors (such as angiostatin and endostatin). Removal of the primary tumor reduces production of growth inhibitors, which accelerates the growth of metastases. Additionally, and perhaps more importantly, healing of surgery-or radiation-related injury is accompanied by a surge in the local and systemic production of various growth and angiogenesis factors that act synergistically with the decrease in the levels of growth and angiogenesis inhibitors. This mechanism suggests that the same metastasis enhancing effect will also manifest for surgery or wounding unrelated to primary tumor. A statistical analysis based on 418 patients with advanced cancer reported in [53] confirmed that this indeed is the case. Similar but weaker effects may result even from biopsy [51]. The mechanisms described above seem to be also relevant to radiation therapy. B3. A significant fraction of cancer patients treated with chemotherapeutic agents develops resistance to treatment. Undoubtedly, this effect is due to many mechanisms; one of them is selection of resistant cells in the target population. This process is mediated and amplified by the formation of spontaneous and chemotherapy-induced mutations, adaptive reactions of cancer cells causing them to evade cytotoxic action of drugs by switching to alternative metabolic and proliferative pathways, and removal of cytotoxic agents from cancer cells by transporting them across cell membrane due to the action of ATP-binding cassette transporters. Finally, chemotherapy, as any other treatment, confers survival advantage on faster proliferating cells, unless the latter are more sensitive to the treatment than slower growing cells. Selection of resistant and fast proliferating cells leads to the decreased efficiency and ultimate failure of chemotherapy.
Taken collectively, hypotheses A1, B2 and B3, if confirmed, would suggest that there is no such thing as "local treatment" per se: any intervention aimed at the primary affects metastatic progression. Therefore, in the present work, we will use the term "local treatment" only as an indication of intent rather than statement about the mechanism or outcome.

General Methodology
The above six hypotheses are formulated in terms of events and processes that are typically unobservable, or only partially observable, in vivo. In the present work, we bring mathematical modeling to bear on estimation of their probabilities, timing and rate characteristics. A distinct advantage of this approach is that mathematical models make it possible to relate unobservable quantities such as the age at disease onset or the start of its metastatic dissemination to clinical variables that are recorded months, years or even decades after the occurrence of the initiating micro-events.
The individual natural history of cancer consists of two important (but unobservable) pieces of information: (1) time to (or age at) critical micro-events such as the emergence of the first malignant clonogenic cell, shedding of metastases by the primary tumor, their seeding at various secondary sites and the start of their irreversible proliferation (termed inception); and (2) rate parameters descriptive of various cancer progression processes including growth of the primary tumor, shedding of metastases, and their seeding and growth at various secondary sites. The patient's observable clinical variables include age, stage and primary tumor volume at diagnosis, various biochemical or genomic markers, and clinical data resulting from follow-up studies such as the age at tumor recurrence or cancer related death and, for our purpose, site-specific number and volumes of detectable metastases. Establishing a quantitative relationship between parameters of the natural history of cancer and remote clinical endpoints and response variables is the apanage of mathematical modeling. Finally, because cancer progression involves substantial variability of the characteristics of primary and secondary tumors and their microenvironments and is impelled through a number of sporadically occurring random events [54], a stochastic approach would be the modeling tool of choice.
In this work, we apply a comprehensive stochastic model of cancer progression developed in [5][6][7] to metastatic prostate cancer and extend the mathematical formalism to the case where patients receive systemic treatment alone. The basic idea behind the model is to view the process of metastasis shedding by the primary tumor as a Poisson process whose intensity is proportional to a certain power of the size of the primary tumor [6,7], see also [5]. The model allows for arbitrary laws of growth of the primary tumor and metastases, and leads to an explicit formula for the conditional distribution of the volumes of detectable metastases in a certain secondary site, given their number, at any time point with an intact, excised or regrowing primary tumor [6]. The basic model developed in [6] was extended to accommodate distinct rates of growth of metastases prior to and after excision the primary tumor [7,55]. In the present work, we further extend it to incorporate two distinct laws of growth of the primary tumor: before and after the start of treatment. In the case of non-recurring, surgically removed primary tumor the model reduces to the one developed in [7]. A parametric version of the extended model based on exponential growth of the primary tumor and metastases, and exponentially distributed metastasis latency times is developed in Section 5. Knowledge of identifiable model parameters allows one to estimate the most important temporal and rate characteristics of the natural history of metastatic cancer and the effects of treatment. Surprisingly (and reassuringly), the model developed in the present work is sensitive enough to correctly predict whether a given patient had surgery.
The model designed in [7] was applied to a breast cancer patient who developed, by age 82 and 8 years after diagnosis and resection of the primary, 31 detectable bone metastases [55]. The model provided an excellent fit to the empirical distribution of the observed volumes of metastases and led to the following patient-specific conclusions [55]: (1) The onset of the disease occurred at age 42, about 32 years prior to primary diagnosis.
(2) Inception of the first metastasis occurred at age 44.5, that is, about 29.5 years prior to the primary diagnosis and 2.5 years after the onset of the disease at which time the primary tumor was extremely small and certainly undetectable. (3) Inception of all detected metastases except one occurred before excision of the primary. (4) The expected metastasis latency time was about 79.5 years (which means that at the time of surveying most metastases were still dormant and undetectable). (5) Resection of the primary tumor was followed by a 32-fold increase in the rate of growth of bone metastases notwithstanding the fact that after surgery the patient was put on tamoxifen that suppresses growth of metastases and has anti-angiogenic effect.
(6) The process of metastasis shedding was essentially homogeneous (i.e., independent of the volume of the primary tumor), which suggests the presence within the tumor of a small selfrenewing subpopulation of relatively constant size consisting of cells with high metastatic potential. This may serve as indirect evidence for the existence of breast cancer stem cells.
In what follows we examine the applicability of these conclusions to prostate cancer patients.

A Mathematical Model of the Individual Natural History of Metastatic Cancer
The processes of metastasis formation, growth and progression are complex, heterogeneous and selective [10,56,57]. To form a micro-metastasis in an organ or tissue, a malignant cell has to detach itself from the primary tumor, degrade extracellular matrix, intravasate, traverse the circulatory network, evade attacks by the immune system, extravasate, invade a secondary site, survive through the dormancy period, start to proliferate and induce angiogenesis. As a result of this multi-stage selection process, only a tiny fraction of cells shed by the primary tumor give rise to actively growing metastases [10,56,57].
The temporal natural history of metastatic cancer is commonly divided into three overlapping periods: disease-free period, primary tumor growth and metastatic progression. These periods and relevant model assumptions are described below and illustrated in Figure 1. Disease-free period begins with the birth of an individual (or start of exposure to a carcinogen) and ends with the appearance of the first malignant clonogenic cell, a development termed the onset of the disease.
Primary tumor dynamics. The size of the primary tumor (that is, the total number of tumor cells) at any time t counted from the age T of disease onset will be denoted by Φ(t). Prior to the start of treatment, the growth of the primary tumor is governed by a function Φ 0 and thereafter by another function Φ 1 , which acts multiplicatively on the size of the primary tumor at the start of treatment (at age V). The function Φ 0 is strictly increasing, continuous and satisfies the initial condition Φ 0 (0) = 1. As to the function Φ 1 , it is continuous but not necessarily increasing. In particular, for a non-recurrent excised tumor, Φ 1 = 0. Functions Φ 0 and Φ 1 may depend on one or several parameters. We denote by φ the inverse function for Φ 0 . It follows from the above assumptions that Metastasis formation. The process of metastasis shedding is governed by a Poisson process with rate μ proportional to the number, N(t), of metastasis-producing cells at time t: μ(t) = α 0 N(t), where α 0 > 0 is the rate of metastasis shedding per cell. Because N(t) is unobservable, we relate it to the primary tumor size Φ(t) through the formula N(t) = α 1 Φ θ (t) with some constants α 1 > 0 and θ ≥ 0. The value θ = 1 means that a constant fraction of cells in a tumor have metastatic potential. It is known that many solid tumors enclose a core of hypoxic, clonogenically sterile cells or even a broth of proteins, while actively proliferating clonogenic cells are concentrated near the tumor surface; in this case one would expect θ = 2/3. Finally, the case θ = 0 corresponds to the existence of a relatively stable, self-renewing subpopulation of metastasis-producing cells within the primary tumor. In summary, the rate of metastasis shedding is: where α = α 0 α 1 . In the case θ = 0 the rate of metastasis shedding μ is constant and the underlying Poisson process is homogeneous.
It is further assumed that metastases shed by the primary tumor give rise to clinically detectable secondary tumors in a given site independently of each other and with the same probability q. Therefore [58, pp. 257-259], inception of metastases in the site in question is governed by a Poisson process with intensity ν = q μ. Each viable metastasis is assumed to spend some random latency time between detachment from the primary tumor and its inception in a secondary site. We assume that latency times for different viable metastases are independent and identically distributed with some probability density function (pdf) f and corresponding cumulative distribution function (cdf) F. Then, see e.g. [59], the resulting process of metastasis inception is again a Poisson process with the rate Q Timeline of the natural history of metastatic cancer and observables. Suppose that the observed primary tumor size at age V is S. Then the patient's age T at the disease onset is given by the formula We will assume that local or systemic treatment was given (or started) at age V, and that at age W, W > V, a certain number, n, of metastases were detected in the same secondary site with the observed volumes X 1 , X 2 ,..., X n , where X 1 < X 2 < ... < X n . Thus, 0 < T < V < W (Figure 1).
Growth of metastases. Prior to the start of treatment, the growth of the size of any viable metastasis in a given secondary site is governed by a function Ψ 0 , while during or after the treatment, the size of the metastasis is growing according to a potentially different function Ψ 1 , which acts multiplicatively on the size of the metastasis at the start of treatment. We assume for simplicity that actively growing metastases start from a single cell. Functions Ψ 0 , Ψ 1 are strictly increasing, differentiable, and satisfying the initial condition Ψ 0 (0) = Ψ 1 (0) = 1. Additionally, they may depend on one or several parameters. It follows from our assumptions that the size Ψ (y) of a viable metastasis at time y from inception is given by: This function is strictly increasing, continuous, piecewise differentiable and satisfies the condition Ψ(0) = 1.
Secondary metastasis. Secondary metastasizing (that is, formation of "metastasis of metastasis") to a given site, both from other sites and from within, is assumed negligible.
Metastasis detection. The volume of a metastasis becomes measurable when it reaches some threshold value m. This value and the accuracy of volume measurement are determined by the sensitivity of imaging technology. In case of PET/CT imaging involved in this study, m = 0.5 cm 3 , and the accuracy of volume determination is one voxel, i.e., approximately 0.065 cm 3 .
Effects of treatment. Because the rate of secondary metastasizing is assumed negligible, the formation of new metastases is stopped at the time of resection of a non-recurrent primary tumor. Any mode of local or systemic treatment (surgery, radiation, chemo-or hormonal therapy) is assumed to affect metastases after their inception in a given secondary site only through the rate of their growth (and not through prolongation of their latency times).

Distribution of the Sizes of Detectable Metastases
Let X be the size of a detectable metastasis with inception time Y (relative to the onset of the disease) that was surveyed at age W. Then: where W -T -Y is the metastasis progression time from inception to detection and function Ψ is given by (3). The maximum value of X is: and the inverse function, ψ:= Ψ −1 , is given by: The distribution of the sizes of metastases in a given secondary site is specified in the following theorem [6,7,55]. Theorem 1. The sizes X 1 < X 2 < ... < X n of metastases in a given secondary site that are detectable at age W are equidistributed, given their number n, with the vector of order statistics for a random sample of size n drawn from the distribution with the following pdf: and p(x) = 0 for x [m, M], where the tumor onset time T is given by (2), function ψ is defined in (5), M is specified in (4), and F is the cdf of the metastasis latency time corresponding to the pdf f. For a proof of Theorem 1, see [6]. Notice that the distribution given by Equations (6)-(7) is free of parameters α, q and sample size n. Observe also that for θ = 0 we have: In this case, the distribution p(x) is independent of the laws of primary tumor dynamics before and after the start of treatment. Setting in Equations (6)-(8) m = c, the volume of a single cancer cell, we obtain the distribution of the sizes of all (both occult and detectable) metastases in a given site. Finally, the site-specific total number of viable metastases at age t > T is Poisson distributed with parameter (expected value) while the probability of developing viable metastases at age t is 1-exp{-E(t)}. In particular, for a tumor resected at age V we have: Due to the non-stationarity of the process of metastasis seeding and the lack of a "natural" order for listing detectable metastases, the sizes (or volumes) of metastases detectable in a certain secondary site at a given time do not form a random sample from a probability distribution. However, it follows from Theorem 1 that the distribution of any rearrangement-invariant statistic based on observations X 1 , X 2 ,..., X n would be identical to the distribution of the same statistic based on a random sample of size n drawn from the pdf p given by Equations (6)- (7). In particular, the joint likelihood of the observations X 1 , X 2 ,..., X n , where X 1 < X 2 < ... < X n , given by the formula L(X 1 , X 2 ,..., X n ) = n! n i i X p 1 ) ( has the same form (apart from the factor n!) it would take should the observations X 1 , X 2 , ..., X n form a random sample from the distribution with pdf p. Therefore, identifiable parameters of a suitably parameterized model of the natural history of metastatic cancer described in Section 3 can be estimated using the method of maximum likelihood.

Distribution of the Sizes of Detectable Metastases for Exponentially Growing Primary Tumor and Metastases
In this section, we introduce a parametric version of the general model of cancer natural history described in Section 3 and compute the distribution p(x) underlying the site-specific sizes of detectable metastases given by Equations (6)- (7).
Suppose that the size of the primary tumor grows exponentially with constant rate β 0 > 0 before treatment and with rate β 1 after the start of treatment: Φ 0 (t) = exp{β 0 t} , 0 ≤ t ≤ V -T, where time t is counted from the age T of tumor onset, and Φ 1 (t) = exp{β 1 t}, where time t is counted from the start of treatment. Note that rate β 1 can be negative. Then for φ, the inverse function for Φ 0 , we have φ(s) = (ln S)/β 0 , so that T = V -(ln S)/β 0 , see Equation (2). Clearly, we must have T > 0, which implies β 0 > V S ln (11) We will assume that before the start of treatment metastases in the site of interest grow exponentially with rate γ 0 > 0, so that Ψ 0 (t) = exp{γ 0 t}. Note that for all 12 patients analyzed in this work their metastases reached considerable sizes at the time of surveying. That is why we are assuming that after the start of treatment the sizes of metastases also grow exponentially with rate γ 1 > 0. Then Ψ 1 (t) = exp{γ 1 t}, and Equation (3) Suppose additionally that metastasis latency times are exponentially distributed with the expected value ρ: f(s) = ρ −1 e −s/ρ , s > 0.
The resulting parametric model of cancer natural history depends on the following 8 parameters: β 0 (the rate of growth of the primary tumor prior to treatment), β 1 (the rate of growth of the primary tumor after the start of treatment), α and θ (two parameters involved in the expression (1) for the rate of metastasis shedding), q (the probability that a metastases shed by the primary tumor will evolve into a viable, potentially detectable secondary tumor), γ 0 (the rate of growth of metastases in the presence of untreated primary tumor), γ 1 (the rate of growth of metastases after the start of treatment) and ρ (the mean metastasis latency time). Recall, however, that the distribution of the site-specific sizes of metastases depends only on 6 parameters: β 0 , β 1 , θ, γ 0 , γ 1 and ρ.
Introduce the following alternative set of 6 model parameters: Note that 0 < A < M.

The Case of Surgery with Non-Recurrent Primary Tumor
The case where the primary tumor was resected at age V and did not recur by age W was considered in [7,55]. It arises from a more general model described in Section 3 if one sets Φ 1 = 0. Accordingly, the corresponding 5-parameteric version of the general 6-parametric parametric model obtains by setting β 1 → − ∞. The pdf, p(x), underlying the site-specific distribution of the sizes of metastases at age W is given by the following expressions computed on the basis of Equations (6)-(7) [7]: (1) If m A d , then ( Recall also that p(x) = 0 for x < m or x > M. The model (13)-(16) will be called the Surgery model. We will also consider a limiting case of the above parametric model where θ = 0. Here the Poisson process of metastasis shedding is homogeneous. Accordingly, this model will be termed the Homogeneous model. Letting a 0 → 0 in the Surgery model (13)-(16) we have: (1) If m A d , then (17) where The Surgery model (13) [7], in the case A > m all parameters of both models are identifiable (for an at length discussion of identifiability of stochastic models, see [60]).

The Case of Non-Resected or Resected Recurrent Tumor
In this case the distribution p(x) obtained from Equations (6)-(7) has the following form.
The Full model (13), (14), (21) In what follows, quantities x, m, A and M will be expressed as volumes assuming the average volume, c, of a cancer cell to be 10 −9 cm 3 . Observe, however, that because the function xp(x) in all cases depends only on the ratios of metastasis sizes x, m, A and M, the likelihood maximizing parameters are independent of c.

Computing Biological Parameters
The Full 6-parametric model of the natural history of cancer is determined by the biological parameters β 0 , β 1 , θ, γ 0 , γ 1 and ρ. Equations (12) enable their expression through model parameters A, M, a 0 , a 1 , b 0 , b 1 . First, observe that Next, the expression for parameter M allows us to compute the disease onset time: (25) In view of the inequality T > 0, model parameters should satisfy the following constraint: Computing the other three biological parameters requires the knowledge of the primary tumor size, S, at age V:  (27) Because the volume, S v , of the primary tumor was estimated by pathologists based on a rough estimate of tumor margins, determination of the primary tumor size S = S v /c, where c is the volume of a single cancer cell, involves a considerable error. Yet another source of error is our assumption that c =10 −9 cm 3 . Furthermore, for eight patients the data on the volume of the primary was unavailable and was ascribed a value of 20 cm 3 , see Table 1. However, these errors result in only minor deviation in the values of biological parameters β 0 , β 1 , θ due to the fact that S appears in the formulas for these parameters under the sign of logarithm. Note that parameters γ 0 , γ 1 , ρ and the age at disease onset T are independent of S.  (27), depend on the logarithm of S and, as a result, vary only slightly with S.

The Data Set
To estimate model parameters, we used a data base of prostate cancer patients diagnosed and treated at MSKCC. To be useful for our analysis, the patients had to satisfy the following requirements: (1) availability of whole body PET/CT scans; (2) the number of metastases in a single secondary site (e.g., the skeletal system) is large enough ( 10); and (3) W > V, where V is the age at which the volume of the primary was measured immediately prior to surgery and/or start of systemic treatment while W is the age at metastasis surveying. Only 12 patients in the data base satisfied these conditions. Information on their gross primary tumor volume was generally unavailable and, when obtainable, quite likely affected by substantial errors. Among the 12 patients, one had surgery (radical prostatectomy) and one received surgery and external beam radiotherapy. Additionally, all 12 patients were given complex combinations and time courses of chemotherapy and adjuvant hormonal therapy. Relevant clinical variables for the 12 patients are given in Table 1. The number, n, of bone metastases detected in these patients varied between 10 and 58. When applying the above model of cancer progression we assumed that parameters β 1 and γ 1 represented the average rates of growth of the primary tumor and bone metastases, respectively, over the entire period from the start of treatment (age V) to metastasis survey (age W).

Results
Parameters A, M, a 0 , a 1 , b 0 , b 1 of the Full model were estimated by maximizing the likelihood function L under constraints (11) and (26). In all cases, the optimal value of parameter A satisfied the condition A > m which prevented the model from degeneration, see Section 5. Discontinuity of the likelihood as function of A at the observed volumes of metastases and the presence of multiple local maxima prohibited the use of gradient methods for likelihood maximization. Instead, we used a genetic global optimization algorithm (Differential Evolution) built into Mathematica package.
Optimal parameter values of the 6-parametric Full model along with the minimal values of -(log L)/n, where n is the number of bone metastases observed in a given patient, are presented in Table 2 while optimal values of biological parameters are given in Table 3. The optimized model provided an excellent fit to the empirical distribution of the volumes of detectable metastases for all patients; see Figure 2 where empirical and theoretical distribution functions are presented for one particular patient. The corresponding pdf, p(x), for this patient is shown in Figure 3. For all patients, the optimal value of parameter A was equal to (more exactly, approached from the right) one of the observed volumes of metastases. The values of parameter θ were small for all patients. Therefore, we also applied the corresponding 4-parametric Homogeneous model (θ = 0), see Section 5. For all patients, the estimates of biological parameters γ 0 , γ 1 , ρ and T of the Homogeneous model and the likelihood agree reasonably well with (and for patients 3 and 9 are very close to) those for the Full model, see Tables 3 and 4. Along with the age at onset T computed through Equation (25), Tables 3 and 4 also give the age at inception of the first bone metastasis (which, according to the model, is the metastasis with the largest observed volume). At that age, primary tumors of all patients were extremely small so that application of any deterministic law of tumor growth for their estimation would be misleading.
As a self-consistency check, we applied the non-surgery model with two growth rates of the primary tumor (β 0 before and β 1 after the start of treatment) to the two surgery patients (patients 1, 2). As expected, the estimated values of β 1 were small negative numbers: β 1 = −2.     3 where p(A+)/p(A-) = γ 1 /γ 0 = 3.5. Parameter A represents the volume of a metastasis whose inception occurred at the time of surgery while M is the maximum possible volume of bone metastases, that is, the volume of a metastasis whose inception occurred at onset of the primary tumor.

Summary of Results
The results presented in Tables 2-4 lead to the following tentative conclusions about the natural history of metastatic prostate cancer and the effects of its treatment, in good agreement with a similar analysis of a breast cancer patient [55], see Section 2.

Onset of the Primary Tumor
The age at cancer onset displayed a substantial inter-patient variability. For several patients, the disease started in the early childhood, for others in the early middle ages, and for the rest of the patients, in their 50s, 60s or 70s. Such variation can be understood in terms of the heterogeneity of the disease: for some patients, the disease may be heritable or resulting from critical mutations occurring during gestation or early childhood while for others, the initiating genomic events could have occurred later or involved a long promotion time between the first genomic event and the emergence of the first malignant cell. The possibility of very early onset of adult cancers is well recognized. As suggested in [61,62], the earlier the occurrence of critical mutations leading to a particular cancer the larger is the number of stem cells carrying these mutations and hence also the risk of early cancer onset. Further evidence for the possibility of early onset of prostate cancer comes from the good agreement between the estimates of the age T at cancer onset provided by the Full and Homogeneous models for the majority of patients, see Tables 3 and 4 and note that T is an independent biological parameter of the Homogeneous model while within the Full model it is computed through Equation (25) on the basis of other model parameters.

Heterogeneity of Tumor Growth Rates
The pre-treatment rate of growth of the primary tumor displayed substantial variability among the 12 cancer patients analyzed. As shown in Table 3, these rates varied between 0.4 and 34.1 year −1 (equivalently, the tumor doubling times varied between 7.4 and 632.5 days). The rates of growth of metastases before and after the start of treatment also varied widely between the patients. Finally, for 9 out of 12 patients, the rate of growth of metastases after the start of treatment exceeded the pretreatment growth rate of the primary.

Effects of the Systemic Treatment on the Primary Tumor
Among 10 patients who received systemic treatment, 7 patients responded to the treatment with a very fast reduction in the size of the primary tumor. The response of patient 6 was much slower: after the start of treatment the volume shrinkage half-life of his primary tumor was about 149 days. Finally, for patients 10 and 12, the half-life of their primary tumor reduction was 6.9 and 7.7 years, respectively. Thus, according to the model, systemic treatment of their primary tumors essentially failed.

Metastasis Shedding
According to Equation (1) the rate of metastasis shedding depends critically on the value of parameter θ. The latter was found to be uniformly small for all patients analyzed. Thus, the rate of metastasis shedding is essentially constant in time as long as primary tumor remains in situ. This also implies that treatment of the primary tumor has only weak effect on the rate of metastasis shedding and consequently on the total number of viable secondary tumors. This is illustrated in Figure 4. Figure 4. The rate of metastasis shedding Φ θ for patient 1 (assuming α =10 −6 , see Equation (1)) as a function of primary tumor volume v corresponding to the tumor size Φ (v = Φ c, where c = 10 −9 cm 3 is taken to represent the average volume of a single tumor cell). The curve labeled "expected (θ = 2/3)" refers to the plausible (yet incorrect) assumption that metastases are generated by actively proliferating cancer cells localized at the tumor boundary. The curve labeled "observed (θ = 4.2×10 −5 )" describes the shedding rate for the value of θ estimated from the data. This, essentially constant, shedding rate implies that, contrary to the traditional belief, treatment of the primary is unlikely to substantially reduce the probability of metastatic failure. The log-plots of Φ θ as functions of time rather than volume for the same values of θ are represented by two lines: one with the slope 2β 0 /3 and the other essentially horizontal. As discussed in the text, this supports the notion of prostate cancer stem cells.

Onset of Metastasis
Our analysis confirms unequivocally that in all patients metastatic dissemination of prostate cancer occurred soon after the onset of the disease and much earlier than the appearance of a clinically detectable primary tumor. In fact, according to the Full model, the time between the onset of the disease and inception of the first metastasis never exceeded 2.3 years, see Table 3 and Figure 5. At that time, the primary tumor was extremely small and certainly undetectable. Shedding of the first viable metastasis occurred even earlier. Thus, for all patients, their disease was systemic at the outset, in agreement with hypothesis A1.

Metastasis Latency
The mean metastasis latency times ρ computed through the Full model ranged from a few days to as long as 16.3 years, see Table 3 and Figure 6. This supports the notion of metastasis dormancy (hypothesis A2). For patients with long latency times, many metastases were still occult at the time of surveying. It can be hypothesized that the significant inter-patient variability of the latency times is due to heterogeneity in the genetic make-up of cancer cells and conditions of host microenvironment. Figure 6. Expected metastasis latency times, ρ, in patients 1-12. ρ represents the average time spent by a viable metastasis between detachment from the primary tumor and onset of irreversible proliferation in a given secondary site (here, bone).

Inception of Metastases
Importantly, parameter A, see Equation (12), represents the size at age W of a metastasis whose inception occurred at the start of treatment (see Figure 1). Comparison of the values of this parameter with the observed volumes of metastases shows that in all patients inception of all or most of the detected metastases occurred prior to the start of treatment. Additionally, these early metastases have the largest volumes. Therefore, resection or irradiation of the primary tumor (or any other form of local treatment) has only minor effect on the number of metastases relevant for patients' survival.

Treatment-induced Acceleration of Metastatic Growth
The effect of treatment of the primary on metastatic growth is characterized by the ratio γ 1 /γ 0 of the rate of growth of bone metastases after the start of treatment to their pre-treatment growth rate. For all patients, this ratio was larger than 1, see Figure 7. According to the Full model, the smallest one, 3.5, occurred for patient 1 whose treatment consisted of surgery, external beam radiation and chemotherapy while the largest, 504, occurred for patient 11 who received chemo-and hormonal therapies alone. For patient 2, the second patient who had surgery and chemotherapy, the value of γ 1 /γ 0 was 35.3. For all other patients, the ratio varied between 4 and 126. Thus, all modes of treatment of the primary lead to a dramatic irreversible exacerbation of the disease. At the same time, this suggests that primary tumor in situ has a strong inhibiting effect on the growth of secondary tumors. In fact, for all patients analyzed, the pre-treatment rate of growth of metastases was smaller than the rate of growth of the primary by an order of magnitude, see Table 3. Figure 7. The ratios γ 1 /γ 0 of the rates of growth of bone metastases after the start of treatment and before treatment for patients 1-12 are all significantly larger than 1. In case of chemotherapy, this unexpected feature may be the result of a treatment-related selection of most resistant and fastest growing metastases while in the case of surgery and possibly radiation it is likely due to the treatment-induced acceleration of the growth and vascularization of dormant or slowly growing secondary tumors, see Sections 1 and 9.
As discussed in Section 3, the mechanism is likely to involve treatment-related weakening or total abrogation of the metastasis-inhibiting action of the primary tumor, as well as direct accelerating impact of treatment on metastatic growth. For surgery and radiation, the latter is caused by local and systemic production of growth and angiogenesis factors due to wound healing processes while in the case of chemo-and hormonal therapy it may be due to selection of fastest growing and most resistant cells. Finally, it is well-known [28] that avascular tumors cannot reach sizes exceeding a few mm in diameter. This suggests that at the start of treatment bone metastases in all patients were avascular, while at the time of their detection they reached considerable sizes unattainable for avascular tumors. Therefore, it is very likely that one of the most notable effects of treatment is turning on the angiogenic switch. The question as to how chemo-and/or hormonal therapies bring about this effect calls for further study. Thus, our analysis unequivocally supports hypotheses B1 and B2. Finally, accelerated growth of metastases after the start of treatment brings about a sharp increase in the total metastatic burden; see Figure 8, where the total volume of all detectable metastases is compared to the volume of the primary tumor. Figure 8. The volume of the primary tumor and the total volume of all detectable metastases represented as functions of age for patient 1. It is notable that as time progresses, the total metastatic volume by far exceeds the volume of the primary at the time of surgery. The total volume of metastases is dominated by the volume of the largest metastasis (that is, the metastasis with the earliest inception time, see also Figure 5).

Discussion and Conclusions
The field of oncology has so far been dominated by the tacit notion that treating the primary tumor reduces the chance of distant metastatic failure and is thereby beneficial with respect to the patient's life expectancy. This belief is based on four premises. First, the probability of developing metastases increases in direct relation to the volume of the primary tumor. Second, the longer the primary tumor is in situ, the larger the number of metastases it produces. Third, the number of metastases at presentation is insignificant, i.e., the onset of metastasis is typically a late event in relation to the development of a clinically observable primary tumor. Fourth, it is taken as axiomatic that treatment of the primaryeven if suboptimal in terms of controlling the distant spread of the disease -does not cause a detrimental outcome. If correct, the sine qua non of the above notion is that the earlier and more aggressive the treatment the better the outcome.
Our results directly challenge the validity of these premises. It follows from Equations (9) and (10) that the probability of developing metastases and their expected number increase with the time the primary tumor is in situ. However, since θ is invariably very small, this dependence is not critical. For example, in our analysis the expected number of viable metastases grows essentially linearly in time. Equally important, the third premise appears to be in error: by the time of primary tumor diagnosis there are already a large number of viable metastases dormant or slowly growing at various secondary sites, which makes the reduction of the number of subsequently formed secondary tumors resulting from the treatment of the primary tumor practically irrelevant. Finally, and most importantly, based on our (admittedly limited) analysis, surgery and chemotherapy of metastatic prostate cancer accelerate the progression of the disease.
While in the case of surgery the mechanisms underlying this effect have been known for decades and are fairly well understood, this is not the case for systemic treatment. A number of important questions arise and call for further study: (1) what is the role of selection of resistant and/or fast proliferating cancer cells in the sharp increase in the average rate of growth of metastases predicted by the model? (2) Does systemic treatment cause a drop in the number of metastases or only elimination of slowly growing cells within each metastasis followed by its repopulation by fast growing resistant cells? (3) Does elimination of metastases or reduction in their size have the same boosting effect on other metastases as would the elimination of the primary tumor?
If confirmed, the results of this work potentially could have profound effects on the strategies of cancer treatment. They would place more emphasis on the control of cancer progression through maintaining the state of homeostasis including dormancy, in particular by encouraging watchful waiting. Although such conservative strategies carry significant risks of their own, they should be balanced against the impending risks of aggressive treatment demonstrated in this work.