Quantitative mathematical modeling of clinical brain metastasis dynamics in non-small cell lung cancer.

Brain metastases (BMs) are associated with poor prognosis in non-small cell lung cancer (NSCLC), but are only visible when large enough. Therapeutic decisions such as whole brain radiation therapy would benefit from patient-specific predictions of radiologically undetectable BMs. Here, we propose a mathematical modeling approach and use it to analyze clinical data of BM from NSCLC. Primary tumor growth was best described by a gompertzian model for the pre-diagnosis history, followed by a tumor growth inhibition model during treatment. Growth parameters were estimated only from the size at diagnosis and histology, but predicted plausible individual estimates of the tumor age (2.1-5.3 years). Multiple metastatic models were further assessed from fitting either literature data of BM probability (n = 183 patients) or longitudinal measurements of visible BMs in two patients. Among the tested models, the one featuring dormancy was best able to describe the data. It predicted latency phases of 4.4-5.7 months and onset of BMs 14-19 months before diagnosis. This quantitative model paves the way for a computational tool of potential help during therapeutic management.

Introduction

Lung cancer is the leading cause of cancer-related death worldwide[1]. Nearly 80% of lung cancers are non-small cell lung cancer (NSCLC) and 60% of them are diagnosed at the metastatic stage[2]. Brain metastases (BMs) affect more than 20% of patients with NSCLC[3,4]. Despite recent advances in this field, BMs remain a major challenge as they are associated with poor prognosis[5-8]. In addition, BMs are responsible for disabling symptoms, with impaired patients’ reported outcomes.

Several biological and clinical open questions can be formulated: 1) What was the pre-diagnosis course of a patient presenting with NSCLC? 2) For a patient developing overt BMs, when did cerebral invasion first occur? 3) Were most of the BMs spread by the primary tumor (PT) or did BMs themselves represent a significant source of secondary spread[9-11]? 4) How did growth kinetics compare between the PT and the BMs? Were there dormancy periods during the BMs history[12,13]? 5) For a patient with no or few (≤5) BMs at a given time point, what is the risk and extent of the occult disseminated disease in the brain (if any)?

Answers to these questions would have direct clinical utility. The clinical follow-up and planning of cerebral magnetic resonance imaging (MRI) could indeed highly benefit from individualized predictions of the probability of BM relapse or development. Moreover, the prediction of the risk and timing of BM may help clinical decision between localized (e.g. stereotactic radiotherapy), pan-cerebral (whole-brain radiotherapy (WBRT)) and/or systemic treatment. As of today, prophylactic cranial irradiation in the management of NSCLC is not recommended despite a reduction in the risk of BMs, due to the lack of overall survival benefit and important neuro-cognitive toxicities[6,14,15]. Despite recent improvement in radiotherapy techniques, with the development of hypocampal-sparing WBRT[16], whole brain volume reduction and neurocognitive function decline were described[17]. This highlights the need for a better prediction of patients with highest risk of BM development or relapse, in order to select patients with the higher benefit/risk ratio for WBRT. Several phase III trials were conducted but no firm conclusion applicable for the entire patient population could be drawn[18,19]. The prediction of the risk of BM development and indication for WBRT to prevent brain dissemination is particularly interesting for specific populations such as EGFR or ALK positive NSCLC patients, for whose incidence of BMs is high[20-22]. This points to the need of rational tools to decide therapeutic action in a patient-specific way.

For such aims, quantitative mathematical modeling may be of considerable help, by inferring underlying processes from observed data or providing useful numerical tools in the era of personalized medicine[23,24]. However, despite numerous studies in mathematical and computational oncology, the majority of the efforts have to date focused on mathematical models at the scale of a single tumor[25]. Moreover, the comparison of the models to empirical data remains infrequent. Historically, modeling efforts in the field of metastasis were first initiated by statistical models phenomenologically describing relapse hazards[26-29]. In the 1970’s, Liotta et al. pioneered the development of biologically-based and experimentally-validated models for all the main steps of the metastatic process[30]. Since then, only a small number of studies addressed this topic[31-40]. Of specific relevance to the current work, the Iwata model[31] introduced a size-structured population approach to capture the time development of a population of regional metastases from a hepatocellular carcinoma. It was further extended and studied from the mathematical and numerical viewpoints[41-43], and was found able to reproduce experimental in vivo data[37-39].

