Introduction Neoplastic growth involves a large number of complex biological processes, including regulation of proliferation and control of the cell cycle, stromal recruitment, angiogenesis and escape from immune surveillance. In combination, these cooperate to produce a macroscopic expansion of the tumor volume, raising the prospect of a possible general law for the global dynamics of neoplasia. Quantitative and qualitative aspects of the temporal development of tumor growth can be studied in a variety of experimental settings, including in vitro proliferation assays, three-dimensional in vitro spheroids, in vivo syngeneic or xenograft implants (injected ectopically or orthotopically), transgenic mouse models or longitudinal studies of clinical images. Each scale has its own advantages and drawbacks, with increasing relevance tending to coincide with decreasing measurement precision. The data used in the current study are from two different in vivo systems. The first is a syngeneic Lewis lung carcinoma (LLC) mouse model, exploiting a well-established tumor model adopted by the National Cancer Institute in 1972 [1]. The second is an orthotopic human breast cancer xenografted in severe combined immunodeficient (SCID) mice [2]. Tumor growth kinetics has been an object of biological study for more than 60 years (see e.g. [3] as one of the premiere studies) and has been experimentally investigated extensively (see [4] for a thorough review and [5]–[8] for more recent work). One of the most common findings for animal [9] and human [10]–[12] tumors alike is that their relative growth rates decrease with time [13]; or equivalently, that their doubling times increase. These observations suggest that principles of tumor growth might result from general growth laws, often amenable to expression as ordinary differential equations [14]. The utility of these models can be twofold: 1) testing growth hypotheses or theories by assessing their descriptive power against experimental data and 2) estimating the prior or future course of tumor progression [9], [15] either as a personalized prognostic tool in a clinical context [16]–[20], or in order to determine the efficacy of a therapy in preclinical drug development [21], [22]. Cancer modeling offers a wide range of mathematical formalisms that can be classified according to their scale, approach (bottom-up versus top-down) or integration of spatial structure. At the cellular scale, agent-based models [23], [24] are well-suited for studies of interacting cells and implications on population-scale development, but computational capabilities often limit such studies to small maximal volumes (on the order of the mm3). The tissue scale is better described by continuous partial differential equations like reaction-diffusion models [19], [25] or continuum-mechanics based models [26], [27], when spatial characteristics of the tumor are of interest. When focusing on scalar data of longitudinal tumor volume (which is the case here), models based on ordinary differential equations are more adapted. A plethora of such models exist, starting from proliferation of a constant fraction of the tumor volume, an assumption that leads to exponential growth. This model is challenged by the aforementioned observations of non-constant tumor doubling time. Consequently, investigators considered more elaborate models; the most widely accepted of which is the Gompertz model. It has been used in numerous studies involving animal [9], [28]–[31] or human [12], [15], [30], [32] data. Other models include logistic [30], [33] or generalized logistic [11], [31] formalisms. Inspired by quantitative theories of metabolism and its impact on biological growth, von Bertalanffy [34] derived a growth model based on balance equations of metabolic processes. These considerations were recently developed into a general law of biological growth [35] and brought to the field of tumor growth [36], [37]. When the loss term is neglected, the von Bertalanffy model reduces to a power law (see [5], [38] for applications to tumor growth). An alternative, purely phenomenological approach led others [39] to simply consider tumor growth as divided into two phases: an initial exponential phase then followed by a linear regimen. Recently, influences of the microenvironment have been incorporated into the modeling, an example being the inclusion of tumor neo-angiogenesis by way of a dynamic carrying capacity [40], [41]. Although several studies have been conducted using specific mathematical models for describing tumor growth kinetics, comprehensive work comparing broad ranges of mathematical models for their descriptive power against in vivo experimental data is lacking (with the notable exception of [30] and a few studies for in vitro tumor spheroids [33], [42]–[44]). Moreover, predictive power is very rarely considered (see [42] for an exception, examining growth of tumor spheroids), despite its clear relevance to clinical utility. The aim of the present study is to provide a rational, quantitative and extensive study of the descriptive and predictive power of a broad class of mathematical models, based on an adapted quantification of the measurement error (uncertainty) in our data. As observed by others [45], specific data sets should be used rather than average curves, and this is the approach we adopted here. In the following sections, we first describe the experimental procedures that generated the data and define the mathematical models. Then we introduce our methodology to fit the models to the data and assess their descriptive and predictive powers. We conclude by presenting the results of our analysis, consisting of: 1) analysis of the measurement error and derivation of an appropriate error model, subsequently used in the parameters estimation procedure, 2) comparison of the descriptive power of the mathematical models against our two datasets, and 3) determination of the predictive abilities of the most descriptive models, with or without adjunction of a priori information in the estimation procedure. Materials and Methods Ethics statement Animal tumor model studies were performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. Protocols used were approved by the Institutional Animal Care and Use Committee (IACUC) at Tufts University School of Medicine for studies using murine Lewis lung carcinoma (LLC) cells (Protocol: #P11-324) and at Roswell Park Cancer Institute (RPCI) for studies using human LM2-4LUC+ breast carcinoma cells (Protocol: 1227M). Institutions are AAALAC accredited and every effort was made to minimize animal distress. Mice experiments Cell culture Murine Lewis lung carcinoma (LLC) cells, originally derived from a spontaneous tumor in a C57BL/6 mouse [46], were obtained from American Type Culture Collection (Manassas, VA). Human LM2-4LUC+ breast carcinoma cells are a metastatic variant originally derived from MDA-MD-231 cells and then transfected with firefly luciferase [47]. All cells were cultured in high glucose DMEM (obtained from Gibco Invitrogen Cell Culture, Carlsbad, CA or Mediatech, Manassas, VA) with 10%FBS (Gibco Invitrogen Cell Culture) and 5% CO2. Tumor injections For the subcutaneous mouse syngeneic lung tumor model, C57BL/6 male mice with an average lifespan of 878 days were used [48]. At time of injection mice were 6 to 8 weeks old (Jackson Laboratory, Bar Harbor, Maine). Subcutaneous injections of 106 LLC cells in 0.2 ml phosphate-buffered saline (PBS) were performed on the caudal half of the back in anesthetized mice. For the orthotopic human xenograft breast tumor model, LM2-4LUC+ cells (1×106 cells) were orthotopically implanted into the right inguinal mammary fat pads of 6- to 8-week-old female severe combined immunodeficient (SCID) mice obtained from the Laboratory Animal Resource at RPCI, as previously described [2]. Tumor measurements Tumor size was measured regularly with calipers to a maximum of 1.5 cm3 for the lung data set and 2 cm3 for the breast data set. Largest (L) and smallest (w) diameters were measured subcutaneously using calipers and the formula was then used to compute the volume (ellipsoid). Volumes ranged 14–1492 mm3 over time spans from 4 to 22 days for the lung tumor model (two experiments of 10 animals each) and 202–1902 mm3 over time spans from 18 to 38 days for the breast tumor data (five experiments conducted with a total of 34 animals). Plots of individual growth curves for both data sets are reported in Figure S1. Mathematical models For all the models, the descriptive variable is the total tumor volume, denoted by V, as a function of time t. It is assumed to be proportional to the total number of cells in the tumor. To reduce the number of degrees of freedom, all the models (except the exponential V 0) had a fixed initial volume condition. Although the number of cells that actually remain in the established tumor is probably lower than the number of injected cells (∼60–80%), we considered 1 mm3 ( cells [49], i.e. the number of injected cells) as a reasonable approximation for V(t = 0). Exponential-linear models The simplest theory of tumor growth presumes all cells proliferate with constant cell cycle duration TC . This leads to exponential growth, which is also valid in the extended cases where either a constant fraction of the volume is proliferating or the cell cycle length is a random variable with exponential distribution (assuming that the individual cell cycle length distributions are independent and identically distributed). As one modification, initial exponential phase can be assumed to be followed by a linear growth phase [39], giving the following Cauchy problem for the volume rate of change (growth rate): (1) Here, the coefficient a 0 is the fraction of proliferative cells times ln 2/TC where TC is either the constant cell cycle length or the mean cell cycle length (under the assumption of exponentially distributed cell cycle lengths). The coefficient a 1 drives the linear phase. Assuming that the solution of the problem (1) is continuously differentiable uniquely determines the value of τ as . The coefficient V 0 denotes the initial volume. From this formula, three models were considered: a) initial volume fixed to 1 mm3 and no linear phase (a 1 = +∞), referred to hereafter as exponential 1, b) free initial volume and no linear phase, referred to as exponential V 0 and c) equation (1) with fixed initial volume of 1 mm3, referred to as the exponential-linear model. Logistic and Gompertz models A general class of models used for quantification of tumor growth kinetics have a sigmoid shape, i.e. an increasing curve with one inflection point that asymptotically converges to a maximal volume, the carrying capacity, denoted here by K. This qualitatively reproduces the experimentally observed growth slowdown [9]–[12] and is consistent with general patterns of organ and organismal growth. The logistic model is defined by a linear decrease of the relative growth rate in proportion to the volume: (2) where a is a coefficient related to proliferation kinetics. This model can be interpreted as mutual competition between the cells (for nutrients or space, for instance), by noticing that under this model the instantaneous probability for a cell to proliferate is proportional to . The logistic model has been used for description of tumor growth, for instance, in [30]. Others (such as [11]) have considered a generalization of the logistic equation, defined by (3) that will be referred to as the generalized logistic model. Equation (3) has the explicit solution which also provides an analytic solution to model (2) when ν = 1. When a different parameterization is employed, this model converges when ν→0 to the Gompertz model, defined by (4) Coefficient a is the initial proliferation rate (at V = 1 mm3) and β is the rate of exponential decay of this proliferation rate. Although first introduced in [50] for a different purpose – the description of human mortality for actuarial applications – the Gompertz model became a widely-accepted representation of growth processes in general [51] and of tumor growth in particular. It was first successfully used in this regard [28] before its applicability was confirmed on large animal data sets [9], [29] and for human breast data [32]. The essential characteristic of the Gompertz model is that it exhibits exponential decay of the relative growth rate. An analytic formula can be derived for the solution of (4): where we can see that asymptotically, the volume converges to a carrying capacity given by . A unified model deriving these three sigmoidal models from specific biophysical assumptions about different types of cellular interactions can be found in [52]. Dynamic carrying capacity Taking the next step up in complexity brings us to a model that assumes a dynamic (time-dependent) carrying capacity (CC) [40], [41] that can be taken, for example, to represent the tumor vasculature. If one assumes that stimulation of the carrying capacity is proportional to the tumor surface, and neglects angiogenesis inhibition, this model can be formulated in terms of two coupled equations: (5) and will be referred to as the dynamic CC model. It should be noted that this model was first developed with the intent of modeling the effect of anti-angiogenic therapies on tumor growth and not strictly for describing or predicting the behavior of V alone. However, we integrated it into our analysis in order to investigate and quantify whether consideration of a dynamic carrying capacity could benefit these tasks. Von Bertalanffy and power law Von Bertalanffy [34], followed later on by others [35], proposed to derive general laws of organic growth from basic energetics principles. Stating that the net growth rate should result from the balance of synthesis and destruction, observing that metabolic rates very often follow the law of allometry (i.e. that they scale with a power of the total size) [34] and assuming that catabolic rates are in proportion to the total volume, he derived the following model for growth of biological processes (6) Employing our usual assumption that V(t = 0) = 1 mm3, we will refer to this model as the von Bertalanffy model (note that others [14], [30] often identify this model as the specific case γ = 2/3, termed “second type growth” in [34]). It has already been successfully applied to describe tumor growth [36], [37]. More elaborate considerations linking tumor growth, metabolic rate and vascularization leading to equation (6) can be found in [37]. That work also provides expressions of the coefficients in terms of measurable energetic quantities. Explicit solution of the model is given by From the observation that our data does not exhibit a clear saturation phase, a qualitative feature of equation (6), we also considered another model, derived from (6), by neglecting the loss term, i.e. taking b = 0. This model will be termed the power law model. Pushing further the reasoning of [34] and arguing that the rate of synthesis of new material, in the context of tumor growth, should be proportional to the number of proliferative cells (under the assumption of a constant cell cycle length), this model suggests that the proliferative tissue is proportional to Vγ . This could be further interpreted as a possible fractional Hausdorff dimension of the proliferative tissue, when viewed as a metric subspace of the full tumor volume (viewed itself as a three-dimensional subset of the three-dimensional Euclidean space). This dimension would be equal to 3γ and could be less than 3 when γ 0.2 by Student's t-test). Among these simulation replicates, we only considered as significant the cases where , for the same reasons as explained earlier. For the lung tumor data set, this did not lead to any exclusion for most of the situations, the only exceptions being for S 3,5 and S 3,6 where only 89/100 and 72/100 replicates were eligible, respectively. In contrast, for the breast tumor data and depths 1 to 10, respectively 99, 16, 76, 3, 100, 3, 100, 0, 34 and 77 replicates were eligible. Therefore, results of S 3,2, S 3,4, S 3,6, S 3,8 and S 3,9 were considered non-significant and were not reported. Results Measurement error The following method was used for analysis of the error made when measuring tumor volume with calipers. One volume per time point per cage was measured twice within a few minutes interval. This gave a total of 133 measurements over a wide range of volumes (20.7–1429 mm3). These were subsequently analyzed by considering the following statistical representation where Y is a random variable whose realizations are the measured volumes, is the true volume, ε is a reduced centered Gaussian random variable, and σE is the error standard deviation. The two measures, termed y 1 and y 2, were, as expected, strongly correlated (Figure 1A, ). Statistical analysis rejected variance independent of volume, i.e. constant E ( , χ 2 test) and a proportional error model (E = Y) was found only weakly significant ( , χ 2 test, see Figure 1B). We therefore introduced a dedicated error model, defined by (22) Two main rationales guided this formulation. First, we argued that error should be larger when volume is larger, a fact that is corroborated by larger error bars for larger volumes on growth data reported in the literature (see Figure 4 in [2] for an example among many others). This was also supported by several publications using a proportional error model when fitting growth data (such as [42], [60]). Since here such a description of the error was only weakly significant, we added a power to account for lower-than-proportional uncertainty in large measurements. Second, based on our own practical experience of measuring tumor volumes with calipers, for very small tumors, the measurement error should stop being a decreasing function of the volume because of detectability limits. This motivated the introduction of the threshold Vm . After exploration of several values of Vm and α, we found to be able to accurately describe dispersion of the error in our data ( , χ 2 test, see Figure 1C). This yielded an empirical value of 10.1371/journal.pcbi.1003800.g001 Figure 1 Volume measurement error. A. First measured volume y 1 against second one y 2. Also plotted is the regression line (correlation coefficient R = 0.98, slope of the regression = 0.96). B. Error against approximation of the volume given by the average of the two measurement . The χ 2 test rejected Gaussian distribution of constant variance ( ) C. Histogram of the normalized error applying the error model given by with α = 0.84 and Vm = 83 mm3. It shows Gaussian distribution (χ 2 test not rejected, p = 0.196) with standard deviation . We did not dispose of double measurements for the breast tumor data and the error analysis was performed using the lung tumor data set only. However, the same error model was applied to the breast tumor data, as both relied upon the same measurement technique. This result allowed quantification of the measurement error inherent to our data and was an important step in the assessment of each model's descriptive power. Descriptive power We tested all the models for their descriptive power and quantified their respective goodness of fit, according to various criteria. Two distinct estimation procedures were employed. The first fitted each animal's growth curve individually (minimization of weighted least squares, with weights defined from the error model of the previous section, see Material and Methods). The second method used a population approach and fitted all the growth curves together. Results are reported in Figure 2 and Tables 1 and 2. Parameter values resulting from the fits are reported in Tables 3 and 4. 10.1371/journal.pcbi.1003800.g002 Figure 2 Descriptive power of the models for lung and breast tumor data. A. Representative examples of all growth models fitting the same growth curve (animal 10 for lung, animal 14 for breast). Error bars correspond to the standard deviation of the a priori estimate of measurement error. In the lung setting, curves of the Gompertz, power law, dynamic CC and von Bertalanffy models are visually indistinguishable. B. Corresponding relative growth rate curves. Curves for von Bertalanffy and power law are identical in the lung setting. C. Residuals distributions, in ascending order of mean RMSE (13) over all animals. Residuals (see formula (15) for their definition) include fits over all the animals and all the time points. Exp1 = exponential 1, Exp-L = exponential-linear, Exp V 0 = exponential V 0, Log = logistic, GLog = generalized logistic, PL = power law, Gomp = Gompertz, VonBert = von Bertalanffy, DynCC = dynamic CC. 10.1371/journal.pcbi.1003800.g003 Figure 3 Examples of predictive power. Representative examples of the forecast performances of the models for the lung data set (mouse number 2). Five data points were used to estimate the animal parameters and predict future growth. Prediction success of the models are reported for the next day data point (OK1) or global future curve (OKglob), based on the criterion of a normalized error smaller than 3 (meaning that the model prediction is within 3 standard deviations of the measurement error) for OK1 and the median of this metric over the future curve for OKglob (see Materials and Methods for details). 10.1371/journal.pcbi.1003800.g004 Figure 4 Prediction depth and number of data points. Predictive power of some representative models depending on the number of data points used for estimation of the parameters (n) and the prediction depth in the future (d). Top: at position (n,d) the color represents the percentage of successfully predicted animals when using n data points and forecasting the time point , as quantified by the score (multiplied by 100), defined in (17). This proportion only includes animals having measurements at these two time points, thus values at different row d on the same column n or reverse might represent predictions in different animals. White squares correspond to situations where this number was too low ( 0.05 # Generalized logistic 0.12 (0.019–0.42) [1] −13 (−30–0.98) [3] 0.7 [6] 2114 [6] 0.4 (0.17–0.82) [1] 0.98 (0.94–1) 100 3 Gompertz 0.155 (0.019–0.67) [4] −13.4 (−32–2.4) [1] −7.62 [1] 2108 [5] 0.41 (0.16–0.93) [2] 0.97 (0.82–1) 100 2 Power law 0.155 (0.016–0.71) [5] −13.4 (−34–2.9) [2] −7.59 [2] 2091 [2] 0.41 (0.15–0.96) [3] 0.96 (0.78–1) 100 2 Dynamic CC 0.136 (0.013–0.61) [2] −12.5 (−32–3.5) [4] 1.19 [7] 2063 [1] 0.42 (0.14–0.96) [4] 0.97 (0.82–1) 100 3 Von Bertalanffy 0.14 (0.016–0.67) [3] −12.5 (−32–4.4) [5] 1.19 [8] 2096 [3] 0.42 (0.16–1) [5] 0.97 (0.81–1) 100 3 Exponential V0 0.217 (0.0069–0.91) [6] −10.7 (−34–5.1) [6] −4.85 [3] 2099 [4] 0.49 (0.096–1.1) [6] 0.93 (0.68–1) 100 2 Exponential-linear 0.22 (0.048–0.76) [7] −8.51 (−17–3.8) [7] −2.7 [4] 2174 [7] 0.51 (0.27–1) [7] 0.96 (0.91–0.99) 100 2 Logistic 0.232 (0.05–0.73) [8] −8.34 (−18–3.4) [8] −2.52 [5] 2214 [8] 0.52 (0.27–0.98) [8] 0.96 (0.92–0.99) 100 2 Exponential 1 1.36 (0.31–2.4) [9] 6.01 (−5.4–13) [9] 8.31 [9] 2442 [9] 1.2 (0.59–1.6) [9] 0.64 (0.28–0.94) 15 1 Models were ranked in ascending order of the RMSE, defined by expression (13). For each metric, indicated are the mean value (among all animals) and in parenthesis the minimal and maximal values (not reported for AICc as they were redundant with the range of AIC). When reported, value inside brackets is the rank of the model for the underlying metric. The model ranking first is highlighted in bold. For animal j, is the minimal value of the objective that was minimized in the individual fits approach (see (8)), divided by the number of time points Ij , and represents the variance of the weighted residuals. AIC and AICc are defined in (11) and (12), AICcpop is the AICc resulting from the mixed-effect estimation (see Materials and Methods) and R 2 is defined in (14). Values reported in the p column are percentages of animals for which Kolmogorov-Smirnov test for normality of residuals was not rejected at the significance level of 0.05. # = number of parameters. J = 20 animals. 10.1371/journal.pcbi.1003800.t002 Table 2 Fit performances of growth models for the breast data. Model 1/I χ2 AIC AICc AICcpop RMSE R2 p>0.05 # Exponential-linear 0.0919 (0.016–0.49) [2] −11.7 (−25–1) [1] 0.798 [1] 2832 [1] 0.34 (0.16–0.83) [1] 0.92 (0.66–0.99) 100 2 Gompertz 0.0976 (0.015–0.33) [4] −11.3 (−28–−0.85) [2] 1.21 [2] 2866 [3] 0.35 (0.14–0.68) [2] 0.92 (0.67–0.99) 100 2 Generalized logistic 0.0814 (0.0037–0.33) [1] −10.7 (−26–0.19) [4] 11.9 [6] 2870 [5] 0.36 (0.096–0.76) [3] 0.94 (0.8–0.99) 100 3 Power law 0.102 (0.016–0.32) [5] −10.9 (−22–−0.017) [3] 1.53 [3] 2913 [7] 0.36 (0.16–0.71) [4] 0.92 (0.61–0.99) 100 2 Exponential V0 0.118 (0.011–0.37) [7] −9.86 (−20–1) [5] 2.61 [4] 2870 [4] 0.39 (0.13–0.78) [5] 0.9 (0.56–0.99) 100 2 Von Bertalanffy 0.0928 (0.015–0.32) [3] −9.77 (−26–1.2) [6] 11.9 [7] 2876 [6] 0.39 (0.15–0.8) [6] 0.93 (0.67–0.99) 100 3 Dynamic CC 0.11 (0.018–0.5) [6] −8.58 (−20–3.2) [7] 13 [9] 2862 [2] 0.42 (0.21–0.94) [7] 0.91 (0.58–0.99) 100 3 Logistic 0.145 (0.0037–0.42) [8] −8.1 (−22–−0.13) [8] 4.37 [5] 2921 [8] 0.43 (0.078–0.76) [8] 0.86 (0.65–0.99) 100 2 Exponential 1 2.19 (0.62–3.4) [9] 8.77 (0.62–13) [9] 12.6 [8] 3518 [9] 1.6 (0.85–2) [9] −0.91 (−5.9–0.88) 53 1 Models were ranked in ascending order of the RMSE, defined by expression (13). For each metric, indicated are the mean value (among all animals) and in parenthesis the minimal and maximal values (not reported for AICc as they were redundant with the range of AIC). When reported, value inside brackets is the rank of the model for the underlying metric. The model ranking first is highlighted in bold. For animal j, is the minimal value of the objective that was minimized in the individual fits approach (see (8)), divided by the number of time points Ij , and represents the variance of the weighted residuals. AIC and AICc are defined in (11) and (12), AICcpop is the AICc resulting from the mixed-effect estimation (see Materials and Methods) and R 2 is defined in (14). Values reported in the p column are percentages of animals for which Kolmogorov-Smirnov test for normality of residuals was not rejected at the significance level of 0.05. # = number of parameters. J = 34 animals. 10.1371/journal.pcbi.1003800.t003 Table 3 Parameter values estimated from the fits: Lung data. Model Par. Unit Median value (CV) Mean normalized std error (CV) Power law a [mm3(1-γ)· day−1] 0.921 (38.9) 11.9 (48.7) γ - 0.788 (9.41) 4 (53.4) Gompertz a [day−1] 0.743 (25.3) 6.02 (51.3) β [day−1] 0.0792 (42.4) 13.7 (65.4) Exponential-linear a0 [day−1] 0.49 (19.3) 3.08 (41.5) a1 [mm3· day−1] 115.6 (22.6) 15.7 (40.7) Dynamic CC a [day−1] 0.399 (106) 447 (89.8) b [mm−2· day−1] 2.66 (241) 395 (176) K0 [mm3] 2.6 (322) 6.5e+04 (345) Von Bertalanffy a [mm3(1-γ)· day−1] 7.72 (112) 1.43e+04 (155) γ - 0.947 (13.5) 40.9 (73) b [day−1] 6.75 (118) 2.98e+07 (222) Generalized logistic a [day−1] 2555 (148) 2.36e+05 (137) K [mm3] 4378 (307) 165 (220) ν - 0.00014(199) 2.36e+05 (137) Exponential V0 V0 [mm3] 13.2 (47.9) 28.9 (55) a [day−1] 0.257 (15.4) 7.49 (48.3) Logistic a [day−1] 0.502 (17.5) 3.03 (48.9) K [mm3] 1297 (23.1) 17.2 (43.8) Exponential 1 a [day−1] 0.399 (13.8) 2.87 (24.5) Shown are the median values within the population and in parenthesis the coefficient of variation (CV, expressed in percent and defined as the standard deviation within the population divided by mean and multiplied by 100) that quantifies inter-animal variability. Last column represents the normalized standard errors (nse) of the maximum likelihood estimator, defined in (11). 10.1371/journal.pcbi.1003800.t004 Table 4 Parameter values estimated from the fits: Breast data. Model Par. Unit Median value (CV) Mean normalized std error (CV) Power law a [mm3(1-γ)· day−1] 1.32 (74.1) 31.2 (48.6) γ - 0.58 (23) 12.1 (62.2) Gompertz a [day−1] 0.56 (18.4) 7.52 (43) β [day−1] 0.0719 (26.4) 12.5 (65.5) Exponential-linear a0 [day−1] 0.31 (16.8) 6.22 (65.9) a1 [mm3· day−1] 67.8 (33.2) 12.9 (45.4) Dynamic CC a [day−1] 2.63 (81.3) 597 (339) b [mm−2· day−1] 0.829 (399) 1.33e+03 (571) K0 [mm3] 12.7 (525) 6.48e+03 (361) Von Bertalanffy a [mm3(1-γ)· day−1] 2.32 (113) 1.17e+04 (181) γ - 0.918 (22.5) 128 (65.5) b [day−1] 0.808 (132) 1.48e+08 (300) Generalized logistic a [day−1] 2753 (131) 7.41e+05 (160) K [mm3] 1964 (557) 232 (433) ν - 2.68e-05 (166) 7.41e+05 (160) Exponential V0 V0 [mm3] 68.2 (57.2) 34.5 (50.8) a [day−1] 0.0846 (27.7) 13.7 (44.1) Logistic a [day−1] 0.305 (10.2) 3.17 (34.9) K [mm3] 1221 (31.4) 11.8 (73.8) Exponential 1 a [day−1] 0.223 (5.9) 3.72 (21.3) Shown are the median values within the population and in parenthesis the coefficient of variation (CV, expressed in percent and defined as the standard deviation within the population divided by mean and multiplied by 100) that quantifies inter-animal variability. Last column represents the normalized standard errors (nse) of the maximum likelihood estimator, defined in (11). Figure 2A depicts the representative fit of a given animal's growth curve for each data set using the individual approach. From visual examination, the exponential 1 (1), logistic (2) and exponential-linear (1) models did not well explain lung tumor growth and the exponential 1 (1) and logistic (2) models did not satisfactorily fit the breast tumor growth data. The other models seemed able to describe tumor growth in a reasonably accurate fashion. These results were further confirmed by global quantifications over the total population, such as by residuals analysis (Figure 2C) and global metrics reported in Tables 1 and 2. When considering goodness-of-fit only, i.e. looking at the minimal least squared errors possibly reached by a model to fit the data (metric in Tables 1 and 2), the generalized logistic model (3) exhibited the best results for both data sets (first column in Tables 1 and 2). This indicated a high structural flexibility that allowed this model to adapt to each growth curve and provided accurate fits. On the other hand, the exponential 1 (1) and logistic (2) models clearly exhibited poor fits to the data, a result confirmed by almost all the metrics (with the exception of the AICc). Influence of the goodness-of-fit metric Being able to closely match the data is not the only relevant criterion to quantify the descriptive power of a model since parameter parsimony of the model should also be taken into account. Other metrics were employed that balanced pure goodness-of-fit and the number of parameters (see Materials and Methods for their definitions). Among them, AICc exhibited the strongest penalization for a large number of parameters. However, this metric was in multiple instances in disagreement with the other metrics dealing with parsimony. For this reason, we also reported the values of AIC. These were found globally in accordance with the RMSE. The AICcpop gave a weaker importance to the number of parameters, due to the large number of data points in the setting of the population approach, since all the animals were pooled together. For the same reason, values of AICcpop were almost identical to values of AICpop and only the former were reported. Other structural and numerical differences (for instance, the individual approach used a deterministic optimizer while the population approach was based on a stochastic algorithm) also explained the discrepancies between the two approaches. When comparing the results generated by the two approaches, better individual fits were obtained using the individual approach (see Tables S2). Indeed, the population approach is better designed for settings where the number of data points is too low to individually estimate the parameters, which was not our case. Taking all these considerations into account, we deemed the RMSE metric to be a good compromise and used this criterion for ranking the models in Tables 1 and 2. Descriptive power and identifiability of the models for each data set For the lung data, five models (generalized logistic (3), Gompertz (4), power law (6), dynamic CC (5) and von Bertalanffy (6)) were found to have similar RMSE (Table 1), suggesting an identical descriptive power among them. However, having one less parameter, the Gompertz and power law models had smaller AIC (and much smaller AICc) and should thus be preferred for parsimonious description of subcutaneous tumor growth of LLC cells. Having an additional degree of freedom translated into poor identifiability of the parameters for the generalized logistic (3), dynamic CC (5) and von Bertalanffy (6) models, as indicated by high standard errors on the parameter estimates (last column of Table 3) and low robustness of these estimates with regard to the initialization of the parameters (see the study of practical identifiability of the models in supplementary text S2 and Table S3). The Gompertz model (4) was also supported by the observation that the median value of ν estimated by the generalized logistic model was close to zero. For the breast data, superior fitting power was obtained by the exponential-linear model (1), for all but one of the metrics considered (R 2, see Table 2). For all the animals, the fits were in the linear phase of the model indicating linear tumor growth dynamics in the range of volumes observed. The Gompertz (4), generalized logistic (3) and power law (6) models still had high descriptive power, with mean RMSE and AIC similar to the exponential-linear model (Table 2). Again, as a consequence of their larger number of parameters, the dynamic CC (5), von Bertalanffy (6) and generalized logistic (3) models exhibited very large standard errors of the parameter estimates as well as large inter-animal variability (Table 3). Consequently, moderate confidence should be attributed to the specific values of the parameters estimated by the fits, although this did not affect their descriptive power. As a general result for both data sets, all the models with two parameters were found to be identifiable (Tables 2 and 3). This was confirmed by a study of practical identifiability performed by systematically varying the initial condition of the minimization algorithm (see supporting text S2 and Table S3). For the theories that were able to fit (power law (6) and Gompertz (4) models for the lung tumor data and additionally the exponential-linear model (1) for the breast tumor data), the values of the parameters and their coefficient of variability provided a fairly good characterization of the tumor growth curves dynamics and inter-animal variability. In particular, the power γ of the power law model identified in the lung tumor data set seemed to accurately represent the growth of the LLC experimental model (low standard errors and coefficient of variation). Results of inter-animal variability suggested a larger heterogeneity of growth curves in the breast tumor data than in the lung tumor data set, which could be explained by the different growth locations (orthotopic versus ectopic). Taken together, our results show that, despite the complexity of internal cell populations and tissue organization, at the macroscopic scale tumor growth exhibits relatively simple dynamics that can be captured through mathematical models. Models with three parameters, and more specifically the generalized logistic model (3), were found highly descriptive but not identifiable. For description of subcutaneous in vivo tumor growth of LLC cells, the Gompertz (4) and power law (6) models were found to exhibit the best compromise between number of parameters and descriptive power. Orthotopic growth of LM2-4LUC+ cells showed a clear linear trend in the range of observed volumes, well captured by the exponential-linear (1), power law (6) and Gompertz (4) models. Forecasting tumor growth: Individual curves The two models that were shown unable to describe our data in the previous section, namely the exponential 1 (1) and logistic (3) models, were excluded from further analysis. The remaining ones were assessed for their predictive power. The challenge considered was to predict future growth based on parameter estimation performed on a subset of the data containing only n data points (with n