Longitudinal nonlinear mixed effects modeling of EGFR mutations in ctDNA as predictor of disease progression in treatment of EGFR-mutant non-small cell lung cancer.

Correlations between increasing concentrations of circulating tumor DNA (ctDNA) in plasma and disease progression have been shown. A nonlinear mixed effects model to describe the dynamics of epidermal growth factor receptor (EGFR) ctDNA data from patients with non-small cell lung cancer (NSCLC) combined with a parametric survival model were developed to evaluate the ability of these modeling techniques to describe ctDNA data. Repeated ctDNA measurements on L858R, exon19del, and T790M mutants were available from 54 patients with EGFR mutated NSCLC treated with erlotinib or gefitinib. Different dynamic models were tested to describe the longitudinal ctDNA concentrations of the driver and resistance mutations. Subsequently, a parametric time-to-event model for progression-free survival (PFS) was developed. Predicted L858R, exon19del, and T790M concentrations were used to evaluate their value as predictor for disease progression. The ctDNA dynamics were best described by a model consisting of a zero-order increase and first-order elimination (19.7/day, 95% confidence interval [CI] 14.9-23.6/day) of ctDNA concentrations. In addition, time-dependent development of resistance (5.0 × 10-4 , 95% CI 2.0 × 10-4 -7.0 × 10-4 /day) was included in the final model. Relative change in L858R and exon19del concentrations from baseline was identified as most significant predictor of disease progression (p = 0.001). The dynamic model for L858R, exon19del, and T790M concentrations in ctDNA and time-to-event model adequately described the observed concentrations and PFS data in our clinical cohort. In addition, it was shown that nonlinear mixed effects modeling is a valuable method for the analysis of longitudinal and heterogeneous biomarker datasets obtained from clinical practice.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Correlations between increasing concentrations of circulating tumor DNA (ctDNA) in human plasma and progression of disease have been shown, suggesting that longitudinally collected ctDNA can be used to guide treatment decision making. However, the analysis of heterogeneous longitudinally collected data is challenging.

WHAT QUESTION DID THIS STUDY ADDRESS?

A nonlinear mixed effects model to describe the dynamics of EGFR driver and resistance mutation ctDNA data in combination with a parametric survival model was developed to evaluate the ability of these modeling techniques to describe ctDNA data from a clinical practice.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

A dynamic model for the driver and resistance mutation concentrations in ctDNA and time‐to‐event model adequately described the observed concentrations and progression‐free survival data in our clinical cohort.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY OR TRANSLATIONAL SCIENCE?

Nonlinear mixed effects modeling is a valuable method for the analysis of longitudinal and heterogeneous biomarker datasets obtained from clinical practice and provides a framework for prediction of progression based on observed ctDNA concentrations.

INTRODUCTION

Over the past two decades, new targeted therapies have been introduced in the treatment of several types of cancer. The pharmacology of these drugs is characterized by targeting specific biological processes that have been related to tumor development. As a result, new approaches to identify and to quantify the response to these drugs are needed. The identification of drug‐specific biomarkers can provide guidance in monitoring the response to treatment and can enable individualization of therapy, both based on the change in biomarker levels. To this end, mutation analyses in liquid biopsies has gained interest. These patient friendly methods are being used to determine the genotype of tumors in human blood or plasma. These methods allow for a less invasive way of taking repeated measurements compared to the traditional method in which a biopsy from the tumor is obtained. This may facilitate the determination of the quantitative change in biomarker response. , During treatment, selected drug‐specific genes can be quantified and monitored to predict either response or progression of the tumor.

