Can evolutionary therapy be applied in non-small cell lung cancer?

Introduction
Resistance to therapy remains one of the most formidable challenges in the treatment of metastatic cancer. While many tumours initially respond to therapy, sustained responses are rare due to the evolution of drug-resistant cell populations. This is particularly problematic in metastatic disease, where tumour heterogeneity and evolutionary pressures are prevalent1. Traditional treatment protocols that rely on maximum tolerated dose (MTD) often accelerate resistance and compromising long-term outcomes2–9. Moreover, MTD-based regimens frequently reduce patient quality of life due to high drug toxicity.
In metastatic non-small cell lung cancer (NSCLC), the leading cause of cancer-related death worldwide2, the introduction of tyrosine kinase inhibitors (TKIs) has markedly improved survival in patients with actionable driver mutations10. Erlotinib, a first-generation EGFR-TKI, is commonly used in this context, but resistance to therapy eventually develops in nearly all patients11,12. This challenge has motivated the exploration of alternative dosing strategies that can delay resistance evolution while reducing toxicity. Intermittent and metronomic regimens have shown promise in some cancers, such as prostate and melanoma13–15, and improved quality of life in elderly NSCLC patients treated with chemotherapy16. Preclinical studies in NSCLC have also shown that resistance to erlotinib may decline in drug-free conditions17, and pulsed high-dose regimens do not increase toxicity in phase I trials18.
One promising alternative to standard of care is evolutionary therapy, which applies principles from evolutionary biology to adjust treatment in response to tumour behaviour19–25. Rather than aiming for maximal cell-kill, evolutionary therapy typically seeks to exploit ecological interactions among cancer cell populations, such as competition between drug-sensitive and resistant cells, to delay resistance and prolong tumour control25,26. A pioneering example of this approach is the trial by Zhang et al. in metastatic castrate-resistant prostate cancer, where adaptive on–off abiraterone dosing significantly prolonged time to progression27–30. This success has sparked interest in extending evolutionary protocols to other cancer types28,31,32. Its applicability to metastatic NSCLC remains unexplored, in part due to challenges in biomarker availability26,33, and the aggressive nature of NSCLC. In this context, mathematical models provide a valuable tool to explore and optimise evolutionary therapies in silico before clinical testing3,20,24,34–36. They enable comparison of different dosing strategies and how they affect clinical outcomes.
A recent application of mathematical modelling to NSCLC was in a study by Yin et al., where a variant of the polymorphic exponential model was fitted to longitudinal tumour burden data from NSCLC patients treated with erlotinib37. Their model includes drug-sensitive and drug-resistant cancer cell populations and assumes a constant, treatment-independent mutation rate from sensitive to resistant cells. It does not, however, incorporate density-dependent or frequency-dependent interactions, which have been shown to be key factors in cancer growth and evolution38–41.
In this study, we investigate whether a specific ET protocol could theoretically work in patients with metastatic NSCLC. We choose the protocol of Zhang et al. as it is the best known ET protocol and has already been successfully trialled28. Unlike Yin et al.’s exponential model, which assumes no ecological interactions, we investigate a broader class of two-population models that explicitly include density- and frequency-dependent competition. To do this, we systematically validate 24 alternative two-population models based on classical cancer dynamics42,43, using the same dataset analysed by Yin et al., and identify which models best capture the dynamics of metastatic NSCLC under erlotinib.
We then use the best-fitting models to simulate the evolutionary therapy protocol of Zhang et al.27, which uses on–off dosing to maintain competition between the two cell types, and we evaluate its predicted impact on time to progression (TTP), defined as the time until tumour burden exceeds 1.1 its initial value-a common proxy in ET studies.

Results

Mathematical models

Population growth models
We consider three base models for the population growth of sensitive cells (with population at time t) and resistant cells (with population size at time t): the Logistic model, the Gompertzian model, and the von Bertalanffy model42,44. In all of these models, and refer to the growth rate of the sensitive and resistant cancer populations, respectively, refers to the competition effect of sensitive cells on the resistant population, refers to the competition effect of the resistant cells on the sensitive population, and K serves as a carrying capacity. These parameters are fitted to longitudinal volumetric measurement data as described in the section titled ‘Data Sources and Processing’ in the Methods. For all models, the sensitive population grows as long as , and the resistant population grows as long as . If , then growth depends only on the total cancer population. However, when these parameters are not fixed to 1,  frequency dependence arises, which means that the growth of each cell type depends on the relative frequencies (or ratios) of the cell types within the population. The Logistic two-population model is given byThe Gompertzian two-population model takes the formThe third model that we consider, the two-population von Bertalanffy model45,46, is given by

