# Metamats: A mechanistic software for the simulation, inference and prediction of clinical metastasis

Célestin Bigarré, Alice Daumas, Laurent Greillier, Xavier Muracciole, Laeticia Padovani, Sébastien Benzekry

► To cite this version:

Célestin Bigarré, Alice Daumas, Laurent Greillier, Xavier Muracciole, Laeticia Padovani, et al.. *Metamats: A mechanistic software for the simulation, inference and prediction of clinical metastasis*. 2025. *hal-04935804v2*

HAL Id: hal-04935804

https://inria.hal.science/hal-04935804v2

Preprint submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# METAMATS: A mechanistic software for the simulation, inference and prediction of clinical metastasis

Célestin Bigarré¹✉, Alice Daumas¹, ², Laurent Greiller¹, ³, Xavier Muracciole¹, ², Laeticia Padovani², ⁴, Sébastien Benzekry¹

¹COMPO, Inria Méditerranée, Cancer Research Center of Marseille, Inserm UMR1068, CNRS UMR7258, UM105, Aix-Marseille Université, 13273 Marseille, France.; ²Oncology Radiotherapy Department, Assistance Publique des Hôpitaux de Marseille, Aix-Marseille Université, 13005 Marseille, France.; ³Multidisciplinary Oncology and Therapeutic Innovations Department, Assistance Publique des Hôpitaux de Marseille, Aix-Marseille Université, 13005 Marseille, France.; ⁴Scientific Research National Center (CNRS), Institute of Neurophysiopathology, Aix-Marseille University, 13005 Marseille, France.

## Abstract

The development of metastases is a complex process that can be better understood and predicted using mathematical mechanistic models. We propose "METAMATS", a ready-to-use modeling framework implementing a semi-mechanistic model for simulating the time dynamics of metastatic development, including a primary tumor and a population of metastatic tumors. It relies on a reduced set of mathematical parameters: $\alpha$ (tumor proliferation rate), $\mu$ (dissemination rate), and $\gamma$ (fractal scale parameter). "METAMATS" supports both individual and population-level analyses.

At the individual level, it can simulate metastatic dynamics, infer parameters from longitudinal metastasis size data, and predict the natural history of cancer both retrospectively and prospectively. At the population level, leveraging nonlinear mixed effects modeling, "METAMATS" performs inference of parameter distributions and assessment of biologically interpretable covariate effects on proliferation and/or dissemination, from distant metastasis-free survival data. "METAMATS" also serves as a generative model for simulating virtual patients or populations.

We demonstrate its applicability for modeling, inference, prediction, and simulation in a clinical setting: the dynamics of brain metastases (BM) in small-cell lung cancer patients and the impact of prophylactic cranial irradiation (PCI). Data included two cohorts: 103 patients from Assistance publique-Hôpitaux de Marseille with longitudinal individual measurements of BM sizes, and 100 patients from the CONVERT study (NCT00433563). PCI was found to have an impact on BM by significantly reducing metastatic appearance (parameter $\mu$), rather than metastatic growth (parameter $\alpha$), a biological finding impossible to obtain by means of classical survival analysis only.
Introduction

The development of metastases is the leading cause of death in cancer patients (Dillekås et al., 2019) and most patients with metastatic cancer die from the consequences of the cancerous disease (Mani et al., 2024). Some metastatic sites are of further clinical importance, such as brain metastasis (BM) for lung cancer, where 20% of the patients present BM at the time of diagnosis (Achrol et al., 2019; Barnholtz-Sloan et al., 2004; Nayak et al., 2012), often leading to disability symptoms. Treatment of metastases is a major challenge and is based on a combination of systemic and local treatments (Rittberg et al., 2021). Understanding the onset and development of metastatic disease is thus of utmost clinical importance.

Mechanistic mathematical models can help in this endeavor by translating biomedical knowl-
edge into equations of dynamical systems. Subsequently implemented into computational tools,
they have relevance for testing biological hypotheses (Baratchart et al., 2015; Benzekry et al., 2016)
or predicting outcomes, at the individual or population level (Nicolò et al., 2020). However, only
a few mathematical models of metastases have been proposed, and most of them focus on de-
scribing the intrinsic properties of the metastatic process at the cellular level (Bocci et al., 2019;
Gallaher et al., 2023; Haeno et al., 2012; McDonald et al., 2023). One of the main challenges is
comparing the model outputs with clinical and experimental data. This arises from the difficulty
of accessing high-quality clinical data, the large complexity of some models, and the lack of a well-
defined mathematical basis for the model fitting procedure. In fact, most studies rely on simple
least-squares criteria without a dedicated statistical inference framework (Strobl et al., 2021).

A seminal paper by Iwata et al. (Iwata et al., 2000) proposed a model of the population of metastases based on a transport partial differential equation. This model describes the growth and appearance of metastases and was initially applied to a single patient with hepatocellular carcinoma and liver metastases. It has since been extended and used in different settings (Benzekry et al., 2016) (Bilous et al., 2019) (Nicolò et al., 2020) (Bigarré et al., 2023) (Benzekry et al., 2024) (Schlicke et al., 2021).

With the rise of computational methods, it is now widespread, especially in the machine learn-
ing community, to use software packages (Comets et al., 2017; Masotti et al., 2023; Nicol et al.,
2022) to share implementations of new methods. These allow the scientific community to build
on previously published models. Currently, when metastasis modeling studies have released the
code for their analysis, it is often an experimental research code that is difficult to reuse (Benzekry
et al., 2024).

The value of semi-mechanistic models of metastatic progression lies in the ability of researchers to apply such models to a large variety of data sets varying biological assumptions, clinical settings, and mathematical formulation. This can only be achieved if the models are publicly released in a form that allows others to adapt the model to fit their use-case. Evaluation and calibration of computational models of metastatic dissemination necessarily involves comparison with real-world data. Two types of data can be used to this effect: individual longitudinal follow-up data (e.g., the sizes of the primary tumor and metastases), and standard clinical endpoints such as overall and metastases-free survival. For the former, identification of the model parameters requires simulation of the model; inference is classically performed by minimization of a likelihood-like criterion (e.g., least squares). However, model fitting is only rarely embedded in a robust statistical framework. Population data requires development of a novel statistical framework because 1)
the presence of right-censoring and 2) the necessity to define a population-level statistical model. Such a framework has been initially proposed in (Nicolò et al., 2020) and further developed in (Bilous et al., 2019), relying on the mixed-effects approach (Lavielle, 2014). Finally, calibration of models, even when based on numeric criteria or metrics that are not always available, requires visualization tools to ensure the correspondence between the data and the model outputs.

Here, we developed a modeling framework encompassing modeling, inference, visualization, and simulation of metastasis, at the individual and population level. This framework was implemented in a series of R packages called METAMATS that followed high-level standards of software development. As an example of the capabilities of the METAMATS software we studied the impact of prophylactic cranial irradiation (PCI) in the brain-metastasis free survival of patients with small cell lung cancer (SCLC). Using clinical data from the Public hospitals of Marseille as well as as data from the CONVERT study (NCT00433563) we demonstrated both the individual and population level capabilities of the METAMATS Software.

# Materials and methods

## A mechanistic model of metastasis

### Derivation of the model

We considered that cancerous diseases begin with a primary tumor (PT) that starts at an initial size of one cell and will grow according to a specified growth law. As the primary tumor grows, each cell has a small probability of forming a new distant metastasis daily. When seeded, metastases begin independent growth (which may or may not follow the same growth law as for the PT) and may begin to seed new metastases.

The following modeling framework can be used to simulate this situation with significant control over the different parts of the model.

Time is denoted by $t$ and $t_0$ is the time of initiation of the disease. The size of the PT at time $t$ is denoted by $V_p(t)$, and is assumed to follow the differential equation:

$$ \begin{cases} dV_p/dt = g_p(V_p) \\ V_p(t_0) = V_0 \end{cases}, $$

with $V_0$ the initial tumor size. We assumed that $V_p(t) = 0$ for all $t < t_0$.

Similarly, a metastasis of age $\tau$ is of size $V(\tau)$, given by the solution of the differential equation:

$$ \begin{cases} dV/dt = g(V) \\ V(t_0) = V_0 \end{cases}. $$

The daily number of new metastases seeding from the PT (respectively, metastases) as a function of tumor size is denoted by $d_p(V)$ (resp $d(V)$), where $V \ge 0$.

The population of metastases is described with a tumor size density function $\rho$, such that the integral of $\rho(t, V)$ between $V_1$ and $V_2$ represents the number of metastases with a size between $V_1$ and $V_2$ at time $t$.

The general model, given $t_0$, $V_0$, $g_p$, $g$, $d_p$, and $d$, is given by the partial differential equation (PDE):

$$ \partial_t \rho(V,t) + \partial_V(g(V)\rho(t,V)) = 0, $$

with the initial condition,

$$ \rho(V, t_0) = 0, $$
ensuring that there are no metastases initially. The boundary condition

$$g(V_0) \rho(V_0, t) = d_p(V_p(t)) + \int_{V_0}^{\infty} d(V) \rho(V, t) dV,$$

defines the apparition of new metastases, with the first term the metastasis originating from the PT and the second term the ones originating from the metastases themselves (Iwata et al., 2000).

The total number of metastases at time t is given by

$$N(t) = \int_{V_0}^{\infty} \rho(V, t) dV.$$

## Discrete model

In computational settings, it is easier to use a discrete model version. This implementation was defined on a time grid $(t_k)_{k \in [0..1]}$, but only for models without secondary dissemination, that is, $d(v) = 0$ for all v. Instead of solving the PDE, we computed the number of metastases using a midpoint integration method:

