# Stability Analysis of Mathematical Model including Pathogen-Specific Immune System Response with Fractional-Order Differential Equations

Authors: Bahatdin Daşbaşı
Journal: Computational and Mathematical Methods in Medicine
 (2018)

---

## Abstract

In this study, the mathematical model examined the dynamics between pathogen and specific immune system cells (memory T cells) for diseases such as chronic infection and cancer in which nonspecific immune system cells are inadequate to destroy the pathogen and has been suggested by using a system of the fractional-order differential equation with multi-orders. Qualitative analysis of the proposed model reveals the equilibrium points giving important ideas about the proliferation of the pathogen and memory T cells. According to the results of this analysis, the possible scenarios are as follows: the absence of both pathogen and memory T cells, only the existence of pathogen, and the existence of both pathogen and memory T cells. The qualitative analysis of the proposed model has expressed the persistent situations of the disease where the memory T cells either do not be able to respond to the pathogen or continue to exist with the disease-causing pathogen in the host. Results of this analysis are supported by numerical simulations. In the simulations, the time-dependent size of the tumor population under the pressure of the memory T cells was tried to be estimated.


---

## Body

## 1. Introduction

For three centuries, the theory of fractional derivatives was developed as a pure theoretical field of mathematics, useful only for mathematicians. But, the use of fractional-orders differential and integral operators in mathematical models has become increasingly common of late years. Therefore, various forms of fractional-order differential equations are suggested for standard models. In this sense, the fractional-order calculus plays an important role in physics [1], thermodynamics [2], viscoelasticity [3], electrical circuits theory [4], fractances [5], mechatronics systems [6], signal processing [7], chemical mixing [8], chaos theory [9], engineering [10], biological system [11], and other applications [12]. Also, a large number of literatures on the application of fractional-order differential equations (FODEs) in nonlinear dynamics have been improved. Especially, when the biological applications of FODEs have considered, it is a rich source for mathematical ideas [13].The mathematical modeling of diseases in biological applications is a subject discussed in the literature. Such models are considered under two main headings as by modeling the size of the spread of infected individuals in a population (SIR) and modeling the population size of the pathogens such as the tumor in an individual as it is here. The word tumor simply refers to a mass. This is a general term that can refer to benign (generally harmless) or malignant (cancerous) growths. Many types of tumors are considered to be a major factor in many fatal diseases in human history. Fundamentally, it is said that this disease is a complex process for both tumor and host. Although different treatment strategies are proposed for tumors, the first and foremost role in disease progression belongs to the immune system of the individual (or host) [14]. The immune system is stated as a system of biological structures and processes in an organism that protects the body from the possible hazardous organism by recognizing and responding to antigens. In more detail, the immune system cells such as T cells are generally described in terms of two different types. These are the effector and the memory of T cells. The ordinary behaviour of the immune system is generally an acute infection, controlled initially by effector T cells (aspecific response or the innate immune system response), later by memory T cells (specific response or the adaptive immune response), until complete clearance of the pathogen. The adaptive immune response is much slower to respond to threats and infections than the innate immune response, which is primed and ready to fight at all times [15]. Especially, T cells are a crucial component of the adaptive immune response against malignancies. Antigen-experienced T cells specific for tumor antigens can be recovered from the blood, lymphoid organs, and tumors of both cancer patients and tumor-bearing mice [16]. Concordantly, the reactions of different hosts in case of the same disease may be different because of the immune system response given by host, and so, the disease progression varies from person to person. Within this context, dynamics of relevances between immune systems cells (T cells) and tumor are significant to find out the nature of the disease. The problem is to try to obtain the known biological features without making the mathematics too complicated.The basic of a most useful explanation of fractional calculus is memory concept. If the output of a system at each timet depends only on the input at time t, then such systems are said to be memoryless systems. On the contrary, if the system has to remember previous values of the input in order to determine the current value of the output, then such systems are said memory systems [17, 18]. Accordingly, the behaviour of most biological systems has memory or after-effects. The modeling of these systems by FODEs has more advantages than classical integer-order modeling, in which such effects are neglected. Also, FODEs are, at least, as stable as their integer order counterpart [11]. In the process of modeling real-life situations, the created mathematical models by using the fractional-order differential operations allow to display the some extra cases regarding the stability region of the equilibrium point of the mathematical model caused by parameters such as derivative orders. For this reason, the mathematical models formed by FODEs are more realistic and feasible [19]. Additionally, stability analysis of equilibrium points for mathematical models consisting of FODEs with multi‐orders and its systems is more general than those of the same-orders too.Although there are many studies that examined the dynamics between tumor and immune system response, the proposed model in this study differs from them in terms of both mathematical structure such as the use of Holling function type-2 (functional and numerical responses) in the model consisting of the FODE system with multi-orders and examination of qualitative analysis of the proposed model. In this sense, it was tried to bring a different perspective from the previous studies.In this study, a FODE model with multi-orders considering the basic mechanisms of tumor and the memory T cells having functional and numerical responses, respectively, has been constructed, and so, the qualitative analysis of the proposed model was performed. The reason for using the Holling function type-2 is to show the limit cycle behaviour of system [20]. The certain conditions dependent on the development of the tumor population under the pressure of memory T cells was obtained. In this respect, all of the possible scenarios related to the tumor size were tried to be explained as parameter-dependent. Additionally, numerical analysis of the model was given as to be compatible with the qualitative analysis.