Density and frequency dependence
We include two types of competition in our models: density-dependent and frequency-dependent competition. Density-dependent competition arises when cells compete for space and/or limited resources. This is modelled using the carrying capacity K. All our models include density-dependent competition. Frequency-dependent competition, on the other hand, depends on the relative abundance of different populations. This means the impact of one population on the other’s growth depends on their relative proportions. Frequency-dependent competition is included in our models by using the competition coefficients and to quantify how strongly the populations influence each other’s growth. For each growth model, we fit two different versions: a purely density-dependent model, where we set , and a frequency-dependent model, where in addition to the density dependence, we include frequency dependence by modelling and as free parameters which are estimated from the data.

Pharmacokinetics (PK) and pharmacodynamics (PD)
Following Yin et al.37, we choose a linear PD model, i.e., the effect of the medication increases linearly with the treatment dose47. Treatment-induced cell death depends on how much treatment reaches cancer cells. Plasma drug concentration varies between treatment doses. Generally, in PK modelling, an exposure metric, , such as area under the concentration-time curve (AUC), maximum or minimum observed concentration ( or ), lowest concentration observed before the next dose (), or the dose (m(t)) is used48,49. These are examples of time-aggregated metrics where the dynamics of the drug concentration are ignored and information is condensed into a single number. We refer to such models as aggregated-PK (agPK) models.
Alternatively, a time-varying metric may be used, such as the concentration predicted by a PK model fitted to PK data48,49. In this case, the concentration dynamics are accounted for and affect the tumour growth. We term such models dynamic-PK (dyPK) models. Here, we used a two-compartment model with first-order absorption and linear elimination similar to Yin et al. 37 (see Appendix 2 for detailed explanation).
In our models, we define a medication function c(t) that captures both the above-mentioned cases:with the maximum tolerable dose defined in milligrams (mg). The physician’s decision is captured by the medication ratio m(t) at time t and lies in the interval [0, 1],  with and corresponding to no treatment and maximum tolerable dose, respectively, and C(t) gives the plasma concentration normalised to the treatment dose at time t (see Appendix 2 for its definition). The parameters and are scaling factors with units and respectively. In the rest of this paper, comparisons of the agPK and dyPK models are made. It is important to note that the results for the agPK models are independent of the choice of exposure metric . In other words, if we replace the treatment dose () at time t by, for instance, AUC or , the results will not change.

Treatment-induced death rate
We consider two standard assumptions regarding treatment-induced death of treatment-sensitive cancer cells: the log-kill (LK) model and the Norton-Simon (NS) model50–53. With the log-kill treatment-induced death rate, the population of sensitive cancer cells killed at time t is with denoting the treatment effect, a parameter to be fitted. Hence, for example, the Logistic model with the LK treatment-induced death rate becomesA key principle of the NS hypothesis is that the treatment-induced death rate is proportional to the sensitive cells’ growth rate. The effect of treatment is thus incorporated by multiplying the right-hand side of the evolution equation for by where represents the treatment effect and the medication function defined in (4).
As an example, applying this modification to the Logistic growth model yields the following system:We consider three growth models, two alternative assumptions regarding frequency dependence, two alternative assumptions regarding the pharmacokinetics, and two different treatment-induced death rate functions, giving 24 two-population models. We compare the ability of these models to capture NSCLC dynamics under erlotinib treatment to that of the polymorphic exponential model37 given byThe polymorphic exponential model assumes a single growth rate r, i.e., . The expression denotes treatment-induced death rate of sensitive cancer cells, with parameter and medication function c(t),  defined by (4). It also assumes a mutation rate of sensitive cells into resistant cells, independent of medication. For this model we also consider two alternative assumptions for pharmacokinetics. This model has been previously fitted37 to the same dataset that we use for our simulations.
Table 1 summarises the key variables and parameters of the 24 models we propose, as well as those of the polymorphic exponential model.