$$\tilde{N}(t_k) = \lfloor \sum_{l \le l \le k} (t_l - t_{l-1}) \times \frac{d(V_p(t_{l-1})) + d(V_p(t_l))}{2} \rfloor,$$

where $\lfloor x \rfloor$ represents the integer part of x.

Since all tumors grow independently and $\tilde{N}$ is always an integer, the number of metastases of size greater than $V$ at time $t$ (the discrete cumulative size distribution) is given by

$$c \tilde{s} d(V, t_k) = \tilde{N}(t_k - \tau(V)),$$

where $\tau(V) = \arg \inf_t V(t) \ge V$ is the time it takes for a metastasis to reach size $V$.

Where necessary for clarity, the model is expressed with explicit dependency on the complete parameter vector $\phi$, which contains all free parameters for the growth and dissemination function, as

$$c \tilde{s} d(V, t_k | \phi)$$

## Time to first metastasis

If we set a minimum tumor size for detection at imaging at $V_{threshold}$, the minimum time required to detect a metastasis is given by

$$T_{meta}(\phi) = \min\{t \ge 0 | \tilde{\rho}(V_{threshold}, t | \phi) \ge 1.\}$$

## Metastatic burden

A key element to describe the extent of metastatic disease is the metastatic burden (MB), that is, the total number of metastatic cells present in the patient. It is given by:

$$M(t) = \int_{1}^{+\infty} V \rho(t, V) dV,$$

and for the discrete model:

$$\tilde{M}(t_k) = \sum_k V_m(t_k),$$

where $V_m(t)$ is the size of the m-th metastasis at time $t$.
# Data

## AP-HM dataset

The AP-HM dataset is composed of data from 103 patients treated for small cell lung cancer (SCLC) at the Assistance publique - Hôpitaux de Marseille (AP-HM) institution. All patients in the dataset were diagnosed with non-brain-metastatic SCLC. We used this dataset to demonstrate the individual-level functionalities of the METAMATS package. The complete clinical characteristics of this population are presented in Table 1. In the AP-HM population, 19% of the patients had a Tumor Node Metastasis (TNM) T staging of 1, 20% a T staging of 2, 32% a T staging of 3, and 28% a T staging of 4 Rami-Porta et al., 2015. Fifty-one percent of the patients had disseminated disease at the time of diagnosis. Most of the patients had a good performance status (PS) at the time of diagnosis (44% with PS of 0, 46% with PS of 1), some had degraded PS (9.7% with PS of 2) and one patient (1%) was severely affected and unable to fully take care of himself (PS of 3). Primary treatment consisted of a combination of surgery, radiation therapy, and chemotherapy. In addition, some patients (17%) received adjuvant immunotherapy.

In this population, around 46% of the patients received prophylactic cranial irradiation (PCI) after the initial course of treatment as prescribed by the French national guidelines; the remaining 54% did not receive PCI by choice or due to contraindications. A comparison between the PCI and non-PCI groups failed to show statistical differences for the clinical variables reported Table 1.

Patients were followed with regular brain imaging to assess the presence and size of brain metastases. Follow-up time ranged from 5 months to 8 years, with a median of 2 years. Out of the 103 patients, 41 developed brain metastases and were used to demonstrate the individual-level functionalities of METAMATS. Brain metastases were individually measured at multiple times (the number of time points ranged from 1 to 11, with a median of 4). Along the course of their disease, these patients developed between 1 and 28 brain metastases (median 4 brain metastases).

## BMFS dataset

Data on TNM tumor size Rami-Porta et al., 2015, brain metastasis-free survival, and PCI status for patients from the AP-HM dataset were combined with data from 87 patients who participated in the CONVERT study (Faivre-Finn et al., 2017) (NCT00433563) of which 94% received PCI and 6% did not, for a total of 190 patients of which 52% experienced brain metastases. Altogether with the time to first metastasis from the AP-HM dataset, these data form the Brain Metastasis-Free Survival (BMFS) dataset (see Supplementary Table S1 and Supplementary Figure S2).

To compensate for the structural absence of patients with brain metastases at the time of diagnosis in our dataset, we supplemented the dataset with a proportion of "phantom" patients, left-censored at time 0. We added around 18% of phantom patients to match the ratio for the Surveillance, Epidemiology, and End Results Program data Li et al., 2021, TNM T status for phantom patients was generated randomly using the proportions described in the same article. This second dataset was used to estimate the parameters at the population level of the model using the statistical framework described in the section "Population survival model".

## TNM to tumor size conversion

As the model needs the quantitative size of the PT as an input parameter, we had to convert the TNM tumor sizes (TNM T classification) to a number of cells in the primary tumor. We chose to consider PT as isotropic spheres with a density of $1 \times 10^6 \text{ cell mm}^{-3}$. The diameter of the PT was inferred from the TNM classification using the lower limits of the eighth edition of the TNM Staging classification Rami-Porta et al., 2015 (see Table 2).
<table><thead><tr><th>Variable</th><th>PCI, N = 47<sup>1</sup></th><th>No PCI, N = 56<sup>1</sup></th><th>Overall, N = 103<sup>1</sup></th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>Sex</td><td></td><td></td><td></td><td>0.7</td></tr><tr><td>F</td><td>20 (43%)</td><td>22 (39%)</td><td>42 (41%)</td><td></td></tr><tr><td>M</td><td>27 (57%)</td><td>34 (61%)</td><td>61 (59%)</td><td></td></tr><tr><td>Age</td><td>64 (57, 68)</td><td>65 (57, 72)</td><td>64 (57, 70)</td><td>0.6</td></tr><tr><td>Smoking</td><td></td><td></td><td></td><td>0.8</td></tr><tr><td>Non smoker</td><td>0 (0%)</td><td>1 (1.8%)</td><td>1 (1.0%)</td><td></td></tr><tr><td>Former smoker</td><td>32 (68%)</td><td>35 (63%)</td><td>67 (65%)</td><td></td></tr><tr><td>Active smoker</td><td>15 (32%)</td><td>20 (36%)</td><td>35 (34%)</td><td></td></tr><tr><td>TNM T</td><td></td><td></td><td></td><td>0.2</td></tr><tr><td>1</td><td>5 (11%)</td><td>15 (27%)</td><td>20 (19%)</td><td></td></tr><tr><td>2</td><td>9 (19%)</td><td>12 (21%)</td><td>21 (20%)</td><td></td></tr><tr><td>3</td><td>17 (36%)</td><td>16 (29%)</td><td>33 (32%)</td><td></td></tr><tr><td>4</td><td>16 (34%)</td><td>13 (23%)</td><td>29 (28%)</td><td></td></tr><tr><td>TNM N</td><td></td><td></td><td></td><td>0.023</td></tr><tr><td>0</td><td>2 (4.3%)</td><td>9 (16%)</td><td>11 (11%)</td><td></td></tr><tr><td>1</td><td>5 (11%)</td><td>6 (11%)</td><td>11 (11%)</td><td></td></tr><tr><td>2</td><td>20 (43%)</td><td>31 (55%)</td><td>51 (50%)</td><td></td></tr><tr><td>3</td><td>20 (43%)</td><td>10 (18%)</td><td>30 (29%)</td><td></td></tr><tr><td>TNM M</td><td>19 (40%)</td><td>27 (48%)</td><td>46 (45%)</td><td>0.4</td></tr><tr><td>Dissemination status</td><td>23 (49%)</td><td>30 (54%)</td><td>53 (51%)</td><td>0.6</td></tr><tr><td>Performance status</td><td></td><td></td><td></td><td>0.2</td></tr><tr><td>0</td><td>25 (53%)</td><td>20 (36%)</td><td>45 (44%)</td><td></td></tr><tr><td>1</td><td>17 (36%)</td><td>30 (54%)</td><td>47 (46%)</td><td></td></tr><tr><td>2</td><td>5 (11%)</td><td>5 (8.9%)</td><td>10 (9.7%)</td><td></td></tr><tr><td>3</td><td>0 (0%)</td><td>1 (1.8%)</td><td>1 (1.0%)</td><td></td></tr><tr><td>Surgery</td><td>1 (2.1%)</td><td>3 (5.4%)</td><td>4 (3.9%)</td><td>0.6</td></tr><tr><td>Radiotherapy</td><td>21 (45%)</td><td>27 (48%)</td><td>48 (47%)</td><td>0.7</td></tr><tr><td>Chemotherapy</td><td>26 (55%)</td><td>29 (52%)</td><td>55 (53%)</td><td>0.7</td></tr><tr><td>Radiochemotherapy</td><td></td><td></td><td></td><td>0.9</td></tr><tr><td>none</td><td>26 (55%)</td><td>29 (52%)</td><td>55 (53%)</td><td></td></tr><tr><td>sequential</td><td>8 (17%)</td><td>9 (16%)</td><td>17 (17%)</td><td></td></tr><tr><td>concurrent</td><td>13 (28%)</td><td>18 (32%)</td><td>31 (30%)</td><td></td></tr><tr><td>Adjuvant immunotherapy</td><td>5 (11%)</td><td>13 (23%)</td><td>18 (17%)</td><td>0.094</td></tr></tbody></table>

<caption>**Table 1.** BM-SCLC population characteristics</caption>

<footnote n="1">¹n (%); Median (Q1, Q3)</footnote>

<footnote n="2">²Pearson's Chi-squared test; Wilcoxon rank-sum test; Fisher's exact test</footnote>
<table><caption>**Table 2.** Conversion rules to convert TNM data to diameters and number of cells based on the eight edition of the TNM staging Rami-Porta et al., 2015.</caption><thead><tr><th>TNM T</th><th>Diameter (mm)</th><th>Number of cells</th></tr></thead><tbody><tr><td>T0</td><td>2</td><td>4188790</td></tr><tr><td>T1a</td><td>10</td><td>523598776</td></tr><tr><td>T1, T1b</td><td>25</td><td>8181230869</td></tr><tr><td>T2a</td><td>30</td><td>14137166941</td></tr><tr><td>T2, T2b</td><td>40</td><td>33510321638</td></tr><tr><td>T3</td><td>50</td><td>65449846950</td></tr><tr><td>T4</td><td>70</td><td>179594380030</td></tr></tbody></table>