## 2. Preliminaries and Definitions

In here, the main definitions and properties of fractional derivative operators have been expressed. Also, the FODE systems with multi-orders have been introduced, and the properties such as stability and existence of the equilibrium points of such systems are given.


### 2.1. Fractional Differential Operators

There are various definitions of a fractional derivative with the orderα>0. The definitions of Riemann–Liouville and Caputo are used most widely. The Caputo sense was used in this study. Taking into account the definition of Caputo sense, the fractional derivative of the function ft is identified as(1)Dαft=Jm−αDmft=1Γm−α∫0tfmτt−τα−m+1dτ,for m−1<α≤m, m∈N, t>0 [21].


### 2.2. The FODE System with Multi-orders

Let us consider thatt is the time parameter. We have assumed that the system of FODE with multi-orders is given as the following equation:(2)DtαXt=Ft,X,X0=X0,where the variable X=x11t,x21t,…,xn1tT∈Rn, the initial conditions by X0=x100,x200,…,xn00T∈Rn, the functions by F=f1,f2,…,fnT∈Rn and fi:0,+∞xRn⟶R for i=1,2,…,n, and the derivative orders by α=α1,α2,…,αnT.Also, when it is considered asDtα=Dtα1,Dtα2,…,DtαnT, Dtαi indicates αi th-order fractional derivative in the Caputo sense. In this sense, it is DtαXt=Dtα1x11t,Dtα2x21t,…,Dtαnxn1tT. The multi-orders can be mathematically any real or complex vector. In this study, the real case was only taken into account. Throughout the paper, we restrict αi to a rational number in the interval 0,1 [22].Remark 1.
From (2), we have assumed that(3)Ft,X=FX,where the independent variable t is not clearly seen in the function F. The equilibrium point of (3) is the point X¯=x1¯,x2¯,…,xn¯ obtained from the equations FX¯=0.Remark 2.
For each equilibrium pointX¯ of the autonomous system in (3), the eigenvalues λ obtain from the following equation:(4)detdiagλmα1,λmα2,…,λmαn−JX¯=0,where JX¯ is the Jacobian matrix evaluated at the equilibrium point and m is the smallest of the common multiples of the denominators of the rational numbers α1,α2,…,αn [23].Theorem 1.
For each equilibrium pointX¯ of system (3), X¯ is locally asymptotically stable (LAS), if the eigenvalues λi obtained from (4) satisfy Routh–Hurwitz Stability Criteria or the inequalities argλi>π/2m for i=1,2,…,mα1+α2. Here, m has been defined in Remark 2 [22]. Because the 2-dimensional of system (3) is used in the proposed model in this study, the stability analysis of such systems are described in detail below.Remark 3.
Let us assume that the autonomous system of FODE with multi-orders is as following:(5)Dα1x1t=f1x1,x2,Dα2x2t=f2x1,x2,with the nonnegative initial conditions(6)x10=xo1andx20=xo2,where the derivative orders α1 and α2 are rational numbers in the interval 0,1. The equilibrium point of system (5) is the point X¯=x1¯,x2¯ obtained from the equations Dαixit=0 for =1,2. To evaluate locally asymptotically stability (LAS) of equilibrium point, the Jacobian matrix, J=∂f1/∂x1∂f1/∂x2∂f2/∂x1∂f2/∂x2=f1x1f1x2f2x1f2x2, is used. Considering Remark 2, the eigenvalues λi for i=1,2,…,mα1+α2 are obtained from the following equation:(7)detdiagλmα1,λmα2−JX¯=λmα1−f1x1x1¯,x2¯−f1x2x1¯,x2¯−f2x1x1¯,x2¯λmα2−f2x2x1¯,x2¯=0.
Therefore, the characteristic equation for eigenvalues is(8)λmα1+α2−λmα1f2x2x1¯,x2¯−λmα2f1x1x1¯,x2¯+f1x1x1¯,x2¯f2x2x1¯,x2¯−f1x2x1¯,x2¯f2x1x1¯,x2¯=0.
If the eigenvaluesλi for i=1,2,…,mα1+α2 satisfy Routh–Hurwitz stability criteria or the conditions,(9)argλi>π2m,then the equilibrium point x1¯,x2¯ is the LAS point for system (5).
For the system of FODE with multi-ordersα1 and α2, the stability region is as shown in Figure 1 (where σ and ω are the real and imaginary parts of the eigenvalues, respectively, and j=−1). By Figure 1, we openly see that the stability region of the equilibrium point of the FODE with multi-orders is greater than the stability regions of the integer-order case and the same fractional-order case [24].Figure 1
Stability region of the equilibrium pointX¯ of system (5).Remark 4.
Letα1=α2=α in system (5). In this case, we have the system(10)Dαx1t=f1x1,x2,Dαx2t=f2x1,x2,with the nonnegative initial conditions x10=xo1andx20=xo2. From the equations Dαxit=0 for i=1,2, we have presumed that the equilibrium point of system (10) is X¯=x1¯,x2¯. If the eigenvalues λ1 and λ2 obtained from the equation(11)DetJx1,x2=x1¯,x2¯−λI2=0provide the conditions(12)argλ1>απ2,argλ2>απ2,then the equilibrium point x1¯,x2¯ is the LAS point for system (10).
Conditions expressed in (12) can be detailed as the followings. Characteristic equation of (11) is the following generalized polynomial:(13)pλ=λ2+a1λ+a2=0.
When both the conditions (12) and the polynomial (13) are considered together, the conditions for LAS of the equilibrium point x1¯,x2¯ are either Routh–Hurwitz conditions [25, 26]:(14)a1,a2>0,or(15)a1<0,4a2>a12,tan−14a2−a12a1>απ2.


