Prediction of Microscopic Metastases in Patients with Metachronous Oligo-Metastases after Curative Treatment of Non-Small Cell Lung Cancer: A Microsimulation Study

Simple Summary Many patients with metachronous oligo-metastases in non-small cell lung cancer have their recurrences surgically removed, although the 5-year recurrence-free survival of this group is 16%. This does not provide any benefit for patients with additional undetected metastases. Therefore, we aim to find patient characteristics that are predictive for having additional undetected microscopic metastases. Based on a theoretical approach, we identified the size and number of detected oligo-metastases, as well as the presence of symptoms that are the most important risk predictors. Abstract Metachronous oligo-metastatic disease is variably defined as one to five metastases detected after a disease-free interval and treatment of the primary tumour with curative intent. Oligo-metastases in non-small cell lung cancer (NSCLC) are often treated with curative intent. However additional metastases are often detected later in time, and the 5-year survival is low. Burdensome surgical treatment in patients with undetected metastases may be avoided if patients with a high versus low risk of undetected metastases can be separated. Because there is no clinical data on undetected metastases available, a microsimulation model of the development and detection of metastases in 100,000 hypothetical stage I NSCLC patients with a controlled primary tumour was constructed. The model uses data from the literature as well as patient-level data. Calibration was used for the unobservable model parameters. Metastases can be detected by a scheduled scan, or an unplanned scan when the patient develops symptoms. The observable information at time of detection is used to identify subgroups of patients with a different risk of undetectable metastases. We identified the size and number of detected oligo-metastases, as well as the presence of symptoms that are the most important risk predictors. Based on these predictors, patients could be divided into a low-risk and a high-risk group, having a model-based predicted probability of 8.1% and 89.3% to have undetected metastases, respectively. Currently, the model is based on a synthesis of the literature data and individual patient-level data that were not collected for the purpose of this study. Optimization and validation of the model is necessary to allow clinical usability. We describe the type of data that needs to be collected to update our model, as well as the design of such a validation study.


Introduction
Metachronous oligo-metastatic disease is variably defined as one to three or one to five metastases detected after a progression-free interval and a controlled primary tumour [1][2][3][4]. Oligo-metastases in non-small cell lung cancer (NSCLC) are often treated with curative intent, but additional metastases are commonly detected later in time, and the 5-year survival is low in this group of patients [5]. These patients are assumed to have had poly-metastatic disease at time of the oligo-recurrence, with most metastases below the detection threshold of the scan. In this paper, we will refer to oligo-metastatic disease with additional undetected metastases as "oligo+", and without undetected metastases as "oligo−". Currently, clinicians have no tools available that can distinguish "oligo+" from "oligo−" patients [6][7][8].
Curative therapy is likely to provide limited benefit in oligo+ patients [9]. Furthermore, such therapies can result in adverse events or even treatment-related mortality [10,11]. Therefore, it is vital to develop tools to distinguish between oligo+ and oligo− patients; that is, identify the predictors that would allow for such a distinction.
There are two important obstacles for the development of such a tool. Firstly, the true outcome of interest is the presence of unobserved metastases, an outcome that is by definition unknown (Box 1). Therefore, most studies use surrogate endpoints, such as prognostic factors, for survival [12][13][14]. These studies include factors such as the patients' age and performance score into their models, which may select a subgroup of both curable and incurable patients that live relatively longer, or administration of chemotherapy. Those factors are good predictors when the focus is on prognosis in a mixed group of patients with oligo+ and oligo−.
A second obstacle is that there are difficulties with patient accrual in randomized controlled trials [15]. Currently, most studies on oligo-metastases are retrospective studies, with small groups of patients. As a result, these studies are susceptible to selection bias, and the levels of evidence are often weak [16][17][18][19].
To tackle both obstacles, we have chosen to construct a microsimulation model of the development and detection of metastases in stage I NSCLC patients. Microsimulation models allow simulation of detailed disease trajectories of a large number of individual patients based on a mathematical model describing how patients transition between health states. Microsimulation models have the advantage that all available evidence can be synthesized to calibrate unobservable parameters, such as the presence of undetectable metastases. In our model, curatively treated patients may develop a number of metastatic lesions that grow exponentially over time. The microsimulation model keeps track of the growth all metastases within one patient, the number of metastases detected either by surveillance or by symptoms, and the recurrence-free survival.
We used a microsimulation model to construct a simulated patient-level dataset, in which individuals were classified as having detected oligo-metastases or poly-metastases, with the former being subdivided into an oligo− and oligo+ group. Subsequently, we used the simulated dataset to identify the clinically observable patient characteristics that can predict the presence of undetected metastases.
We describe the development of the microsimulation model and the underlying evidence. In addition, we present the simulated dataset and the identification of the observable patient characteristics that are related to the risk of having undetected recurrences. Finally, we discuss the requirements for future validation of our findings and identification of additional predictors, such that the identification of patients that benefit from curative treatment of oligo-metastases may be improved.