It has been shown that 5–20% of patients with non‐small cell lung cancer (NSCLC) present with activating mutations in the EGFR gene. , Of these mutations, the L858R point mutation on exon 21 (L858R) and deletions on exon 19 (exon19del) of the EGFR gene are the most common. , Erlotinib and gefitinib are orally available tyrosine kinase inhibitors (TKIs) targeting these EGFR mutations. Both TKIs are ATP competitors at the ATP‐binding pocket in the intracellular kinase domain of EGFR and it has been shown that mutant kinases are more sensitive to inhibition by erlotinib and gefitinib. Significant improvements in progression‐free survival (PFS) and overall survival (OS) have been shown in patients that were treated with either erlotinib or gefitinib. However, after an initial response, the tumor is likely to acquire resistance mechanisms, resulting in relapse of the disease. , Resistance to EGFR inhibitors in NSCLC has been related to the development of a point mutation on exon 20 (T790M) in approximately 50% of patients. , , Several studies have shown correlations between increasing concentrations of EGFR mutants (L858R, exon19del, and T790M) that were measured by droplet digital polymer chain reaction (ddPCR) and progression of disease during treatment with erlotinib or gefitinib. , , , In addition, low concentrations of these mutants have been observed during response to treatment. , In addition, we previously reported on the dynamics of the L858R, exon19del, and T790 mutations in a NSCLC cohort and showed that an increase in EGFR driver mutations may predict clinical progression. However, a significant relationship between the increase in EGFR driver mutation and time to progression could not be identified in this cohort.

The analysis of heterogenous longitudinal data collected in clinical practice is highly challenging. Conventional statistical techniques are impractical because most of these require homogenous sampling. Therefore, more sophisticated techniques are required. The use of nonlinear mixed‐effects modeling enables the analysis of longitudinal and heterogeneous datasets as has extensively been described in the field of pharmacokinetics (PKs) and pharmacodynamics (PDs). In this study, the application of these modeling techniques to circulating tumor DNA (ctDNA) biomarker and survival data was explored.

METHODS

Patients, study design, and data

Data were available from an observational study that was performed in the outpatient clinic of the Netherlands Cancer Institute, Amsterdam, The Netherlands. In this study, ctDNA from patients with NSCLC who received erlotinib or gefitinib as first‐line therapy was analyzed in surplus plasma samples. Plasma samples were collected at random timepoints during treatment during routine hospital visits and put into EDTA tubes and stored at (at least) −20°C until DNA isolation. In addition, erlotinib and/or gefitinib dosing information, plasma sampling time, and date and PFS details were collected retrospectively from medical records. The period from first dosage of erlotinib or gefitinib to the date of disease progression, based on computed tomography (CT) scan, was defined as PFS. For this retrospective observational study, no informed consent was required in accordance with the code of conduct for responsible use of human tissue and medical research.

Cell‐free DNA was isolated from 1 ml plasma for mutation analysis. The EGFR mutations L858R, exon19del, and T790M were a priori selected and hence quantified in purified DNA using ddPCR assays. Quantification of the ddPCR results was based on the quantity of positive droplets for mutation and wildtype, double positive droplets, and total accepted droplets in the assay. Allele concentrations were calculated by using the initial plasma sample volume from which the ctDNA was obtained and presented as copies per ml plasma.

ctDNA model development

Different dynamic models, such as baseline, turnover, and first‐order growth models, were tested to describe the longitudinal data of ctDNA concentrations of L858R, exon19del, and T790M. Initially, models for the three mutations were developed separately. In order to explore correlations among the mutations, the final model was a joint model including the complete dataset. Treatment with erlotinib or gefitinib was tested as a binary covariate in the final ctDNA model. We also used previously developed population PK models to derive individual parameter estimates based on relevant available covariates, as only dosing information was available for our patient population. These individual PK parameters were used to predict individual concentration‐time profiles of erlotinib or gefitinib and derive the area under the concentration‐time curve (AUC) as a measure of drug exposure. , The ctDNA concentrations were related to individual drug exposure, dose, and cumulative dose, to explore exposure‐response relationships. Linear, maximum effect (Emax), and sigmoid Emax relationships were investigated.

Interindividual variability (IIV) was evaluated for all parameters using an exponential error model (Equation 1): where the typical population parameter estimate and the individual parameter estimate for individual i are represented by P pop and P , respectively.