To understand biological dynamics of metastasis, animal studies allowing to capture the natural course of the disease[37,38] and the effect of therapeutic interventions in controlled environments[39] provide valuable data for quantitative analysis, although with limited clinical relevance. Experimental procedures are tedious and often only provide access to the total metastatic burden, neglecting its distribution into distinct metastatic lesions when no anatomical imaging modalities are employed[38,39].

Individual-level clinical data of metastatic development are also challenging to obtain because precise number and size of existing lesions are not routinely collected. Apart from the landmark work of Iwata et al.[31] (which only included one untreated patient with hepatocellular carcinoma), no study has yet provided a quantitative modeling analysis of longitudinal data of individual-level number and size of metastases. The purpose of the current work is to develop this approach for BMs dynamics in NSCLC, in order to address the biological and clinical questions stated above.

Results

Pre-diagnosis natural course of lung primary tumors

We investigated two possible growth models for the natural course of lung PTs: exponential and Gompertz. Exponential growth is the simplest model expressing uncontrolled proliferation and is often well adapted to describe tumor growth kinetics during limited observation periods[44]. However, it has been shown that on longer timeframes (typically under volumes increases of 100 to 1000-fold), the specific growth rate of tumors decreases[45], which is well captured by the Gompertz model[45-47]. For calibration of the parameters, we used the data of the PT size and mean doubling time at diagnosis (see Methods). The latter was assumed to be histology-dependent and retrieved from a meta-analysis comprising a total of 388 adenocarcinomas and 377 squamous cell carcinomas[44,48] (see Supplementary Table S1). For the Gompertz model (two parameters), we additionally assumed a carrying capacity of 1012 cells[49]. Under such assumptions, we obtained predictions of pre-diagnosis phases of 19 years versus 5.4 years in the exponential and Gompertz cases, respectively, for an adenocarcinoma with diameter 35 mm (median value from ref.[50]), see Fig. S1. The first prediction seems highly unrealistic compared with previous reports estimating the age of lung primary tumors to be 3–4 years old[51], based on a different method using time to recurrence[52]. The Gompertz prediction on the other hand is rather consistent with the literature range. Moreover, the obtained value of the cellular proliferation rate was realistic, corresponding to a length of the cell cycle of 24.4 hours[53].

We therefore concluded that the Gompertz model was best adapted to describe the pre-diagnosis natural history of the PT.

Probability of brain metastasis as a function of primary tumor size can be described by a mechanistic model of metastasis

PT size was reported by Mujoomdar et al.[50] as a major predictive factor of BM in NSCLC. We investigated whether our mechanistic model of BM was able to fit the data from these authors, at the population level (see Methods for the definition of the model). To do so, we assumed a statistical distribution associated with the inter-individual variability for the metastatic dissemination parameter μ (see Methods). The growth parameters were determined as above, i.e. from histology only. This approach yielded size-dependent population probabilities of BM occurrence, which could be compared to the data. The model fits were reasonably accurate, for both histological types (Fig. 1). Importantly, this demonstrates the descriptive power of our model, which is based on mechanistic simulations of the underlying biological process. In this regard, this is superior to a mere fit of a biologically agnostic statistical model such as a linear regression. Notably, this result was obtained with a minimal number of parameters to describe the inter-patient heterogeneity in metastatic potential. Namely, only two parameters defining the lognormal statistical distributions of μ in the populations. The coefficients of variation were very large (>6,000%). The median value μ gives a quantitative way to measure median BM aggressiveness. It was found two log10 orders of magnitude larger for adenocarcinomas than for squamous cell carcinomas ( vs ).

Figure 1

Probability of brain metastasis as a function of the primary tumor size. Fit of the mechanistic model for probability of BM to data from Mujoomdar et al.[50]. (A) Adenocarcinoma data (n = 136). Inferred values for the distribution of μ: and . (B) Squamous cell carcinoma data (n = 47). Inferred values for the distribution of μ: and .

Model-based comparison of biological scenarios with individual longitudinal data suggests dormancy

Data

These consisted of 10 and 11 PT size measurements, and 47 and 16 BM size measurements in two patients, respectively (one with an adenocarcinoma, the other with a squamous cell carcinoma). We first analyzed the adenocarcinoma patient (Fig. 2) to develop the model and used the second patient to investigate the generalizability of our findings. In the adenocarcinoma patient, the PT first responded to systemic therapy before slowly regrowing (Fig. 2A). A first distant BM was detected 20 months after diagnosis, and kept growing uncontrolled afterwards (Fig. 2B). Other BMs appeared during follow-up (Fig. 2C), reaching a total of 20 BMs at 48 months, date of last examination (Fig. 2D).

