Comparing the Efficacy of Two Generations of EGFR-TKIs: An Integrated Drug-Disease Mechanistic Model Approach in Egfr-Mutant Lung Adenocarcinoma

1. Introduction

Lung cancer remains the leading cause of cancer-related deaths globally, with over 1.8 million deaths and 2.2 million new cases yearly [1]. The high mortality rate is largely attributed to late-stage diagnosis despite advances in screening [1]. Non-small cell lung cancer (NSCLC) constitutes 80% of cases, with about 40% being lung adenocarcinoma (LUAD) [2]. Driver mutations, like the mutation of the epidermal growth factor receptor (EGFR) gene, prevalent in 48% of Asian and 19% of Western LUAD patients, play a crucial role in tumor growth and invasiveness [3]. EGFR mutations, commonly either the exon 21 L858R insertion or deletion on exon 19, trigger uncontrolled activation of the EGFR receptor, promoting cancer cell growth and survival [4]. This understanding has led to targeted therapies focusing on these driver mutations, offering an alternative to traditional chemotherapy. Targeted therapies specifically target driver mutations and disrupt their ability to drive tumor proliferation and metastasis’ apparition. Since 2015, EGFR tyrosine kinase inhibitors, which compete with ATP at the intracellular kinase domain of EGFR, inhibiting its auto-phosphorylation and constitutive downstream pathways’ activation, have been approved by the Food and Drug Administration (FDA) to treat LUAD patients with an EGFR driver mutation. The first generation TKIs, namely erlotinib and gefitinib, are reversible inhibitors meaning that once the drug unbinds its target, the effects of the drug disappear, in contrary to second generation TKIs afatinib and dacomitinib that are irreversible inhibitors. Osimertinib, a third generation TKI, was later developed to overcome resistance to 1st and 2nd generation EGFR TKIs, as it targets both common EGFR mutations and the T790M resistance mutation, an EGFR mutation on exon 20 that frequently appears during treatment with 1st or 2nd generation TKIs and that causes resistance to these treatments [5] . Despite these advancements in targeted therapies, predicting patient response to specific EGFR-TKIs remains a challenge , especially because resistance to treatment is almost inevitable [6]. This lack of response predictability underscores the need for more refined decision-making tools in clinical practice. Here, we propose the application of a Knowledge-Based Model (KBM) as a potential solution [7]. The KBM can simulate the responses of a specific patient population to different TKIs targeting the same receptor, thereby enabling comparative analysis of various treatment strategies and sequences and offering insights into optimal patient-specific therapeutic approaches. To achieve this, we refined a KBM named ISELA-V1, that we previously developed [8] and which mechanistically represents the main biological phenomena related to EGFR-mutant LUAD disease, and combined it with physiologically-based pharmacokinetics models of osimertinib and gefitinib that incorporate data on these two EGFR TKIs and their mode of action. Using this updated model, named ISELA-V2, we simulated responses to both EGFR TKI in the same patient population and compared the outcomes and efficacy of both treatments. This article details the methodology behind this approach and demonstrates how it can be used to inform treatment choices, filling the gaps in the current understanding of differential treatment responses and potentially guiding more effective clinical decision-making for LUAD patients.

2. Materials and Methods

2.1. Construction of the ISELA-V2 Model

In order to evaluate the impact of a treatment on LUAD evolution, we first need a model of this disease compatible with treatment models. To develop this disease model, we used the ISELA-V1 model presented by L’Hostis et al. [8] and whose credibility has already been assessed. This disease model was then improved in order to represent in more detail biological phenomena that are essential to better grasp the complexity of tumoral characteristics, namely the neoangiogenesis, the cell cycle and death, and the formation of metastases. On top of this disease model, we added PBPK treatment models combined with the implementation of their mechanism of action for the two EGFR-TKI gefitinib and osimertinib. All these additions are presented in the supplementary materials 1. The combination of the improved disease model with gefitinib and osimertinib treatment options has been named ISELA-V2. Its outputs are notably composed of the time to progression as defined by the RECIST criteria version 1.1 [9].

2.2. PBPK Model of Osimertinib and Gefitinib

