Modelling and investigation of the CD4+ T cells-macrophages paradox in melanoma immunotherapies

It is accepted that tumour cells can be eliminated by M1 anti-tumour macrophages and CD8 + T cells. However, experimental results the 10 – have shown that B16 mouse melanoma cells can be eliminated by the CD4 + T cells alone (either Th1 or Th2 sub-types), in the absence of CD8 + T cells. In some studies, elimination of B16 melanoma was associated with a Th1 immune response (i.e., elimination occurred in the presence of cytokines produced by Th1 cells), while in other studies melanoma elimination was associated with a Th2 immune response (i.e., elimination occurred in the presence of cytokines produced by Th2 cells). Moreover, macrophages have been shown to be present inside the tumours, during both Th1 and Th2 immune responses. To investigate the possible biological mechanisms behind these apparently contradictory results, we develop a class of mathematical models for the dynamics of Th1 and Th2 cells, and M1 and M2 macrophages in the presence/absence of tumour cells. Using this mathematical model, we show that depending on the re-polarisation rates between M1 and M2 macrophages, we obtain tumour elimination in the presence of a type-I immune response (i.e., more Th1 and M1 cells, compared to the Th2 and M2 cells), or in the presence of a type-II immune response (i.e., more Th2 and M2 cells). Moreover, tumour elimination is also possible in the presence of a mixed type-I/type-II immune response. Tumour growth always occurs in the presence of a type-II immune response, as observed experimentally. Finally, tumour dormancy is the result of a delicate balance between the pro-tumour e ﬀ ects of M2 cells and the anti-tumour e ﬀ ects of M1 and Th1 cells.