Figure 2

Data of primary tumor and brain metastases in a patient with non-small cell lung cancer. (A) Post-diagnosis kinetics of the primary tumor largest diameter, measured on follow-up computerized tomography images (inlet). This EGFR mutated patient was treated first with a tyrosine kinase inhibitor (erlotinib) then with several rounds of additional systemic treatments upon relapse (cytotoxic chemotherapy (docetaxel), re-challenge with erlotinib and a second cytotoxic agent (gemcitabine)), as indicated by arrows and dashed lines in the figure. The solid line corresponds to the model fit during treatment (tumor growth inhibition (TGI) model). (B) Growth kinetics of the brain metastases. The solid line corresponds to Gompertz growth predictions based on parameters estimated from the primary tumor size at diagnosis and histological type, using as initial condition the first observation of each BM. (C) Time course of the apparition of visible metastases. (D) Cerebral magnetic resonance image showing two brain metastases 48 months post-diagnosis (other brain lesions not visible on this slice). EGFR = Epithelial Growth Factor Receptor, TKI = Tyrosine Kinase Inhibitor, mTKI = maintenance Tyrosine Kinase Inhibitor, CT = (cytotoxic) chemotherapy. Time is in months from diagnosis.

Growth kinetics

To model the effect of systemic treatment on the PT, we found that a tumor growth inhibition model[54] (equation (1)) was able to adequately fit the data (Fig. 2A).

Interestingly, for the BM kinetics, we found that the untreated PT Gompertz model predicted surprisingly well the observed data (Figs 2B and S2). Specifically, when employing the Gompertz parameters used for the pre-diagnosis PT history, and using initial conditions on each BM, the future growth curves were accurately predicted. This indicated that: 1) for this patient, the BMs did not respond to systemic therapy, at least from the time of first observation, 2) it is reasonable to assume that all BMs grew at the same growth rate, 3) the BMs growth rate might be similar to the PT growth rate, at least during the clinically overt phase.

Quantitative assessment of five biologically-based models of metastatic dissemination and colonization

However, this mere description of the BMs growth is not satisfactory as a model of systemic disease, since the dissemination part is absent. In particular, no model is given for the initiation times of the BMs. To include the dissemination component of the metastatic process, we relied on a modeling framework first initiated by Iwata et al.[31] (see Methods).

Five models, corresponding to five possible biological scenarios, were investigated (Fig. 3). These included: (A) a “basic” model with same growth model and parameters for the PT and the BMs and no secondary dissemination (i.e. BMs spread only by the PT); (B) a model allowing different growth parameters for the PT and the BMs; (C) a model with secondary dissemination, i.e. the ability of BMs to spread BMs themselves[31]; (D) a model allowing a delay before initiation of metastatic ability of the PT (the so-called linear model in which dissemination occurs at a late stage, opposed to the parallel model where dissemination is an early event[49]) and (E) a model allowing BM dormancy phases, i.e. time periods with null growth rate[12,13,55]. Each of these scenario was translated into corresponding mathematical equations (see Methods).

Figure 3

Schematic of the models and investigated hypotheses. Several biological hypotheses about metastatic dynamics were formalized as mathematical models and tested against the data. These included: same (A) or different (B) growth rates between primary and secondary tumors; only primary dissemination from the primary tumor or (C) secondary dissemination from brain metastases themselves; (D) possibility of a delay between cancer initiation and onset of metastatic ability; and (E) the possibility of a dormancy lag time between metastatic spread and growth initiation.

The models exhibited different descriptive power, as quantified by the best-fit value of the objective function (Table 1). The best description of the data was obtained by the dormancy model (Fig. 4). Importantly, this model achieved accurate goodness-of-fit for 47 BM size measurements over 6 time points (Fig. 4A,B) with only 3 parameters fitted from the BM sizes data (μ, and γ). Consequently, parameter identifiability was found excellent (relative standard errors <5%, Table 2). Apparition of the first BM above the visibility threshold was also accurately described (19 vs 20 months, model vs data, Fig. 4C). The model overestimated the number of small visible BMs at 40 months, possibly because a lot of these were present but too small to be visible at MRI. Indeed, a lot of BMs were present at 43 months – in accordance with the model (Fig. 4A) – which thus emphasizes its predictive potential.

Table 1