## Individual fits

The model can be used to describe longitudinal individual patient data. In this setting, we identified the individual model parameters to reproduce the observed data, which consisted of information on the patient's disease at diagnosis and several follow-up visits (time points). The necessary data included the size of the PT at the time of diagnosis and, for each visit, the list of the sizes of the observed metastases.

### Growth model

In the rest of the article, the growths of the PT and the metastases were supposed to be identical. We used exponential growth defined as:

$$g(V) = a \cdot V,$$

and the Gompertz growth (Norton, 1988),

$$g(V) = V \cdot (\alpha - \beta V),$$

where the parameter $\alpha$ is the initial proliferation rate and parameter $\beta$ is a coefficient driving the exponential decrease of the specific growth rate $\frac{1}{V} \cdot \frac{dV}{dt}$, characteristic of the Gompertz model. The solution of this equation asymptotically tends to a finite value (carrying capacity) given by $K = \exp(\frac{\alpha}{\beta})$ and has the following explicit expression:

$$V(t) = \exp\left(\log K \left(1 - \exp\left(-\frac{\alpha t}{\log K}\right)\right)\right).$$

To limit the number of free parameters in the growth model, $K$ was fixed at $1 \times 10^{12}$ (corresponding to a tumor mass of 1 kg).

The effect of PT treatment was accounted for by modifying the PT growth after diagnosis. Specifically, disease control was modeled by fixing the PT size to its size at diagnosis, 50% of the diagnosis size in the case of partial response, and zero in the case of complete surgical resection or complete response.

### Dissemination model

The dissemination was assumed to follow a power law of the form:

$$d(v) = \mu \cdot v^{\gamma}$$

and we defined $r = \mu \cdot v_r^{\gamma}$ as the instant dissemination rate for a tumor with a diameter of $v_r = 10$ mm as a way to compare models with different values of $\mu$ and $\gamma$.
## Cumulative size distribution

We converted the observed individual BM size data — given as tumor diameters — to number of cells, assuming that the tumors were spherical and isotropic with a cell density of 10⁶ cells per mm³.

For each observation time $t_k$, the patient's data was represented by the cumulative observed size distribution of the metastases, i.e., the number of metastases larger than a given size. Formally, at time $t_k$, if the patient had observed metastases with respective sizes $V_j$, the cumulative size distribution function is:

$$csd(t_k, s) = \#\{V_j \le s\}.$$

Given the observed data for a patient, the model parameters for this patient $\varphi$ were determined by minimizing the mean squared error on the cumulative size distribution:

$$\varphi = \arg \min_{\varphi} \sum_{1 \le k \le K_{\max}} \int_{0}^{\infty} (csd(t_k, s) - c\tilde{sd}(t_k, s | \varphi))^{2} ds.$$

Computation of the integral above is explicit since both $csd$ and $c\tilde{sd}$ are decreasing step functions with a finite number of distinct values. The minimization procedure was performed over a discrete set of values for $\gamma$ (0.1, 0.3, 0.5, and 1) and with numerical bounds to restrain the values of $\alpha$ ([0.01, 1]) and $\mu$ ([10⁻⁷, 1] for $\gamma = 0.1$, [10⁻¹⁵, 10⁻³] for $\gamma = 0.3$, and [10⁻¹⁵, 10⁻⁶] for $\gamma = 0.5$ and $\gamma = 1$). We used the SUBPLEX algorithm (Rowan, 1990) implemented in the NLOPT library (Johnson, 2008) (using the nloptr R package version 2.2.1).

## Model training and diagnostics

The determination of individual parameters for the patient population from the AP-HM dataset used the data available at the first metastatic event $T_{BM}$ (the first follow-up time with observed brain metastases).

With an individual parameter vector $\varphi$, we plotted the model simulation, which included the simulated pre-diagnostic history and the observed data.

For a fixed time point, we confronted the model prediction and observed data for the cumulative size distribution of metastases, the number of metastases, and the metastatic burden defined as the total number of metastatic cells in the patient in all metastases.

The graphical validation of the model fit was assessed by comparing the model prediction with the data at $T_{BM}$.

## Population survival model

The mechanistic model can also be used with brain metastasis-free survival (BMFS) data. In this context, we used the (potentially censored) time between diagnosis and the first metastasis and the initial size of the PT.

## Non-linear mixed-effect model

The survival model was defined for Gompertz growth with the same parameters for PT and metastases ($g_p(V) = g(V)$).

To handle the BMFS data, we embedded the model within a non-linear mixed-effects statistical framework, where:

