Estimating latent, dynamic processes of breast cancer tumour growth and distant metastatic spread from mammography screening data

We propose a framework for jointly modelling tumour size at diagnosis and time to distant metastatic spread, from diagnosis, based on latent dynamic sub-models of growth of the primary tumour and of distant metastatic detection. The framework also includes a sub-model for screening sensitivity as a function of latent tumour size. Our approach connects post-diagnosis events to the natural history of cancer and, once refined, may prove useful for evaluating new interventions, such as personalised screening regimes. We evaluate our model-fitting procedure using Monte Carlo simulation, showing that the estimation algorithm can retrieve the correct model parameters, that key patterns in the data can be captured by the model even with misspecification of some structural assumptions, and that, still, with enough data it should be possible to detect strong misspecifications. Furthermore, we fit our model to observational data from an extension of a case-control study of post-menopausal breast cancer in Sweden, providing model-based estimates of the probability of being free from detected distant metastasis as a function of tumour size, mode of detection (of the primary tumour), and screening history. For women with screen-detected cancer and two previous negative screens, the probabilities of being free from detected distant metastases 5 years after detection and removal of the primary tumour are 0.97, 0.89 and 0.59 for tumours of diameter 5, 15 and 35 mm, respectively. We also study the probability of having latent/dormant metastases at detection of the primary tumour, estimating that 33% of patients in our study had such metastases.