Box 1. Definitions used.
Curative treatment: Treatment focused on removing or destroying tumours locally, with the aim to make the patient disease free. Metastases: All tumours originating from the primary tumour or another metastasis. Recurrences: All metastases detected on a scan. Oligo-metastasis: Several definitions have been used in the literature. In this paper, oligo-metastasis is defined as three or less detected recurrences within one patient, with a controlled primary tumour. Poly-metastasis: Several definitions have been used in the literature. In this paper, poly-metastasis is defined as four or more metastases detected within one patient. Metachronous oligo-metastasis: The oligo-metastasis is detected after curative treatment of the primary tumour, contrary to synchronous oligo-metastasis. Undetectable metastases: In this paper, all (microscopic) metastases of a size below the detection threshold. RFS: Recurrence-free survival-the proportion of patients that remain recurrence-free after curative treatment of the primary tumour. PFS: Progression-free survival-the proportion of patients that have no additional recurrences detected after curative treatment of oligo-metastases. oligo+: A patient with oligo-metastasis(es) with additional undetectable metastases, who would be classified as poly-metastatic if the true number of metastases would be known. oligo−: A patient with oligo-metastasis(es) without additional undetectable metastases. VDT: Volume doubling time-the time required for the tumour to double in volume.

Concept of the Microsimulation Model
Curative treatment of oligo-metastases has a pooled 16% 5-year progression-free survival (PFS) [6][7][8]20,21]. This low PFS is assumed to be caused by undetected metastases that already existed at time of treatment of the primary tumour ( Figure 1). The purpose of the microsimulation model is to test which observable patient characteristics during detection of oligo-metastases may be predictive for the presence of undetected metastases, by synthesizing the available evidence on the growth of tumours, recurrence patterns, and methods of detection using the current theory and available data. Metastases are considered undetectable when their size is below the detection threshold (dashed lines), and thus become detectable above the detection threshold (solid lines). If a patient is scanned at time "a", all metastases are invisible, and this patient would be considered to be "recurrence-free". At time "b", three recurrences would be visible on the scan. This is defined as an oligo-metastasis, even though there are several metastases under the detection threshold (oligo+). At time "c", 8 recurrences are visible on the scan, which is defined as poly-metastatic disease.
The microsimulation model was developed in C++ and describes the growth and detection of metastases in individual patients. The model stores all patient-specific features and outcomes in comma-separated values files for further analysis. Microsoft Excel Professional plus 2016 was used to create the figures. Certain assumptions and simplifications had to be made, to allow modelling of cancer progression and metastases, without making the model overly complex for the research question that we aim to address (Box 2). These model assumptions have been tested in the sensitivity analyses. The model contains both directly observable parameters as well as unobservable parameters. The observable parameters were directly estimated based on patient data and the literature (Table 1). Calibration against the observable output variables was used to estimate the unobservable (hidden) parameters, such as the number of undetectable metastases (Table 1).
The most important assumptions of the microsimulation model are listed below, as well as the reasoning behind these assumptions. 1.
An exponential tumour growth model using volume doubling time was chosen, because volume doubling time is the most commonly used statistic in the literature. All metastases within one patient are assumed to have the same growth rate, and are formed consecutively with fixed time intervals. This assumption has been tested in Scenario 2 of the sensitivity analyses.

2.
Each patient is assumed to have a fixed number of metastases after their primary tumour had been curatively treated. The proportion of patients that have zero metastases after curative treatment equals the proportion of patients that are recurrence-free after 5 years. In all other patients, the number of metastases is randomly drawn from a rounded truncated normal distribution.

3.
We assume that metastases below the minimum detectable size for the Computed Tomography scan will always be missed.

4.
Metastases of detectable size are either found during surveillance or on an unscheduled scan because of symptoms, whichever happens first. Other scenarios are not considered. Surveillance CT scans are able to detect recurrences to the lung, liver, and adrenal glands. Bone and brain metastases are highly symptomatic. Less than 3% of NSCLC metastases are found in other organs, and these metastases are often also symptomatic [22,23]. Therefore, we assume that the com-bination of symptomatic and surveillance detection sufficiently describes the detection patterns.

5.
Once one recurrence is detected, a more rigorous examination, that is, Positron Emission Tomography-Computed Tomography, follows, resulting in detection of all metastases above the minimum detectable size. 6.
The proportion of patients with microscopic metastases within those with detected oligorecurrent disease is assumed to be equal to the proportion of patients with a 5-year PFS after treatment of oligo-recurrent disease. This may lead to an underestimation of the proportion of oligo-metastases. Therefore, this assumption was further investigated in Scenario 1 of the sensitivity analyses.