Models using the log-kill death rate fit better than those using the Norton-Simon death rate
To fit model parameters to individual patient data, we used differential evolution, a population-based heuristic optimisation algorithm for non-linear problems. The method iteratively improves candidate solutions through mutation, crossover, and selection, and is often considered less sensitive to premature convergence to local minima than gradient-based approaches54,55.
We fit the model parameters both individually, where each patient has unique fitted parameters, and population-wide, where the same parameter values are fitted for all patients. For individual fits, we minimised the relative percentage error, whereas for population-wide fits, we minimised the average relative percentage error across all patients.
There is a considerable difference in the goodness of fit of models with different treatment-induced death rates. Figure 11 in Appendix 4 shows that the LK death rate decreases the relative percentage error of individual fits for all models compared to the NS death rate. A similar effect is observed for population-wide fits, where acceptable goodness of fit is achieved only by the Gompertzian and Logistic models with a log-kill treatment-induced death rate (Table 2). Given these results, we continue analysing only models with the log-kill death rate in the rest of the study.

Pharmacokinetic type does not impact model goodness of fit
We observed that the choice of pharmacokinetic modelling type, aggregated pharmacokinetics (agPK) or dynamic pharmacokinetics (dyPK), did not significantly affect the goodness of fit. Among the five best-fitting models, only the poorest-fitting one (the Logistic model with log-kill dynamics and without frequency dependence, LG_LK) showed a statistically significant difference in paired t-tests (see Table 7 and Fig. 1). The measured effect was less than 0.1%, suggesting that the difference in goodness of fit is negligible in practical terms. This indicates that, from the perspective of goodness of fit, both pharmacokinetic types are equally suitable. Likewise, on a population-wide level, the best-fitting models show little to no difference in relative percentage error between the two pharmacokinetic types (Tables 2 and 7). Boxplots illustrating the distribution of relative percentage errors for the five models with the best individual-level fits are shown in Fig. 1, while boxplots for all models are provided in Appendix 4. The polymorphic exponential model no longer provides an adequate fit, whereas the other four top-performing models maintain good and consistent fits (Table 2).

The polymorphic exponential model is not suitable for predicting disease progression or guiding evolutionary therapy
The populations of sensitive and resistant cancer cells in the polymorphic exponential model at any time t can be calculated analytically (see Appendix 3):Here, r is the growth rate for the sensitive and resistant populations (denoted and , respectively), is the effect of medication on the growth of the sensitive cell population, is the mutation rate, and k is the ratio between the initial resistant and sensitive populations; refer to Table 1 for parameter definitions. It can clearly be seen from the equations in Appendix 3 that the size of the resistant population is entirely independent of whether any treatment has been given. Referring back to the model given by (7), we see that the mutation rate is simply directly proportional to the sensitive population.
Moreover, in the same appendix we show that the time t at which the tumour burden reaches any proportion p of the initial tumour burden satisfiesa quantity independent of the initial tumour burden. This implies that the time before exceeding any given progression threshold of the polymorphic exponential model is a function of the parameters and independent of the initial tumour burden. For population-wide fits, the parameters are identical for all patients, meaning all patients progress at the same time. Since the mutation rate is proportional to the sensitive population and independent of the medication, when medication is paused and the sensitive population increases, this in turn causes the growth rate of the resistant population to increase. The absence of a competitive effect and the medication-independent mutation rate mean that the model may predict worse results under evolutionary therapy compared to MTD. Although the polymorphic exponential model can capture tumour burden data relatively well for many patients (Appendix 4), the mathematical properties discussed above disqualify its further usage for treatment optimisation.

The Gompertzian model with frequency dependence most accurately describes NSCLC treatment dynamics
For both individual and population-wide fits, the Gompertzian model with frequency dependence has the best relative percentage error (Fig. 1 and Table 2). Moreover, of all the models tested, this model is best able to capture the U-shaped trajectory typical for treatment-induced resistance (Fig. 2).
We, therefore, proceed with the Gompertzian model with frequency dependence for simulations of evolutionary therapy. We select Zhang et. al’s therapeutic protocol, which aims to keep the tumour burden below the progression threshold by strategically switching the treatment on and off, allowing the sensitive cell population to regrow in order to keep the resistant population in check27. For comparison, we also simulate evolutionary therapy on the other best-fitting models and show this alongside the polymorphic exponential model (Fig. 1).

