Abstract
Defining tumor stage at diagnosis is a pivotal point for clinical decisions about patient treatment strategies. In this respect, early detection of occult metastasis invisible to current imaging methods would have a major impact on best care and long-term survival. Mathematical models that describe metastatic spreading might estimate the risk of metastasis when no clinical evidence is available. In this study, we adapted a top-down model to make such estimates. The model was constituted by a transport equation describing metastatic growth and endowed with a boundary condition for metastatic emission. Model predictions were compared with experimental results from orthotopic breast tumor xenograft experiments conducted in Nod/Scidγ mice. Primary tumor growth, metastatic spread and growth were monitored by 3D bioluminescence tomography. A tailored computational approach allowed the use of Monolix software for mixed-effects modeling with a partial differential equation model. Primary tumor growth was described best by Bertalanffy, West, and Gompertz models, which involve an initial exponential growth phase. All other tested models were rejected. The best metastatic model involved two parameters describing metastatic spreading and growth, respectively. Visual predictive check, analysis of residuals, and a bootstrap study validated the model. Coefficients of determination were |$R^2 = 0.94$| for primary tumor growth and |$R^2 = 0.57$| for metastatic growth. The data-based model development revealed several biologically significant findings. First, information on both growth and spreading can be obtained from measures of total metastatic burden. Second, the postulated link between primary tumor size and emission rate is validated. Finally, fast growing peritoneal metastases can only be described by such a complex partial differential equation model and not by ordinary differential equation models. This work advances efforts to predict metastatic spreading during the earliest stages of cancer. Cancer Res; 74(22); 6397–407. ©2014 AACR.
This work advances efforts to predict metastatic spreading during the earliest stages of cancer, at points that could help clinicians make the best decisions on treatment strategies.
Primary tumor size |$x_p$| is described by the Gomp-Ex growth equation
The Gomp-Ex adaptation of the Gompertz model constrains the tumor doubling rate to the in vitro doubling rate.
The metastatic burden is given by
where |$\beta$| is the metastatic emission rate from the primary tumor
and |$x_m$| is the size of a metastasis emitted at |$t = 0$|, again described by a Gomp-Ex model:
A metastasis emitted at time|$s$| reaches the size |$x_m \left({t - s} \right)$| at time |$t$|, which yields the structure of the integral of the total metastatic burden |$M\left(t \right)$|. This formulation is derived from a PDE model as shown in the Supplementary Materials.
|$a_p$| and |$a_m$| are the Gompertzian growth rates of the primary tumor and of the metastases, respectively;
|$a_{vitro}$| is the in vitro proliferation rate, which is determined experimentally;
|$x_0$| denotes the primary tumor size at time |$t = 0$|;
|$b$| is the maximal tumor size of both primary and metastatic tumors;
|$\mu$| is the rate of metastatic emission.
In the model, size is measured in number of cells, which is in turn linked to the signal from bioluminescence imaging (measured in photons emitted per second) by an experimentally determined ratio.
The in vivo growth rates |$a_p$| and |$a_m$| cannot exceed the in vitro growth rate |$a_{vitro}$|;
Metastatic emission is proportional to a |$2/3$| power of primary tumor size, which corresponds to surficial metastatic emission;
Secondary metastatic emission is negligible in the experimental framework;
All metastases grow at the same rate;
Newborn metastases consist of a single cell;
The number of viable tumor cells is proportional to the signal measured by bioluminescence imaging.
Interindividual variation is incorporated in a mixed-effects model. Each of the five structural parameters |$a_p, a_m, x_0, b,\mu$| is log-normally distributed among the individuals in the population. The measurements are log-normally distributed around the predicted primary tumor size and metastatic burden.
For each structural parameter, mean and variance of the random variable's logarithm are estimated.
|$\sigma_{p}$| and |$\sigma_m$| are the residual error SDs for the logarithmized primary tumor size and metastatic burden, respectively.
Each individual's structural parameters are realizations of a random variable characterizing the population.
Introduction
Tumor classification as localized or advanced disease is a major issue in clinical oncology because subsequent therapeutic strategies (i.e., adjuvant chemo- and/or radiotherapy) depend on the initial disease staging. In this respect, occult metastasis (i.e., lesions smaller than 2–3 mm not detectable by medical imaging) impedes decisions on appropriate adjuvant therapies aiming at limiting the risk of recurrence. A tool that predicts the metastatic risk for a given patient, despite the lack of imaging-based evidence for secondary lesions, would therefore be clinically relevant by orienting the oncologist toward the most appropriate therapeutic options.
Along with the development of a personalized medicine (1), mathematical modeling has been constantly on the rise in the field of medical sciences. Tumor growth dynamics have received considerable attention and numerous models have been proposed to describe growth on different scales. On a purely macroscopic level, the Gompertz model, exhibiting a sigmoidal growth curve, has been used first (2). Many other sigmoidal growth curves have also been proposed (3, 4). Some of these models are derived on a mechanistic basis, such as the von Bertalanffy (5) and West (6, 7) models, which assume an auto-similar tissue organization. A model correction limiting the growth rate of small tumors also exists (3).
To describe dynamic processes such as the spatial evolution of tumors, more complex strategies are used, for example, models based on reaction-diffusion equations (8). Most of the recent approaches aim at integrating newly discovered factors and determinants implicated into tumor growth such as stem cells (9), microenvironment (10), or tumor-induced angiogenesis (11). In addition, models have been developed to describe the impact of cancer therapy onto tumor growth so as to predict efficacy (12, 13).
Because metastasis is a highly complex, multistep process depending on many determinants, modeling metastatic spreading is a particularly challenging task. Compared with the numerous models dedicated to tumor growth, relatively few models that describe metastasis have been proposed so far. The existing approaches range from simple phenomenological to complex mechanistic models and contain deterministic and/or stochastic components (14).
Modeling complexity should depend on the purpose of the model. Phenomenological approaches are particularly suitable for models conceived as prognostic tools. These models intrinsically allow for an easier confrontation with preclinical and clinical data, which is essential for model validation.
The metastatic model introduced by Iwata and colleagues accurately described the metastatic growth dynamics visible on CT scans of a single patient with a metastatic hepatocellular carcinoma (15). It was further fitted to clinical data collected over a 20-year follow-up study for metastatic relapse in 3,500 patients with breast cancer and accurately predicted the metastatic risk given the primary tumor size (16). Despite these promising results, a more thorough validation and a model adaptation to more general settings seem necessary. Indeed, the former study includes only one patient and growth and emission take place in a single organ. In the latter study, the model is fitted to the percentage of patients with a metastatic disease, but no information on growth dynamics is available.
Recording the necessary clinical data for a thorough model validation would be challenging, due to ethical constraints and the long duration of metastatic development. Before this, a suitable preclinical approach, such as the MDA-MB-231 orthotopic xenograft in mice (17, 18), could allow the record of a large dataset within a reasonable delay.
Nonlinear mixed-effects (NLME) modeling is a statistical framework for the description of repeated measurements from a population of individuals. Characterizing both typical population values and individual variation around these values, this approach efficiently uses information from sparse longitudinal data. NLME modeling has been widely used in clinical and preclinical settings and in particular for the description of tumor growth data (19). Several algorithms for maximum likelihood estimation of NLME models have been developed. Among them, the SAEM algorithm (20, 21), first implemented in the Monolix software, has recently received considerable attention, especially in the pharmacology community. In practice, NLME modeling is almost exclusively used for models written as systems of ordinary differential equations (ODE). To our knowledge, the only approach stretching this framework to partial differential equations (PDE) is due to Grenier and colleagues (22).
In this work, we obtained detailed experimental data on tumor growth and metastatic spreading from tumor-bearing mice using noninvasive 3D bioluminescence imaging. The metastatic model described above was adapted to this nonclinical setting and fitted to the collected experimental data within an NLME modeling framework. To this end, a tailored computational approach for the resolution of the PDE model was used in the Monolix software.
Materials and Methods
Cell culture
The human breast cancer cell line MDA-MB-231-luc-D3H2LN was purchased from Caliper Life Sciences (Villebon sur Yvette; product number 119 261). Cells were used within 6 months upon reception from the supplier. This BioWare light producing cell line is derived from the MDA-231 human adenocarcinoma by stable transfection of the North American firefly gene expressed from SV40 promoter. Tissue origin was a mammary gland adenocarcinoma from the MD Anderson Cancer Center (Houston, TX). Upon reception, cells were gently thawed then quickly cultured following supplier's recommendation until grafted as described below. Cells were immediately amplified in T150 culture flasks and frozen into nitrogen when in exponential phase. Next, cells were used as soon as a sufficient number was reached for grafting with never more than two passages prior engraftment. The cell line was regularly authenticated on the basis of viability, recovery, growth, morphology, and bioluminescence. Cells in exponential phase were grown following standard recommendations in RPMI-1640 medium (Invitrogen) supplemented with 2 mmol/L l-glutamine (Invitrogen), 5 IU/mL penicillin/streptomycin (Eurobio), 5 IU/mL fungizon (Eurobio), and 10% of FCS (Eurobio) at 37°C in a humidified atmosphere with 5% CO2.
Animal experiments
All experiments were approved by an ethical committee and were conducted in compliance with European regulations, based on the UKCCCR guidelines for the welfare of animals in experimental neoplasia (23). Pathogen-free 6-week-old female Nod Scid γ (NSG) mice (Charles River laboratories) were acclimatized in a sterile environment for 2 weeks before starting the experiments. Mice were maintained in sterilized filter-stopped cages throughout the experiments, kept in a sterile and thermostated cabinet. They were daily monitored for signs of distress, decreased physical activity and weight. Animals were euthanized under anesthesia when showing signs of distress, cachexia (i.e., loss of 10% of body weight), or when tumor reached an apparent mass of 2 g (i.e., 2 cm3). Bioluminescence measurements were performed twice a week. Carcass bodyweight was monitored twice a week as a standard marker for quality of life in animals.
Preparation of cell suspensions and orthotopic cell implantation
Before implantation, the MDA-MB-231 cells were trypsinized, counted, and centrifuged (5 minutes, 2,500 rpm) and washed twice in sterile PBS. Cells were resuspended in RPMI-1640 with 60% of Matrigel (BD Sciences) and maintained in ice-cooled conditions until engraftment. A volume of 50 μL containing |$1.5 \times 10^5$| cells was injected in the mammary inguinal right fat pad through the nipple under gas anesthesia [2% isofluran, (Abbott France) in standard oxygen-nitrogen protoxide mix].
Bioluminescence imaging
Imaging and growth monitoring for both primary tumors and metastases started on day 6 postinjection. The experiment was carried out until mice euthanasia for ethical reasons. Intraperitoneal injection of firefly D-Luciferin (Caliper Life Sciences) at a dose of 150 mg/kg body weight was performed in mice before starting imaging. Acquisitions started 12 minutes after Luciferine injection, the delay required to reach a plateau in light emission. For 3D bioluminescence, images were captured at six different wavelengths ranging from 560 nm to 660 nm (step: 20 nm). Imaging and data processing were performed in an IVIS Spectrum imager equipped with the Living Image 4.2 software (PerkinElmer). Signal classification as primary or metastatic tumor was supported by the entire observation time history and, in unclear cases, by autopsy.
Photons emitted per second by one cell
The photon-to-cell ratio was determined before starting the in vivo experiments by a first calibration step of the imager, following standard procedure. MDA-MB-231-luc-D3H2LN cells ranging from |$10^2$| to |$10^6$| cells were seeded in a 24-well plate and were immediately imaged. The cell-to-signal ratio was calculated by linear regression.
Models for tumor growth and metastatic spreading
For the description of primary tumor growth, a set of phenomenological candidate models reported in the literature was considered (see Table 1), having the general form
Phenomenological tumor growth models
| Growth model | Growth rate | Reference |
|---|---|---|
| Exponential | |$g\left(x \right) = ax$| | |
| Gompertz | |$g\left(x \right) = ax \ {\rm log}\left({\frac{b}{x}} \right)$|, | (2) |
| Logistic | |$g\left(x \right) = ax\left(1 - {\frac{x}{b}}\right)$| | (4) |
| Power growth | |$g\left(x \right) = ax^c$| | (4) |
| Von Bertalanffy | |$g\left(x \right) = ax^{2/3} - bx$| | (4) |
| West | |$g\left(x \right) = ax^{3/4} - bx$| | (6) |
| Models including the in vitro doubling rate | ||
| Gomp-Ex model | |$g\left(x \right) = \min \left({a_{vitro} \ x, ax \ {\rm log}\left({\frac{b}{x}} \right)} \right)$| | (3) |
| West-Ex model | |$g\left(x \right) = \min \left({a_{vitro} \ x, \ ax^{3/4} - bx} \right)$| | Analogy to (3) |
| Berta-Ex model | |$g\left(x \right) = \min \left(a_{vitro} \ x, \ ax^{2/3} - bx\right)$| | Analogy to (3) |
| Growth model | Growth rate | Reference |
|---|---|---|
| Exponential | |$g\left(x \right) = ax$| | |
| Gompertz | |$g\left(x \right) = ax \ {\rm log}\left({\frac{b}{x}} \right)$|, | (2) |
| Logistic | |$g\left(x \right) = ax\left(1 - {\frac{x}{b}}\right)$| | (4) |
| Power growth | |$g\left(x \right) = ax^c$| | (4) |
| Von Bertalanffy | |$g\left(x \right) = ax^{2/3} - bx$| | (4) |
| West | |$g\left(x \right) = ax^{3/4} - bx$| | (6) |
| Models including the in vitro doubling rate | ||
| Gomp-Ex model | |$g\left(x \right) = \min \left({a_{vitro} \ x, ax \ {\rm log}\left({\frac{b}{x}} \right)} \right)$| | (3) |
| West-Ex model | |$g\left(x \right) = \min \left({a_{vitro} \ x, \ ax^{3/4} - bx} \right)$| | Analogy to (3) |
| Berta-Ex model | |$g\left(x \right) = \min \left(a_{vitro} \ x, \ ax^{2/3} - bx\right)$| | Analogy to (3) |
In several models, the tumor growth rate is limited by an experimentally determined in vitro growth rate |$a_{vitro}$|. This has been described for the Gompertz model (3), and the other models were taken by analogy.
The metastatic process was modeled by a transport equation describing the evolution of the size-structured metastatic density ρ
endowed with a boundary condition for metastatic emission:
and the initial condition |$\rho \left({x,0} \right) = 0$| for |$x \in \left(1,{\rm b} \right)$|. In this model, from now on called size-structured model, it is assumed that
newly created metastases consist of a single cell,
both primary tumor and metastases create new metastases at the size-dependent rate |$\beta \left(x \right) = \mu x^\alpha$|,
primary tumor size |$x_p$| is described by the primary tumor growth function |$g_p$|,
metastases grow at a rate |$g_m$|, which could be different from that for primary tumor growth.
For small |$\Delta x$|, the expression |$\rho \left({x, t} \right)\Delta x$| can be interpreted as the number of metastases with size between |$x$| and |$x \ + \ \Delta x$| cells, observed at time |$t$|. The maximal size |$b$| metastasis can attain depends on the metastatic growth rate |$g_m$|. The parameter |$\alpha$| corresponds to an emission proportional to a fractal dimension of the emitting tumor and was estimated close to 2/3 by Iwata and colleagues, which they interpreted as surficial emission (15). Mathematical details on the model are given as Supplementary Material.
The total number of metastases at time |$t$| is given by |$N\left(t \right) = \mathop \smallint \int_1^b \rho \left({x, t} \right)dx$| and the total metastatic burden by |$M\left(t \right) = \mathop \smallint \int_1^b x\rho \left({x, t} \right)dx$|. The primary tumor size |$x_p$| and the total metastatic burden are observed by bioluminescence.
This model was first introduced by Iwata and colleagues (15), with |$x_0 = 1$| and Gompertzian growth |$g\left(x \right) = ax\log \left(b/x\right)$| for primary and secondary tumors. Subsequently, different methods for the numerical resolution of the model were proposed in the literature (24, 25). Here, we used an approach consisting of a model reformulation into a Volterra integral equation, which can be solved efficiently (26). The equivalence of the two formulations is shown in the Supplementary Material. In this setting, the metastatic burden is given by
The first right-hand term describes metastases emitted from the primary tumor and the second term describes secondary metastatic emission. The function |$x_m$| denotes the size of a metastasis emitted at |$t = 0$|:
Statistical model and software
The typical population values and the interindividual variability of each model parameter were the parameters of the NLME model; they were estimated with the Monolix software. During model building, data on primary tumor size and metastatic burden were analyzed sequentially. In the final model, primary tumor and metastatic parameters were estimated simultaneously.
Model selection was primarily based on Akaike's information criterion (AIC). The condition number of the Fisher Information Matrix [denoted by cond(FIM)] was considered as a measure of parameter identifiability. In addition to that, two biologic constraints were also considered, namely the positivity of |$\alpha$| and a tumor doubling time superior to an experimentally determined in vitro doubling time.
The parameters |$a_p, b, x_0, a_m, \mu$| of the final model were assumed to be log-normally distributed among the individuals in the population. For implementation, the parameter |$\mu$| was log transformed and a normal distribution was chosen, yielding a linear dependence of the final structural model on this parameter. For the residual error, a normal distribution with constant variance of the log-transformed data was chosen for primary and secondary tumors.
Model validation
For the final model, model validation consisted of analysis of individual-weighted residuals (IWRES) for normality, graphical analysis of visual predictive checks (VPC), and bootstrapping. We generated 1,000 bootstrap samples, took the 95% confidence intervals of the parameter estimates in the bootstrapped samples and checked whether the estimated typical values were within this interval. Data were stratified to account for different degrees of data richness (27). The |$R^2$| coefficient of determination was calculated for the individual primary tumor and metastatic fits.
Results
Setting optimal experimental conditions
In vitro doubling time was estimated to 34 hours (data not shown). In preliminary experiments, signal-to-cell ratio was evaluated by measuring the number of photons emitted from the MDA-MB-231 line. The determined ratio was |$88.37$| photons/cell/second. Next, in vivo acquisitions were performed at different times after 150 mg/kg luciferine injection on satellite tumor-bearing mice to identify the plateau ensuring a maximum signal-to-noise ratio. The emitted signal was maximal 10 minutes after injection, with a plateau lasting up to 15 minutes (data not shown). Consequently, it was decided to perform all subsequent acquisitions 12 minutes after luciferine injection for optimal sensitivity. Finally, accuracy of the 3D reconstruction and localization of the signal-emitting sources were checked on a group of satellite mice by comparing in silico reconstructed images with necropsy. Excellent accuracy was obtained, both for localization and signal-to-mass ratio (Fig. 1). The site of each signal could be localized consistently and even metastatic signals very close to the primary tumor were recognized as distinct sources. Repeated reconstructions led to consistent estimates of the signal intensity per metastatic site (error rate inferior to 5%).
Correlation of localization on reconstructed images and necropsy. Red (big ellipse), primary tumor and peritoneal metastases. Yellow (rectangle), hepatic metastases. Orange (polygon), pulmonary metastases. Blue (small ellipse), axillary metastases.
Correlation of localization on reconstructed images and necropsy. Red (big ellipse), primary tumor and peritoneal metastases. Yellow (rectangle), hepatic metastases. Orange (polygon), pulmonary metastases. Blue (small ellipse), axillary metastases.
3D monitoring of metastatic spreading in tumor-bearing mice
Tumor injection yielded a 100% engraftment rate in NSG mice. Figure 2 displays typical 3D reconstructions of a representative mouse for the whole observation period. Diffuse light imaging tomography allowed signal detection from |$7.5 \times 10^5$| cells. The peritoneum was first colonized by metastatic cells, beginning on day 23. All mice eventually metastasized in this area. Lung metastases were first observed on day 30 and 50% of the mice eventually developed secondary pulmonary lesions. First, lymph node metastases were observed on day 34 and 25% of the mice showed such lesions.
Measurements on primary tumor size ranged from |$6.6 \times 10^7$| photons/second to |$8.3 \times 10^{10}$| photons/second (3.1 orders of magnitude) and measurements on the metastatic burden ranged from |$1.5 \times 10^8$| to |$4.5 \times 10^{10}$| photons/second (2.5 orders of magnitude).
Primary tumor model
The selection criteria for the primary tumor models supplied by Monolix are shown in Table 2. According to AIC, the Berta-Ex, West-Ex, West, Gomp-Ex, Gompertz, and Bertalanffy models were the best, whereas the other models (Exponential, Logistic, Power growth) were clearly rejected.
AIC, model estimates of the initial number of cells |$x_0$| and model-derived doubling times for one cell and for |$x_0$|
| Model | AIC | Estimated |${\bi x}_{\bf 0}$| | DT (one cell) | DT (|${\bi x}_{\bf 0}|) |
|---|---|---|---|---|
| Exponential | 387.1 | 1.7 × 106 | 109 h | 109 h |
| Gompertz | 262.0 | 6.2 × 104 | 11 h | 25 h |
| Logistic | 303.1 | 5.0 × 105 | 65 h | 65 h |
| Power growth | 318.2 | 3.3 × 105 | 4 h | 41 h |
| Bertalanffy | 263.6 | 3.6 × 10−9 | 16 min | 2 sec |
| West | 258.8 | 26.3 | 43 min | 2 h |
| Gomp-Ex | 260.3 | 9.0 × 104 | 34 h (imposed) | 34 h (imposed) |
| West-Ex | 256.5 | 8.0 × 104 | 34 h (imposed) | 34 h (imposed) |
| Berta-Ex | 255.5 | 7.6 × 104 | 34 h (imposed) | 34 h (imposed) |
| Model | AIC | Estimated |${\bi x}_{\bf 0}$| | DT (one cell) | DT (|${\bi x}_{\bf 0}|) |
|---|---|---|---|---|
| Exponential | 387.1 | 1.7 × 106 | 109 h | 109 h |
| Gompertz | 262.0 | 6.2 × 104 | 11 h | 25 h |
| Logistic | 303.1 | 5.0 × 105 | 65 h | 65 h |
| Power growth | 318.2 | 3.3 × 105 | 4 h | 41 h |
| Bertalanffy | 263.6 | 3.6 × 10−9 | 16 min | 2 sec |
| West | 258.8 | 26.3 | 43 min | 2 h |
| Gomp-Ex | 260.3 | 9.0 × 104 | 34 h (imposed) | 34 h (imposed) |
| West-Ex | 256.5 | 8.0 × 104 | 34 h (imposed) | 34 h (imposed) |
| Berta-Ex | 255.5 | 7.6 × 104 | 34 h (imposed) | 34 h (imposed) |
NOTE: The Gomp-Ex, West-Ex, and Berta-Ex doubling times are limited by the in vitro doubling time of 34 hours.
Abbreviation: DT, doubling time.
At the beginning of the experiment (outside of the observation range), the Bertalanffy and West models predicted very low initial tumor sizes and very short doubling times. To a lesser degree, this was also the case of the Gompertz model with a doubling time below the in vitro value (24 hours vs. 34 hours). The model versions limiting the doubling rate by the in vitro value (Berta-Ex, West-Ex, and Gomp-Ex) extrapolated to small sizes correctly and even improved the fit. The Berta-Ex and West-Ex fits were equivalent and slightly better than the Gomp-Ex fit.
Metastatic model
First of all, the inclusion of secondary metastatic emission in the size-structured model did not alter the results. For faster simulations, it was thus neglected.
Several combinations of free parameters in the primary and metastatic Gomp-Ex growth models and the emission model |$\beta$| were then compared. When the same parameters were assumed for |$g_p$| and |$g_m$|, the emission parameter |$\alpha$| estimate was unreliable (negative value) and the model was rejected. Equally, AIC clearly rejected all models where |$g_p = g_m$|. It was thus necessary to consider metastatic growth rates different from the primary tumor growth rate.
The best structure of the metastatic model with Gomp-Ex growth was then determined (Table 3, top). Models with free |$\alpha$| or with only free |$\mu$| were mis-specified. The model with free|$\mu, a_m, b_m$| led to inconsistent estimation. Models with one free growth parameter and one emission parameter (|$\mu, a_m$| or |$\mu, b_m$|) described the data best. Because of lacking parameter identifiability of the latter and model inconsistency in the bootstrap procedure (data not shown), the model |$\mu, a_m$| was selected.
Determination of the optimal structure of the size-structured model and maximum likelihood estimates of independent metastatic growth models
| Metastatic model (size-structured) | ||||
|---|---|---|---|---|
| Number of parameters | Estimated parameters | Fixed parameters | AIC | cond(FIM) |
| 1 | |$\mu$| | |$\alpha = 2/3, a_m = a_p, b_m = b_p$| | 246.0 | 2.5 |
| 2 | |$\mu, \alpha$| | |$a_m = a_p, b_m = b_p$| | 241.7 | 4,800 |
| 2 | |$\mu, a_m$| | |$\alpha = 2/3,b_m = b_p$| | 183.1 | 80 |
| 2 | |$\mu, b_m$| | |$\alpha = 2/3,a_m = a_p$| | 180.2 | 490 |
| 3 | |$\mu, a_m, b_m$| | |$\alpha = 2/3$| | 189.9 | 690 |
| Independent growth of a metastasis emitted at |${\bi t}_{\bf 0}$| with size 1 | ||||
| Model name | Estimated | Fixed | AIC | DT 1 cell |
| Gompertz | |$a_m = 0.04,t_0 = - 3.7$| | |$b_m = b_p$| | 187.6 | 17h |
| Gomp-Ex | |$a_m = 0.048,t_0 = - 10.0$| | |$b_m = b_p$| | 187.7 | 34h |
| Metastatic model (size-structured) | ||||
|---|---|---|---|---|
| Number of parameters | Estimated parameters | Fixed parameters | AIC | cond(FIM) |
| 1 | |$\mu$| | |$\alpha = 2/3, a_m = a_p, b_m = b_p$| | 246.0 | 2.5 |
| 2 | |$\mu, \alpha$| | |$a_m = a_p, b_m = b_p$| | 241.7 | 4,800 |
| 2 | |$\mu, a_m$| | |$\alpha = 2/3,b_m = b_p$| | 183.1 | 80 |
| 2 | |$\mu, b_m$| | |$\alpha = 2/3,a_m = a_p$| | 180.2 | 490 |
| 3 | |$\mu, a_m, b_m$| | |$\alpha = 2/3$| | 189.9 | 690 |
| Independent growth of a metastasis emitted at |${\bi t}_{\bf 0}$| with size 1 | ||||
| Model name | Estimated | Fixed | AIC | DT 1 cell |
| Gompertz | |$a_m = 0.04,t_0 = - 3.7$| | |$b_m = b_p$| | 187.6 | 17h |
| Gomp-Ex | |$a_m = 0.048,t_0 = - 10.0$| | |$b_m = b_p$| | 187.7 | 34h |
NOTE: Top, determination of the optimal structure of the size-structured model, supposing both primary tumor and metastatic Gomp-Ex growth; bottom, maximum likelihood estimates (estimated), AIC, and doubling time of one cell (DT 1 cell) of independent metastatic growth models.
Abbreviation: DT, doubling time.
Next, the three best-fitting primary tumor growth functions (Berta-Ex, West-Ex, and Gomp-Ex) were compared as metastatic growth functions in the size-structured model. The Gomp-Ex model fitted the data on metastatic burden slightly better than the other models (West-Ex: |${\rm \Delta}AIC = 2.6$| and Berta-Ex: |${\rm \Delta}AIC = 3.5$|). This model was selected for primary tumor and metastatic growth, preserving the same parametric structure for both growth processes. The final structure of the metastatic model was
with a primary tumor regression parameter |$b$| and two free parameters |$\mu$| and |$a_m$|.
To validate the direct link between primary tumor size and metastatic emission, parameters were estimated in two settings. In both, the primary tumor parameters were first estimated. Then, the metastatic parameters were estimated using either the post hoc individual primary tumor parameters or the typical population values. The comparison of the two settings showed that individual differences in primary tumor size highly impact on the goodness of fit|$\left({\rm \Delta\ AIC = 3.7} \right)$|.
This size-structured model was then compared with ODE models describing the metastatic burden as one tumor emitted by the primary tumor at a given time point|$t_0$|and growing independently (Table 3, bottom). Gompertzian growth yielded low doubling times and a metastatic inception before primary tumor xenograft (|$t_0 = - 3.7$| day). A Gomp-Ex model placed the metastatic inception even earlier (|$t_0 = - 10.0$| day). In addition, goodness of fit of both models was inferior to that of the size-structured model (|${\rm \Delta}AIC = 4.5$| and |${\rm \Delta}AIC = 4.6$|).
Population parameter estimates with FIM and bootstrap-computed confidence intervals are shown in Table 4. Parameter estimates were close to the bootstrap estimates and inside the 95% confidence intervals computed by bootstrapping. Residuals of the primary tumor and metastatic models were approximately normally distributed (Fig. 3). A VPC validated both models because all empirical percentiles (10%, 50%, and 90%) were within the theoretical 95% confidence intervals (Fig. 3).
Monolix outputs. Left, primary tumor fit (Gomp-Ex model). Right, metastatic model fit (Gomp-Ex growth). Top, IWRES assess whether the difference between individual fits and observations is consistent with random error. Residuals are approximately normally distributed. Bottom, VPC compares observations and simulations of model predictions visually. All empirical percentiles (10%, 50%, and 90%) are within the corresponding 95% VPC confidence intervals.
Monolix outputs. Left, primary tumor fit (Gomp-Ex model). Right, metastatic model fit (Gomp-Ex growth). Top, IWRES assess whether the difference between individual fits and observations is consistent with random error. Residuals are approximately normally distributed. Bottom, VPC compares observations and simulations of model predictions visually. All empirical percentiles (10%, 50%, and 90%) are within the corresponding 95% VPC confidence intervals.
Population parameter estimates with FIM and bootstrap-computed confidence intervals
| Primary tumor parameters (U) | Typical population value (|$ \pm$|SE) | Bootstrap estimate (95% CI) | Interindividual variability (|$ \pm$|SE) |
|---|---|---|---|
| |$a_p$| (1/day) | |$0.075 \ \left({\pm 0.006} \right)$| | |$0.075 \ \left({0.063,0.088} \right)$| | |$0.11 \ \left({\pm 0.048} \right)$| |
| |$b$| (cells) | |$5.4 \times 10^8 \ \left({\pm 8.8 \times 10^7} \right)$| | |$5.5 \times 10^8 {\rm} \ \left({4.0 \times 10^8, 7.6 \times 10^8} \right)$| | |$0.23 \ \left({\pm 0.14} \right)$| |
| |$x_0$| (cells) | |$9.0 \times 10^4 \ \left({\pm 1.3 \times 10^4} \right)$| | |$9.4 \times 10^4 \ \left({6.9 \times 10^4, 1.3 \times 10^5} \right)$| | |$0.18 \ \left({\pm 0.21} \right)$| |
| Residual error |$\sigma _p$| | |$0.47\left \ ({\pm 0.03} \right)$| | — | — |
| Metastatic parameters (U) | |||
| |$\log \left(\mu \right)\ \left(\mu: 1/({\rm cells}^{2/3}\ {\rm day}\right)$|: | −0.46 (±0.7) | −0.45 (−1.98, +0,80) | 0.21 (±0.74) |
| |$a_m$| (1/day) | |$7.9 \times 10^{- 3} \ {\rm}\left({\pm 2.5 \times 10^{- 3}} \right)$| | |$7.6 \times 10^{- 3} \ {\rm}\left({3.5 \times 10^{- 3}, 1.3 \times 10^{- 2}} \right)$| | |$0.11 \ \left({\pm 0.28} \right)$| |
| Residual error |$\sigma _m$| | |$0.90 \ \left({\pm 0.09} \right)$| | − | − |
| Primary tumor parameters (U) | Typical population value (|$ \pm$|SE) | Bootstrap estimate (95% CI) | Interindividual variability (|$ \pm$|SE) |
|---|---|---|---|
| |$a_p$| (1/day) | |$0.075 \ \left({\pm 0.006} \right)$| | |$0.075 \ \left({0.063,0.088} \right)$| | |$0.11 \ \left({\pm 0.048} \right)$| |
| |$b$| (cells) | |$5.4 \times 10^8 \ \left({\pm 8.8 \times 10^7} \right)$| | |$5.5 \times 10^8 {\rm} \ \left({4.0 \times 10^8, 7.6 \times 10^8} \right)$| | |$0.23 \ \left({\pm 0.14} \right)$| |
| |$x_0$| (cells) | |$9.0 \times 10^4 \ \left({\pm 1.3 \times 10^4} \right)$| | |$9.4 \times 10^4 \ \left({6.9 \times 10^4, 1.3 \times 10^5} \right)$| | |$0.18 \ \left({\pm 0.21} \right)$| |
| Residual error |$\sigma _p$| | |$0.47\left \ ({\pm 0.03} \right)$| | — | — |
| Metastatic parameters (U) | |||
| |$\log \left(\mu \right)\ \left(\mu: 1/({\rm cells}^{2/3}\ {\rm day}\right)$|: | −0.46 (±0.7) | −0.45 (−1.98, +0,80) | 0.21 (±0.74) |
| |$a_m$| (1/day) | |$7.9 \times 10^{- 3} \ {\rm}\left({\pm 2.5 \times 10^{- 3}} \right)$| | |$7.6 \times 10^{- 3} \ {\rm}\left({3.5 \times 10^{- 3}, 1.3 \times 10^{- 2}} \right)$| | |$0.11 \ \left({\pm 0.28} \right)$| |
| Residual error |$\sigma _m$| | |$0.90 \ \left({\pm 0.09} \right)$| | − | − |
NOTE: First column, typical population values are the median values of the parameter distributions. SE are given in parenthesis; second column, bootstrap estimates of the typical population values with 95% bootstrap confidence intervals (95% CI); third column, depending on the parameter distribution, the inter-individual variability is the parameter variance (normally distributed parameter) or the variance of the logarithmized parameters (log-normally distributed). SEs are given in parenthesis.
The |$R^2$| coefficient of determination was calculated for logarithmized primary tumor size and logarithmized metastatic burden. The individual primary tumor fits yielded |$R^2 = 0.94$| and the individual metastatic burden fits |$R^2 = 0.57$|.
Discussion
An operational diagnostic model for metastasis could assess the risk of occult metastasis and predict growth when no information is available upon medical imaging. Here, we have validated a phenomenological model of the metastatic process in a combined experimental and modeling approach.
Tumor heterogeneity and intertumor diversity result in an enhanced complexity (28), which the most sophisticated mechanistic models fail to match. For an adequate description of sparse experimental and clinical data on tumor growth and metastatic emission, only phenomenological approaches summing up most underlying mechanisms are conceivable.
Human MDA-MB-231 breast cancer cells were selected as a substrate widely used in experimental therapeutics with breast cancer, especially for addressing cell migration and distant metastasis (29). Bioluminescence tomography was chosen as a noninvasive method for monitoring accurately both primary tumor growth and metastatic spreading (30) and for detecting deep lesions as small as |$7.5 \times 10^5$| cells. The primary tumor was not resected during the experiment because of the potential impact on spreading (31). NSG mice were chosen because of their high level of immunodeficiency, ensuring a maximal metastatic spreading. Preliminary attempts to generate distant metastasis with other animals (e.g., NMRI or Swiss nu/nu mice) grafted with the MDA-MB-231 line had a high failure rate. Here, all animals generated secondary lesions within 3 to 4 weeks postengraftment, mainly in the peritoneal cavity. Secondary sites were the lungs and lymph nodes. Of note, no bone metastases were observed in our experiments, although this lesion is frequently observed in patients with breast cancer.
With 166 data records on 16 mice, ranging over three orders of magnitude, and a moderate noise level, our primary tumor data were an excellent benchmark for tumor growth models. Among six growth models reported in the literature, three were found to describe our data significantly better than the others, namely the Bertalanffy, West, and Gompertz models. Because of the close resemblance of three-parameter tumor growth models, some authors have cast doubts on their distinguishability (4, 32). It is thus remarkable that our data permitted this distinction.
As a subsequent model validation step, doubling times were estimated for the initial predicted size|$x_0$|. For all three models, doubling times were estimated below the experimentally determined in vitro doubling time, which is biologically aberrant. This has been reported previously for the Gompertz model (11), but we found that the effect was much more pronounced for the Bertalanffy and West models, which missed the doubling time and the initial size |$x_0$| by several orders of magnitude. To overcome this biologically unrealistic model behavior, hybrid models (Berta-Ex, West-Ex, and Gomp-Ex) featuring a growth speed limitation by the in vitro growth speed were considered. These hybrid models also improved the model fits.
In the three final primary tumor growth models, the initial tumor size was estimated to |$9.0 \times 10^4$|, |$8.0 \times 10^4$|, and |$7.6 \times 10^4$|, respectively. Because approximately |$1.5 \times 10^5$| cells were initially injected, this means that approximately 40% to 50% of cells were lost during xenograft, a finding consistent with the steps (washing, trypsinization, washing, suspension in Matrigel-based matrix, injection) the cells underwent before being grafted. This finding corroborates the in vivo validity of the signal-to-cell ratio, which was determined in vitro.
The original metastatic model by Iwata and colleagues had to be adapted and different growth functions were needed to describe primary and secondary tumor growth. The data on total metastatic burden were best described by a Gomp-Ex growth model with two free parameters, one describing metastatic emission and the other metastatic growth.
The final metastatic model structure was corroborated by VPC, analysis of residuals, and bootstrapping. The |$R^2$| values were excellent for the main tumor and satisfactory for the metastases. Several factors could explain why |$R^2$| values were lower for the total metastatic burden than for the primary tumor size. The signal of early metastases was lower than primary tumor signal, effectively making the primary tumor a source of noise for close peritoneal metastases. The anatomy of the lungs also degraded the signal localization. Furthermore, the metastatic development follows different dynamics in each organ, which was lost when summing up the site-specific metastatic burdens to a global quantity. Nevertheless, in contrast with the number of metastases and site-specific burdens, the total metastatic burden had the advantage of robustness against misclassification between different metastases.
Our population-based modeling approach revealed that individual differences in metastatic development cannot be explained by differences in metastatic growth and emission rates only. Individual differences in primary tumor size also have a considerable impact on spreading.
Because most of the data on metastatic burden come from one or two large metastases in the peritoneum, the question arose whether the size-structured model describes this situation better than a simple ODE growth model for one peritoneal metastasis. Most interestingly, ODE models fitted to the peritoneal data yield both unreasonably short estimates of small-size doubling times and a metastatic inception before primary tumor xenograft. In other words, the peritoneal metastases grow too fast to be described by ODE models with biologically reasonable parameters. The size-structured model provides an explanation to fast peritoneal growth, predicting rapid metastatic emission and large numbers of metastases in this site. Because of the high degree of immunodeficiency of the NSG mice, rapid metastatic emission was certainly plausible in our preclinical experiments. The simultaneous growth at biologically reasonable rates could have led to a very quick development of metastases that need not be distinguishable macroscopically.
The case of the peritoneum shows that the apparent number of metastases reveals little about the underlying growth and emission processes and that little insight would be gained by confronting these predictions to the number of metastases visible on images or at necropsy. The modeled number of metastases should rather be understood as the number of metastatic foci, meaning that the fast growth of the peritoneal metastasis is explained by a large number of foci.
Conclusion
Developing in silico tools to better forecast the evolution of cancer is a rising trend in experimental and clinical oncology. In this work, the confrontation of a model for tumor growth and metastatic spreading to experimental data yielded biologically significant findings. First, parameters describing both metastatic growth and metastatic spreading can be identified from measures of the total metastatic burden. Second, our data validate a major modeling hypothesis, the structural link between metastatic emission and primary tumor size. Finally, fast peritoneal growth is accurately described by our metastatic model but not by simpler phenomenological models. Although the best parametric structure is data dependent and could be different if more data were available, the predictions confirm the validity of the general modeling approach.
In a clinical setting, metastatic growth and emission parameters would have be inferred from information available at diagnosis, which could consist of oncogenic biomarkers and measurements of circulating tumor cells or tumor DNA (33). This type of information could be naturally included in the metastatic model via the metastatic growth and emission rates. Kinetic information on metastatic spreading would also facilitate the experimental validation of biologic hypotheses. It has been previously described how measurements of circulating tumor cells would permit to determine the role of secondary metastatic emission (34). Beyond that, a multisite metastatic model could map kinetic information onto the heterogeneity of cancer in different organs.
Optimization of cancer therapy is another important field of use of mathematical modeling. This includes optimizing dosing regimen through better understanding of PK/PD relationships (35). Several extensions of the metastatic model including a treatment have been proposed (36–38). The next step of this work should be to confront model predictions with data in a therapeutic framework. Furthermore, once parameter values are determined and validated, the impact of a variety of treatments could be tested in silico so as to identify the best design.
A model able to predict metastatic growth could give a more precise description of the metastatic state of a patient than the TXNXMX-staging score commonly used in the clinical practice. On the basis of model predictions of the number of metastases and metastatic burden, the metastatic model discussed here could yield a score indicating the risk of occult metastasis during diagnosis or metastatic relapse after surgery. This supplementary information could help the clinician with decisions on treatment strategies.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: S. Mollard, D. Barbolosi, A. Benabdallah, J. Ciccolini, F. Hubert
Development of methodology: N. Hartung, S. Mollard, D. Barbolosi, A. Benabdallah, G. Chapuisat, J. Ciccolini, C. Faivre, F. Hubert
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Mollard, S. Giacommetti, J. Ciccolini
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): N. Hartung, S. Mollard, D. Barbolosi, G. Chapuisat, S. Giacommetti, A. Iliadis, J. Ciccolini, C. Faivre, F. Hubert
Writing, review, and/or revision of the manuscript: N. Hartung, S. Mollard, D. Barbolosi, G. Chapuisat, A. Iliadis, J. Ciccolini, C. Faivre, F. Hubert
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Mollard, G. Henry, F. Hubert
Study supervision: D. Barbolosi, F. Hubert
Other (analysis of the mathematical model): A. Benabdallah
Grant Support
The authors were partially supported by the French Agence Nationale de la Recherche under grant MEMOREX-PK ANR-09-BLAN-0217-01, by the Cancéropôle PACA (PROCAN project ModeMeta), Project PhysCancer 2009-2013 (grant no. 101194), GEFLUC Marseille-Provence, and the French Association pour la Recherche contre le Cancer (ARC, grant no. 5009).
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
