A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression

The heterogeneity and plasticity of macrophages have become a topic of great inter- est, due to their role in various diseases ranging from cancer to bacterial infections. While initial experimental studies assumed an extreme polarisation situation, with the (anti-tumour) M1 and (pro-tumour) M2 macrophages representing the two extreme cell phenotypes, more recent studies showed a continuum of macrophages polarisation phe- notypes. Here, we focus on tumour-macrophage interactions and develop a mathematical model based on kinetic equations for active particles to describe (i) the dynamics of macrophages with a continuum of diverse functional states, ranging from pro-tumour to anti-tumour states; and (ii) the dynamics of tumour cells with a variety of progression (i.e. mutation) states. With the help of this model we show that the growth of solid tumours is associated with an increased clonal heterogeneity, as well as with an increased macrophages phenotypic heterogeneity (caused by a shift from an initial anti-tumour M1-like phenotype to a mixed M1–M2 phenotype). Moreover, we show that the assumption of exponential tumour/immune cell growth leads to an unbounded macrophages growth, which is biologically unrealistic. In contrast, the assumption of logistic tumour/immune cell growth can lead Finally, we show that tumour dormancy is associated with an increase in the clonal heterogeneity of tumour cells and in the phenotypic heterogeneity of macrophages.


Introduction
Macrophages are a very important component of the innate immunity, being present in almost all tissues. 50,55 Cells belonging to the monocyte/macrophage lineage are characterised by heterogeneity and plasticity. 42 The heterogeneity and plasticity of macrophages have become a topic of great interest in immunology, due to the implications to the treatment of various diseases including cancer, 43,55 bacterial infections, 13,48 and even liver fibrosis. 58 In solid tumours, between 5-40% of tumour mass is formed of macrophages, and these macrophages are involved in tumour progression and collective migration. 26,57 The two extreme classes of macrophages are the M1 (or classically activated) cells and the M2 (or alternatively activated) cells, and their nomenclature mirrors the Th1-Th2 paradigm -since the M1 and M2 cell reciprocally engage with the Th1 and Th2 cells. 49,58 The M1 cells are proinflammatory, anti-tumour, have anti-bacterial and anti-viral activity, and can cause tissue damage (see also Fig. 1). The M2 cells are anti-inflammatory, pro-tumour, have anti-parasitic activity and are involved in tissue remodelling (see Fig. 1). Although the view of two extreme classes of macrophages simplifies the understanding of these complex cells, many macrophages express markers of both M1 and M2 differentiation (and therefore can exhibit mixed M1 and M2 activities). Because the majority of experimental studies in the literature do not investigate the differentiation markers (these studies usually look at the total number of macrophages), we are still lacking a full understanding of the role of various types macrophages in different diseases.
Understanding the balance between M1/M2 cells (with their varying polarisation phenotypes) could eventually help predicting the outcome of diseases associated with M1 or M2 immune responses. 49 For example, advanced metastatic tumours have been associated with tumour-infiltrating macrophages that have a M2 phenotype, 28 while early tumours have been associated with macrophages that have a M1 phenotype 50 (see also Fig. 1). In fact, it has been shown that a high M1/M2 ratio was associated with extended survival of some cancer patients. 66 Moreover, the M2 macrophages are involved in the collective migration of tumour cells during tissue invasion and metastasis, through their direct communication with tumour cells via different signalling pathways. 20,26,37 Recent studies have suggested immunotherapeutic approaches (e.g. via cytokines such as IFN-γ and receptormediated activation signals) to re-polarise the M2 macrophages towards an M1 phenotype with the purpose of treating cancer. 21,28 The re-polarisation of M1 and M2 cells to a different phenotype can be influenced by the cytokine environment (cytokines produced by the macrophages them-

selves, or by other cells such as Th1 and Th2 cells, NK cells or even tumour cells).
For example, the exposure of murine macrophages to type-I cytokines (e.g. IFN-γ, IL-12) triggers their polarisation towards M1 cells. 28 In contrast, the exposure of macrophages to type-II cytokines (e.g. IL-4, IL-10, IL-13) triggers their polarisation towards M2 cells.
The outcome of tumour (immuno-)therapies does not depend only on the heterogeneity of tumour-infiltrating macrophages, but also on the clonal heterogeneity of the tumour itself. 17 This intra-tumour heterogeneity is caused not only by genetic and epigenetic instability, 15 but can be caused also by tumour dormancy. 56 Various studies have shown that tumour heterogeneity impacts not only tumour evolution and metastasis, but also drug resistance and immune response. 17 Therefore, Therefore, in order to understand and control tumour evolution, it is important to first understand the interplay between the heterogeneity of solid tumours and the heterogeneity of the immune responses (which can be either pro-tumour or anti-tumour 52 ).

Mathematical approaches.
Over the past decade, various mathematical models have been derived to investigate the interplay between the M1 and M2 macrophages phenotypes in various diseases: from tumour progression, 23,24,39,41 to heart failure, 61 lung infections, 27,44 healing of bone fractures 60 or some general inflammatory processes. 59 However, none of these models have investigated the continuum of different polarisation states between the two extremes, M1 and M2, and their role on disease progression.
The kinetic theory of active particles 6,7 can be used to model and investigate this continuum distribution of different macrophages polarisation states, and the interactions between macrophages and tumour cells. In fact, kinetic theory models have been used quite often to investigate the competition between tumour and immune cells; see, for example Refs. 3,5,8,9,10,11,12,34 and the references therein. Some of these models investigated also the interactions between tumour cells and anti-tumour macrophages, 3 but the focus was exclusively on the antitumour properties of macrophages.
The goal of this study is to develop and assess a model for the kinetic theory of active particles, which considers the dynamics of a macrophage population with a continuum of polarisation states (from the anti-tumour M1 phenotypes to the pro-tumour M2 phenotypes), and its interaction with a clonally-heterogeneous population of tumour cells. To this end, we will consider the following sequential steps that describe the modelling approach of the kinetic theory of active particles (for details see, for example, Refs. 9 and 10): (1) The cells that form the system, which are referred to as active particles, are partitioned into functional subsystems labelled with the subscripts i = 1, . . . , n.
The cell functionality is heterogeneously distributed among active particles, and corresponds to an individual state defined as activity. Hence the state of each functional subsystem is defined by a probability distribution function over the activity variable. The macroscopic quantities, such as density, mean cells' functionality, and so on, are obtained by moments of the probability distribution functions. (2) Active particles interact within the same functional subsystem as well as with particles of other subsystems. Interactions are modelled by theoretical tools of stochastic games: each player/particle interacts with the collectivity; the outcome of a single interaction event can be known only in probability; the rules of interaction can evolve in time. (3) The evolution of the probability distribution is obtained by a balance of particles (cells) within elementary volumes of the space of microscopic states, the inflow and outflow of particles (cells) being related to the aforementioned interactions.

Derivation of the Mathematical Model
To derive the mathematical model, we will go through the details of the steps defined in Sec. 1. For simplicity, in this study, we will ignore the spatial collective movement of macrophages and tumour cells (which will be the subject of a future work), and instead focus only on their temporal and phenotypic dynamics. Functional subsystems. We first define the following two subsystems: i = 1: macrophages; i = 2: tumour cells. For simplicity, we use the same activity variable, u, to describe the functionality of each of the two subsystems: • For macrophages (i = 1), u ∈ [−1, 1] describes the polarisation state of cells.
Positive values (u > 0) describe the ability of macrophages to destroy tumour  Thus, the state of the two functional subsystems is described by the following probability distribution functions: Tumour cells: Note that we choose finite domains for the variables u since we assume: (i) a limit for macrophages cellular polarisation (i.e. u = −1 for fully polarised M2 macrophages and u = 1 for fully polarised M1 macrophages); (ii) a limit for tumour progression (i.e. u = 1 corresponds to the most aggressive tumour cell types, with the highest proliferative and metastatic potential).
Interactions. It has been shown experimentally that macrophages can have both anti-tumour and pro-tumour effects, depending on their phenotype. 55 Therefore, we will assume that positive values of u lead to tumour elimination (due to the M1 phenotype of macrophages), while negative values of u lead to tumour growth (due to the M2 phenotype of macrophages). Moreover, the tumour cells have the ability to reduce the anti-tumour activity of macrophages, by shifting their phenotype from u > 0 to u < 0 (note that this is in contrast with the apoptosis effect of tumour cells on some of the adaptive immune cells, such as the T cells 29 ). Also, the increase in the number of macrophages is the result of their recruitment to the tissue in response to tumour chemoattractants. 38 We model the interactions between the macrophages and tumour cells via the following terms: progression probability density (B ij (u * → u)), proliferative rates (α, β f ) and destructive rates (γ).
Based on the above description, we can formulate the mathematical model as follows: The following parameters appear in the above equations: • λ 1 = re-polarisation rate of f 1 macrophages. Although this rate depends on the cytokine environment, 42 for simplicity, here we assume that this rate is constant; • λ 2 = progression rate of f 2 tumour cells. Throughout this study we assume that this rate is constant; • α = proliferation/recruitment rate of macrophages, in response to chemokines secreted by the tumour cells. 38 For simplicity (and since we do not hold enough information on the possible differences in the proliferation of M1 and M2 macrophages), we assume that α is constant; • μ d = half-life of macrophages. Since during steady state conditions, circulatory monocytes from which macrophages originate have a half-life of one to three days, 64 for simplicity in this study we assume that all macrophages have a similar averaged constant half-life μ d ; • β = proliferation rate of tumour cells, due to interactions between tumours and some external nutrient (β n ), as well as interactions between tumours and M2 macrophages, i.e. macrophages with activity u < 0 (β f ). Again, for simplicity (and since we could not find information about a possible correlation between the tumour proliferation rate and the progression stage), we assume that these tumour proliferation rates are constant; • γ = destruction rate of tumour cells following interaction with macrophages with M1-like phenotype (i.e. macrophages with activity variable u * > 0). Again, for simplicity, we assume that the tumour killing rate is constant. • B 12 = macrophages polarisation probability density, which is the result of interactions between macrophages and tumour cells. Experimental studies have shown that tumours can switch the macrophages polarisation phenotype: from an initial M1-like phenotype (with u > 0) to a M2-like phenotype (with u < 0). 18,28 In this study, we take the approach in Ref. 25 and consider a normalised step function Here, parameter σ 1 describes the uncertainty of changing macrophages' polarisation (i.e. phenotype): large σ 1 means less exact change in polarisation, while small σ 1 means more exact change in polarisation. We note that σ 1 depends implicitly on the level of cytokines in the environment. Hence, the function above models the probability for a macrophage with phenotype u * to change its phenotype to u = u * − k 1 (u * + u * ). Parameter k 1 is a measure of the magnitude of polarisation: k 1 = 0 no biased polarisation, while k 1 = 1 strong biased polarisation following interactions with tumour cells u * > 0 (and tumour-produced cytokines). Throughout this study, we consider k 1 ≈ 1, and thus the new phenotype will be u ≈ −u * ∈ [−1, 0] (provided that σ 1 is small). Remark 2.1. As observed experimentally, 18 the presence of the tumour always changes the phenotype of macrophages, from an initial M1 phenotype to a later M2 phenotype. Thus, it makes sense to have u ≈ −u * < 0 in the B 12 term. This means that more aggressive tumours (i.e. u * ≈ 1) will lead to macrophages with a full M2 phenotype (i.e. u ≈ −1) that support tumour progression. On the other hand, less aggressive tumours (i.e. u * ≈ 0) will lead to macrophages with a less dominant M2 phenotype (i.e. u ≈ 0).
• B 21 = tumour progression probability density, which is the result of interactions between tumours and M2-like macrophages. Experimental studies have shown that tumour-associated macrophages (with an M2-like phenotype) enhance tumour progression through the inflammatory environment they create. 54 We assume again that B 21 is given by Here, σ 2 describes the uncertainty of changes in tumour progression. This function models the probability of tumour cells with progression state u * to evolve into cells with progression state u = u * − k 2 (u * + u * ) upon interactions with macrophages with phenotype u * < 0 (i.e. M2 cells). As before, k 2 is a measure of the progression of tumour cells, and we take k 2 ≈ 1 which then leads to u ≈ −u * ∈ [0, 1] (provided that σ 2 is small). Thus, macrophages with an M2-distinct phenotype (i.e. u * ≈ −1) will lead to aggressive tumour cells (with u ≈ 1).

Remark 2.2.
The transition probability densities B 12 and B 21 satisfy the following normalisation conditions (for any u * , u * ): It is also worth noting that in this paper, attention is restricted to the so-called "linear" interactions, namely interactions which only depend on the microscopic state of the active particles. This assumption leads to a model with quadratic nonlinearities in contrast with the "structural" nonlinearities that arise when the transition probability densities also depends on the probability distribution functions of the interacting active particles. Models with "nonlinear" interactions have been extensively treated in Refs. 7 and 12.

Remark 2.3.
These interactions involve the probabilistic assessment of three types of cells (particles): test cells with microscopic (i.e. molecular) state u, candidate cells with microscopic state u * , and field cells with microscopic state u * . The candidate cells can acquire, with a certain probability, the microscopic state of the test cells, following an interaction with the field cells. Note that for cells, these proliferativedestructive-progressive interactions usually occur through a variety of molecules (i.e. cytokines, chemokines, growth factors) secreted by the cells themselves. Nevertheless, to keep the model simple, we will assume that these interactions are direct.

Results
To start investigating the dynamics of model (2.1a), in Sec. 3.1, we will focus on the steady states and note that this model cannot support nonzero equilibrium distribution functions that are constant over the activity variable. Then, in Sec. 3.2, we will consider a particular case where λ 1 = 0, and show that the model can be reduced to a classical Lotka-Volterra-type model whose dynamics is known. In Sec. 3.3, we will perform numerical simulations to support the analytical predictions. Finally, in Sec. 3.4, we will investigate numerically the effect of a logistic tumour and macrophages growth on overall system dynamics. Throughout this study, the numerical simulations are performed for two different initial conditions: • Homogeneous macrophages and tumour cells distributions across the phenotype space: The parameter values used for the numerical simulations are shown in Table 1. To obtain some of these values we used the following mouse experimental knowledge: In this study, we consider a killing rate γ ∈ (0.17, 0.75). • Regarding the proliferation of murine solid tumours (or human solid tumours implanted into mice), which are all infiltrated by macrophages, we note that there is a variety of proliferation rates: from the highly-aggressive B16F10 melanoma that has a doubling time of only 17.2 h, 22 to the less aggressive non-small cell lung cancer (NSCLC) that has a median doubling time between ≈2.1-7.9 days. 62 This corresponds to a tumour proliferation rate β ∈ (0.087, 0.967). Since in our study we split the tumour proliferation rate into β n and β f , we will assume that each of these rates are half the value of β: β n,f ∈ (0.043, 0.483).

Steady states
The search for equilibrium and local non-equilibrium states is a challenging problem of the kinetic theory of active particles. This topic has been thoroughly investigated, especially in the field of econophysics, sociodynamics, 2 and, more recently,  in biological sciences. 1 Here, we look for uniform equilibrium distribution functions of the form Substituting these values into model (2.1a), we obtain It is not difficult to prove that (f * 1 , f * 2 ) = (0, 0) is a fixed point for the steady state equations (3.2a). In regard to the existence of nonzero fixed points, it is worth noting thatB 12 (u) andB 21 (u) are not constant functions; see also Fig. 4. Therefore, equilibrium distribution functions which are constant with respect to the activity variable u do not exist for model (2.1a).

Numerical simulations
Numerical solutions of model (2.1a) are obtained by discretising the distribution functions on uniform grids. This yields a system of ordinary differential equations which are numerically solved by a first-order explicit Euler scheme. Integrals on the right-hand side of Eqs. (2.1a) are computed by means of the midpoint quadrature formula. The solution provides the distribution functions f i (t n , u j ), i = 1, 2, and related moments at different times t n and activity discrete points u j .
We start the numerical investigation into the dynamics of the kinetic model (2.1a) by comparing the cases (a) λ 1 = 0 (when the model can be reduced to Eqs. (3.5a)) and (b) λ 1 > 0. Figure 5    for λ 1 = 0. In fact, the solution trajectories eventually approach infinity -as it will be seen more clearly in the next figures.
Next, we investigate the time-dynamics of model (2.1a) in the activity space, for the baseline parameter values summarised in Table 1. In Fig. 6(a), we consider initial  conditions that are homogeneous distributions (as in Fig. 3(a)), while in Fig. 6(b), we assume initial conditions that are heterogeneous distributions (as in Fig. 3(b)). Both types of initial conditions lead to a time-periodicity in tumour-immune interactions, and for very large times the macrophage population will eventually grow unbounded (see Figs. 6(i)). Moreover, the evolution of the tumour leads to a particular cell distribution over the activity space, which is independent of the initial conditions: all macrophages eventually acquire phenotypes u < 0.5 (see Figs. 6(ii)), and all tumour cells acquire mutations u ∈ [0, 1] with slightly more cells having mutations u = 0.5 (see Figs. 6(iii)).
Note that the large growth in macrophage population, which is triggered by a relatively high proliferation rate, is not biologically realistic in the context of solid tumours. Next, we investigate the effect of decreasing macrophages' proliferation rate α and increasing their death rate μ d . Figure 7 shows the dynamics of model (2.1a) for α = 0.5 and α = 0.05, while μ d = 0.9. Even in this case the macrophage population eventually grows unbounded, triggered by a significant growth in the tumour population.
We should emphasise that changes in other parameter values lead to tumourimmune behaviours similar to those in Fig. 7 -not shown here. For example, increasing the tumour-killing rate γ (e.g. to γ = 0.75) leads to macrophage oscillations with increased amplitude, as in Fig. 7(i). Decreasing the transition rates λ 1,2 (e.g. to λ 1 = 0.001, λ 2 = 0.0001) leads to a very slow growth in the amplitude of macrophage oscillations; but eventually the amplitude of oscillations will approach infinity.

Remark 3.2.
These numerical results suggest that the non-negative solutions for model (2.1a) with exponential tumour/macrophages growth and with the initial conditions f 1,2 (0, u) = f * 1,2 (u) described at the beginning of Sec. 3 exists only locally. Note that the local existence of non-negative solutions for the kinetic model (2.1a) can be proved analytically by taking the same approach as in Ref. 4. Alternatively, if we re-write the averaged density model (3.3a) as (3.6a), then, local existence and uniqueness of non-negative solutions for this ODE can be shown using the classical Picard-Lindelöf existence theorem and Gronwall's inequality. 53 Finally, we note in Eqs. (3.6a)-(3.6b) that as soon as n 2 (t) > μ d /α, the macrophages will grow unbounded (i.e. n ± 1 → ∞ as t → ∞), which explains the above numerical results. On the other hand, if n ± 1 ≈ 0 (by choosing a very small α, such that n 2 < μ d /α), then the tumour cells will grow unbounded (i.e. n 2 (t) → ∞ as t → ∞). Therefore, the blow-up in the solutions of model (3.3a) (as well as the original model (2.1a) -see the numerical simulations in Figs. 6(i) and 7) can be triggered by both the uncontrolled growth of macrophages as well as the uncontrolled growth of tumour cells.

Logistic tumour and macrophages growth
Next, we modify model (2.1a) to incorporate a logistic growth of tumour cells caused by the availability of resources, as observed in many experimental and clinical studies that emphasise the slow-down in tumour growth at large volumes. 14,35, 36 We also assume a logistic growth of macrophages, since experimental studies have shown that macrophages growth curves displays an initial exponential growth phase followed by a stationary phase. 19 Moreover, very large macrophage populations secrete large amounts of cytokines that could eventually kill the patient (e.g. by triggering the cytokine storm syndrome 40 ). The carrying capacities for the f 1 and f 2 cell populations, namely K 1 and K 2 , are assumed to satisfy the relation K 1 = 0.4K 2 (since up to 40% of tumour mass can be formed of macrophages 57 ). Thus, the new model with logistic cell growth is given by the following equations: Remark 3.3. The above kinetic system can be reduced to the following equations for the averaged cell population densities: where n ± 1 , n 1 and n 2 are defined by Eqs. (3.6a). If we further split Eq. (3.8a) into two separate equations for n + 1 and n − 1 (with n + 1 + n − 1 = n 1 ) we have This model admits five fixed points: the trivial state (n +, * 1 , n −, * 1 , n 2 ) = (0, 0, 0), a semi-trivial tumour-only state (n +, * 1 , n −, * 1 , n 2 ) = (0, 0, K 2 ), a semi-trivial state with M2 cells infiltrating the tumour (n +, * 1 , n −, * 1 , n 2 ) = 0, that exists for αK 2 > μ d , a semi-trivial state with M1 cells infiltrating the tumour (n +, * 1 , n −, * 1 , n 2 ) = N +, * 1 , 0, and a tumour-macrophages coexistent state (n +, * 1 > 0, n −, * 1 > 0, n 2 > 0). The existence of so many steady states suggests that model (3.9a) can exhibit more complex dynamics, with solutions approaching or moving away from these states, depending on their stability (note that for the parameter values discussed in this study -see Table 1 -the coexistence state is stable while the other states are  unstable).
Model (3.7a) is discretised using the same numerical method described in the previous section. Here, we focus only on the case of heterogeneous initial cell distributions (as in Fig. 3(b)), since this case describes the more realistic situation where the initial tumour does not contain too many mutations. Figure 8 shows a few examples of tumour dynamics as we vary different model parameters. We note in particular the dormant tumour behaviour displayed in panels (a)-(c): the tumour size does not grow significantly while the tumour cells increase their clonal heterogeneity; see the n 2 curve for (a) t ∈ (300, 650), (b) t ∈ (400, 1000) and (c) t ∈ (150, 300). This result is consistent with various experimental observations of tumour phenotype plasticity during dormancy, where the accumulation of tumour mutations leads to immune escape. 47,63 The tumourmacrophages dynamics in Figs. 8(a)-8(c) shows that the tumour is controlled at very small sizes (i.e. n 2 ≈ 1 in panels (a) and (b) and n 2 ≈ 2 in panel (c)) by a much larger macrophage population. This macrophage population is formed of a mixed M1-M2 phenotype, although the macrophages f 1 (t, u) with an extreme M2-phenotype, i.e. with u ∈ [−1, −0.75], have a much lower density compared with those macrophages with u > −0.75 (see Fig. 8(a)(ii) for t = 400 and t = 600). The tumour relapse in panel (a) is associated with a large density of macrophages with extreme M2 phenotype (see Fig. 8(a)(ii) for t = 800). Moreover, at this stage, the tumour density is almost homogeneously distributed over the activity interval [0, 1], meaning that the tumour is formed of equal densities of cells with all possible clonal mutations. However, if we compare the tumour dynamics in panels (a) and (b) we note that even when the tumour is controlled at smaller sizes (as in panel (b)) there is a mixed phenotype of M1 and M2 macrophages infiltrating the tumours. Nevertheless, the level of f 1 (800, u) in panel (a)(ii) is much higher than the level of f 1 (800, u) in panel (b)(ii). This, together with the higher γ values in sub-panels (b), explains the lower tumour sizes in these sub-panels.
Another aspect that we investigated in Fig. 8(e) is the effect of decreasing the transition rates λ 1 (for macrophage re-polarisation) and λ 2 (for tumour clonal mutation). We note that lower λ 1,2 rates could allow the M1 macrophages to keep the tumour under control at very small sizes for a very large time: n 2 ≈ 0.6 for t ∈ (600, 1600).
Note that all simulations in Fig. 8 were performed under the assumption that the tumour carrying capacity is small, i.e. K 2 = 40 (and correspondingly the carrying capacity of tumour-infiltrating macrophages is small K 1 = 40%, K 2 = 10). This describes the growth of small avascular tumours, 33 up to the point where they deplete the environment from nutrients and to continue growing they need to induce  Table 1. Sub-panels (ii) and (iii) show the probability distribution of macrophage (f 1 ) and tumour (f 2 ) sub-populations over the activity space u, at five different time steps. (b) Dynamics of model (3.7a) when we increase the rate at which M1-macrophages kill tumour cells: γ = 0.5 → γ = 0.75; Sub-panel (i) shows the time-evolution of the activity-averaged solutions n 1 and n 2 , while sub-panels (ii) and (iii) show the probability distribution of macrophage (f 1 ) and tumour (f 2 ) sub-populations over the activity space u, at five different time steps. (c) Time-evolution of n 1 and n 2 as we decrease the proliferation rate of macrophages to α = 0.5 and increase the half-life of macrophages to μ d = 0.9. (d) Time-evolution of n 1 and n 2 as we decrease α even further to α = 0.05. (e) Time-evolution of n 1 and n 2 as we decrease the transition rates λ 1 = 0.01 → λ 1 = 0.001 and λ 2 = 0.01 → 0.0001. For all these simulations we consider K 1 = 10, K 2 = 40 (corresponding to some small tumour volumes) and all other parameter values as in Table 1  angiogenesis. However, increasing the tumour carrying capacity to K 2 = 400 leads to dynamics similar to the one in Fig. 8 -not shown here. All these numerical results confirm that in the presence of a logistic growth, the non-negative solutions for model (3.7a) exist globally in time, as mentioned at the end of Remark 3.3.

Summary and Discussion
Understanding the collective dynamics of solid tumours, many of which being infiltrated with large macrophage populations that have different phenotypes and functions, is an important step in the future development of new anti-cancer therapies. To shed some light on the evolution of such solid tumours, in this study, we derived a new kinetic model for active particles that focuses on the interactions between a heterogeneous population of macrophages and a heterogeneous population of tumour cells. While there are many other kinetic models for tumour-immune interactions (see, for example, Refs. 12, 16, 30 and 34), the majority of those models are quite generic, being formulated as general competition models. In this study, we aimed to derive a new model based on latest experimental studies regarding the dynamics of macrophages with different phenotypes (anti-tumour M1 and pro-tumour M2) and their specific interactions with the tumour cells. The main results have been summarised in Table 2.
We started our investigation into the dynamics of this new kinetic model by showing that it cannot display phenotypically/clonally homogeneous steady states (Fig. 4). The model could be reduced to the classical Lotka-Volterra model when assuming an initial uniform macrophage phenotype distribution and no macrophages re-polarisation (i.e. λ 1 = 0). And as expected, the reduced model exhibited an oscillatory dynamics (Fig. 5). The introduction of macrophages re-polarisation (λ 1 > 0) combined with an exponential cell growth (for macrophages and tumour cells) led to an unbounded growth of macrophages, which then triggered the growth of tumour population (see Fig. 6). For very low macrophage proliferation rates (α = 0.05 in Fig. 7) one could see an initial unbounded tumor growth, which triggers an unbounded macrophages growth.
To avoid this unrealistic unbounded growth of macrophages, and to incorporate the phenomenological observation that tumour growth slows down at large sizes, 35 dynamics (see Figs. 8(a)-8(c)), during which the tumour cells increased their clonal heterogeneity (see Fig. 8). This simulation result is consistent with various experimental studies, which emphasise the divergence in tumour clones during dormancy. 47,56,63 Understanding tumour dormancy is important because it can lead to further tumour relapse and possible patient death (in Ref. 56, the authors noted that more than 60% of all deaths from breast cancer take place after the fiveyear survival mark, which implies that dormancy might play an important role on patient survival). In our study, the tumour-dormant behaviour was mediated by the immune cells (i.e. M1-like cells), and was also characterised by a re-polarisation of M1 macrophages towards an M2 phenotype. And when the density of macrophages with the extreme M2-phenotype (i.e. u ≈ −1) increased significantly, this was associated with a tumour relapse.
Note that for large phagocytosis rates γ, it is possible that the tumour is controlled at relatively small values also in the presence of macrophages with M2phenotype (see Fig. 8(b)). This result could explain some of the contradictory experimental studies in the literature regarding the contribution of M1 versus M2 macrophages to patient survival time. For example, in Ref. 46, the authors observed that increased M1 macrophage densities inside tumour islets and stroma were associated with improved survival, but the M2 macrophage densities (even if they were higher than the M1 densities inside tumour islets and stroma) were not associated with patients survival time. In Ref. 51, the authors observed increased patient survival following the infiltration of tumour islets with M1 cells (as well as relatively large numbers of M2 cells). Extended survival was also observed for large densities of M1 and M2 cells in the tumour stroma. On the other hand, in Ref. 32, the authors suggested that while high infiltration of M1 macrophages inside tumour islets was indeed associated with increased overall survival, high infiltration of M2 macrophages (inside islets and stroma) was associated with reduced overall survival. Therefore, there is still quite a lot of debate on the pro-tumour/anti-tumour roles of M1 and M2 cells infiltrating the tumours, and our study (Fig. 8) suggests that the rates at which these cells eliminate or support tumour growth might lead to improved or reduced patient survival.
While some of the model parameters listed in Table 1 were obtained from the published biological literature, we need to emphasise that there is a lot of variability in these parameters. For example, the proliferation (α) and death rates of macrophages (μ d ) overlap over a large interval, which leads to the possibility of having α ≤ μ d or α ≥ μ d . Also the proliferation rates of the tumour cells (β n,f ) or the rates at which the M1-like macrophages can kill tumour cells (γ) can vary over large intervals, depending on the tumour cell lines. For the baseline values listed in Table 1 we chose arbitrary values from the middle of the parameter ranges. A more rigorous approach for parameter identification would have been to fit the models (2.1a) and (3.7a) to a published dataset for tumour growth in the presence/absence of macrophages with different phenotypes. Unfortunately, we are not aware of any experimental/clinical studies that quantify (through time series data) not only tumour volume growth but also the evolution of macrophage populations with intermediate phenotypes (although there are some studies that quantify tumour growth in the presence of macrophages with extreme M1 and M2 phenotypes, also quantified in terms of their percentages on days 7 and 14 after tumour injection 18 ). This lack of experimental data makes it difficult to parametrise kinetic mathematical models for active particles (as opposed to simpler ODE models that were parametrised with the tumour/macrophages data from Ref. 18; see Ref. 23). For this reason, the purpose of this current study was not to make quantitative predictions, but to try to obtain a better understanding of the qualitative dynamics exhibited by models (2.1a) and (3.7a).
Previous experimental studies have shown a co-migration of tumour cells and macrophages during the invasion and metastasis process. 20,26 Therefore, future work will focus on generalising these kinetic models to incorporate also the collective spatial movement of both macrophages and tumour cells, as a results of cell-cell communication via different signalling pathways.