$$\log(Y^i) = \log(T_{\text{meta}}(V_{\text{diag}}^i; \varphi^i)) + \varepsilon^i, \quad \varepsilon^i \sim \mathcal{N}(0, \sigma^2)$$
$$
\left\{
\begin{array}{ll}
\alpha^i = \alpha_{pop} + \beta_{\alpha} \cdot C^i + \eta_{\alpha}^i, & \eta_{\alpha}^i \sim \mathcal{N}(0, \omega_{\alpha}^2) \\
\log(\mu^i) = \log(\mu_{pop}) + \beta_{\mu} \cdot C^i + \eta_{\mu}^i, & \eta_{\mu}^i \sim \mathcal{N}(0, \omega_{\mu}^2)
\end{array}
\right.
$$

with $Y^i$ being the observed (censored) time of metastasis-free survival for patient $i$, $V_{diag}^i$ the size of the PT at diagnosis, $C^i$ the vector of patient covariates, $\varphi^i = (\alpha^i, \mu^i)$ the patient's individual parameters, $\eta^i$ the individual random effects. Overall, the parameters of the population model are $\sigma, \alpha_{pop}, \mu_{pop}, \beta_{\alpha}, \beta_{\mu}, \omega_{\alpha},$ and $\omega_{\mu}$.

Population fit

We used a modified likelihood expression to account for censoring (Bigarré et al., 2023) and iden-
tified the parameters for the population model with the Stochastic approximation of expectation
maximization algorithm (SAEM) implemented in the saemix R package (Comets et al., 2017). The
estimation procedure was performed with σ fixed at 0.02 and γ fixed at 0.01. We used a normal
distribution for α and a log-normal distribution for μ.

Plausible values were chosen as initial values of the population parameters: $\alpha_{pop,0} = 0.0185$ and $\mu_{pop,0} = 0.000078$. Initial values of the inter-patient variability parameters $\omega_{\alpha,0}$ and $\omega_{\mu,0}$ were set to 0.01 and 1, respectively. The population model was trained using the entire BMFS dataset, running 4 chains in parallel for 950 iterations of the SAEM algorithm.

The variance of the estimation was evaluated with a bootstrap study with 200 bootstrap repli-
cates using the same initial values and hyperparameters for the SAEM algorithm as for the esti-
mation step.

Survival function prediction

For a patient with a given PT size at diagnosis $V_{diag}^i$ and covariate vector $C^i$, we predicted the brain
metastasis-free survival function with Monte Carlo simulation:

$$
S^i(t) = \frac{1}{p} \sum_{k=1}^{p} P\left(\exp\left(\log T_{\text{meta}}(V_{\text{diag}}^i; \varphi_k) + \varepsilon_k\right) \geq t\right)
$$

with $\eta_k$ and $\varepsilon_k$ being sampled from the population distributions.

$$
\text{The mean survival in the population was defined as } S(t) = \frac{1}{N} \sum_{i=1}^{N} S^{i}(t).
$$

Identifiability simulation study

The identifiability of the NLME model was assessed with a simulation study that generated syn-
thetic patient populations (N = 50, 100, 250, and 1000). Tumor sizes at diagnosis were sampled
from the BMFS dataset. All synthetic patients were drawn from a population with $\gamma = 0.1$ for
all individuals. The parameter $\alpha$ was normally distributed with mean $\alpha_{pop} = 0.03$ and variance
$\omega_{\alpha}^2 = 0.0001$ and $\mu$ was log-normally distributed (with median $\mu_{pop} = 0.00012$ and logarithmic vari-
ance $\omega_{\mu}^2 = 1.75$). Simulated patients with negative $\alpha$ values were excluded from the population,
effectively truncating the distribution. Brain metastasis-free survival was generated for each pa-
tient from the simulation of the model with a log-normal residual error model ($\sigma = 0.15$). These
values of the population parameters were chosen to be consistent with the data used in this study.

The censoring time was also generated from a log-normal model with a median of 36.4 months
and a standard deviation of the logarithm of 0.71.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Overview of the modeling framework. **A.** Schematic representation of the mechanistic model. The red curve represents the growth of the primary tumor. After diagnosis, the dashed lines illustrate different potential scenarios for the primary tumor depending on the treatment assumptions. The blue curves depict the size evolution of metastases: solid lines correspond to metastases emitted before diagnosis, while dashed lines represent those emitted after diagnosis assuming unchanged growth of the primary tumor. **B–G.** Illustration of the impact of modifying key modeling assumptions. All simulations are performed considering a primary tumor size of 10 mm at diagnosis. The red curve represents the growth of the primary tumor, while the blue curves indicate the growth of individual metastases. **B.** Simulation based on an exponential growth model with parameters $α = 0.12$, $μ = 0.001$, and $γ = 0.1$. **C.** Simulation using a Gompertz growth model with parameters $α = 0.044$, $μ = 0.001$, and $γ = 0.1$. **D.** Simulation using a Gompertz growth model with a reduced proliferation rate ($α = 0.02$, $μ = 0.001$, $γ = 0.1$), leading to slower tumor and metastasis growth dynamics. **E.** Simulation using a Gompertz growth model with an increased metastatic potential ($α = 0.044$, $μ = 0.005$, $γ = 0.1$), resulting in a higher number of emitted metastases. **F.** Simulation using a Gompertz growth model with increased dissemination efficiency ($α = 0.044$, $μ = 1.4 \times 10^{-11}$, $γ = 1$), with a reduced baseline emission rate to preserve the overall rate of metastasis emission ($r = 0.0074$). **G.** Simulation using a Gompertz growth model ($α = 0.044$, $μ = 0.001$, and $γ = 0.1$) with complete resection of the primary tumor at the time of diagnosis, illustrating the resulting arrest of post-diagnosis metastatic emission. dt: doubling time</figcaption></figure>
# Results

## Modeling metastasis growth and dissemination

Our proposed mathematical model of metastatic development describes the temporal evolution of the disease ranging from the initiation of the PT, early infra-clinical metastases, diagnosis, and post-diagnosis dynamics of the minimal residual disease. The flexibility of our framework allows to model different scenarios following diagnosis that correspond to different responses to treatment (Figure 1A). For example, the absence of treatment (or the absence of response to treatment) is modeled by continuing with the same growth law for the PT. Control of the tumor size through chemotherapy or immunotherapy during the early stages can be modeled as a PT keeping a constant size during the treatment, and successfully removing the primary disease by setting the tumor size to be very small (or effectively null).

Using Gompertz growth with fixed carrying capacity and power law for metastatic dissemination yields three free parameters: $\alpha$ (in $d^{-1}$) controlling growth dynamics, $\mu$ (in $d^{-1}$) the instantaneous probability for a cell in the PT to seed a new metastasis, and $\gamma$ (unitless) the parameter of the fractal scale that defines the proportion of the PT contributing to the formation of metastases.

The exponential growth (Figure 1B) does not saturate and continues to grow rapidly outside the range of possible tumor sizes. With an exponential parameter $\alpha_{exp} = 0.012$ that gives the same growth speed at a diameter of 10 mm (a doubling time of 60 days) as the Gompertz growth (Figure 1C), the PT takes longer to reach this size. The first metastasis also takes more time to seed; after that, the interval between successive metastases quickly increases.

Figure 1D-F shows the qualitative effects of each parameter for the model with Gompertz growth and power law dissemination. For a Gompertz growth with a lower value of $\alpha$ (Figure 1D, $\alpha = 0.02$, the doubling time when the diameter reaches 10 mm is 133 days), the PT takes longer to reach the same size (60 months to reach 10 mm with $\alpha = 0.02$, 27 months to reach 10 mm for $\alpha = 0.044$), leaving more time for metastases to appear. On the other hand, for a higher value of $\alpha$ (Supplementary Figure S3A, $\alpha = 0.1$) the tumors grow faster and fewer metastases appear when the PT reaches a size of 10 mm.

Changing the value of $\mu$ (Supplementary Figure S3B $\mu = 0.0004$, Figure 1E $\mu = 0.005$) while keeping $\gamma$ constant impacts the scale of the metastasis number function $N$, but not its shape. Changing $\gamma$ and $\mu$ simultaneously to maintain a constant instantaneous dissemination rate of 0.0074 cell/d for a tumor diameter of 10 mm shows that the value of $\gamma$ strongly impacts the shape of $N$. A small value of $\gamma$ (Supplementary Figure S3C, $\gamma = 0.01$) induces an almost constant time between the birth of two successive metastases, while a large value (Figure 1F, $\gamma = 1$) produces a long time before the first metastasis and then an exponentially decreasing time between the birth of the metastases.

Figure 1G shows the effect of a complete resection of the primary tumor (PT) when it reaches a diameter of 10 mm, while Supplementary Figure S3D illustrates the impact of a 90% reduction, simulating alternative treatment outcomes for the PT. After partial resection, new metastases appear more slowly and at a constant rate. Total resection prevents the birth of new metastases.

## Brain metastasis-free survival

Population-level parameters of the NLME model were estimated using the SAEM algorithm, with resulting estimates presented in Table 3. The training dataset was augmented with "phantom patients" representing left-censored observations at the time of diagnosis to account for patients presenting with brain metastases at initial diagnosis.

The fitted population model demonstrated adequate descriptive performance successfully re-
<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Nonlinear Mixed-Effects Population Model.* Fit of the population model to the BMFS dataset without any covariates, assuming a normal distribution for α and a log-normal distribution for μ. The residual error model is log-normal and adjusted for censoring when applicable. **A.** Comparison between the Kaplan-Meier estimate of the full BMFS dataset (including phantom left-censored patients) and the model-predicted population survival function (including negative relapse times). **B.** Comparison between the Kaplan-Meier estimate of the true BMFS dataset (real patients only) and the model-predicted population survival function (excluding non-positive relapse times). **C-F.** Bootstrap analysis (200 replicates) assessing parameter identifiability. Each panel shows the bootstrap density of a population parameter, with the solid vertical line representing the estimate from the full dataset and the dashed lines indicating the 90% bootstrap confidence interval. **C.** Bootstrap density of α<sub>pop</sub> **D.** Bootstrap density of μ<sub>pop</sub> **E.** Bootstrap density of ω<sub>α</sub> **F.** Bootstrap density of ω<sub>μ</sub></figcaption></figure>
<table><caption>**Table 3.** Estimands of the population parameters
BMFS dataset</caption><thead><tr><th>Parameters</th><th>value</th><th>R.S.E.</th></tr></thead><tbody><tr><td>α<sub>pop</sub> (per day)</td><td>4.60 × 10<sup>-2</sup></td><td>13.99%</td></tr><tr><td>ω<sup>2</sup><sub>α</sub></td><td>3.95 × 10<sup>-4</sup></td><td>45.00%</td></tr><tr><td>μ<sub>pop</sub> (per day)</td><td>7.03 × 10<sup>-4</sup></td><td>12.63%</td></tr><tr><td>ω<sup>2</sup><sub>μ</sub></td><td>0.68</td><td>70.66%</td></tr><tr><td>γ</td><td>0.01</td><td>(Fixed)</td></tr><tr><td>σ</td><td>0.02</td><td>(Fixed)</td></tr></tbody></table>

producing the overall population survival function (Figure 2 A) while remaining within the 90% confidence bounds of the Kaplan-Meier estimator. However, systematic deviations were observed, with underestimation of survival probability during the initial phase (<1 year) and modest overestimation of long-term survival rates. To evaluate model performance specifically for patients without brain disease at diagnosis, survival functions were also conditioned on disease-free status at baseline (Figure 2 B). While the model maintained reasonable overall concordance with the observed brain disease-free survival patterns, it failed to adequately capture the initial phase dynamics, characterized by slower progression rates in the early post-diagnosis period.

The precision of the parameter estimation was assessed through 200 bootstrap replicates (Figure 2 C-F). The fixed-effects parameters α<sub>pop</sub> (Figure 2 C) and μ<sub>pop</sub> (Figure 2 E) demonstrated robust estimation properties with low relative standard errors (RSE), all values falling below the conventional 30% threshold for acceptable fixed-effects parameter estimation (Table 3). Between-subject variance components exhibited greater estimation uncertainty: ω<sup>2</sup><sub>α</sub> (Figure 2 D) achieved acceptable precision with RSE values below the established 50% threshold for variance parameter estimation, while ω<sup>2</sup><sub>μ</sub> (Figure 2 F) showed considerably higher estimation variability (RSE = 70.66%), exceeding conventional acceptability criteria and indicating reduced confidence in this parameter estimate.

## Identification of individual parameters

METAMATS allows the identification of individual parameters from PT size data at diagnosis and metastasis sizes at the time t<sub>BM</sub>. We used METAMATS on the AP-HM datasets for patients who have had a brain metastasis relapse (N = 40 patients). Simultaneous estimation of the three parameters α, μ, and γ was not computationally feasible as no stable solution could be found using the SUBPLEX algorithm, which led us to estimate α and μ for each patient while keeping γ fixed. We compared four possible values for γ (0.1, 0.3, 0.5, and 1) corresponding to different hypotheses in the subset of PT cells that can disseminate. A value of γ = 0.1 gave significantly better fits for the whole population; values of γ = 0.3, 0.5, or 1 led to a similar higher error level (Figure 4 A).

At the individual level, we found that for a large number of patients γ = 0.1 was indeed the best fit (17 patients out of 40, Supplementary Figure S4), but still more than half of the population was estimated to have a different value of γ (6 patients with γ = 0.3, 10 patients with γ = 0.5 and 7 patients with γ = 1). This result emphasizes the potential inter-individual variability in terms of metastatic aggressiveness.

The model was able to fit data from patients with different metastasis size profiles, to run simu-
<figure><img src="image_4.png" /><figcaption>**Figure 3. Individual patient fits and calibration of the metastatic size distribution model.** Model parameters were estimated through minimization of squared error between observed and predicted cumulative size distributions for each patient. **A-C.** Model simulations using patient-specific parameter estimates demonstrating predicted metastatic growth dynamics for representative cases: A. Patient 8, B. Patient 23, and C. Patient 42. **D-F.** Comparative analysis of size distributions for observable metastases (>3 mm) between empirical observations and model predictions at the time of metastatic relapse: D. Patient 8, E. Patient 23, and F. Patient 42. **G-I.** Direct comparison of individual metastatic lesion sizes between observed and model-predicted observable metastases at metastatic relapse: G. Patient 8, H. Patient 23, and I. Patient 42. **J.** Model calibration assessment for the number of visible metastases, displaying predicted versus observed metastatic counts. Data points represent discrete integer values with random positioning applied to enhance visualization of the underlying distribution. **K.** Model calibration assessment for visible metastatic burden, comparing predicted total metastatic cell number against observed values.</figcaption></figure>
lations of the natural history of the disease for visible and occult metastases, and to predict the age of the disease at diagnosis, as well as metastatic aggressiveness metrics such as the dissemination rate for a 10 mm PT. Patient 8 (Figure 3 A) had a tumor of 40 mm diameter at diagnosis and a brain metastasis of approximately 15 mm at the first brain metastatic event (12 months after diagnosis); the estimated parameters show a small $\alpha = 0.016$, leading to an estimated age of PT at diagnosis of 10 years and a doubling time at diagnosis of 394 days, $\gamma = 0.1$, $\mu = 0.00024$ and $r = 0.0024$. This set of parameters predicts that 2 tumors, the first of size around 15 mm and the second of size around 5 mm, should be present at the time of the brain metastatic event (Figure 3 B–C). Patient 23 (Figure 3 D) had a smaller 25 mm PT, but developed four BMs 7 months after diagnosis, with 2 between 8 mm and 5 mm and 2 smaller than 5 mm. For this patient, $\alpha = 0.017$ was similar to that of patient 8, giving a similar growth of PT (estimated age around 8 years, doubling time at diagnosis of 253 days) but $\gamma = 0.3$ and $\mu = 5 \times 10^{-5}$ ($r = 0.048$) led to a different metastatic scenario with 3 tumors predicted by the model (Figure 3 E–F). Patient 42 (Figure 3 G) had the largest PT (50 mm) and developed, 12 months after diagnosis, three 9 mm BM and six BM smaller than 8 mm. Their PT was also the youngest (3.7 years) with a large $\alpha = 0.047$ and a doubling time at diagnosis of 172 days. Metastatic parameters $\gamma = 1$ and $\mu = 1.7 \times 10^{-12}$ ($r = 0.014$) induce rapid metastatic dissemination with the prediction of 5 small visible BM (less than 5 mm), 3 BM in the 5 mm to 8 mm range, and 2 BM larger than 9 mm (Figure 3 H–I).

Due to differences in estimated values of $\alpha$ and PT at diagnosis, predicted ages of the PT ranged from 2 to over 12 years, with a mean of 6.4 years and a standard deviation of 3.0 years (Supplementary Figure S5).

Beyond the comprehensive size distribution analysis, METAMATS enables the prediction of total metastatic burden and metastasis count at any given time. The model’s prediction of visible metastatic burden was compared against observed metastatic burden at $T_{BM}$ in the AP-HM patient cohort, demonstrating excellent concordance between predicted and actual visible metastatic burden (Figure 3J). Additionally, METAMATS facilitated comparison between predicted and observed numbers of visible metastases (Figure 3K), revealing accurate model performance for this parameter, albeit with a systematic tendency toward overestimation, particularly pronounced at higher metastasis counts.

## Clinical application: impact of prophylactic cranial irradiation

We compared the distribution of the individual parameters in the PCI and non-PCI groups to see how the model would reflect this therapeutic intervention. No significant differences were found in the tumor growth parameter $\alpha$ between the two groups ($p = 0.29$ Wilcoxon rank sum test, Figure 4 A), while the fractal parameter $\gamma$ was found to be greater in the PCI group than in the non-PCI group ($p = 0.044$ Fisher’s exact test, Figure 4 B).

Direct comparison of $\mu$ between patients is difficult because the scale of the $\mu$ parameter depends on the value of $\gamma$ (Supplementary Figure S6). We thus used the distribution of the parameter $r$ that incorporates both $\gamma$ and $\mu$ to compare all patients on the same scale, and $r$ was found to be significantly lower in the PCI versus non-PCI group (Wilcoxon rank sum test, p-value 0.011). The joint distribution of $\alpha$ and $r$ (Figure 4D) did not reveal a clear separation of the two groups.

The NLME model framework facilitates the investigation of covariate effects on BMFS. This modeling approach enables independent covariate influence on the distribution of individual model parameters. We incorporated the PCI status as a covariate on both $\alpha$ and $\mu$ distributions and trained this NLME model on the BMFS dataset yielding the parameter estimates presented in Table 4. Integration of PCI status yielded superior model performance, with improved log-likelihood (-2LL
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Influence of prophylactic cranial irradiation (PCI) status on model parameter estimation and predictive performance. Individual parameter distributions are stratified by PCI status, with patients without PCI represented in red and patients with PCI represented in blue. **A–C.** Distribution analysis of patient-specific parameter estimates stratified by PCI status: A. growth parameter α, B. dissemination rate r, and C. fractal scale parameter γ. **D.** Bivariate distribution of growth parameter α and dissemination rate r with color-coding by PCI status (red: no PCI; blue: PCI present). **E.** Kaplan–Meier visual predictive check for the nonlinear mixed-effects (NLME) population model incorporating PCI status as a covariate for the distribution of μ. Solid lines represent model-predicted means, dotted lines indicate 90% prediction intervals, with stratification by PCI status (red: no PCI group; blue: PCI group). The visual predictive check assesses model adequacy through comparison of observed survival patterns against model-generated predictions under varying PCI conditions.</figcaption></figure>
<table><caption>**Table 4.** Population parameter estimates for the NLME model incorporating PCI status effects on α and μ parameters.</caption><thead><tr><th>Parameters</th><th>value</th><th>R.S.E.</th></tr></thead><tbody><tr><td>α<sub>pop</sub> (per day)</td><td>2.50 × 10<sup>-2</sup></td><td>34.12%</td></tr><tr><td>ω<sup>2</sup><sub>α</sub></td><td>1.49 × 10<sup>-4</sup></td><td>79.05%</td></tr><tr><td>μ<sub>pop</sub> (per day)</td><td>7.90 × 10<sup>-4</sup></td><td>22.91%</td></tr><tr><td>ω<sup>2</sup><sub>μ</sub></td><td>0.56</td><td>32.42%</td></tr><tr><td>β<sub>PCI,α</sub></td><td>0.02</td><td>27.57%</td></tr><tr><td>β<sub>PCI,μ</sub></td><td>-0.38</td><td>62.65%</td></tr><tr><td>γ</td><td>0.01</td><td></td></tr><tr><td>σ</td><td>0.02</td><td></td></tr></tbody></table>

BMFS dataset

= 1560) compared to the covariate-free model (-2LL = 1599). The Bayesian Information Criterion was similarly enhanced in the complete model (BIC = 1598 vs 1626), at the same time the inter-individual variance was reduced both for α ($\omega_{\alpha}^2 = 3.95 \times 10^{-4}$ without PCI and $\omega_{\alpha}^2 = 1.49 \times 10^{-4}$ with PCI as a covariate) and for μ ($\omega_{\mu}^2 = 0.68$ without PCI and $\omega_{\mu}^2 = 0.56$ with PCI as a covariate).

Parameter estimation for this enhanced model remained robust, with $\alpha_{pop}$ exhibiting a relative standard error (RSE) marginally exceeding 30% (Supplementary Figure S7 A) and $\mu_{pop}$ demonstrating adequate precision (Supplementary Figure S7 B). Inter-individual variability parameters maintained the previously observed estimation imbalance, with $\omega_{\alpha}$ displaying substantial RSE exceeding 70% (Supplementary Figure S7 C), while $\omega_{\mu}$ achieved satisfactory estimation precision (Supplementary Figure S7 D). The PCI status effect on α was well-characterized, with $\beta_{PCI,\alpha}$ RSE below 30% (Supplementary Figure S7 E). Conversely, the PCI effect on μ proved more challenging to estimate accurately within our dataset, resulting in $\beta_{PCI,\mu}$ RSE exceeding 70% (Supplementary Figure S7 F). Despite estimation uncertainty, the bootstrap distribution (Supplementary Figure S7) provided strong evidence for significant PCI status effects on μ, with fewer than 5% of bootstrap estimates for $\beta_{PCI,\mu}$ yielding strictly negative values ($p < 0.00025$, Wald test).

## Simulations and identifiability

Investigation of the impact of sample size on the model identifiability of the NLME time-to-relapse model was assessed through simulation studies using synthetic patient populations. Model parameters were estimated across populations of varying sizes (n = 50, 100, 250, 500, and 1000) to evaluate the relationship between sample size and estimation precision. Simulation parameters were set to $\alpha = 0.04$, $\mu = 0.0007$, $\omega_{\alpha}^2 = 0.0004$, $\omega_{\mu}^2 = 0.7$, $\sigma = 0.02$, and $\gamma = 0.01$ to replicate conditions analogous to those observed in the BMFS population estimation. Censoring was simulated using a lognormal distribution with parameters derived from the observed censoring distribution in the BMFS dataset (median censoring time = 1109 d, standard error of log-censoring time = 0.71).

The population growth parameter $\alpha_{pop}$ exhibited minimal bias across all sample sizes (Figure 5 A), with estimates converging to approximately 0.046 (relative bias: 15%). Conversely, the population dissemination parameter $\mu_{pop}$ demonstrated negligible bias (relative bias: 7%), with the true value consistently contained within the 90% confidence interval (Figure 5 B). The inter-individual variability parameter $\omega_{\alpha}$ proved to be the least identifiable parameter (Figure 5 C), exhibiting con-
A

<figure><img src="image_10.png" /></figure>

B

<figure><img src="image_11.png" /></figure>

C

<figure><img src="image_12.png" /></figure>

D

<figure><img src="image_13.png" /></figure>

Simulation parameter median 90% interval

<caption>**Figure 5.** Identifiability assessment of the nonlinear mixed-effects (NLME) population model through simulation-based parameter recovery analysis. Synthetic datasets were generated using the established model framework across varying sample sizes (N = 50, 100, 250, 500, 1000 subjects). For each sample size configuration, 100 independent simulations were conducted, with subsequent re-estimation of population model parameters performed to assess parameter identifiability and estimation precision. Solid lines represent the true parameter values employed in the data generation process, dashed lines indicate the median of parameter estimates across 100 replicates for each sample size, and dotted lines demarcate the 90th percentile boundaries of the estimated parameter distributions. **A.** Population-level growth parameter $\alpha_{pop}$. **B.** Population-level location parameter $\mu_{pop}$. **C.** Between-subject variance component for the growth parameter $\omega_{\alpha}^2$. **D.** Between-subject variance component for the location parameter $\omega_{\mu}^2$. The convergence patterns demonstrate the asymptotic properties of maximum likelihood estimators and quantify the sample size requirements for reliable parameter identification within the proposed modeling framework.</caption>
sistent substantial bias with estimates converging to 0.0002 (relative bias: 48%). In contrast, the
inter-individual variability parameter for μ demonstrated excellent identifiability (*Figure 5 D*) with
negligible bias.

As anticipated, estimation variance decreased with increasing sample size. A sample size of N = 250 appeared sufficient to achieve low estimator dispersion, with no substantial variance improvement observed when increasing from N = 500 to N = 1000. Given that the dataset comprises censored time-to-event data, individual contributions to information content vary, with only non-censored patients providing complete likelihood information.

Discussion

This article presents a comprehensive framework for modeling metastatic size distributions, with a demonstration of its clinical application to brain metastases in small-cell lung cancer (SCLC) and the impact of prophylactic cranial irradiation. The framework is implemented in the R package METAMATS, which provides extensibility for incorporating novel modeling hypotheses and diverse clinical scenarios. The METAMATS package integrates computational tools for parameter inference from both individual patient data and population-level time-to-event data, alongside capabilities for simulation and visualization.

Modeling

METAMATS uses a discrete formulation of a partial differential equation to model the temporal evolution of metastatic size distributions, based on an original paper by Iwata et al. Iwata et al., 2000. This modeling framework encompasses a broad family of models capable of representing diverse biological hypotheses regarding the metastatic process Benzekry et al., 2016; Bigarré et al., 2023; Bilous et al., 2019; Nicolò et al., 2020. The framework provides substantial flexibility in defining the fundamental model functions that govern tumor growth dynamics and instantaneous dissemination from the primary tumor. Users can implement custom functions tailored to their specific modeling assumptions. The default implementation incorporates either Gompertz or exponential growth laws for both primary tumors and metastases, with optional incorporation of a fixed dormancy period. The dissemination function is implemented by default as a power-law relationship with primary tumor size. Treatment effects on primary tumors and metastases are incorporated through time-dependent modifications of the growth dynamics. The package can accommodate surgical interventions, including partial or complete resection of the primary tumor at diagnosis. Other types of non-constant treatment responses could also be incorporated into the model.

The current implementation assumes independence between individual tumor growth dynam-
ics and enforces uniform growth kinetics for all metastases as a function of time since seeding.
Future developments could include the incorporation of time-dependent differential equations
to define growth dynamics, thereby expanding the range of biological phenomena that can be
modeled. This extension would enable the representation of host-tumor interactions, including
the influence of the local tumor microenvironment and systemic immune system. An additional
prospective enhancement involves the incorporation of metastatic site specificity, as the current
model remains organ-agnostic with respect to metastatic localization.

Inference and predictions

In metastatic disease studies, two distinct data types are available, necessitating different infer-
ence procedures. The primary model output consists of the complete size distribution of metas-
tases as a function of time. The most direct data type comprises longitudinal records of metastatic count and size for individual patients. However, this presents two challenges: such data are difficult to acquire in practice, typically requiring specialized radiological studies with systematic enumeration and measurement of individual tumors; additionally, no established statistical framework exists for these data, as model predictions and observations constitute non-normalized integer-valued functions lacking classical probability distributions.

To address these challenges, we implemented a minimization approach for estimating individual patient parameters using a criterion based on squared differences of cumulative distribution functions. While this method enables individual parameter estimation, the absence of a rigorous statistical framework precludes population-level model assessment. Furthermore, simultaneous estimation of all three free parameters (α, μ, and γ) proved computationally intractable, necessitating a modified approach wherein γ was fixed at discrete values (0.1, 0.3, 0.5, 1.0) with separate estimation of α and μ for each γ value.

A second data type applicable to this modeling framework is metastasis-free survival, which provides limited information at the individual level but becomes informative at the population level. Within our modeling framework, metastatic relapse is defined as the time at which the first metastasis exceeds a specified size threshold $V_{vis}$. This definition enabled the construction of a robust non-linear mixed-effects (NLME) model for population-level parameter inference.

Due to technical constraints, the current NLME implementation was restricted to a specific model configuration wherein primary tumor and metastatic growth follow identical Gompertz kinetics with fixed carrying capacity, and dissemination was governed by a power law function of primary tumor size. The model additionally assumed complete resection of the primary tumor at diagnosis, thereby precluding post-diagnostic metastatic seeding. Future developments could encompass NLME implementation for longitudinal metastatic size data and extension of the time-to-event framework to accommodate alternative growth and dissemination functions, including the impact of pre- and post-surgical therapy Benzekry et al., 2024.

## Simulations

As a mechanistic model, METAMATS enables direct simulation of patient trajectories given individual parameters. This capability allows reconstruction of pre-diagnostic disease history and estimation of the occult metastatic burden at diagnosis. Forward simulations can predict disease progression for individual patients and evaluate different therapeutic response scenarios and treatment strategies. This functionality positions the model as a potential tool for personalized metastatic risk assessment. We demonstrated that the model generates realistic patient trajectories suitable for *in silico* studies investigating treatment effects and informing clinical trial design.

Simulation studies provided the foundation for comprehensive identifiability analysis of the NLME population model. Synthetic populations of varying sizes were generated using population parameters, initial tumor sizes, and censoring patterns comparable to those observed in clinical data. Parameter recovery was assessed using the estimation procedure with $\gamma$ and $\sigma$ (observation error) fixed at values used in the clinical dataset analysis.

For the growth-related population parameters $\alpha_{pop}$ and $\omega_{\alpha}$, similar sample size effects on variance were observed; however, both parameters exhibited systematic bias. Two potential mechanisms may explain this phenomenon: (i) the likelihood landscape may contain multiple local minima for the parameter pair $(\alpha_{pop}, \omega_{\alpha})$, with a non-global minimum attracting a substantial proportion of estimates; (ii) approximately 2% of simulated patients exhibited negative $\alpha$ values under the simulation parameters, necessitating their exclusion and effectively truncating the population
distribution. This truncation likely contributed to the observed bias in the identifiability analysis. Despite the systematic nature of this bias, the magnitude remained relatively modest (15% relative error for $\alpha_{pop}$ and 50% relative error for $\omega_\alpha$).

## Application: clinical insights for the management of brain metastases in SCLC patients

Prophylactic cranial irradiation represents the current standard of care for patients without detectable brain metastases at diagnosis; however, recent concerns regarding the benefit-risk profile of this intervention have emerged due to its substantial impact on quality of life (Heredia et al., 2025; Zeng et al., 2020). The study population comprised patients who either received or did not receive PCI following successful primary tumor treatment.

For a subset of the population (N = 41), records detailing brain metastasis number and size at relapse were available. Individual parameters were estimated for this subpopulation using cumulative size distribution criteria, enabling per-patient parameter estimation describing both time to first metastatic event and brain metastasis size distributions. The $\gamma$ parameter was estimated consistently closer to zero than to one, with most patients demonstrating optimal fit at $\gamma = 0.1$; however, substantial variability was observed in $\gamma$ distribution, with over half of patients achieving superior fit with $\gamma \neq 0.1$ (ranging from 0.3 to 1.0). This suggests that in most cases, only a small fraction of the primary tumor contributes to brain metastasis formation. Such low fractal dimension of the metastasis-competent primary tumor subpopulation could be attributed to limited tumor vascularization (assuming uniform metastatic potential among primary tumor cells proximal to blood vessels) or restricted development of metastatic phenotype within only a small primary tumor fraction, potentially involving a combination of both mechanisms.

The individual parameter $\alpha$ was estimated to range from 0.014 to 0.091, with a mean of 0.034, consistent with the wide range of doubling rates for SCLC primary tumors reported in the literature (Brigham et al., 1978). From this information, we estimated tumor age since initiation to range from 2 to 12 years (mean 6.7 years), which is biologically plausible.

For the dissemination component of the model, parameter $\mu$ was estimated to have a very large inter-patient variability, ranging from $1.9 \times 10^{-13}$ to $6.6 \times 10^{-4}$ (mean $1.1 \times 10^{-4}$), but with strong dependence on the $\gamma$ values. When $\gamma = 0.1$ the median $\mu$ was $2.38 \times 10^{-4}$; when $\gamma = 0.3$ the median $\mu$ was $6.27 \times 10^{-6}$; when $\gamma = 0.5$ the median $\mu$ was $7.68 \times 10^{-8}$; and when $\gamma = 1$ the median $\mu$ was $3.32 \times 10^{-13}$ (Supplementary Figure S6). Better comparison was achieved using $r$, the dissemination rate of successfully established metastatic cells at a fixed diameter of 10 mm. In our population, $r$ ranged from $9.8 \times 10^{-5}$ cell/d to $2.12 \times 10^{-2}$ cell/d (median $= 1.77 \times 10^{-3}$ cell/d) and was much more homogeneous across $\gamma$ values (except potentially for $\gamma = 1$). This means that the most aggressive tumors were estimated to seed a new metastasis every 47 days (at constant size) and the least aggressive every 28 years (median 565 days).

Distribution of parameters gamma and $r$ were found to be significantly different for patients who received PCI, while values of alpha were comparable in both groups. This suggests that PCI does not much alter the growth of the brain metastases but rather acts as if the dissemination from the PT was lower, which could correspond to the destruction of part of the brain tumor population. One of the drawbacks of the individual fitting approach is that a lot of data per individual is needed to obtain a precise estimation of the individual parameters, and a large number of patients is needed to see small or medium effects of PCI on the parameters.

Using the metastatic-free survival data, population model parameters could be reliably estimated only when incorporating information regarding the proportion of patients presenting with
metastatic disease at diagnosis. By design, the database excluded patients with brain metastases at diagnosis, resulting in poor convergence and unstable parameter estimates when using this dataset alone. To address this limitation, the original cohort was augmented with phantom patients exhibiting left-censored relapse times at diagnosis. The proportion of supplementary patients was calibrated to match metastatic prevalence reported in SEER database analyses (Morgensztern et al., 2012). Following the inclusion of left-censored patients, population parameters achieved acceptable estimation precision, with the exception of $ω_{μ''}$, which exhibited relative standard error exceeding 70%. The necessity of population augmentation supports our choice of population distribution for individual parameters, as fitting without left-censored patients corresponds to estimating a truncated (log-)normal distribution from data originating from a complete distribution.

The NLME approach confirmed that a very small value for $γ$ was best suited to describe brain metastasis relapse in this population, consistent with the cancer stem cell-like hypothesis where only a small pool of cells develops metastatic capabilities and the size of this pool is largely independent of PT size. The estimated distribution for $α$ matched the values obtained in the individual approach and values of $μ$ were of the same order of magnitude as for the individual approach when $γ = 0.1$, showing that both methods are compatible and can provide complementary information.

The population model also brought quantitative estimates of the impact of PCI. Including PCI status as a linear covariate on the log-distribution of $μ$ led to a significant reduction of the dissemination rate for patients who received PCI by a factor of $exp(-1.22) \approx 0.30$. At the population level, the inclusion of the PCI status covariate discriminated the two sub-populations in terms of BMFS. From a biological point of view, there is little evidence that brain irradiation could directly impact the metastatic capability of PT cells. The impact on $μ$ should probably be viewed as a proxy for an effect on part of the metastatic process not described by our model, such as a delay between metastatic seeding and growth.

This study not only introduces the METAMATS software but also represents a proof of concept for the practical application of our mechanistic modeling framework. It is the first mechanistic modeling study of BMFS in SCLC patients and demonstrates interesting clinical results on the impact of PCI in this context. Upon further validation, METAMATS could provide a valuable tool to discriminate patients and personalize treatments based on predicted metastatic risks.

## References

Achrol, A. S., R. C. Rennert, C. Anders, R. Soffietti, M. S. Ahluwalia, L. Nayak, S. Peters, N. D. Arvold, G. R. Harsh, P. S. Steeg, and S. D. Chang (2019). "Brain metastases". In: *Nature Reviews Disease Primers* 5.1, p. 5. DOI: 10.1038/s41572-018-0055-y.

Baratchart, E., S. B. *, A. B. *, T. C. *, L. S. Cooley, R. Pineau, E. J. Ribot, O. Saut, and W. Souleyreau (2015). "Computational Modelling of Metastasis Development in Renal Cell Carcinoma". In: *PLoS Computational Biology* 11.11. Ed. by N. L. Komarova, e1004626. DOI: 10.1371/journal.pcbi.1004626. URL: http://dx.plos.org/10.1371/journal.pcbi.1004626.

Barnholtz-Sloan, J. S., A. E. Sloan, F. G. Davis, F. D. Vigneau, P. Lai, and R. E. Sawaya (2004). "Incidence Proportions of Brain Metastases in Patients Diagnosed (1973 to 2001) in the Metropolitan Detroit Cancer Surveillance System". In: *Journal of Clinical Oncology* 22.14. PMID: 15254054, pp. 2865-2872. DOI: 10.1200/JCO.2004.12.149. eprint: https://ascopubs.org/doi/pdf/10.1200/JCO.2004.12.149. URL: https://ascopubs.org/doi/abs/10.1200/JCO.2004.12.149.
Benzekry, S., M. Mastri, C. Nicolò, and J. M. L. Ebos (2024). "Machine-learning and mechanistic modeling of metastatic breast cancer after neoadjuvant treatment". In: *PLOS Computational Biology* 20.5, pp. 1-20. DOI: 10.1371/journal.pcbi.1012088. URL: https://doi.org/10.1371/journal.pcbi.1012088.

Benzekry, S., A. Tracz, M. Mastri, R. Corbelli, D. Barbolosi, and J. M. Ebos (2016). "Modeling Spon-
taneous Metastasis following Surgery: An In Vivo-In Silico Approach". In: *Cancer Research* 76.3,
pp. 535-547. ISSN: 0008-5472. DOI: 10.1158/0008-5472.CAN-15-1389. eprint: https://aacrjournals.
org/cancerres/article-pdf/76/3/535/2743651/535.pdf. URL: https://doi.org/10.1158/0008-5472.CAN-
15-1389.

Bigarré, C., F. Bertucci, P. Finetti, G. Macgrogan, X. Muracciole, and S. Benzekry (2023). "Mechanistic modeling of metastatic relapse in early breast cancer to investigate the biological impact of prognostic biomarkers". In: *Computer Methods and Programs in Biomedicine* 231, p. 107401. ISSN: 0169-2607. DOI: https://doi.org/10.1016/j.cmpb.2023.107401. URL: https://www.sciencedirect.com/science/article/pii/S0169260723000688.

Bilous, M., C. Serdjebi, A. Boyer, P. Tomasini, C. Pouypoudat, D. Barbolosi, F. Barlesi, F. Chomy, and S. Benzekry (2019). "Quantitative mathematical modeling of clinical brain metastasis dynamics in non-small cell lung cancer". In: *Scientific Reports* 9.1, p. 13018. ISSN: 2045-2322. DOI: 10.1038/s41598-019-49407-3. URL: https://doi.org/10.1038/s41598-019-49407-3.

Bocci, F., M. Kumar Jolly, and J. N. Onuchic (2019). "A Biophysical Model Uncovers the Size Distribution of Migrating Cell Clusters across Cancer Types". In: *Cancer Research* 79.21, pp. 5527-5535. ISSN: 0008-5472. DOI: 10.1158/0008-5472.CAN-19-1726. eprint: https://aacrjournals.org/cancerres/article-pdf/79/21/5527/2871952/5527.pdf. URL: https://doi.org/10.1158/0008-5472.CAN-19-1726.

Brigham, B. A., P. A. Bunn, J. D. Minna, M. H. Cohen, D. C. Ihde, and S. E. Shackney (1978). "Growth rates of small cell bronchogenic carcinomas". In: *Cancer* 42.6, pp. 2880-2886. ISSN: 0008-543X. DOI: 10.1002/1097-0142(197812)42:6<2880::aid-cncr2820420650>3.0.co;2-x.

Comets, E., A. Lavenu, and M. Lavielle (2017). "Parameter Estimation in Nonlinear Mixed Effect Models Using saemix, an R Implementation of the SAEM Algorithm". In: *Journal of Statistical Software* 80.3. DOI: 10.18637/jss.v080.i03.

Dillekås, H., M. S. Rogers, and O. Straume (2019). "Are 90% of deaths from cancer caused by metastases?" In: *Cancer Medicine* 8.12, pp. 5574-5576. DOI: https://doi.org/10.1002/cam4.2474. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cam4.2474. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cam4.2474.

Faivre-Finn, C., M. Snee, L. Ashcroft, W. Appel, F. Barlesi, A. Bhatnagar, A. Bezjak, F. Cardenal, P. Fournel, S. Harden, C. L. Pechoux, R. McMenemin, N. Mohammed, M. O'Brien, J. Pantarotto, V. Surmont, J. P. V. Meerbeeck, P. J. Woll, P. Lorigan, F. Blackhall, and C. S. Team (2017). "Concurrent once-daily versus twice-daily chemoradiotherapy in patients with limited-stage small-cell lung cancer (CONVERT): an open-label, phase 3, randomised, superiority trial". In: *The Lancet Oncology* 18.8, pp. 1116-1125. ISSN: 1470-2045. DOI: 10.1016/s1470-2045(17)30318-2.

Gallaher, J., M. Strobl, J. West, R. Gatenby, J. Zhang, M. Robertson-Tessi, and A. R. Anderson (2023). "Intermetastatic and Intrametastatic Heterogeneity Shapes Adaptive Therapy Cycling Dynamics". In: *Cancer Research* 83.16, pp. 2775-2789. ISSN: 0008-5472. DOI: 10.1158/0008-5472.CAN-22-2558. eprint: https://aacrjournals.org/cancerres/article-pdf/83/16/2775/3355407/2775.pdf. URL: https://doi.org/10.1158/0008-5472.CAN-22-2558.

Haeno, H., M. Gonen, M. B. Davis, J. M. Herman, C. A. Iacobuzio-Donahue, and F. Michor (2012). "Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Opti-
mum Treatment Strategies". In: *Cell* 148.1-2, pp. 362-375. ISSN: 0092-8674. DOI: 10.1016/j.cell.2011.11.060.

Heredia, L., M. Murcia-Mejía, and M. Torrente (2025). "Dataset on the quality of life and neurocognitive effects of prophylactic cranial irradiation, with and without hippocampal avoidance, in small-cell lung cancer patients". In: *Data in Brief* 58, p. 111197. ISSN: 2352-3409. DOI: 10.1016/j.dib.2024.111197.

Iwata, K., K. Kawasaki, and N. Shigesada (2000). "A Dynamical Model for the Growth and Size Distribution of Multiple Metastatic Tumors". In: *Journal of Theoretical Biology* 203.2, pp. 177-186. ISSN: 0022-5193. DOI: https://doi.org/10.1006/jtbi.2000.1075. URL: https://www.sciencedirect.com/science/article/pii/S0022519300910750.

Johnson, S. G. (2008). *The NLopt nonlinear-optimization package*. URL: https://github.com/stevengj/nopt.

Lavielle, M. (2014). *Mixed Effects Models for the Population Approach*. Chapman and Hall/CRC.

Li, N., Y. Chu, and Q. Song (2021). "Brain Metastasis in Patients with Small Cell Lung Cancer". In: *Int J Gen Med* 14, pp. 10131-10139.

Mani, K., D. Deng, C. Lin, M. Wang, M. L. Hsu, and N. G. Zaorsky (2024). "Causes of death among people living with metastatic cancer". In: *Nature Communications* 15.1, p. 1519. ISSN: 2041-1723. DOI: 10.1038/s41467-024-45307-x. URL: https://doi.org/10.1038/s41467-024-45307-x.

Masotti, M., N. Osher, J. Eliason, A. Rao, and V. Baladandayuthapani (2023). "DIMPLE: An R package to quantify, visualize, and model spatial cellular interactions from multiplex imaging with distance matrices". In: *Patterns* 4.12, p. 100879. ISSN: 2666-3899. DOI: 10.1016/j.patter.2023.100879.

McDonald, T. O., Y.-C. Cheng, C. Graser, P. B. Nicol, D. Temko, and F. Michor (2023). "Computational approaches to modelling and optimizing cancer treatment". In: *Nature Reviews Bioengineering* 1.10, pp. 695-711. DOI: 10.1038/s44222-023-00089-7.

Morgensztern, D., S. Waqar, J. Subramanian, F. Gao, K. Trinkaus, and R. Govindan (2012). "Prognostic Significance of Tumor Size in Patients with Stage III Non-Small-Cell Lung Cancer: A Surveillance, Epidemiology, and End Results (SEER) Survey from 1998 to 2003". In: *Journal of Thoracic Oncology* 7.10, pp. 1479-1484. ISSN: 1556-0864. DOI: 10.1097/jto.0b013e318267d032.

Nayak, L., E. Q. Lee, and P. Y. Wen (2012). "Epidemiology of Brain Metastases". In: *Current Oncology Reports* 14.1, pp. 48-54. ISSN: 1523-3790. DOI: 10.1007/s11912-011-0203-y.

Nicol, P. B., D. L. Barabási, K. R. Coombes, and A. Asiaee (2022). "SITH: An R package for visualizing and analyzing a spatial model of intratumor heterogeneity". In: *Computational and Systems Oncology* 2.2. ISSN: 2689-9655. DOI: 10.1002/cso2.1033.

Nicolò, C., C. Périer, M. Prague, C. Bellera, G. MacGrogan, O. Saut, and S. Benzekry (2020). "Machine Learning and Mechanistic Modeling for Prediction of Metastatic Relapse in Early-Stage Breast Cancer". In: *JCO Clinical Cancer Informatics* 4. PMID: 32213092, pp. 259-274. DOI: 10.1200/CCI.19.00133. eprint: https://ascopubs.org/doi/pdf/10.1200/CCI.19.00133. URL: https://ascopubs.org/doi/abs/10.1200/CCI.19.00133.

Norton, L. (1988). "A Gompertzian model of human breast cancer growth." In: *Cancer research* 48.24 Pt 1, pp. 7067-71. ISSN: 0008-5472.

Rami-Porta, R., V. Bolejack, J. Crowley, D. Ball, J. Kim, G. Lyons, T. Rice, K. Suzuki, C. F. Thomas, W. D. Travis, and Y.-L. Wu (2015). "The IASLC Lung Cancer Staging Project: Proposals for the Revisions of the T Descriptors in the Forthcoming Eighth Edition of the TNM Classification for Lung Cancer". In: *Journal of Thoracic Oncology* 10.7, pp. 990-1003. ISSN: 1556-0864. DOI: https://doi.org/10.1097/JTO.0000000000000559. URL: https://www.sciencedirect.com/science/article/pii/S1556086415334869.
Rittberg, R., S. Banerji, J. O. Kim, S. Rathod, and D. E. Dawe (2021). "Treatment and Prevention of Brain Metastases in Small Cell Lung Cancer". In: *American Journal of Clinical Oncology* 44.12, pp. 629-638. ISSN: 0277-3732. DOI: 10.1097/coc.0000000000000867.

Rowan, T. H. (1990). "Functional stability analysis of numerical algorithms". PhD thesis. University of Texas at Austin.

Schlicke, P., C. Kuttler, and C. Schumann (2021). "How mathematical modeling could contribute to the quantification of metastatic tumor burden under therapy: insights in immunotherapeutic treatment of non-small cell lung cancer". In: *Theoretical Biology and Medical Modelling* 18.1, p. 11. ISSN: 1742-4682. DOI: 10.1186/s12976-021-00142-1. URL: https://doi.org/10.1186/s12976-021-00142-1.

Strobl, M. A., J. West, Y. Viossat, M. Damaghi, M. Robertson-Tessi, J. S. Brown, R. A. Gatenby, P. K. Maini, and A. R. Anderson (2021). "Turnover Modulates the Need for a Cost of Resistance in Adaptive Therapy". In: *Cancer Research* 81.4, pp. 1135-1147. ISSN: 0008-5472. DOI: 10.1158/0008-5472.can-20-0806.

Zeng, H., L. E. Hendriks, W. H. v. Geffen, W. J. Witlox, D. B. Eekers, and D. K. D. Ruysscher (2020). "Risk factors for neurocognitive decline in lung cancer patients treated with prophylactic cranial irradiation: A systematic review". In: *Cancer Treatment Reviews* 88, p. 102025. ISSN: 0305-7372. DOI: 10.1016/j.ctrv.2020.102025.
Supplementary Figures

<figure><img src="image_14.png" /><figcaption>Supplementary Figure S1. Disease history of the patients in the BM-SCLC dataset</figcaption></figure>

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>**Supplementary Figure S3.** Impact of the Model Parameters. All simulations assume a primary tumor diameter of 10 mm at diagnosis. The red curve represents the growth of the primary tumor, while the blue curves correspond to the growth of individual metastases. **A.** Simulation using a Gompertz growth model with parameters $\alpha = 0.1$, $\mu = 0.001$, and $\gamma = 0.1$. **B.** Simulation using a Gompertz growth model with parameters $\alpha = 0.044$, $\mu = 5 \times 10^{-4}$, and $\gamma = 0.1$. **C.** Simulation using a Gompertz growth model with parameters $\alpha = 0.044$, $\mu = 0.001$, and $\gamma = 0.01$. **D.** Simulation using a Gompertz growth model with a 90% reduction in primary tumor diameter after diagnosis, with parameters $\alpha = 0.044$, $\mu = 0.001$, and $\gamma = 0.1$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Supplementary Figure S4.** Distribution of best-fitting *γ* values across the AP-HM patient population. The bar chart shows the number of patients for whom each value of *γ* (0.1, 0.3, 0.5, or 1) provided the best fit to their observed metastatic size distribution data.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Supplementary Figure S5.** Distribution of the predicted age of PT at the time of diagnosis in the AP-HM population</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Supplementary Figure S6.** Distribution of estimated individual *μ* parameter as a function of the estimated *γ* values in the AP-HM population.</figcaption></figure>
A

<figure><img src="image_20.png" /></figure>

B

<figure><img src="image_21.png" /></figure>

C

<figure><img src="image_22.png" /></figure>

D

<figure><img src="image_23.png" /></figure>

E

<figure><img src="image_24.png" /></figure>

F

<figure><img src="image_25.png" /></figure>

Supplementary Figure S7. Bootstrap distributions of model parameters incorporating PCI status as covariate.

A. Distribution of α parameter estimates. B. Distribution of μ parameter estimates. C. Distribution of ω₂_α parameter estimates. D. Distribution of ω₂_μ parameter estimates. E. Distribution of β_{PCI,α} coefficient estimates. F. Distribution of β_{PCI,μ} coefficient estimates. Vertical solid lines indicate maximum likelihood estimates, and vertical dashed lines indicate 95% confidence intervals for each parameter.
Supplementary Tables

<table><thead><tr><th>Variable</th><th>No PCI, N = 61<sup>1</sup></th><th>PCI, N = 129<sup>1</sup></th></tr></thead><tbody><tr><td>Number of events</td><td>26 (43%)</td><td>57 (44%)</td></tr><tr><td>TNM T</td><td></td><td></td></tr><tr><td>T0</td><td>1 (1.6%)</td><td>3 (2.3%)</td></tr><tr><td>T1</td><td>15 (25%)</td><td>15 (12%)</td></tr><tr><td>T1b</td><td>1 (1.6%)</td><td>1 (0.8%)</td></tr><tr><td>T2</td><td>10 (16%)</td><td>32 (25%)</td></tr><tr><td>T2a</td><td>1 (1.6%)</td><td>0 (0%)</td></tr><tr><td>T2b</td><td>1 (1.6%)</td><td>3 (2.3%)</td></tr><tr><td>T3</td><td>17 (28%)</td><td>34 (26%)</td></tr><tr><td>T4</td><td>15 (25%)</td><td>41 (32%)</td></tr></tbody></table>

¹n (%)

**Supplementary Table S1.** BMFS population characteristics