Minimal value of the objective function obtained when fitting each of the models to the data.

Model	Patient 1	Patient 2
Basic	5.51	2.53
Secondary	5.43	2.3
Delay	5.23	1.53
Dormancy	4.93	1.71
Diff. growth	4.95	1.79

The objective function was the sum of squared residuals (see expression (7)).

Figure 4

Fit of the dormancy model. (A) Time course of the visible brain metastases (BMs) size distributions during follow-up. Comparison between model calibration and data. (B) Time course of the number of visible BMs. T = time of diagnosis. T = detection time of the first brain metastasis. (C) Comparison of the BM size distribution between the model fit and the data at last follow-up (48 months post-diagnosis). Time is in months from diagnosis.

Table 2

Patient-specific parameters.

Parameter	Meaning	Unit	Patient 1	Patient 2	Determination
S d	PT size at diagnosis	mm	36.0	53.7	Data
Histology			Adenocarcinoma	Squamous cell carcinoma	Data
T d	Age of the PT at diagnosis	year	5.3	2.1	Predicted
Visibility threshold		mm	3	3	Fixed
S 0, p	Initial size of the PT	cell	1	1	Fixed
S 0	Initial size of a BM	cell	1	1	Fixed
αp	Proliferation rate at one cell (PT)	day−1	0.0284	0.0858	Computed
βp	Exponential decrease of growth rate (PT)	day−1	1.03 × 10−3	3.10 × 10−3	Computed
DT PT	PT doubling time at 3 cm	day	173	57	Predicted
α1	PT growth rate during relapse	day−1	5.72 × 10−4 (23.6)	0.0101 (5.28 × 10−3)	Fit PT
κ	PT log-kill effect	day−1	4.46 × 10−3 (20.6)	0.0439 (1803)	Fit PT
t res	Half-life of treatment response	day	149 (1.05 × 10−3)	35.4 (0.998)	Fit PT
α0	Proliferation rate at one cell (BM)	day−1	0.0284	0.0858	Computed
β	Exponential decrease of growth rate (BM)	day−1	1.03 × 10−3	3.10 × 10−3	Computed
τ	Dormancy duration	day	133 (4.23)	171 (33.3)	Fit BM
DT BM	BM doubling time at 3 cm	day	173	57	Predicted
μ	Cellular metastatic potential	cell−1 × day−1	2.00 × 10−12 (2.83)	1.02 × 10−12 (71.9)	Fit BM
γ	Fractal (Hausdorff) dimension of the vasculature (x1/3)	—	1	1	Fit BM*
d	Dissemination rate at 3 cm	day−1	0.0282	0.0144	Predicted

Values of the parameters from either the data (Data), the model assumptions (Fixed), estimated from direct computations from diagnosis data (Computed), fit of the TGI model to the PT data (Fit PT) or fit of the metastatic dormancy model to the brain metastases data (Fit BM), or computed from these estimations (Predicted). Values in parenthesis for fitted parameters correspond to the relative standard errors expressed in percent, i.e. where se is the standard error (square root of the covariance diagonal entry, see Materials an methods) and est is the parameter estimate. *For value of the parameter , due to identifiability issues, only five possible values were considered during the fit procedure (see Methods).

The basic (A), secondary (C) and delay models (D) exhibited substantially worse goodness-of-fit (Table 1 and Figs S3–5), suggesting that their underlying hypotheses are insufficient to provide valid explanations of BM dynamics. In the delay model, the fit was improved compared to models (A) and (B) but the value of the delay t was found unrealistically large (4.8 years after the first cancer cell, corresponding to 6 months before diagnosis).

Finally, model (D) yielded a significant improvement of the goodness-of-fit for both the dynamics of the number of visible BMs and cumulative size distributions, while not deteriorating the practical identifiability of the parameters (Fig. S6). Of note, the estimated values of the BM growth parameters were close to the PT ones. However, given the previous observation that during the clinical phase the BMs grew at a growth rate consistent with the one predicted for the PT, combined with the facts that the optimal objective value was achieved for the dormancy model and that it is more parsimonious (one parameter less), we selected the dormancy model for further analysis.

Clinically relevant simulations of the disease course reveal times of BM initiation

From the quantitative calibration of the dormancy model to the data, several predictions of clinical interest can be made. The value of γ that generated the best-fit was 1, suggesting equal probability for all the cells within the PT to disseminate. In turn, this could indicate of a well-vascularized tumor, which might be prone to anti-angiogenic therapy. Interestingly, the value of μ inferred from this patient-level data was in the same range as the one inferred from the above population analysis (μ = 2 × 10−12 cell−1.day−1 versus μ = 1.39 × 10−11 cell−1.day−1), giving further support to our approach.