Introduction
The anti-tumour role of the immune system has been documented for at least a century, with one of the earliest studies on the role of immune surveillance against transformed cells being published by Ehrich (1909). The last 20-30 years have seen a very rapid increase in the number of experimental studies that investigate the molecular and cellular mechanisms behind the tumour-immune interactions. However, in many cases, the experimental results are contradictory. For example, Mattes et al. (2003) investigated the anti-tumour effects of two types of CD4 + T cells (Th1 and Th2 cells) on B16 melanoma, and concluded that contrary to the generally accepted idea that the CD4 + T cells have only a helper role, they can actually eliminate tumours on their own via the cytokines they produce. Moreover, the authors showed that while the Th1-tumour interactions led to temporary tumour control followed by tumour escape and growth (see Fig. 1(a)), the Th2-tumour interactions led in the long term to tumour elimination (see Fig. 1(a)). In fact, Mattes et al. (2003) suggested that tumour elimination in the presence of Th2 cells is helped by the influx of eosinophils to the tumour site. In addition to eosinophils, the authors also showed the presence of tumour-infiltrating macrophages (see Fig. 1(b)), which seemed to be associated with tumour growth (but the authors did not investigate the possible anti-tumour/pro-tumour action of these macrophages). In a later study, Xie et al. (2010) showed that the Th1 cells can actually eliminate B16 melanoma cells (see Fig. 2(a)). Kobayashi et al. (1998) showed that the growth of B16F10 cells is associated with a large number of Th2 cells and a high concentration of IL-4 cytokines (see Fig. 2(b)). Moreover, Chen et al. (2011) showed that the growth of B16 melanoma cells is associated with a shift from anti-tumour M1 macrophages to pro-tumour M2 macrophages (see Fig. 2(c)). (Note that the classification of macrophages into M1 and M2 phenotypes mirrors the Th1 and Th2 nomenclature , and despite this strict classification there is actually a continuum of phenotypes between the M1 and M2 extremes).
The anti-tumour effects of Th1 and Th2 cells are exerted by the cytokines they produce: (i) the Th1 cells produce type-I cytokines, such as IFN γ − , IL−2, TNF α − and TNF β − (Magombedze et al., 2014;Lucey Fig. 1. Data approximated and re-drawn from Mattes et al. (2003), where the authors transfer Th1 cells or Th2 cells into C57BL/6 mice that were previously injected with B16-OVA melanoma cells. a) Number of tumour metastases after the adoptive transfer of Th1 cells, Th2 cells and for the control case (i.e. no treatment with immune cells). b) Number of tumourinfiltrating macrophages following the adoptive transfer of Th1 cells and Th2 cells, and comparison with the number of macrophages in control tumours (with no adoptive transfer of Th1/Th2 cells). Cytokine level (pg/ml) in INF−gamma levels (pg/ml) on day 7 Nbr of tumour metastases on day 14 Tumour size (mm2)  Xie et al. (2010), where the authors inject RAG −/− mice (which do not have any CD8 + T cells, B cells or NKT cells) with B16F10 melanoma cells. Panel (i) shows tumour size on day 20 for mice injected with CD4 + T cells and for control mice (with no injection of CD4 + T cells); Panel (ii) shows the level of IFN-γ in mice injected with CD4 + T cells and in control mice, suggesting that the CD4 + T cells that reduce the size of the tumour are actually Th1 cells (which produce high levels of IFN-γ). (b) Data approximated and re-drawn from Kobayashi et al. (1998), where the authors inject C57BL/6 mice with B16F10 melanoma cells. Panel (i) shows the number of metastatic colonies on day 14 after injection; Panel (ii) shows the level of IL −2 IFN-γ and IL −4 cytokines produced by naive CD4 + T cells in normal mice and in mice injected with B16F10 cells. (c) Data approximated and re-drawn from Chen et al. (2011), where the authors inject C57BL/6 mice with B16F10 melanoma cells. Panel (i) shows tumour volume on days 7 and 12 after transfer of tumour cells; Panel (ii) shows the percentage of M1 and M2 macrophages inside the tumour, on days 7 and 14. et al., 1996); (ii) the Th2 cells produce type-II cytokines, such as IL−4, IL−5, IL−6, IL−10 and IL−13 Romagnani (1999); Lucey et al. (1996). It is usually thought that the type-I cytokines (e.g., IFN-γ, IL−2) have an anti-tumour role (Lucey et al., 1996), while the type-II cytokines (e.g., IL−10) are generally associated with tumour growth (Lucey et al., 1996). These cytokines are not only produced by the Th1/Th2 cells, but also by other cells in the environment: e.g., macrophages, neutrophils, eosinophils, etc. (Lucey et al., 1996). In particular, the macrophages can produce, and respond to, both type-I and type-II cytokines. Classically activated M1 macrophages are induced by cytokines such as IFN γ − or TNF α − . Alternatively activated M2 macrophages are induced by cytokines such as IL−4 and IL−13 . Moreover, the M1 cells are associated with Th1 responses, being involved in resistance against tumours . On the other hand, the M2 cells are associated with Th2 responses, being involved in tumour progression, tissue repair and remodelling . We emphasise here the crosstalk between the Th cells and macrophages via the type-I and type-II cytokines, which might influence the tumour microenvironment (see also Fig. 3).
The goal of this study is to derive a class of mathematical models that can propose hypotheses regarding the apparent paradoxical results in the anti-tumour effects of Th1 and Th2 cells, and M1 and M2 macrophages. We note that in the mathematical literature there are various models investigating different aspects of the interactions between Th1 and Th2 cells, and between M1 and M2 macrophages. For example, the Th1-Th2 dynamics was investigated in the context of cell differentiation and cross-regulation (Yates et al., 2000;Bergmann et al., 2001;Fishman and Perelson, 1999), during the immune response to allergens (Gross et al., 2011) and asthma development (Kim et al., 2013), during autoimmune diseases (Louzoun et al., 2001), following T cell vaccination (Severins et al., 2008), during bacterial infection in ruminants (Magombedze et al., 2014), or in the rejection of cancers such as melanoma (Eftimie et al., 2010;Kogan et al., 2013). The M1-M2 dynamics was investigated during macrophage activation post-myocardian infarction (Wang et al., 2012), during wound healing (Yu, 2014), or in the rejection of pancreatic cancer (Louzoun et al., 2014). However, very few mathematical models investigate the interplay between M1/M2 macrophages and Th1/Th2 cells during cancer evolution (den Breems and . For example, the study in den Breems and  investigated (numerically and with the help of sensitivity analysis) the influence of the ratio of M1 and M2 macrophages on early and advanced tumour growth, for normal and mutated tumour cells. The authors showed that their model can only exhibit tumour growth (i.e., no tumour elimination). Moreover, they showed that while a ratio of M2:M1>1 can always predict growth towards tumour carrying capacity, a ratio of M2:M1<1 can lead to either growth towards carrying capacity or growth towards a lower tumour size.
In this study, we will investigate the possible mechanisms that could explain the elimination of B16 melanoma by Th2 cells in Mattes et al. (2003) and by Th1 cells in Xie et al. (2010), and the role played by M1 and M2 macrophages in tumour growth and elimination (given the crosstalk between Th1/Th2 cells and M1/M2 cells via the cytokines they produce; see Fig. 3). To this end we develop two mathematical models: (i) a model for the interactions between the Th cells and macrophages alone, which is used to investigate the type-I and type-II immune responses they generate (where we define a type-I immune response to be the response dominated by Th1 and M1 cells, and a type-II immune response to be the response dominated by Th2 and M2 cells); (ii) a model for the interactions between tumour cells, Th cells and macrophages. We show that tumour can be eliminated both in the presence of a type-I immune response and a type-II immune response. Tumour growth is always associated with the presence of a type-II immune response.
The structure of this article is as follows. In Section 2 we introduce a mathematical model for the Th cells-macrophages interactions and discuss the long-term behaviour of the model by investigating the number and stability of the steady states. We also investigate numerically the dynamics of this model, and discuss the conditions under which the model displays a type-I or a type-II immune response. In Section 3 we generalise the previous model to incorporate also tumour dynamics. Again, we calculate the steady states and their stability to emphasise the complexity of the new model. We also investigate numerically the short-term and long-term dynamics of the model for tumour-immune interactions, and discuss the parameter values for which we see tumour elimination in the presence of a type-I immune response and in the presence of a type-II immune response. We conclude in Section 3.3 with a summary and discussion of the results.

Modelling the Th1 & Th2 and M1 & M2 interactions
We first ignore the presence of the tumour, and investigate the dynamics of the interactions between the Th cells and macrophages, following their cross-talk (via cytokines, which we consider implicitly). Thus we define four variables: the density of Th1 cells (H 1 ), the density of Th2 cells (H 2 ), the density of M1 macrophages (M 1 ) and the density of M2 macrophages (M 2 ). The time-evolution of these variables is given by  2012). These cells grow at a rate p H 1 in the presence of type-I cytokines such as IL −2 (Taylor-Robinson, 1997) or IL −12 (His et al., 1993) (which can be also produced by M1 macrophages), up to maximum carrying capacity m 1 . The growth term also incorporates the competition between the Th1 and Th2 cells for antigens (Magombedze et al., 2014). Note that high Th2 responses lead to a suppression of Th1 responses and vice-versa, as observed experimentally (Magombedze et al., 2014). The natural death rate of Th1 cells is e H 1 (Magombedze et al., 2014).
• The Th2 cells are activated at a rate a H 2 in the presence of IL −4 and IL −13 cytokines that can be produced by M2 macrophages (Romagnani, 1999). Moreover, the Th2 cells grow at a rate p H 2 in the presence of IL −4 (Zhu et al., 2002), up to maximum carrying capacity m 1 . The natural death rate of Th2 cells is e H 2 (Magombedze et al., 2014).
• The M1 macrophages are activated at a rate a M 1 in the presence of IFN γ − cytokine, produced also by Th1 cells (Preuße et al., 2012;Weisser et al., 2013). Also, the M1 cells grow at a rate p M 1 via a self renewal process (Helming, 2011), up to a maximum carrying capacity m 2 . The apoptosis rate of M 1 cells is e M 1 (Gauthier et al., 2013). Note that M1 macrophages can become M2 macrophages, in the presence of type-II cytokines (Allavena and Mantovani, 2012). We denote by r M 1 the re-polarisation rate from M1 to M2 macrophages (Wang et al., 2012).
• The M2 macrophages are activated at a rate a M 2 in the presence of IL −4, IL −13 (which can be produced by Th2 cells) (Weisser et al., 2013). Moreover, the M2 cells proliferate in the presence of IL −4 cytokines characteristic to a Th2-environment (Jenkins et al., 2011) (hence the proliferation rate p H M 2 2 ), up to a maximum carrying capacity of m 2 cells. (Note that, in contrast to the M2 cells, the M1 cells proliferate via self-renewal (Helming, 2011), and thus we do not multiply the p M 1 rate with the H 1 variable.) The apoptosis rate of M2 cells is e M 2 (Gauthier et al., 2013). Finally, since the M2 macrophages can change their phenotype and become M1 macrophages in the presence of type-I cytokines (Allavena and Mantovani, 2012), we denote by r M 2 the re-polarisation rate from M2 to M1 cells (Wang et al., 2012).
We note here that there are a few studies that suggest the possibility of Th1 ↔ Th2 re-polarisation based on the environment (Panzer et al., 2011). However, since this concept of Th re-polarisation is still new, we will not investigate it in this study.
A non-dimensionalised version of the model (1) is shown in Appendix C. However, throughout this study we prefer to work with this dimensional model since in the next two sections we will discuss some of the results in the context of dimensional experimental studies. Moreover, the non-dimensionalisation approach does not reduce significantly the number of model parameters. and thus M M * > * 1 2 (a mixed type-I/type-II immune response, as predicted by the bifurcation diagram in Fig. 5(a)). In regard to the transient immune dynamics: during the first 19 days the Th2 response is lower than the Th1 response, but after day 19 the Th1 response becomes lower than the  (1) and (3), and their values used throughout the numerical simulations. References marked by "*" correspond to parameter values that were approximated based on experimental studies. Some of the elements in column "Value" show not only the specific values used for the simulations, but also the parameter ranges (in parentheses) over which we varied those parameters.  (which leads to H H * < * 1 2 ). (a) Short-term dynamics (panel (i)) and long-term dynamics (panel (ii)) obtained when . For the rest of parameter values see Table A.1.

