# Response of Immunotherapy to Tumour-TICLs Interactions: A Travelling Wave Analysis

Authors: Joseph Malinzi; Precious Sibanda; Hermane Mambili-Mamboundou
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

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.


---

## Body

## 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 [2–9] 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, [12–14], 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 E, tumour cells T, interleukin 2 concentration IL2, tumour-immune cell complex C, a chemokine α, and resting cells R. The class IL2 represents a population of cultured immune cells that have antitumour reactivity with the tumour host. We assume that IL2 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 IL2 increases the rate of conversion of resting TICLs to primed TICLs (see [17, 18]). R 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.Figure 1
Schematic diagram of the local kinetic cell interactions [17].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 k1 which later can either lead to tumour cell death with probability p at a rate k2p or lead to inactivation of TICLs at a rate k2(1-p). In case of the latter, the tumour-immune complex is dissociated at a rate k-1. k2 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 k1ET. Similarly the TICLs and tumour cell densities, respectively, increase by (k-1+k2p)C and (k-1+k2(1-p))C, 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 χ∇·(E∇α) 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 ρR, where ρ is the supply rate. We consider the immune cells proliferation term to be fC/(g1+T) and similarly the chemokine production term to be fC/(g3+T), where f, g1, and g3 are constant parameters derived from experimental results. fC/(g1+T) 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]):
(1)∂E∂t=D1∇2E-χ∇·(E∇α)+ρR+fCg1+TWWI-d1E-k1ET+(k-1+k2p)CWWI+ωIL2·R+θ2E·IL2g2+IL2+eT,∂T∂t=D2∇2T+a1T(1-b1T)-k1ETWWI+(k-1+k2(1-p))C,∂C∂t=k1ET-(k-1+k2)C,∂α∂t=D3∇2α+fCg3+T-d3α,∂IL2dt=D4∇2IL2+s2-d2IL2,∂R∂t=D5∇2R+s3+a2R(1-b2R)WWI-ωIL2·R-ρR,
where Di, i=1,2,…,5, are diffusion coefficients of primed TICLs, tumour, IL2, α, 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 IL2. The capacity of IL2 to stimulate the production of antibodies is denoted by eT and θ2E·IL2/(g2+IL2) is a proliferation term also considered by Kirschener and Panetta [3]. It models the stimulation of TICLs by IL2 and is of the Michaelis-Menten form (see [3]). a1T(1-b1T) and a2R(1-b2R) are logistic growth terms, respectively, modelling tumour and resting cells’ growth, where ai and bi-1, i=1,2, are, respectively, the growth rates and carrying capacities, s2 is the IL2 supply, and d2, d3 are, respectively, the deactivation rates of IL2 and α. s3,   a2, and b2-1 are, respectively, the resting cells supply rate, growth rate, and carrying capacity.We consider a one-dimensional spatial domain on the interval[0,x0] 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 [0,L], where L=0.2x0 [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
(2)n·∇E=n·∇T=n·∇IL2=n·∇R=n·∇C=n·∇α=0atx=0,x=x0,E(x,0)={0,0≤x≤L,E0[1-exp⁡(-1000(x-L)2)],L≤x≤x0,R(x,0)={0,0≤x≤L,R0[1-exp⁡(-1000(x-L)2)],L≤x≤x0,IL2(x,0)=IL20,∀x∈[0,x0],C(x,0)=C0,∀x∈[0,x0],α(x,0)=0,∀x∈[0,x0],T(x,0)={T0[1-exp⁡(-1000(x-L)2)],0≤x≤L,0,L≤x≤x0.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 E, T, IL2, and R as fractions of their initial concentrations with t0=x0/D1 and x0=1 cm give
(3)∂E∂t=∇2E+ϕ1¯R+θ1¯ETη1+T-ψE-νETWWI+ω1¯IL2·R+θ2¯E·IL2η2+IL2+e¯T,∂T∂t=ϕ∇2T+β1T(1-β2T)-μ1¯ET,∂IL2∂t=ξ∇2IL2+σ2-μ2¯IL2,∂R∂t=ζ∇2R+σ3+α1R(1-α2R)WWI-ω2¯IL2·R-ϕ2¯R,
where
(4)θ1¯=θ1t0,ψ=d1t0,θ2¯=θ2t0,e¯=eT0t0E0,ω1¯=ωR0IL20t0E0,η1=gT0,η2=g2t0,β1=a1t0,β2=b1T0,μ1¯=mE0t0,σ3=s3t0R0,ν=lT0t0,ϕ=D2t0,ξ=D4t0,η3¯=g3t0,ϕ1¯=ρR0t0E0,σ2=σ2t0IL20,ζ=D5t0,μ2¯=μ2t0,ω2¯=ωIL20t0,α1=a2t0,α2=b2R0,ϕ2¯=ρt0,l=Kk2(1-p),θ1¯=fK,m=Kk2p,K=k1(k-1+k2).The boundary and initial conditions are, respectively,(5)∂E∂x(0,x)=∂IL2∂x(0,t)=∂R∂t(0,x)=∂T∂x(0,t)=0,∂E∂x(1,t)=∂R∂x(1,t)=∂IL2∂x(1,t)=∂T∂x(1,t)=0,E(x,0)={0,0≤x≤L,[1-exp⁡(-1000(x-L)2)],L≤x≤1,R(x,0)={0,0≤x≤L,[1-exp⁡(-1000(x-L)2)],L≤x≤1,IL2(x,0)=IL20,∀x∈[0,1],T(x,0)={[1-exp⁡(-1000(x-L)2)],0≤x≤L,0,L≤x≤1.


## 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 z=x-ct, where c the travelling wave speed is greater than zero (c>0), and let E~(z)=E(x,t),  T~(z)=T(x,t), IL2~(z)=IL2(x,t), and R~(z)=R(x,t). For simplicity, we drop the tildes and the system (3) is transformed into
(6)-cdEdz=d2Edz2+ϕ1¯R+θ1¯ETη1+T-ψE-νETWWWW+ω1¯IL2R+θ2¯EIL2η2+IL2+e¯T,-cdTdz=ϕd2Tdz2+β1T(1-β2T)-μ1¯ET,-cdIL2dz=ξd2IL2dz2+σ2-μ2¯IL2,-cdRdz=ζd2Rdz2+σ3+α1R(1-α2R)WWWW-ω2¯IL2R-ϕ2¯R.For simple phase space analysis, we define variables(7)E1=dEdz,T1=dTdz,IL21=dIL2dz,R1=dRdz,
and (6) are transformed into a system of autonomous first order differential equations as follows:
(8)dXdz=f(X),whereX=(E1ET1TIL21IL2R1R)∈R8,(9)f(X)=[-cE1-ϕ1¯R-θ1¯ETη1+T+ψE+νET-ω1¯IL2R-θ2¯EIL2η2+IL2-e¯TE11ϕ(-cT1-β1(1-β2T)+μ1ET)T11ξ(-cIL21-σ2+μ2IL2)IL211ζ(-cR1-σ3-α1R(1-α2R)+ω2¯IL2R+ϕ2¯R)R1],
with boundary conditions
(10)lim⁡z→-∞(E1,E,T1,T,IL21,IL2,R1,R)WW=X0=(0,E,0,0,0,IL2,0,R),lim⁡z→+∞(E1,E,T1,T,IL21,IL2,R1,R)WW=X1=(0,E,0,T,0,IL2,0,R),
where X0 and X1 correspond to the equilibrium points of the system (8).The system (8) can be regarded as an eigenvalue problem because the wave velocity c is unknown. We take c=20, 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 c 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 X0 and X1, where, after substituting parameter values in Table 1,
(11)X0≈(00.00010000.796800.001),X1≈(05.893400.798600.796800.0002).
Here, X0 and X1 are equilibrium points of the system (8). Our interest is to establish the existence of an orbit Xcon(z) that satisfies
(12)lim⁡z→-∞Xcon(z)=X0,lim⁡z→+∞Xcon(z)=X1.
The existence of such an orbit would imply that travelling wave solutions do exist [1].Table 1
Dimensional parameter values for model (1).
Parameter
Estimated value
Units
Source
a
0.18
day-1
[2]
k
1
1.3
×
1
0
-
7
day cells-1 cm
[2]
k
2
7.2
day-1
[2]
d
1
0.0412
day-1
[2]
g
2.02
×
1
0
7
cells cm-1
[2]
b
2.0
×
1
0
-
9
cells-1 cm
[2]
k
-
1
24
day-1
[2]
p
0.9997
Dimensionless
[2]
f
0.2988
×
1
0
8
day-1 cells cm-1
[2]
s
1.36
×
1
0
4
day-1 cells cm-1
[2]
D
1
1
0
-
6
cm2 day-1
[2]
D
2
1
0
-
6
cm2 day-1
[2]
θ
2
0.1245
day-1
[3]
e
0
≤
c
≤
0.005
day-1
[3]
g
2
1
0
7
cm3
[3]
d
2
10
day-1
[18]
a
2
0.0245
day-1
[18]
b
2
1
1
0
7
cell-1
[18]
ρ
6.4
×
1
0
-
6
cells-1 day-1
[18]We consider the linearization(13)dXdz=Df(X0)X,dXdz=Df(X1)X,
of the vector field f at the equilibrium points X0 and X1, respectively. From the Jacobian
(14)Df(x)=(-cA10A20A30-Iω1-ϕ1100000000Tμ1ϕ-cϕA40000001000000000-cξμ2000000100000000ω2ξ-cξA500000010),
where
(15)A1=Tν-Tθ1T+η1-Iθ2I+η2+ψ,A2=Ee¯+Eν-Eη1θ1(T+η1)2,A3=Eη2θ2(I+η2)2,A4=2Tβ1β2ϕ+Eμ1ϕ-β1ϕ-Rω1,A5=2Rα1α2ξ+Iω2+ϕ2ξ-α1ξ,
we determine the spectrum of the matrices Df(X0) and Df(X1). For parameter values in Table 1, Df(X0) has eight real eigenvalues (213.22,27.11,20,15.68,  -5.35×10-9,-7.11,-15.48,-193.22), four positive and four negative. The four positive eigenvalues imply the existence of a 4-dimensional unstable manifold Wu(X0). Similarly, Df(X1) has eight eigenvalues (213.2,27.11,20,-15.68,-0.0002,-7.11,-15.48,-193.22), three positive and five negative, implying the existence of a 5-dimensional stable manifold Ws(X1). From this result, we note that
(16)dim⁡(Wu(X0))+dim⁡(Ws(X1))=dim⁡R8+1.Equation (16) suggests that Wu(X0) and Ws(X1) 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 (IL2≠0) and without (IL2=0) treatment connecting the tumour-free equilibrium point to the cancer dormant equilibrium point. In this section we seek the minimum wave speed c. We apply the same technique used by Chahrazed [22] and Maidana and Yang [23] in determining c. 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 thatX0, the tumour-free equilibrium point of (8), must not oscillate. In other words, the eigenvalues λi corresponding to this equilibrium point must have real values; that is, λi∈R. We seek the travelling wave speed both with and without immunotherapy.


### 4.1. No Treatment Case

WithIL2=0, the tumour-free equilibrium point of the system (8) is
(17)X0=(0,E*,0,0,0,R*),whereE*=(α1-ϕ2+4α1α2σ3+(α1-ϕ2)2)ϕ2α1α2ψ,R*=α1-ϕ2+4α1α2σ3+(α1-ϕ2)22α1α2.
For the equilibrium point X0 to be biologically meaningful, E* and R* must be positive. E* and R* are positive provided that
(18)α1+4α1α2σ3+(α1-ϕ2)2≥ϕ2.
The eigenvalues λi, i=1,2,…,6, corresponding to X0 are
(19)λ1=-12c+12c2+4ψ,(20)λ2=-12c-12c2+4ψ,(21)λ3=-c+c2+4(4α1α2σ3+(α1-ϕ2)2)ξ2ζ,(22)λ4=-c-c2+4(4α1α2σ3+(α1-ϕ2)2)ξ2ζ,(23)λ5=-(α1α2cψ2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+DWW+2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D)W×(2α1α2ϕψ)-1,(24)λ6=-(α1α2cψ2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+DWW-2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D)W×(2α1α2ϕψ)-1,
where
(25)C=(4α12α22β1ϕ-α12α22c2)ψ2,D=2(α12α2μ1ϕ2-α1α2μ1ϕ2ϕ2)ψ.The first four eigenvalues (19)–(22) are real. Therefore (23) or (24) should determine the minimum wave speed which we obtain by setting
(26)2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D=0,
since we require λ5,6 to be real. Solving for c in (26) gives
(27)c=(4β1ϕ+2μ1ϕ2ϕ2α1α2ψ-2μ1ϕ2α2ψ24α1α2σ3+α12-2α1ϕ2+ϕ22μ1ϕ2α1α2ψW-24α1α2σ3+α12-2α1ϕ2+ϕ22μ1ϕ2α1α2ψ)1/2.
Substituting parameter values from Table 1 into (27) gives c≥4.176. This indicates that the minimum wave speed cmin⁡ for the tumour-immune interaction model without immunotherapy is approximately 4.176.


### 4.2. Treatment Case

WithIL2≠0, the tumour-free equilibrium points of the system (8) are
(28)X0=(0,E¯,0,0,0,IL2¯,0,R¯),whereE¯=(μ2ϕ1+ω1σ2)(η2μ2+σ2)(p1-p2+(A+P)-(B+Q))2(η2μ2ψ+ψσ2-σ2θ2)α1α2μ22≥0,IL2¯=σ2μ2,R¯=p1-p2+(A+P)-(B+Q)2α1α2μ2≥0,provided(A+P)≥(B+Q),(η2μ2+σ2)ψ>σ2θ2,p1+(A+P)-(B+Q)≥p2,whereA=4α1α2μ22σ3+α12μ22,B=2α1μ22ϕ2,P=μ22ϕ22+ω22σ22,Q=2(α1μ2ω2-μ2ω2ϕ2)σ2,p1=α1μ2,p2=μ2ϕ2+ω2σ2.The eigenvaluesλi, i=1,2,…,8, corresponding to X0 (with immunotherapy) are
(29)λ1=-c+4μ2ξ2+c22ξ,(30)λ2=-c-4μ2ξ2+c22ξ,(31)λ3=-cμ2+c2μ22+4(A+P)-(B+Q)μ2ξ2μ2ξ,(32)λ4=-cμ2-c2μ22+4(A+P)-(B+Q)μ2ξ2μ2ξ,(33)λ5=-(cη2μ2+cσ2(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2WW+((η2μ2+σ2)(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2)WWW×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2)W×2(η2μ2+σ2)-1,(34)λ6=-(cη2μ2+cσ2×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))-1/2WW-((η2μ2+σ2)(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2)WWW×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2)W×2(η2μ2+σ2)-1,(35)λ7=--Γ1+(A+P)-(B+Q)-C+D+Γ2α1α22(α1α2η2μ22ϕψ+α1α2μ2ϕψσ2-α1α2μ2ϕσ2θ2),(36)λ8=--Γ1-(A+P)-(B+Q)-C+D+Γ2α1α22(α1α2η2μ22ϕψ+α1α2μ2ϕψσ2-α1α2μ2ϕσ2θ2),
where
(37)Γ1=-α1α2cη2μ22ψ+α1α2cμ2ψσ2-α1α2cμ2σ2θ2,Γ2=J(4A+α12μ22-2α1μ22ϕ2+μ22ϕ22+ω22σ22WWWIω22σ22-2(α1μ2ω2-μ2ω2ϕ2)σ2)1/2,
where J=2(η2μ1μ22ϕϕ1+μ1ω1ϕσ22+(η2μ1μ2ω1ϕ+μ1μ2ϕϕ1)σ2).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 Γ1-(A+P)-(B+Q)-C+D+Γ2α1α2=0 and substituted parameter values in Table 1. The result gave the value c≥4.176 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 c≥2 [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.


## 4.1. No Treatment Case

WithIL2=0, the tumour-free equilibrium point of the system (8) is
(17)X0=(0,E*,0,0,0,R*),whereE*=(α1-ϕ2+4α1α2σ3+(α1-ϕ2)2)ϕ2α1α2ψ,R*=α1-ϕ2+4α1α2σ3+(α1-ϕ2)22α1α2.
For the equilibrium point X0 to be biologically meaningful, E* and R* must be positive. E* and R* are positive provided that
(18)α1+4α1α2σ3+(α1-ϕ2)2≥ϕ2.
The eigenvalues λi, i=1,2,…,6, corresponding to X0 are
(19)λ1=-12c+12c2+4ψ,(20)λ2=-12c-12c2+4ψ,(21)λ3=-c+c2+4(4α1α2σ3+(α1-ϕ2)2)ξ2ζ,(22)λ4=-c-c2+4(4α1α2σ3+(α1-ϕ2)2)ξ2ζ,(23)λ5=-(α1α2cψ2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+DWW+2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D)W×(2α1α2ϕψ)-1,(24)λ6=-(α1α2cψ2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+DWW-2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D)W×(2α1α2ϕψ)-1,
where
(25)C=(4α12α22β1ϕ-α12α22c2)ψ2,D=2(α12α2μ1ϕ2-α1α2μ1ϕ2ϕ2)ψ.The first four eigenvalues (19)–(22) are real. Therefore (23) or (24) should determine the minimum wave speed which we obtain by setting
(26)2(4α1α2σ3+(α1-ϕ2)2)α1α2μ1ϕ2ψ-C+D=0,
since we require λ5,6 to be real. Solving for c in (26) gives
(27)c=(4β1ϕ+2μ1ϕ2ϕ2α1α2ψ-2μ1ϕ2α2ψ24α1α2σ3+α12-2α1ϕ2+ϕ22μ1ϕ2α1α2ψW-24α1α2σ3+α12-2α1ϕ2+ϕ22μ1ϕ2α1α2ψ)1/2.
Substituting parameter values from Table 1 into (27) gives c≥4.176. This indicates that the minimum wave speed cmin⁡ for the tumour-immune interaction model without immunotherapy is approximately 4.176.


## 4.2. Treatment Case

WithIL2≠0, the tumour-free equilibrium points of the system (8) are
(28)X0=(0,E¯,0,0,0,IL2¯,0,R¯),whereE¯=(μ2ϕ1+ω1σ2)(η2μ2+σ2)(p1-p2+(A+P)-(B+Q))2(η2μ2ψ+ψσ2-σ2θ2)α1α2μ22≥0,IL2¯=σ2μ2,R¯=p1-p2+(A+P)-(B+Q)2α1α2μ2≥0,provided(A+P)≥(B+Q),(η2μ2+σ2)ψ>σ2θ2,p1+(A+P)-(B+Q)≥p2,whereA=4α1α2μ22σ3+α12μ22,B=2α1μ22ϕ2,P=μ22ϕ22+ω22σ22,Q=2(α1μ2ω2-μ2ω2ϕ2)σ2,p1=α1μ2,p2=μ2ϕ2+ω2σ2.The eigenvaluesλi, i=1,2,…,8, corresponding to X0 (with immunotherapy) are
(29)λ1=-c+4μ2ξ2+c22ξ,(30)λ2=-c-4μ2ξ2+c22ξ,(31)λ3=-cμ2+c2μ22+4(A+P)-(B+Q)μ2ξ2μ2ξ,(32)λ4=-cμ2-c2μ22+4(A+P)-(B+Q)μ2ξ2μ2ξ,(33)λ5=-(cη2μ2+cσ2(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2WW+((η2μ2+σ2)(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2)WWW×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2)W×2(η2μ2+σ2)-1,(34)λ6=-(cη2μ2+cσ2×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))-1/2WW-((η2μ2+σ2)(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2)WWW×(c2η2μ2+4η2μ2ψ+(c2+4ψ)σ2+4σ2θ2))1/2)W×2(η2μ2+σ2)-1,(35)λ7=--Γ1+(A+P)-(B+Q)-C+D+Γ2α1α22(α1α2η2μ22ϕψ+α1α2μ2ϕψσ2-α1α2μ2ϕσ2θ2),(36)λ8=--Γ1-(A+P)-(B+Q)-C+D+Γ2α1α22(α1α2η2μ22ϕψ+α1α2μ2ϕψσ2-α1α2μ2ϕσ2θ2),
where
(37)Γ1=-α1α2cη2μ22ψ+α1α2cμ2ψσ2-α1α2cμ2σ2θ2,Γ2=J(4A+α12μ22-2α1μ22ϕ2+μ22ϕ22+ω22σ22WWWIω22σ22-2(α1μ2ω2-μ2ω2ϕ2)σ2)1/2,
where J=2(η2μ1μ22ϕϕ1+μ1ω1ϕσ22+(η2μ1μ2ω1ϕ+μ1μ2ϕϕ1)σ2).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 Γ1-(A+P)-(B+Q)-C+D+Γ2α1α2=0 and substituted parameter values in Table 1. The result gave the value c≥4.176 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 c≥2 [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 Table1, 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 R and IL2 diffuse at the same rate as TICLs (i.e., D1=D4=D5=10-6) and used the diffusivity value 10-6 for immune cells by Matzavinos et al. [2]. We took the travelling wave speed to be c=20 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 c≥4. This is consistent with our numerical simulations for which we used c=20.Travelling wave solutions of the system (3) for different travelling wave coordinates without treatment.
(a)
(b)
(c)
(d)
(e)
(f)Travelling wave solutions of the system (3) for different travelling wave coordinates with immunotherapy.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 (β1,α1,α2,σ3,ϕ,ϕ1,μ1) 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.



---

*Source: 137015-2014-07-15.xml*