A mathematical model of cytotoxic and helper T cell interactions in a tumour microenvironment

Abstract We develop a mathematical model to examine the role of helper and cytotoxic T cells in an anti-tumour immune response. The model comprises three ordinary differential equations describing the dynamics of the tumour cells, the helper and the cytotoxic T cells, and implicitly accounts for immunosuppressive effects. The aim is to investigate how the anti-tumour immune response varies with the level of infiltrating helper and cytotoxic T cells. Through a combination of analytical studies and numerical simulations, our model exemplifies the three Es of immunoediting: elimination, equilibrium and escape. Specifically, it reveals that the three Es of immunoediting depend highly on the infiltration rates of the helper and cytotoxic T cells. The model’s results indicate that both the helper and cytotoxic T cells play a key role in tumour elimination. They also show that combination therapies that boost the immune system and block tumour-induced immunosuppression may have a synergistic effect in reducing tumour growth.


Introduction
Cancer is a disease that affects millions of people worldwide and whose impact continues to increase (Siegel, Miller, & Jemal, 2015). Understanding how and why cancers develop, and how they may be effectively treated, is one of the major societal challenges of the twenty-first century. Cancer is distinguished from other diseases by its ability to hijack normal cell functions, to initiate uncontrolled cell proliferation and, of particular interest here, to evade immune detection and elimination (Hanahan & Weinberg, 2000. T cells are one of the most important components of the immune system in the fight against cancer. They originate from pluripotent haematopoietic stem cells in the bone marrow which migrate to the thymus where they mature into naive T cells. The naive T cells move to the lymph nodes where they become activated on contact with their cognate antigens. Activated T cells proliferate rapidly to produce a substantial army of antigenspecific T cells. These short lived T cells are then transported through the blood vessels to the tumour where they bind to and kill infected cells and also produce cytokines that recruit other immune cells to the tumour. This process continues until either the tumour has been removed or the tumour adapts to, and evades, targetting by the T cells (Chen & Mellman, 2013;Janeway, Travers, Walport, & Shlomchik, 2001).
There are two populations of T cells, helper and cytotoxic, which are distinguished by their expression of CD4 and CD8 proteins, respectively. Naive helper T cells are activated by antigen presented with a major histocompatibility complex (MHC) class II molecule. (A glossary of the terminology used in this article is provided in Appendix 1.) Naive cytotoxic T cells are activated by antigen presented with an MHC class I molecule. Once activated the helper and cytotoxic T cells perform complementary functions to eliminate the tumour. Helper T cells further differentiate into subpopulations classified by the specific cytokines that they produce. In this way, they regulate multiple aspects of an immune response. For example, they promote the proliferation of cytotoxic T cells, they recruit and promote cells of the innate immune response, and they control levels of inflammation at the tumour site (Hung et al., 1998;Magombedze et al., 2013). (In this work, we do not distinguish these subpopulations.) Cytotoxic T cells scan the body for transformed cells (cancer cells) to which they bind before inducing cell killing.
Most experimental and clinical studies have focused on the role played by cytotoxic T cells in tumour elimination (Fridman, Pages, Sautes-Fridman, & Galon, 2012) and neglected the role of helper T cells (Hiraoka et al., 2006;Rad, Ajdary, Omranipour, Alimohammadian, & Hassan, 2015). Recently, several authors have argued that helper T cells form a vital component of an anti-tumour immune response (Badoual et al., 2006;Magombedze et al., 2013). In this work, we focus on the complementary roles played by helper and cytotoxic T cells as cytokine producers and tumour cell killers, respectively. For simplicity, we do not consider naive and memory helper and cytotoxic T cell populations as in Robertson-Tessi, El-Kareh, and Goriely (2012). Our focus is to identify which T cell population should be targeted to boost the anti-tumour actions of helper and cytotoxic T cell populations in the tumour micro-environment.
As first hypothesised by Ehrlich (1909), and subsequently confirmed experimentally, the immune system can detect and eliminate cancer before it reaches clinically detectable sizes. During 'immunosurveillance', the immune system detects and destroys cancerous cells. In some cases, the tumour evades the immune system and develops into a large mass, which may eventually prove fatal. The theory of 'immunoediting' proposed by Dunn, Old, & Schreiber (2004) represents the possible outcomes of tumour-immune interactions, via the three E's of immunoediting -elimination, equilibrium and escape.
While immunosurveillance should, in theory, rid the body of pre-cancerous and/or cancerous cells, in practice this does not always occur. If, for example, immunogeneic cells are removed by the immune system then a poorly immunogeneic and antigenic population of tumour cells remains (DuPage, Mazumdar, Schmidt, Cheung, & Jacks, 2012;Pardoll, 2003;Vicari, Caux, & Trinchieri, 2002). As these tumour cells proliferate, they mutate, creating cells that present low levels of tumour antigen and, consequently, evade recognition by T cells (Vesely & Schreiber, 2013). Cancer cells also manipulate the local microenvironment by producing immunosuppressive cytokines such as TGF-β and IL-10 that promote differentiation of helper T cells to a regulatory phenotype (Kawamura, Bahar, Natsume, Sakiyama, & Tagawa, 2002). Cancer cells may also produce cytolytic molecules such as Granzyme B that induce apoptosis of immune cells (Igney & Krammer, 2002).
A number of novel immunotherapies have been designed to either enhance/boost the immune response or reverse/block the immunosuppressive effects of the tumour. In this work, we focus on adoptive T cell transfer therapy. This therapy acts to boost the army of T cells specific for a chosen tumour antigen (Rosenberg, Restifo, Yang, Morgan, & Dudley, 2008;Ruella & Kalos, 2014). Here circulating T cells are extracted from the host and those that are specific to the targeted tumour antigen are selected. This pool of antigen specific T cells are then stimulated to proliferate without the constraints of an immunosuppressive environment. Once a substantial population of these T cells have amassed, they are re-infused into the host. With the identification of an appropriate tumour antigen to target, this therapy holds the potential to elicit a potent anti-tumour immune response. Most therapies have been designed to target the cytotoxic T cells which are able to recognise and directly kill tumour cells. Work needs to be done to explore the potential of targeting the helper T cell population instead of or in addition to the cytotoxic T cell population (Zanetti, 2015). In this work, we explore which T cell population should be the target of immunotherapies designed to boost the T cell populations in the tumour micro-environment.
To understand how the immune system works to eradicate cancer, it is necessary to identify its key components. Several mathematical models have been proposed, typically involving systems of time-dependent ordinary differential equations (ODEs). For example, Eftimie, Bramson, and Earn (2011) provides a comprehensive review, covering models of tumour growth (single equations) through to models consisting of many cell types (e.g. innate and adaptive immune cells) and molecules (e.g. cytokines). Many existing ODE models view tumour-immune interactions as a predator-prey system, where the tumour cells are the prey and the immune cells are the predators (Frascoli, Kim, Hughes, & Landman, 2014;Kuznetsov, Makalkin, Taylor, & Perelson, 1994;López, Seoane, & Sanjuán, 2014;Piotrowska, Bodnar, Poleszczuk, & Foryś, 2013;Wilkie & Hahnfeldt, 2013). With two dependent variables it is possible to characterise both the local and global stability of steady states through phase-plane analysis, and the bifurcation structure of the system as model parameters vary. In this way, it is possible to identify those parameters that have the strongest influence on tumour escape.
To explore the role of helper T cells, in addition to cytotoxic T cells and tumour cells, further variables must be considered. This population has been included both implicitly and explicitly. Kirschner and Panetta (1998) proposed a system of three ODEs for the cytokine IL-2 produced by the helper T cells, cytotoxic T cells and tumour cells. The inclusion of the cytokine IL-2 can be thought of as implicitly accounting for the helper T cell population. They found that the tumour antigenicity is a critical parameter determining whether tumour escape, equilibrium or elimination occurs. They also identified stable periodic solutions, mirroring clinical observations of tumour suppression and regrowth. On the other hand, Robertson-Tessi, El-Kareh, and Goriely (2012) proposed a detailed mathematical framework (twelve ODEs) describing interactions between multiple cell types and cytokines involved in an immune response to cancer, capturing the three E's of Immunoediting (Dunn, Old, & Schreiber, 2004). Helper T cells were explicitly included and further subdivided into three distinct subpopulations based on their maturity status (e.g. naive, effector, memory). Their aim was to determine which aspect of the immune response is dominant at different stages of tumour development, and which immunosuppressive effects facilitate tumour escape. However due to its complexity, investigation of the model was limited to numerical simulations rather than analytical identification of the steady states and their linear stability.
Delay differential equation models have also been used to account for the time delay between initial encounter of a naive T cell with its matching antigen and the effector response (Dong, Huang, Miyazaki, & Takeuchi, 2015;Rihan, Rahman, Lakshmanan, & Alkhajeh, 2014). In this paper, we follow Kuznetsov et al. (1994), and adopt an ODE framework to study the dynamics associated with interactions between tumour cells, helper T cells and cytotoxic T cells.
In this paper, a simplified ODE model accounting for the interactions of tumour cells with helper and cytotoxic T cells is developed. This model extends the pioneering study of Kuznetsov et al. (1994) by including helper T cells and accounting implicitly for immunosuppressive effects. The proposed model also draws on key components of the more complex model developed by Robertson-Tessi et al. (2012) that included both helper and cytotoxic T cells in an anti-tumour immune response, and examined the impact of different effects of immunosuppression at various stages of tumour development. Due to the complexity of their model, the analysis was limited to a numerical investigation. The simpler model proposed here has the advantage of being analytically tractable while retaining both a helper T cell population and implicit forms of immunosuppression. This permits the identification and characterisation of the stability of steady states as a function of model parameters. This also allows for the derivation of an amplitude equation to describe time-dependent oscillatory solutions observed from numerical integration of the model equations. As clinical trials have shown that the level of tumour infiltrating T cells correlates with patient outcomes (Lee et al., 1989;Kawai et al., 2008;Hiraoka et al., 2006), we focus on understanding how changes in the rates of infiltration of cytotoxic and helper T cells affect the model dynamics. This model reveals that the three Es of immunoediting depend highly on the infiltration rates of the helper and cytotoxic T cells. The results identifying where in parameter space tumour elimination, equilibrium and escape occur are used to identify conditions under which cytotoxic and/or helper T cells should be targeted by adoptive T cell therapy. We find that the parameters associated with immunosuppression play a decisive role in tumour-elimination.
The remainder of this paper is organised as follows. We formulate our mathematical model in Section 2, while in Section 3, we present simulation results that show how our model captures the three E's of immunoediting (elimination, equilibrium and escape) and that it may exhibit oscillatory solutions which we view as equilibrium solutions since the tumour mass remains bounded at all times. In Section 4, we identify the steady state solutions, and indicate those regions of parameter space in which multiple steady state solutions exist. In Section 5, we characterise the linear stability of the steady state solutions and identify a Hopf bifurcation from which a branch of oscillatory solutions emerges from a steady-state solution branch as the rate of infiltration of the cytotoxic T cells is varied. In Section 6, we confirm the existence and stability of the associated limit cycle and determine its amplitude by performing a weakly nonlinear analysis in a neighbourhood of the Hopf bifurcation. In Section 7, we perform a parameter sensitivity analysis, focusing on how the tumour's growth dynamics change as parameters controlling the immunosuppressive effects of the tumour are varied. The paper concludes in Section 8 where we summarise the key results, discuss their therapeutic implications and outline directions that merit further investigation.

Model development
We propose a system of three, time-dependent, ordinary differential equations (ODEs) to model the interactions between cytotoxic and helper T cells and tumour cells in a wellmixed (i.e. spatially uniform) tumour microenvironment, see Figure 1. With t denoting time, the dependent variables N(t), T H (t) and T C (t), represent numbers of tumour cells, and helper and cytotoxic T cells respectively.
We assume that the evolution of the tumour cells is dominated by proliferation and cell death due to interactions with cytotoxic T cells. For simplicity, we assume that, in the absence of an immune response, the tumour undergoes logistic growth, with carrying capacity κ (units: number of cells) and growth rate γ (units: day −1 ). We interpret solutions for which the system approaches its carrying capacity (N ≈ κ) as tumour escape. We assume further that, once an immune response has been stimulated, the tumour cells interact with cytotoxic T cells at a rate which is proportional to the product of their cell numbers, with constant of proportionality k (units: number of cells −1 day −1 ), and that such interactions lead to tumour cell death, with probability p, or inactivation of the cytotoxic T cell, with probability (1 − p). Combining these processes, we propose the following ODE for the time evolution of the tumour cells: We assume that the evolution of the helper and cytotoxic T cells is dominated by infiltration from the lymph nodes, proliferation in response to the tumour cells and natural cell death. We assume that the two T cell populations infiltrate the tumour from the lymph nodes at constant rates σ H and σ C (units: number of cells day −1 ), and die at rates δ H and δ C (units: day −1 ). Naive T cells are regularly produced by hematopoietic stem cells and therefore a constant population of these naive T cells are present in the circulating blood stream. We can assume that the flow of blood is constant and therefore, in addition to direct stimulation of T cells in the presence of tumour antigen, there is a constant influx of primed T cells from the circulating blood. Furthermore, this assumption enables us to implicitly account for the presence of a constant pool of memory T cells present in the absence of antigen stimulation (Farber, Yudanin, & Restifo, 2014). This approach has been widely adopted in similar models to account for a background level of circulating T cells (de Pillis et al., 2005;Eftimie et al., 2011;Kuznetsov et al., 1994).
We suppose that the tumour suppresses the proliferation of the helper T cells and, thereby, indirectly inhibits the cytotoxic T cells. Specifically, we assume that the helper T cells proliferate at a rate which has biphasic dependence on the number of tumour cells, N(t), with constant of proportionality α (units: day −1 ). The parameterÑ (units: number of cells) is the number of tumour cells at which helper T cell proliferation is half-maximal. The proliferation rate (per T H cell) attains a maximum value of α/2 when N =Ñ and decays to zero as N → ∞. A similar biphasic term was adopted by Serre et al. (2016). Through their production of cytokines (e.g. IL-2, IFN-γ ), helper T cells promote the proliferation of cytotoxic T cells (Akhmetzyanova et al., 2013;Rad et al., 2015). Rather than introducing additional dependent variables to represent these cytokines, for simplicity, we account for their effects by assuming that cytotoxic T cells proliferate at a rate proportional to the number of helper T cells, with constant of proportionality β (units: number of cells −1 day −1 ). Combining these assumptions, we propose the following ODEs to describe the time evolution of the helper and cytotoxic T cells respectively: We close equations (1)-(3) by prescribing the following initial conditions: wherein the non-negative constants N 0 , T H0 and T C0 denote the initial concentrations of the tumour cells, helper T cells and cytotoxic T cells respectively.

Parameter values
In the absence of suitable experimental data, we base our estimates of the parameters values in Equations (1)-(3) on those available in the modelling literature (see Table  1). Specifically, we choose dimensional parameter values for which the model exhibits qualitative behaviour consistent with the three E's of Immunoediting (Dunn et al., 2004) (see Sections 5 and 7).
For ease of presentation, we henceforth omit hats in Equations (6)-(9). In Table 2, we list the default dimensionless parameters appearing in these equations. We note that there may be considerable variation in these values due to uncertainty in their dimensional values.

Qualitative behaviour
We solve Equations (6)-(9) numerically using the Dormand-Prince explicit adaptive timestepping method in Python (dopri5). This numerical method is based on a combination of fourth and fifth order Runge-Kutta schemes (Dormand & Prince, 1980) and can accommodate stiff equations, such as ours, where parameter values range over several orders of magnitude (see Table 2).
In Figure 2, we present simulation results for three choices of σ C , the rate at which cytotoxic T cells infiltrate the tumour, which reveal that our model reproduces the three E's of immunoediting (Dunn et al., 2004). For moderate values of σ C (see Figure 2 (a,b,c)), we observe bistability between a stable limit cycle which oscillates with small amplitude about an unstable steady state corresponding to a small tumour burden, and a stable steady state with a large tumour burden corresponding to tumour escape. Thus, depending on the choice of initial conditions, at long times the system will approach either the stable limit cycle or the stable steady state with a large tumour burden. For the periodic solution, the tumour, cytotoxic and helper T cell populations oscillate out of phase: growth of the tumour stimulates the proliferation of the helper T cells which, in turn, promote proliferation of the cytotoxic T cells, driving down the tumour population. When the tumour escapes, the cytotoxic and helper T cell populations evolve to constant values that are smaller than their baseline values of T H = σ H /δ H and T C = σ C attained when there are no tumour cells present. We remark that the tumour cell population has been nondimensionalised with its carrying capacity, κ, and therefore a large tumour population corresponds to a dimensionless tumour population that is close to one.
For a small increase of 12.5% in the value of σ C (from 2 to 2.25) (see Figure 2 (d,e,f)), we retain bistability between a tumour equilibrium and tumour escape. In this case however, we observe that the tumour equilibrium solution is no longer oscillitory but settles to a stable steady state with a small tumour burden. For a large increase of 60% in the value of σ C (from 2 to 3.2) (see Figure 2 (g,h,i)), we observe monostability of a tumour-free steady state solution which corresponds to tumour elimination: in this case, the immune response is strong enough to eliminate the tumour. We conclude that our model may exhibit multi-stability and that the nature and multiplicity of stable solutions change as σ C , the infiltration rate of the cytotoxic T cells, varies.

Steady state analysis
The results presented in Section 3 reveal that Equations (6)-(9) possesses several different types of stable steady state solutions and a stable limit cycle solution. As a first step to understanding the solution structure of our model, in this section we identify and characterise its steady state solutions.
Setting the time-derivatives to zero in Equations (6)-(8) yields three algebraic equations that define the steady state solutions (N, We note, from Equation (11), that for physically realistic solutions N * ∈ [0, 1). The solutions depend on five parameter groupings, First of all, Equations (11)-(13) admit a tumour-free solution,  Table 2). In each plot three initial conditions are considered and distinguished by linestyle: This solution is physically realistic provided 0 ≤σ H < 1. All other steady state solutions have 0 < N * < 1 and T * C =γ (1 − N * ). Furthermore, by eliminating T * C from Equations (12)-(13), it is straightforward to show that these steady state solutions lie at the intersection of the following two curves: In Figure 3, we use Equations (16)- (17) to show how the qualitative behaviour of the curves of T * H = H a (N * ) and T * H = H b (N * ) changes as the system parameters vary. For T * H = H a (N * ) three cases may arise: 1.ᾱ < 2Ñ: H a (N * ) > 0 for N * ∈ (0, 1) and attains a finite maximum in the feasible range; 2. 2Ñ ≤ᾱ < 1+Ñ 2 : H a (N * ) has two asymptotes N * ± in the range N * ∈ (0, 1); between the asymptotes, H a (N * ) < 0 and otherwise it is positive; 3.ᾱ ≥ 1 +Ñ 2 : H a (N * ) > 0 to the left of the asymptote at N * − ∈ (0, 1) and H a (N * ) is negative otherwise. (17) reveals that three cases may arise as the parameters are varied:

Inspection of Equation
Since the dependence of H b (N * ) on the model parameters is more involved than for H a (N * ), further details are left to Appendix 2. Notably H b (N * ) always becomes negative for N * sufficiently close to 1. The full range of behaviours can only be realised by varying bothk andσ C . It is not enough to fixk to a single value and varyσ C .
Recall that the physically realistic steady state solutions for which N * > 0 lie at those intersections of Equations (16) Guided by the results presented in Figure 3, in Figure 4 we show how the number and nature of these steady state solutions change as we varyᾱ andσ C for a fixed value ofk = 2 close to the default value (see Table 2).
For completeness, also indicated on these figures is the tumour-free steady state solution defined by Equation (15). In practice, up to four such steady states may arise (see Appendix 2 for details), although these cases are limited to small region of parameter space.
It is straightforward to show that if 2Ñ <ᾱ < 1+Ñ 2 (see Figure 4(b)) then Equation (16) has two asymptotes at N * = N * ± say. In addition, if 0 <σ C < 1 −σ H then H b (0) > H a (0) and there is a physically realistic steady state solution with 0 < N * < N * − . In this case, for small values ofσ C (i.e. low rates of infiltration of the cytotoxic T cells), there are three physically realistic steady state solutions with N * > 0, while for large values ofσ C the immune response is so strong that the tumour is eliminated and only the tumour-free steady state solution persists. Asσ C increases between these extremes, first the steady state with small N * is lost (at a transcritical bifurcation; see Section 5) and then the other two steady state solutions with 0 < N * collide at a fold bifurcation. (The caseᾱ < 2Ñ illustrated in Figure 4(a), where H a (N * ) attains a finite maximum, is more complicated and detailed in Appendix 2.) For the default parameters used in this study 2Ñ ≤ᾱ ≤ 1 +Ñ 2 (see Figure 4 (b)). In this case, for sufficiently smallσ C , three intersections occur: one at small N * and two at moderate N * . Asσ C increases, the two solutions at moderate N * are lost followed by, for largerσ C , the solution at small N * .

Linear stability analysis
In Section 3, we showed that, for different choices of the parameter values, our model may exhibit tumour elimination, equilibrium and/or escape. Tumour elimination and escape correspond to stable steady states, whereas equilibrium corresponds to either a stable steady state of intermediate size or a stable limit cycle. In this section, we decompose parameter space into distinct regions according to the number of steady states that exist and their linear stability, as well as identify regions of parameter space in which stable limit cycles exist.
We start by linearising about the steady state solutions (N * , T * H , T * C ), seeking solutions of the form  (11) and (13) change asᾱ andσ C vary, or, equivalently, the intersections of Equations (16) and (17) whereŇ,Ť H ,Ť C and λ are constants, and 0 < 1 is a small parameter. Substituting (18) into Equations (6)-(8) and equating to zero terms of O( ), we deduce that the eigenvalues λ = λ n (n = 1, 2, 3) satisfy |J − λI| = 0 where the Jacobian J (N * , T * H , T * C ) is given by Instability of a steady state solution is predicted when max n=1,2,3 (λ n ) > 0. Conversely, the steady state is said to be linearly stable if (λ n ) < 0 for n = 1, 2, 3. Note that in general the eigenvalues depend not only on the parameter groupings in Equation (14), but also on the unscaled parameters γ , p, k, α and δ H . For ease of presentation, it is convenient to use both scaled and unscaled parameters in subsequent analyses. It is straightforward to show that the eigenvalues associated with the tumour-free steady We deduce that where it exists (i.e. where 0 <σ H < 1) the tumour-free steady state is linearly stable if max 0, 1 −σ C <σ H < 1 .
We remark that the region of parameter space in which the small tumour steady state exists corresponds to the region in which the tumour-free steady state is linearly unstable. Characterisation of the other steady state solutions must be performed numerically. Figure 5 summarises how the multiplicity and linear stability of the steady state solutions of Equations (6)-(8) change as we varyσ C , the (scaled) basel rate at which cytotoxic T cells infiltrate the tumour (qualitatively similar behaviour is observed whenσ H is varied; results not presented). For the default parameter values (see Table 2), we identify five distinct regions separated by four bifurcation pointsσ F C ,σ T C ,σ Hopf C andσ Hom C . The qualitative behaviour in each region can be summarised as follows, •σ C >σ F C : the system is monostable; all initial conditions evolve to the tumour-free steady state (tumour elimination). •σ C =σ F C : fold bifurcation at which two non-trivial steady states are created. •σ T C <σ C <σ F C : the system is bistable; an unstable steady state with intermediate tumour burden separates the stable, tumour-free steady state from a stable solution with a large tumour burden (tumour escape); asσ C decreases, the tumour burden on the upper branch of solutions increases towards its (dimensionless) carrying capacity. •σ C =σ T C : transcritical bifurcation at which a new branch of stable steady state solutions emerges from the tumour-free branch.

•σ
Hopf C <σ C <σ T C : bistability is preserved; the tumour-free steady state exchanges stability with the new branch of solutions characterised by a very small tumour burden.

•σ C =σ
Hopf C : Hopf bifurcation at which the branch of solutions characterised by a small tumour burden loses stability to oscillatory solutions.

C <σ C <σ
Hopf C : bistability between a stable limit cycle and the steady state with a large tumour burden; all other steady state solutions are unstable. •σ C =σ Hom C : the limit cycle disappears when it collides with the intermediate tumour steady state (see Figure 7 and below for details).
In Figure 6, we present numerical results which confirm the predicted stability of the steady state solutions depicted in Figure 5 for small values ofσ C . All trajectories are attracted to the large tumour steady state (with N * ≈ 1). Some trajectories start within its basin of attraction while others first approach the tumour-free steady state solution along its stable manifold, before being repelled along its unstable manifold and attracted to the large tumour steady state. This behaviour is replicated for trajectories that start near the small and intermediate steady state solutions, both of which are unstable.
The bifurcation diagrams in Figure 5 indicate those values ofσ C for which oscillatory solutions occur. The limit cycle emerges whenσ C decreases through the critical valuẽ σ Hopf C and disappears asσ C →σ Hom C . Numerical results which confirm these dynamics are presented in Figure 7. Forσ C σ Hom C , we observe periodic solutions for which tumour growth stimulates infiltration by helper T cells which, in turn, recruit large numbers of cytotoxic T cells. As the tumour increases in size, immunosuppression of the helper T cells becomes important and their rate of proliferation falls while levels of cytotoxic T cells continue to rise, causing the tumour burden to fall. As levels of helper T cells fall, the proliferation rate of the cytotoxic T cells falls, causing their levels also to decrease. The immunosuppressed tumour then resumes its growth until it reaches a size at which an effective immune response is again stimulated, levels of helper T cells rise, and the cycle repeats. As the base rate at which cytotoxic T cells infiltrate the tumour approaches (from above)σ Hom C , the period of the limit cycles increase, with the tumour spending longer time periods during which the immune response is low.

Asymptotic behaviour forÑ 2 1
WhenÑ 2 1, immunosuppression occurs at small tumour sizes, and acts by reducing the proliferation rate of the helper T cells. We show in Section 6.1 how, in this case, it is possible to construct approximate expressions for the small tumour steady state solutions from which oscillatory solutions emerge and also to determine an expression for the critical value of σ C at which the Hopf bifurcation occurs. We also perform a weakly nonlinear analysis to determine the amplitude of the oscillations in the neighbourhood of the Hopf bifurcation.

The small tumour steady state
ConsiderÑ 2 ᾱ, as holds for the default parameters. Then the asymptote occurs at N * − =Ñ 2 /ᾱ = δ 1, to leading order in δ. We showed in the Section 4 that a tumour equilibrium steady state solution occurs with N * < N * − for certain parameter groupings. We can explicitly calculate the form of the equilibrium to leading order in δ by assuming  Table 2). Solution branches are distinguished by colour: the blue branch corresponds to the large tumour (escape) state; the red branch corresponds to the intermediate tumour state; the green branch corresponds to the small tumour state and the magenta branch corresponds to the tumour-free state. The local stability of the steady states is distinguished by linestyle: solid and dashed curves correspond to stable and unstable solutions respectively. The cyan curves correspond to the maximum and minimum values of the various cell populations for the periodic solutions. The bifurcations that occur are indicated on each plot: F=fold bifurcation; T=transcritical bifurcation; Hopf=Hopf bifurcation and Hom=homoclinic bifurcation.  Table 2). A series of typical trajectories starting in different regions of the phase plane are presented. While all trajectories eventually evolve to the large tumour steady state, in some cases they first approach close to the tumour-free steady state. Red and blue circles denote unstable and stable steady states respectively. The continuous and dashed lines represent different initial conditions. the expansions Inserting these expansions into Equations (16) and (17) gives and As anticipated from the analyses presented in sections 4 and 5, existence of the small tumour steady state occurs where 0 <σ C < 1 −σ H which coincides with instability of the tumour-free steady state solution.