R. Eftimie, H. Hamam
Journal of Theoretical Biology 420 (2017)   as predicted by the bifurcation diagram in Fig. 5(a)). Note in Figs. 6 and 7 that there are points where the curves have non-continuous derivatives. This is likely a numerical artefact, the result of the number of points used to plot the curves and the scale of the plot.
We conclude that the dynamics of model (1) can be dominated by a type-I, a type-II or a mixed type-I/type-II immune responses, depending on the ratio r r / M M 1 2 and the activation rate of immune cells. Note that for these simulations, we also varied the macrophages activation rates (a a ,

M M
1 2 ) within the interval (10 , 10 ) −4 −2 , but the overall dynamics did not change. We acknowledge that model dynamics might change if we would vary some of the fixed parameters (i.e., those parameters for which we found values in the literature; see Table A.1).

Modelling the Th1 & Th2 and M1 & M2 interactions with tumour cells
Next, we investigate the anti-tumour and pro-tumour effects of M1/ M2 macrophages and Th1/Th2 cells. Thus, we consider five variables: the density of tumour cells (T), the density of Th1 cells (H 1 ), the density of Th2 cells (H 2 ), the density of M1 macrophages (M 1 ) and the density of M2 macrophages (M 2 ). The time-evolution of these variables is given by In addition to the assumptions incorporated in model (1), for model (3) we make also the following assumptions: • Tumour cells grow at a rate α, up to a carrying capacity β (which is chosen to correspond to the maximum tumour size allowed for experimental protocols in mice (NIH, 1996)). To model the phenomenological observation that tumour growth slows down as tumour becomes very large and depletes the available nutrients (Laird, 1964), we choose logistic growth. Tumour cells have a very low natural death (i.e., apoptosis) rate f (Wong, 2011). The Th1 cells kill the cancer cells at a rate g H 1 (via IL −2 and IFN γ − ); see Knutson and Disis (2005). Moreover, the tumour cells can be killed by the Th2 cells at a rate g H 2 (via IL −4 & IL −13 cytokines that attract eosinophils (Mattes et al., 2003)). Also, M1 macrophages kill tumour cells at a rate g M 1 (through the release of tumouricidal products such as NO Lamagna et al., 2006). Finally, the presence of M2 macrophages increases the proliferation of cancer cells (Sica et al., 2008). We denote by g M 2 the proliferation rate of cancer cells in the presence of M2 cells. For simplicity, we . For the rest of parameter values see Table A.1.