## 2.1. Fractional Differential Operators

There are various definitions of a fractional derivative with the orderα>0. The definitions of Riemann–Liouville and Caputo are used most widely. The Caputo sense was used in this study. Taking into account the definition of Caputo sense, the fractional derivative of the function ft is identified as(1)Dαft=Jm−αDmft=1Γm−α∫0tfmτt−τα−m+1dτ,for m−1<α≤m, m∈N, t>0 [21].


## 2.2. The FODE System with Multi-orders

Let us consider thatt is the time parameter. We have assumed that the system of FODE with multi-orders is given as the following equation:(2)DtαXt=Ft,X,X0=X0,where the variable X=x11t,x21t,…,xn1tT∈Rn, the initial conditions by X0=x100,x200,…,xn00T∈Rn, the functions by F=f1,f2,…,fnT∈Rn and fi:0,+∞xRn⟶R for i=1,2,…,n, and the derivative orders by α=α1,α2,…,αnT.Also, when it is considered asDtα=Dtα1,Dtα2,…,DtαnT, Dtαi indicates αi th-order fractional derivative in the Caputo sense. In this sense, it is DtαXt=Dtα1x11t,Dtα2x21t,…,Dtαnxn1tT. The multi-orders can be mathematically any real or complex vector. In this study, the real case was only taken into account. Throughout the paper, we restrict αi to a rational number in the interval 0,1 [22].Remark 1.
From (2), we have assumed that(3)Ft,X=FX,where the independent variable t is not clearly seen in the function F. The equilibrium point of (3) is the point X¯=x1¯,x2¯,…,xn¯ obtained from the equations FX¯=0.Remark 2.
For each equilibrium pointX¯ of the autonomous system in (3), the eigenvalues λ obtain from the following equation:(4)detdiagλmα1,λmα2,…,λmαn−JX¯=0,where JX¯ is the Jacobian matrix evaluated at the equilibrium point and m is the smallest of the common multiples of the denominators of the rational numbers α1,α2,…,αn [23].Theorem 1.
For each equilibrium pointX¯ of system (3), X¯ is locally asymptotically stable (LAS), if the eigenvalues λi obtained from (4) satisfy Routh–Hurwitz Stability Criteria or the inequalities argλi>π/2m for i=1,2,…,mα1+α2. Here, m has been defined in Remark 2 [22]. Because the 2-dimensional of system (3) is used in the proposed model in this study, the stability analysis of such systems are described in detail below.Remark 3.
Let us assume that the autonomous system of FODE with multi-orders is as following:(5)Dα1x1t=f1x1,x2,Dα2x2t=f2x1,x2,with the nonnegative initial conditions(6)x10=xo1andx20=xo2,where the derivative orders α1 and α2 are rational numbers in the interval 0,1. The equilibrium point of system (5) is the point X¯=x1¯,x2¯ obtained from the equations Dαixit=0 for =1,2. To evaluate locally asymptotically stability (LAS) of equilibrium point, the Jacobian matrix, J=∂f1/∂x1∂f1/∂x2∂f2/∂x1∂f2/∂x2=f1x1f1x2f2x1f2x2, is used. Considering Remark 2, the eigenvalues λi for i=1,2,…,mα1+α2 are obtained from the following equation:(7)detdiagλmα1,λmα2−JX¯=λmα1−f1x1x1¯,x2¯−f1x2x1¯,x2¯−f2x1x1¯,x2¯λmα2−f2x2x1¯,x2¯=0.
Therefore, the characteristic equation for eigenvalues is(8)λmα1+α2−λmα1f2x2x1¯,x2¯−λmα2f1x1x1¯,x2¯+f1x1x1¯,x2¯f2x2x1¯,x2¯−f1x2x1¯,x2¯f2x1x1¯,x2¯=0.
If the eigenvaluesλi for i=1,2,…,mα1+α2 satisfy Routh–Hurwitz stability criteria or the conditions,(9)argλi>π2m,then the equilibrium point x1¯,x2¯ is the LAS point for system (5).
For the system of FODE with multi-ordersα1 and α2, the stability region is as shown in Figure 1 (where σ and ω are the real and imaginary parts of the eigenvalues, respectively, and j=−1). By Figure 1, we openly see that the stability region of the equilibrium point of the FODE with multi-orders is greater than the stability regions of the integer-order case and the same fractional-order case [24].Figure 1
Stability region of the equilibrium pointX¯ of system (5).Remark 4.
Letα1=α2=α in system (5). In this case, we have the system(10)Dαx1t=f1x1,x2,Dαx2t=f2x1,x2,with the nonnegative initial conditions x10=xo1andx20=xo2. From the equations Dαxit=0 for i=1,2, we have presumed that the equilibrium point of system (10) is X¯=x1¯,x2¯. If the eigenvalues λ1 and λ2 obtained from the equation(11)DetJx1,x2=x1¯,x2¯−λI2=0provide the conditions(12)argλ1>απ2,argλ2>απ2,then the equilibrium point x1¯,x2¯ is the LAS point for system (10).
Conditions expressed in (12) can be detailed as the followings. Characteristic equation of (11) is the following generalized polynomial:(13)pλ=λ2+a1λ+a2=0.
When both the conditions (12) and the polynomial (13) are considered together, the conditions for LAS of the equilibrium point x1¯,x2¯ are either Routh–Hurwitz conditions [25, 26]:(14)a1,a2>0,or(15)a1<0,4a2>a12,tan−14a2−a12a1>απ2.