The IIV for subject individual i is represented by η , which was assumed to be normally distributed following N (0, ω2). Residual unexplained variability was described by a proportional error (Equation 2): where C obs, represents the observed concentration for individual i and observation j, C pred, represents the individual predicted concentration, and ε represents the proportional error distributed following N (0, σ 2). Separate residual unexplained variability was estimated for T790M concentrations and for L858R and exon19del concentrations combined, in order to account for observed differences in variability between the various measurements.

Survival model

To investigate the relationship between ctDNA concentrations and disease progression, a parametric time‐to‐event model was developed. First, the exponential, Gompertz and Weibull hazard functions were explored. From the final developed ctDNA model, individual ctDNA concentrations for the driver and T790M resistance mutations were predicted until the recorded time of progression/censoring. These predicted concentrations were then used to evaluate their value as predictors for disease progression. Absolute concentrations and relative change from baseline over time were tested for significance and biological plausibility. Harboring the exon19del or L858R mutation and treatment with erlotinib or gefitinib were tested as (binary) covariates in the survival model.

Model evaluation

Discrimination between models was guided by physiological and scientific plausibility, general goodness‐of‐fit (GOF) plots, precision of parameter estimates, and change in objective function value (dOFV). A drop of greater than or equal to 3.84, corresponding to a p < 0.05 (χ2‐distribution with 1 degree of freedom), was considered a significant improvement of the fit for hierarchical nested models. Kaplan–Meier plots were generated to evaluate the survival model. Prediction intervals were derived from 500 simulations and visually compared with the observed data. Parameter precision was assessed by the sampling importance resampling procedure.

Software

Nonlinear mixed‐effects modeling was performed using NONMEM (version 7.3; ICON Development Solutions) and Perl‐speaks‐NONMEM (PsN, version 4.4.8) with First‐Order Conditional Estimation with interaction as the estimation method for the ctDNA model. For the survival model, the Laplacian estimation method was used. Pirana (version 2.9.7) was used to compare and interpret model output. , , R (version 3.4.3) was used for data management and graphical diagnostics.

RESULTS

Patients and data

In total, 198 samples from 54 patients with NSCLC were available for this analysis. The L858R and exon19del driver mutations were detected in ctDNA of 19 (35.2%) and 35 (64.8%) patients, respectively. Observed T790M concentrations were lower compared with L858R and exon19del concentrations. In terms of treatment, 32 (59.3%) patients were treated with erlotinib, 16 (29.6%) with gefitinib, and 6 (11.1%) patients switched from erlotinib to gefitinib or vice versa.

ctDNA model

The observed L858R, exon19del, and T790M concentrations were best described by a zero‐order growth model. The data was best described when the L858R and exon19del concentrations were combined and thus described by shared parameters. , Equations 3 and 4 describe the model for the L858R and exon19del driver mutations:

where y(t) is the change in either L858R or exon19del over time, k in represents the zero‐order increase in ctDNA concentrations, k out represents the drug‐driven decrease in ctDNA concentrations, and R(t) accounts for the time‐dependent development of resistance where λ is the progression rate. The addition of this time‐dependent resistance term resulted in a significant increase of the model fit (dOFV = −32.5) and was thus included in the model. T790M concentrations were also best described by a zero‐order growth model, similar to Equation 3. In contrast to the time‐dependent resistance development of the driver mutations, the T790M concentrations were better described (dOFV = −21.4) by a model in which the resistance is depending on an increase in driver concentrations (y(t)) (Equation 5):

Plasma samples were collected after the start of treatment with erlotinib or gefitinib and, as a consequence, the baseline before start of treatment was unavailable. An additional parameter was estimated in the final model to account for the baseline concentrations. Parameter estimates are depicted in Table 1. High IIV was observed for all parameters for which IIV could be estimated (451%, 782%, 185%, and 1245% for driver baseline, driver progression factor, T790M baseline, and T790M progression factor, respectively). Mixture models and models incorporating eta‐transformations were tested in order to explore this high variability but were unable to improve the model. The GOF (Figure 1) and a selection of individual plots for three representative patients (Figure 2) show an adequate description of the L858R, exon19del, and T790M concentrations over time.

