On the growth and dissemination laws in a mathematical model of metastatic growth

Metastasis represents one of the main clinical challenge in cancer treatment since it is associated with the majority of deaths. Recent technological advances allow quantification of the dynamics of the process by means of noninvasive techniques such as longitudinal tracking of bioluminescent cells. The metastatic process was simplified here into two essential components -- dissemination and colonization -- which were mathematically formalized in terms of simple quantitative laws. The resulting mathematical model was confronted to in vivo experimental data of spontaneous metastasis after primary tumor resection. We discuss how much information can be inferred from confrontation of theories to the data with emphasis on identifiability issues. It is shown that two mutually exclusive assumptions for the secondary growth law (namely same or different from the primary tumor growth law) could fit equally well the data. Similarly, the fractal dimension coefficient in the dissemination law could not be uniquely determined from data on total metastatic burden only. Together, these results delimitate the range of information that can be recovered from fitting data of metastatic growth to already simplified mathematical models.


Introduction
Metastasis (from the greek µετά = change and στάσιζ = place) is a process by which secondary tumors emerge at distant sites from the location of the primary disease [Weiss, 2000]. This is the result of cells detaching the primary tumor, invading the surrounding tissue, traveling to distant sites (mostly through hematogenous ways) before establishing secondary colonies in distant organs, such as the lungs or the liver (the two organs most often affected nowadays) [Talmadge andFidler, 2010, Weiss, 2000]. Despite very complex processes happening at various biological scales, the phenomenon can be summarized into two major phases: dissemination and colonization [Chaffer andWeinberg, 2011, Steeg, 2006].
In recent years, novel experimental techniques shed new lights on the process of metastasis [Sahai, 2007]. On the one hand, the advent of molecular biology revealed the existence and importance of several processes, such as the epithelial-to-mesenchymal transition and its regulation [Chaffer and Weinberg, 2011, Valastyan and Weinberg, 2011 or the association between genetic signatures and the probability of metastatic relapse [van de Vijver et al., 2002], of crucial importance considering that metastases are the ones that ultimately kill cancer patients. On the other hand, other processes were discovered that happen at the system scale of the organism. These include the distant suppression of metastatic growth by endogenous inhibitors of angiogenesis [O'Reilly et al., 1994] or the self-seeding phenomenon [Norton andMassagué, 2006, Kim et al., 2009]. Quantification of the dynamics of these processes is now conceivable thanks to non-invasive experimental techniques that allow to longitudinally track the evolution of the disease [Francia et al., 2011]. One of these consist in transfecting cancer cells with an excitable protein such as the luciferin. After orthotopic implantation of these cells in mice, when the enzyme luciferase is subsequently injected to the animals, the cells emit photons, and the emission intensity can be externally recorded.
The advent of these new experimental techniques generates new data that call for quantitative models in order to analyze them and infer general laws of the disease development. Maybe due to the relative scarcity of data about the dynamics of metastasis, relatively few mathematical models have been proposed that could be used to test hypotheses and quantify aspects of the process [Michor et al., 2006, Bartoszyński, 1987, Hanin, 2013, Retsky et al., 1997, Yorke et al., 1993, Hartung et al., 2014, Haeno et al., 2012. Of note however, in 2000, Iwata, Kawasaki and Shigesada proposed a model particularly well suited for description of the time dynamics of a population of secondary tumors [Iwata et al., 2000]. This model can be easily adapted to longitudinal data of total metastatic burden and has been recently shown able to fit such data in several animal models [Hartung et al., 2014. Based on one of this study, we will focus here on the identifiability aspects of particular coefficients of the model, respectively related to the growth and dissemination laws. In other words, our aim is to determine what can and what cannot be discriminated based on the data that we dispose on one hand, and the formalism of [Iwata et al., 2000] on the other.

Model
The models that were confronted in this study were largely inspired by the initial model of Iwata, Kawasaki and Shigesada [Iwata et al., 2000]. We refer to [Benzekry et al., 2015] for a detailed description. Briefly, dynamics of the primary tumor volume V p (t) at time t is defined by the following Cauchy problem: where g p is the primary tumor (PT) growth law (here, either Gompertz or exponential) and V i is the number of cells at injection, converted into the appropriate unit (either photons per seconds for bioluminescence data or mm 3 for caliper-measured volumes). Metastatic development is then defined by two components: the dissemination law d(V p ) and the colonization (or growth) law g(v) for a tumor of volume v. These shall be discussed in details below. Specification of these two laws allows to write a partial differential equation for the time evolution of the density of metastases ρ(t, v) structured by the size v of the lesions [Iwata et al., 2000]: where V 0 is the size at which metastasis are born (here assumed to be the size of one cell). For calibration of the conversion ratio from number of cells to bioluminescence, we refer the reader to [Benzekry et al., 2015]. In (2), the first equation expresses conservation of the number of metastases when growing in size, the second equates the entering flux of new tumors with the rate of (successful) dissemination from the primary tumor and the last one is the initial condition. From (2), the main quantity of interest for confrontation to bioluminescence data is the total metastatic burden:

Results
We investigated growth laws being either Gompertz (or Gomp-Exp) or exponential. The dissemination law was assumed to have the following expression (see below for modeling details) d(V p ) = µV γ p and identifiability of γ was considered. When not otherwise specified, γ = 1.

Growth law
Material and methods generating the data employed here are extensively described in [Benzekry et al., 2015]. For the growth of LM2-4 luc+ cells orthotopically xenografted in the mammary fat pad of severe combined immune-deficient mice, as employed here, it had been previously demonstrated that the growth kinetics can be accurately described using the Gompertz model [Benzekry et al., 2014] Arguing that for small volumes (such as the one of metastases at initiation), the specific growth rate (i.e. g(v) v ) is bounded, and choosing as a reasonable bound the in vitro specific growth rate, we rather considered the Gomp-Exp model [Wheldon, 1988]: where λ is the in vitro proliferation rate, retrieved from preliminary experiments (see [Benzekry et al., 2015]). We then investigated what structural and parametrical shape of the secondary growth coefficient g(v) would better fit the data, while ensuring reasonable identifiability of its parameters. Four scenarii were considered: A) same growth between the PT and the metastases, B) Gomp-Exp for the PT and exponential for the metastases, C) Gomp-Exp for both the PT and the metastases, with α constrained to be identical for the two and D) Gomp-Exp for both the PT and the metastases, with both parameters α and β allowed to differ between the PT and the metastases. To use the information at the population level and increase robustness of the parameters estimation (which was found very poor on individual growth curves alone), we used the nonlinear mixed-effects statistical framework to fit the data [Lavielle, 2014]. Briefly, it consists in maximizing the likelihood of the entire dataset over all individuals, under the assumption of a distribution of the parameters within the population (which was assumed here to be lognormal). Plain and dashed lines in Figure 1 are respectively the medians and 10th and 90th percentiles of the outputs of simulations of the model under this distribution.
Selection of the best model was made based on the following considerations: i) visual accuracy of the fit (A-D in Figure 1) and statistical criteria of goodness-of-fit such as the Akaike Information Criterion (F in Figure 1), but also ii) practical identifiability of the parameters, as quantified by the normal standard errors reported in Figure 1.E. As can be observed, model B) has to be rejected for inaccuracy, while model D is clearly over-parameterized (see the values of criteria in Figure 1.F and more importantly the uncertainty on the parameters estimation, especially the normalized standard error on µ in Figure 1.E). With similar visual accuracy of the population fits (also observable in individual fits of particular mice, results not shown here), the model C) nevertheless generated substantially higher uncertainty on the parameter estimation, especially parameter µ, with respective normalized standard error of 26.5% for the "same growth" model and 72.1% for the "different growth with alpha fixed" model. Additionally, probably due to sharper estimation of the parameters, predictive performances were improved with the "same growth" model (results not shown). On the other hand, results of the AIC suggest that the addition of one degree of freedom in model C) yields significant improvement of the goodness-of-fit.
However, under the model C), the difference inferred for parameter β translated into minor differences between primary and secondary growth. Together, these results suggest that, although definitive determination of the relationship between the primary and secondary tumor growth law cannot be assessed based only on data on PT and MB growth, the theory assuming equality of the growth rates is conceivable and the most parsimonious.