Model Functions
Each simulated patient starts with a number of metastases (M total ) drawn from a truncated normal distribution, N(σ,µ) ≥ 1, rounded up to integer numbers. All metastases in the model grow with a patient-specific volume doubling time (VDT) in days −1 [35].
Here, V m (t) is the volume (in cm 3 ) of metastasis m, m = 1 . . . M total , at time t (in days). The size of the largest metastasis is used to calculate the sizes of the other metastases. The relative sizes of the metastases are defined by the size ratio R [36]: In the model, metastases can either be detected during a routine surveillance scan, or when the metastases become symptomatic. Metastases become detectable on a CT scan once their size is equal or larger than the minimum detectable volume, V det . An exponential hazard function for becoming symptomatic is assumed: with λ symptom the hazard to develop symptoms. Time t det m is defined such that V m (t detm ) = V det . Metastases smaller than V det are assumed to be too small to cause symptoms.

Parameter Estimation
The model parameters have been directly estimated from patient data and the literature when possible. All the unobservable parameters were estimated using calibration in combination with other observable target parameters.
Data of stage I NSCLC patients curatively treated with video-assisted thoracoscopic surgery and stereotactic body radiation therapy were obtained from two studies performed in the Netherlands between 2003 and 2013 [24,25]. Patients were excluded if they had ≥stage II disease, an ECOG performance score ≥2, a second primary tumour, or history with previous cancer. Patients of both studies were pooled and their 1:1 propensity score matched on treatment using the Matching R package, version 4.9-2 [37]. Recurrence-free survival (RFS) was analysed with a Kaplan-Meier survival curve using Statistical Package for Social Sciences version 22.0 [38].
The VDTs were obtained from papers that reported the distribution of the growth rates measured with CT equipment, in NSCLC with mixed histologies and a non-specific tumour morphology [26][27][28][29]. Some papers classified growth rates into groups (fast, moderate, slow, and no growth). For these groups, the average of the minimum and maximum growth rates was assumed to be the average VDT of the reported group. A cumulative exponential distribution was fitted to the reported data using a least squared estimate function ( Figure 2). This distribution was used to sample the VDTs for patients in the microsimulation model.
Smaller metastases are more likely to be missed on a scan [39,40]. As such, the reported scan sensitivities and specificities for tumours of any size as found in the literature are not suitable for our continuous tumour growth microsimulation model [41]. Instead, we have therefore chosen to use a detection threshold. We assumed a detection threshold of 5 mm diameter on the CT scans [31,32], corresponding to a spherical volume of 0.07 cm 3 .  [26][27][28][29]. A cumulative distribution function is fitted using a least square function to the data shown in this scatterplot to obtain a quantile function, which was subsequently used to determine the random tumour growth rates per patient.

Model Calibration
Four (sets of) parameters needed to be calibrated: (1) the detection rate parameter λ det , which is used to randomly draw the time when the largest metastasis passes the detection threshold; (2) the parameters used to determine the total number of metastases per patient M total (µ,σ); (3) the parameters used to determine the size ratio of the metastases R(α,β); and (4) the hazard of a metastasis becoming symptomatic (λ symptom ). Each of these model parameters is strongly associated with a specific model output parameter that is also available from the literature or patient-level data (see below). The mean squared error of 1000 simulations was used in combination with a univariate or grid search algorithm (for M total (µ,σ), and R(α,β)) and repeated until the calibration target was matched up to 3 significant figure places [42,43].
Calibrated parameters only interact in one direction: the number of metastases influences the chance that a patient becomes symptomatic, but the hazard of a single metastasis becoming symptomatic does not influence the number of metastases per patient. Therefore, parameters could be calibrated consecutively if done in the correct order-the most dominant parameter first.
Firstly, the distribution of the time to detectability (that is, reaching a 5 mm diameter) of the largest metastasis was assumed to be a negative exponential distribution . λ det was calibrated against the 1 2 yearly RFS, weighted for the number of patients at risk. Patient-level data of the two retrospective Dutch cohort studies described above were pooled to obtain the RFS statistics [24,25]. Both the simulations and the patient-level data used the same surveillance schedule of the Dutch guidelines [44].
Secondly, M total represents the number of both detected and undetectable metastases per patient. There are no data available on this number, nor on the distribution of the detected metastases. Therefore, M total is assumed to be normally distributed, but is truncated such that M total ≥ 1. M total (µ,σ) was calibrated against the proportion of oligo-metastases without undetected metastases. To determine this proportion, it was assumed that patients who are curatively treated for oligo-metastases and are progression-free at 5 years do not have undetected metastases. A weighted average of the estimates reported in the medical scientific literature was calculated, resulting in an estimated 16% of oligo-metastases that have no underlying undetected metastases [6][7][8]20,21]. To allow a unique solution for µ and σ, we added to following constraint: the solution for µ and σ that minimizes the maximum number of detected metastases under the current surveillance schedule.
Thirdly, the size ratio of the metastases, R, is calculated as the volume of the second largest metastasis divided by the volume of the largest metastasis. Therefore, R has a value between 0 and 1 by definition. We used a Beta distribution to draw the R for each patient. The shape parameters α,β were defined as integer values. R(α,β) was calibrated against the pooled average proportion of oligo− [7,20,21]. As a constraint, the solution with the smallest lower tail was selected.
The last calibrated variable was λ symptom , the hazard of a single metastasis becoming symptomatic. λ symptom was calibrated against the ratio of the symptomatically (or unscheduled) detected metastases to metastases detected by a scheduled scan [33,34].

