Mathematical analysis identifies the optimal treatment strategy for epidermal growth factor receptor-mutated non-small cell lung cancer

Introduction In Asians, more than half of non-small cell lung cancers (NSCLC) are induced by epidermal growth factor receptor (EGFR) mutations. Although patients carrying EGFR driver mutations display a good initial response to EGFR-Tyrosine Kinase Inhibitors (EGFR-TKIs), additional mutations provoke drug resistance. Hence, predicting tumor dynamics before treatment initiation and formulating a reasonable treatment schedule is an urgent challenge. Methods To overcome this problem, we constructed a mathematical model based on clinical observations and investigated the optimal schedules for EGFR-TKI therapy. Results Based on published data on cell growth rates under different drugs, we found that using osimertinib that are efficient for secondary resistant cells as the first-line drug is beneficial in monotherapy, which is consistent with published clinical statistical data. Moreover, we identified the existence of a suitable drug-switching time; that is, changing drugs too early or too late was not helpful. Furthermore, we demonstrate that osimertinib combined with erlotinib or gefitinib as first-line treatment, has the potential for clinical application. Finally, we examined the relationship between the initial ratio of resistant cells and final cell number under different treatment conditions, and summarized it into a therapy suggestion map. By performing parameter sensitivity analysis, we identified the condition where　osimertinib-first therapy was recommended as the optimal treatment option. Discussion This study for the first time theoretically showed the optimal treatment strategies based on the known information in NSCLC. Our framework can be applied to other types of cancer in the future.


Introduction
Among all the cancer types, lung cancer causes the highest number of cancer-related deaths.26% of cancer-related deaths in males and 25% in females are induced by lung cancer (1).In Asians, 53% of non-small cell lung cancer (NSCLC) progression is induced by epidermal growth factor receptor (EGFR) mutations, such as the L858R mutation, exon 19 deletion, and exon 20 insertion (2).Besides, EGFR has been recognized as an oncogenic driver of NSCLC, making it increasingly important in the era of precision medicine for lung cancer (3).
EGFR belongs to the receptor tyrosine kinase (RTK) family.EGFR is activated by various ligands in the extracellular environment and transmits cellular responses to mediate many cellular activities, including cell proliferation, survival, growth, and development.It is expressed in many organs, with its abnormal expression associated with a variety of cancers.EGFR has an extracellular ligand-binding domain, hydrophobic transmembrane domain, and cytoplasmic tyrosine kinase domain.The driver mutations in EGFR associated with cancers are concentrated in the tyrosine kinase domain, forming exons 18-21 (4)(5)(6)(7).More than 200 types of EGFR mutations have been identified, but the most common types are exon-19 deletion and the L858R mutation in exon 21 (8,9).Approximately 44% of EGFR-mutated patients harbor exon-19 deletion, and 31% have the L858R mutation (10).
Although EGFR was first identified in 1977, EGFR-targeted antitumor drugs were first reported in 1994 (11).After the first report of EGFR-targeted therapy, first-generation EGFR-Tyrosine Kinase Inhibitors (EGFR-TKIs) were not approved until 2004 (12).Subsequently, the second-generation EGFR-TKI, afatinib, was approved in 2014.First-and second-generation EGFR-TKIs are effective in most cases of lung cancer harboring EGFR driver mutations (13)(14)(15)(16).However, acquiring mutations, such as the T790M mutation, causes drug resistance and induces treatment failure (17,18).In 2015, the third-generation EGFR-TKI (osimertinib), which inhibits both driver mutations and the T790M mutation, was approved as a second-line drug for patients with EGFR mutations (19)(20)(21)(22).Although osimertinib is clinically effective, the emergence of additional mutations, such as the C797S mutation, induces resistance to osimertinib (23)(24)(25).Clinical observations suggest that optimized treatment schedules can help patients achieve better therapeutic effects (26)(27)(28).Thus, predicting resistance evolution and making reasonable treatment schedules in advance are necessary to delay the appearance of drug resistance and prolong the time of recurrence.However, even with knowledge of medical and genetic information in the early stage, such as tumor size and the proportion of different genotypes, it is still difficult to simulate the future development of tumors using traditional biological techniques alone.
Mathematical modeling is an approach for simulating realistic problems using mathematics and computational algorithms.This can offer a better understanding of how tumors evolve during treatment, which can be visualized in vivo.Thus, it can help us predict tumor dynamics under certain treatment schedules, compare different treatments, and even suggest optimal treatment strategies.Many studies have demonstrated the capability of mathematical modeling in cancer-related research (29)(30)(31)(32)(33).For example, Diaz Jr. et al. developed a mathematical model of cell kinetics during chemotherapy to predict the emergence of resistant genotypes in colorectal cancer (30).Castagnino et al. established a mathematical model of a genetic network to identify novel molecular targets for the treatment of colorectal cancer (34).In this way, we decided to employ mathematical modeling to predict tumor evolution and direct reasonable treatment schedules for lung cancer patients harboring EGFR mutations.
In this study, we establish a novel mathematical model of lung cancer evolution under EGFR-targeted therapy based on clinical observations.Parameter values in the model are estimated from published literature, and the results are confirmed using clinical observations.Moreover, we examine the relationship between the timing of switching drugs and the final number of cells in the tumor.Furthermore, we compare the combinatorial use of EGFR-TKIs to their sequential use to test their potential for clinical application.Finally, we investigate how intratumoral heterogeneity at the initial time of therapy affects treatment outcomes.The simulation results are comprehensively tested by parameter sensitivity analysis in order to identify the condition where each treatment strategy becomes the best option.Our framework is expected to be capable of suggesting reasonable individualized medicine for EGFR-mutated NSCLC.