TABLE 1

Population estimates for the final ctDNA model

	Units	L858R and exon19del		T790M
		Estimate (95% CI a )	IIV b (95% CI a )	Estimate (95% CI a )	IIV b (95% CI a )
Population parameter
Baseline	Copies/ml	12.5 (7.73–18.6)	451% (279%–887%)	7.96 (5.76–10.8)	185% (139%–290%)
k in	Copies/ml/day	101 (76.7–121)	–	88.5 (70.2–109)	–
k out	/day	19.7 (14.9–23.6)		16.3 (12.9–20.0)
Progression rate (λ)	/day	5.0 × 10−4 (2.0 × 10−4–7.0 × 10−4)	782% (320%–986%)	1.6 × 10−3 (7.0 × 10−4–2.2 × 10−3)	1245% (324%–1542%)
Residual variability
Proportional residual error	CV%	30.9% (27.6%–34.9%)	–	80.5% (69.0%–94.5%)	–

Abbreviations: CI, confidence interval; CV%, coefficient of variation; IIV, interindividual variability.

The 95% CI values were obtained from sampling importance resampling.

IIV expressed as CV%, calculated as sqrt(exp[variance]‐1).

FIGURE 1

Goodness‐of‐fit plots of the final ctDNA model for the driver (blue dots) and T790M (red dots) concentrations. Including individual predictions (IPREDs) and population predictions (PREDs) vs. observed values and conditional weighted residuals (CWRES) vs. time after dose and PRED. Blue dots represent the driver concentrations and red dots are the T790M concentrations.

FIGURE 2

Individual plots for observed (solid line) and individual predicted (dashed line) driver (a) and T790M (b) concentrations over time for three representative patients (1–3). One patient with an increase in driver concentrations and stable T790M concentrations after ~ 1 year of treatment. Two patients with an initial decrease in both driver and T790M concentrations. Three patients with stable driver and T790M concentrations during ~ 1.5 years of treatment.

A continuous drug‐driven effect of erlotinib or gefitinib dosing or exposure (AUC) on the dynamics (k out) of the ctDNA concentrations could not be identified, most likely explained by the fact that all data was collected during treatment. In addition, no difference in parameters was identified for patients that were treated with gefitinib compared with patients that were treated with erlotinib.

Exponential, Gompertz, and Weibull hazard functions were explored in order to describe the disease progression data. The Weibull distribution with increasing hazard over time best described the baseline hazard for disease progression (dOFV = +5.55 and dOFV = +5.37 for the exponential and Gompertz models, respectively; Table 2). The visual predictive check of the survival model showed a good predictive performance of the model (Figure 3).

TABLE 2

Population estimates for the final survival model for progression‐free survival

	Units	Estimate (95% CI a )
Population parameter
Hazard coefficient, λ	/day	0.0013 (0.0011–0.0017)
Shape, α	–	1.3338 (1.0494–1.6130)
Covariate effect, β	–	0.0387 (0.0229–0.0503)

Abbreviation: CI, confidence interval.

The 95% CI values were obtained from sampling importance resampling.

FIGURE 3

Kaplan–Meier plot of progression‐free survival data. Blue line is the observed Kaplan–Meier curve and the shaded area represents the 95% prediction intervals (n = 500 simulations).

Relative change in driver concentration from the estimated baseline was a statistically significant predictor of disease progression, as shown by the diagnostic stratified Kaplan–Meier curve (Figure 4, stratified by the median relative change of driver concentration of 0.034) and was therefore included in the model (dOFV = −10.9, p = 0.001). Absolute driver and T790M concentrations did not show an improvement of the model fit compared to the relative change in driver concentrations, as well as the relative change in T790M concentrations. The final Weibull model is described by Equation 6:

FIGURE 4

Kaplan–Meier plot of observed progression‐free‐survival data stratified by median relative change in driver concentration from baseline at time of progression (p = 0.005, log rank).