The best-fitting models predict increased time to progression with Zhang et al’s protocol
Our simulations using the fitted parameters show that Zhang et al.’s protocol increases time to progression (defined as 1.1 times the initial tumour burden) when compared to the standard of care. We employ simulations to demonstrate the size of the potential benefit of Zhang et al.’s protocol for patients treated with erlotinib. As an illustrative example, we show the Gompertzian model with frequency dependence fitted to one patient in Fig. 3.
Figures 3 and 4 demonstrate expected time to progression (defined as 1.1 times the initial tumour volume) with MTD and Zhang et al.’s evolutionary therapy with the Gompertzian and Yin et al.’s models, respectively, for one of the erlotinib patients. While both models predict similar growth dynamics under MTD, they predict very different outcomes under an evolutionary therapy protocol. Similarly to the Gompertzian model, all best-fitting models presented in Fig. 1 predict a qualitatively similar benefit of Zhang et al.’s evolutionary therapy when compared to standard of care, as shown in Fig. 5. The simulations in Fig. 5 were continued beyond the initial span of the data. If the simulated trajectories had not reached the progression threshold during the simulated time, the data were right censored. In the Gompertzian and Logistic models, censoring occurred only during the simulations of evolutionary therapy as the dynamics for some patients allowed multiple repeated cycles. In contrast, when simulating with the polymorphic exponential model, censoring was also required for predictions under MTD, as the model predicted very slow long-term growth.
The median predicted TTPs for both MTD and Zhang et al.’s evolutionary therapy in models with frequency dependence (FD) are lower than those in models without FD (Fig. 5). Models with FD also exhibit lower RMSE values and fit the data more accurately (Fig. 1). This is expected since these models include two additional parameters compared to models without FD. For all patients, the optimised parameters in both the Gompertzian and logistic models with FD consistently satisfied , indicating asymmetric competitive effects. In both models, the change in predicted TTP between models with and without FD depends on the ratio of the competition coefficients . For instance, in the Gompertzian models GM_LK_FD (with FD) and GM_LK (without FD), all patients with exhibited higher predicted TTP under the Zhang et al.’s protocol in GM_LK_FD than in GM_LK. For patients with , the opposite was observed. The median TTP decreases in both Gompertzian and logistic models when FD is included because the majority of patients in our dataset (approximately 77% for Gompertzian and 69% for logistic) had .
Additionally, the TTP gains achieved with Zhang et al.’s evolutionary therapy compared to MTD were greater for Gompertzian models than for logistic models (Fig. 5). This is consistent with previous theoretical studies showing that the reduction in growth rate for a marginal increase in total tumour burden is greater in the Gompertzian than in the logistic model of tumour growth34,56. Consequently, competitive effects are stronger in Gompertzian models, resulting in larger TTP gains than in logistic models.
Figure 6 presents Kaplan–Meier plots of the expected time to progression under Zhang et al.’s evolutionary therapy and MTD for all patients, for the Gompertzian and polymorphic exponential models, respectively. The Gompertzian model with frequency-dependent competition predicts a substantial benefit of evolutionary therapy, with the median time to progression increasing from 24.8 months under MTD to 42.3 months, an improvement of 17.5 months. All well-fitting models that include competition show a similar benefit. Kaplan–Meier plots for the other well-fitting models, and are provided in Appendix 7. In contrast, the polymorphic exponential model, which lacks both density- and frequency-dependent competition, predicts the opposite trend. It also suggests that some patients do not progress at all under MTD during the simulated time frame, which occurs when a high fitted mutation rate is offset by a very low growth rate of resistant cells. Figure 4 alsor illustrates that the absence of competition in the polymorphic exponential model leads to worse outcomes under Zhang et al.'s evolutionary therapy, reinforcing its unsuitability for ET protocol design (see Appendix 3 for analytical justification).