Mathematical model
Based on clinical observations (35, 36), we established a mathematical model that describes the dynamics of the four types of EGFR-mutated cells under two types of EGFR-TKIs (Figure 1).There are two different types of EGFR-TKIs in the model: one is "DrugA," representing the firstor second-generation EGFR-TKIs named gefitinib, erlotinib or afatinib; and the other is "DrugB," representing osimertinib.Four cancer cell types are denoted by "Type-W," "Type-X," "Type-Y," and "Type-Z".Type-W is sensitive to both DrugA and DrugB, indicating cancer cells with driver EGFR mutations, such as L858R mutations or exon-19 deletion.Type-X cells are resistant to DrugA but sensitive to DrugB, indicating cells with T790M mutations.Type-Y is sensitive to DrugA but resistant to DrugB, indicating cells with C797S mutations.Type-Z is resistant to both DrugA and DrugB.Summarizing the relationship between drugs and cells, under DrugA treatment, Type-W and Type-Y will decrease, but Type-X and Type-Z will increase, whereas under DrugB treatment, Type-W and Type-X will decrease, but Type-Y and Type-Z will increase.According to published clinical studies (37-39), when DrugA was administered as first-line treatment, Type-X eventually became dominant in the tumor.After switching from DrugA to DrugB, the frequency of Type-X decreased, and only Type-Z continued to grow and dominate the tumor.However, when using DrugB as the first-line treatment, Type-Y will replace Type-W as the major population.After switching from DrugB to DrugA, only Type-Z became the donor population in the tumor.
In this mathematical model, we assumed that each cell type itself increases in number by cell division and mutates into a resistant type at a low mutation rate.We did not consider back mutations that were resistant to the sensitive cells.Moreover, according to the purpose of this study, we only focused on mutations related to drug resistance and assumed that other mutations are neutral and do not affect the growth kinetics.Then the dynamics of Type-W, Type-X, Type-Y, and Type-Z are given by Eqs.
Here, the variables w, x, y, and z represent the cell numbers of Type-W, -X, -Y, and -Z, respectively.Parameters a, b, c, and f are the growth rates of Type-W, Type-X, Type-Y, and Type-Z, respectively, and g, h, k, p, and q are the mutation rates from type-W to Type-X, Type-W to Type-Y, Type-W to Type-Z, Type-X to Type-Z, and Type-Y to Type-Z, respectively.Because no other cell type can mutate into Type-W, the number of Type-W cells is affected by its kinetics only.However, Type-W will mutate into Type-X, Type-Y, and Type-Z.