Model Simulations
The life histories of 100,000 patients with stage I NSCLC and one or more undetected metastases were simulated from treatment of the primary tumour until the detection of a metastasis (RFS). Patients without micro-metastatic disease after curative treatment of the primary tumour were not simulated but were added to the cohort post-simulation in accordance with the proportion of patients that were recurrence-free after 5 years.
For the simulation, each individual was randomly generated from the fitted distributions: the total number of metastases M total , the size ratio R of the metastases, the volume doubling time, the time that the metastases would grow past the detection threshold, t det , and the time that the patient would develop symptoms, t symp . Subsequently, the model determines if symptomatic detection occurs before or after the time that the tumour would be detected on a surveillance scan. Once the time of detection is determined, the number of metastases visible on a scan m det can be calculated with Function (4). Function (4) can be derived from Functions (1) and (2) (Appendix A.1). The point in time that the largest metastasis becomes detectable on a scan is denoted by t det = 0.

Prognostic Groups
A binary logistic regression model was used for the purpose of identifying the clinically observable predictors that can distinguish oligo− from oligo+. The parameters that were considered clinically observable were the number of recurrences detected; the size category of the largest detected recurrence: small (<6 mm), medium (6-8 mm), or large (>8 mm) [32], symptomatic or surveillance detection; and the RFI (years) [12]. All predictive covariates were tested for multi-collinearity by calculating the variance inflation factor (VIF) for each variable. A residual and a binned residual plot for the logistic regression analysis was made to test if the assumptions of the logistic regression model were met.
Predictors with a Wald statistic p < 0.05 were used to divide the simulated patients into prognostic sub-groups and determine the proportion of oligo+ for each of these groups. These sub-groups were pooled into a low-risk group with less than 30% oligo+, and a high-risk group with more than 30% oligo+. The 30% risk threshold was considered an acceptable level to consider curative therapy, based on a discussion with our clinical experts (LA, EAK, SYS, and FMNHS). Still, this 30% threshold is suggested tentatively, as it is not obtained from any formal consensus procedure. However, the threshold allows for accessible presentation of the performance and potential use of our modelling approach, and can be seen as a first step towards decision support in this area.

Sensitivity Analyses
To test the robustness of our results, we have designed a series of scenarios and tested how much the predictors used to determine the risk groups, proportion of oligo+, and size of the risk groups are affected by these scenarios.

1.
Recalibration of model parameters to the upper and lower confidence interval of their targets.

2.
Random variation in VDT of metastases within one patient and in the detection threshold.

3.
Correlation between the volume doubling time and the total number of metastases per patient.
The ability of the metastases to produce new metastases.
All sensitivity analysis scenarios are described in detail in Appendix A.

Calibration
An overview of all the calibrated parameters can be found in Table 1. Figure 3 shows the resulting model-based progression-free survival curve compared to the two Dutch cohorts [24,25]. The simulated progression-free survival is within the 95% CI range of the cohort data with the exception of the first and third month, but the 95% CI was never exceeded by more than 1.1% after calibration of λ detectable . Recurrence-free survival after curative treatment of the primary tumour of the simulated patients is calibrated to match two Dutch cohorts [24,25]. Recurrences are detected with a surveillance scan schedule and unplanned additional scans using a symptom hazard model. A total of 62.5% of the Dutch cohorts remained recurrence-free. Patients without recurrences were not simulated but were added after the simulation for the purpose of showing the differences between the simulated patients and the Dutch cohorts in this figure. Number of detected recurrences per month. Asymptomatic recurrences were only detected during surveillance scans. These scans were planned according to a fixed schedule, with a small variation in scan date N (µ = FU schedule, σ = 15.5). Symptomatic recurrences were detected throughout the year, and detection times were determined with the hazard model (Function (3)).
Calibration of M total resulted in a maximum of 70, an average of 9, and a median of 6 metastases detected. This solution was considered realistic, based on discussion with our clinical experts LA, EAK, SYS, and FMNHS.
Calibration of the size ratio R resulted in α = 9 and β = 1, giving an average R of 0.93. The proportion of symptomatically detected metastases was 33.5% after calibration of the parameter λ symptom , similar to the proportion observed in clinical practice [33,34].

Simulation Results
The microsimulation model was used to produce a patient-level database of 100,000 patients. The first row of Table 2 shows the characteristics for all the simulated patients. These patients are subdivided into three subgroups in row two to four: poly-metastases (patients with 4 or more detected metastases), oligo+ (patients with 3 or less detected metastases and additional undetectable metastases), and oligo− (patients with 3 or less detected metastases and no additional undetectable metastases). Between these three subgroups, the ranges for all the reported characteristics are largely overlapping, although the skewness of the distributions is different. The median oligo+ patient has a large size ratio between the metastases, slow-growing metastases, small diameters at the time of detection, and a low RFI. Table 2. Description of the generated patient-level data (the median, with the 2.5% and 97.5% percentiles in brackets).