Discussion
This study evaluated the potential of evolutionary therapy for non-small cell lung cancer (NSCLC) using mathematical models fitted to tumour burden data from patients treated with erlotinib. We first evaluated the goodness of fit of multiple two-population models, at both the individual and population levels, to identify those that best described the observed tumour dynamics. The models that provided the best fits were then used to simulate Zhang et al.’s evolutionary therapy protocol27. Time to progression, defined as the time when tumour burden exceeds 110% of its initial value, was compared between this protocol and standard-of-care treatment.
Our analysis showed that although several models can achieve a satisfactory goodness of fit to patient data, their ability to describe tumour dynamics and resistance evolution varies substantially. Among the fitted models, the two-population Gompertzian model with frequency dependence provided the best overall goodness of fit, capturing tumour-burden trajectories in metastatic NSCLC patients treated with erlotinib more accurately than alternative models. This emphasises the importance of including parameters that represent density-dependent interactions as well as frequency-dependent interactions between cancer cell types, when developing models to optimise systemic therapies. This is consistent with previous in-vitro studies38,39 demonstrating frequency-dependent competition in the tumour microenvironment.
We further demonstrated, using a mathematical proof, that in Yin et al.’s polymorphic exponential model37, the predicted time required for the tumour burden to exceed any fixed multiple of its initial value is independent of the initial size and determined solely by the fitted parameter values (Appendix 3). This structural property implies that, although the model can achieve reasonable goodness of fit to the observed data, it cannot capture treatment-dependent differences in resistance dynamics. Consequently, the polymorphic exponential model is not suitable for guiding evolutionary therapy in NSCLC. In contrast, for all two-population models deemed suitable for therapy guidance, simulations of the Zhang et al.’s evolutionary protocol27 predicted substantially longer times to progression compared with standard treatment.
The primary objective of evolutionary therapies is to forestall or delay the development of treatment-induced resistance. This is often done through predicting and controlling cancer cells’ response to therapy26,32,57. This study demonstrates that not all models are capable of capturing this. Although many models, including the Yin et al.’s polymorphic exponential model, fit the data well when the tumour burden shows a monotonic decrease, not all models can predict the U-shaped patterns typically observed when treatment-induced resistance develops. To combat resistance, being able to capture its development is vital. Therefore, though they may fit data well, models that predict long-term tumour suppression at the MTD should be approached with caution, as they may not be able to guide evolutionary therapy.
Even among the models that we found to be suitable for describing tumour dynamics in NSCLC, there are differences in the predicted time to progression (TTP). We find that the predicted TTP depends strongly on the ratio of the competition coefficients, . Importantly, implies that resistant cells affect the growth of sensitive cells less than sensitive cells affect resistant cells, and hence retaining a large sensitive population yields a higher TTP under Zhang et al.’s protocol and a lower TTP under MTD. Conversely, implies that resistant cells have a competitive advantage over sensitive cells, resulting in a lower TTP under Zhang et al.’s protocol and a higher TTP under MTD. In models without FD, by contrast, the two populations are assumed to interact symmetrically.
Furthermore, the increase or decrease in predicted TTP can depend on the existence, reachability, and stability of interior equilibria, which are determined by the intraspecific competition coefficients58. In our dataset, most patients had fitted parameters such that . Although the median decrease in TTP predicted by models with and without frequency dependence was small, for some individuals the difference was substantial. Notably, Zhang et al.’s evolutionary therapy consistently led to longer TTP than MTD across all models, regardless of FD, underscoring the importance of accounting for asymmetric competition between sensitive and resistant cancer-cell populations.
Our study provides a methodological framework for evaluating evolutionary therapy protocols using longitudinal clinical data. By systematically comparing 26 models with varying assumptions about competition, pharmacokinetics, and treatment-induced death rate, we demonstrate how model structure shapes both goodness of fit and mechanistic interpretation. We deliberately focused on qualitative two-population models that capture population-level resistance dynamics while remaining interpretable and clinically relevant. Although the Gompertzian model with frequency-dependent competition provided the best fits for this NSCLC dataset, other cancer types, drug combinations, and treatment regimens may yield different optimal model structures. For instance, in contrast to other research59, we found that models fitted to this dataset did not achieve good fits under a Norton–Simon treatment-induced death rate, whereas using a log-kill assumption significantly improved goodness of fit. This is likely due to the specific drug–disease combination of NSCLC and erlotinib. We also found that including dynamic PK effects did not significantly improve the model’s goodness of fit in this dataset compared with aggregate PK.
This study is positioned within a broader landscape of alternative treatment strategies that aim to personalise cancer care and contributes to the literature on the model-informed precision dosing 60,61. While therapeutic drug monitoring62 focuses on adjusting doses to achieve target drug exposures, model-informed precision dosing uses models to adjust dosing based on patient characteristics.
Evolutionary therapies consider the ecological and evolutionary dynamics of the cancer cells21,27,28,32,63. The protocol simulated in this study is evolutionary as it is based on models of intra-tumoural evolution. Other approaches also based on such models include dose stabilisation14,23,64, intermittent dosing16,65, extinction therapy66 and the double-bind therapy67,68. The models used in our study may also be extended to simulate and develop these alternate strategies and/or find the best evolutionary therapy through optimisation, similarly to69,70. We believe that targeting resistant populations is the necessary next step in personalised therapy and will lead to patient- and treatment-specific treatment regimens which will improve patients’ quality of life and survival26.
Models for guiding evolutionary therapy must be based on robust patient data and incorporate key features of cancer biology, such as competition between distinct cancer cell subpopulations23,32,58,71–73. Bringing evolutionary therapy to clinic necessitates a collaborative approach, where pharmacologists provide insights into drug action and kinetics, mathematicians develop and refine the models, and clinical oncologists ensure biological plausibility and translate findings into patient benefit26. Such interdisciplinary efforts are crucial for creating models and developing therapies that will benefit patients.
Our findings emphasise the shortcomings of overtly simplistic models, such as the polymorphic exponential model, which assumes a constant mutation rate independent of treatment and neglects competition. In contrast, models with frequency-dependent selection can better capture the interplay between drug-sensitive and drug-resistant cells, improving predictive accuracy74. The consistent prediction across the best-fitting models that Zhang et al.’s evolutionary therapy outperforms the standard of care for all studied patients is promising. These findings provide strong theoretical support for the use of evolutionary therapy to improve care for NSCLC patients.
Our analysis indicates that, in the context of evolutionary therapy as applied in the protocol of Zhang et al. for stage IV NSCLC, the choice between dynamic and aggregate PK representations does not significantly affect the goodness of fit. Dynamic PK may, however, be important for predictions in alternative evolutionary therapy protocols, such as dose titration14.
In this work, we have demonstrated that evolutionary therapy can theoretically be applied in NSCLC. We showed that two-population models, which include density- and frequency-dependent competition, can capture resistance dynamics and predict treatment outcomes. These findings provide a theoretical foundation for extending evolutionary therapy to NSCLC and underscore the need for empirical validation before clinical implementation. This will require larger datasets to ensure that the data used are representative of the broader NSCLC population treated with erlotinib.
Notably, in the dataset employed here, not all patients had recorded measurements indicating tumour regrowth. Although all patients are palliative and known to progress, this event often occurs outside of the observation window. As a result, some of the models fitted to this dataset failed to predict progression. In future work, the absence of progression data could be addressed by incorporating population-level parameter estimates from progressing patients in order to constrain the range of parameter values for these patients.
The future clinical implementation of evolutionary therapy for NSCLC will require overcoming several practical challenges26. Validation through clinical trials will be essential to confirm the efficacy and safety of the proposed protocols. A key challenge for fast-growing cancers such as NSCLC is the need for frequent and detailed measurements of tumour burden. One promising avenue is the use of circulating tumour DNA (ctDNA) as a biomarker75. Collected through a simple blood sample, ctDNA is less invasive and more cost-effective than imaging, allowing for more frequent follow-up. Moreover, it can provide higher-resolution information on tumour composition, leading to more accurate and responsive model predictions. In the current study, we considered a single resistant population, which may limit predictive accuracy at the individual level, as intra-tumour heterogeneity can cause divergent responses across cancer cell types1,3,76. Incorporating time-series ctDNA measurements in future work will enable more dynamic, patient-specific model updates and further improve clinical relevance.
Our findings provide theoretical support for extending evolutionary therapy to NSCLC. By fitting mathematical models directly to patient data and identifying key mechanisms such as frequency-dependent selection, we offer a modelling framework to evaluate and optimise evolutionary treatment strategies. With further validation and clinical integration, these model-informed approaches may contribute to personalised regimens that improve long-term tumour control, limit resistance, and ultimately enhance patient survival and quality of life.