Linear stability analysis
Following the general linear stability analysis in Section 5, we next perturb the above solution by an infinitesimal time-dependent function proportional to e λt (see Equation (18)).  Table 2). The trajectories start near the intermediate tumour steady state (see red branch in Figure 5). Panels (a) and (b) show how the cytotoxic T cells T C (t) evolve whenσ C σ Hom C and σ C <σ Hom We deduce that the eigenvalues λ solve |J − λI| = 0 where the Jacobian is given by (see (19)):

LETTERS IN BIOMATHEMATICS S53
Thus, the eigenvalues λ solve the following cubic Periodic solutions emerge when two roots of Equation (25) are purely imaginary, of the form λ 1,2 = ±iω. Hence, we locate the Hopf bifurcation point by seeking solutions to Equation (25) of the form where Eliminating λ 3 and ω 2 in (27) yields the following algebraic equation forσ C in terms of σ H and other model parameters at the Hopf bifurcation point: Viewing Equation (28) as a quadratic for σ H rather than a cubic forσ C , we can show where the small tumour steady state exists (i.e. where 0 <σ C < 1 −σ H ) there is a unique, positive root for σ H : where χ = (1 −σ C )(σ 2 C + γ (1 −σ C )). The eigenvalue λ 3 ≤ 0 over this range ofσ C , indicating stability of the oscillatory solutions where they exist.
The results of this analysis are verified in Figure 8(a), showingσ H versusσ C for both the asymptotic approximation given in Equation (29) and the numerical solution of the eigenvalue equation at the Hopf bifurcation point. There is close agreement for all values ofσ H . Figure 8(b) illustrates how the frequency at the point of emergence of the limit cycle varies withσ C . These results underpin the weakly nonlinear analysis presented in the following subsection.

