How combination therapies shape drug resistance in heterogeneous tumoral populations

Treatment of cancer relies increasingly on combination therapies to overcome cancer resistance, but the design of successful combined protocols is still an open problem. In order to provide some in...


Introduction
Cancer appears in organisms not as a single cancer clone but as a variety of clonal populations which undergo selection by processes of competition and predation: tumour cells compete for resources and are attacked by immune cells. The fitness of a tumour type, i.e. the ability to adapt and grow, depends on how effectively it outcompetes the others and how successfully develops mechanisms to prevent detection and elimination by the immune system. In turn, natural selection promotes cell clones that have acquired advantageous heritable characteristics (Gerlinger & Swanton, 2010). Thus, cancer types can CONTACT M. Delitala marcello.delitala@polito.it be considered coexisting populations, embedded in an environment comprising normal and immune cells (Marusyk, Almendro, & Polyak, 2012): the growth, or otherwise, of different tumour types depends, clearly, on a variety of environmental factors, so it makes sense to talk of an ecology of cancer (Hillen & Lewis, 2014;Pacheco, Santos, & Dingli, 2014). Cancer heterogeneity has been invoked to explain one of the major aspects of cancer development, namely acquired drug resistance, by which phases of remission are often followed by a rapid growth of tumour cells. Cancer cells can acquire resistance to chemotherapy by a range of mechanisms but the results, however, are the emergence of a resistant population, (Dean, Fojo, & Bates, 2005;Pisco et al., 2013). Thus, even though therapies exist that can decimate a given cancer type, if one or more variants are present in the tumour population which are resistant, the resurgence of treatment-refractory disease may occur, Gerlinger and Swanton (2010). The coexistence of cancer types makes the design of optimal protocols a challenging problem, and specifically the optimal dosing and timing of combination of chemotherapy and immunotherapy are still an open issue, Slovin (2012).
Obviously, given the relevance of the problem, it is not surprising that a variety of mathematical tools has been used to model cancer dynamics and the effects of therapies. For a review see for instance (Bellomo, Li, & Maini, 2008;Eftimie, Bramson, & Earn, 2011;Eladdadi, Kim, & Mallet, 2014;Wilkie, 2013) and references therein.
The work presented here concerns the analysis of different therapies on two competing cancer populations with particular concern on the role of competition on the outcome of medical treatments.
For this purpose, a model based on population theory (Murray, 2002;Hofbauer &Sigmund, 1998) andpresented in Piretto, Delitala, andFerraro (2018) is used as an in silico laboratory to simulate the outcomes of protocols differing in doses and timing of the treatments and to provide an analysis of the results.
Two kinds of treatment are simulated: a therapy reducing the rate of growth of cancer used alone or in combination with drugs increasing the effectiveness of the immune system. Time courses of medications are, throughout the paper, a metronomic (André, Carré, & Pasquier, 2014) and a pulsed therapy; in last section, the switch to pulse protocol, (Foo, Chmielecki, Pao, & Michor, 2012), has been simulated.
In the sequel, first some general simulations are presented, centred on the study of the therapies effects for different cases of inter-clonal competition, and next a specific type of cancer is considered, namely the non-small cell lung cancer (NSCLC).