Solution of equations
The Eqs. (1) can be solved and given by Eqs. ( 2) (2:1) W 0 e (at−bt)e bt (2:2) The equations describe the cell number of each type over time (t) during therapy.W 0 , X 0 , Y 0 , and Z 0 represent the initial cell numbers of Type-W, -X, -Y, and Type-Z in the tumor.Please refer to Table 1 for the meaning of each letter in the model.

Parameter evaluation
The parameter values were obtained from the published data (Table 2) (40)(41)(42).Since we obtained growth parameters under erlotinib and osimertinib treatments, we regarded these drugs as representative of DrugA and DrugB, respectively.Since Starrett et al. (41) reported that first-line therapy with erlotinib and osimertinib delayed the emergence of secondary mutations in untreated EGFR-mutated NSCLC, thus, for combination therapy, we defined a combinatorial regimen of erlotinib plus osimertinib as DrugC.Based on genome-editing cell line experiments (40), we adopted the growth rate of EGFR-L858R mutated cells for Type-W as -0.17 [/day] under DrugA (a A ) and -0.32 [/day] under DrugB (a B ).Note that the subscript of each growth rate represents the condition of drugs, i.e., a A represents the growth rate of Type-W under DrugA.From Starrett et al. (41), we adopted the growth rate of EGFR-L858R/C797S mutated cells for Type-Y as -0.13 [/day] under DrugA (c A ) and 0.024 [/day] under DrugB (c B ).The growth Illustration of the Model.Blue cells represent EGFR-TKI sensitive genotypes (for example, EGFR L858R and EGFR del-19 ), orange cells represent osimertinib sensitive genotypes (such as EGFR L858R-T790M or EGFR del-19-T790M ), yellow cells represent osimertinib-only resistant genotypes (for instance, EGFR L858R-C797S or EGFR del-19-C797S ), green cells represent all EGFR-TKI resistant genotypes (like EGFR L858R-T790M-C797S and EGFR del-19- T790M-C797S ).DrugA involved the firstand second-generation EGFR-TKIs (erlotinib, gefitinib, and afatinib), and DrugB is osimertinib.
Based on the adopted parameters, we assumed all the other parameter values.Because the growth rate of Type-Y under DrugC is approximately 26% of that under DrugA (-0.0335/-0.13),we calculated the growth rate of Type-W in DrugC as 26% of that under DrugA, which is -0.064 [/day] (a C ).We assume the growth rate of Type-X under DrugC is same as Type-Y, which is -0.0335 [/day] (b C ).Moreover, we assumed the growth rate of Type-X under DrugB as -0.15 [/day] (b B ), which is smaller than that of Type-Y under DrugA (c A ) based on clinical observation (20,21) where the first-line treatment by DrugB showed better prognosis than that by DrugA.Based on the same reason, we assumed the growth rate of Type-X cell under DrugA as 0.045 [/day] (b A ). Finally, since Type-Z is resistant to both DrugA and DrugB, we assume its growth rates under DrugA and DrugC are same as that under DrugB effect, which is -0.022 [/day] (f A and f C ).
As for the initial condition of simulations, the initial total cell number of the tumor is set to be 10 9 , and the standard initial cell number of Type-X (X 0 ), Type-Y (Y 0 ), Type-Z (Z 0 ), and Type-W (W 0 ), is 10 4 , 10 4 , 10, and the rest component, respectively.The initial total cell number is set to be 10 9 because the diameter of a tumor at this point is about 1cm and a detectable size clinically.About the drug switching time in monotherapy, we set day-307 (sta=307) under DrugA-first therapy and day-567 (stb=567) under DrugB-first therapy based on clinical statistic data of the median Progression-Free Survival (mPFS) (20).The whole treatment time is assumed to be 1000 days (T=1000) in our simulation because 1000 days is long enough to compare the treatment options.

Computational simulation
We used Python (version 3.8.8) to simulate our model.We did time course simulation of different therapies for checking whether our model can express the progression of tumor as clinical observation.Then, we simulated the relationship between drug switching time and final cell number to theoretically figured out the possible affects that could influence therapy effects.Moreover, we simulated how parameters affected final cell number in different therapies.All the codes for simulations can be found in our GitHub open repository: https://github.com/yuqianxigua/EGFR-TKIs-therapy.