Dissemination law
With the growth law fixed to the "same growth" model, we investigated determination of the shape of the dissemination law d(V p ). It has been proposed to have the following shape [Iwata et al., 2000]: with µ a parameter quantifying the per cell and per day probability of generating a new metastasis, and γ possibly related to the fractal dimension of the primary tumor vasculature. Parameter γ thus controls the number of cells that are susceptible to leave the tumor while µ can be linked to a more intrinsic metastatic potential of the disease, aggregating chance of success to multiple steps of the metastatic cascade (including genetic mutations, epithelial-to-mesenchymal transition, blood vessel intravasation, survival in transit, extravasation at the distant site and colonization of the new organ's parenchyma).
We wondered whether the value of γ could be determined given the data we dispose. As appears in the results reported in Figure 2, several values generated relatively similar fits. Shown are only three representative values between 0 and 1 but every value of γ investigated gave similar curves. Although the value of γ = 0.5 generated a slightly better fit than γ = 0 or γ = 1, we did not consider this strong enough to support rejection of either of the possibilities for γ. We   Figure 1: Population fits of the breast xenograft data under different growth theories. A: Same growth = same Gompertz growth parameters (α and β) for primary and secondary tumors.B: Different growth Exp = exponential growth law for the metastases. C: Different growth alpha fixed = for each animal, same value of parameter α was imposed while value of β was allowed to vary between the PT and the secondary tumors. D: Different growth = the two Gompertz parameters α and β were allowed to vary between the PT and the metastases.

Discussion
Any data-based modeling inference is limited by two major aspects: 1) the richness of the data and 2) the complexity of the model. Only when taking into account the balance between these two aspects can we reject hypotheses and thus improve our knowledge about natural mechanisms.
Metastasis is a complex process that remains poorly understood, especially during the colonization phase (after extravasation at the distant site). Using a mathematical formalism of basic laws of the process (dissemination and growth), we evaluated here the limits of what can be learned from longitudinal noninvasive data of post-surgery spontaneous metastasis. Our results demonstrated that similar growth between the primary tumor and the metastases was a conceivable theory for metastatic development in a xenograft breast cancer animal model. Exponential growth of the metastases had to be rejected while, on the other hand, a model with different growth of the metastases could also explain the data. Further quantification of metastatic growth are required to discriminate between the "same growth" and "different growth alpha fixed" models, although growth of secondary tumors was not substantially different when inferred under the second model. To this end, longitudinal follow up of individual metastatic lesions in mice using magnetic resonance imaging might be of great help [Baratchart et al., 2015]. Non-uniqueness was also observed for the values of the fractal dimension of the vasculature in the dissemination coefficient able to fit the data. This justifies the use of γ = 1, i.e. the simplest theory (all cells within the primary tumor have equal probability of generating a metastasis) to describe the data [Benzekry et al., 2015]. Additional data are required to discriminate if a particular value of γ is better adapted for specification of the dissemination law.
Together, our results demonstrate that, even with a simple mathematical framework (3 de-grees of freedom in the most adapted version) and relatively rich data, several mutually exclusive models can fit data of total metastatic burden dynamics, leaving open the question of definitive determination of the growth and dissemination laws of metastatic development.