In order to have access to the concentration of the TKI present in the tumor tissues that will contribute to reducing growth signals, we developed two PBPK models, one for gefitinib and one for osimertinib. PBPK models rely on the known physiology to predict how a drug will be absorbed, distributed, metabolized and eliminated (ADME) by the body depending on its physico-chemical properties. Those models are particularly useful in our context because they estimate the concentration of the drug in all the main organs of the body and therefore in all the organs that carry the primary tumor or a metastasis. Moreover, even if their structure is complex, it stays the same between two drugs that are similar in terms of physicochemical properties. For this reason, the PBPK models of gefitinib and osimertinib were built using the same equations. Nonetheless, a particularity of osimertinib is that one of its metabolites, AZ5104, also acts as an EGFR TKI [10]. The metabolization of osimertinib to AZ5104 in the liver has therefore been modeled, as well as the equations of distribution of AZ5104 from the liver to the different tissues and tumors. The structure of the PBPK models is presented in more detail in the supplementary materials (Section 2). The equations of the PBPK models regulate the pharmacokinetics of the drugs and grant access to the concentration of the drug in each of the main organs of the body. For the distribution of the drug in the tumors (primary tumor or metastasis), it was assumed that the concentration of drug in the tumor is directly proportional to its concentration in the organ the tumor is located in. Because the quantitative public data on osimertinib and gefitinib mechanism of action were found mainly in mice, a mice setting was included in the model, keeping all the equations but changing all the physiological parameters to correspond to a mouse. The parameters to calibrate related to the pharmacokinetics of the drugs were calibrated separately in mice and humans to reproduce the plasmatic profile of the drug after administration in both species.

2.3. Modeling TKI Mechanism of Action

The disease model has been developed in the context of EGFR mutated LUAD, and the modeling of the growth factors signaling, thus including the signals arising from EGFR activation, follows the activation pathways depicted in Figure 1.

Once activated, the phosphorylation of EGFR activates MAPK and PI3K/AKT downstream pathways, eventually leading to cell proliferation and survival. Not all intermediaries have been modeled, the first modeled actors being RAS and PI3K. Of note, there exist several other tyrosine kinase receptors that would activate the same downstream pathways as EGFR [11]. Among them, we chose to model MET and VEGFR receptors, as they would play specific roles in phenomena of interest for the context of use of the model: the resistance of a tumor to EGFR TKI treatment can be caused by an amplification of the MET receptor [6] and VEGFR is a receptor playing a role in the neoangiogenesis of the tumor [12]. The activation of RAS and PI3K by EGFR, VEGFR and MET are respectively depicted in Equation (1) and (2).

$$\begin{array}{cc}\hfill \frac{dTumorCell.RasActive}{dt}=\frac{TumorCell.RasInactive}{vol\left(TumorCell\right)}& (kaRasByMet\times TumorCell.MetActive\hfill \\ \hfill & +kaRasByVegfr\times TumorCell.VegfrActive\hfill \\ \hfill & +kaRasByEgfr\times TumorCell.EgfrActive\times propPhosphoEgfrActive)\hfill \\ \hfill & -TumorCell.RasActive\times kdaRas\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{dTumorCell.Pi3kActive}{dt}=\frac{TumorCell.Pi3kInactive}{vol\left(TumorCell\right)}\times & (kaPi3kByMet\times TumorCell.MetActive\hfill \\ \hfill & +kaPi3kByVegfr\times TumorCell.VegfrActive\hfill \\ \hfill & +kaPi3kByEgfr\times TumorCell.EgfrActive\times propPhosphoEgfrActive)\hfill \\ \hfill & +kaPi3kByRas\times TumorCell.RasActive)\hfill \\ \hfill & -TumorCell.Pi3kActive\times kdaPi3k\hfill \end{array}$$

(2)

The model represents the different states of EGFR: Inactivated, activated and phosphorylated: The parameter propPhosphoEgfrActive corresponds to the proportion of activated EGFR that are phosphorylated. It multiplies the concentration of activated EGFR to model the concentration of EGFR activating downstream pathways. It includes the mechanism of action of tyrosine kinase inhibitors that will impact the strength of the signal by competing with ATP in the intracellular domain and preventing the phosphorylation of the receptor. It is defined by:

$$propPhosphoEgfrActive=\frac{\left[ATP\right]}{\left[ATP\right]+atpAffinity\times (1+\frac{1+\left[TKI\right]}{K{i}_{TKI}})}$$

(3)

The value of the ATP affinity with EGFR intracellular domain can vary among the existing EGFR mutations, as well as the inhibition coefficient of the tyrosine kinase inhibitor [13,14,15]. The values of the inhibition coefficient are given in Table 1:

As expected, the T790M mutation that causes resistance to gefitinib treatment is linked to a high $K{i}_{TKI}/atpAffinity$ ratio for gefitinib compared to other sensitive mutations. The exon 20 resistant mutation that causes resistance to both gefitinib and osimertinib also presents a high $K{i}_{TKI}/atpAffinity$ ratio for gefitinib and osimertinib.

2.4. Mechanisms of Resistance

The two EGFR TKI models built are using the same model structure, whether it be for the mechanism of action of the drug or the pharmacokinetic model. Some elements still allow distinguishing those two models, which are the inhibition constant for ATP binding in the tumor site (Table 1) and the PK data the models are calibrated on (Section 2.4). But what is more important to differentiate osimertinib and gefitinib is the mechanisms that will grant the tumor resistance to the treatment. In most cases, the development of a resistance mechanism is marked by a regrow of the tumor after a phase of decline. Whether those mechanisms of resistance appear during the treatment or are already in place but only for a limited number of cells which are then selected by the treatment is uncertain and may depend on each patient and each mechanism of resistance [16,17]. In our approach, for simplification purposes, we chose the second option and assumed that a limited number of cell bearing resistance mechanisms are already present in the tumor when starting the treatment. Gefitinib and osimertinib, even though they have an identical mechanism of action, have similar but also distinct mechanisms of resistance [18]. What also differentiates the two treatments is the frequency of apparition of the mechanisms of resistance, which can also be impacted by whether the treatment is given in first or second line [19,20]. For each treatment, we have implemented the most prevalent associated mechanisms of resistance namely EGFR alteration (amplification, T790M mutation loss and C797S mutation), MET amplification, KRAS mutation or amplification and PIK3CA mutation or amplification for osimertinib [19] and T790M EGFR mutation and MET amplification for gefitinib (listed in decreasing order of frequency) [21]. Note that some mechanisms of resistance to osimertinib would also provide resistance to gefitinib treatment but they were not implemented in the model as their prevalence is negligible compared to the other resistance mechanisms to gefitinib.

In order to be able to represent both the tumor heterogeneity and the dynamic of resistance to a treatment, each tumor is divided into separate subclones that are distinct subpopulations of tumor cells each carrying their own set of mutations. Among those subclones is a resistant subclone that can carry the resistance mechanisms described above. This subclone is initialized with a negligible size compared to other subclones but if it carries a mechanism of resistance, it will grow to the point that it will predominate over the other subclones, which in contrast are declining in response to the treatment. To reproduce the variability of response to EGFR TKI, more than 70 patient descriptors were introduced with the complete list available in the supplementary materials (spreadsheet). As an example, regarding the resistance to treatment, for each virtual patient the presence or absence of each resistance mechanism in the resistant subclone is drawn from a Bernoulli distribution, whose parameter has been informed by the literature [19,21,22,23]. The presence of a resistance mechanism in the primary tumor does not necessarily mean that all the metastases will carry this resistance. Therefore, we introduced a parameter that either ensures the same mechanism of resistance in all the tumors and metastases (with a 70% probability) or allows for independent mechanisms of resistance between the primary tumor and the metastases (with a 30% probability). The assumption that in 70% of EGFR mutated LUAD cases, the mechanisms of resistance are the same in tumor and metastases has been made based on discussions with clinicians.

2.5. Data for Model Calibration

The Table 2 describes the dataset used to inform the model on the drug pharmacokinetics and mechanism of action. The methodology used to optimize the parameters of the model using various datasets has been described by Palgen et al. [24].

2.6. Comparison of Model Prediction with a Retrospective Study