The model
Basic elements determining the evolution of the cell populations are proliferation, predation and competition for resources. Two competing types of cancer are considered together with immune system cells, and the model is formulated as a system of three differential equations: (1) The model is proposed and described in greater details in Piretto et al. (2018). Here we briefly summarize the main assumptions: • Tumour clone x 1 . The first two terms in the RHS of the equation represent the growth of x 1 in isolation, i.e. in absence of other cancer types, of the immune system and of medical treatments. In this case, x 1 undergoes a logistic growth and the parameter r 1 is the reproduction rate whereas K 1 corresponds to the carrying capacity of the environment, i.e. the maximum value that x 1 can take. Development of x 1 is constrained by the competition with clone x 2 (measured by the parameter b 12 ) and by the interaction with the immune system (parameter c 1 ). • Tumour clone x 2 . The equation for x 2 is analogous to the one for x 1 and the same considerations apply. • Immune system z. It grows with a net rate β, and in absence of tumour, it is limited by H. In presence of cancer, z undergoes a clonal expansion weighted, respectively, by parameters α 1 , α 2 which measure the ability of immune cells to detect and recognize cancer cells. • Chemotherapy. The effects of chemotherapy on the population x i are represented in the model by the term g i (t)x i where g i takes into account the drugs kinetics in the organism, see (De Pillis & Radunskaya, 2001). This is equivalent to rewrite the growth The action of the immune system can be enhanced by immunotherapy, whose effect is expressed by the term h(t)/K i x i z, h(t) > 0 which is equivalent to define a new parameter κ i (t) = c i + h(t). More specifically, the immunotherapy modelled here is a vaccination (DCs transduced with adenovirus containing fulllength mouse wild-type p53 (Ad-p53)) that results in the generation of immune system cells (CTLs specific for p53-derived peptide) inducing a specific antitumoral immune response (Nikitina et al., 2001). The effect of the immunotherapy is then an increased ability by the immune system to recognize and kill cancer cells as proposed, for instance, in Wilson and Levy (2012). In the following, x = x 1 + x 2 denotes the total cancer load, and for simplicity's sake, The model considered here, albeit, admittedly, a drastic simplification of the processes underlying cancer evolution, includes some of the important factors involved in cancer development, i.e. the activity of immune system, the competition between two cancer clones and the effects of medical treatments.
Obviously, a cancer population may contain several different clonal types; however, restricting to just two is justified by the fact that experimental literature usually distinguishes between two main populations, i.e. sensitive or resistant to therapy (Foo et al., 2012).
Also the representation of the immune system proposed here is simplified, as it does not consider the several types of immune cells, with different effects on the tumour population, the system is composed of, but, rather, one population acting as predator on cancer cells (Kuznetsov, Makalkin, Taylor, & Perelson, 1994). The rationale for this simplification is that here the main focus is the development of cancer populations, under different selective pressures, of which the action of the immune system is a component and, therefore, has been modelled in a simple way.
Finally, we are aware that chemotherapy not just affects cancer cells reproduction rate, but it has at least a twofold effect (Zitvogel, Apetoh, Ghiringhelli, & Kroemer, 2008) on the immune system and on the healthy host tissue: on one hand, it weakens the proliferation of immune cells, and on the other, it elicits or reactivates anticancer immune response, enhancing the immunogenic characteristics of tumour cells. Here the focus is on the primary aspect of chemotherapy, namely the reduction of cancer cells proliferation.