where λ is the hazard coefficient and α is the shape parameter, DriverREL is the relative change in driver concentrations from baseline, and the covariate effect is represented by β. The final estimate for β was 0.0387 (95% confidence interval 0.0229–0.0503). No differences in PFS were observed for patients with the exon19del compared to the L858R mutation and for patients treated with erlotinib or gefitinib.

DISCUSSION

We successfully developed a nonlinear mixed effects model that describes the dynamics of L858R, exon19del, and T790M concentrations in plasma of patients with NSCLC who were treated with erlotinib or gefitinib. In addition, by using a parametric time‐to‐event model, the relative change in L858R and exon19del concentrations was identified as significant predictors of disease progression in this patient population. To the best of our knowledge, we are the first in presenting a nonlinear mixed effects biomarker model that describes ctDNA concentrations that were longitudinally collected from patients with NSCLC. To this end, we used a similar model as has been used previously to describe tumor size dynamics following anticancer treatment.

It has previously been shown that the trend of ctDNA concentrations measured in the plasma of patients correlated with the mutations detected in patient tumor tissue and the change in tumor volume as measured by CT scan. In this study, patients with ctDNA concentrations below the detection limit (i.e., 0 copies/ml) showed longer PFS and OS compared with patients with detectable ctDNA concentrations. Furthermore, both the EGFR driver and resistance mutation concentrations in plasma showed a progressive increase during disease progression. Measurements of ctDNA allowed earlier detection of NSCLC relapse compared to standard CT scan diagnostics with a median interval of 70 days (range 10–346 days). , In addition to an earlier detection of progression, measuring ctDNA concentrations in plasma samples is much less invasive for patients compared to conventional diagnostics, such as taking tumor tissue biopsies or making CT scans. This permits repeated determinations of ctDNA concentrations, thereby allowing for the monitoring of the change in these levels over time during treatment with EGFR inhibitors. Our survival model was significantly improved by the addition of the relative change in EGFR driver concentrations from baseline as a predictor for disease progression. This finding further confirms that a relationship exists between EGFR driver concentrations in ctDNA and progression of NSCLC in patients treated with the EGFR inhibitors erlotinib and gefitinib.

Our longitudinal model of driver and resistance mutations consisted of a zero‐order growth function (k in) in combination with an elimination function (k out) and resistance term (R(t)). The systemic measurements of ctDNA are derived from tumor cells, either in apoptosis or otherwise. Occurrence of driver and resistance mutations in plasma may be partly the result of the same underlying mechanism of systemic DNA shedding, probably depending on, for example, the level of drug exposure. For T790M, the resistance term appeared to be depending on an increase in driver mutation concentration. An abundance of driver mutation ctDNA might therefore reflect an overall increase in tumor load, as such leading to increased T790M shedding in circulation.

Nonlinear mixed effects modeling supports the longitudinal and mechanism‐based prediction of the full time course of biomarker concentrations in response to anticancer treatment. By coupling dynamic biomarker models to time‐to‐event models, the relationship between time‐varying biomarker concentrations and survival can be assessed. fIn addition, the effects of drug administration and thus development of resistance under drug treatment can be explored. This model can also be used for the analysis of ctDNA concentrations obtained from future clinical trials with targeted anticancer drugs.

Previous work has been focusing on fixed timepoint descriptors of drug exposure, biomarker concentrations, or tumor size as predictor for disease progression or OS. Modeling of continuous biomarker concentrations uses all available data and allows for simulations of the response to future treatment regimens. This has also been shown by a PK/PD model that describes the tumor size dynamics in patients with gastrointestinal stromal tumors. This model describes OS and tumor size dynamics as a function of the relative change in stem cell factor receptor, daily sunitinib exposure and relative change in the soluble vascular endothelial growth factor receptor 3.