Once calibrated from the data, our model allowed to simulate the predicted natural history of the disease. The Supplementary Movie S1 shows a simulation of the PT growth together with the apparition of the entire population of BMs (visible + invisible). Stars indicate tumors that are invisible (<5 mm), and the BMs size distribution time course is also simulated. BMs represented in gray were born before the diagnosis, while BMs in white are the ones born after.

Apart from the age of the PT, prediction of birth times and invisible BM burden at any time point can be performed. Interestingly, we found that the first BM – which clinically appeared 19 months after diagnosis – was already present 14 months before clinical detection of the PT (Fig. 5A). In fact, at the time of diagnosis, our model predicted the presence of 10 occult BMs, representing a total burden of 1,167 cells mostly distributed into the largest (first) BM (size 1,046 cells  0.126 mm), see Fig. 5B. Notably, at the time the PT had reached the imaging detection limit, no metastasis had occurred yet. This suggests that if the disease had been detected earlier (e.g., by systematic screening), BM occurrence might have been prevented provided the tumor had been operable at this time. The amount of predicted BMs present at T was 28 tumors, with an overall BM burden of 1.7 × 107 cells (Fig. 5C). Therefore, provided that neuro-cognitive risks would be acceptable, the model would recommend WBRT or systemic treatment rather than localized intervention only.

Figure 5

Model simulated predictions of the pre- and post-diagnosis time course of the primary and cerebral disease. (A) Model-inferred growth kinetics of the primary tumor (in blue) and the brain metastases (in red), compared with data (circles). T = time of first cancer cell, T = time of diagnosis, T = detection time of the first brain metastasis. Only brain metastases that will become visible are shown. (B) Predicted size distribution of the brain metastases at diagnosis. (C) Predicted size distribution of the brain metastases at the time of clinical occurrence of the first one. Simulations performed with the discrete version of the model. PT = primary tumor, BM = brain metastasis. Time is in months from diagnosis.

Together, these results demonstrate the potential clinical utility of the model for prediction of the invisible BM, in order to inform therapeutic decision.

Generalization to a second patient

To test whether the dormancy model was generalizable, we used data from a second patient, which was not employed during the model development phase. Given the different histological type of the lung PT (squamous cell carcinoma), we adapted the doubling time accordingly (Table S1) and found a younger age of the PT of 2.1 years. Estimation of the log-kill parameter from the TGI PT model during therapy suffered from lack of identifiability, due to an estimated short duration of treatment effectiveness (see parameter in Table 2). Response of the PT to therapy was indeed characterized by a faster relapse (parameter α1), compared with the first patient. The qualitative structure of the model confirmed its descriptive power by being able to give an accurate description the BMs size dynamics (Figs S7 and 8), while the “basic” and “secondary dissemination” models were not (data not shown). Several parameters appeared to be patient-dependent, such as the dormancy duration , estimated to 171 days versus 133 days for the first patient. Due to the lower number of data points available for this patient, parameters identifiability was found worse, but still acceptable (Table 2).

Resulting clinical predictions were distinct (Fig. S8), emphasizing the patient-specific nature of BM dynamics. The first BM was clinically detected 14 months after diagnosis, but was predicted to have been disseminated 45 months prior to diagnosis. While for both patients cerebral dissemination had already occurred at the time of diagnosis, its extent was different with a much larger mass (>1,000 cells) for the first patient than for the second (8 cells). This is due to the long period of dormancy for patient 2, resulting in all 8 predicted BMs being still dormant at the time of diagnosis.

These results demonstrate the flexibility of our model, which is able to describe and predict distinct situations in terms of repartition of the metastatic burden into individual lesions.

Discussion

Using both population- and individual-level data of BM development in NSCLC, we have developed a general method based on biologically grounded computational models that allowed: 1) to infer the disease age from data on PT size and histology only, 2) to test multiple scenarios of metastatic dissemination and colonization against macroscopic data available in the clinic and 3) to infer times of BM initiations and number and size of invisible lesions during the clinical course of the disease.

Estimation of the duration of the pre-diagnosis PT growth phase has important implications, since BMs are more likely to have occurred for an old PT compared to a young one. Our results showed a significant difference whether considering exponential or Gompertz growth, which was consistent with previous findings[48]. Of note, our age estimates of 5.3 and 2.1 years-old are in relative agreement with the 3–4 years range found by others relying on a different method[51].