Stability analysis
A detailed stability analysis of system (1) and the stationary points has been presented in Piretto et al. (2018), and it can be carried out with the standard methods of analysis, see (Murray, 2002;Hofbauer & Sigmund, 1998). The main results are reported here in view of the discussion of Section 3.
In the positive orthant R 3 + , the system (1) has 8 stationary points, whose components will be denoted by x * 1 , x * 2 , z * . Four of these points are of the form P = (x * 1 , x * 2 , 0) with x 1 ≥ 0, x 2 ≥ 0, and they are clearly unstable as dz/dt > 0 if z < H: thus, they are of no biological interest.
The other points correspond to the case of tumour eradication, survival of a single cancer clone (competitive exclusion) and coexistence of both cancer clones. In more detail: • Tumour eradication. Tumour eradication corresponding to P te = (0, 0, H), occurs when: in this case P te is the only stationary point, and it is globally asymptotically stable, corresponding to the well known fact that the immune system routinely suppresses tumoral cells in the organism.
• Competitive exclusion: Only one tumour type survives, and then two cases are possible: EC1 corresponding to the survival of just the clone x 1 (stationary point Permanence of a single type can happen in two ways: either one of conditions (2) does not hold, and the corresponding population survives, or both cancer clones escape eradication by the immune system, but one dominates the other and the principle of competitive exclusion applies. In this case, there exists a single globally stable stationary point, namely P ce1 , or P ce2 , depending on which clone prevails.
• Coexistence. In this case, called COE for short, both tumour types are present and the stationary point is P coe = (x * 1 , x * 2 , z * ), with all components greater than zero. Coexistence occurs when both types survive, the action of the immune system and competitive interaction are not so strong to cause competitive exclusion, the corresponding stable point is P coe .
In conclusion, stability analysis allows to make explicit which conditions should be satisfied by parameters for the occurrence of the different trends. Thus, what is of interest are not the stationary points, per se, but rather the correspondence between parameters and composition of the cancer population.
In the following, it will be assumed that, before therapy, both cancer clones escape eradication by the immune system and hence the system will be in a configuration leading, asymptotically, to EC1 (resp. EC2), i.e. survival of just x 1 (resp. x 2 ) or to COE, i.e. coexistence of the two cancer clones.

Simulations
The model outlined in the previous sections is the basis for in silico experiments. A variety of simulations have been carried out, corresponding to different cancer types under the effect of different protocols (with respect to the dose and the time schedule), involving chemotherapy, immunotherapy and molecular target therapy.

Therapies
As mentioned before the effect of chemotherapy is simulated by decreasing the rate of growth r i to f i , whereas simulation of immunotherapy consists in increasing the efficiency of the immune system, from c i to κ i .
Thus, administered doses are expressed in terms of differences between new and old values of these parameters, where [T i , T f ] is the time interval of administration of the treatments. Any measure of the efficiency of a treatment needs to take into account not just the decrease of the total cancer load x = x 1 +x 2 but also the trend of resistant cells numerosity, a critical factor in assessing how well a therapy works. A parameter depending on variation (or otherwise) of the total load and of the number of resistant cell can be defined as follows.
Let parameters , and δ be given by: here x pre , x post are respectively the total cancer loads one day before and one day after the administration of the therapy, whereas x 2,pre and x 2,post are, respectively, the values of x 2 (the resistant clones) before and after therapy. Parameters and δ are not independent, as x 2,post ≤ x post ; indeed, it is straightforward to show that δ is a linear function of and, in particular, δ = 1 if = 1. Note that and δ can take negative values, whereas, to combine them in a global parameter, they need to be non-negative. To this aim, a transformation is applied to , δ by means of the sigmoid function .
where the division by g (1) is just to ensure that sup y = sup y δ = 1. Variables y , y δ define a feature space where outcomes of different therapies can then be mapped. An efficiency index that depends on the decrease of both x and x 2 can be defined in a natural way as γ = y y δ .
When = δ = 0, γ = 0.467 and thus one can choose γ = 0.5 as a threshold for the effectiveness of therapies (i.e. γ > γ corresponds to a positive treatment). Note that δ is usually very small in our simulations; however, what matters here it is not the absolute value of the parameter but, rather, the sign and the relative values in the different cases.
Obviously no index can fully capture the complexity involved in the results of medical treatments and parameters introduced here are meant to give just some insight on therapeutic outcomes.

Parameters and simulation settings
Parameters used in the simulations have been evaluated from data reported in the literature (De Pillis, Gu, and Radunskaya, 2006;Kuznetsov, Makalkin, Taylor, & Perelson, 1994, or, as in case of b 12 , b 21 , varied in an exploratory way. They are reported in Table 1. A sensitivity analysis has been carried out in Piretto et al. (2018) to evaluate the effects of an error in the parameter estimation. The parameters with the greater impact on the output of the model (i.e. r i , K i and H) are set to fit the data of (Ramakrishnan et al., 2010).
The different parameters and simulation settings used in the sequel can be summarized as follows: • Cancer type. It is assumed that x 1 has a larger proliferation rate, r 1 > r 2 , and it is susceptible to chemotherapy, whereas x 2 is resistant (f 2 = r 2 ). According to  (Foo et al., 2012), the slower grow rate of drug-resistant cells is due to the evolutionary cost to be resistant. Immunotherapy affects both cancer clones in the same way.
In addition, different competitive efficiencies, expressed by b 12 , b 21 , are considered, corresponding to clones competition for resources, as, for instance, space, glucose or oxygen.
We focus on three types of cancer populations, characterized by different asymptotic behaviours, as highlighted in Section 2.1: -Coexistence (COE). Before therapy, x 1 , the susceptible clone, is also the fittest, as it has the larger rate of growth and also a competitive advantage (b 12 = 0.02 and b 21 = 0.07) even though x 2 is not totally eliminated. -Competitive exclusion (EC1). x 1 is the dominant clone. In this configuration, x 1 increases its competitive advantage (b 12 = 0.1 and b 21 = 0.5) so that, asymptotically, x 2 is wiped out. -Competitive exclusion (EC2). The resistant clone x 2 dominates x 1 (b 12 = 0.5 and b 21 = 0.09), so that, asymptotically, x 1 disappears.
• Immune system. The immune threshold H is also derived from the literature (Kuznetsov, Makalkin, Taylor, & Perelson, 1994), corresponding to a patient with standard immune system. Obviously, H can vary with age and in particular is lower for older patients. • Initial conditions. Cancer evolution starts with resistant population, x 2 , small with respect to the susceptible, x 1 , i.e. x 1 (0) = 9 · 10 4 , x 2 (0) = 4 · 10 4 , and the immune cells are z(0) = 5 · 10 2 . • Therapy dose. In the simulated medical protocols, two therapies (chemo-and immuno-) are administered simultaneously. The effects of two doses of chemotherapy, standard and low, denoted, respectively, by Ch + and Ch − , corresponding to different values of integral (3a), are evaluated in combination with the immunotherapy, denoted by I. More precisely the combinations considered here are Ch + (standard dose of chemotherapic drug alone), Ch − I (low dose of chemotherapy in combination with immunotherapy) and Ch + I (standard chemotherapic dose combined with immunotherapy). Total doses are Ch − = 2.41, Ch + = 9.66, I = 971. Standard and low doses are set according to the data of (Ramakrishnan et al., 2010). • Time protocols. Two time courses are examined here: a metronomic therapy (André, Carré, & Pasquier, 2014), i.e. an administration of the two drugs every day guaran-teeing a quasi stationary level of drug in the bloodstream and a pulsed protocol (Foo et al., 2012) in which drugs are administered in a series of "pulses" each followed by an interval of no therapy (one pulse every two weeks).

Some in silico experiments
In general, we have found that outcomes from the pulsed and metronomic (continuous) therapies are similar, but in some cases the pulsed therapy yields better results and therefore the graphs refer to it. Moreover, we are interested here on the evolution of cancer and thus, for simplicity's sake, in figures, the time course of the immune system z is not included. Results, in the different cases of competition, are summarized here and they refer, unless otherwise indicated, to the pulsed therapy: • Coexistence. Standard doses of chemotherapy alone Ch + are inefficient, γ = 0.15; the negative value of ( = −0.54) indicates that the total cancer load x keeps growing during the treatment (see the purple curve in Panel (a) of Figure 1 and the red circle in Figure 4). The population of resistant clone x 2 survives and increases unchecked by x 1 (δ = −1.26), and thus, after an initial drop, the cancer rebounds, as a consequence of a clonal inversion where x 2 become the larger clonal population, as shown in Panel (a) of Figure 1. A combination of weak chemotherapy and immunotherapy Ch − I is more efficient, γ = 0.61, in reducing the total load ( = 0.58) and preventing the growth of the resistant clone (δ = 0.04): x 1 is kept under control since it is sensitive to both therapies, while the growth of x 2 is curtailed by immunotherapy and by x 1 . The results are presented with a red square in Figure 4 and in Panel (b) of Figure 1.
Results with Ch + I, a combination of standard chemotherapy combined with immunotherapy, are similar, γ = 0.62: obviously, there is a large drop of x ( = 0.82), but now x 2 can grow without interference by x 1 , as indicated by the negative value of δ = −0.09, and indeed the total population is composed mostly by resistant cells (see Panel (c) of Figure 1 and Panel (a) of Figure 4). When the combination Ch − I is administered the metronomic protocol is slightly less efficient (γ = 0.57), whereas in the other two cases there is no difference. Results do not vary if the values of parameters b ij are exchanged within the COE configuration.
• Competitive exclusion: clone x 1 dominates. In this case, Ch + is effective, γ = 0.62 ( = 0.66); given the superiority in fitness of x 1 even a reduced population of this clone it is sufficient to keep x 2 small (see Figure 2 Panel (a)). The combination Ch − I is less efficient than in COE configuration, γ = 0.60 ( Figure 2 Panel (b)); here what matters is to control x 1 and thus the introduction of immunotherapy cannot, in this case, compensate the reduction of the chemotherapy dose. Finally Ch + I is an effective treatment (γ = 0.68) and the number of cancer cells becomes very low, almost to the level of eradication ( = 0.98), see Figure 2 Panel (c). The metronomic protocol is slightly less effective for Ch + (γ = 0.60) and Ch − I (γ = 0.57).
• Competitive exclusion: clone x 2 dominates. In this case, one should expect to find therapy outcomes not much different from COE, but with the effect of clonal inversion somehow intensified by the fact that x 2 is now dominant. For instance Ch + (γ = 0.08) is useless, as it merely accelerates the disappearance of x 1 , already doomed by the  competition. The result is an increased total load ( = −0.84), composed almost entirely by x 2 (δ = −1.75). This effect is clearly shown in Figure 3: the Panel (a) shows the time courses of x 1 , x 2 and x during the administration of Ch + with a pulsed therapy and Panel (b) displays the fraction x 2 /x of the resistant clone during the treatment's cycle. It is interesting to note that in this case Ch − I has the same effectiveness as Ch + I (γ = 0.61 and γ = 0.60, respectively), since, even though the latter yields a large reduction of x ( = 0.61, = 0.79, resp.), the combination Ch − I gives a better performance in controlling x 2 (δ = 0.03, compared with δ = −0.09).
In conclusion, competitive interactions contribute to shape the outcome of therapies. In cases COE and EC2, combination therapies have comparable efficiency, Ch + I being more effective against the growth of x and Ch − I giving a better performance in controlling x 2 . In other words, there exists a trade off between the decrease of the total load and the control of resistant clones; of course, Ch − I has the advantage of less dangerous side effect on the patient and this could make it preferable.
On the contrary, in EC1, when the sensitive clone is strongly dominant, Ch + yields the best results, since, as mentioned before, the growth of x 2 is curtailed by x 1 competitive advantage.

Drug holiday
Evolution of drug resistance is a dynamic process (Dhawan et al., 2017;Sharma et al., 2010); thus in clinical protocols, repeated courses of treatment are often separated by a periods of "drug holiday" during which no therapy is administered (Carrere, 2017;Pouchol, Clairambault, Lorz, & Trélat, 2017). This method enables a regrowth of sensitive cells and a corresponding reduction of resistant ones thus reestablishing the effectiveness of a therapy; moreover, it allows patients to recover from the therapy side effects. Its results depend on a variety of factors, e.g. the effect of the sensitive clone on the growth of resistant cells and the timing and duration of the holiday.
In the framework of the present model, one should expect to find the greater effect when populations x 1 (sensitive) and x 2 (resistant) are in a EC1 case. Simulations reported below compare the outcomes of the same treatment with and without drug holiday; in both cases, values of x pre , x 2,pre are taken at the beginning of the therapy (t = 50 days) and x post , x 2,post just after the end of the second part of the treatment (t = 200 days).  The administration of Ch + , in the EC1 case, clearly shows the effect of a drug holiday and, in particular, the decrease of the resistant clone (see Figure 5), besides, obviously, confirming that chemotherapy alone is inefficient (γ ≤ γ = 0.5 in both cases).
When a combination therapy is adopted the action of x 1 over x 2 is less relevant since immunotherapy affects both cancer types and then the growth of x 2 is limited in any case, see Figure 6 where the results for different time spans of the drug-free interval are compared. Values of γ are, respectively, γ a = 0.48, γ b = 0.57, γ c = 0.67, with indices referring to the corresponding panel. Note that, even though the best result is obtained in Panel (c), during the period of no treatment there is a large increase of the total load, an effect that must be obviously taken into account, when planning a drug holiday.
Finally if a Ch + I therapy is used, the population of resistant cells x 2 can become dominant during the holiday, since x 1 becomes so small that its action on x 2 is negligible.
Outcomes of the therapies, for the EC1 case, are shown in Figure 7, and it is apparent that results for Ch + I (γ = 0.68) and Ch − I (γ = 0.67) are very similar.
In the COE case, γ = 0.69 for both combination therapies, but Ch − I is slightly better in checking the growth of the resistant cells, with δ = 0.07 vs. δ = 0.04 for Ch + I.  (6)) for the EC1 configuration, comparing treatments with and without the drug holiday.

Application to lung cancer
So far our treatment has been quite general; here, we aim to provide some more specific conjecture on the action of therapy. For this reason, we focus our study on non-small cell lung cancer (NSCLC).
This type of cancer shows resistance to several therapies. One of the most clinically effective treatment involves the use of a small-molecule tyrosine kinase inhibitor (TKIs) Erlotinib (Pao et al., 2004). Erlotinib is an orally active selective inhibitor of the epidermal growth factor receptor (EGFR) tyrosine kinase (Shepherd et al., 2005), whose mutation appears in 70% of patients with NSCLC. This drug is part of the more general class of molecular target therapies.
Recent experiments show that Erlotinib improves both survival and quality of life, compared with placebo in patients after first or second-line chemotherapy (Shepherd et al., 2005). These results are nevertheless controversial (Garassino et al., 2013) because one of the most serious problems is the existence of a second mutation resulting in a new cancer type resistant to Erlotinib (Foo et al., 2012).
Erlotinib has the effect of reduce the rate of growth of a population of cancer cell TKI-sensitive. For this reason, we modelled the effect of this drug as done before with the chemotherapy namely rewriting the tumour growth term as f i (t)x i = (r i − g i (t))x i where g i takes into account the drugs kinetics in the organism. However, in order to perform the simulations, we must make sure that the parameters fixed before apply also at NSCLC. To this aim, we compare a simplified version of our model (consisting only of two populations of cancer cells, without immune system population) with data from in vitro experiments, presented in [ Figure 3B] (Chmielecki et al., 2011), during which a mixture of sensitive and resistant cells was treated for 3 days with increasing concentrations of Erlotinib and afterwards the ratio x post /x c was measured, where x post is the number of cells after the treatment and x c refers to control. In Figure 8, experimental data are compared with simulation results. Data are represented by squares on coloured lines, each colour corresponding to a different initial mixture of resistant and susceptible cells, whereas coloured circles show the simulation results for the same mixtures. We have considered different growth kinetics for the two populations, in particular the drug-resistant cells have a slower growth as proposed in Foo et al. (2012). The comparison shows good agreement with data for most mixtures of the two types.

Clinical protocols
In order to test some of the clinical protocols used in the treatment of NSCLC, the schedule proposed in Foo et al. (2012) has been simulated.
• Figure 9, Panel (a). A pulsed protocol has been considered, where a dose of 1600 mg of Erlotinib is administrated one time a week for a cycle of 30 days. Each pulse of the cycle reduces the rate of growth of x 1 to f 1 = −1.3 for 36 hours, corresponding to Erlotonib half-life in patients. The reduction has been calculated in a way that, with 5 pulses of 1.5 days the total dose is Ch = 9.66 (compare Equation (3a)). • Figure 9, Panel (b). A low-dose concentration in a metronomic schedule has been taken into account (150mg of Erlotinib every day). The parameter f become f 1 (150mg)= 0.014 in the exact proportions with respect to the standard dose. • Figure 9, Panel (c). A mixed schedule of pulsed and continuous administration (a switch to pulse protocol) has been simulated by adding to the pulse therapy a concentration of Erlotinib administered in a metronomic way (1600mg/week +150mg/day).
We have found that, irrespective of the type of cancer, starting with 5% of resistant cells both switch to pulse and pulsed protocol are very efficient in term of reduction of total cancer load, compare Panels (a) and (c) of Figure 9, where Panel (b) presents, for comparison, the effects of a low-dose metronomic schedule. Obviously in all cases, the decrease is temporary because of the growth of the resistant population. The last panel in Figure 9, representing the growth of the fraction x 2 /x of resistant cells for all analysed schedules and levels of dose, clearly demonstrates that a protocol involving an high selective pressure on the system accelerates the growth of the resistant population.
In order to increase the efficiency of the treatments, and motivated by the results obtained in the first part of this work, another treatment, affecting both cancer types, is modelled, namely the p53 cancer vaccine, a standard clinical in the treatment of NSCLC (Antonia et al., 2006).
In the Figure 10, a low-dose (i.e. approximately 1/3 of the one used in simulations of Section 3.2) continuous immunotherapy schedule is considered. The outcomes show that the response at the end of treatments is better than with the Erlotinib alone, in every schedule considered (compare Figure 9).
The improvement is particularly marked in case of metronomic schedule where a low-dose combination therapy yields results comparable with high-intensity single TKI schedules (pulsed and switch to pulse). Panel (d) of Figure 10 shows how the fraction of resistant cells is affected by the combination therapy; comparison with the same panel of Figure 9 show clearly a slower increase of the fraction of resistant cells for the metronomic schedule.

Conclusions
Results presented here clearly suggest that inter-specific competition among cancer cells has a powerful influence on the outcomes of therapies.
In our study, a standard chemotherapy works effectively when the sensitive clone has a strong competitive advantage (case EC1) whereas is useless in the other cases (COE and EC2). Combination therapies are more effective: the use of standard doses of chemotherapy combined with immunotherapy (Ch + I) is the most efficient treatment in the EC1 case, but it gives, for COE and EC2, results that are similar to a combination with a lower chemotherapic dose (Ch − I) and it is, in all cases, less effective than Ch − I in reducing the number of resistant cells.
Simulations of drug holidays confirm the effectiveness of this method, in that both combination therapies give better results than without interruption of the treatment; moreover, in this case values of the efficiency index γ are very close for Ch + I and Ch − I, also in the EC1 case. Thus, a relatively low dose of chemotherapic drug combined with immunotherapy seems to be preferable in most instances, considering also the obvious advantages in terms of tolerance by the organism. Comparison of the temporal protocols shows that metronomic application of the treatment yields, in general, comparable, albeit in some cases slightly worse, results.
In the case of non-small cell lung cancer (NSCLC), the simulations clearly show that a combination therapy with low doses of both immunotherapy and molecular target therapy (Erlotinib) administered every day can give results similar to the monotherapy of Erlotinib with higher doses. An additional benefit of a low-dose combination therapy is the slower increase of resistance (measured as the relative frequency of resistant cells), besides the obvious advantage in term of tolerance by the organism.
Findings presented here agree with the growing experimental literature showing the superiority in most cases of the combination therapies (Ledford, 2016). However, this cannot be an absolute rule, as it depends on the interactions among cancer populations and hence there may exist cases in which a single therapy performs better, as in case of a strong competitive advantage of the susceptible clone. Thus, a knowledge of the evolutionary interactions among cancer types appears to be necessary for a successful strategy against cancer.
In order to contrast the insurgence of resistant populations, exploiting the competition between the cancer populations seems to be the best recipe. A combination of low-doses therapies could yield a substantial reduction on the total tumour population without impose a massive selective pressure that would suppress the more competitive, but therapysensitive population, which controls clonal types with the evolutionary advantage of resistance to therapies. In this context, drug holiday offers a way to prevent the increase of the resistant clone and our simulations highlight the role of timing and duration in agreement with results reported in Dhawan et al. (2017).
In conclusion the best procedure seems to be, in general, not a treatment aiming to cancer eradication, that may lead the emergence of a resistant clone, but rather protocols that can control cancer, also taking advantage of the inter-clonal competition, as also suggested in Gatenby (2009).

Disclosure statement
No potential conflict of interest was reported by the authors.