## 3. Model Formulation

The proposed model is particularly well suited for describing diseases such as chronic infection and cancer in which the nonspecific immune system cells are inadequate to destroy the pathogen. Consequently, it has been proposed, and another extension of the models in [20, 27–33] has been analyzed.It has been identified mathematically as a pathogen load, specifically tumor population, and level of memory T cells, namely, the adaptive immune response or specific response, in an individual. In this sense, we have assumed that the population densities of pathogen and memory T cells at timet are denoted by Pt and Tt, respectively. Additionally, the memory T cells predate the tumor cells by a Holling function type-2.Under the assumptions aforementioned, we have proposed the following system of FODE with multi-ordersα1 and α2:(16)Dtα1P=βPP1−PΛ−cP1+aPT,Dtα2T=μP1+νPT−δT,0<α1,α2≤1,where Dtαi for i=1,2 indicates αith-order fractional derivatives in the Caputo sense, it is P=Pt and T=Tt, and the parameters have the following properties:(17)βP,Λ,c,a,μ,ν,δ∈R+.In addition that, system (16) has to be finished with positive initial conditions Pt0=P0 and Tt0=T0. The parameters used in the model are defined as follows.It is presumed that the pathogen follows a logistic growth rule with the carrying capacityΛ and the growth rate βP. The memory T cells proliferate proportionally to the pathogen load by the Holling function type-2. Since the pathogen capture rate of memory T cells is assumed to be proportional to the per capita growth rate of memory T cells, the constant μ represents the maximum growth rate for memory T cells and the constant ν is the pathogen population size at which the growth rate of memory T cells in half of its maximum. These situations are very suitable for the growth of memory T cells especially in case of chronic infection or tumor. Memory T cells have per capita natural death rate δ. Moreover, the pathogen die due to the action of the memory T cells, and we have presumed that the effect of these cells on pathogen is modeled using a saturating response, cP/1+aP, subject to a maximum killing rate c and the level of memory T cells required for the half maximum effect, a.Remark 5.
Rate of replication of the immune system cells is higher than its death rate, at least every time the pathogen load is very high [20]. In this case, we have limt⟶∞μPt/1+νPt=μ/ν by (16). Therefore, the following inequality is obtained:(18)μν>δ,by this limit.Proposition 1.
System (16) provides the followings. The free-disease equilibrium point E00,0 and the equilibrium point E1Λ,0, where only the pathogen exists and always exist. In addition to E0 and E1, there exists a third equilibrium point as E2=P∗,T∗ for(19)P∗=δμ−νδ,T∗=βP1−P∗/Λc/1+aP∗,when δ/μ−νδ<Λ.Proof.
The steady states of the model (16) are again the intersection of null clines Dα1P=0 and Dα2T=0 in (16). We have accepted that the solutions of theses equations consist of the pairs P¯,T¯. Then, we have the following system:(20)P¯βP1−P¯Λ−c1+aP¯T¯=0,T¯μP¯1+νP¯−δ=0.
From the first equation of (20), it is P¯=0 or βP1−P¯/Λ−c/1+aP¯T¯=0. Let P¯=0, and then T¯=0. Therefore, the system (16) has the free-disease equilibrium point E00,0. On the contrary, let βP1−P¯/Λ−c/1+aP¯T¯=0, that is, T¯=βP1−P¯/Λ/c/1+aP¯. If the value T¯ is rewritten in the second equation of system (20), then we have found the equilibrium points E1Λ,0 and E2δ/μ−νδ,βP1−δ/μ−νδ/Λ/c/1+aδ/μ−νδ. Considering (19), if E2 is rewritten, the point E2P∗,T∗ is obtained. E1 always exists due to (17). Let us consider E2. P∗ is positive due to (17) and (18). On the contrary, T∗ is positive due to (17), when(21)δμ−νδ<Λ.
Therefore, we have a positive equilibrium pointE2P∗,T∗ where P∗ and T∗ are in (19).
In Table1, biological existence conditions of equilibrium points of system (16) are showed.Table 1
Biological existence conditions for the equilibria of system (16).
Equilibrium points
Biological existence conditions
E
0
0,0
Always exists
E
1
Λ
,
0
E
2
P
∗
,
T
∗
δ
/
μ
−
ν
δ
<
ΛProposition 2.
In system (16), let us consider derivative orders as(22)α1=k1m1,α2=k2m2andk1,k2,m1,m2∈Z+,where the smallest common multiple of m1 and m2 is m. System (16) satisfies the following:(a)
E00,0 is a unstable point(b)
E1Λ,0 is LAS, when δ/μ−νδ>Λ. Also, if δ/μ−νδ≤Λ, then this point is a unstable point(c)
E2=P∗,T∗ where P∗ and T∗ defined in (19) is LAS, when all roots λi for i=1,2,…,mα1+α2 found from the equation(23)λmα1+α2−λmα2βPA1A2+1−1+βPδ2A1μΛ1−A1=0satisfy Routh–Hurwitz stability criteria or the condition argλi>1/mπ/2. In here, it is(24)1−P∗Λ=A1,aP∗aP∗+1=A2.Proof.
For the stability analysis of the equilibrium points, the functions in system (16) are assigned as(25)fP,T=PβP1−PΛ−c1+aPT,gP,T=TμP1+νP−δ.
In this respect, the Jacobian matrix evaluated at each equilibrium point showed in Table1 is(26)JEiP¯,T¯=βP1−2P¯Λ+T¯caP¯+1aP¯aP¯+1−1−cP¯aP¯+1T¯μ1+νP¯2μP¯1+νP¯−δ,for i=0,1,2.(a)
ForE0, the Jacobian matrix in (26) is(27)JE00,0=βP00−δ.
Now, we will investigate whether the inequality (9) has been achieved. By (27), the eigenvalues have been found from the following determinant:(28)detdiagλmα1,λmα2−JE00,0=detλmα1−βP00λmα2+δ=0.
In this respect, the characteristic equation obtained from (28) is(29)λmα1−βPλmα2+δ=0.
Thus, we haveλmα1=βP and λmα2=−δ. λmα1 is real positive due to (17). Moreover, it is obtained as(30)λi=βP1/mα1∈R+fori=1,2,…,mα1.
For the eigenvalues in (30), it is argλi=0 for i=1,2,…,mα1. These eigenvalues are positive real number on the right side of the complex plane, and so, it is argλi=0<π/2m. Since the stability condition is not supplied, the equilibrium point E00,0 is a unstable point for system (16).(b)
From (24), the Jacobian matrix related to E1Λ,0 is(31)JE1Λ,0=−βP−cΛaΛ+10μΛ1+νΛ−δ.
From the equationdetdiagλmα1,λmα2−JE1Λ,0=0, the characteristic equation of eigenvalues is(32)λmα1+βPλmα2−μΛ1+νΛ−δ=0.
Therefore, it isλmα1=−βP and λmα2=μΛ/1+νΛ−δ. These equations are examined as the following:
(1)
λ
m
α
1 is a negative real number due to (17). By De Moivre’s formula, we have λmα1=βPcisπ⇒λi=βP1/mα1cisπ/mα1 fori=1,2,…,mα1, such that cisπ=cosπ+isinπ,i=−1. Considering (9), the stability condition for E1Λ,0 is α1<2 due to argλ=π/mα1>π/2m. This condition has been always provided since 0<α1,α2≤1 in (16).
(2)
On the other hand, we have considered the equationλmα2=μΛ/1+νΛ−δ. If μΛ/1+νΛ<δ, then λi=δ−μΛ/1+νΛ1/mα2cisπ/mα2 for i=1,2,…,mα2 is obtained from λmα2=δ−μΛ/1+νΛcisπ by De Moivre formulas. In this respect, the stability condition is α2<2 by argλ=π/mα2>π/2m. This condition has been always provided from (16). Additionally, If μΛ/1+νΛ≥δ, then the eigenvalues are positive real number due to λi=μΛ/1+νΛ−δ1/mα2cis0 for i=1,2,…,mα2. In this sense, we have argλi=0. The stability condition is not provided due to argλ=0<π/2m. Therefore, the equilibrium point E1Λ,0 is an unstable point for system (16).
Consequently, if the inequality(33)μΛ1+νΛ<δis provided, then the equilibrium point E1Λ,0 is LAS, and if μΛ/1+νΛ−δ≥0, then this point is an unstable point for system (16). When (33) is rearranged, the stability condition of E1Λ,0 is(34)δμ−νδ>Λ.(c)
Jacobian matrix evaluated atE2P∗=δ/μ−νδ,T∗=βP1−P∗/Λ/c/1+aP∗ is(35)JE2=βP1−2P∗Λ+T∗caP∗+1aP∗aP∗+1−1−cP∗aP∗+1T∗μ1+νP∗20,that is,(36)JE2=βPA1A2+1−1−caA2aβPδ2A1cμΛ1−A1A20,where A1 and A2 are defined in (24).
Let us give more details forA1 and A2. Accordingly, it is(37)0<A1<1,since the value P∗, the pathogen size may take, is less than or equal to its carrying capacity Λ. In addition, it is(38)0<A2<1,due to (17), and the components of the equilibrium point in E2 are positive.
From (36), we have the characteristical equation as follows:(39)λmα1+α2−λmα2βPA1A2+1−1+βPδ2A1μΛ1−A1=0.
To be LAS ofE2, it should be that all roots λi for i=1,2,…,mα1+α2 found from the (39) satisfy the inequalities argλi>π/2m in (9).
Proposition is proved.Corollary 1.
Equation (39) can be examined in more detail as shown below. This equation can be rewritten by De Moivre formulas such that(40)λ=rcisθ=cosθ+isinθ,λmα1=rmα1cismα1θ,λmα2=rmα2cismα2θ,where r∈R+, angle θ∈0,2π, and i=−1. By (40), (39) transforms to(41)rmα1+α2cismα1+α2θ−βPA1A2+1−1rmα2cismα2θ+βPδ2A1μΛ1−A1=0,and so,(42)rmα1+α2cosmα1+α2θ−βPA1A2+1−1rmα2cosmα2θ+βPδ2A1μΛ1−A1+irmα1+α2sinmα1+α2θ−βPA1A2+1−1rmα2sinmα2θ=0.
By arranging (42), there is the following system:(43)rmα1+α2sinmα1+α2θ−βPA1A2+1−1rmα2sinmα2θ=0,rmα1+α2cosmα1+α2θ−βPA1A2+1−1rmα2cosmα2θ+βPδ2A1μΛ1−A1=0,and so,(44)sinmα1+α2θ−βPA1A2+1−1r−mα1sinmα2θ=0,cosmα1+α2θ−βPA1A2+1−1r−mα1cosmα2θ+βPδ2A1μΛ1−A1r−mα1+α2=0.
From the first equation in system (44), we have found(45)r=βPA1A2+1−1sinmα2θsinmα1+α2θ1/mα1.
By substituting (45) in the second equation in (44), it is found(46)sinmα2θα2/α1sinmα1θsinmα1+α2θ1+α2/α1=βPδ2A1/μΛ1−A1βPA1A2+1−11+α2/α1.
Consequently, if the anglesθi for i=1,2,…,mα1+α2 obtained from (46) satisfy Routh–Hurwitz stability criteria π/2<θ<π or the condition (9) θ>π/2m, then E2 is LAS.
For equilibria of system (16), the conditions found for LAS and biological existence are summarized in Table 2.Table 2
LAS and biological existence conditions for the equilibria of system (16).
Equilibrium points
Biological existence conditions
LAS conditions
E
0
0,0
Always exists
Unstable point
E
1
Λ
,
0
Always exists
δ
/
μ
−
ν
δ
>
Λ
E
2
P
∗
,
T
∗
δ
/
μ
−
ν
δ
>
Λ
Fori=1,2,…,mα1+α2, either the eigenvalues λi obtained from the characteristic equation λmα1+α2−λmα2βPA1A2+1−1+βPδ2A1/μΛ1−A1=0 satisfy argλ>1/mπ/2 or the angles θi obtained from the equation sinmα2θα2/α1sinmα1θ/sinmα1+α2θ1+α2/α1=βPδ2A1/μΛ1−A1/βPA1A2+1−11+α2/α1 satisfy θi>π/2m
where (P∗,T∗) and (A1,A2) are defined in (19) and (24), respectively.Corollary 2.
E
1 is an unstable point, when E2 exists biologically. Therefore, these equilibrium points cannot be stable when together. Similarly, E2 is biologically meaningless, when E1 is LAS. These circumstances appeared are also seen in Table 2.Corollary 3.
Let us consider the special case ofα1=α2=α for E2. In this case, we have Remark 4. The characteristical equation obtained from DetJE2−λI2=0 is(47)λ2−λβPA1A2+1−1+βPδ2A1μΛ1−A1=0.
The conditions for LAS of the equilibrium pointE2 are either Routh–Hurwitz stability conditions.(48)A1A2+1−1<0,due to (37) and (38), or the conditions(49)A1A2+1−1>0,4βPδ2A1μΛ1−A1>A1A2+1−12,tan−14βPδ2A1/μΛ1−A1−A1A2+1−12−A1A2+1−1>απ2,.
In addition that, by considering equation (46), there may be found a similar condition for stability of this point. In this sense, it is(50)θ=12mαcos−1βPA1A2+1−122δ2A1/μΛ1−A1−1,by formulas of half angle. Thereby, the LAS condition of E2 is(51)cos−1βPA1A2+1−122δ2A1/μΛ1−A1−1>πα,where −1≤βPA1A2+1−12/2δ2A1/μΛ1−A1−1≤1. Corollary 3 is summarized in Table 3.Table 3
The LAS conditions forE2P∗,T∗, in case of α1=α2=α.
LAS conditions ofE2
A
1
A
2
+
1
−
1
<
0 (from Routh–Hurwitz criteria in (14))
Or
A
1
A
2
+
1
−
1
>
0
,
4
β
P
δ
2
A
1
/
μ
Λ
1
−
A
1
>
A
1
A
2
+
1
−
1
2
,
tan
−
1
4
β
P
δ
2
A
1
/
μ
Λ
1
−
A
1
−
A
1
A
2
+
1
−
1
2
/
−
A
1
A
2
+
1
−
1
>
α
π
/
2
, (from argλi>απ/2 in (15))
Or
cos
−
1
β
P
A
1
A
2
+
1
−
1
2
/
2
δ
2
A
1
/
μ
Λ
1
−
A
1
−
1
>
π
α (from (46))