The retrospective study presented by Li et al. [33] compares the efficacy of osimertinib and gefitinib and was not used for calibration purposes. Therefore, we were able to use this study to confront the model predictions with real-life data . This study included 102 patients with stage III-B or IV NSCLC, of whom 49 were treated with osimertinib and 53 with gefitinib. The PFS was measured and will serve as a comparative endpoint for our model. Table 3 summarizes the characteristics of the patients that participated in this study and that were used to create the associated virtual population.

2.7. Effect Model

Lastly, we performed an in silico clinical trial to compare the efficacy of gefitinib and osimertinib in the same virtual population. By doing so, we can assess for each patient his time to progression whether he is treated with osimertinib or with gefitinib. The population chosen to perform this trial is summarized in Table 4 and contains the default distribution for each of the patient’s descriptors to represent a generic EGFR mutated LUAD population. The number of virtual patients has been set to 10 times the number of patients in the retrospective study in order to have enough variability in the virtual population as well as enough statistical units to use bootstrapped approaches. This method has been described by Boissel et al. as a way to give insights on treatment efficacy and best responders characteristics [34].

3. Results

3.1. Reproducing Pharmacokinetic Data of Gefitinib and Osimertinib

After calibration, the PBPK models of gefitinib and osimertinib were able to reproduce the plasmatic concentrations of the drugs both in mice and in humans (Figure 2). All the simulated data points are within the range defined by the standard deviation.

3.2. Reproducing the Mechanism of Action of Gefitinib and Osimertinib

After calibration, the model was also able to reproduce the tumor volume evolution after gefitinib or osimertinib treatment in mice transplanted with xenografts bearing an EGFR exon 19 deletion (Figure 3) but also with xenograft bearing both an EGFR exon 19 deletion and a PIK3CA mutation (Palgen et al. [24]). The model also reproduced the dynamics of the tumor without treatment (Palgen et al. [24]), ensuring that the treatment is indeed responsible for the shrinkage of the tumor when administered. The Table 5 summarizes the calibration results in mice.

Using population characteristics composed of generic information such as age and gender as well as disease specific information such as initial cancer stage and EGFR mutation, the model was able to reproduce the probability of tumor progression as measured in the FLAURA and NEJ002 clinical trials. In order to achieve these results, several parameters of the model that could not be informed by in vitro experiments or via animal testing were calibrated, such as the initial size of the micrometastases, the initial size of the resistant subclone, the reduced impact of the immune system in the metastases compared to the primary tumor [35] and the scaling factor regulating the distribution of the treatments in the tumors. To compare the goodness of fit of the simulated Kaplan-Meier curve with the one from the associated clinical trial, bootstrapped combined weighted log-rank tests were performed using 5000 bootstraps with subsamples containing a tenth of the global population. This statistical methodology has been described by Jacob et al. [32]. The threshold of the ratio of non-significant tests required to consider that the two curves are similar was set to 80% and was largely exceeded for the two simulations (99.86% for the population treated with osimertinib and 98.42% for the population treated with gefitinib).

3.3. Reproducing the Results from a Retrospective Study

To compare the efficacy of gefitinib and osimertinib in treating EGFR mutated LUAD, we simulated a clinical trial with virtual patients representing the participants of the retrospective study from Li et al. [33]. In contrast to the previous results, no parameter was changed to fit the outcome of the retrospective study, only the patients’ characteristics were updated to match the ones of the real population described in this study. As presented in Table 6, the results of the in silico clinical trial show that osimertinib is associated with longer TTP than gefitinib for the designed population with a median TTP of 20 months (95%CI: 15-24) for osimertinib and 11 months (95% CI: 7.5-12) for gefitinib. The retrospective study reported a median PFS of 18.1 months (95% CI: 15.4-20.7) for osimertinib and 10.7 months (95% CI: 9.9-11.4) for gefitinib, in coherence with the TTP results of our in silico clinical trial, although the endpoints are not identical as the model predicts the TTP while the retrospective study reports the PFS. This highlights the robustness of the model to predict the efficacy of osimertinib and gefitinib in EGFR mutated LUAD patients.