R. Eftimie, H. Hamam
Journal of Theoretical Biology 420 (2017)   assumed that all immune cells interact with tumour cells in a linear manner. Under this assumption, the term modelling tumour proliferation can be written as T α g M αT β ( + − / ) M 2 2 , suggesting that the presence of M 2 cells can increase the maximum tumour size. This seems to be confirmed by experimental studies showing that tumours co-inoculated with M2 macrophages grow much larger than control tumours (see, for example, Fig. 5 in Yamaguchi et al. (2016)).
• The Th1 cells can be inactivated by the tumour cells at a rate n H 1 (Magombedze et al., 2014;Eftimie et al., 2010). All other rates that control the dynamics of Th1 cells are as described in Section 2.
• The Th2 cells can be inactivated by the tumour cells at a rate n H 2 (Magombedze et al., 2014). All other rates that control the dynamics of Th2 cells are as described in Section 2.
• The anti-tumour M1 cell population can be reduced, at a rate n M 1 , by the tumour cells that secrete pro-tumour cytokines (e.g., IL −10, TGF-β) . All other rates that control the dynamics of M1 macrophages are as described in Section 2.
• The recruitment of M2 cells at the tumour site is helped by cytokines (e.g., IL −10) and chemokines (e.g., CCL2) produced by the tumour cells (Solinas et al., 2009). We denote this recruitment rate by n M 2 . For simplicity, throughout this study we consider n n = . All other rates that control the dynamics of M2 macrophages are as described in Section 2.
We emphasise that in model (3), we incorporated only an example of tumour-macrophage-Th cell interactions. Continuous development of this research area, will likely reveal more types of interactions among these cells. However, it is not the goal of this article to model detailed dynamics of tumour-immune interactions. Rather, we plan to investigate whether the assumptions incorporated in (3) can explain the paradoxical anti-tumour and pro-tumour immune dynamics observed experimentally in B16 melanoma cells (as discussed in Section 1).
We also note that while there are many other types of tumour growth laws (e.g., exponential, power, von Bertalanffy, Gompertz or sub-linear) that can fit various experimental data sets, recent studies suggest that the most appropriate growth laws seem to be dependent on the details of the experiments and on the particular tumour cell lines (Murphy et al., 2016;Sarapata and de Pillis, 2014;Benzekry et al., 2014;Talkington and Durrett, 2015). Since the goal of this study is not to compare in detail our results to various experimental data sets, we decided to focus only on one law, the logistic growth, and to investigate whether this assumption on tumour growth can help propose some generic biological mechanisms that can explain the apparent paradox in the observed anti-tumour immune responses.
Before investigating the dynamics of system (3), we note that (3) has non-negative solutions (see the discussion in Appendix B).

Steady states and stability
Next, we study the long-term behaviour of model (3), when the system is at equilibrium. The existence of four possible equilibrium points (listed below) emphasises the complexity of (3).