A.2 Data-Generating Mechanisms
We simulated independent datasets of 1500 individuals each, assuming tumours originate from a single spherical cell of volume cell corresponding to a diameter of Cell = 0.01 mm; we also assumed that tumours are detectable only when when they reach a given volume 0 corresponding to a diameter 0 = 0.5 mm. We assumed the parameter to be fixed at = 4.
We simulated inverse growth rates from a GammaΓ distribution with equal shape and rate parameters 1 = 2 = (for identifiability purposes), and we fixed = exp(0.2938482). This corresponds to a GammaΓ distribution with a mean of 1 and a variance of approximately 0.75).
Then, we simulated the volume of the tumour at symptomatic detection at the time from inception to detection from the model for the rate of symptomatic detection and the density for tumour volume at symptomatic detection; we fixed the parameter ′ = − log( ) = 9.3811683.
The number of metastases being seeded per each simulated individual was drawn from a Poisson distribution with subject-specific intensity depending on = det (simulated at the previous step) and assuming = exp (−15.5223994).
The timing of the first seeding was simulated using the inversion method as described in Bender et al. (2005).
The true values of the parameters described above were selected by fitting the model to a subset of data from CAHRES consisting of women that were detected symptomatically and with no history of previous negative screens. CAHRES is described in more detail in the body of the manuscript and elsewhere (Eriksson et al. 2013;Rosenberg et al. 2006Rosenberg et al. , 2008.
We therefore calculated time to metastasis met as the sum where 1s is the time to the first seeding, 0 is the time it takes for each tumour to growth from volume Cell to volume 0 (conditional on each subject-specific inverse growth rate ), and det is the time to symptomatic detection (simulated above).
Finally, we applied administrative censoring after 20 years of follow-up.

A.3 Estimands
The estimands of interest are the parameters of the model for tumour growth and metastatic spread in the absence of screening: = { , , }. In practice, we modelled and on the log-scale and ′ = − log( ).

A.4 Methods
We fitted the true data-generating model in the absence of screening (as described in the body of the manuscript) to each independent simulated dataset. Standard errors of the model parameters were computed by inverting the Hessian matrix at the optimum, which was calculated using Richardson extrapolation as implemented in the R package numDeriv (Fornberg and Sloan 1994;Gilbert and Varadhan 2019;Lindfield and Penny 1989).

A.5 Performance Measures and Number of Replications
Given the aim of this simulation study, the performance measures of interest are: • Bias, to quantify whether the estimation procedure can retrieve the true data-generating parameters on average. This is the key performance measure. We also report average and median estimated values of each estimand; • Empirical and model-based standard errors, to quantify the precision of the point estimates; • Coverage probability, quantifying the probability that a confidence interval contains the true value of each estimand.
Further details on each performance measure are described elsewhere (Morris et al. 2019); we report Monte Carlo standard errors of each performance measure as well, quantifying the uncertainty in their estimation.
To estimate the key performance measure with a given precision, we first ran 10 replications of this simulation study. Allowing a Monte Carlo standard error for bias of 0.01, we would require sim = Var/MCSE 2 ≈ 31, where Var = 0.0031 is estimated from the largest value among all empirical and model-based standard errors forô btained from the above-mentioned 10 replications. To be conservative we rounded up that value to the nearest 100 (with a minimum of 200 replications), yielding a total of sim = 200.

A.6 Software
This simulation study was coded and run using R version 4.0.2 and analysed using the rsimsum package (Gasparini 2018; R Core Team 2020). The likelihood function was optimised using the optim function in R and the limitedmemory modification of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method (Byrd et al. 1995).

A.7 Results
Descriptive characteristics of the metastatic process across all simulated datasets are presented in Table S1. On average, 5.07% of individuals had detectable metastases at the time of primary tumour detection, 64.04% developed a detectable metastasis during follow-up, and 30.89% were event-free at the end of follow-up. The proportion of individuals developing detectable metastases during follow-up was rather high, but we did not simulate any competing event (such as death) here.
All models converged to a solution, with no replications showing issues with the convergence of the optimisation algorithm or calculation of the Hessian matrix. The comparison of model-based and empirical standard errors showed that our procedure (inverting the numerical Hessian at the optimum) slightly overestimated the standard errors for log( ) and log( ), with a relative percentage error of 7.93% and 9.16%, respectively; model-based standard errors for − log( ) were underestimated, with a relative percentage error of -1.69%.

B Supplementary Figure S1
Supplementary Figure

D.1 Aim
This simulation study aims to assess the robustness of the proposed model to alternative data-generating mechanisms. We focus on the model in the absence of screening, for (1) simplicity and (2) to keep the focus on the core of the new model introduced in this paper.

D.2 Data-Generating Mechanisms
We simulate data under five different data-generating mechanisms. Specifically, we define three scenarios where we relax the assumption of common growth rates between the primary tumour and seeded metastases (scenarios 1-3), where we relax the assumption of Gamma-distributed growth rates (scenario 4), and where we allow for heterogeneity in the spread parameter (scenario 5).
In scenarios 1-3, we simulate growth rates for the primary tumour and each seeded metastasis from a bivariate distribution that we constructed using a bivariate Gaussian copula (Meyer 2013) with correlation parameter values = 1.00, 0.90, and 0.75, respectively. The margins of the two growth rate distributions are assumed to be Gamma-distributed, with equal shape and rate parameters 1 = 2 = 2.3. Note that scenario 1 corresponds to perfectly correlated growth rates, and therefore represents the model introduced in the main body of the paper; it was however included as a check on the data generation procedure based on the copula model.
In scenario 4, we relax the assumption that growth rates follow a Gamma distribution. Specifically, we simulate growth rates from a log-Normal distribution with mean and standard deviation (on the log-scale) of −0.25 and 0.75, respectively.
In scenario 5 we allow for heterogeneity in the spread parameter . Specifically, we draw subject-specific values of the spread parameter from a Gamma distribution with scale and rate parameters of 0 and 16 , respectively.
For all scenarios, we simulate data assuming ′ = 9.0; the remaining parameter values and assumptions are the same as those used in the simulation study of Appendix A.

D.3 Estimands
We consider four model-based predictions under each data-generating scenario:

D.4 Performance Measures
Predictions under each simulation scenario are presented/evaluated graphically. Specifically: 1. We compare the distribution of predicted median tumour doubling times using box plots; 2. We compare the probabilities of detected and latent metastasis at diagnosis (i.e. estimands 2 and 3, above) with empirical probabilities calculated from a large simulated dataset (with 100,000 subjects) under each data-generating mechanism. Empirical probabilities are computed across categories of tumour sizes, defined as 5 mm wide intervals. Exact confidence intervals for the empirical probabilities are additionally included; 3. We compare the predicted survival probabilities to survival curves obtained using the Kaplan-Meier estimator on the same large simulated datasets used in point 2, above, with size categories defined using ±2.5 mm from the central point; 4. We compare the fitted distribution of tumour volume at detection with a histogram of tumour volumes on raw data; 5. Finally, to provide a comprehensive view of the results, we also include box plots for the distribution of fitted model parameters across scenarios and repetitions.
Furthermore, as an additional comparison, we include another scenario where we simulate data from the correct model (as in Appendix A) but using the parameters described in Appendix D.2 instead. This is referred to below as the Reference.
We run 200 repetitions per simulation scenario (as in Appendix A).

D.5 Software
This simulation study was coded and run using R version 4.0.4; once again, the model likelihood was optimised using the optim function and the limited-memory modification of the BFGS algorithm (Byrd et al. 1995; R Core Team 2020). The distribution of median tumour doubling times across scenarios is depicted in Figure S3. The distributions largely overlap, with increasingly larger values in scenarios 2 and 3 as the correlation between growth rates of the primary tumour and metastases diminishes. Conversely, in the scenario with heterogeneity in the spread parameter, the median of the distribution was slightly lower. Nevertheless, differences between all scenarios are still relatively small (with a range of approximately 10 days between the largest and smallest median of each distribution).  where growth rates of the primary tumour and distant metastases were not strongly correlated. Figure S5 depicts the fitted probability of having latent, yet undetected metastases at diagnosis as a function of tumour size; once again, empirical probabilities are included for comparison. All scenarios showed a decent fit between predicted and empirical probabilities, with the model only slightly overestimating the probability of having latent metastases.  Finally, in Figure S7 we compare the parameter estimates across the simulation scenarios. Generating data based on non-perfectly-correlated growth rates between the primary tumour and metastases (scenarios 2 and 3) resulted in log( 1 ) and − log( ) being (on average) overestimated, with there being increasingly larger differences as the correlation diminished. This corresponds to overestimating the median inverse growth rate (which is in accordance with Figure S3) and underestimating the rate of symptomatic detection. Conversely, log( ) was estimated with little bias, across all scenarios, except for scenario 5 (data-generating mechanism with heterogeneity in ) where log( ) was slightly overestimated. Nevertheless, the bias was marginal and (as seen in the previous plots) did not substantially affect model-based predictions.