Tumor evolution under different treatments
We simulated Eqs.(2) to predict tumor progression under different treatments, including monotherapy and combination therapy.When DrugA was used as first-line treatment (Figures 2A, D), Type-W and Type-Y decreased, whereas Type-X and Type-Z increased.Once Type-X became the dominant population, the tumor started to grow again and would no longer be sensitive to the first treatment.We then changed this drug to DrugB.In this study, we set the drug-switching time at day 307 (t= 307) based on clinical observations of the median Progression-Free Survival (mPFS) of erlotinib treatment (21).Under the second-line medication of DrugB, the growth of Type-W and Type-X was suppressed, but that of Type-Y and Type-Z increased.Finally, Type-Z became the major population.The simulation results demonstrated a trend of tumor evolution under erlotinib-first treatment.At the end of our simulation period, set as day 1000 (t=1000), the cell number of the tumor was 2.12× 10 12 .Next, we examined the DrugB-first treatment.The main population changed from Type-W to Type-Y and Type-Z (Figures 2B, E).During the treatment period, Type-W and Type-X decreased, whereas Type-Y and Type-Z increased.Once Type-Y became the main population, the tumor started to grow again and was no longer sensitive to DrugB.Herein, we switched drugs at day 567 (t=567) because the mPFS was approximately 567 days under the osimertinib treatment (21).When DrugA was used as the second-line treatment, Type-Y was suppressed, and Type-Z continued to grow and dominated the tumor.Compared with the presumed evolution (Figure 2E), our model profitably reflected the tumor progression of osimertinibfirst treatment.In this treatment schedule, the tumor recurred at day 490, which was longer than that of erlotinib-first therapy.Additionally, at the end of our simulation period (t=1000), the total cell number was 1.09× 10 12 , which was less than that of the DrugA-first treatment.Furthermore, we investigated the outcomes of combination therapy (DrugC) by using DrugA and DrugB at the same time as first-line treatment (Figures 2C, F).When DrugC was applied as the Time course simulation results of monotherapy and combination therapy.The results of the simulations are depicted in (A-C).The x-axis is time, and the y-axis is the cell number.The blue, orange, yellow and green curves represent the dynamics of Type-W, -X, -Y, and -Z, respectively.The purple curve represents the total cell number.The expected tumor progression tendencies are depicted in (D, E, F).The blue, orange, yellow, and green cells are Type-W, -X, -Y, and -Z, respectively.In the simulation of erlotinib-first (A), the main population changed from Type-W to Type-X for a while.After changing erlotinib to osimertinib at day 307, Type-X decreased, and Type-Z became the dominant population in the end.This simulation result represents the tumor evolution tendency shown in (D).The simulation product of osimertinib-first is shown in (B).The tumor response to osimertinib increased in the beginning, but as Type-Y became the main population, osimertinib-resistance appeared.After the change to erlotinib at day 567, Type-Y decreased, and the tumor response to treatment increased again.However, Type-Z became the central population causing drug resistance.This simulation result represents the tumor evolution tendency shown in (E).In the combination therapy (C), the main population changed from Type-W to Type-Z.This result represents the tumor evolution tendency shown in (F).

Yu et al.
10.3389/fonc.2023.1137966 Frontiers in Oncology frontiersin.orgfirst-line treatment, Type-W, Type-X, and Type-Y decreased, and only Type-Z continued to increase because it was resistant to both DrugA and DrugB.Type-W is the main population in the early phase of treatment, and eventually, Type-Z replaced Type-W to become the main population in the tumor.During medication, neither X nor Type-Y dominated the population.At day 1000, the total cell number was 4.2 × 10 12 .