No tumour cells and no immune cells:
T H H ( * , * , * , For the parameter values shown in Table A.1, there are three such steady states that are real and positive (see Fig. E.18 in Appendix E), and their stability is illustrated in Fig. E.19(e)-(g). As in Section 2.1, we are now interested in investigating the parameter space where tumour growth and elimination occurs in the presence of a type-I immune response, a type-II response or a mixed response. Thus, we focus on the two steady states with non-zero immune responses. Since the steady state H H M M (0, * , * , * , * ) 1 2 1 2 is similar to the state investigated in Fig. 5, we can conclude that tumour elimination can occur in the presence of a type-I response, a type-II response, or a mixed type-I/type-II response. We will return to this aspect in Section 3.2, when we will investigate numerically the long-term dynamics of system (3).
For the tumour-immune coexistence state T H H M M ( * , * , * , * , * ) 1 2 1 2 , let us first investigate the parameter values for which T * > 0, which is equivalent (from (4a)) with solving the following equation for T * > 0: values investigated in this study (panels (a), (a')), whenever tumours grow they are accompanied by a type-II immune response. However, for very large g M 2 values, tumours can exist also for M > 1 and H > 1 (see panels (b), (b')). This result suggests that there could be fewer M2 cells compared to M1 cells, but if these cells secrete large amounts of type-II cytokines, they can skew the tumour microenvironment in favour of tumour sustenance and growth. (We will return to this hypothesis in the Discussion section.) We also need to emphasise here that an increase in g H 1 (from 4.2 × 10 −9 in panels (a), (b), to 1.26 × 10 −7 in panels (a'), (b')) reduces the parameter space over which we can expect tumour-growth in the presence of a type-I response.
In Fig. 9 we notice that the 10-fold increase in g M 1 (from g = 6 × 10 (i) forces T * > 0 to exist mainly during a type-II response, and (ii) induces the requirement for much higher g M 2 values for tumour persistence in the presence of a type-I response (i.e., at least a 150fold increase in g M 2 ; see panels (b), (b'), where only a mixed type-I/ type-II response was obtained after a 126-fold increase in g M 2 ). We emphasise here that small changes in g H 2 do not have a significant effect on tumour growth (also supported by the sensitivity analysis in Fig. 14). To observe a difference between the diagrams in panels (a), (b) and those in panels (a'), (b') we had to increase g H 2 by more than 40-fold (shown in panels (a'), (b') is the effect of a 60-fold increase in g H 2 ). In this case, the increase in g H 2 affected mainly the region where (see the right figures in panels (b), (b')). Overall, Figs. 8 and 9 suggest that the parameters most likely to impact tumour growth/decay are g H 1 , g M 2 and g M 1 . We will return to this aspect in Section 3.3, when we will perform a sensitivity analysis for the transient dynamics of model (3).

Short-term and long-term dynamics
To investigate numerically the long-term dynamics of immune cells and cancer cells, we use the parameter values described in Table A Tumour elimination. First, we focus on the parameter ranges for r M1 and r M2 that ensure tumour elimination in the presence of a type-I immune response, a type-II immune response, or a combination of both type-I/type-II immune responses. In this case, the dynamics will approach the stable steady state H H M M (0, * , * , * , * ) 1 2 1 2 , and the dominant immune responses are consistent with those in the bifurcation diagram shown in Fig. 5. We emphasise this aspect by discussing separately the following two cases involving the activation rates a a , H H 1 2 for the Th1 and Th2 cells: (1) Case a a < H H 1 2 . Fig. 10 illustrates the short-term dynamics (panels (i); t < 30 days) and long-term dynamics (panels (ii); t ≤ 100 days) . In panel (a)(i) we observe a double switch between the M1 and M2 cells that dominate the dynamics (and this is associated with only one switch in the Th1-Th2 dynamics). In panel (b)(i) we observe a single switch the dynamics of both M1 and M2 cells, and Th1 and Th2 cells. In all cases the tumour is eliminated, and the results are consistent with the bifurcation diagrams in Fig. 5(a).
(2) Case a a > H H 1 2 . Fig. 11 illustrates the short-term dynamics (panels (i)) and long-term dynamics (panels (ii)) of model (3)  . In panel (a)(i) we observe a double switch between the M1 and M2 cells that dominate the dynamics (but this is not associated with any switch in the Th1-Th2 dynamics). In panel (b)(i) we observe a double switch between the Th1 and Th2 cells that dominate the dynamics (associated with a single switch in the M1-M2 dynamics).
Tumour persistence. Fig. 12 2 ). By investigating the short-term dynamics of model (3) (see panel a(i)) we observe a switch in both the Th1-Th2 and M1-M2 dynamics, from an initial type-I response to a later type-II response. This is consistent with the bifurcation results in Fig. 8( ). Moreover, we would like to emphasise that the dormant behaviour exhibited by the tumour for t ∈ (5, 15) is mainly the result of a very large M1 population that keeps the tumour under control. As soon as this M1 population is reduced, the tumour grows fast towards its carrying capacity.
The difference between tumour dormancy/growth in Fig. 12 and tumour elimination in Figs. 10 -11 is the result of (a) a small change in the rate at which tumour cells are eliminated by the Th1 cells via the cytokines they produce (from g = 4.4 × 10 for tumour growth). However, different other combinations of parameter changes can lead to similar tumour dormant behaviours (which seem to be controlled by relatively high levels of M1 cells). To investigate the effect of small changes in parameter values on the level of tumour and immune cells during dormancy (not only M1 but also M2, Th1 and Th2 cells), in Section 3.3 we will perform a sensitivity analysis.
In Section 1 we mentioned the experimental results in Chen et al. (2011) (see also Fig. 2(c)), which showed tumour growth being associated with a shift in the ratio of M1 and M2 cells: from M1:M2 ≈ 90:10 on day 7, to M1:M2 ≈ 20:80 on day 14. To compare these experimental results with our numerical results, in Fig. 12(b) we show the percentage of Th cells and macrophages on day 4.5 (when tumour is small), day 14 (when tumour is dormant) and day 19 (when tumour approaches its carrying capacity). We see that tumour growth is not only associated with an increase in the percentage of M2 cells (as shown experimentally in Chen et al. (2011)), but also with an increase in the percentage of Th2 cells (as shown experimentally in Protti and Monte (2012)). Note that this is one possible outcome of the model. Changes in parameter values could lead to different ratios of Th2:Th1 cells and M2:M1 cells as tumour progresses.
Finally, we recall that the results in Fig. 8(b) suggested that by increasing g M 2 one could observe tumour existence also in the case of a  Table A.1.

R. Eftimie, H. Hamam
Journal of Theoretical Biology 420 (2017) 82-104 type-I immune response with M H , > 1 (in addition to a type-II response, with M H , < 1). We show in Fig. 13 the short-term and long-term dynamics of model (3), characterised by the persistence of tumour cells at lower values (with a maximum of about 5×10 7 cells). This persistence is the result of a type-I immune response, which alternates for short periods of time with a type-II response. We emphasise that these oscillations in tumour growth/decay (triggered by oscillations in the type-I/type-II immune responses) might not be always observable in a clinical setting. Friberg and Mattson (1997) showed that in humans, the tumour diagnostic level is between 10 7 and 10 9 cells. Therefore, 5×10 7 cells might not be always detected clinically.

Sensitivity analysis
Since the majority of parameter values could not be approximated from the literature, in the following we perform a sensitivity analysis to investigate the effect of changes in these parameters on the growth of the tumour. To this end, we vary each parameter P by ±10% or ±90% at a time (i.e., P P ± ▵ , with P P ▵ = 0.1 or P P ▵ = 0.9 ), and investigate the impact of this change on tumour size on day 10 (an arbitrarily-chosen day, when the tumour has not reached its maximum size yet). The relative change in tumour size on day 10 (i.e., T ▵ (10)) is used in Fig. 14 to plot the ratio of relative changes: ( )/( ) |▵ | | | . Fig. 14 illustrates tumour sensitivity to changes in the parameter values: (a) by ±10% and (b) by ±90%. The parameters that have the most significant effect on tumour size when varied by ±10% are: the tumour growth rate (α), the proliferation of Th1 cells (p H 1 ), the elimination rate of tumour cells by the Th1 cells (g H 1 ) and by M1 macrophages (g M 1 ), the carrying capacity of Th cells (m 1 ), the carrying capacity of macrophages (m 2 ), the transition rate from M1 to M2 cells (r M 1 ), the activation rate of M1 cells (a M 1 ) and the proliferation of M2 cells in the presence of type-II cytokines (p M 2 ). It is likely that p M 1 might also have higher impact on tumour if we would consider higher selfproliferation rates for M1 cells. The parameters that have the most significant impact on tumour size when varied by ±90% are p H 1 and α (similar to case (a)). Also a decrease in m 1 , g H 1 , g M 1 and r M 1 leads to a significant increases in tumour size (see the inset in the right panel of Fig. 14(b)). (Note that, in Fig. 14(b) is difficult to see the reduction in tumour size as we vary the parameter valuesbecause of the very large increases in tumour size.) We also need to emphasise that g H 2 and g M 2 (both associated with a type-II immune response) do not have a significant impact on tumour reduction. This is a particularly interesting result that might be of biological interest, since at least g H2 has the same order of magnitudesee Table A.1as parameters g H 1 and g M 1 (which have a significant effect on tumour reduction/growth). Moreover, this result supports the idea that the elimination of tumour cells by the Th2 cells in Mattes et al. (2003) was not the result of direct Th2-tumour interactions (via Th2-cytokines), but the combined effect of different anti-tumour cells.
To gain a better understanding on tumour dormancy (and on the role of immune response in controlling tumour growth), next we perform a tumour and immune sensitivity to small changes in four parameter values associated with anti-tumour/pro-tumour immune responses: g H 1 , g H 2 , g M 1 , g M 2 . To this end, we start with the baseline parameters that lead to tumour dormancy/growth in Fig. 12(a), and we vary them by ±10% to investigate the changes in tumour and immune sizes at day t=10 (when dormancy occurs). First, we note that during tumour dormancy, changes in parameter g H 1 have a slightly bigger impact on tumour at day t=10 (T (10)) compared to changes in parameter g M 1 -see Fig. 15(a). This is in contrast to the case of tumour elimination (see Fig. 14(a), left panel) where g M 1 has a bigger impact on T (10) compared to g H 1 . Second, we note that during tumour dormancy g M 2 has a stronger impact on T (10) (see Fig. 15(a)) compared to the case of tumour elimination where g M 2 barely affects T (10) (see Fig. 14(a)). In fact, we observe that ±10% changes in the three parameters g H 1 , g M 1 and g M 2 , lead to changes of relatively similar magnitudes in tumour cells ( Fig. 15(a)), and in each of the four types of immune cells ). This suggest that tumour dormancy is the result of a delicate balance between the anti-tumour effect of Th1 and M1 cells, and the pro-tumour effect of M2 cells. Moreover, by looking at panels (b)-(e) we observe that the effects of g M 1 and g M 2 do not balance perfectly during dormancy: g M 2 causes slightly larger effects in both tumour and immune responses compared to g M 1 (and this imbalance eventually translates into tumour relapse).
To conclude the discussion on the effects of parameters g H 1 , g M 1 and g M 2 on the immune responses during tumour dormancy, we stress that while it was expected that an increase in g M 1 and g H 1 would be associated with an increase in M 1 and H 1 (through the direct reduction of tumour), it was however unexpected that g H 1 would have an effect on H 2 and M 2 cells (stronger than the effects of parameters g H 2 and g M 2 ).  For the rest parameters values see Table A.1. In this case, the tumour persists being controlled alternatively by a type-I and a type-II immune response.