Importantly, we were then able to include this description of PT past growth into a mechanistic model to describe the probability of BM occurrence. We found that the cell-scale dissemination parameter μ was significantly larger in adenocarcinomas than in squamous cell carcinomas, which offers a quantitative theory for the reported differences in BM occurrence between these two histological types[50].

More generally, our computational platform provided a way to translate biological findings into clinically useful numerical tools. However, in order to provide robust inference, the complexity of the models had to remain commensurate with the available data. Thus, several higher order phenomena relevant to the metastatic process were ignored or aggregated into mesoscopic parameters. For instance, the metastatic potential μ is the product of several cell-scale probabilities relating to the multiple steps of the metastatic cascade[56]. The median values inferred from population analyses based on probabilities of BM ( and ) matched relatively well the ones that were found from analysis of patient-specific data ( and ), given the standard deviations. The larger value found in patient 1 could be due to histology or to the EGFR mutational status, known to impact on the BM aggressiveness[3,22,57,58]. Interestingly, while relying on distinct modeling techniques, the value of μ that we obtained was in the same range as obtained by others using stochastic evolutionary modeling[34,59].

When evaluating the descriptive power of multiple biologically-based models against longitudinal data of BM sizes, our results suggested possible phases of dormancy. This finding is supported by several preclinical and clinical reports[13,60]. For instance, Nguyen et al. mention possible dormancy phases for cerebral metastases from lung cancer[13]. Durations of these latent periods are mentioned to be of the order of months, which is consistent with our quantitative results (4.4 and 5.7 months for patients 1 and 2, respectively). BM dormancy could possibly be explained by the fact that colonizing lung carcinoma cells face a foreign micro-environment when reaching the brain, which hampers their establishment and growth[61]. This could cause single cells to remain quiescent or tumors to remain at a small, pre-angiogenic size[12,62]. Indeed, vascular endothelial growth factor-A (VEGF-A) has been shown to mediate BM dormancy in a mouse model of lung cancer[60].

Several second-order phenomena were ignored here for the sake of identifiability, but could nevertheless impact BM dynamics and dormancy. These include tumor-tumor interactions, either through soluble circulating factors[63] distantly inhibiting growth and possibly maintaining dormancy[62], or by exchange of tumor cells between established lesions[9,64]. We have recently proposed a model for distant interactions that was validated in a two-tumors experimental system and could be incorporated into the current modeling platform[65,66].

While not problematic for the two patients that were investigated here since BMs did not respond to the systemic treatment – possibly because of the blood-brain barrier hampering delivery of the anti-cancer agents – a major limitation of our model is that the effect of therapy on the BMs was neglected. We intend to address this in future work, in particular for optimizing and personalizing combination therapies[24,67,68]. Moreover, the current study needs to be extended to a larger number of patients.

Our study demonstrates the potential clinical utility of our computational platform as a personalized predictive tool of BM that could be helpful in the decision of WBRT versus localized treatment only. For instance for patient 1, at the time of detection of the first BM or at 40 months when only 6 BMs were visible, large amounts of invisible BMs were predicted by the model, which would advise in favor of WBRT. In contrast, a smaller metastatic burden was predicted in patient 2, thus advising for localized intervention only. Thus, our quantitative modeling approach could complement existing predictive metrics such as the brain metastasis velocity index[69,70], for patients undergoing stereotactic radiosurgery (SRS), i.e. patients presenting with a small (≤3–5) number of BMs. While the latter requires data at a second time point after a first SRS, our model could predict the invisible burden of BM already at the first SRS, and consequently inform the risk and timing of future relapse.

In order to translate our findings into a clinically usable tool, further methods need to be developed to calibrate the parameters of the model from data already available at diagnosis or at the first BM occurrence. By relying on a small number (3) of free parameters for the BM dynamics, our model represents a realistic candidate for such a purpose. To this respect, in addition to routine clinico-demographic features, molecular gene expression signatures[71] (for parameter μ, for instance) as well as radiomics features[72] (for parameter γ) might represent valuable resources.

Materials and Methods

Patient data