Drug-switching time and final cell number
To examine the relationship between drug-switching time and the development of the total cell number of the tumor, we simulated the tumor dynamics and measured the total cell number at day 1000 with various drug-switching times (Figure 3).In the case of DrugA-first treatment (Figure 3A), the lowest final total number of cells was 2.0× 10 12 , while it was about 1.0 × 10 12 in the case of DrugB-first treatment (Figure 3B).This implied that first-line treatment with DrugB displayed better treatment outcomes than DrugA-first treatment.Moreover, the total number of cells at day 1000 remained essentially the same in an appropriate range of drug-switching times under both DrugA-and DrugB-first treatments.This suggested the existence of an optimal drug-switch period, and it was not advisable to switch drugs too early or too late.Furthermore, comparing the suitable drugswitching time period for these two treatments, DrugB-first therapy had a broader range than DrugA-first.In the DrugB-first treatment, switching drugs from days 200 to 900 was acceptable (Figure 3B).However, in DrugA-first therapy, the suitable drug-switching time ranged from day 100 to day 450 (Figure 3A).

Cell initial proportion dependence
To investigate the effect of the initial proportion of different mutant cells on the final cell number, we simulated how the final cell number changes with the increase of mutant cell proportion in different treatment strategies (Figure 4).We explored the effect of only one resistant cell type at one time, keeping other conditions constant as the standard condition.For Type-X and Type-Y, we tested the change in initial proportion from 10 -8 to 10 -1 , and for Type-Z from 10 -9 to 10 -5 .With the increase of Type-X cell (Figures 4A-C), the final cell number did not change under DrugB-first therapy (Figure 4B) and combination therapy (DrugC) (Figure 4C) but increased in DrugA-first therapy (Figure 4A) once the initial proportion of Type-X exceeded 10 -4 .Similarly, in Type-Y dependence simulations (Figures 4D-F), the final cell number increased only in DrugB-first therapy (Figure 4E) when the initial proportion of Type-Y became larger than 10 -3 .In addition, in Type-Z dependence simulations (Figures 4G-I), the final cell number increased once the initial proportion of Type-Z cell became larger than 10 -7 under all treatments.Within the range of initial cell proportion that did not cause an increase in the final cell number, DrugB-first therapy always showed the smallest number of total cells at day 1000.

Therapy selection map
In order to identify which treatment strategy is optimal in a given case, we compared the final total cell number in different treatments with the change of the initial Type-X and Type-Y cell proportion and summarized the results in a therapy selection map (Figure 5).In this simulation, we kept the initial number of Type-Z constant as 10.By comparing the final total cell number under these three treatment strategies in the different initial proportions of Type-X and Type-Y cells, we determined the best strategy by realizing the smallest cell number at day 1000.The simulations were performed in the same method as used in Figure 2. From this map, we noticed that DrugB-first therapy was the optimal choice when tumors harbored a low initial proportion of Type-Y cells.However, DrugA-first therapy could still be advisable if the initial proportion of Type-Y cells was more significant in the tumor cluster.Furthermore, this map indicated that when both Type-X and Y cells had a high initial proportion in the tumor cluster, combination therapy (DrugC) was the optimal choice.

Parameter sensitivity analysis
To investigate the parameter sensitivity, we analyzed how the total cell number at day 1000 changed with parameters under those Drug switch time and final total cell number.The x-axis is drug-switch time, and the y-axis is the total cell number at day 1000.In panel (A), the simulation result using erlotinib as a first-line treatment is shown.The lowest total cell number at day 1000 is approximately 2 × 10 12 .In panel (B), the case of first-line Osimertinib treatment is shown.The lowest total cell number at day 1000 is approximately 1 × 10 12 .three treatment strategies (Figure 6).In the analysis of the growth rate of Type-W cell (a), the final total cell number increased with the increase of a A under DrugA-first therapy (Figure 6A1); the increase of a B under DrugB-first therapy (Figure 6A2); and the increase of a C under the combination therapy (Figure 6A3).As for the growth rate of Type-X cell (b) and Type-Y cell (c), they did not affect the final total cell number significantly in our simulated value range (Figures 6B, C).Moreover, about the growth rate of Type-Z cell (f), the final total cell number increased with f C under the combination therapy (DrugC) (Figure 6D).Concerning the effect of mutation rates (Figure 6E-I), their influence was different based on therapy strategies.In DrugA-first therapy, the increase of g A , h A , k A , p A and q A increased the final total cell number.Meanwhile, the increase of g B , h B , k B , p B and q B increased it in DrugB-first therapy.In the combination therapy (DrugC), only k C increased it (Figure 6).