## 4. Numerical Analysis for Model

The proposed model in this study summarizes the general dynamics of the pathogen-immune system. The immune system cells described herein are memory T cells specially produced by the host against the pathogen. Hence, the proposed model is suitable for modeling of diseases such as chronic infections or tumors in which the nonspecific immune system cells of the host at the beginning of the disease have failed to destroy the pathogen.The conditions found in Table2 have been supported by numerical studies shown below. In this section, the time-dependent sizes of the tumor and memory T cells for cancer tried to be estimated by giving the different values to the parameters in the proposed model. The reason for this is to be able to obtain different scenarios and to better demonstrate the results of qualitative analysis. The values of parameters used in system (16) are shown in Table 4.Table 4
The interpretation and considered values of the parameters in the proposed model.
Parameters
Descriptions
Units
Values
For Figure2
For Figure3
For Figure4
β
P
Growth rate of the tumor
Day−1
2.4
2.4
2.4
Λ
Carrying capacity of the tumor
Cells
1
10
5
c
Maximum killing rate of the tumor by immune cells
Day−1
4
4
4
a
Immune cells for half maximum effect on the tumor
Cell−1·day−1
0.2
4
1
μ
The effect of capture rate of immune cells
Day−1
3.98
3.9
3.9
ν
The tumor population size at which the growth rate of immune cells is half its maximum
Cell−1·day−1
1.9
1.9
1.9
δ
Natural death rate of immune cells
Day−1
1.99
1
1
α
1
Fractional-order of the first equation in (16)
A rational number in the interval0,1
0.9
0.8
0.8
α
2
Fractional-order of the second equation in (16)
A rational number in the interval0,1
0.75
0.6
0.6The values calculated by Table4 are given in Table 5.Table 5
The values calculated from Table4 according to Table 2.
Expressions
Terms
Values
For Figure2
For Figure3
For Figure4
Equilibrium pointE1
E
1
Λ
,
0
E
1
1,0
E
1
10,0
E
1
5,0
Stability condition ofE1
δ
/
μ
−
ν
δ
>
Λ
10
>
1 (E1 is LAS)
0.50
<
10 (E1 is unstable)
0.50
<
5 (E1 is unstable)
Equilibrium pointE2
E
2
P
∗
=
δ
/
μ
−
ν
δ
,
T
∗
=
β
P
1
−
P
∗
/
Λ
/
c
/
1
+
a
P
∗
E
2
10
,
−
16.20(biologically meaningless)
E
2
0.50
,
1.71
E
2
0.50
,
0.81
ParameterA1
1
−
P
∗
/
Λ
—
0.95
0.90
ParameterA2
a
P
∗
/
a
P
∗
+
1
—
0.66667
0.3333
Least common multiple of order's denominator
m
—
5
5
Characteristical equation of eigenvalues forE2
λ
m
α
1
+
α
2
−
λ
m
α
2
β
P
A
1
A
2
+
1
−
1
+
β
P
δ
2
A
1
/
μ
Λ
1
−
A
1
=
0
—
λ
7
−
1.40
λ
3
+
1.1692
=
0
λ
7
−
0.48
λ
3
+
1.1077
=
0
The eigenvalues forE2
—
—
λ
1
≈
0.9614
+
0.2454
i
λ
2
≈
0.9614
−
0.2454
i
λ
3
≈
0.1199
+
1.1495
i
λ
4
≈
0.1199
−
1.1495
i
λ
5
≈
−
1.2004
λ
6
≈
−
0.4812
+
0.7135
i
λ
7
≈
−
0.4812
−
0.7135
i
λ
1
≈
0.9287
+
0.3757
i
λ
2
≈
0.9287
−
0.3757
i
λ
3
≈
0.1852
+
1.0410
i
λ
4
≈
0.1852
−
1.0410
i
λ
5
≈
−
1.0798
λ
6
≈
−
0.5741
+
0.7646
i
λ
7
≈
−
0.5741
−
0.7646
i
Angle of eigenvalues forE2
θ
=
arg
λ
—
θ
1
≈
14.3191
∘
,
θ
2
≈
−
14.3191
∘
,
θ
3
≈
84.0452
∘
,
θ
4
≈
−
84.0452
∘
,
θ
5
≈
180
∘
,
θ
6
≈
123.997
∘
,
θ
7
≈
−
123.997
∘
θ
1
≈
22.0255
∘
,
θ
2
≈
−
22.0255
∘
,
θ
3
≈
79.9123
∘
,
θ
4
≈
−
79.9123
∘
,
θ
5
≈
180
∘
,
θ
6
≈
126.901
∘
,
θ
7
≈
−
126.901
∘
Stability condition ofE2
θ
i
>
π
/
2
m
—
E
2 is an unstable point, since θ1,θ2≈14.31∘<18∘.
E
2 is LAS, since θ1,θ2,…,θ7>18∘.
Initial conditions
P
0
,
T
0
0.3
,
0.01
0.3
,
0.01
0.3
,
0.01Through the values in the first columns of Table4, the stability of the equilibrium point E11,0, where the tumor exists and it approaches its carrying capacity, is obtained as shown in Figure 2. This happens within at least 200 days.Figure 2
In case of the first column data in Table4, temporary course of population size of the tumor and memory T cells.When the values in the second columns of Table4 are used, a situation where the system behaviour is limit cycle and E2P∗,T∗ is an unstable point is obtained as shown in Figure 3.Figure 3
In case of the second column data in Table4, temporary course of population size of the tumor and memory T cells.Finally, let us consider the third column in Table4. In here, P∗ for the tumor has a value of 0.5 and T∗ for memory T cells has a value of 0.81. This occurs in at least 100 days as seen in Figure 4. Therefore, these two types of cells stay permanently in the host.Figure 4
In case of the third column data in Table4, temporary course of population sizes of the tumor and memory T cells.