Methods

Data sources and processing
We performed our experiments using the two data sets containing patient data from the START-TKI trial (NCT05221372)77, previously analysed in 37. All patients are stage IV patients with Non-Small Cell Lung Cancer treated with erlotinib. The data includes: Longest diameters of 2-7 tumours for 18 patients measured approximately every 6 weeks, with 2-18 records per patient. All 18 patients were initially given the MTD dose of 150 mg of erlotinib. Of these, 4 patients had their doses reduced to 100 mg or 50 mg due to the effects of drug-induced toxicity. For 16 of these patients, pharmacokinetic (PK) data (blood concentration levels) of erlotinib measured between 0.25 and 30.5 h after the last dose was administered. We refer to this as ‘sparsely sampled’ PK data.

PK data from 30 patients treated with 50–150 mg erlotinib daily with blood samples collected before treatment initiation and at 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 6, 8, 12, and 24 h after drug administration. We refer to this as ‘intensively sampled’ PK data.

From the tumour measurement data (18 patients, 2-7 tumours, 2-18 time points), we removed measurements from dates with some of the tumours unmeasurable. For the purpose of fitting the data to ordinary differential equation models, we also excluded patients with fewer than 6 measurements in total. We converted the longest tumour diameter measurements to volumes by assuming the tumour to be spherical. For each patient, we summed up the volumes of all their tumours to obtain their total tumour burden at each measurement date. This left a final tumour volume data set consisting 13 patients with total tumour burden information for 6-18 time points. PK data was processed by binning the sparsely sampled data into intervals of 0.5–3 h before merging with the intensively sampled PK dataset. This larger dataset was then used to create a population-level PK model.