Figure 5. Exploratory in silico clinical trial: comparison of gefitinib and osimertinib performed with the LUAD model and with the same patient baseline characteristics as the retrospective study from Li et al. The uncertainty intervals stand for the 95% confidence intervals of the progression curves. The dotted lines represent the median TTP. Patients were orally treated with 250mg of gefitinib daily in the gefitinib arm (in orange) and with 80 mg of osimertinib daily in the osimertinib arm (in blue). Both treatments are starting at t=0.

Figure 5. Exploratory in silico clinical trial: comparison of gefitinib and osimertinib performed with the LUAD model and with the same patient baseline characteristics as the retrospective study from Li et al. The uncertainty intervals stand for the 95% confidence intervals of the progression curves. The dotted lines represent the median TTP. Patients were orally treated with 250mg of gefitinib daily in the gefitinib arm (in orange) and with 80 mg of osimertinib daily in the osimertinib arm (in blue). Both treatments are starting at t=0.

3.4. Effect Model

By simulating the time to progression for each patient of the virtual population whether they are treated with osimertinib or gefitinib, we can directly compare the efficacy of the two treatments by plotting an effect model (Figure 6). This effect model suggests that treatment with osimertinib is linked to a longer TTP for the majority of patients. This result is even more accentuated for patients carrying an exon 20 mutation, which should however be interpreted with caution as the data that was used to calibrate the model did not contain any patients carrying exon 20 mutations. Only patients who progressed when treated with both gefitinib and osimertinib were included in the effect model but it is interesting to note that out of the 1020 virtual patients, 329 did not progress during the 24 months of simulation when treated with osimertinib while there were only 119 for gefitinib. This is an additional prediction towards the superior efficacy of osimertinib.

4. Discussion and Conclusions

The success of the calibration to reproduce the pharmacokinetics and drug effect of gefitinib and osimertinib in mice and humans constitutes evidence of the credibility of the model [36]. The model credibility is also increased by the reproduction of the retrospective study, containing patient data not used during the calibration process, using only the patient’s characteristics at the population level. This constitutes evidence that knowledge based mechanistic models are able to grasp the complexity of biological phenomena to predict the result of experiments or clinical trials by using previously gathered knowledge and data. Such models could play a major role in the development of new therapies [37]. They could be used to identify patients that would respond the best to a given treatment by analyzing the correlations of patient’s characteristics with treatment outcome. This would help guide the recruitment of patients for clinical trials, a recruitment whose quality can dictate the outcome of the trial and that can be challenging depending on the rarity of the disease [38]. Another attractive aspect of such models is the ability to generate control arms for clinical trials. Indeed, the model we developed was able to predict the response to osimertinib and gefitinib, two standards of care, in an EGFR mutated NSCLC population. By generating virtual patients representative of the target population, the model could create a control arm with simulated patients receiving either osimertinib or gefitinib. This virtual control arm could then be compared to an arm of real patients receiving an investigational treatment, allowing for a more accurate assessment of the investigational treatment’s efficacy. While the model introduced in this article holds promise for EGFR mutated LUAD patients, its applicability beyond this specific context remains limited. Any simulations performed outside of this intended use, while potentially yielding intriguing perspectives, should be approached with caution due to the inherent credibility gap compared to predictions within its domain of application [36]. In addition, the model still presents some shortcomings that could be solved with more development time. As an example, the PBPK models are based on the physiology of healthy patients while it is known that patients carrying a tumor can suffer from many symptoms such as renal impairment that can affect the pharmacokinetics of the drug and therefore its effect [39]. Another limiting point is that the model is not able to predict death, as the exact mechanisms leading to the death of a patient during the follow-up period are unknown and depend on a variety of factors that are difficult to assess (adverse effect of a treatment, exhaustion, metastases). The model only reproduces the time to progression of patients, which has to be compared to the more standard clinical endpoints that are the overall survival and the progression free survival with caution. Nevertheless, the model presented here is on the path to become an invaluable asset in the search of new therapeutic strategies in EGFR mutated lung cancer and has even recently been used to predict with success the results of a Phase III clinical trial whose results were not published yet, a first in the in silico industry [40,41].