Summary and discussion
In this article, we derived two mathematical models for the dynamics of immune responses involving Th1 & Th2 and M1 & M2 cells, in the absence and in the presence of tumour cells. We then used these models to propose mechanistic hypotheses that could explain the contradictory results in the experimental data for the immune response against melanoma B16 cells.
We started with a model that considered only the interplay between M1 and M2 macrophages, and Th1 and Th2 cells in response to some external pathogen that first triggered an M1 response (i.e., M (0) > 0 1 ). To shed light on the complexity of model dynamics, we first calculated the steady states (to study the long-term behaviour of the model) and then we performed numerical simulations for the short-term and longterm model dynamics. By focusing on the ratio r r / M M 1 2 (of macrophages re-polarisation rates), and the activation rates of Th cells (a H1 , a H2 ) in the presence of signals received from macrophages, we were able to classify the immune responses into: a type-I dominated response  Barros et al. (2013) (Table 1), the authors showed that about 60.7% of Th1 disease cases investigated (in the context of infectious mononucleosis and Crohn's disease) have M1>M2, and about 72.5% of Th2 disease cases investigated (in the context of allergic nasal polyps, oxyuriasis, wound healing and foreign body granulomas) have M2 > M1. Thus their results suggest that there are Th1 diseases with a higher level of M2 cells, and Th2 diseases with a higher level of M1 cells (consistent with our numerical results).
Next, we generalised the mathematical model to consider also tumour dynamics. We showed numerically that tumour elimination can occur both in the presence of a type-I dominated immune response, as well as in the presence of a type-II dominated response (as observed experimentally in Mattes et al. (2003), Xie et al. (2010), Kpbayashi et al. (1998; see also Figs. 1 and 2). We need to emphasise that tumour elimination also required a relatively large tumour lysis rates g H 1 and g M 1 and a low g M 2 . As before, the type of immune response that dominated the dynamics was decided by the ratio r r / M M 1 2 and the activation level of immune cells (a H1 , a H2 ).
Tumour growth towards carrying capacity (or some very large size) was always associated in our study with a long-term type-II immune response, i.e., H H > for t ≤ 10 days) was always replaced in the long-term by a type-II immune response. This shift from a type-I to a type-II response was observed also experimentally in the context of cancer growth. For example, Chen et al. (2011) showed a 90:10 ratio of M1:M2 macrophages in B16F10 melanoma tumours around day 7, and a 20:80 ratio of M1:M2 macrophages around day 14 (see Fig. 2(c)). Other experimental studies have described a shift from a Th1 response to a Th2 response during the first 14-20 days of progression of malignant tumours (see Tatsumi et al., 2002 for human melanoma). These experimental studies also suggested that one could improve cancer outcome by re-polarising the macrophages and Th cells from a type-II response associated with tumour growth to a type-I response associated with tumour decay (Heusinkveld and van der Burg, 2011). Our theoretical results are in agreement with the experimental suggestion that a type-I response improves long-term cancer outcome. Moreover, our results also emphasise the complexity of the tumour-immune system, in which a type-I immune response might alternate with a type-II immune response (for short-term or long-term), thus leading only to tumour control but not tumour elimination.
We stress that the interaction between the pro-tumour/anti-tumour effects of macrophages and Th cells affects tumour dynamics in a nonlinear manner. For example, a 10-fold increase in the rate of tumour clearance by M1 macrophages (g M 1 ) caused tumour persistence only in the presence of a type-II immune response (i.e., a type-I immune response would be associated to tumour clearance). To ensure tumour persistence also in the presence of a type-I response, the 10fold increase in g M 1 needed to be counter-balanced by at least a 150-fold increase in the tumour growth rate in the presence of M2 cells, g M 2 (see Figs. 8 and 9). This nonlinearity in the anti-tumour response is likely the result of the interplay between the macrophages and the Th cells, an aspect not very well studied at experimental level. Although there are some studies on the interactions between macrophages and CD4 + T cells, for example, in the context of breast and lung cancer (DeNardo et al., 2009;Almatroodi et al., 2016), or in the context of rheumatoid arthritis (Roberts et al., 2015), such studies do not shed much light on the nonlinear interactions between these different types of immune cells.
In the context of the anti-tumour effect of macrophages, the sensitivity analysis in Fig. 14(a) suggested that tumour elimination was mainly the effect of M1 macrophages (and to a lesser extent the effect of Th1 cells). This is an interesting hypothesis generated by the model, which, if validated experimentally, could influence the current anti-tumour immune therapies that focus mainly on T cell responses Voena and Chiarle, 2016). In contrast, the sensitivity analysis in Fig. 15 suggested that the transient decrease in tumour size on day 10 during tumour dormancy was mainly the effect of Th1 cells (and to a lesser extent the effect of M1 cells). In fact, the tumour dormant behaviour was the result of a delicate balance between the anti-tumour responses of Th1 and M1 cells, and the pro-tumour responses of M2 cells. In addition, the results in Figs. 8 and 9 suggested that the three parameters, g H 1 , g M 1 and g M 2 , influenced also the asymptotic behaviour of model (3). This is in support of the idea that anti-cancer immunotherapies should focus on the combined effect of T cells and M1 macrophages.
The results in Fig. 8 suggested that there could be very few M2 cells (and many M1 and Th1 cells), but if these M2 cells secrete large amounts of type-II cytokines (i.e., large g M 2 ), they can skew the tumour microenvironment in favour of tumour sustenance and growth. This would support the experimental results in Mattes et al. (2003), where a type-I environment was not enough to eliminate B16F10 melanoma cells. The authors in Mattes et al. (2003) recognised that the inability of Th1 cells to eradicate tumours might have been influenced by the presence of pro-angiogenic tumour-infiltrating macrophages (i.e., M2 cells), but they did not measure the levels of M2 and M1 macrophages, nor the levels of Th1 and Th2 cells. In fact, Mattes et al. (2003) identified the Th1 and Th2 immune responses by the levels of type-I and type-II cytokines produced by these cells: high IL-5, IL-13 and IL-4 for a Th2-dominated response, and high IFN-γ, TNF-α and IL-13 for a Th1-dominated response (note here the relatively high levels of IL-13 observed during both Th1 and Th2 responses; and the fact that IL-13 is also involved in the alternative activation of M2 macrophages (Martinez and Gordon, 2014)). Since many experimental studies focus on the levels of cytokines as a proxy for the number of immune cells corresponding to a type-I or type-II response (Mattes et al., 2003;Almatroodi et al., 2016), to be able to test our hypothesis regarding the role of g H 2 and M2 cells on tumour persistence during type-I responses, we need to extend model (3) by incorporating explicitly the effects of type-I and type-II cytokines on tumour-immune interactions (i.e., an approach similar to Eftimie et al. (2010), where a mathematical models incorporated the effects of type-I, type-II, tumour-promoting and tumour-suppressing cytokines). In this study, to keep the models relatively simple, we ignored deliberately the microenvironment which can alter the immune response against cancer (Hanahan and Weinberg, 2011). However, the incorporation of the explicit effects of type-I and type-II cytokines (which can be further altered by the tumour cells (Burkholder et al., 2014)) would allow us not only to compare our results with available experimental cytokine data, but also to gain a better understanding of how to control cell-cell communication (by controlling cytokine signalling) with the ultimate goal of improving cancer immunotherapies.
Note from Table A.1 that models (1) and (3) contain both fast and slow variables. One could have used a quasi-steady state analysis to simplify the models. However, such an analysis might lead to limitations in our understanding of the transient dynamics of the Th1-Th2 and M1-M2 cells (see, for example the study in Flach and Schnell (2006)). This type of transient dynamics was observed in experimental studies on early tumour behaviours, which suggested that the ratios of Th1/Th2 cells or M1/M2 cells can be used as independent predictive markers of patient survival (Monte et al., 2011;Protti and Monte, 2012;Chen et al., 2011). In this theoretical study we showed that these ratios of immune cells can change once or twice before they stabilise towards a steady state (and they stabilise when the tumour reaches either a very large size or is eliminated; see Figs. 10-12). The changes in the dominating Th or macrophages dynamics are not always correlated with each other. Moreover, we showed the possibility of having a long-term oscillatory tumour-immune dynamics characterised by low tumour values and periodic changes between type-I and type-II immune responses; see Fig. 13. While sustained periodic tumour oscillations are not very often observed in clinical studies (although see Gliozzi et al., 2010), we emphasise that model (3) exhibits such oscillations for tumour sizes around the detection threshold (of about 10 − 10 7 8 cells (Friberg and Mattson, 1997)). This suggest that oscillations between type-I and type-II immune responses (in the presence of tumour) might be more common in clinical/experimental settings but they might not be measured since the tumour cannot be detected. Overall, we hypothesise that trying to predict the long-term outcome of the tumour while the ratios Th1/Th2 and M1/M2 are still varying due to the cross-talk with the tumour environment, might not always offer accurate predictions on patient survival.
At a more theoretical level, it would be interesting to investigate the differences between the double feedback in tumour-immune dynamics modelled in this study, and a single feedback for tumour-immune interactions. Such an investigation (to be the subject of a future study) would allow us to uncover the minimal biological mechanisms that need to be incorporated into a model to explain the dominant type-I and/or type-II immune responses associated with cancer immunotherapies.
Finally, these numerical results for systems (1) and (3) have generated two new mathematical questions that will be answered analytically in future studies: (i) analytical investigation of fast and slow parameters that control transient and long-term tumour-immune behaviours, and how the simplified dynamics in the slow/fast models matches the original dynamics; (ii) analytical investigation of the Hopf bifurcation that generated the limit cycle shown in Fig. 13.
Biological realism of the parameter values and overall results. The results of this study depend on the parameter values described in Table  A.1. Some of these values were taken from the literature, others were approximated based on published experimental results, and the remaining values were varied within some estimated ranges (see Appendix A). This approach is very common in the mathematical immunology literature, due to a lack of quantitative results regarding the immune responses following various antigen stimulations. In addition to the fact that very few labs measure and estimate kinetic parameters (the majority of such studies focusing on lymphocyte kinetics following pathogen stimulation; see for example Borghans and Boer, 2007;Asquith et al., 2009;Boer and Perelson, 2013), there is also the difficulty of interpreting kinetic data; see the review in Boer and Perelson (2013). Moreover, the few rigorously estimated kinetic parameters in the mathematical immunology literature depend on the estimation method used, as emphasised in Laydon et al. (1675). A more detailed discussion on model validation and parameter estimation in mathematical immunology can be found in Eftimie et al. (2016).
Based on these facts, we acknowledge that the majority of models in the mathematical immunology literature, including this particular study, can have at this moment only a theoretical value. In particular, the model presented here can only propose hypotheses regarding the possible outcomes of the interactions between the Th1-Th2 and M1-M2 immune responses, in the absence/presence of tumour cells.
We showed that small variations in the values of parameters that control tumour cells lysis via anti-tumour cytokines (e.g., g H 1 , g M 1 , g M 2 ), or the parameters for the activation of Th cells (a a , H H 1 2 ), or the macrophages re-polarisation rates (r r ,

M M
1 2 ) could explain the variety of tumour-immune dynamics observed in the experimental literature. To obtain a better understanding of immune responses to specific diseases, the next step would be to quantify the rates that control various type-I and type-II immune responses. Therefore, for a better mechanistic understanding of the in vivo immune responses, which can be obtained with a more realistic in silico model, mathematicians (and immunologists) need to have access to relevant experimental data that could then be used to parametrise the mathematical models. The goal of our present study was not to parametrise the models to specific diseases, but to propose some general hypotheses regarding the processes involved in different immune responses. (but we note that we could have chosen any other parameter). The four symbols in Fig. E.16 show the real parts of the four eigenvalues corresponding to the Jacobian matrix (E.1). Numerical calculations of the eigenvalues corresponding to the steady state (0, 0, 0, 0) show that this state is stable for the parameter values shown in Table A