Ethics statement
This study involved secondary analyses of individual-level data from the START-TKI trial (ClinicalTrials.gov ID: NCT05221372). The data themselves have not been published, although analyses based on them have been previously reported (Yin et al., CPT: Pharmacometrics & Systems Pharmacology, 2024; 13(4): 612–623. https://doi.org/10.1002/psp4.13105). The data were obtained under a data use agreement with the original investigators at Erasmus Medical Center.
The START-TKI trial was conducted in accordance with the Declaration of Helsinki and all relevant guidelines and regulations. Ethical approval for the trial was granted by the Medical Ethics Review Committee (METC) of Erasmus Medical Center, in accordance with the Dutch Medical Research Involving Human Subjects Act (WMO), as registered on Onderzoekmetmensen.nl and ClinicalTrials.gov (NCT05221372). Informed consent was obtained from all participants prior to enrolment in the trial.
All data used in this study were de-identified prior to access, and no additional ethical approval was required for this secondary analysis.

Pharmacokinetics
Pharmacokinetics (PK) of a drug describes how its concentration in the blood changes over the drug dosing interval. In the data we use in this study, patients typically took erlotinib once a day, resulting in a 24-h cycle in the medication dose dynamics. To study any effects drug concentration dynamics may have on tumour dynamics and our predictions, we investigate versions of each model with agPK and dyPK, explained in detail in the section titled ‘Pharmacokinetics (PK) and pharmacodynamics (PD)’ in the Results. For dyPK, the plasma concentration of the drug (in ) at time t given a dose D(t) (in mg) is given by , where C(t) describes the dynamics of the plasma concentration according to a two-compartment PK model with first-order absorption and lag time, similar to Yin et al.’s study37 (see Appendix 2 for details). The dose D(t) is multiplied by 1000 to convert it to from mg. This can also be written as where m(t) is the medication dose ratio and is the maximum tolerable dose in mg.
For agPK models, we assume that the entire dose is absorbed and the concentration stays constant, and thus .
The fitting procedure used for the PK model is explained in detail in Appendix 2. Of the six parameters of the model, four parameters () were fitted and two parameters were fixed (), based on the findings of Yin et al.37, due to their insensitivity to small changes.

Model fitting and parameter estimation
To investigate suitable models for use in evolutionary therapy, we fitted the data of 13 patients to two-population mathematical models. We considered three standard models, with two treatment-induced death functions from the mathematical oncology literature, two pharmacodynamic options, and two options regarding frequency dependence, bringing us to 24 models (see section titled ‘Mathematical Models’ in the Results). We compared the fits of these models to the model of Yin et al., with both agPK and dyPK.

Parameter optimisation using differential evolution
We fitted the tumour volume data to all model combinations (see section titled ‘Mathematical Models’ in the Results) as well as to the polymorphic exponential model of Yin et al. The total initial tumour burden was defined as , with the initial population set to and .
The polymorphic exponential model as described by Yin et al., contains no initial resistant population. Instead, all resistance occurs through mutation. As mutation rate is constant and independent of medication, the structure of this polymorphic exponential model means that the resistant population is initiated immediately and always grows from the start of the simulation. Consequently, the dynamics of the model are not substantially changed by our including a small initial population of resistant cells.
The remaining model parameters were estimated using Python 3.11.4 with the differential evolution algorithm54 implemented in the scipy.optimize library. Differential evolution (DE) is a global optimisation algorithm designed for continuous, non-linear problems. It evolves a population of candidate parameter vectors through mutation, recombination, and selection54. In each generation, parameter vectors are randomly perturbed and recombined to explore the parameter space. The recombined vector’s cost (or goodness of fit in the context of this paper) is compared to the target vector’s cost, and the vector with the lower cost is retained for the next generation. The optimization process terminates when the standard deviation of the population’s costs falls below , where is the relative tolerance. This convergence criterion ensures the population stabilises and no significant improvements are achievable. The model was fitted to the data by solving the differential equations with the Runge–Kutte 45 method from the solve_ivp module in the scipy.integrate library. Although DE is stochastic and explores a large portion of the parameter space on a single run, we ran the fitting procedure five times per patient and selected the result with the lowest error. We used the default settings for differential evolution, except for decreasing relative tolerance to .

Parameter constraints
We fit the model parameters through differential evolution, applying bounds summarised in Table 3. The bounds for the growth rates and based on the literature. For this patient data set, Yin et al reported upper and lower bounds for the growth in tumour’s longest diameter. Converting diameter to tumour volume (assuming tumour is spherical) gave upper and lower bounds of 0.00255 and 0.0921, respectively. We also considered growth rates reported from in vitro studies, such as 39. This gave us growth rates being bounded between 0.01 and 0.05 for the Gompertzian and polymorphic exponential models. In the von Bertalanffy model   (Eq. 3), and are multiplied by the populations and raised to fractional powers. This results in smaller net growth rates compared to the other models. This means that the values of and need to be higher to achieve the same net growth rates of the population. For this reason, the upper bounds for and for the von Bertalanffy model were set to 0.5.
Individual bounds for the carrying capacity K were set between 1.5 and 5 times the initial tumour burden.
As where mg of erlotinib, we set different upper bounds for , specifically 0.001 for models with LK, and 0.01 for models with NS.
The lower and upper bounds for the competition coefficients and were set at 0.1 and 10,  respectively. We chose to explore the high upper bound due to the findings of relatively high competition coefficients in recent studies, such as those of Freischel et al.38.

Individual and population-level parameter estimation
Each model was fitted for each patient using differential evolution to minimise the following cost function:Here, is the total tumour volume of the ith patient at the start of treatment, and RMSE is the root mean squared error. We converted RMSE to relative percentage error in order to simplify comparison of fits between patients with different total tumour burdens. By minimizing this cost function, we determined the individual best-fit parameters for each patient.
Although we have confirmation that all patients in the dataset progressed, this is often not reflected in the data. We have also observed that the earlier and later data points are most important for capturing the trend. When, as occurs in our data, the last data points do not indicate progression, it is difficult for an individual model to predict the tumour growth trajectory accurately, even though many models may fit the data well.
To account for this, we additionally sought for each model combination a single set of parameters that would fit best to the all individual patients, given the unique initial tumour burden of each patient. We fitted a cost function for population-level relative percentage error which calculates the error as a percentage of the patient’s initial tumour burden for each patient, then averages this over the population (Eq. 9).In this equation, n is the number of patients, and are the initial tumour volume and root mean squared error of the ith patient respectively.

Generalizing model fits across patients
In the population fitting we sought a single set of parameters that would give the lowest error across all patients, given their individual starting size. The results from the population fitting are shown in Table 2, and it can be seen that the Gompertzian model with fitted competition coefficients fits best to all patients with a relative percentage error of with and without PK. The fits of this model to all patients in the population are shown in Fig. 12.

Simulation of evolutionary therapy protocols
We used the models selected in the previous step to compare the predicted time before the tumour first exceeds 1.1 times its original size (time to progression) under the MTD treatment protocol with that under the evolutionary protocol of Zhang et al.78. The Zhang et al.’s protocol consists of two phases: with and without medication. With medication, the tumour growth is simulated using the optimised parameters for the selected model, with the dose ratio m(t) set to 1. Without medication, the medication dose is set to 0, keeping all other parameters the same.
Using the selected ordinary differential equations population dynamics, the simulation is initialised with medication on phase, starting with the provided initial total tumour burden. When the total predicted tumour burden decreases to of its initial value, the medication off phase commences. When the total predicted tumour burden returns to its initial value, then the medication on phase resumes. This cycle continues for the desired length of the simulation. To measure time to progression we defined progression as 1.1 times the initial estimated sum of tumour volumes. For each patient, we used the parameters found in the fitting process to simulate future tumour growth and calculate the predicted time to progression under the MTD and Zhang et al.’s protocols.
The impact of evolutionary therapy on all models was assessed using Kaplan-Meier survival plots, where the event of interest is TTP. A Kaplan-Meier plot is a non-parametric approach to survival analysis which allows the comparison of time-to-event data between different groups. Both MTD and evolutionary therapy scenarios were used to simulate tumour volume over a period 10 times the original duration. Patients who did not progress by the end of this period were assigned to the longest time to progression. We employed the log-rank test to compare survival distributions of the MTD and evolutionary therapy groups. As the log-rank test accounts for censored data, the log-rank test statistic, which follows a chi-squared distribution, provides an unbiased comparison of observed versus expected outcomes between groups.