Prognostic Groups
The VIF for the predictors in the logistic regression model ranged from 1.016 to 1.097, and no transformation of variables were assumed to be needed based on the residual and binned residual plots ( Figure A1). Although there is a large difference in the median RFI of the oligo+ and oligo− groups, as shown in Table 2, RFI is not a significant predictor in the logistic regression model (Table 3), when the size of the largest metastasis is included in the model. All the other predictors (size, number of detected metastases, and symptomatic detection) are significant, and are therefore used to construct Table 4. A logistic regression model was used as an explorative search for the parameters suitable to determine the risk groups. Odds Ratios are calculated as exp(β). Year of detection was removed from the regression model because it was not a significant predictor (Wald test p > 0.05). All other covariates were significant (Wald test p < 0.0001). Note that the significance is influenced by the large number of simulated patients.  Table 4 shows the proportion of oligo+ at time of detection of recurrence within each of the subgroups defined by the three identified predictors. The 30% oligo+ threshold resulted in two distinct groups of patients. Patients with a large-sized metastasis or a single medium-sized metastasis fall in the low-risk group, and patients with small or 2 or 3 medium-sized metastases fall in the high-risk group. Although the symptoms do have an effect on the proportion of oligo+, this effect is not large enough to differentiate between low and high risk.
The pooled low-risk group has a proportion of 8.1% oligo+, and the high-risk group has a proportion of 89.3% oligo+. As shown in Table 2, the differences between the highand the low-risk groups are larger than the differences between oligo+ and oligo−. This is a result of the selection procedure for the risk groups. The high-and low-risk groups were compared to the "treat-all" and "treat-none" strategies on their strategy performance (Table 5). In total, 86% of the patients with an oligo-recurrence are in the oligo+ group, which is the cause of the difference in accuracy between "treat-all" and "treat-none". Table 5. Performance of the risk-group-based treatment selection compared to "treat-all" and "treat-none" strategies.

Strategy
Low  Additional analysis of the simulation output by extrapolation of the growth model of metastases showed that 73.8% of the oligo+ patients would switch to detectable polyrecurrent disease within 3 months, and 98.2% would become detectable within one year ( Figure A2).

Sensitivity Analyses
A summary of all the sensitivity analyses is given in Table 6. Details of the proportion of oligo+ per subgroup are shown in Tables A1-A17. The model output, such as the diameters of the metastases, or RFI, remained within realistic ranges under all scenarios. Table 6. Proportion of detected oligo-metastases that fall into the risk groups and proportion of oligo+ within the risk groups per scenario.

Sensitivity Analysis Scenario
All The classification into two risk groups and the proportion of oligo+ in the low-versus high-risk group were similar to our initial analysis in those scenarios in which the model parameters were recalibrated (Scenarios 1-8). A change in the proportion of oligo+ in all simulated patients was associated with a change in the proportion of oligo+ in the risk groups. This was the most extreme for scenarios in which the M total and R varied ( Figure A1).
Random variation in the detection size and growth rates of the metastases within one patient only caused a small increase in oligo+ patients in the low-risk group (Scenarios 9 and 10).
In Scenarios 11 and 12, in which the model parameters were correlated, the proportion of oligo+ in the low-risk group decreased even though the size of the low-risk group increased. Simultaneously, the proportion of oligo+ in the high-risk group increased. Thus, the ability to separate low-from high-risk patients increased with the correlated model parameters.
When oligo metastases are defined as 1 metastasis or as 1 to 5 metastases, this change in definition directly affects the total number of patients considered for curative treatment of their oligo-metastases. Furthermore, the number of detected metastases is a predictor that determines if a patient is high risk or not (Tables A13 and A14). Using the definitions of oligo-metastasis as a single detected metastasis results in a risk model with overall better accuracy, but also a smaller number of patients that is considered to have an oligo-recurrence.
Metastatic metastases had a large direct effect on the proportion of oligo+, resulting in a significant increase in the proportion of oligo+ in all subgroups. As such, more subgroups passed the 30% oligo+ threshold, resulting in fewer subgroups to be pooled into the lowrisk group. There are no patients with an acceptable risk to undergo curative therapy left, if 80% of all undetectable metastases are metastasizing at the time of curative treatment. Because a significant amount of patients do survive 5 years without additional recurrences after curative treatment of oligo-metastases in clinical practice, it is likely that the number of metastasizing metastases is low [6][7][8]20,21].
Over the different sensitivity analyses, the proportion of oligo+ at the time of detection of an oligo-metastasis remained relatively stable and resulted in acceptable average risks for curative therapy in the low-risk group. Furthermore, the combinations of predictors that were pooled into the low-risk group remained stable with the exception of the metastatic metastases scenarios.