Parameter dependence on the therapy selection map
Since several parameter values were set by our own assumptions, we investigated how these values affected the optimal choice of treatment in detail (Figure 7).In this analysis, we changed one focused parameter value, made a therapy-selection Initial Proportions of Type-X and -Y cells and treatment strategy selection.The x-axis is the initial ratio of Type-X, and the y-axis is the initial ratio of Type-Y.The yellow region means that osimertinibfirst therapy is the optimal therapy, the blue region means erlotinibfirst therapy is the optimal therapy, and the green region means combination therapy (erlotinib+osimertinib) is the optimal choice.map as described in Figure 5, calculated the area of each strategy on the map and showed the area composition of each strategy at each parameter value.Especially, we investigated the dependence of growth rates of Type-W (a) and Type-Z (f) cell, and mutation rate from Type-W to Type-Z cell (k) under the three treatment strategies.As a result, the area where DrugB-first therapy exhibited superiority was large in the a A , a B , and a C ,-dependent analyses (Figures 7A-C).When the DrugB effect was weak against Type-W cell, the DrugA-first therapy became superior (Figure 7B).Moreover, when we changed growth rates of Type-Z under the three strategies, the DrugB-first therapy was the best option again except the cases where the growth rate of Type-Z under DrugA and DrugB was large (Figures 7D, E), and the growth rate of Type-Z under DrugB and DrugC was small (Figures 7E, F).Finally, changing the mutation rate under the three treatment strategies, DrugB-first therapy was the best option in most cases (Figures 7G-I).When the mutation rate (k) was small under DrugA and DrugC, and large under DrugB, DrugA-first or DrugC therapy became the best option (Figures 7G-I).