## 5. Conclusions

In this study, the mathematical model examining the changes in the pathogen population size under pressure of specific immune system response in case of cancer or chronic infection has been constructed by the FODE system with multi-orders. According to the results of analysis of model, the pathogen causing disease never disappears in host unless an additional treatment is provided, since the disease is continued by a pathogen, and the free-disease equilibrium pointE00,0 is an unstable point. This case is very suitable for the presumed diseases in the proposed model.When the existence condition ofE2P∗,T∗ and the stability condition of E1Λ,0 in Table 2 is rearranged, then(52)δ<μΛ1+νΛ,δ>μΛ1+νΛ,is obtained, respectively. The term μΛ/1+νΛ is the rate of growth of the memory T cells when the tumor reaches its carrying capacity Λ. The parameter δ is the natural death rate of memory T cells. Accordingly, the value μΛ/1+νΛ−δ can be interpreted as the proliferative power of the memory T cells against the tumor. When this value is negative as seen in Figure 2, the memory T cells cannot reproduced as enough, and consequently, it can be mentioned from the stability of the equilibrium point E1, that only the tumor exists, and the tumor approaches its carrying capacity. Let us consider that this value is positive, and some additional conditions in Table 2 are met. As can be seen in Figure 4, it can be said that the positive equilibrium point E2, in which both tumor and memory T cells exist, is stable.In numerical studies, we tried to estimate the timing and magnitude of the development of the tumor. If treatment procedure for the individual has not been applied, then the results obtained from the proposed model emphasize the fact that either the tumor reaches its maximum size and the memory T cells collapse or the tumor and memory T cells continue to stay together in the host. In the last case mentioned above that the memory T cells of the individual does not collapse, the tumor maintains its presence in the host in a limited manner. The results obtained from analysis are quite consistent with the scenarios of real situations related to the tumor.



---

*Source: 7930603-2018-12-04.xml*