Key Findings
A microsimulation model of the growth and detection of metastases in individual patients was built based on what we considered to be realistic assumptions and an appropriate level of complexity. The model is able to reproduce detection patterns of metastases as seen in the clinical setting, based on the growth rates and other model parameters obtained from the literature. The model output was used to identify the characteristics that allow prediction of the presence of microscopic disease; that is, the presence of additional undetectable metastases in patients with oligo-metastases. Based on three selected risk predictors, the population of patients with an oligo-recurrence could be grouped into a small low-risk group with an 8.1% risk and a high-risk group with an 89.3% risk of microscopic disease. Sensitivity analyses showed that the model is robust to changes in parameters.

Clinical Implications
Currently, many stage I NSCLC patients in whom oligo-metastatic disease is detected are curatively treated for these oligo-metastases, although only 16% of those patients remain progression-free for 5 years [6][7][8]20,21]. For these patients, curative treatment provides no benefit but does provide additional risks and side-effects. Multiple prognostic factors for patients with oligo-metastases have been identified, but most of these predictors, such as usage of adjuvant chemotherapy, are not suitable to predict the presence of undetected metastases and do not allow for selecting those patients that potentially benefit from curative treatment [12]. The key novelty in this paper is that we study the underlying tumour growth to predict which patients are most likely to have additional undetected metastases, rather than identifying the prognostic factors for survival.
Current guidelines suggest that curative treatment for oligo-metastases can be considered [44]. As a result, there is considerable variation in treatment of oligo-metastases in clinical practice. As an alternative to a "treat-all" or a "treat-none" strategy, patients can be classified into high-risk and low-risk groups on the basis of prognostic-factors for additional undetected metastases. Low-risk patients could be curatively treated, while high-risk groups could be treated as patients with poly-recurrences. Validation of our model is essential before using our model for clinical decision making. Such validation effort seems worthwhile, as our results show that the use of prognostic subgroups for clinical decision making could potentially reduce incorrect assignment of curative or systemic treatment for patients without or with undetected metastases, compared to the "treat-all" or "treat-none" strategies ( Table 5). The use of risk groups is likely to be both cost saving and beneficial for patients, as it may reduce unnecessary harmful treatments. A full cost-benefit analysis, including life expectancy, quality of life, and costs of the treatments, would be needed to confirm this assumption.
There are some patient subgroups in the model that have a relatively high chance to be assigned the wrong treatment. Ninety percent of the simulated patients had a very high or low risk of microscopic disease (above 80% or below 10%). A potential management strategy in the other ten percent could be to offer these patients more frequent surveillance, since the model predicted that 74% of patients with additional microscopic disease would be diagnosed as poly-recurrence within three months ( Figure A2).
One feature not added to the model is death due to other causes than metastatic disease. Inclusion of death due to another cause would have reduced the benefit of curative treatment in oligo-recurrent disease. Other factors that should be taken into consideration for treatment decision making are the eligibility criteria of a patient for a certain therapy and the effects of a therapy on symptoms. Therefore, multidisciplinary teams and shared decision making are important to provide the best treatment for a specific patient. This simulation model is meant as a tool to support such a decision-making process.
The microsimulation model and prognostic groups are not based on assumptions nor data of second primary lung cancers, synchronous oligo-metastases, or other types of cancer than NSCLC. Predictions of risk of undetectable metastases are therefore not necessarily applicable for these patients, although it is possible to collect data on the possible predictors, such as the sizes of the metastases in the clinical setting, and test if it is possible to create risk groups for those patients in a similar fashion. It can be informative and of great value for patients to collect clinical data to see if these conclusions are valid in those scenarios. Finally, once our prediction model for NSCLC is validated, the criteria of Martini and Melamed [45] should be applied to distinguish patients with potential oligo recurrences from patients with second primary lung cancers.

Microsimulation Model
Simulation modelling is a useful tool to synthesize the available evidence and extrapolate beyond currently observed data, especially in circumstances when relevant information and processes are either not available from a single source or unobservable. With respect to the growth of micro-metastatic lesions, there is a paucity in information that could support medical decision making. However, evidence from cell cultures [46], mouse models [47], and observed growth rates in larger tumours [26] provide indications about the tumour growth of micro-metastases in real patients. This information can be combined in a simulation model to allow inference on the disease aspects important for clinical or policy decisions. Although this microsimulation model is useful for improving insight into the metastatic process, and the selection of patients that may benefit from curative therapy, we need to keep in mind that all models are based on various assumptions. We have made an effort to be fully transparent by making all important model assumptions explicit and presenting the reasoning behind many of the choices made. Some model assumptions could be tested in the sensitivity analyses. However, validation of the model outcomes in clinical practice is still needed to determine if the model and its assumptions are correct.