Discussion
In this study, we proposed a new mathematical model of EGFRmutated NSCLC.First, our model successfully reproduced the process of tumor evolution under different treatment schedules, including monotherapy and combination therapy (Figure 2).In the erlotinib-first treatment (i.e., DrugA-first treatment), the drug was switched at day 307, while at day 576 it was switched in the osimertinib-first treatment (i.e., DrugB-first treatment).Next, we compared the effects of the two therapies.Our simulation results indicated that first-line osimertinib therapy was better than erlotinib.Within the same time period, for example, 1000 days in our study, osimertinib-first therapy resulted in a lower total cell number.Furthermore, the tumor recurred at nearly 500 days in osimertinib-first therapy compared to approximately 300 days in erlotinib-first therapy.This implied that first-line osimertinib therapy could suppress the growth of tumors more effectively than first-line erlotinib therapy, and could prolong the time of tumor recurrence.In the FLAURA project, clinical statistical data also revealed that EGFR-mutated NSCLC patients treated with osimertinib-first therapy had longer median mPFS (20,21).This statistical study indicated the validity of our proposed model.Additionally, we noticed that in monotherapy, the total cell count was relatively low over a range of drug-switching times (Figure 3).This finding describes the existence of a suitable drugswitching phase, which suggests that it is not advisable to change the drug at a very early or late stage.In the suitable range of drug switching times, our simulation results showed that osimertinibfirst therapy had a relatively lower total cell number than erlotinibfirst therapy at day 1000.This result also indicates the potential of osimertinib as a first-line therapy in clinical applications.In Parameter sensitivity analysis.Parameter dependence on the total total cell number at day 1000 under the three therapy strategies was analyzed.In monotherapy related analysis, the x-and y-axis are the parameters in the effect of DrugA and DrugB, and the color bar presented the final total cell number.In combination therapy, the x-axis is the parameter, while the y-axis is the final total cell number.The analysis of growth rates, a, b, c, f is showed in (A-D), the mutation rates analysis is showed in (E-I).addition, the appropriate drug-switching time range in osimertinibfirst treatment was broader than that in erlotinib-first.Furthermore, we explored the influence of tumor heterogeneity on therapeutic effects (Figure 4).By analyzing the relationship between the initial cell number and total cell number at the end of the tested time, we learned that the therapeutical effects depended on the initial ratio of resistant types in untreated tumors, and sensitive type did not affect it.According to the simulation results, when the initial ratio of Type-X exceeded the threshold, only the total cell number in the erlotinib-first therapy became large (Figure 4A).In the case of Type-Y, only osimertinib-first therapy resulted in large number of cells (Figure 4E).As for Type-Z, when its number became sufficiently large, the final total cell number developed rapidly in all the tested treatment schedules (Figures 4G-I).These results indicated that a high proportion of drug-resistant cells is associated with poor treatment efficacy.This conclusion suggests that if the tumor harbors a high ratio of Type-X, osimertinib-first is better than erlotinib-first.However, with a high initial ratio of Type-Y, erlotinib-first was better.Importantly, by combining this information, we for the first time theoretically revealed the relationship between the choice of treatment strategy and the initial proportion of Type-X and -Y cell (Figure 5).These findings indicated the advantage of first-line osimertinib treatment and revealed the influencing factors when determining treatment plans.Parameter sensitivity analysis about the total cell number and the best treatment choice confirmed the region where osimertinibfirst therapy was superior to other options (Figures 6 and 7).Especially, we noticed that among the parameters, growth rate of Type-W and Type-Z cell and mutation rate from Type-W to Type-Z made a significantly change in the therapy selection map (Figure 7).Thess findings indicated the importance of suppressing all-drug-sensitive (Type-W) and all-drug-resistant (Type-Z) cells.This implied that during the treatment, not only the emergence of secondary resistant cells, but also the response of all-drug-resistant and -sensitive cells to drugs should be considered.
Furthermore, the simulation results showed that the combination of two types of EGFR-TKIs (erlotinib + osimertinib) as first-line therapy has the potential for clinical applications.Recently, several clinical studies have combined two types of EGFR-TKIs (first-and third-generation TKIs) as first-line treatment.In 2017, Wang et al. first reported the combination of erlotinib and osimertinib in patients with EGFR-mutated NSCLC patients (37).They illustrated the expediency of this type of treatment strategy.In addition, Rotow et al. applied gefitinib plus osimertinib as the first-line treatment for untreated patients with EGFR-mutated NSCLC (43).Their results showed the feasibility of conducting EGFR-TKI combination therapy, and survival analysis is in progress.In this study, we explored the combination of erlotinib and osimertinib as first-line therapy and explained the advantages of this method from a theoretical level.Especially, the drug response time with combination therapy was longer than that with monotherapy.Based on our simulation results, the recurrence time under combination therapy was longer than 500 days, whereas it was approximately 450 days in osimertinib-first treatment and 300 days in erlotinib-first therapy.This finding implied that the ability of combination therapy to prevent the emergence of acquired mutations and prolong the drug response time was even better than osimertinibfirst treatment, which suggested its potential in clinical applications.
Based on the above, the versatility exhibited by the simulation results suggests that our model has the potential to be applied to simulate other similar cases in different cancer types.For further study, some clinical information about patients, such as age, sex, and the degree of malignancy of the tumor, may be considered in the parameter estimation.Thus, this model can be used to develop individual treatment schedules in the future.Parameter dependence on the area compositions of the three strategies in the optimal strategy map.We analyzed the composition of optimal therapies among DrugA-first therapy, DrugB-first therapy, and DrugC therapy in therapy selection map with the change of parameters.The x-axis is the parameter to be focused, the y-axis is the percentage of the area of the optimal therapies in the therapy selection map.The dependence of growth rates of Type-W (A-C), and Type-Z cell (D-F), and mutation rate from Type-W to Type-Z cell (G-I) under the three treatment strategies were tested.The yellow, blue, and red lines are DrugA-first therapy, DrugB-first therapy, and DrugC therapy, respectively.

4
FIGURE 4 Relationship between the initial proportion and final cell number.The total cell number at day 1000 of each therapy with the different initial cell numbers of Type-X cells are shown in (A-C).The final total cell number of each therapy with the different initial cell numbers of Type-Y cells are shown in (D-F).The total cell number at day-1000 of each therapy with the different initial cell numbers of Type-Z cells are shown in (G-I).The xaxis is initial proportion of mutation cells, the y-axis is the final total cell number.

TABLE 2
Parameter values with different therapies.

TABLE 1
Parameter notation in the mathematical model.