Weakly nonlinear analysis
In this section, we perform a weakly nonlinear analysis in a neighbourhood of σ Hopf C , specifically σ C = σ Hopf C − 2 with 1, in order to determine the amplitude of the limit cycle and confirm its local stability. Since the analysis is involved, results are summarised in the main text (see Appendix 3 for details).
For ease of presentation, we let x = N, y = T C and z = T H . We introduce the long-time τ = 2 t and proceed by expanding x = (x, y, z) T as follows: where the long time scale τ = 2 t determines the limit cycle amplitude. Here x 0 is the steady state solution with N = O(Ñ 2 ) at σ C = σ Hopf C . To obtain the amplitude equation  Frequency ω associated with the eigenvalues λ = ±iω at the emergence of the limit cycle asσ C varies. Note ω → 0 asσ H → 0 and 1.
for the limit cycle we substitute (30) into Equations (6)-(8) and equate to zero terms of like orders of magnitude (full details are provided in Appendix 3). At O( ), we recover the linear dynamics ∂x 1 ∂t where J 0 = J (x 0 ). From the linear stability analysis, The most general solution of Equation (31) is Here, φ(τ ) is the amplitude of the first-order (linear) mode. The amplitude of the decaying mode, χ(τ ), plays no role in the subsequent analysis because it is non-secular. The equation for φ(τ ) is determined at O( 3 ) by eliminating secular terms as discussed below.
At O( 2 ) we find that x 2 is given by where x 20 = J 0 −1 (0, 1, 0) T and R a 2 through R d 2 are constant vectors (see Appendix 3). Notably, x 20 arises from taking σ C = σ Hopf C − 2 . At O( 3 ), the secularity condition supplies the following equation for the amplitude of the limit cycle dφ dτ where μ = μ r + iμ i and ν = ν r + iν i depend on the governing parameters.
Substituting φ(τ ) = r(τ ) exp (iϑ(τ )) into Equation (34) and equating real and imaginary parts yields dr dτ = (μ r + ν r r 2 )r and The equation for r is independent of ϑ and has solution If μ r > 0 > ν r , as is typical for the parameter values considered (see Table 2), then In such cases the limit cycle is a stable attractor.
We confirmed our weakly nonlinear analysis by computing the mean square amplitudes of ( Here ω nl is the nonlinear correction to the frequency, given by ω nl = μ i − ν i μ r /ν r . Samplying the solution 3600 times over the last 20T time units yields estimates of the amplitude A which agree to within 0.4% when = 0.1 and to within 0.016% when = 0.01.

The effect of immunosuppression
Having described the behaviours exhibited by our model for default parameter values, we now consider the effects of varying parameters associated with the immunosuppressive term in Equation (2). The steady state analysis in Section 4 showed that the values of the two immunosuppressive parameters,ᾱ andÑ, significantly affect the system equilibria.
We start by considering how changes inᾱ, the (scaled) proliferation rate of the helper T cells upon encounter with tumour antigen, affect the bifurcation structure in (σ C ,σ H ) space. We focus on varying the (scaled) rates of T cell infiltration,σ H andσ C , as experiments have indicated large patient-to-patient variability in these parameters (Oelkrug and Ramage, 2014;Lee et al., 1989;Kawai et al., 2008;Hiraoka et al., 2006). Furthermore, high levels of infiltrating T cells correlate with good prognoses (Lee et al., 1989;Kawai et al., 2008;Hiraoka et al., 2006).
In Figure 9, we present results which show how the bifurcations that occur in (σ C ,σ H ) parameter space change asᾱ varies. We observe seven distinct regions: 1. Black shaded region: monostability of the tumour-free steady state (i.e. tumour elimination). Here no other steady states exist. 2. Region shaded with vertical lines: bistability between the tumour-free steady state and large tumour steady state solutions (i.e. depending on the initial conditions either the tumour is eliminated or it escapes). 3. Dark grey shaded region: bistability between the small and large tumour steady state solutions (i.e. depending on the initial conditions either tumour escape or equilibrium will occur). 4. Light grey shaded region: existence and stability of limit cycle solutions. 5. White shaded region: multiple unstable steady state solutions and stability of the large tumour steady state (i.e. tumour escape). 6. Region shaded by horizontal lines: instability of the tumour-free steady state and stability of the large tumour steady state (i.e. tumour escape). 7. Region shaded by crossings: instability of the tumour-free steady state and stability of the small tumour steady state (i.e. tumour equilibrium).
Features thats are common for all values ofᾱ can be summarised as follows: • A stable tumour-free steady state exists for large values ofσ C andσ H (top right corner); • An unstable tumour-free steady state and a stable large tumour steady state exists for small values ofσ C andσ H (lower left corner); • The series of bifurcations that separate regions in which the tumour-free steady state is stable (large values ofσ C ) from regions in which the large tumour state is stable (small values ofσ C ) depend upon the value ofᾱ.
We now discuss the changes in the bifurcation structure that occur whenᾱ is varied about its default value (ᾱ = 0.19). For smaller values ofᾱ (e.g.ᾱ = 0.06, the small and intermediate steady state solutions coalesce at a fold bifurcation at a small value ofσ C (see Figure 9(a,d)). Increasingᾱ increases the size of the region in which the tumour-free state is the only stable attractor (black region) and reduces the size of the region in which the tumour-free and large tumour states are stable (region shaded with vertical lines), with bistability disappearing altogether ifᾱ is sufficiently large (see Figure 9 (c,f)). The size of the region characterised by bistability between the large and small tumour steady states disappears for large values ofᾱ (dark grey region). Comparison of Figure 9 (b,c) and (e,f) reveals that asᾱ increases the region in (σ C ,σ H ) parameter space in which only the large tumour steady state is stable decreases in size (white regions). Increasingᾱ increases the size of the region in which limit cycles exist (light grey region). We remark further that asᾱ increases the fold point at which the large and intermediate tumour steady states coalesce occurs at smaller values ofσ C . In particular, forᾱ > 1 +Ñ 2 , the intermediate and large tumour states disappear and only the small tumour and tumour-free steady state solutions exist (results not shown). In this case the tumour cells stimulate such a strong increase in the proliferation rate of the helper T cells that they increase to sufficient numbers that they can control or eliminate the tumour (see right hand panel of Figure 4(c)).
The results in Figure 9 provide a guide for potential treatment protocols for varying rates of helper T cell proliferation. In all cases, increasingσ C andσ H in equal amounts is most effective in driving the system away from tumour escape. However, the increase in either population must be limited to avoid an overzealous immune response. Targeting both populations in combination appears to be the most effective treatment. If only one population can be targeted alone, then targeting the cytotoxic T cell population appears to be most effective. The only exception occurs for an effective immune system (relatively high helper T cell proliferation, as in Figure 9(c)); in this case targeting either population is comparably effective.
In this section, we have examined how the bifurcation structure in (σ C ,σ H ) space changes as we varyᾱ, a parameter that reflects the extent to which the helper T cells increase their rate of cell proliferation on contact with tumour cells. We identified six distinct regions, distinguished by the existence and stability of steady states and the limit cycle. Variations inÑ yield similar results and are, therefore, omitted. Variation of other  Table 2). Key: black region (tumour-free steady state is the only stable attractor, i.e. tumour elimination); region shaded with vertical lines (bistability between the tumour-free and large tumour steady state solutions, i.e. depending on the choice of the initial conditions, the tumour will either escape or be eliminated); dark grey shaded region (bistability of large tumour steady state and small tumour steady state solution, i.e. depending on the choice of the initial conditions, either tumour escape or control will occur); light grey region (existence and stability of limit cycle solutions, i.e. tumour control); white region (multiple unstable steady state solutions and stability of large tumour steady state, i.e. tumour escape); region shaded with crossings (instability of tumour-free steady state solution and stability of the small tumour steady state i.e. tumour control); region shaded with horizontal lines (instability of tumour-free steady state solution and stability of the large tumour steady state, i.e. tumour escape). model parameters (e.g. γ , k and p) affects the location of the fold point relative to the transcritical bifurcation, and the Hopf bifurcation. Our results show that changing the parameters (ᾱ,Ñ) has a significant effect on the system's bifurcation structure in (σ C ,σ H ) parameter space. Taken together, our results indicate the extent to which combination therapies which not only boost the immune system (by increasingσ H andσ C ) but also block immunosuppressive effects may outperform monotherapies, with a single mode of action.

Discussion
We have proposed a new mathematical model of tumour-immune interactions in which helper and cytotoxic T cells interact with tumour cells. Our model captures the three Es of immunoediting -elimination, equilibrium and escape (Dunn et al., 2004). We have examined how the number and nature of the attractors (stable steady states and limit cycles) change as the infiltration rates of cytotoxic and helper T cells, σ C and σ H , are varied. Our focus on σ C and σ H is motivated by experimental observations which show high inter-patient variability in T cell infiltration rates and corresponding variability in observed outcomes (Oelkrug and Ramage, 2014;Fridman et al., 2012;Halama et al., 2011;Ali et al., 2016). Specifically, high levels of tumour infiltrating T cells have been correlated with long-term progression free survival in a variety of different cancer types (Rad et al., 2015;Preston et al., 2013;Hiraoka et al., 2006).
Our model exhibits three distinct behaviours as the infiltration rates of cytotoxic and helper T cells σ C and σ H vary. When σ C and σ H are both small the tumour always escapes. For intermediate values of σ C and σ H , two types of bistability can occur: either (a) bistability between tumour escape and elimination, or (b) bistability between tumour escape and tumour coexistence with the immune system, the co-existence state representing either a finite equilibrium state or a time-dependent periodic solution. When σ C and σ H are both large, the tumour is always eliminated. Based upon these results, we propose that patients may be categorised into three main groups: immunocompromised patients have low rates of T cell infiltration and their tumours will escape; healthy patients have large rates of T cell infiltration and any nascent tumours will be eliminated; the tumours of patients with intermediate levels of T cell infiltration will either escape or be controlled, the outcome depending on their initial size.
Stable limit cycles, characterised by periodic growth and suppression of the immune and tumour populations, have been observed in other mathematical models of tumourimmune dynamics (Kirschner and Panetta, 1998) and experimentally (Coventry et al., 2009). We confirmed that our model admits a stable limit cycle via a weakly nonlinear analysis which also yielded an equation for the amplitude of the limit cycle, valid near the Hopf bifurcation point.
To date, segregation of the T cell pool into helper and cytotoxic T cell populations, that can be targeted separately by immunotherapies, has received little attention. Which T cell population may elicit a stronger immune response is still widely debated (Hanson et al., 2000;Dudley et al., 2002). In this paper, we have developed a mathematical model which distinguishes their distinct roles as immune promoter and tumour killer respectively. We found that both T cell populations are important for tumour elimination. Given a sufficiently high rate of infiltration of either population, tumour elimination occurs. If, however, the rate of infiltration of helper T cells is low, then a large rate of infiltration of cytotoxic T cells is needed to eliminate the tumour. On the other hand, if the cytotoxic T cell infiltration rate is low, only a moderate increase in the helper T cell pool is required to eliminate the tumour. These findings suggest that targeting the helper T cell population may be more effective, and result in fewer immune-related adverse effects (Gangadhar and Vonderheide, 2014) than targeting the cytotoxic T cells directly.
Our model has also revealed that the immunosuppressive parameters α andÑ play a major role in determining response outcomes. Specifically, increasing the proliferation rate of helper T cells, α, increases the size of the region occupied by the stable tumour-free state. IncreasingÑ, the size of the tumour at which immunosuppressive effects come into play, decreases the size of the region in which tumour escape occurs. These results suggest that manipulating α andÑ may be effective therapies in cancer treatment. This could be achieved by blocking the PD-1/PD-L1 axis on T cells infiltrating the tumour (Hirano et al., 2005), a therapy which has already shown significant promise (Powles et al., 2014).
In future work, it would be interesting to explore the effect of different treatment protocols for the different patient groups identified in this study. Immunotherapies that could be explored include: (a) treatment with an immune promoter cytokine such as IL-2, which stimulates the proliferation of cytotoxic T cells (Rosenberg, 2014); (b) adoptive T cell therapy which targets the infiltration of cytotoxic and helper T cells (Maus et al., 2014;Restifo et al., 2012); (c) antibodies that block the PD-1/PD-L1 axis expressed by tumour and immune cells (Pardoll, 2012) affecting the size of the tumour at which immunosuppressive effects come into play; and (d) bi-specific antibodies that promote interactions between tumour and T cells (Schreiner et al., 2016). Comparing different treatment strategies, both individually and in combination, may reveal effective treatment protocols that could be used to treat particular patient groups.
A limitation of the current model is the assumption of spatial homogeneity. Future developments include extending the model to account for cell movement. This would allow us to study the spatial dynamics of interacting cytotoxic T cells, helper T cells and tumour cells. This can be achieved by extending our model as a system of partial differential equations where cell movement is assumed to be diffusive (Matzavinos et al., 2004;Mallet and De Pillis, 2006). Despite this limitation, our simplified model is able to capture the three E's of immunoediting and reveals distinct patient groups depending on levels of T cell infiltration into the tumour. Furthermore, examination of the parameters controlling the immunosuppressive affects included in the model reveals how these parameters and the rates of infiltration of the T cells might be manipulated to strengthen a weakened or ineffective immune system so that it can successfully eliminate any tumour.
To simplify yet retain key aspects of the complex model proposed by Robertson-Tessi et al. (Robertson-Tessi et al., 2012), we developed a system of three ODEs in the form of a predator-prey model like that of Kuznetsov et al. (Kuznetsov et al., 1994) but with the addition of a helper T cell population. Our model enabled us to study, through a combination of analytical and numerical methods, the effectiveness of targeting helper and cytotoxic T cells in an anti-tumour immune response. Our model does not distinguish the specific immunosuppression that impact the tumour at various stages of development as in Robertson-Tessi et al. (Robertson-Tessi et al., 2012), but was nonetheless able to illustrate how the anti-tumour immune response varies with the level of infiltrating helper and cytotoxic T cells. Our results suggest that a combined therapy targeting the infiltration of both T cell populations may result in the best outcome for patients having a broad range of hindered immune systems. Nonetheless there is a preference for targeting the cytotoxic T cell population for highly hindered immune systems.

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

Funding
This work was supported by the EPSRC [EP/G037280/1] and Hoffmann-La Roche.

S60
H. DRITSCHEL ET AL. Table A1. Glossary of terms used to describe tumour-immune interactions and immunotherapy.

Term Description
Antibody A molecule, produced by B cells, that binds a foreign body, labelling it for recognition by the immune system Antigen A molecule that provokes a specific immune response (T and B cells) against a foreign body Antigen presenting cell (APC) A cell of the innate immune response which presents antigen to naive T and B cells to provoke a specific immune response Effector T cells Short-lived, functional T cells (e.g. cytotoxic or cytokine-producing, helper T cells) Elimination The stage of immunoediting where the immune system recognises and eliminates the cancer cells before they grow to clinically detectable sizes Equilibrium The stage of immunoediting where the rate of tumour cell death caused by the immune response balances the tumour cell proliferation rate Escape The stage of immunoediting where tumour cell proliferation outstrips the suppressive effects of the immune system, enabling the tumour to grow to a substantial size Memory T cells Long-lived T cells that have already encountered their cognate antigen which are primed to respond rapidly if they re-encounter the pathogen Major Histocompatibility Complex (MHC) Proteins found on the surface of cells that present antigen for recognition and activation of the cognate immune cells Granzyme B A protein produced by cytotoxic T cells and natural killer cells which perforates the target cell Humoral response Target cells are eliminated indirectly via production of antibodies by B cells Interferon-γ (IFN-γ ) A cytokine which is critical for both innate and adaptive immunity; it activates macrophages of the innate immune response induces expression of MHC class I complexes on cell surfaces Interleukins (ILs) Proteins produced by cells of the immune system that regulate (e.g. promote or suppress) the response to the foreign body Transforming Growth Factor β (TGF-β) A cytokine which regulates cell growth, proliferation, differentiation and death Programmed Cell Death Protein 1 (PD-1) A transmembrane protein primarily expressed by T cells which plays a role in suppressing the immune response Programmed Cell Death Ligand 1 (PD-L1) A transmembrane protein primarily expressed by tumour cells which plays a role in suppressing the immune response

Appendix 2. Conditions for tumour equilibria
In Section 4, we introduced the functions T H = H a (N ) and T H = H b (N) (see Equations (16)- (17)) and showed that their positive intersections give rise to tumour equilibria. Figure 4 illustrates the ways in which these curves may intersect. In this appendix, we detail how the qualitative behaviour of the functions T H = H a (N) and T H = H b (N) vary with model parameters and identify a region of parameter space in which four (positive) equilibria may occur. For the curve T H = H a (N), three cases may arise, depending on the values ofᾱ andÑ: (1)ᾱ < 2Ñ: T H = H a (N) > 0 for N ∈ (0, 1) and attains a finite (positive) maximum: at N =Ñ.
(2)σ C /k ≤ 1: T H = H b (N) attains a maximum H b,max = 1 +k − 2 kσ C at N = 1 − σ C /k ∈ (0, 1). For physically realistic equilibria we require H b,max > 0 or, equivalently, These results are summarised in Figure ?? where we indicate the region of (k,σ C ) parameter space in which physically realistic tumour equilibria may occur.
where f (t, τ ) contains time dependencies other than e iωt and is therefore non-secular. The components of the vectors B l and B c are given by and The particular solution to Equation (C13) is x 3 (t, τ ) = X 3 (t, τ ) + u(τ )e iωt + c.c.
where X 3 consists of non-secular terms coming from f (t, τ ) + c.c., while u accounts for all nonsecular terms in B l and B c . This means u must be a linear combination of v 2 and v 3 : where r j = (iωI − J 0 )v j . Using this in Equation (C13) above, upon equating all potentially secular e iωt terms we obtain dφ dτ v 1 + C 2 r 2 + C 3 r 3 = (B l + B c |φ| 2 )φ.
This can be regarded as a system of equations to determine the vector (dφ/dτ , C 2 , C 3 ). Define Therefore, the amplitude equation satisfied by the limit cycle is where μ = s lx and ν = s cx .