Prognostic Groups
Prognostic factors have been identified and used to pool patients into risk groups. The 30% threshold value for low and high risk is arbitrary. Other thresholds may be used, and will influence the proportion of undetected metastases in the high-and low-risk groups. This threshold should therefore be in balance with the expected effects on quality of life and survival. Additional cost-effectiveness research can be used to determine the optimal threshold.
The diameter of the metastases was found to be the most important prognostic factor in this model. When building the microsimulation model, the focus was to make the model realistic with an appropriate level of complexity. The diameter predictor was an emergent property of the simulation model, which makes sense theoretically. The fact that the metastasis was invisible on the previous scan in combination with the size of the metastasis once detected, is informative about the growth rate of the metastases. A large metastasis therefore also has a high growth rate. Patients with fast-growing metastases have a relatively small chance of being detected as oligo+.
The tumour volume can also have a negative impact on survival. For instance, Oh et al. reported a 1.04 death hazard ratio for each cm 3 rise in volume of a brain metastasis [48]. However, in our model all metastases larger than 0.27 cm 3 were considered to be large, but these metastases may still be classified as small in the model of Oh et al. Furthermore, the model outcomes are different: Oh et al. predict the chance to die from a specific metastasis, while our model predicts the risk of having additional metastases. The findings from the two models are, however, not contradictory, as patients with larger brain oligo-metastases may both have a greater hazard to die from their oligo-metastases and have a lower risk of additional metastases. Both models suggest that those patients could benefit from curative treatment. Sensitivity analyses showed a high level of model robustness. Extreme parameter values resulted in minor effects on the ability to predict underlying metastases, and the outcomes still fall within a clinically observable range. Random variation of the parameters had little effect on model outcomes. Correlation between the growth rate and the numbers of metastases per patient improved the ability to distinguish between patients with and without additional undetected metastases. We expect that these features are associated; however, evidence for this hypothesis is lacking.
Using a different definition of oligo-metastatic disease, that is, detection of only one recurrence, affected the proportion of undetected metastases in the low-risk group. Extension of the definition to five recurrences did not increase the size of the low-risk group, but would have a negative effect on the relative number of unnecessary curative treatments if the "treat all" strategy was used. These findings match the literature in which the number of oligo-metastases per patient is also a prognostic factor [12].
When increasing the ability of metastases to metastasize, the size of the low-risk group diminishes accordingly. In the case that 80% of the patients have metastatic metastases, no low-risk patients with an indication for curative therapy remain identifiable. Although this model prediction is to be expected, it does not reflect the current clinical practice since there are patients that have long progression-free survival after curative treatment of oligo-metastases. Therefore, we argue that metastatic metastases may either be rare [49] or need more than 5 years to become detectable from the single-cell state.

Future Developments
There have been intriguing developments in molecular diagnostics to identify microscopic tumour load, such as circulating tumour cells (CTCs) and circulating tumour DNA (ctDNA). Although insight into the true potential of these techniques requires continued research efforts, preliminary modelling studies have investigated its potential in the clinical setting. For example, the study by Coemans et al. [50] studied whether CTCs could be used to detect distant metastases before surgery of the primary tumour in breast cancer. Tang et al. show that ctDNA may be predictive for distinguishing "true" oligometastatic disease from poly recurrence, although they have a small sample size [51]. As soon as clinical studies are available that elucidate the relationship between ctDNA or CTCs and oligo recurrence, it will be interesting to include this information into our simulation model.
Symptom emergence is largely dependent on the location and size of the metastasis. The data on the location of the metastases is currently not of high enough quality to be included in this model. Future versions of this model could be made that includes the locations of the metastases, using location-specific hazard rates to give symptoms.

Recommendations on Data Collection
Most studies on curative treatment for oligo-metastases in NSCLC are either prospective, single-arm studies, or retrospective studies, which have a small sample size of patients or report below the number of parameters that are required for our research question. These studies are often based on highly selected patients with favourable inclusion criteria, using varying definitions of oligo-metastases. These studies are therefore susceptible to selection bias [17,52]. Two randomized studies have been reported with 29 and 49 enrolled patients [53,54]. Although the results show increased progression-free survival, strong evidence for a benefit of curative therapy for oligo-metastatic NSCLC is still absent.
In the absence of high-quality patient-level data, we have used microsimulation to provide guidance and generate new hypotheses. To allow validation and improvement or extension of the model, structural reporting of the tumour and patient features in a prospective manner would be highly valuable, preferably within the context of a randomized controlled trial investigating treatment for oligo-metastatic disease.
Currently, three phase 2 trials in oligo-metastatic disease, including lung cancer (NCT03905317, NCT03349203, NCT03965468), are ongoing. These studies, finishing in 2021, focus on the treatment of oligo-metastases and compare two treatments for that purpose. It is not clear whether the identification of prognostic and predictive factors is an additional goal of these studies, although it is possible that the studies' data may be used for this purpose. Alternatively, the hypothesis of this model could also be tested in different types of cancer, such as colon cancer. To validate the model we present here, reporting on the diameters of the detected recurrences, the affected organs, the number of recurrences, and mode of detection (by symptoms or not, and which symptoms) is needed. To allow potential further improvement of the model, factors such as genetic markers, ctDNA and CTCs, histology of the primary tumour, or other prognostic factors [12] may be investigated for this specific purpose. These factors are not included in the current model.   (4) m det is the number of detected metastases and is per definition an integer number between 0 and M total . The number of detected metastases is equal to the number of metastases of detectable size (V det ) at the time that a scan occurs (t scan ). Metastasis i reaches detectable size (V det ) at t det i ; thus, i metastases are detectable when t deti+1 ≥ t scan ≥ t deti . Furthermore, the detection threshold should also fall in between V i and V i+1 at the time of this scan: V i+1 (t scan ) ≤ V det ≤ V i (t scan ). mdet can be obtained by solving i.
Applying Function (1) on the formula above results in: And according to Function (2): This can be simplified to: m det can be obtained by rounding i to an integer number: VDT×log 2 (R) = m det .