The data used in this study concerned patients with non-small cell lung cancer (NSCLC) and were of two distinct natures: 1) population data of probability of BM as a function of PT size retrieved from the literature[50] and 2) longitudinal measurements of PT and BM diameters in two patients with NSCLC retrieved from imaging data (CT scans for lung lesions, MRI for brain tumors). Both patients had unresectable PT at diagnosis. The first patient (used for model development) was extracted from an EGFR mutated cohort from Institut Bergonié (Bordeaux, France). The second patient (used for model validation) had an EGFR wild type squamous cell lung carcinoma and was issued from routinely treated patients in the thoracic oncology service of the University Hospital of Marseille. The data comprised 10 PT sizes and 47 BM sizes (spanning 6 time points) for the first patient and 11 PT sizes and 16 BM sizes (spanning 4 time points) for second patient.

Data use from the Marseille patient was approved by a national ethics committee (CEPRO, Comité d’Evaluation des Protocoles de Recherche Observationnelle, reference number 2015-041), according to French law. Data use from the Marseille patient was approved by the institutional ethics committee of Institut Bergonié (Comité de Protection des Personnes Sud-Ouest et Outre Mer III) and IRB approval was obtained for use of the images. The need for written informed consent of data collection was waived for this patient, in accordance with the related policy of Institut Bergonié. The study was performed in accordance with the Declaration of Helsinki, Good Clinical Practices, and local ethical and legal requirements.

Mathematical modeling of primary tumor growth and metastatic development

Primary tumor growth

The pre-diagnosis natural history of the primary tumor size S(t) for times (diagnosis time) – expressed in number of cells in the model computations – was assumed to follow the Gompertz growth model[46,47], i.e.where time t = 0 corresponds to the first cancer cell, parameter α0, is the specific growth rate (i.e. ) at this time and β is the exponential rate of decrease of the specific growth rate. Conversions from diameter measurements to number of cells were performed assuming spherical shape and the classical assumption 1 mm3 = 106 cells[73]. After treatment start (T), the primary tumor size was assumed to follow a tumor growth inhibition model[54] consisting of: 1) exponential growth (rate α1), 2) log-kill effect of the therapy (efficacy parameter )[74] and 3) exponential decrease of the treatment effect due to resistance, with half-life t. The equation is:

Metastatic development

Basic and different growth models. The general modeling framework we employed was derived from the work of Iwata et al.[31]. It consists in modeling the population of metastases by means of a size-structured density , of use to distinguish between visible and invisible tumors. Metastatic development of the disease is reduced to two main phases: dissemination and colonization[75]. The multiple steps of the metastatic cascade[56] are aggregated into a dissemination rate with expression:

which corresponds to the number of successfully born BM per unit of time. In this expression, the geometric parameter γ corresponds to the intra tumor repartition of the cells susceptible to yield metastasis and μ is the per day per cell probability for such cells to overcome all the steps of the metastatic cascade (acquisition of metastatic-specific mutations, epithelial-to-mesenchymal transition, invasion of the surrounding tissue, intravasation, survival in transit, extravasation and survival in the brain). For γ = 1 all cells in the PT have equal probability to give a BM whereas a value of γ = 0 indicates a constant pool of cells having metastatic ability (cancer stem cells). Intermediate values 0 < γ < 1 can be interpreted as the geometric disposition of the metastatic-able cells, including the surface of the tumor (γ = 2/3) or a fractional dimension linked to the fractal nature of the tumor vasculature[76]. Assuming further that the growth of the metastasis follows a gompertzian growth ratewith growth parameters α0 and β possibly equal (basic model (A)) or distinct (different growth model (B)) compared to the PT ones, the density satisfies the following transport partial differential equation[31]:where S0 is the size of a BM at birth (here assumed to be one cell). From the solution to this equation, the main quantity of interest for comparison to the empirical data is the number of metastasis larger than a given size s (cumulative size distribution):

The total number of metastases – denoted N(t) – is obtained by using s = V0 above and its expression can be directly computed without solving the entire problem (3) as it is given by:

Using the method of characteristics, one can derive the following relationship between N and f:where t(s) is the time for a tumor growing at rate g to reach the size s. In the case of Gompertz growth one has:

Of particular interest is the number of visible BMs with S the minimal visible size at CT scan taken here to be 5 mm in diameter.

Secondary dissemination. In the previous model formulation, all BMs were assumed to have been seeded by the primary tumor. When BMs are also allowed to spread metastases themselves, this results in a second term in the boundary condition of (3) and the model becomes[31]:

In this case, formula (4) is not valid anymore, which complicates substantially the computation of the cumulative size distribution. A dedicated scheme based on the method of characteristics was employed[43].