Our study does, however, harbor some limitations. As this was a retrospective analysis of data that was collected in a cohort study from routine clinical care patients, only limited data were available. The dynamics of the ctDNA concentrations, in particular, the time‐dependent development of resistance during treatment, could be captured by our model but was driven by a limited number of patients in our dataset. As a result, very high IIV was estimated for the progression factors. Mixture models and eta‐transformations were tested but could not improve the eta‐distributions. In addition, patient samples were collected as part of routine care. The number of samples collected directly after the start of treatment was therefore limited. As a result, the final model estimates for k in and k out were overestimated, complicating simulations based on this model early after initiation of treatment. In addition, repeated tumor size and PK measurements were unavailable. Therefore, we could not relate the model for ctDNA to tumor dynamics or assess the effects of PK on the ctDNA dynamics. It has, however, been shown that ctDNA concentrations in plasma correlated with tumor size measurements in patients with NSCLC. Therefore, our model probably provides a good approximation of the tumor size dynamics in this patient population. We used previously developed PK models for erlotinib and gefitinib and individual dosing records to incorporate the effect of PK on the ctDNA dynamics but this did not improve the model. Most likely as a result of the very small dosing range and patients being at steady‐state.

We analyzed the combined concentrations of patients harboring the exon19del or L858R mutations. Separate modeling showed similar parameter estimates for both mutations, but with poor parameter precision. In order to improve the power of this analysis, we decided to describe both driver concentrations by the same dynamics and parameters. It has previously been suggested that patients with the exon19del mutation have a higher response rate and longer median OS compared with patients with the L858R mutation. We performed a covariate analysis to explore this relationship, but could not identify a difference between these groups in our data. In addition, no differences were found in either the ctDNA model or the survival model for patients who were treated with erlotinib or gefitinib. These results are in line with previous studies that compared erlotinib with gefitinib for differences in PFS, OS and biomarker tumor size dynamics. , ,

CONCLUSION

In conclusion, we successfully developed a model that describes the longitudinal dynamics in ctDNA in a clinical cohort and identified the relative change in EGFR driver mutations from baseline as a predictor for disease progression following erlotinib or gefitinib treatment. In addition, this model shows that nonlinear mixed effects modeling is a valuable method for the analysis of longitudinal and heterogeneous biomarker datasets obtained from clinical practice and provides a framework for prediction of disease progression based on observed ctDNA concentrations.

CONFLICT OF INTEREST

As of November 2017, R.B.V. is an employee of AstraZeneca Cambridge, UK. N.S. provided consultation or attended advisory boards for Boehringer Ingelheim, Ellipses Pharma, and received research grants for the institute from AB Science, Abbvie, Actuate Therapeutics, ADCtherapeutics, Amgen, Array, Ascendis Pharma, Astex, AstraZeneca, Bayer, Blueprint Medicines, Boehringer Ingelheim, BridgeBio, Bristol‐Myers Squibb, Cantargia, Celgene, CellCentric, Cresecendo, Cytovation, Deciphera, Eli Lilly, Exelixis, Genentech, Genmab, Gilead, GlaxoSmithKline, Incyte, InteRNA, Janssen/Johnson&Johnson, Kinate, Merck, Merck Sharp & Dohme, Merus, Molecular Partners, Novartis, Numab, Pfizer, Pierre Fabre, Regeneron, Roche, Sanofi, Seattle Genetics, Servier, Taiho, and Takeda, outside of the submitted work. J.H.B. is a part‐time employee and (indirect) shareholder of Modra Pharmaceuticals BV. He is (partly) a patent holder of oral taxane formulations, which are clinically developed by Modra Pharmaceuticals BV (a spin‐off company of the Netherlands Cancer Institute, not related to this work). All other authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

J.M.J., L.L., M.M.H., J.H.B., N.S., D.B., A.D.R.H., and T.P.C.D. wrote the manuscript. A.D.R.H. and T.P.C.D. designed the research. J.M.J., R.B.V., T.T.D., D.B., and T.P.C.D. performed the research. J.M.J. analyzed the data. J.H.B. contributed new reagents/analytical tools.