Appendix A.2. Sensitivity Analyses: Alternative Calibration Targets
The accuracy of the risk groups was tested by generating a new patient database using the regular prediction model with extreme parameters.
For each of the four model parameters that were obtained by calibration, one by one, the upper and lower limit of the 95% confidence intervals around the calibration targets were used to recalibrate the model. For the pooled averages of the rates (µ) from the literature (for oligo-metastases detected, oligo-, and symptomatic detections) the 95% confidence interval around this rate was estimated to be µ ± 1.96 µ (1−µ) n . These recalibrated microsimulation models were run again, and the risk sub-group analyses were repeated. All the other model parameters were left unchanged.

Appendix A.3. Sensitivity Analyses: Random Normal Variation around the Model Parameters
In the second set of sensitivity analyses, the parameters VDT and V det were separately varied around the point estimates by randomly drawing a new value from a normal distribution per patient. The scan detection threshold was altered around 5 mm with a standard deviation of 1 mm. The individually drawn volume doubling time was altered for the smallest detectable metastases with a standard deviation of 10 days. The microsimulation model was run again, and the risk sub-group analyses were repeated. All the other model parameters were left unchanged.

Appendix A.4. Sensitivity Analyses: Correlation between Volume Doubling Time and the Total Number of Metastases
In the Base Case model, no interactions were assumed between model parameters, because proof for such interactions was lacking. However, it is a credible assumption to make that patients with higher VDTs also have higher numbers of metastases. We explored the effect of such an interaction by adding a correlation between VDT and M total using the formula: with ρ as the "correlation coefficient", and r as a random uniform value that is drawn to determine the VDT of a patient's metastases. Two simulations were used for the sensitivity analyses, with ρ = 0.5 and ρ = 1.0. The microsimulation model was run again, and the risk sub-group analyses were repeated. All the other model parameters were left unchanged.

Appendix A.5. Sensitivity Analyses: Adaptation of the Definition of Oligo-Metastases
For these sensitivity analyses, the simulated patient-level database was unaltered, but the prediction model was modified by setting the maximum number of detected recurrences that is considered to be an oligo-metastasis Moligo max to 1 and to 5. This change in definition alters the number of patients with detected oligo-metastases and the proportions of oligo+ and oligo− within the risk sub-groups.
Appendix A.6. Sensitivity Analyses: Metastasizing Metastases A sensitivity analysis was included that estimates the impact of the unverified assumption that both oligo and poly metastases are able to produce additional metastases, without the presence of the primary tumour. This increases the proportion of oligo+. For this sensitivity analysis, we assumed that the mutations required for the metastasis to metastasize are related to the number of DNA duplication steps and cell divisions (A4). Therefore, we assume that the hazard of metastases acquiring the ability to metastasize is related to the total tumour volume in time: Here, p mm is the probability that one metastatic cell within the total tumour volume will acquire the ability to disseminate metastases, and V total (t) is the total tumour volume in time. The hazard was used to estimate how many of the patients with oligo-metastases in the model have developed additional metastases after removal of the primary tumour. p mm was calibrated to create scenarios with 20%, 50%, and 80% metastatic metastases ( Table 1). The proportions of oligo+ per risk sub-group was subsequently adjusted for patients that switched from oligo− to oligo+. Figure A1. Detection of recurrences in sensitivity analyses. Sensitivity analysis with t detectable calibrated to upper 95% CI of the progression-free survival (a) and the lower 95% CI of the progression-free survival (b). Figure A2. Time required for oligo+ to become distinguishable as poly-recurrent disease. This is estimated as the time that the 4th metastasis passes the detection threshold. 73.8% of oligo+ metastases become detectable after 3 months, and 98.2% of the oligo+ becomes detectable after one year. Data is based on the extrapolation of the microsimulation-model.            Proportions of oligo+ per risk sub-group. Grey fields represent the low-risk group. This scenario has fewer sub-groups than the Base Case. Table A14. Proportion of oligo+ within patients with detected oligo-metastases in the redefinition of oligo metastases to 1-5 metastases scenario.

Asymptomatic Symptomatic
Metastases detected: Proportions of oligo+ per risk sub-group. Grey fields represent the low-risk group. This scenario has more sub-groups than the Base Case.