Activation of the immune response by cytokines and its effect on tumour cells: a mathematical model

Abstract In this paper, we present a mathematical model at the cellular level of the tumour–immune competition mediated by the cytokines. The model consists of a system of nonlinear differential equations describing the intracellular interactions between the tumour and the immune cells in the presence of the cytokines. A detailed phenomenological description of the model based on the kinetic theory for active particle approach is carried out to formulate the model. Well-posedness is presented to establish local and global existence. Numerical simulations are addressed to show how initial conditions and model parameters influence the output of the model. Under a suitable choice of the model’s key parameters and the cytokines’ initial activation levels, the simulation results show that the activated immune system is able to achieve a total elimination of the cancer cells.


Introduction
The significance of the immune system in tumour development, progression and destruction, as well as the complexity of the tumour-immune interactions have been extensively studied in recent decades (see reviews and references therein (Zamarron & Chen, 2011;Kim, Emi, & Tanabe, 2007;Cal, Molon, & Viola, 2017;Finn, 2012;Wargo, Reddy, Reuben, & Sharma, 2016). The evolution of tumour cells in a host is governed by more than just the activities of cancer cells alone. It is largely due to a plethora of inter-and intracellular interactions within the tumour microenvironment (Whiteside, 2010;Wargo et al., 2016). The immune cells play a key role in the tumour microenvironment, which contains heterogeneous cell populations such as cancer stem cells, stromal cells, fibroblasts and endothelial cells (Hanahan, Lisa, & Coussens, 2012).
Complex cellular interactions and crosstalks take place within the tumour microenvironment and are governed by an intricate network of biological pathways (Bor-Ching et al., 2018). The various cell types in the tumour microenvironment upon activation, secrete multiple signals also known as cytokines (Chen & Mellman, 2017). These signals control the strength and timing of the anticancer response and play an important role in tumour pathogenesis (Burkholder, 2014). The cytokines are low molecular weight proteins that bind on a responsive target cell and act as intracellular messenger invoking particular biological activity (Landskron, De la Fuente, Thuwajit, Thuwajit, & Hermoso, 2014). The cytokines produced by the cells in the tumour microenvironmentare are key factors in modulating immune response either against or in favour of tumorigenesis (Showalter et al., 2017 September). Depending on the cytokines and other signals in the tumour microenvironment, recruited immune cells will either form a pro-tumour immunity or an anti-tumour immunity (Burkholder et al., 2014). Cytokines are also considered for potential therapeutic and preventive targets as well as prognostic factors. They have been recently explored in cancer immunotherapy (i.e. the use of the immune system to treat cancer) (Hong-Mei, 2014;Sylvia & Margolin, 2011).
Understanding such a complex and continuously changing network of multiple signals that control cancer-immune interactions not only does it require ongoing and sophisticated experimental and clinical research, but it also calls for mathematical modelling. Several mathematical models of tumour-immune system interactions have been developed over the past decades. They have included a broad range of mathematical methods such as differential equations, spatial and nonspatial multiscale models and agent-based models to name just a few. These models have examined different characteristics of the cancerimmune complex interactions both in general and in other specific cases such as immune surveillance, suppression and escape, and during therapies. For a comprehensive survey of these models, the reader is referred to excellent review articles in (Altrock, Liu, & Michor, 2015;Adam & Bellomo, 2012;De Pillis, Eladdadi, & Radunskaya, 2014;Eladdadi, Kim, & Mallet, 2014).
Our mathematical model describes the interactions between single cells of the tumourimmune system and is developed according to the kinetic theory for active particle approach. This mathematical framework applied to cancer modelling was pioneered by (Bellomo, Preziosi, & Forni, 1996) and developed further in numerous publications as documented in the review paper (Bellomo, Elaiw, Althiabi, & Alghamdi, 2015). The main advantage of the kinetic theory is that it provides a deeper insight into the interactions between the tumour and immune cells at a cellular level at an early stage of tumour development, that is before the tumour becomes a macroscopically observable spatial structure. This early stage is significant since the competition between the tumour and immune cells can still lead to a complete elimination of cancer cells by the activated immune system (Kim et al., 2007). In fact, the tumour immuno-surveillance hypothesis, first formulated by Burnet back in 1954, indicates that the immune system is capable of inhibiting the growth of very small tumours and eliminating them before they become clinically recognized (Burnet, 1954).
The tumour-immune competition involves several interacting cell populations within the tumour microenvironment each characterized by a unique microscopic internal biological function. The immune system consists of various cell subpopulations each with different biological functions. They are recruited into the tumour microenvironment as a result of the immune system activation by different cytokines (Nagarsheth et al., 2017). Similarly, there is a large network of active cytokines with different biological activities and attributes in the tumour microenvironment (Bor-Ching et al., 2018). To simplify the complexity generated by a large number of cell subpopulations and to concentrate on capturing the overall evolution of the immune system activation by the cytokines, we consider the immune cells as one whole population. Similarly for the cytokines, we consider all type of cytokines as one population.
In this paper, we build upon our previous model in (Bellouquid & CH-Chaoui, 2014), where we presented a nonlinear integro-differential model of the tumour-immune competition mediated by the cytokines activity. The model describes the conservative interactions between the tumour and immune cells, that is, the competition where no cell proliferation/destruction occurs. In this work, we extend this model to include the activation density of the cytokines by considering both the conservative and proliferative/destructive interactions of the tumour and immune cells. Particularly, we investigate how and under which conditions the immune activation by cytokines leads to a complete elimination of cancer cells in their early stage of development.
This paper is organized as follows. The formulation of the mathematical model is detailed in Section 2. The well-posedness of the initial value problem is presented in Section 3. Numerical study is given in Section 4. The conclusion and discussion are presented in Section 5.

Model formulation
In this section, we present the mathematical model describing the interactions between the tumour cells and the cytokines-activated immune cells at the cellular level. We first start by briefly reviewing the kinetic approach in modelling a system of a large number of interacting cells (Sec 2.1), followed by the main assumptions (in Section 2.2) needed for the formulation of the mathematical model which is detailed in (Section 2.3).

The kinetic theory for active particles approach
The kinetic theory for active particles approach refers to the classical models of kinetic theory defined by the Boltzmann equations which are also known as the generalized kinetic Boltzmann models. They are the fundamental models of nonequilibrium statistical mechanics that provide the conceptual framework for generalizing the methods of kinetic theory to various fields of applied sciences (Cercignani, Illner, & Pulvirenti, 1994). The kinetic theory is praised for its ability to help understand phenomena of nonequilibrium statistical mechanics which are not described by the traditional macroscopic approach (Bellomo & Pulvirenti, 2000). The models based on the kinetic theory describe the evolution of the one-particle distribution function over the physical state of a large systems of interacting particles/cells. Each particle/cell is characterized by a certain microscopic state that describes its functional state such as the position, velocity and the biological activity. The overall system can be divided into finite subsystems, where each subsystem is a collection of particles/cells showing common activities (or biological expressions). These subsystems usually consist of interacting populations with the possibility of shifting from one population to the other. The state of each subsystem is described by a distribution function over the microscopic state of the interacting cells, whereas the evolution of the system is defined by local cells interactions. For a comprehensive derivation of the kinetic theory for active particles mathematical frameworks, the interested reader is referred to the books in (Bellouquid & Delitala, 2006;Bellomo, 2008;Bellomo & De Angelis, 2008). Readers interested in the kinetic approach applied to cancer modelling is referred to these excellent papers and the references therein (Bellomo, Bellouquid, & Delitala, 2004;Bellomo & Delitala, 2008;Bellouquid, De Angelis, & Knopof, 2013;De Angelis, 2014;Bellouquid & De Angelis, 2011;De Angelis & Lods, 2008;De Angelis & Jabin, 2003;De Angelis & Jabin, 2005).

Model description
In this section, we present the assumptions for our model based on the discussion above and the scientific literature on the interaction and competition at the cellular level between the tumour and immune cells (Zamarron & Chen, 2011;Kim et al., 2007;Cal et al., 2017;Finn, 2012;Wargo et al., 2016). The early stage of cancer development includes activation and inhibition of the immune system through cytokine signals which regulate cellular activities. We start by detailing the first four steps of the kinetic theory that will help in formulating our mathematical model.

Cell populations
The biological system considered in this model consists of three interacting populations in the tumour microenvironment: cancer cells, immune cells and cytokine signals labelled, respectively, by the indices 1, 2 and 3.

The microscopic state
The tumour and immune cells as well as the cytokine signals are viewed as active particles with a microscopic state describing a unique internal biological function or activity. From a mathematical viewpoint, this biological activity is represented by a real variable u ∈ D u ⊂ R. The value of u is a measure of the cell's ability to prevail in the cell-cell or binary interactions. For the cancer cells, u represents the progression from the normal to cancerous state, whereas for the immune cells and the cytokines, u describes activation. Cells with a positive activation state value are called active, whereas cells with a negative value are called passive (or inhibited). Specifically, in the case of cancer cells, negative values of u correspond to normal endothelial cells, while positive values of u correspond to pathological state or cancer cells. For the immune cells and the cytokines, positive values of u means activation, while negative values means inhibition.
Based on the generalized kinetic Boltzmann equations (Bellouquid and Delitala 2006), the microscopic state of the system of each population is described by the normalized density functions f i (where i = 1, 2, 3) over the activity u as follows: where the densities N i (t, u) are such that dn i = N i (t, u)du denotes the number of cells per unit volume whose state is, at time t, in the interval [u, u + du], and that ν 0 is the total number of cells at t = 0.

The macroscopic quantities
The macroscopic quantities portray the behaviour of the system and are usually computed using the normalized density functions f i (t, u). The two main macroscopic quantities that we focus on in our mathematical model are the size (or number density) of the cell populations denoted by n(t) and the activity of each population denoted by A(t), which are respectively defined as follows: • Size or number density of each of the populations at time t is recovered as zerothorder momenta of the distribution function f i = f i (t, u) (i = 1, 2, 3) and is given by: while the total number of cells at time t is given by , which may depend on time due to the birth (proliferation) and death (destruction) processes related to the nature of the cell interactions, as well as to the flux of cells through the boundaries of the volume.
• The activation of each cell population at time t is given by the first-order momenta: The sizes and activation of each cell type and their meanings are summarized in Tables 1 and 2, respectively: where the total endothelial cell density is given by,

Microscopic interactions
Below, we detail the binary cellular interactions in our proposed model both biologically and mathematically. We first start by describing the conservative interactions followed by the proliferating/destructive interactions.
• Conservative interactions: they modify the progression of the endothelial cells, the activation of the immune cells and the cytokines. Mass conservative encounter rates between two population, say i and j, with microscopic states u * and u * , respectively, are quantitatively represented by the transition probability density m ij (u * , u * ). It is assumed to be a Gaussian distribution function with the most probable output defined by the mean value m ij (u * , u * ) which are described by the δ function over m ij . The conservative interactions (C.1 -C.7) are listed below.  Table 2. Activation and their meaning. The inner tendency of endothelial cells to degenerate and progress towards a pathological state c 12

Activations Description
The ability of the active immune cells to reduce the progression of progressing cells c 21 The ability of tumour cells to inhibit immune cells c 23 The activation of immune cells by cytokine signals c 32 The progressive decay of the cytokine signals activity p 11 The proliferation rate for progressing cells due to the interaction with normal endothelial cells p 12 The destructive ability of the immune system p 21 The proliferation of immune cells due to interaction with progressing cells (C.1) Interactions between the endothelial cells: endothelial cells show an inner tendency to degenerate and progress towards a pathological state with most probable output given by: u * ∈ R : m 11 = u * (1 + c 11 ), where c 11 is the ability of endothelial cells to degenerate and progress towards a pathological state. (C.2) Interactions between the tumour and immune cells: Tumour cells means that u * > 0, and in this case there are two scenarios depending on the activity of the immune cell: • if the immune cell is not active, i.e. u * ≤ 0, then the state of the tumour cell does not change and: u * ≤ 0, u * ∈ R, u * ≥ 0, u * ≤ 0 : m 12 = u * • if the immune cell is active, i.e. u * ≥ 0, then the most probable output is given as follows: where c 12 refers to the ability of the active immune cells to reduce the evolution of the progressing cells (C.3) Interactions between the immune and endothelial cells: • if the endothelial cells is not progressing (normal), i.e. u * ≤ 0, the state of an immune cell (u * ∈ R) does not change.
• if the endothelial cell is progressing (i.e. u * > 0), then the state of the immune cell does not change after interaction with normal the endothelial cell: u * ∈ R, u * ≤ 0, u * < 0, u * > 0 : m 21 = u * • if the endothelial cell is progressing (i.e. u * ≥ 0), after interaction with an active immune cell, the progressing cells will inhibit the activation of the immune cell, and the mean value is given as follows: where c 21 is a parameter which indicates the ability of tumour cells to inhibit immune cells:

4) Interactions between the immune cells: interactions between immune cells do
not change their respective states: u * , u * ∈ R : m 22 = u * (C.5) Interactions between the immune cell and cytokines: Interaction between an immune cell (either activated or inhibited) with molecules of cytokines increases the activation of the immune cell. The mean value is then given by: u * ∈ R, u * > 0, m 23 = u * + c 23 , where c 23 indicates the activation of immune cells by cytokine signals. (C.6) Interactions between the cytokines and the immune cell: Interaction between a molecule of cytokines and an immune cell decreases the activation of cytokine signals, and the mean value is given by: where c 32 depicts the progressive decay of the cytokine signals activity (C.7) Interactions between the cytokine and endothelial cells: interaction between endothelial cells and cytokine signals do not lead to any modification of their states. u * , u * ∈ R : m 13 = m 31 = m 33 = u * • Modelling microscopic proliferating/destructive interactions: the nonconservative interactions involving the three populations are modelled via the proliferating/ destructive rate μ ij (u * , u * ), referring to every possible nonconservative interactions between the cell pairs and are listed below ((P.1)-(P.6)): (P.1) Cancer cells undergo uncontrolled mitosis stimulated by encounters with normal cells due to their angiogenic ability: The term μ 11 (u * , u * ) models the net proliferation into progressing endothelial cells due to interaction with normal endothelial cells • p 11 is a parameter which characterizes the proliferating ability of tumour cells. (P.2) The proliferation rate of normal endothelial cells, u * ≤ 0, due to encounters with other endothelial cells, is equal to zero. Encounters between tumour cells do not lead to any proliferation or destruction. (P.3) Tumour cells are partially destroyed due to encounters with active immune cells (u * ∈ [0, ∞)), where p 12 is a parameter that characterizes the destructive ability of active immune cells. (P.4) The proliferation rate of nonprogressing cells u * ≤ 0 due to encounters with immune cells is equal to zero.
(P.5) Active immune cells (u * ∈ [0, ∞)) proliferate due to encounters with progressing cells, although some of them may be inhibited: where p 21 is a parameter that describes the proliferation ability of active immune cells. (P.6) The proliferation rate of inhibited immune cells (u * ≤ 0) due to encounters with cells of the first population is equal to zero.

Mathematical model
Based on the phenomenological assumptions ((C.1)-(C.7)) and ((P.1)-(P.6)) listed above, we now proceed to formulate the mathematical model by putting all of them together. This model which describes the tumour-immune competition mediated by the cytokines at the cellular level consists of a system of three coupled nonlinear integro-differential equations. The dynamics of the system refers to an early stage of the competition with both conservative and nonconservative cellular interactions. The evolution equations are derived by suitable mass balance equations between the rate of change in the number of cells per unit volume and the gain (positive flow) and loss (negative flow) of cells entering or leaving this volume due to binary interactions as described in (Arlotti, De Angelis, Fermo, Lachowicz, & Bellomo, 2012;Bellouquid & Delitala, 2006) and are given by the following mathematical framework: Inlet and outlet flux due to conservative interactions for i = 1, 2, 3, where C i and P i are suitable smooth operators acting over the distribution functions [f ](t, u), which lead to the calculation of the net flow of active particles in the elementary volume of the micro-states as follows: • C i [f ](t, u) is the net flow at time t due to conservative interactions.
• P i [f ](t, u) is the net flow at time t due binary interactions that are either proliferative or destructive.
Technical calculations (Bellouquid & Delitala, 2006) of the mathematical structure defined in Equation (2) lead to the following evolution system of equations: Self progression toward pathological state Reduction of the progression of progressing cells by active immune cells Destruction of progressing cells by active immune cells Inhibition of immune cells by progressing endothelial cells

Activation of immune cells by cytokine signals
Utilization and decay of cytokine signals The corresponding macroscopic model related to the number densities n T is obtained by suitable integrations of Equation (3) over the biological activity variable u ∈ (0, +∞): The macroscopic model related to the activities A T 1 , A A 2 and A A 3 corresponding respectively to the tumour cells, activated immune cells and cytokines is obtained using Equation (3): System of Equation (5) is not in a closed form because f 1 (t), f 2 (t) and f 3 (t) are unknown functions. This macroscopic model can, however, describe the evolution of the activations of cells expressing common biological expressions. Specifically, it describes the effects of the cytokine signals on the interaction between the tumour cells and activated immune cells. The asymptotic analysis of this model is not possible because of its dependence on both f i (t) and zeroth-order moment. In this case, we restrict the qualitative analysis of this model to numerical simulations which give an insight into the behaviour of the system with different activation values.
The model is characterized by eight phenomenological parameters describing the transitions and proliferation/destructive rates of the cellular interactions. Transitions parameters in conservative encounters are denoted by c ij , while proliferation/destructive encounters are denoted by p ij . The first index denotes the population that undergoes a change, while the second index indicates the population that causes the change. The model parameters and their meaning are summarized in Table 3.

Well-posedness of the initial value problem
In this section, we present the well-posedness of the Initial value problem for the mathematical model Equation (3). The problem is stated by linking the evolution equations to suitable initial conditions as follows: where h = (f 1 , f 2 , f 3 ) and f i0 = (f 10 , f 20 , f 30 ).
The subscripts i = 1, 2, 3 correspond to tumour cells, immune cells and cytokines, respectively, while H[h] represents the right-hand side terms of the system of Equation (3).
Let L 1 (1 +|u|) m du denote the (1 +|u|) m −weighted Lebesgue L 1 space, with the norm is denoted by . 1 , with m ≥ 0, and let X = L 1 ((1 + |u|) m du) 3 is the Banach space equipped with the norm Let X + be a positive cone: and Y = C([0, T], X + ) the space of the functions continuous on [0, T] with values in a Banach space X + , endowed with the uniform norm The well-posedness of (6) is given by the following two theorems: Theorem 3.1: Well-posedness of the problem and nonnegativity of the solution.
For any initial data f 0 ∈ X + , there exist two positive constants T and a 0 such that Problem (6) has a unique local in time solution f ∈ C([0, T], X + ), which satisfies the following estimate: Proof: The problem (6) can be written as follows: where the operator J , defined by Equation 3, is written as follow: Briefly, the proof of the local existence and uniqueness is based on standard fixed point arguments. Let f and g in X, then J [f ] ∈ X, and (10) Exploiting (9) and (10), one has and This implies, there exist two constant T and a 0 determinate only by C, and f 0 , such that which implies that is a contraction on a ball in Y of radius a 0 f 0 if T < 1 4C f 0 .
Thus, there exists a unique local solution f (t) of (6) on [0, T]. Finally, in order to prove the nonnegativity of the solution, let us write the componentwise definition of Problem (6) in the following equivalent form:

S188
M. CH-CHAOUI ET AL. where and ( 1 4 ) which after integrating in t one obtain This conclude the proof, as long as f exist, due to nonnegativity of f i0 , and nonnegativity The existence of solution f of Problem (6) can be extended over the whole real positive axis R + , by the following: where C T constant depending on T and on the initial data.
Proof: Global existence can be proved for Equation (6) as follows. Bearing in mind the results of Theorem 1, it remains to find a priori estimates for the solution. Integrating the first equation of (3) with respect to u in ( − ∞, 0) yields then by integrating the first equation with respect to u in R one obtain

Hence the total number of progressing cells is bounded on each finite interval [0, T].
Integrating the second equation of (3) with respect to u yields It follows, from Equations (17) and (18), that which gives that n 2 (t) is bounded on each finite time interval [0, T]. Now integrating the third equation of (3) with respect to u in R, one gets easily which finally gives (17) with C T given by:

Numerical study
Simulations are obtained using the so-called generalized collocation method (Bellomo, 1997). The activity variable u is discretized by suitable set of collocation points, then the distributions functions f 1 (t, u), f 2 (t, u)and f 3 (t, u) are interpolated. Sinc functions is used for interpolations in assuming the natural trend to zero at infinity (Bellomo & Ridolfi, 1995). The integral terms are approximated by means of weighted sums in each collocation points. The integro-differential initial value problem is transformed into an initial value problem for ordinary differential equations, describing the evolution of f i (t, u), i = 1, 2, 3 in the node of discretization (Bellouquid & Delitala, 2006). We are interested in the time evolution of the activations A T 1 and A A 2 corresponding to the aggressiveness of the cancer cells and the efficiency of the activated immune cells with and without the effect of the cytokines, respectively. We are also interested in the time evolution of f i (t, u), for i = 1, 2, 3 represented in the evolution equations Equation (3). The aim of the numerical simulation is to show some dynamics of the tumour-immune cell competition in the presence of the cytokines. We are also interested in determining the conditions under which the activated immune system can win the competition by achieving a total regression of the tumour's activity. In particular, we focus our numerical study on the two main aspects of the immune response against the onset of cancer cells: • The ability of the tumour cells to inhibit the activated immune cells.
• The role of the cytokines' activation of the immune cells and the progressive decay of the cytokine signals activity.
Our model involves eight phenomenological parameters (Table 3), each one describes a well-defined biological function. Of interest to our study are the following model parameters: c 12 (the ability of the active immune cells to reduce the progression of the tumour cells), c 21 (the ability of tumour cells to inhibit immune cells), and c 23 (the activation of immune cells by cytokine signals), c 23 (the progressive decay of the cytokine signals activity). To investigate the effects of the cytokines activation on the immune response to cancer cell progression, we first start by varying the values of c 12 , c 21 , c 23 and c 32 while keeping the rest of the parameters' values (in Table 3) fixed.
Particularly, we focus our simulations on a biological situation where endothelial cells are characterized by a nonnegligible rate of progression toward a pathological state (c 11 = 0.3), a low proliferation rate of cancer cells (p 11 = 0.1) and an intermediate level of the immune response (p 12 = 0.3). We use a low proliferation rate of the immune cells (p 21 = 0.1) in the simulation. These parameters values were used in previous studies as well (Afraites, Atlas, Bellouquid, CH-Chaoui, 2012;Bellouquid & Delitala, 2006;Bellomo & Delitala, 2008;Bellouquid, CH-Chaoui, & De Angelis, 2015). All simulations are performed with the following nonzero initial conditions, unless otherwise noted: The results of these numerical simulations are shown in both 2D and 3D figures. The vertical axis show the population under consideration over the state variable and time. The population size is rescaled by dividing the population densities N i (t, u), i = 1, 2, 3 by the constant size of the total cell population, as explained in Equation (1). We describe the numerical results in the following two sections.

The ability of the tumour cells to inhibit the activated immune cells
The first objective of the numerical study is to examine the ability of the tumour cells to inhibit the activated immune system, in the absence of the immune activator. For that, we focus on two-model parameters c 21 and c 12 , which describe the ability of the tumour cells to inhibit immune cells and the ability of the active immune cells to reduce the progression of the tumour cells, respectively. To this end, we vary the model parameters c 21 with respect to c 12 , in order to investigate their effects on the dynamics of the competition between the tumour and immune cells. In this first scenario, we consider the onset of the tumour cells and their competition with the immune system for increasing values of c 21 . Specifically, we set c 12 = 0.1 which represents a weaker active immune system and vary the numerical values of c 21 . It is expected that increasing c 21 will result in a increase in the ability of the immune system to control cancer cells progression.
• For a low value of c 21 = 0.1, the activation of the immune cells decreases more rapidly than the activation of the tumour cells as depicted in Figure 1(a). This simulation shows that the immune cells are not able to completely reduce the aggressiveness of the tumour cells as a consequence of a weaker immune system. Figure 2  • For a high value of c 21 = 0.7, we observe that the immune cells activation is able to reduce the activation of the tumour cells as represented by the blue line in Figure 1, and consequently leads to a total reduction of the tumour cells' density as depicted in Figure 2(d).
In summary, Figures 1 and 2 show a comparison between the activation's trends of the competing tumour-immune cells. Specifically, Figures 1(b) and 2(d) highlight the crucial role of an efficient immune system that led to a total regression of the tumour cells.
In the second scenario, we set c 12 = 0.6 (higher immune defence ability) and vary c 21 (the ability of the tumour cells to inhibit immune cells) accordingly.
• For a low value of c 21 = 0.1, simulation results in Figure 3 show a higher defence level of the immune system. Specifically, panels Figure 3(a) & (c) shows a complete depletion of the tumour cells while panel Figure 3(b) illustrates an active immune state. This higher state of activation is observed for an early stage of the competition (Figure 3(c)). In addition, we observe that conservative interactions shift the immune cells to a higher value of activation (Figure 3(b)). These simulation results show that the tumour cells are asymptotically destroyed and completely go to extinction, suggesting that the tumour cells are not able to inhibit the immune cells under a higher immune defence ability (represented by c 12 = 0.6).
• For a higher value of c 21 = 0.7, the simulation results show a decrease in the immune cells' activation. Indeed, as shown in Figure 4(b) the immune system is able to reduce the activation of tumour cells. On the other hand, this behaviour induce a progressive weakness in the immune system and the final outcome depends on which of the two competing populations is the strongest in inhibiting the other.

The role of the cytokines' activation of the immune cells and the progressive decay of the cytokine signals activity
The second objective of the numerical study is to investigate how the choice of the initial activation A A 30 affect the dynamics of tumour cells contrasted by activated immune cells, especially how the two model parameters c 23 and c 32 modify the output of the competition with respect to the initial values of cytokines activation.
• Case of higher initial cytokines activation A A 30 = 0.071 As a first step, considering the case expressed in Figure 4(a). Introducing the cytokines with the value A A 30 = 0.071. This value is higher than both the tumour cells' activation, A T 10 = 0.070, and the immune cells' activation A A 20 = 0.043. We set c 23 = 0.4 and c 32 = 0.5. Figure 4(c)shows a decrease in the tumour cells activation, and a slow decrease in the immune cells activation. Figure 5(b) shows the distribution of the immune cells at different instant of time. Nevertheless, increasing the value of c 32 = 0.7 induces immune cells towards a higher level of defence against progression of tumour cells, while these one decay asymptotically to zero in time. This is due to immune action reinforced by the cytokines. • Case of small initial cytokines activation A A 30 = 0.007 Considering the same case where c 12 = 0.6 & c 21 = 0.7, as shown in Figure 4(a). We choose a low value of A A 30 = 0.007, while A T 10 and A A 20 kept the same initial values (A T 10 = 0.070 and A A 20 = 0.043). In Figure 6(a), the immune system contrast more efficiently the tumour activation. In fact, the activation of the tumour cells is decreased quickly and kept under control, while the immune system express rapid manifestation toward killing the tumour cells. Figure 5(e) highlights this behaviour of immune cells to overcome and eliminate the pathology. Nevertheless, introducing the cytokines with a low initial activation value induces a rapid decay in the tumour cells as reported in Figure 6 . If we stretch out the time simulation, a rapid initial growth of immune cells is identified up to a maximum value which correspond to an initial phase characterized by a rapid activation of immune cells up to maximum value (corresponding to time 10) that slowly relax and return to it sentinel state (it initial activations). However, in consistency with our previous work (Bellouquid & CH-Chaoui, 2014), considering the ration α = c 23 c 32 between these two parameters, corresponding to the ability of immune cells to prevent the aggressiveness of tumour cells, such that that increasing α from 0.1 to 0.9 (Figure 6(c)) induces a strong immune defence, and the tumour cells are controlled and can no longer return to grow. Through simulations, we can conclude that immune cells contrasts the tumour activation more efficiently when the cytokines is applied with low initial activation, while the immune system displays a faster defence to kill tumour cells.

Conclusion
In this paper, we derived a phenomenological model based on the mathematical kinetic theory describing the response of the cytokines-activated immune system to the onset of tumour cells. The model consists of a system of integro-differential equations representing the time evolution of the distribution functions over the microscopic biological state of the interacting tumour, immune cells and cytokine signals. Comparing to the analysis developed in (Bellouquid & CH-Chaoui, 2014) which dealt only with the conservative intracellular interactions, this model includes both the conservative and proliferative/destructive interactions between the tumour and immune cells in the presence of the cytokines. The addition of the proliferative/destructive interactions resulted in more tumour-immune dynamics that are biologically relevant such as the immune defence need large time to eradicate the pathology. The numerical study aimed at visualizing specific features of the tumour-immune competition dynamics. The simulation study concentrated on two main aspects of the immune response against the evolution of cancer cells: (1) the ability of the tumour cells to inhibit the activated immune cells, and (2) the role of the cytokines' activation of the immune cells, and the progressive decay of the cytokine signals. Numerical study also helped in defining the conditions under which the activated immune system was able to win the competition by achieving a total regression of the tumour's activity as shown in Figures 5 and 6. The numerical study showed that the tumour cells are completely depleted when the ability of the active immune cells to reduce the progression of the tumour cells was (six times) higher than the ability of tumour cells to inhibit immune cells. Other cases where the tumour was completely eliminated is when the cytokines is introduced with a low initial activation, an initial growth is observed up to a maximum value corresponding to a faster activation of immune cells. In the main time, the tumour is suppressed and the immune activation still being activated.
Identifying the phenomenological parameters of the model is a challenge due to the lack of the experimental data at the microscopic level. Nonetheless, developing such single-cell mathematical models of the tumour-immune competition may lead to the characterization of these parameters by theoretical methods based on methods of immunology. In this paper, we modelled the immune system as one whole population. Our model can be easily extended to include other subpopulations of the immune system with the aim of specializing the biological functions within each population. Additionally, in this model, we only considered the cytokines within the tumour microenvironment, and not as external therapeutic drug. Cytokines as natural mediators of the immune response and therefore, considered for potential therapeutic and preventive targets as well as prognostic factors. They have been recently explored in cancer immunotherapy (i.e. the use of the immune system to treat cancer) (Hong-Mei, 2014;Landskron et al., 2014). An extension of this model could be to consider the cytokines as therapeutic drug by adding an external source to the model or making the key model parameters time dependent to account for the external actions modifying the immune response.