Despite internal complexity, tumor growth kinetics follow relatively simple laws that can be expressed as mathematical models. To explore this further, quantitative analysis of the most classical of these were performed. The models were assessed against data from two in vivo experimental systems: an ectopic syngeneic tumor (Lewis lung carcinoma) and an orthotopically xenografted human breast carcinoma. The goals were threefold: 1) to determine a statistical model for description of the measurement error, 2) to establish the descriptive power of each model, using several goodness-of-fit metrics and a study of parametric identifiability, and 3) to assess the models' ability to forecast future tumor growth. The models included in the study comprised the exponential, exponential-linear, power law, Gompertz, logistic, generalized logistic, von Bertalanffy and a model with dynamic carrying capacity. For the breast data, the dynamics were best captured by the Gompertz and exponential-linear models. The latter also exhibited the highest predictive power, with excellent prediction scores (≥80%) extending out as far as 12 days in the future. For the lung data, the Gompertz and power law models provided the most parsimonious and parametrically identifiable description. However, not one of the models was able to achieve a substantial prediction rate (≥70%) beyond the next day data point. In this context, adjunction of a priori information on the parameter distribution led to considerable improvement. For instance, forecast success rates went from 14. 9% to 62. 7% when using the power law model to predict the full future tumor growth curves, using just three data points. These results not only have important implications for biological theories of tumor growth and the use of mathematical modeling in preclinical anti-cancer drug investigations, but also may assist in defining how mathematical models could serve as potential prognostic tools in the clinic.
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.
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. 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 V0) 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). For a given animal j and model M, the general setting considered for prediction was to estimate the model' s parameter set using only the first n data points and to use these to predict at a depth d, i. e. to predict the value at time, provided that a measurement exists at this day (in which case it will be denoted by). The resulting best-fit parameter set will be denoted.
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 representationwhere 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 y1 and y2, 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 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. 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. 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). 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<Ij for a given j). We refer to the Materials and Methods section for the definitions of prediction metrics and success scores. When relatively fewer data points were used, for example with only three, individual predictions based on individual fits were shown to be globally limited for the lung tumor data, especially over a large time frame (Figure 4. A, Table 5). However, this situation is likely to be the clinically relevant since few clinical examinations are performed before the beginning of therapy. On the other hand, large databases might be available from previous examinations of other patients and this information could be useful to predict future tumor growth in a particular patient. In a preclinical setting of drug investigation, tumor growth curves of animals from a control group could be available and usable when inferring information on the individual time course of one particular treated animal. An interesting statistical method that could potentiate this a priori information consists in learning the population distribution of the model parameters from a given database and to combine it with the individual parameter estimation from the available restricted data points on a given animal. We investigated this method in order to determine if it could improve the predictive performances of the models. Each dataset was randomly divided into two groups. One was used to learn the parameter distribution (based on the full time curves), while the other was dedicated to predictions (limited number of data points). For a given animal of this last group, no information from his growth curve was used to estimate the a priori distributions. The full procedure was replicated 100 times to ensure statistical significance, resulting in respectively 2000 and 3400 fits performed for each model. We refer to the Materials and Methods for more technical details. Results are reported in Figure 5. Predictions obtained using this technique were significantly improved for the lung tumors, going from an average success score of 14. 9%±8. 35% to 62. 7%±11. 9% (means ± standard deviations) for prediction of the total future curve with the power law model (6) (see Figure 5. A). Prediction success rates were improved even at large future depths. For instance, predictions 7 days in the future reached an average success rate of 50. 6%, power law model (6), see Figure 5. C, while their success rate was very low with direct individual prediction (6. 07%). Prediction successes reached 90% (power law model (6) ) at the closer horizon of the next day data point (), while success rate was only 57. 1% using an individual approach (Figure 5. B). Other small horizon depths also reached excellent prediction scores (Figure 5. C). The largest improvement of success rates for the power law model was observed for that went from an average score of 6. 86% (with standard deviation 7. 47) to an average score of 75. 2% (with standard deviation 12. 9), representing more than an 11-fold increase. We report in Figure S5 the details of predictions with and without a priori information for all the animals within a given forecast group from the lung tumor data set (power law model (6) ). It can be appreciated how additional information on the parameter distribution in the estimation procedure significantly improved global prediction of the tumor growth curves. The impact of the addition of the a priori information was however less important when using more data points for the estimation (results not shown). For the breast data, due to its already high prediction score without adjunction of a priori information, the exponential-linear model did not benefit from the method. For the next day data point of the breast tumor growth curves, predictability was already almost maximal without adjunction of a priori information and thus no important impact was observed. For both data sets, not all the models equally benefited from the addition of a priori information (Figure 5). Models having the lowest parameter inter-animal variability, such as the power law (6), Gompertz (4), exponential-linear (1), and exponential V0 (1) models (Table 3), which also had better practical identifiability (Tables 2 and S3), exhibited great benefit. In contrast, the models with three parameters showed only modest benefit or even decrease of their success rates (see and for the von Bertalanffy model (6) on the breast tumor data in Figure 5. B), with the exception of the generalized logistic model (3) on the breast tumor data. In these cases, adjunction of a priori information translated into poor enhancement of predictive power because the mean population parameters did not properly capture the average behavior within the population and were therefore not very informative. On the other hand, models such as the power law model (6) on the lung tumor data set, whose coefficient γ characterized particularly well the growth pattern (Table 3), had a more informative a priori distribution that translated into the highest improvement of predictive power. For the generalized logistic model (3) on the breast data, the mean parameters were able to inform the linear regimen of the growth phase and thus protected the model from too early saturation. These results demonstrated that addition of a priori information in the fit procedure considerably improved the forecast performances of the models, in particular when using a small number of data points and low-parameterized models for data with low predictability, such as the power law model for the lung tumor data set.
In our analysis, constant variance of the error was clearly rejected and although a proportional error (used by others [33]) was not strictly rejected by statistical analysis (p = 0. 08), a more adequate error model to our data was developed. However, using a proportional or even constant error model did not significantly affect conclusions as to the descriptive power of the models, identifying the same models (Tables 1,2) as most adequate for description of tumor growth (results not shown). Nevertheless, the use of an appropriate error model could have important implications in the quantitative assessment of a model' s descriptive performance and rejection of inaccurate tumor growth theories. For instance, using the same human tumor growth data, Bajzer et al. [60] found the assumption of proportional error variance to favor the Gompertz model for descriptive ability, whereas Vaidya and Alexandro [30] had observed the logistic model to be favored, under a constant-variance assumption. The error model used might additionally have important implications on predictions. Although detailed analysis of the impact of the error model on prediction power is beyond the scope of the present study, we performed a prospective study of predictive properties when using a constant error model on the lung tumor data and found changes in the ranking of the models (results not shown). As expected, our results confirmed previous observations [9]–[11], [13], [29], [32] that tumor growth is not continuously exponential (constant doubling time) in the range of the tumor volumes studied, ruling out the prospect of a constant proliferating fraction. A less expected finding was that the logistic model (linear decay in volume of the relative growth rate) was also unable to describe our data, although similar results have been observed in other experimental systems [31], [33]. On the other hand, the Gompertz and power law models could give an accurate and identifiable description of the growth slowdown, for both data sets. More elaborate models such as the generalized logistic, von Bertalanffy, and dynamic CC models could describe them as well. However, their parameters were found not to be identifiable from only tumor growth curves, in the ranges of the observed volumes. Additional data could improve identifiability, such as related to later growth and saturation details. It should be noted in this case that the dynamic CC model was not designed with the intent to quantify tumor growth, but rather to describe the effects of anti-angiogenic agents on global tumor dynamics. Because the model carries angiogenic parameters that are not directly measureable, or even inferable, from the experimental systems we used, it stands to reason that they would not be easily identifiable from the data. Kinetics under the influence of antiangiogenic therapy might thus provide useful additional information that could render this model identifiable. For the breast tumor experimental system, the slowdown was characterized by linear dynamics and was most accurately fitted by the exponential-linear model. Observed was exponential growth from the number of injected cells (during the unobserved phase) that switched smoothly to a linear phase (exponential-linear model). It should be noted that in the breast tumor data set, no data were available during the initiation phase (below 200 mm3) and only the linear part of a putative exponential-linear growth was observed. Explorations of the kinetics of growth during the initial phase (at volumes below the mm3) are needed for further clarification. Despite structural similarities, important differences were noted in the parameter estimates between the two experimental models, in agreement with other studies emphasizing differences between ectopic and orthotopic growth [61], [62]. Our results and methodology may help to identify the impact on kinetics of the site of implantation, although explicit comparisons could not be made here due to the differences in the cell lines used. The Gompertz model (exponential decay in time of the relative growth rate) was able to fit both data sets accurately, consistently with the literature [13], [15], [29], [31], [33]. One of the main criticisms of the Gompertz model is that the relative tumor growth rate becomes arbitrarily large (or equivalently, the tumor doubling time gets arbitrarily small) for small tumor volumes. Without invoking a threshold this becomes biologically unrealistic. This consideration led investigators [13], [63] to introduce the Gomp-exp model that consists in an initial exponential phase followed by Gompertzian growth when the associated doubling time becomes realistic. This approach could also be applied to any decreasing relative growth rate model. We did not consider it in our analysis due to the already large initial volume and the lack of data on the initiation phase where the issue is most relevant. The power law model was also able to describe the experimental data and appeared as a simple, robust, descriptive and predictive mathematical model for murine tumor growth kinetics. It suggests a general law of macroscopic in vivo tumor growth (in the range of the volumes observed): only a subset of the tumor cells proliferate and this subset is characterized by a constant, possibly fractional, Hausdorff dimension. In our results, this dimension (equal to 3γ) was found to be significantly different from two or three (p<0. 05 by Student' s t-test) in 14/20 mice for the lung tumor data set and 13/34 mice for the breast tumor data, effectively suggesting a fractional dimension. A possible explanation of this feature could come from the fractal nature of the tumor vasculature [64], [65], an argument supported by others who have investigated the link between tumor dynamics and vascular architecture [37]. More precisely, the branching nature of the vascularization generates a fractal organization [37], [64], [65] that could in turn produce a contact surface of fractional Hausdorff dimension. Considering further that the fraction of proliferative cells is proportional to this contact surface (for instance because proliferative cells are limited to an area at fixed distance from a blood vessel or capillary, due to diffusion limitations), this could make the connection between fractality of the vasculature and proliferative tissue. These considerations could therefore provide a mechanistic explanation for the growth rate decay that naturally happens when the dimension of the proliferative tissue is lower than three. Our results were obtained using two particular experimental systems: an ectopic mouse syngeneic lung tumor and an orthotopic human xenograft breast tumor model. Although consistent with other studies that found the power law model adequate for growth of a murine mammary cell line [38] or for description of human mammography density distribution data [5], these remain to be confirmed by human data. This model should also be taken with caution when dealing with very small volumes (at the scale of several cells for instance) for which the relative growth rate becomes very large. Indeed, the interpretation of a fractional dimension then fails, since the tumor tissue can no longer be considered a continuous medium. In this instance, it may be more appropriate to consider exponential growth in this phase [37]. Our results showed that a highly descriptive model (associated to large flexibility) such as the generalized logistic model, might not be useful for predictions, while well-adapted rigidity – as provided by the exponential-linear model on the breast tumor data – could lead to very good predictive power. Interestingly, our study revealed that models having low identifiability (von Bertalanffy and dynamic CC) could nevertheless exhibit good predictive power. Indeed, over a limited time span, different parameter sets for a given model could generate the same growth curves, which would be equally predictive. For the Gompertz model, predictive power might be improved by using possible correlations between the two parameters of this model, as reported by others [15], [63], [66]–[68] and suggested by our own parameter estimates (R = 0. 99 for both data sets, results not shown). If a backward prediction is desired (for instance for the identification of the inception time of the tumor), the use of exponential growth might be more adapted for the initial, latency phase, e. g. by employment of the Gomp-exp model [13], [63]. Translating our results to the clinical setting raises the possibility of forecasting solid tumor growth using simple macroscopic models. Use of a priori information could then be a powerful method and one might think of the population distribution of parameters being learned from existing databases of previous patient examinations. However, the very strong improvement of prediction success rates that we obtained partly comes from the important homogeneity of our growth data (in particular the LLC data) that generated a narrow and very informative distribution of some parameters (for instance parameter γ of the power law model), which in turn powerfully assisted the fitting procedure. In more practical situations such as with patient data, more heterogeneity of the growth data should be expected that could alter the benefit of the method. For instance, in some situations, growth could stop for arbitrarily long periods of time. These dormancy phases challenge the universal applicability of a generic growth law such as the Gompertz or power law [69]. Description of such dormancy phenomena could be integrated using stochastic models that would elaborate on the deterministic models reviewed here, as was done by others [70] to describe breast cancer growth data using the Gompertz model. Moreover, further information than just tumor volume could be extracted from (functional) imaging devices, feeding more complex mathematical models that could help design more accurate in silico prediction tools [18], [71]. Our analysis also has implications for the use of mathematical models as valuable tools for helping preclinical anti-cancer research. Such models might be used, for instance, to specifically ascertain drug efficacy in a given animal, by estimating how importantly the treated tumor deviates from its natural course, based on a priori information learned from a control group. Another application can be for rational design of dose and scheduling of anti-cancerous drugs [22], [72], [73]. Although integration of therapy remains to be added (and validated) to models such as the power law, more classical models (exponential-linear [39] or dynamic CC [41]) have begun to predict cytotoxic or anti-angiogenic effects of drugs on tumor growth. Our methods have allowed precise quantification of their respective descriptive and predictive powers, which, in combination with the models' intrinsic biological foundations, could be of value when deciding among such models which best captures the observed growth behaviors in relevant preclinical settings.
Tumor growth curves display relatively simple time curves that can be quantified using mathematical models. Herein we exploited two experimental animal systems to assess the descriptive and predictive power of nine classical tumor growth models. Several goodness-of-fit metrics and a dedicated error model were employed to rank the models for their relative descriptive power. We found that the model with the highest descriptive power was not necessarily the most predictive one. The breast growth curves had a linear profile that allowed good predictability. Conversely, not one of the models was able to accurately predict the lung growth curves when using only a few data points. To overcome this issue, we considered a method that uses the parameter population distribution, informed from a priori knowledge, to estimate the individual parameter vector of an independent growth curve. This method was found to considerably improve the prediction success rates. These findings may benefit preclinical cancer research by identifying models most descriptive of fundamental growth characteristics. Clinical perspective is also offered on what can be expected from mathematical modeling in terms of future growth prediction.
oncology
medicine and health sciences
mathematics
theoretical biology
tumor physiology
calculus
biology and life sciences
physical sciences
basic cancer research
differential equations
Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth