There are several cancers for which effective treatment has not yet been identified. Mathematical modelling can nevertheless point out to clinicians tumour invasion properties that should be targeted to mitigate these cancers. We present a travelling wave analysis of a tumour-immune interaction model with immunotherapy. We use the geometric treatment of an apt-phase space to establish the intersection between stable and unstable manifolds. We calculate the minimum wave speed and numerical simulations are performed to support the analytical results.

1. Introduction

In travelling wave analysis, the medium moves in the direction of propagation of the wave. Travelling wave analysis is important in tumour-immune interaction dynamics since if travelling waves exist, then we may estimate the potential with which the tumour cells invade healthy tissue [1]. Tumour-immune interaction studies have revealed a lot of information regarding cancer and cancer treatments [29] including cancer dormancy, when tumour cells remain in a quiescent state for a long period of time without metastasizing. Cancer dormancy has been attributed to tumour-immune interactions, particularly tumour infiltrating cytotoxic lymphocytes (TICLs) [2]. Travelling wave analysis could lead to an understanding of the analytical connection between model parameters and tumour invasion properties.

Most of the standard cell invasion models are related to the Fisher-Kolmogorov equation. The Fisher-Kolmogorov equation [10, 11] is the simplest macroscopic reaction-diffusion evolution equation for modelling cancer invasion just as seen in [12]. Many authors, for example, [1214], have used the Fisher-Kolmogorov equation in modelling diffusive tumours and the evolution of cancer on a macroscopic scale. Several studies have shown that this equation exhibits travelling wave solutions and the minimum wave speed for these models has been estimated (see [15, 16]). The tumour-immune interaction model presented in this paper employs the Fisher-Kolmogorov equation to model the random movement of cells. The aim of this paper is to investigate the existence of travelling wave solutions in a tumour-immune interaction model with and without immunotherapy and to estimate the minimum wave speed with which tumour cells invade healthy tissue. In this way we obtain an estimate of the strength with which a tumour invades immune cells or the ability of tumour cells to resist invasion by immune cells and also identify the tumour invasion properties in the form of parameters that should be targeted to mitigate cancer in body tissue. The work presented in this paper complements the analysis done by Mambili-Mamboundou et al. [17]. They presented similar model equations, analyzed their equilibria, and found numerical solutions. The main objective in Mambili-Mamboundou et al.’s work [17] was to ascertain the cause of cancer dormancy and investigate the effect that immunotherapy has on the response of TICLs to solid tumour invasion.

2. The Model

The model considered here was derived by Mambili-Mamboundou et al. [17]. It subdivides the cell population into local concentrations of primed TICLs , tumour cells , interleukin 2 concentration , tumour-immune cell complex , a chemokine , and resting cells . The class represents a population of cultured immune cells that have antitumour reactivity with the tumour host. We assume that does not necessarily bind with TICLs to form a cell complex but rather stimulates the TICLs to fight cancer through lymphocyte activation, growth, and differentiation. We also assume that increases the rate of conversion of resting TICLs to primed TICLs (see [17, 18]). is a class representing the population of TICLs which have not yet matured or been activated by antigens. During a tumour attack on immune cells or any other body tissue infection, naive or resting TICLs are primed by antigen presenting cells (APCs) in secondary lymphoid organs such as lymph nodes and spleen [19]. Figure 1 shows the cells’ local kinetics.

Following the receptor-ligand kinetics theory in [20], when a tumour cell and an immune cell come into contact, it may lead to the formation of a tumour-immune complex at a binding rate which later can either lead to tumour cell death with probability at a rate or lead to inactivation of TICLs at a rate . In case of the latter, the tumour-immune complex is dissociated at a rate . is a parameter describing the detachment rate of TICLs from tumour cells, resulting in an irreversible programming of the tumour cells for lysis. Complex formation reduces both TICLs and tumour cell densities and increases the complex density by . Similarly the TICLs and tumour cell densities, respectively, increase by and , in case the tumour or immune cell dies. The binding of the primed TICLs to tumour cells leads to the production of a chemokine   . The chemokine gradient defines the migration of the TICLs towards the tumor by a process known as chemotaxis which is represented by in the model, with being the chemotaxis constant. We assume that the rate of supply of immune cells into the region of tumour localization is , where is the supply rate. We consider the immune cells proliferation term to be and similarly the chemokine production term to be , where , , and are constant parameters derived from experimental results. is a function that explains how tumour cells proliferate as a result of interaction with immune cells. We consider that all cell densities diffuse at constant rates. We thus consider the following system of parabolic nonlinear partial differential equations (Mambili-Mamboundou et al. [17]): where , , are diffusion coefficients of primed TICLs, tumour, , , and resting cell densities, respectively, and is the rate of stimulation of resting cells into activated TICLs as a result of injecting a patient with . The capacity of to stimulate the production of antibodies is denoted by and is a proliferation term also considered by Kirschener and Panetta [3]. It models the stimulation of TICLs by and is of the Michaelis-Menten form (see [3]). and are logistic growth terms, respectively, modelling tumour and resting cells’ growth, where and , , are, respectively, the growth rates and carrying capacities, is the supply, and , are, respectively, the deactivation rates of and . ,    , and are, respectively, the resting cells supply rate, growth rate, and carrying capacity.

We consider a one-dimensional spatial domain on the interval and assume that there are two regions in this interval, one fully occupied by tumour cells and the other fully occupied by TICLs (both activated and resting). We propose that the initial interval of tumour localization is , where [2]. In our model, we do not include the Heaviside function since we consider a resting cell class. We further assume that these resting cells can be recruited into the activated cell class. The boundary and initial conditions therefore are

It has been shown that chemotaxis does not influence the existence of travelling wave solutions (see, e.g., [1]). We therefore do the travelling wave analysis without the effect of chemotaxis. Assuming that the formation of cellular conjugates occurs on a time scale of a few hours while that of tumour cells as well as the influx of immune cells into the spleen occurs on a much slower time scale, probably tens of hours, and nondimensionalizing the above system of (1) by taking , , , and as fractions of their initial concentrations with and  cm give where

The boundary and initial conditions are, respectively,

3. Travelling Wave Solutions

In this section we investigate whether model (3) exhibits travelling wave solutions or not. We use the geometric treatment of an apt-phase space to establish the intersection between stable and unstable manifolds, a method also employed by Bellomo et al. [1] in investigating travelling wave solutions. The gist of this method is to establish the presence of a heteroclinic orbit joining two different equilibrium points in the phase space. We specify a travelling coordinate , where the travelling wave speed is greater than zero ( ), and let ,   , , and . For simplicity, we drop the tildes and the system (3) is transformed into

For simple phase space analysis, we define variables and (6) are transformed into a system of autonomous first order differential equations as follows: with boundary conditions where and correspond to the equilibrium points of the system (8).

The system (8) can be regarded as an eigenvalue problem because the wave velocity is unknown. We take , a value that numerically gives rise to travelling wave solutions. We chose this value after simulating the system of (3). There are several other values of that can give rise to travelling wave solutions. In the next section, we calculate a critical wave speed below which travelling wave solutions do not exist. We find a heteroclinic connection between and , where, after substituting parameter values in Table 1, Here, and are equilibrium points of the system (8). Our interest is to establish the existence of an orbit that satisfies The existence of such an orbit would imply that travelling wave solutions do exist [1].

We consider the linearization of the vector field at the equilibrium points and , respectively. From the Jacobian where we determine the spectrum of the matrices and . For parameter values in Table 1, has eight real eigenvalues ( ,   ), four positive and four negative. The four positive eigenvalues imply the existence of a 4-dimensional unstable manifold . Similarly, has eight eigenvalues ( , ), three positive and five negative, implying the existence of a 5-dimensional stable manifold . From this result, we note that

Equation (16) suggests that and intersect transversally along a one-dimensional curve in the eight-dimensional phase space. This is because the solutions of the system (8) lie in eight dimensions (8D) but the summation of the dimension of the stable and unstable manifolds is nine (9D) just as shown in (16) (see [1, 21]). If this is the case, then this curve would define a generic heteroclinic connection [1]. This therefore confirms that the system (1) exhibits travelling wave solutions for certain parameter values.

4. Minimum Wave Speed

In the previous section, we established that (3) exhibits travelling wave solutions. In this section, we calculate the minimum wave speed for model (3) with ( ) and without ( ) treatment connecting the tumour-free equilibrium point to the cancer dormant equilibrium point. In this section we seek the minimum wave speed . We apply the same technique used by Chahrazed [22] and Maidana and Yang [23] in determining . This technique involves analyzing the phase space by characterizing the equilibrium points of the autonomous system. The minimum wave speed corresponds to a change in the eigenvalues of the travelling-wave differential equations at the equilibrium point ahead of the wave.

To calculate the minimum wave speed, we impose a condition that , the tumour-free equilibrium point of (8), must not oscillate. In other words, the eigenvalues corresponding to this equilibrium point must have real values; that is, . We seek the travelling wave speed both with and without immunotherapy.

4.1. No Treatment Case

With , the tumour-free equilibrium point of the system (8) is For the equilibrium point to be biologically meaningful, and must be positive. and are positive provided that The eigenvalues , , corresponding to are where

The first four eigenvalues (19)–(22) are real. Therefore (23) or (24) should determine the minimum wave speed which we obtain by setting since we require to be real. Solving for in (26) gives Substituting parameter values from Table 1 into (27) gives . This indicates that the minimum wave speed for the tumour-immune interaction model without immunotherapy is approximately .

4.2. Treatment Case

With , the tumour-free equilibrium points of the system (8) are

The eigenvalues , , corresponding to (with immunotherapy) are where where .

The first six eigenvalues (29)–(34) are real provided that the conditions we have imposed for positivity of the tumour-free equilibrium point are fulfilled. Equation (35) or (36) should therefore determine the conditions for the existence of a minimum wave speed. We set and substituted parameter values in Table 1. The result gave the value as in the case without treatment. This implies that the minimum wave speed for model (3) for both with and without treatment is the same. In other words, immunotherapy may possibly not influence the strength with which the tumour cells attack immune cells because the minimum wave speed with or without clinical treatment is the same. Fisher’s equation exhibits travelling wave solutions for [24]. The minimum wave velocity which we obtained is greater than two and therefore not a violation of the minimum wave speed for Fisher’s equation.

5. Numerical Simulations

Using the parameter values in Table 1, we simulate model (3). These parameter values were obtained from data where the murine B cell lymphoma was used as an experimental model of tumour dormancy in mice [25]. The kinetic parameter values that were obtained in this experiment are shown in Table 1. We assumed that and diffuse at the same rate as TICLs (i.e., ) and used the diffusivity value for immune cells by Matzavinos et al. [2]. We took the travelling wave speed to be and implemented the simulations in python using a Runge-Kutta numerical method. The numerical simulations (see Figures 2 and 3) indicate that the system of (3) exhibits travelling wave solutions for certain parameter values. This supports the analytical results on the existence of travelling waves in Section 3. Figures 2 and 3, respectively, show the numerical travelling wave solutions for model (3) without and with clinical treatment, and for different travelling wave coordinates. They depict solutions that are periodic and oscillating around a stable equilibrium state. These solutions describe heterogeneous cell distributions with a relatively low tumour cell density. The travelling wave solutions indicate that tumour cells invade immune cells at a high potential. The minimum wave speed obtained in the previous section indicated that the model exhibits travelling wave solutions for . This is consistent with our numerical simulations for which we used .

6. Conclusions

Many biological and physical phenomena can be described by reaction-diffusion equations. However not many nonlinear reaction-diffusion equations are integrable. It is therefore imperative to find other quantitative methods for tackling such nonlinear systems. The objective of this study was to use a quantitative method to investigate travelling wave solutions of a tumour-immune interaction model and also identify the tumour invasion properties in the form of parameters that should be targeted to mitigate cancer by estimating the minimum wave speed. We investigated the existence of travelling wave solutions and estimated the minimum wave speed of the wave solutions by analyzing the model phase space. The existence of travelling wave solutions confirmed that a tumour attacks immune cells at full potential. The expression from which the minimum wave speed was calculated determined the parameters that need to be targeted to eradicate cancer in body tissue. We simulated model (3) and compared the results to analytical results. The numerical travelling wave solutions depicted periodic cell densities with a low tumor level, oscillating about a stable equilibrium state. These solutions depict cancer dormancy which has been observed in several cancers, for example, osteogenic sarcomas, basal-cell carcinoma, and breast cancers, and they also imply that the tumour cells attack the immune cells at their full potential.

Equation (27) highlights the main parameters ( ) involved in tumour invasion corresponding to tumour growth rate, resting TICLs’ growth rate, carrying capacity of the resting TICLs, resting cells’ supply, diffusion rate of the tumour cells, and the local kinetic interaction parameters (tumour cell death and inactivation of TICLs).

The results obtained in this paper are similar to those in Matzavinos and Chaplain [26]. In their work, they performed a travelling wave analysis of a model describing the growth of a tumour in the presence of an immune system response. Their results showed that indeed a tumor attacks immune cells at full potential since their model exhibited travelling wave solutions. In the future, we hope to consider diffusion in higher dimension due to the fact that body tissue geometry is highly intricate.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors’ Contribution

Precious Sibanda and Hermane Mambili-Mamboundou are coauthors.


The authors are grateful for financial support from the University of KwaZulu-Natal.