Delay model. Consideration of a delay t0 before onset of metastatic dissemination in the model can be taken into account by remarking that

Thus

Dormancy model. For inclusion of dormancy in the model – defined as a period of duration during which a newborn metastasis remains at size S0 – the time to reach any given size becomes . The cumulative size distribution is then given by:

Discrete versions of the models. While continuous versions of the models were used for fitting the model to the data because they allow computations to be tractable, discrete versions were implemented for forward calculations, because of the small number of BMs. Briefly, in a discrete simulation, the appearance time of the i-th BM T was defined bywith N(t) the number of metastases as defined above. The size of the i-th BM S(t) was then defined, for t > T, by:

This discrete version corresponds to a Poisson process where the number of events (metastatic births) is replaced by its expectation: . For further links between stochastic (Poisson process) and continuous versions of the Iwata model, the reader is referred to[77]. Simulations reported in Fig. 5B,C and 6 were generated using discrete versions of the models.

Models fitting and estimation of the parameters

Pre-clinical primary tumor growth

To parameterize the Gompertz function defining the PT growth, two parameters needed to be defined (α0, and β). In the absence of longitudinal measurements of the PT size without treatment, these two parameters were determined from two considerations: 1) the maximal reachable size (carrying capacity, equal to ) of a human tumor is 1012 cells[29,49] and 2) the histology-dependent value of the doubling time at diagnosis, retrieved from a meta-analysis of published literature about the natural growth of lung PTs (see Table S1, extended from[44,48]). The latter yielded values of 201 days for an adenocarcinoma and 104 days for an undifferentiated carcinoma. For the Gompertz model, the doubling time is size-dependent and its value for the PT diagnosis size DT(S) is given by:

Using the formula linking α0, and β to the carrying capacity, this nonlinear equation was numerically solved.

Population level: probability of BM apparition

To fit the data from[50] describing the probability of BM in a population of lung adenocarcinoma patients, we employed a previously described methodology[39,78]. Briefly, we considered that the probability of developing BM after diagnosis was the probability of having already BM at diagnosis, i.e. , with N defined by equation (4). Of note, since this quantity only depends on the dissemination rate, this probability is equal for all the model scenarios introduced above. We fixed the value of the PT growth parameters as described above from the cohort histology and set γ = 1 as the simplest dissemination model. The inter-individual variability was then minimally modeled as resulting from a lognormal population distribution of the parameter Uniform distributions of the PT diameters were assuming within each interval of the PT sizes given as data. The probability of developing a metastasis with a PT size writes:

The best-fit of these probabilities – evaluated by Monte Carlo simulations – to the empirical data was then determined by least squares minimization performed using the function fminsearch of Matlab (Nelder-Mead algorithm)[79].

Individual level: description of longitudinal data of number and size of BM growth

Maximum likelihood estimation. Due to the discrete nature of the data at the individual level (diameters of a small number of BMs at discrete time points), a direct comparison between the size distribution solution of the problem (3) was not possible. Instead, we compared the data to the model by means of the cumulative size distribution. Denoting by t the observation times, the sorted BM sizes at time t and the number of metastases larger than size at time t, we considered the following nonlinear regression problem:where regroups all the parameters of the model. Note that, except for the “secondary dissemination” model, all models can be viewed as submodels of a general model including all these parameters (the “basic model” consisting of the case and , for instance). Classical maximum likelihood estimation then leads to the following estimate:

Parameters identifiability. Standard errors can be computed from this statistics’ covariance matrix, given by[80]:where is the jacobian matrix of the model (with respect to the parameter vector θ) at all time points and all sizes, evaluated at the optimal parameter and is the a posteriori estimate of σ given bywith N the total number of data points and P the number of free parameters.

Using the standard errors as an identifiability metric, we repeatedly observed a lack of identifiability of parameters μ and γ when fitted together. Indeed, standard errors for μ and γ were larger than 200% when fitting the basic model to the data. Further investigation of the shape of the objective function confirmed this lack of idenfiability (Fig. S7). To address this issue, we only considered a finite set of relevant possible values for γ and only optimized the value of μ. These values were and corresponding initial conditions for μ were . When more parameters were let free in the model (delay t or dormancy period ), we generated 4 × 4 parameters grids for initial conditions (with or ). Of these 16 optimization problems, the one with the minimal value of l at convergence was selected.

Supplementary Material

Movie S1

Movie S2