Back
 AM  Vol.11 No.7 , July 2020
Mathematical Model of Classical Kaposi’s Sarcoma
Abstract: In this paper, the global properties of a classical Kaposi’s sarcoma model are investigated. Lyapunov functions are constructed to establish the global asymptotic stability of the virus free and virus (or infection) present steady states. The model considers the interaction of B and progenitor cells in the presence of HHV-8 virus. And how this interaction ultimately culminates in the development of this cancer. We have proved that if the basic reproduction number, R0 is less than unity, the virus free equilibrium point, ε0, is globally asymptotically stable (GAS). We further show that if R0 is greater than unity, then both the immune absent and infection persistent steady states are GAS.

1. Introduction

AIDS-related malignancies such as Kaposi’s sarcoma, non-Hodgkin’s lymphoma, Hodgkin’s disease, primary effusion lymphoma, remain a significant burden for people living with HIV virus. Kaposi’s sarcoma (KS) is the most common neoplasm associated with AIDS.

There are four different forms of known Kaposi sarcoma (KS): Classical (sporadic), African (endemic), acquired immune deficiency syndrome (AIDS)-associated (epidemic), and Transplant or immunosuppression-associated (iatrogenic) KS. The development of each of these forms is dependent on prior infection with Kaposi’s sarcoma-associated Herpesvirus (KSHV) also called Human Herpesvirus-8 (HHV-8). However, KSHV/HHV-8 infection alone is insufficient for the development of KS. Some form of immunodeficiency is also necessary for disease progression.

Several studies [1] [2] have indicated that additional factors contribute to the development of KS in asymptomatic as well as symptomatic KSHV/HHV-8 infected persons. It is known that not every AIDS patient develops KS even in the face of profound immunosuppression, only a minority of KSHV/HHV-8 infected transplant recipients develop iatrogenic KS, and people with classical or endemic KS are not typically immunosuppressed [3] [4]. Most individuals not infected with HIV-1 with strong immune responses have remained latently infected with KSHV/HHV-8 throughout their lifetime [5]. The co-factors involved in classic and endemic KS have yet to be definitely elucidated. Different environmental and genetic factors have been implicated, including age, sex and malnutrition. The progression of KSHV/HHV-8 infection to KS disease depends on a complex and as yet incompletely understood interplay between KSHV and the host immune system that allows for the establishment of a tumor-promoting an environment that influences cellular proliferation, survival, migration, angiogenesis and cytokine/chemokine production. Currently, HIV-1 related KS is the fourth largest killer of people living with HIV/AIDS in sub-Saharan Africa [6].

KS is characterised by abnormal neoangiogenesis, inflammation and proliferation of tumor cells (KS spindle cells—SCs). The review by Foreman et al. [1], has suggested that the infection of progenitor cells by Human Herpesvirus-8 (HHV-8/KSHV) is the pre-requisite for the development of KS.

According to Foreman et al. [1], KS SCs are poorly differentiated endothelial cells (ECs), similar to lymphatic ECs, that are mistakenly infected by the KSHV/HHV-8.

This infection is transmitted either sexually or via saliva and its development is known to be enhanced in individuals with suppressed immune systems. Since HIV-1 is an immunosuppressive virus, it promotes the development of KS in individuals dually infected with both HIV-1 and KSHV/HHV-8.

The HIV-1 growth factors are known to stimulate the immune cells including the healthy and infected B-cells in response to signals from the T cells to proliferate. The response of the actively infected B-cells and the activation of latently infected B-cells to the signals to proliferate leads to the production and increase of KSHV/HHV-8.

The paper is structured as follows: In Section 2, we develop a mathematical model for the pathogenesis of classical Kaposi’s sarcoma based on [1]. We then show that our model is biologically well posed. That is to establish that all model solutions are positive and none of them grows unboundedly. The three model equilibria are found and the basic reproduction number is derived.

In Section 3, we establish the local and global stability of the virus free steady state. The global stability is proved via a Lyapunov function. In Section 4, we establish the local and global stability of the immune absent steady state. This fixed point represents a scenario where the cancer is established and the CD8+ cytotoxic T lymphocytes are dysfunctional due to persistent antigenic stimulation. This phenomenon is prevalent in chronic viral infections [7]. In Section 5, we prove the local and global stability of the infection persistent steady.

In an endeavour to understand the pathogenesis of the KS, we develop an in-host model of KS and find the steady states and establish their local and global stability.

2. Model Formulation

In this study we present a mathematical model of the progression of KS in the absence of HIV-1 (classical KS), by including the interactions of B-cells, HHV-8 infected B cells and progenitor cells, effector cells such as the cytotoxic T lymphocytes (CTLs).

Suppose we first assume that the progenitor cells are as equally prone to HHV-8 infection as the B-cells are. And we further presume that once infected, the progenitor cells proceed directly to the infectious compartment, P i . To simply the model, we ignore the latent infection of these cells.

B-cell and HHV-8 Dynamics

B ˙ = S b β 1 B V 8 μ b B (1)

B ˙ i = β 1 B V 8 μ b i B i β 2 B i E (2)

V ˙ 8 = N μ b i B i β 1 B V 8 β 3 P V 8 μ ν 8 V 8 (3)

Equation (1) describes the dynamics of the susceptible B-cells subjected to HHV-8 infection. The first term represents the constant source for these cells. The second term represents the rate at which the B cells are getting infected by HHV-8. We have assumed a constant infection rate β 1 and the third term represents the natural death rate constant, μ b .

Equation (2) represents a class of productively infected B cells. The first term represents infected B cells that become productively infected at a constant rate β 1 . Most of the infected B cells become latently infected, nevertheless, to lessen the complexity of the model, we have ignored the activation dynamics of latently infected B cells. The second term represents the constant lytic death rate, μ b i , of infected B cells due to bursting of HHV-8. The third term accounts for the lysis of the infected B cells by the HHV-8 specific effector cells, E.

Equation (3) represents a class of the concentration of HHV-8 virions. The first term accounts for the production of HHV-8 virions from the infected B cells where N is regarded as the average number of viral particles produced by the infected B cell during its lifespan. The second and third terms account for the loss of HHV-8 virions as they infect the B and progenitor cells and the fourth term denotes the constant clearance rate, μ V 8 .

Progenitor and KS cells Dynamics

P ˙ = S p μ p P β 3 P V 8 (4)

P ˙ i = β 3 P V 8 μ p i P i β 4 P i E (5)

K ˙ = μ ^ p i P i μ k K (6)

Equations (4) to (6) describe the dynamics of the progenitor cells and KS cells. Upon activation by the cytokines secreted by activated immune cells, progenitor cells predispose themselves to HHV-8 infection [for details see Foreman, et al.]. Then, through autocrine and paracrine mechanisms, progenitor cells progress to KS. Equation (4) represents a class of susceptible progenitor cells. The first term denotes the constant source of these precursor cells from the bone marrow. The second term accounts for the natural death rate of progenitor cells and the third term accounts for the rate at which the progenitor cells get infected by the HHV-8.

Equation (5) represents the class of productively infected progenitor cells. The first term accounts for the progenitor cells that get infected by HHV-8. The second term denotes the constant nature death rate of progenitor cells. The third term accounts for the loss of infected progenitor cells due to killing by effector cells, E.

Equation (6) represents the class of KS cells. The first term accounts for the gain from the infected progenitor cells as they transform into KS cells and the second term denotes a loss due to natural death. This natural death rate is quite negligible due to the fact that these cells produce cytokines that interfere with the apoptosis program.

T-cell Dynamics

E ˙ = r B i E μ e E (7)

Equation (7) describes the dynamics of the HHV-8 specific effector cells that kill infected B and progenitor cells. The first term represents the proliferation of these cells due to infected B cells’ stimulation at a constant rate r and the second term is the constant natural death rate.

The above dynamics are illustrated in the following schematic diagram represented in Figure 1. We summarize the parameter definitions and their values in Table 1 and Table 2, respectively.

2.1. Positivity of Solutions and Boundedness of Solutions

2.1.1. Positivity of Solutions

For the model system above to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all times. In other words, solutions of the model above with positive initial data remain positive for all t > 0 .

We thus show that our model is well posed in the sense that the populations do not become negative, and the populations are bounded. Also, the non-negative orthant is positively invariant, in other words, any trajectory that begins in the non-negative orthant remains there for all time.

For simplicity of notations, let

Figure 1. Schematic diagram of classical KS dynamics.

Table 1. Model parameters and their definitions.

Table 2. Estimation of parameters.

x 1 ( t ) : = B ( t ) , x 2 ( t ) : = P ( t ) , x 3 ( t ) : = E ( t ) , x 4 ( t ) : = B i ( t ) , x 5 ( t ) : = P i ( t ) , x 6 ( t ) : = V 8 ( t ) , x 7 ( t ) : = K ( t ) S 1 : = S b , S 2 : = S p , μ 1 : = μ b , μ 2 : = μ p , μ 3 : = μ e , μ 4 : = μ b i , μ 5 : = μ p i , μ 6 : = μ ν 8 , μ ^ 5 : = μ ^ p i , μ 7 : = μ k

Denote by + 7 the set of points

x t = ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) 7 (8)

with positive coordinates and consider the system with initial values

x 0 = ( x 0 1 , x 0 2 , x 0 3 , x 0 4 , x 0 5 , x 0 6 , x 0 7 ) + 7 (9)

Lemma 1 Let x 0 i 0 , i = 1 , 2 , 3 , , 7 . Then the solution x t + 7 for all t > 0 in the region

Γ 7 = { ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) | x i i , i = 1 , 2 , 3 , , 7 } 7 .

The systems (1)-(7) can be rewritten in terms of two sub-systems as follows:

Susceptible States

x ˙ 1 = S 1 μ 1 x 1 β 1 x 1 x 6 (10)

x ˙ 2 = S 2 μ 2 x 2 β 3 x 2 x 6 (11)

x ˙ 3 = r x 3 x 4 μ 3 x 3 (12)

Infected States

x ˙ 4 = β 1 x 1 x 6 μ 4 x 4 β 2 x 3 x 4 (13)

x ˙ 5 = β 3 x 2 x 6 μ 5 x 5 β 4 x 3 x 5 (14)

x ˙ 6 = N μ 4 x 4 β 1 x 1 x 6 β 3 x 2 x 6 μ 6 x 6 (15)

x ˙ 7 = μ ^ 5 x 5 μ 7 x 7 (16)

Proof. The system with Equations (10)-(12) can be written as a system of differential inequalities

d x i d t ( A i j = 1 7 B i j x j ) x i + S i (17)

where B i j 0 , S = ( S 1 S 2 S 3 ( t ) ) T > ( 0 0 0 ) T .

Suppose the assertion x i ( t ) 0 , i = 1 , 2 , 3 is not true. Then there exists a smallest number t 0 , such that

x i ( t ) > 0 f o r 1 i 3,0 t t 0

x i ( t 0 ) = 0 f o r a t l e a s t o n e i , s a y i 0 .

Then, x i 0 is a decreasing function and

d x i 0 ( t 0 ) d t 0.

From the differential inequality (17) for x i 0 ( t ) we get

d x i 0 ( t 0 ) d t π i > 0

which is a contradiction.

Hence, if x i ( 0 ) 0 , i = 1 , 2 , 3 then, x i ( t ) 0 for all t > 0 , i = 1 , 2 , 3.

The system with equations (13)-(16) representing the infected states can be written in matrix form as

Y ˙ ( t ) = M Y (18)

where

Y = [ x 4 x 5 x 6 x 7 ] T (19)

M = [ ( μ 4 + β 2 x 3 ) 0 β 1 x 1 0 0 ( μ 5 + β 4 x 3 ) β 3 x 2 0 N μ 4 0 ( μ 6 + β 1 x 1 + β 3 x 2 ) 0 0 μ ^ 5 0 μ 7 ] (20)

M is a Mertzler matrix. Hence, the infected states x i ( t ) 0 , i = 4 , 5 , 6 , 7 .

2.1.2. Boundedness of Solutions

We already have a lower bound given by the above lemma 1 for the solutions of model (10)-(16) with nonnegative initial values. We now show that the solutions are bounded from above.

We denote by C b ( I ) the set of all continuous and bounded functions defined on the interval I taking values in 7 .

Lemma 2 Let φ : [ 0, ) 7 be a solution of system (10)-(12). If φ ( 0 ) + 7 , then φ C b [ 0, ) .

Proof. Because of Lemma 1, it only remains to prove the existence of an upper bound to the nonnegative solutions of system (10)-(16). Let y 1 ( t ) = x 1 ( t ) + x 4 ( t ) be the total concentration of the B cells. Then,

y ˙ 1 ( t ) = x ˙ 1 ( t ) + x ˙ 4 ( t ) = S 1 β 1 x 1 x 6 μ 1 x 1 + β 1 x 1 x 6 μ 4 x 4 β 2 x 3 x 4 < S 1 μ 1 x 1 μ 4 x 4 S 1 d 1 y 1 ( t ) , w h e r e d 1 = min { μ 1 , μ 4 } y 1 ( t ) y 1 ( 0 ) + S 1 d 1 = M 1 . (21)

Let y 2 ( t ) = x 2 ( t ) + x 5 ( t ) be the total concentration of the progenitor cells. Then,

y ˙ 2 ( t ) = x ˙ 2 ( t ) + x ˙ 5 ( t ) = S 2 β 3 x 2 x 6 μ 2 x 2 + β 3 x 2 x 6 μ 5 x 5 β 4 x 3 x 5 S 2 μ 2 x 2 μ 5 x 5 S 2 d 2 y 2 ( t ) , w h e r e d 2 = min { μ 2 , μ 5 } y 2 ( t ) y 2 ( 0 ) + S 2 d 2 = M 2 . (22)

For the KS cells,

x ˙ 7 = μ ^ 5 x 5 ( t ) μ 7 x 7 μ ^ 5 M 2 μ 7 x 7 x 7 ( t ) x 7 ( 0 ) + μ ^ 5 M 2 μ 7 = M 3 . (23)

We now consider the HHV-8 dynamics equation,

x ˙ 6 = N μ 4 x 4 β 1 x 1 x 6 β 3 x 2 x 6 μ 6 x 6 N μ 4 M 1 μ 6 x 6 x 6 ( t ) M 4 , w h e r e M 4 = x 6 ( 0 ) + N μ 4 M 1 μ 6 . (24)

For the HHV-8 specific effector cells, x 3 , we argue as follows:

x ˙ 3 ( t ) ( r M 1 μ 3 ) x 3 , μ 3 < r M 1 , 0 < x 3 ( t ) < . (25)

2.2. Equilibria

To find the steady states of the model system (10)-(16), we set each differential equation to solve and solve the resulting system simultaneously from which we obtain three equilibria. These are virus free equilibrium, ε 0 , immune absent equilibrium, ε * and the endemic equilibrium ε * * in which are the coordinates are strictly positive.

ε 0 = ( x 1 0 , x 2 0 , x 3 0 , x 4 0 , x 5 0 , x 6 0 , x 7 0 ) , w h e r e x 1 0 = S 1 μ 1 , x 2 0 = S 2 μ 2 , x i 0 = 0 , i = 3 , 4 , 5 , 6 , 7. (26)

2.2.1. Derivation of the Basic Reproduction Number, R 0

The method of Castillo-Chavez et al. [13] is used to compute the reproductive number for the model (10)-(16) and representing the total number of newly infected B cells arising out of one infected B cell.

Heterogeneity is defined using groups defined by fixed characteristics:

( d X d t = f ( X , Y , Z ) , d Y d t = g ( X , Y , Z ) , d Z d t = h ( X , Y , Z ) . (27)

where X 3 , Y 3 , Z , and h ( X ,0,0 ) = 0 . Assuming that g ( X * , Y , Z ) = 0 implicitly determines a function Y = g ˜ ( X * , Z ) . Let A = D Z h ( X * , g ˜ ( X * ,0 ) ,0 ) and further assume that A can be written in the form A = M D , with M 0 (that is m i j 0 ) and D > 0 , a diagonal matrix. The reproduction numbers are then evaluated from the matrix M D 1 .

The cell population subgroups are divided between X uninfected cells, Y infected cells and Z HHV-8. Then, X = ( x 1 , x 2 , x 3 ) , Y = ( x 3 , x 4 , x 7 ) and Z = x 6 .

Let U 0 = ( X * ,0,0 ) denote the virus free equilibrium, that is, f ( X * , 0 , 0 ) = g ( X * , 0 , 0 ) = h ( X * , 0 , 0 ) = 0 and Y = g ˜ ( X * , Z ) , where

g ˜ ( X * , Z ) = ( β 1 S 1 x 6 μ 1 μ 4 μ ^ 5 x 5 μ 7 ) (28)

Computing A = D Z h ( X * , g ˜ ( X * ,0 ) ,0 ) gives

( M = N β 1 S 1 μ 1 , D = μ 1 μ 2 μ 6 + β 1 S 1 μ 2 + β 3 S 2 μ 1 μ 1 μ 2 , (29)

Therefore,

R 0 = N μ 2 β 1 S 1 μ 1 μ 2 μ 6 + β 1 S 1 μ 2 + β 3 S 2 μ 1 . (30)

Let us now determine respectively, the immune absent and infection persistent equilibria

1) ε * = ( x 1 * , x 2 * , x 3 * , x 4 * , x 5 * , x 6 * , x 7 * ) ,

2) ε * * = ( x 1 * * , x 2 * * , x 3 * * , x 4 * * , x 5 * * , x 6 * * , x 7 * * ) ,

where

{ x 1 * = S 1 μ 1 + β 1 x 6 * , x 2 * = S 2 μ 2 + β 3 x 6 * , x 3 * = 0 , x 4 * = β 1 S 1 x 6 * μ 4 ( μ 1 + β 1 x 6 * ) , x 5 * = β 3 S 2 x 6 * μ 5 ( μ 2 + β 3 x 6 * ) , x 7 * = μ ^ 5 x 5 * μ 7 . (31)

where x 6 * is a positive solution to the following quadratic equation

Q 1 ( x 6 ) = A 1 x 6 2 + A 2 x 6 + A 3 = 0 ,

where

A 1 = β 1 β 3 μ 6 < 0 , A 2 = N β 1 β 3 S 1 β 1 β 3 S 2 μ 2 μ 6 β 1 μ 1 μ 6 β 3 β 1 β 3 S 1 , A 3 = ( β 1 S 1 μ 2 + β 3 S 2 μ 1 + μ 1 μ 2 μ 6 ) ( R 0 1 ) . (32)

{ x 1 * * = S 1 μ 1 + β 1 x 6 * * , x 2 * * = S 2 μ 2 + β 3 x 6 * * , x 3 * * = β 1 ( r S 1 μ 3 μ 4 ) x 6 * * μ 1 μ 3 μ 4 μ 3 β 2 ( μ 1 + β 1 x 6 * * ) , x 4 * * = μ 3 r , x 5 * * = β 3 S 2 x 6 * * ( μ 5 + β 4 x 3 * * ) ( μ 2 + β 3 x 6 * * ) . x 7 * * = μ ^ 5 x 5 * * μ 7 . (33)

where x 6 * * is a positive real solution to the following cubic equation

Q 2 ( x 6 ) = B 1 x 6 3 + B 2 x 6 2 + B 3 x 6 + B 4 = 0 ,

where

B 1 = β 1 β 3 r μ 6 < 0 , B 2 = N β 1 β 3 μ 3 μ 4 r β 1 β 3 S 2 r β 1 r μ 2 μ 6 β 2 r μ 1 μ 6 β 1 β 3 S 1 r , B 3 = N β 1 μ 2 μ 3 μ 4 r β 3 S 2 r μ 1 r μ 1 μ 2 μ 6 β 1 S 1 r μ 2 + N β 3 μ 1 μ 3 μ 4 r , B 4 = N μ 1 μ 2 μ 3 μ 4 r > 0. (34)

We need to establish the existence of a positive value of x 6 * * .

Note:

Q 2 ( 0 ) = B 4 > 0 , a n d l i m x 6 Q 2 ( x 6 ) = ,

Due to continuity of Q 2 , it follows that there exists a positive value of x 6 , say, x 6 * * ( 0, ) , such that Q 2 ( x 6 * * ) = 0 .

Note also that the existence of x 3 * * requires that

x 6 * * > μ 1 μ 3 μ 4 β 1 ( r S 1 μ 3 μ 4 ) , a n d i f x 6 * * = μ 1 μ 3 μ 4 β 1 ( r S 1 μ 3 μ 4 ) , t h e n x 3 * * = 0 (35)

and we obtain the immune absent equilibrium point, ε * .

Lemma 3 Consider the system (10)-(16). Then immune absent equilibrium point, ε * , exists if R 0 > 1 , and unstable otherwise.

Proof. Note that

Q 1 ( 0 ) = A 3 > 0 , s i n c e R 0 > 1 a n d l i m x 6 * Q 1 ( x 6 * ) = .

Hence, there exists a positive value x 6 * ( 0, + ) , such that Q 1 ( x 6 * ) = 0 . In particular,

x 6 * = A 2 A 2 2 4 A 1 A 3 2 A 1 > 0.

However, if R 0 = 1 , then x 6 * = 0 is a root and we get the virus free equilibrium, ε 0 .

And if R 0 < 1 , then Q 1 ( x 6 * ) < 0 , x 6 * > 0 .

3. Stability Analysis of Virus Free Equilibrium Point, ε0

3.1. Local Stability of Virus Free Equilibrium Point, ε0

Theorem 1 For system (10)-(16), if R 0 < 1 , then the virus free equilibrium, ε 0 , is locally asymptotically stable and unstable otherwise.

Proof. Consider the Jacobian matrix, J ( ε 0 ) , evaluated at ε 0 , of the system (10)-(16).

J ( ε 0 ) = ( μ 1 0 0 0 0 β 1 x 1 0 0 0 μ 2 0 0 0 β 3 x 2 0 0 0 0 μ 3 0 0 0 0 0 0 0 μ 4 0 β 1 x 1 0 0 0 0 0 0 μ 5 β 3 x 2 0 0 0 0 0 N μ 4 0 ( μ 6 + β 1 x 1 0 + β 3 x 2 0 ) 0 0 0 0 0 μ ^ 5 0 μ 7 ) (36)

The corresponding characteristic equation has the following eigenvalues

λ 1 = μ 1 , λ 2 = μ 2 , λ 3 = μ 3 , λ 4 = μ 7

and the other three eigenvalues are obtained from the following reduced matrix, J 1

J 1 = ( μ 1 0 β 1 x 1 0 0 μ 5 β 3 x 2 0 N μ 4 0 ( μ 6 + β 1 x 1 0 + β 3 x 2 0 ) ) (37)

From (37), we obtain the characteristic equation and investigate the nature of the roots (or eigenvalues).

f ( λ ) = λ 3 + B 1 λ 2 + B 2 λ + B 3 = 0 , (38)

where

B 1 = μ 4 + μ 5 + μ 6 + β 1 S 1 μ 1 + β 3 S 2 μ 2 > 0 ,

B 2 = μ 5 ( μ 4 + μ 6 + β 1 S 1 μ 1 + β 3 S 2 μ 2 ) + μ 4 ( μ 6 + β 1 S 1 μ 1 + β 3 S 2 μ 2 ) [ 1 R 0 ] > 0 , i f R 0 < 1.

B 3 = μ 4 μ 5 μ 6 + β 1 μ 4 μ 5 x 1 0 + β 3 μ 4 μ 5 x 2 0 N β 1 μ 4 μ 5 x 1 0 = μ 4 μ 5 N β 1 μ 2 S 2 R 0 2 [ 1 R 0 1 ] > 0 , s i n c e R 0 < 1.

B 1 B 2 B 3 = ( μ 4 + μ 5 + μ 6 + β 1 x 1 0 + β 3 x 2 0 ) ( μ 4 μ 5 + μ 4 μ 6 + μ 5 μ 6 + β 1 μ 4 x 1 0 + β 1 μ 5 x 1 0 + β 3 μ 4 x 2 0 + β 3 μ 5 x 2 0 N β 1 μ 4 x 1 0 ) μ 4 μ 5 μ 6 β 1 μ 4 μ 5 x 1 0 β 3 μ 4 μ 5 x 2 0 + N β 1 μ 4 μ 5 x 1 0 = ( μ 4 + μ 6 + β 1 S 1 μ 1 + β 3 S 2 μ 2 ) ( μ 4 μ 5 + μ 5 μ 6 + μ 5 2 + β 1 S 1 μ 5 μ 1 + β 3 S 2 μ 5 μ 2 + μ 4 ( μ 6 + β 1 S 1 μ 1 + β 3 S 2 μ 2 ) [ 1 R 0 ] ) > 0. (39)

Based on Hurwtz criterion, all roots of (38) have the negative real parts. Thus, if R 0 < 1 , the virus free equilibrium ε 0 , is locally asymptotically stable.

If

R 0 > 1 , l e t f ( λ ) = λ 3 + B 1 λ 2 + B 2 λ + B 3 , (40)

then it is obvious that

f ( 0 ) = B 3 < 0 , l i m λ + f ( λ ) = . (41)

Hence, there must exist a λ 0 > 0 such that f ( λ 0 ) = 0 . This shows that equation (38) has at least one root with positive real part. Thus, if R 0 > 1 , the virus free equilibrium ε 0 is unstable.

3.2. Global Stability of the Virus Free Equilibrium, ε0

Theorem 2 For the system (10)-(16), if R 0 β 1 S 1 μ 2 β 1 S 1 μ 2 + β 3 S 2 μ 1 < 1 , then ε 0 is globally asymptomatically stable (GAS).

Proof. Define the Lyapunov function

V 0 = x 1 0 ( x 1 x 1 0 1 ln x 1 x 1 0 ) + x 2 0 ( x 2 x 2 0 1 ln x 2 x 2 0 ) + N β 2 r ( N 1 ) x 3 + N N 1 x 4 + N N 1 x 5 + 1 N 1 x 6 + N μ 5 μ ^ 5 ( N 1 ) x 7 , (42)

Differentiating along the solutions of system (10)-(16), we get

V ˙ 0 = ( 1 x 1 0 x 1 ) x ˙ 1 + ( 1 x 2 0 x 2 ) x ˙ 2 + N β 1 r ( N 1 ) x ˙ 3 + N N 1 x ˙ 4 + N N 1 x ˙ 5 + 1 N 1 x ˙ 6 + N μ 5 μ ^ 5 ( N 1 ) x ˙ 7 = ( 1 x 1 0 x 1 ) [ S 1 β 1 x 1 x 6 μ 1 x 1 ] + ( 1 x 2 0 x 2 ) [ S 2 β 3 x 2 x 6 μ 2 x 2 ] + N β 2 r ( N 1 ) [ r x 3 x 4 μ 3 x 3 ] + N N 1 [ β 1 x 1 x 6 β 2 x 3 x 4 μ 4 x 4 ] + N N 1 [ β 3 x 2 x 6 β 4 x 3 x 5 μ 5 x 5 ] + 1 N 1 [ N μ 4 x 4 β 1 x 1 x 6 β 3 x 2 x 6 μ 6 x 6 ] + N μ 5 μ ^ 5 ( N 1 ) [ μ ^ 5 x 5 μ 7 x 7 ] . (43)

Applying S 1 = μ 1 x 1 0 and S 2 = μ 2 x 2 0 , we obtain

V ˙ 0 = μ 1 x 1 0 ( 2 x 1 0 x 1 x 1 x 1 0 ) + μ 2 x 2 0 ( 2 x 2 0 x 2 x 2 x 2 0 ) N β 2 μ 3 r ( N 1 ) x 3 N N 1 β 4 x 3 x 5 N μ 5 μ 7 μ ^ 5 ( N 1 ) x 7 + β 1 S 1 μ 2 + β 3 S 2 μ 1 + μ 1 μ 2 μ 6 μ 1 μ 2 ( N 1 ) [ ( 1 + β 3 S 2 μ 1 β 1 S 1 μ 2 ) R 0 1 ] x 6 0. (44)

Since the arithmetic mean is greater than or equal to the geometric mean, 2 x i 0 x i x i x i 0 0 , i = 1 , 2 And the global stability follows from the LaSalle’s invariance principle [14].

4. Stability Analysis of Immune Absent Equilibrium Point, ε*

4.1. Local Stability of Immune Absent Equilibrium Point, ε*

Theorem 3 For system (10)-(16), if x 1 * μ 6 β 1 ( N 1 ) and x 4 * < μ 3 r , then the immune absent equilibrium point, ε * , is locally asymptotically stable and unstable otherwise.

Proof. Consider the Jacobian matrix, J ( ε * ) , evaluated at ε * , of the system (10)-(16).

J ( ε * ) = ( ( μ 1 + β 1 x 6 * ) 0 0 0 0 β 1 x 1 * 0 0 ( μ 2 + β 3 x 6 * ) 0 0 0 β 3 x 2 * 0 0 0 ( μ 3 r x 4 * ) 0 0 0 0 β 1 x 6 * 0 β 1 x 4 * μ 4 0 β 1 x 1 * 0 0 β 3 x 6 * β 4 x 5 * 0 μ 5 β 3 x 2 * 0 β 1 x 6 * β 3 x 6 * 0 N μ 4 0 ( μ 6 + β 1 x 1 * + β 3 x 2 * ) 0 0 0 0 0 μ ^ 5 0 μ 7 ) (45)

The corresponding characteristic equation has the following eigenvalues

λ 1 = μ 5 , λ 2 = ( μ 3 r x 4 * ) , λ 3 = μ 7

and the other four eigenvalues are obtained from the following reduced matrix, J 2

J 2 = ( ( μ 1 + β 1 x 6 * ) 0 0 β 1 x 1 * 0 ( μ 2 + β 3 x 6 * ) 0 β 3 x 2 * β 1 x 6 * 0 μ 4 β 1 x 1 * β 1 x 6 * β 3 x 6 * N μ 4 ( μ 6 + β 1 x 1 * + β 3 x 2 * ) ) (46)

From (45), we obtain the characteristic equation and investigate the nature of the roots (or eigenvalues).

G ( λ ) = λ 4 + C 1 λ 3 + C 2 λ 2 + C 3 λ + C 4 = 0 , (47)

where

C 1 = μ 1 + μ 2 + μ 4 + μ 6 + β 1 x 1 * + β 3 x 2 * + β 1 x 6 * + β 3 x 6 * > 0 ,

C 2 = μ 1 μ 2 + μ 1 μ 4 + μ 2 μ 4 + μ 1 μ 6 + μ 2 μ 6 + μ 4 μ 6 + β 1 μ 1 x 1 * + β 1 μ 2 x 1 * + β 1 μ 4 x 1 * + β 3 μ 1 x 2 * + β 3 μ 2 x 2 * + β 1 μ 2 x 6 * + β 3 μ 4 x 2 * + β 3 μ 1 x 6 * + β 1 μ 4 x 6 * + β 1 μ 6 x 6 * + β 3 μ 4 x 6 * + β 3 μ 6 x 6 * + β 1 β 3 x 6 * 2 + β 1 β 3 x 1 * x 6 * + β 1 β 3 x 2 * x 6 * N β 1 μ 4 x 1 * = μ 1 μ 2 + μ 1 μ 4 + μ 2 μ 4 + μ 1 μ 6 + μ 2 μ 6 + β 1 μ 1 x 1 * + β 1 μ 2 x 1 * + β 3 μ 1 x 2 * + β 3 μ 2 x 2 * + β 1 μ 2 x 6 * + β 3 μ 1 x 6 * + β 1 μ 4 x 6 * + β 1 μ 6 x 6 * + β 1 μ 4 x 6 * + β 3 μ 6 x 6 * + β 1 β 3 x 6 * 2 + β 1 β 3 x 1 * x 6 8 + β 1 β 3 x 2 * x 6 * > 0.

C 3 = μ 1 μ 2 μ 4 + μ 1 μ 2 μ 6 + μ 1 μ 4 μ 6 + μ 2 μ 4 μ 6 + β 1 μ 1 μ 2 x 1 * + β 1 μ 1 μ 4 x 1 * + β 1 μ 2 μ 4 x 1 * + β 3 μ 1 μ 2 x 2 * + β 3 μ 1 μ 4 x 2 * + β 3 μ 2 μ 4 x 2 * + β 1 μ 2 μ 4 x 6 * + β 3 μ 1 μ 4 x 6 * + β 1 μ 2 μ 6 x 6 * + β 3 μ 1 μ 6 x 6 * + β 1 μ 4 μ 6 x 6 * + β 3 μ 4 μ 6 x 6 * + β 1 β 3 μ 4 x 6 * 2 + β 1 β 3 μ 6 x 6 * 2 N β 1 μ 1 μ 4 x 1 * N β 1 μ 2 μ 4 x 1 * + β 1 β 3 μ 1 x 1 * x 6 * + β 1 β 3 μ 2 x 2 * x 6 * + β 1 β 3 μ 4 x 1 * x 6 * + β 1 β 3 μ 4 x 2 * x 6 * N β 1 β 3 μ 4 x 1 * x 6 * = μ 1 μ 2 μ 4 + μ 1 μ 2 μ 6 + β 1 μ 1 μ 2 x 1 * + β 3 μ 1 μ 2 x 2 * + β 1 μ 2 μ 4 x 6 * + β 1 μ 1 μ 4 x 6 * + β 1 μ 2 μ 6 x 6 * + β 3 μ 1 μ 6 x 6 * + β 1 μ 4 μ 6 x 6 * + β 1 β 3 μ 4 x 6 * 2 + β 1 β 3 μ 6 x 6 * 2 + β 1 β 3 μ 1 x 1 * x 6 * + β 1 β 3 μ 2 x 2 * x 6 * > 0.

C 4 = μ 1 μ 2 μ 4 μ 6 + β 1 μ 1 μ 2 μ 4 x 1 * + β 3 μ 1 μ 2 μ 4 x 2 * + β 1 μ 2 μ 4 μ 6 x 6 * + β 3 μ 1 μ 4 μ 6 x 6 * + β 1 β 3 μ 4 μ 6 x 6 * 2 N β 1 μ 1 μ 2 μ 4 x 1 * + β 1 β 3 μ 1 μ 4 x 1 * x 6 * + β 1 β 3 μ 2 μ 4 x 2 * x 6 * N β 1 β 3 μ 1 μ 4 x 1 * x 6 * = β 3 μ 1 μ 4 x 6 * ( μ 6 + β 1 x 1 * N β 1 x 1 * ) + β 1 β 3 μ 2 μ 4 x 2 * x 6 * + β 1 μ 2 μ 4 μ 6 x 6 * + β 1 β 3 μ 4 μ 6 x 6 * 2 > 0. (48)

Based on Hurwitz criterion, all roots of (47) have the negative real parts. Thus, the immune absent equilibrium ε * , is locally asymptotically stable.

x 1 * > μ 6 β 1 ( N 1 ) , l e t G ( λ ) = λ 4 + C 1 λ 3 + C 2 λ 2 + C 3 λ + C 4 , (49)

then it is obvious that

G ( 0 ) = B 3 < 0 , l i m λ + P ( λ ) = . (50)

Hence, there must exist a λ * > 0 such that G ( λ * ) = 0 . This shows that equation (47) has at least one root with positive real part.

x 1 * > μ 6 β 1 ( N 1 ) , (51)

the immune absent equilibrium point ε * is unstable.

4.2. Global Stability of the Immune Absent Equilibrium Point, ε*

In the ensuing proof we neglect the loss of HHV-8 virions due to the infection process of both the healthy B and progenitor cells.

Theorem 4 For the system (10)-(16), if R 0 > 1 , then the immune absent steady state, ε * , is GAS.

Proof. Consider the Lyapunov function, V 1 , defined by

V 1 = ( x 1 x 1 * x 1 * l n x 1 x 1 * ) + ( x 4 x 4 * x 4 * l n x 4 x 4 * ) + β 1 x 1 * x 6 * N μ 4 x 4 * ( x 6 x 6 * x 6 * l n x 6 x 6 * ) (52)

Differentiating along the solutions of the system (10)-(16), we obtain

V 1 = ( 1 x 1 * x 1 ) x ˙ 1 + ( 1 x 4 * x 4 ) x ˙ 4 + β 1 x 1 * x 6 * N μ 4 x 4 * ( 1 x 6 * x 6 ) x ˙ 6 (53)

Using the relations at equilibrium point and relabeling y i = x i x i * , i = 1 , 2 , , 7 , we get

V 1 = μ 1 x 1 * ( 2 y 1 1 y 1 ) + β 1 x 1 * x 6 * ( 1 y 1 y 6 1 y 1 + y 6 ) + β 1 x 1 * x 6 * ( y 1 y 6 y 4 y 1 y 6 y 4 + 1 ) + β 1 x 1 * x 6 * ( y 4 y 6 y 4 y 6 + 1 ) = μ 1 x 1 * ( 2 y 1 1 y 1 ) + β 1 x 1 * x 6 * ( 3 1 y 1 y 4 y 6 y 1 y 6 y 4 ) < 0. (54)

Since the arithmetic mean is greater than or equal to the geometric mean, 2 y 1 1 y 1 0 and 3 1 y 1 y 4 y 6 y 1 y 6 y 4 0 . And the global stability follows from the LaSalle’s invariance principle [14].

5. Stability Analysis of Infection Persistent Equilibrium Point, ε**

5.1. Local Stability of Infection Persistent Equilibrium Point, ε**

Theorem 5 For system (10)-(16), if x 4 * * μ 3 r and x 1 * * x 3 * * μ 6 N μ 4 , then the infection persistent equilibrium point, ε * * , is locally asymptotically stable and unstable otherwise.

Proof. Consider the Jacobian matrix, J ( ε * * ) , evaluated at ε * * , of the system (10)-(16).

J ( ε * * ) = ( a 11 0 0 0 0 β 1 x 1 * * 0 0 a 22 0 0 0 β 3 x 2 * * 0 0 0 a 33 r x 3 * * 0 0 0 β 1 x 6 * * 0 β 2 x 4 * * a 44 0 β 1 x 1 * * 0 0 β 3 x 6 * * β 4 x 5 * * 0 a 55 β 3 x 2 * * 0 β 1 x 6 * * β 3 x 6 * * 0 N μ 4 0 a 66 0 0 0 0 0 μ ^ 5 0 μ 7 ) (55)

where

a 11 = ( μ 1 + β 1 x 6 * * ) , a 22 = ( μ 2 + β 3 x 6 * * ) , a 33 = ( r x 4 * * μ 3 ) , a 44 = ( μ 4 + β 2 x 3 * * ) , a 55 = ( μ 5 + β 4 x 3 * * ) , a 66 = ( μ 6 + β 1 x 1 * * + β 3 x 2 * * ) (56)

The corresponding characteristic equation has the following eigenvalues

λ 1 = μ 7 , λ 2 = ( μ 5 + β 3 x 3 * * )

and the other five eigenvalues are obtained from the following reduced matrix, J 3

J 3 = ( b 11 0 0 0 β 1 x 1 * * 0 b 22 0 0 β 3 x 2 * * 0 0 b 33 r x 3 * * 0 β 1 x 6 * * 0 β 2 x 4 * * b 44 β 1 x 1 * * β 1 x 6 * * β 3 x 6 * * 0 N μ 4 b 55 ) (57)

where

b 11 = ( μ 1 + β 1 x 6 * * ) , b 22 = ( μ 2 + β 3 x 6 * * ) , b 33 = ( r x 4 * * μ 3 ) , b 44 = ( μ 4 + β 2 x 3 * * ) , b 55 = ( μ 6 + β 1 x 1 * * + β 3 x 2 * * ) (58)

From (55), we obtain the following characteristic equation and investigate the nature of the roots (or eigenvalues).

H ( λ ) = λ 5 + D 1 λ 4 + D 2 λ 3 + D 3 λ 2 + D 4 λ + D 5 = 0 , (59)

where

D 1 = μ 1 + μ 2 + μ 3 + μ 4 + μ 5 + β 1 x 1 * * + β 1 x 3 * * + β 3 x 2 * * + β 1 x 6 * * + β 3 x 6 * * r x 4 * * > 0 ,

D 2 = μ 1 μ 2 + μ 1 μ 3 + μ 1 μ 4 + μ 2 μ 3 + μ 2 μ 4 + μ 1 μ 6 + μ 3 μ 4 + μ 2 μ 6 + μ 3 μ 6 + μ 4 μ 6 + β 1 2 x 1 * * x 3 * * + β 1 2 x 3 * * x 6 * * + β 1 μ 1 x 1 * * + β 1 μ 2 x 1 * * + β 1 μ 1 x 3 * * + β 1 μ 3 x 1 * * + β 1 μ 2 x 3 * * + β 1 μ 4 x 1 * * + β 3 μ 1 x 2 * * + β 1 μ 3 x 3 * * + β 3 μ 2 x 2 * * + β 3 μ 3 x 2 * * + β 1 μ 2 x 6 * * + β 3 μ 4 x 2 * * + β 1 μ 3 x 6 * * + β 1 μ 6 x 3 * * + β 3 μ 1 x 6 * * + β 1 μ 4 x 6 * * + β 3 μ 3 x 6 * * + β 1 μ 6 x 6 * * + β 3 μ 4 x 6 * * + β 3 μ 4 x 6 * * + β 3 μ 6 x 6 * * r μ 1 x 4 * * r μ 2 x 4 * * r μ 4 x 4 * * r μ 6 x 4 * * + β 1 β 3 x 6 * * 2 + β 1 β 3 x 2 * * x 3 * * + β 1 β 3 x 1 * * x 6 * * + β 1 β 3 x 2 * * x 6 * * + β 1 β 3 x 3 * * x 6 * * β 1 r x 1 * * x 4 * * β 3 r x 2 * * x 4 * * β 1 r x 4 * * x 6 * * β 3 r x 4 * * x 6 * * N β 1 μ 4 x 1 * *

= μ 1 μ 2 + μ 1 μ 4 + μ 2 μ 4 + μ 1 μ 6 + μ 2 μ 6 + μ 4 μ 6 + β 1 2 x 1 * * x 3 * * + β 1 2 x 3 * * x 6 * * + β 1 μ 1 x 1 * * + β 1 μ 2 x 1 * * + β 1 μ 3 x 1 * * + β 1 μ 2 x 3 * * + β 1 μ 4 x 1 * * + β 3 μ 1 x 2 * * + β 1 μ 3 x 3 * * + β 3 μ 2 x 2 * * + β 1 μ 2 x 6 * * + β 3 μ 4 x 2 * * + β 3 μ 1 x 6 * * + β 1 μ 4 x 6 * * + β 1 μ 6 x 6 * * + β 3 μ 4 x 6 * * + β 3 μ 4 x 6 * * + β 3 μ 6 x 6 * * + β 1 β 3 x 6 * * 2 + β 1 β 3 x 2 * * x 3 * * + β 1 β 3 x 1 * * x 6 * * + β 1 β 3 x 2 * * x 6 * * + β 1 β 3 x 3 * * x 6 * * > 0. (60)

D 3 = μ 1 μ 2 μ 3 + μ 1 μ 2 μ 4 + μ 1 μ 3 μ 4 + μ 1 μ 2 μ 6 + μ 2 μ 3 μ 4 + μ 1 μ 3 μ 6 + μ 1 μ 4 μ 6 + μ 2 μ 3 μ 6 + μ 2 μ 4 μ 6 + μ 3 μ 4 μ 6 + β 1 μ 1 μ 2 x 1 * * + β 1 μ 1 μ 3 x 1 * * + β 1 μ 1 μ 2 x 3 * * + β 1 μ 1 μ 4 x 1 * * + β 1 μ 2 μ 3 x 1 * * + β 1 μ 1 μ 3 x 3 * * + β 1 μ 2 μ 4 x 1 * * + β 3 μ 1 μ 2 x 2 * * + β 1 μ 2 μ 3 x 3 * * + β 1 μ 3 μ 4 x 1 * * + β 3 μ 1 μ 3 x 2 * * + β 3 μ 1 μ 4 x 2 * * + β 3 μ 2 μ 3 x 2 * * + β 1 μ 1 μ 6 x 3 * * + β 3 μ 2 μ 4 x 2 * * + β 1 μ 2 μ 3 x 6 * * + β 1 μ 2 μ 6 x 3 * * + β 3 μ 3 μ 4 x 2 * * + β 1 μ 2 μ 4 x 6 * * + β 1 μ 3 μ 6 x 3 * * + β 3 μ 1 μ 3 x 6 * * + β 1 μ 3 μ 4 x 6 * * + β 3 μ 1 μ 4 x 6 * * + β 1 μ 2 μ 6 x 6 * * + β 1 μ 3 μ 6 x 6 * * + β 3 μ 1 μ 6 x 6 * * + β 3 μ 3 μ 4 x 6 * * + β 1 μ 4 μ 6 x 6 * * + β 3 μ 3 μ 6 x 6 * * + β 3 μ 4 μ 6 x 6 * * r μ 1 μ 2 x 4 * * r μ 1 μ 4 x 4 * * r μ 2 μ 4 x 4 * *

r μ 1 μ 6 x 4 * * r μ 2 μ 6 x 4 * * r μ 4 μ 6 x 4 * * + β 1 β 3 μ 3 x 6 * * 2 + β 1 β 3 μ 4 x 6 * * 2 + β 1 β 3 μ 6 x 6 * * 2 + β 1 2 μ 1 x 1 * * x 3 * * + β 1 2 μ 2 x 1 * * x 3 * * + β 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 x 3 * * x 6 * * + β 1 2 μ 3 x 3 * * x 6 * * + β 1 2 μ 6 x 3 * * x 6 * * + β 1 2 β 3 x 3 * * x 6 * * 2 N β 1 μ 1 μ 4 x 1 * * N β 1 μ 2 μ 4 x 1 * * N β 1 μ 3 μ 4 x 1 * * + β 1 β 3 μ 1 x 2 * * x 3 * * + β 1 β 3 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 x 1 * * x 6 * * + β 1 β 3 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 x 3 * * x 6 * * + β 1 β 3 μ 2 x 2 * * x 6 * * + β 1 β 3 μ 3 x 1 * * x 6 * * + β 1 β 3 μ 3 x 2 * * x 6 * * + β 1 β 3 μ 4 x 1 * * x 6 * * + β 1 β 3 μ 3 x 3 * * x 6 * * + β 1 β 3 μ 4 x 2 * * x 6 * * + β 1 β 3 μ 6 x 3 * * x 6 * * β 1 r μ 1 x 1 * * x 4 * * β 1 r μ 2 x 1 * * x 4 * * β 1 r μ 4 x 1 * * x 4 * * β 3 r μ 1 x 2 * * x 4 * * β 3 r μ 2 x 2 * * x 4 * * β 1 r μ 2 x 4 * * x 6 * * β 3 r μ 4 x 2 * * x 4 * * β 3 r μ 1 x 4 * * x 6 * *

β 1 r μ 4 x 4 * * x 6 * * β 1 r μ 6 x 4 * * x 6 * * β 3 r μ 4 x 4 * * x 6 * * β 3 r μ 6 x 4 * * x 6 * * β 1 β 3 r x 4 * * x 6 * * 2 + β 1 2 β 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 x 2 * * x 3 * * x 6 * * N β 1 β 3 μ 4 x 1 * * x 6 * * + N β 1 r μ 4 x 1 * * x 4 * * β 1 β 3 r x 1 * * x 4 * * x 6 * * β 1 β 3 r x 2 * * x 4 * * x 6 * * = μ 1 μ 2 μ 4 + μ 1 μ 2 μ 6 + μ 1 μ 4 μ 6 + μ 2 μ 4 μ 6 + β 1 μ 1 μ 2 x 1 * * + β 1 μ 1 μ 2 x 3 * * + β 1 μ 1 μ 4 x 1 * * + β 1 μ 2 μ 3 x 1 * * + β 1 μ 1 μ 3 x 3 * * + β 1 μ 2 μ 4 x 1 * * + β 3 μ 1 μ 2 x 2 * * + β 1 μ 2 μ 3 x 3 * * + β 3 μ 1 μ 4 x 2 * * + β 3 μ 2 μ 4 x 2 * * + β 1 μ 2 μ 4 x 6 * * + β 3 μ 1 μ 4 x 6 * * + β 1 μ 2 μ 6 x 6 * * + β 3 μ 1 μ 6 x 6 * * + β 1 μ 4 μ 6 x 6 * * + β 3 μ 4 μ 6 x 6 * * + β 1 β 3 μ 4 x 6 * * 2 + β 1 β 3 μ 6 x 6 * * 2 + β 1 2 μ 1 x 1 * * x 3 * * + β 1 2 μ 2 x 1 * * x 3 * * + β 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 x 3 * * x 6 * *

+ β 1 2 μ 3 x 3 * * x 6 * * + β 1 2 μ 6 x 3 * * x 6 * * + β 1 2 β 3 x 3 * * x 6 * * 2 + β 1 β 3 μ 1 x 2 * * x 3 * * + β 1 β 3 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 x 1 * * x 6 * * + β 1 β 3 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 x 3 * * x 6 * * + β 1 β 3 μ 2 x 2 * * x 6 * * + β 1 β 3 μ 4 x 1 * * x 6 * * + β 1 β 3 μ 3 x 3 * * x 6 * * + β 1 β 3 μ 4 x 2 * * x 6 * * + β 1 2 β 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 x 2 * * x 3 * * x 6 * * + N β 1 r μ 4 x 1 * * x 4 * * > 0. (61)

D 4 = μ 1 μ 2 μ 3 μ 4 + μ 1 μ 2 μ 3 μ 6 + μ 1 μ 2 μ 4 μ 6 + μ 1 μ 3 μ 4 μ 6 + μ 2 μ 3 μ 4 μ 6 + β 1 2 β 3 μ 3 x 3 * * x 6 * * 2 + β 1 2 β 3 μ 6 x 3 * * x 6 * * 2 + β 1 μ 1 μ 2 μ 3 x 1 * * + β 1 μ 1 μ 2 μ 4 x 1 * * + β 1 μ 1 μ 2 μ 3 x 3 * * + β 1 μ 1 μ 3 μ 4 x 1 * * + β 1 μ 2 μ 3 μ 4 x 1 * * + β 3 μ 1 μ 2 μ 3 x 2 * * + β 3 μ 1 μ 2 μ 4 x 2 * * + β 1 μ 1 μ 2 μ 6 x 3 * * + β 3 μ 1 μ 3 μ 4 x 2 * * + β 1 μ 1 μ 3 μ 6 x 3 * * + β 3 μ 2 μ 3 μ 4 x 2 * * + β 1 μ 2 μ 3 μ 6 x 3 * * + β 1 μ 2 μ 3 μ 4 x 6 * * + β 3 μ 1 μ 3 μ 4 x 6 * * + β 1 μ 2 μ 3 μ 6 x 6 * * + β 1 μ 2 μ 4 μ 6 x 6 * * + β 3 μ 1 μ 3 μ 6 x 6 * * + β 1 μ 3 μ 4 μ 6 x 6 * * + β 3 μ 1 μ 4 μ 6 x 6 * * + β 3 μ 3 μ 4 μ 6 x 6 * * r μ 1 μ 2 μ 4 x 4 * * r μ 1 μ 2 μ 6 x 4 * * r μ 1 μ 4 μ 6 x 4 * * r μ 2 μ 4 μ 6 x 4 * * + β 1 β 3 μ 3 μ 4 x 6 * * 2 + β 1 β 3 μ 3 μ 6 x 6 * * 2 + β 1 β 3 μ 4 μ 6 x 6 * * 2 + β 1 2 μ 1 μ 2 x 1 * * x 3 * * + β 1 2 μ 1 μ 3 x 1 * * x 3 * *

+ β 1 2 μ 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 μ 3 x 3 * * x 6 * * + β 1 2 μ 2 μ 6 x 3 * * x 6 * * + β 1 2 μ 3 μ 6 x 3 * * x 6 * * β 1 r μ 1 μ 2 x 1 * * x 4 * * β 1 r μ 1 μ 4 x 1 * * x 4 * * β 1 r μ 2 μ 4 x 1 * * x 4 * * β 3 r μ 1 μ 2 x 2 * * x 4 * * β 3 r μ 1 μ 4 x 2 * * x 4 * * β 3 r μ 2 μ 4 x 2 * * x 4 * * β 1 r μ 2 μ 4 x 4 * * x 6 * * β 3 r μ 1 μ 4 x 4 * * x 6 * * β 1 r μ 2 μ 6 x 4 x 6 β 3 r μ 1 μ 6 x 4 * * x 6 * * β 1 r μ 4 μ 6 x 4 * * x 6 * * β 3 r μ 4 μ 6 x 4 * * x 6 * * β 1 β 3 r μ 4 x 4 * * x 6 * * 2 β 1 β 3 r μ 6 x 4 * * x 6 * * 2 + β 1 2 β 3 μ 1 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 2 x 2 * * x 3 * * x 6 * * + β 1 2 β 3 μ 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 3 x 2 * * x 3 * * x 6 * * N β 1 μ 1 μ 2 μ 4 x 1 * * N β 1 μ 1 μ 3 μ 4 x 1 * * N β 1 μ 2 μ 3 μ 4 x 1 * * + β 1 β 3 μ 1 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 μ 3 x 2 * * x 3 * * + β 1 β 3 μ 2 μ 3 x 2 * * x 3 * * + β 1 β 3 μ 1 μ 3 x 1 * * x 6 * * + β 1 β 3 μ 1 μ 4 x 1 * * x 6 * *

+ β 1 β 3 μ 1 μ 3 x 3 * * x 6 * * + β 1 β 3 μ 2 μ 3 x 2 * * x 6 * * + β 1 β 3 μ 2 μ 4 x 2 * * x 6 * * + β 1 β 3 μ 3 μ 4 x 1 * * x 6 * * + β 1 β 3 μ 3 μ 4 x 2 * * x 6 * * + β 1 β 3 μ 1 μ 6 x 3 * * x 6 * * + β 1 β 3 μ 3 μ 6 x 3 * * x 6 * * N β 1 β 3 μ 1 μ 4 x 1 * * x 6 * * N β 1 β 3 μ 3 μ 4 x 1 * * x 6 * * + N β 1 r μ 1 μ 4 x 1 * * x 4 * * + N β 1 r μ 2 μ 4 x 1 * * x 4 * * β 1 β 3 μ 1 x 1 * * x 4 * * x 6 * * β 1 β 3 r μ 2 x 2 * * x 4 * * x 6 * * β 1 β 3 r μ 4 x 1 * * x 4 * * x 6 * * β 1 β 3 r μ 4 x 2 * * x 4 * * x 6 * * + N β 1 β 3 r μ 4 x 1 * * x 4 * * x 6 * * = μ 1 μ 2 μ 4 μ 6 + β 1 2 β 3 μ 3 x 3 * * x 6 * * 2 + β 1 2 β 3 μ 6 x 3 * * x 6 * * 2 + β 1 μ 1 μ 2 μ 4 x 1 * * + β 1 μ 1 μ 2 μ 3 x 3 * * + β 3 μ 1 μ 2 μ 4 x 2 * * + β 1 μ 2 μ 4 μ 6 x 6 * * + β 3 μ 1 μ 4 μ 6 x 6 * * + β 1 β 3 μ 4 μ 6 x 6 * * 2 + β 1 2 μ 1 μ 2 x 1 * * x 3 * * + β 1 2 μ 1 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 μ 3 x 3 * * x 6 * *

+ β 1 2 μ 2 μ 6 x 3 * * x 6 * * + β 1 2 μ 3 μ 6 x 3 * * x 6 * * + β 1 2 β 3 μ 1 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 2 x 2 * * x 3 * * x 6 * * + β 1 2 β 3 μ 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 3 x 2 * * x 3 * * x 6 * * + β 1 β 3 μ 1 μ 2 x 2 * * x 3 * * + β 1 β 3 μ 1 μ 3 x 2 * * x 3 * * + β 1 β 3 μ 2 μ 3 x 2 * * x 3 * * + β 1 β 3 μ 1 μ 4 x 1 * * x 6 * * + β 1 β 3 μ 1 μ 3 x 3 * * x 6 * * + β 1 β 3 μ 2 μ 4 x 2 * * x 6 * * + N β 1 r μ 1 μ 4 x 1 * * x 4 * * + N β 1 r μ 2 μ 4 x 1 * * x 4 * * + N β 1 β 3 r μ 4 x 1 * * x 4 * * x 6 * * > 0. (62)

D 5 = μ 1 μ 2 μ 3 μ 4 μ 6 + β 1 μ 1 μ 2 μ 3 μ 4 x 1 * * + β 3 μ 1 μ 2 μ 3 μ 4 x 2 * * + β 1 μ 1 μ 2 μ 3 μ 6 x 3 * * + β 1 μ 2 μ 3 μ 4 μ 6 x 6 * * + β 3 μ 1 μ 3 μ 4 μ 6 x 6 * * r μ 1 μ 2 μ 4 μ 6 x 4 * * + β 1 β 3 μ 3 μ 4 μ 6 x 6 * * 2 + β 1 2 μ 1 μ 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 μ 3 μ 6 x 3 * * x 6 * * + β 1 2 β 3 μ 3 μ 6 x 3 * * x 6 * * 2 N β 1 μ 1 μ 2 μ 3 μ 4 x 1 * * + β 1 β 3 μ 1 μ 2 μ 3 x 2 * * x 3 * * + β 1 β 3 μ 1 μ 3 μ 4 x 1 * * x 6 * * + β 1 β 3 μ 2 μ 3 μ 4 x 2 * * x 6 * * + β 1 β 3 μ 1 μ 3 μ 6 x 3 * * x 6 * * β 1 r μ 1 μ 2 μ 4 x 1 * * x 4 * * β 3 r μ 1 μ 2 μ 4 x 2 * * x 4 * * β 1 r μ 2 μ 4 μ 6 x 4 * * x 6 * * β 3 r μ 1 μ 4 μ 6 x 4 * * x 6 * * β 1 β 3 r μ 4 μ 6 x 4 * * x 6 * * 2 + β 1 2 β 3 μ 1 μ 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 2 μ 3 x 2 * * x 3 * * x 6 * * N β 1 β 3 μ 1 μ 3 μ 4 x 1 * * x 6 * * + N β 1 r μ 1 μ 2 μ 4 x 1 * * x 4 * *

β 1 β 3 r μ 1 μ 4 x 1 * * x 4 * * x 6 * * β 1 β 3 r μ 2 μ 4 x 2 * * x 4 * * x 6 * * + N β 1 β 3 r μ 1 μ 4 x 1 * * x 4 * * x 6 * * . = β 1 2 μ 1 μ 2 μ 3 x 1 * * x 3 * * + β 1 2 μ 2 μ 3 μ 6 x 3 * * x 6 * * + β 1 2 β 3 μ 3 μ 6 x 3 * * x 6 * * 2 + β 1 β 3 μ 1 μ 2 μ 3 x 2 * * x 3 * * + β 1 2 β 3 μ 1 μ 3 x 1 * * x 3 * * x 6 * * + β 1 2 β 3 μ 2 μ 3 x 2 * * x 3 * * x 6 * * + N β 1 r μ 1 μ 2 μ 4 x 1 * * x 4 * * + N β 1 β 3 r μ 1 μ 4 x 1 * * x 4 * * x 6 * * > 0. (63)

Note that D i > 0 , i = 1 , 2 , 3 , 4 , 5 .

By using the Routh-Hurwitz [15], the five roots of the characteristic equation will have negative real parts if and only if D i > 0 , i = 1 , 2 , 3 , 4 , 5 and D 1 D 2 D 3 > D 3 2 + D 1 2 D 4 .

One can also show that D 1 D 2 D 3 > D 3 2 + D 1 2 D 4 , so the Routh-Hurwitz criteria are always satisfied, implying that the eigenvalues are invariably negative. This means that the steady state is stable and that the infection is established chronically in an infected individual. Based on Hurwitz criterion, all roots of (59) have the negative real parts. Thus, the immune absent equilibrium ε * , is locally asymptotically stable.

It is also evident that if x 4 * * > μ 3 r , then b 33 > 0 and the endemic equilibrium point is unstable.

5.2. Global Stability of the Infection Persistent Equilibrium, ε**

In the ensuing proof we neglect the loss of HHV-8 virions due to the infection process of both the healthy B and progenitor cells.

Theorem 6 For the system (10)-(16), if R 0 > 1 , then the infection persistent steady state, ε * * , is GAS.

Proof. Consider the Lyapunov function, V 2 , defined by

V 2 = ( x 1 x 1 * * x 1 * * l n x 1 x 1 * * ) + β 2 r ( x 3 x 3 * * x 3 * * l n x 3 x 3 * * ) + ( x 4 x 4 * x 4 * l n x 4 x 4 * ) + β 1 x 1 * x 6 * N μ 4 x 4 * ( x 6 x 6 * x 6 * l n x 6 x 6 * ) (64)

Differentiating along the solutions of the system (10)-(16), we obtain

V 2 = ( 1 x 1 * * x 1 ) x ˙ 1 + β 2 r ( 1 x 3 * * x 3 ) x ˙ 3 + ( 1 x 4 * * x 4 ) x ˙ 4 + β 1 x 1 * * x 6 * * N μ 4 x 4 * * ( 1 x 6 * * x 6 ) x ˙ 6 (65)

Using the relations at equilibrium point and relabeling y i = x i x i * * , i = 1 , 2 , , 7 , we get

V 1 = ( 1 1 y 1 ) [ μ 1 x 1 * * ( 1 y 1 ) + β 1 x 1 * * x 6 * * ( 1 y 1 y 6 ) ] + β 2 r ( 1 1 y 3 ) [ r x 3 * * x 4 * * ( y 3 y 4 y 3 ) ] + ( 1 1 y 4 ) [ β 1 x 1 * * x 6 * * ( y 1 y 6 y 4 ) + β 2 x 3 * * x 4 * * ( y 4 y 3 y 4 ) ] + β 1 x 1 * * x 6 * * N μ 4 x 4 * * ( 1 1 y 6 ) [ N μ 4 x 4 * * ( y 4 y 6 ) ] = μ 1 x 1 * * ( 2 y 1 1 y 1 ) + β 1 x 1 * * x 6 * * ( 3 1 y 1 y 4 y 6 y 1 y 6 y 4 ) < 0. (66)

Since the arithmetic mean is greater than or equal to the geometric mean, 2 y 1 1 y 1 0 and 3 1 y 1 y 4 y 6 y 1 y 6 y 4 0 . And the global stability follows from the LaSalle’s invariance principle [14].

6. Numerical Simulation

We now demonstrate the existence of the virus free steady state that we observed in our mathematical analysis. The figure below depicts the dynamics when the reproduction number R 0 = 0.8403 < 1 .

7. Discussion and Conclusion

The purpose of this paper is to investigate the qualitative behaviour of a classical KS model such as positive invariance, boundedness, and global stability of steady states. The global stability of the virus free and infected steady states are established by direct Lyapunov method. The basic reproduction number R 0 is obtained, and it determines the dynamics of the classical KS model. We have proved that if the basic reproduction number, R 0 is greater than unity, then both the immune and infection persistent steady states are globally asymptotically stable (GAS). The immune absent equilibrium mirrors a scenario where an individual has a classical Kaposi’s sarcoma but the cell mediated immunity through the CD8+ cytotoxic T lymphocytes is dysfunctional. This is a popular phenomenon among individuals suffering from chronic viral diseases.

The simulations in Figure 2 illustrate the development of classical KS. In Figure 2(a) we depict the dynamics of uninfected and infected progenitor cells for a period of 800 days. After about 250 days the infected progenitor cell population rises rapidly before it attains a steady state value after about 600 days. The progenitor cells are erroneously infected by HHV-8 and these infected cells are ones that change their morphology and become KS cells. The KS cells begin to secrete inflammatory cytokines that act on these cells in either an autocrine or paracrine manner. In Figure 2(b) we illustrate the relation between the infected progenitor cells and the KS cells. Figure 2(c) depicts the dynamics of the B cells as targets of HHV-8 infection. The uninfected B cell concentration declines to low levels as HHV-8 infects them. The HHV-8 virion population is sustained by the infected B cell population. Most infected HHV-8 infected B cells remain latently infected and hence do not contribute to the growth of HHV-8. However, after activation, these cells undergo a lytic program in which HHV-8 is replicated which in turn contributes to the development of KS.

The simulations in Figure 3 illustrate the infection free steady state for the basic reproduction number R 0 = 0.8403 , less than a unity. In this particular case the disease fails to establish itself. Note that the infected populations decrease to zero while the healthy populations return to their original steady states, their levels, before the infection was introduced in the host.

We are convinced that if the lytic program could be disrupted then infected people could harbour the virus without ever developing the disease.

Figure 2. Population dynamics during the first 800 days with R 0 = 1.1765 . (a) Uninfected and infected progenitor cells; (b) Infected progenitor and KS cells; (c) Uninfected and infected B cells; (d) HHV-8 virions.

Figure 3. Population dynamics during the first 800 days with R 0 = 0.8403 . (a) Uninfected and infected progenitor cells; (b) Infected progenitor and KS cells; (c) Uninfected and infected B cells; (d) HHV-8 virions.

Funding

The author gratefully acknowledges the funding received from Simons Foundation based at Botswana International University of Science and Technology (BIUST).

Cite this paper: Chimbola, O. (2020) Mathematical Model of Classical Kaposi’s Sarcoma. Applied Mathematics, 11, 579-600. doi: 10.4236/am.2020.117040.
References

[1]   Foreman, K.E. (2001) Kaposi’s Sarcoma: The Role of HHV8 and HIV1 in Pathogenesis. Expert Reviews in Molecular Medicine, 3, 1-17.
https://doi.org/10.1017/S1462399401002733

[2]   Lungu, E.M., Massaro, T.J., Ndelwa, E., Ainea, N., Chibaya, S. and Malunguza, N.J. (2013) Mathematical Modeling of the HIV/Kaposi’s Sarcoma Coinfection Dynamics in Areas of High HIV Prevalence. Computational and Mathematical Methods in Medicine, 2013, Article ID: 753424, 12 p.

[3]   Martin, J.N., Ganem, D.E., Osmond, D.H., Page-Shafer, K.A., Macrae, D., et al. (1998) Sexual Transmission and the Natural History of Human Herpesvirus 8 Infection. The New England Journal of Medicine, 338, 948-954.
https://doi.org/10.1056/NEJM199804023381403

[4]   Kestens, L., Melbye, M., Biggar, R.J., Stevens, W.J., Piot, P., et al. (1985) Endemic African Kaposi’s Sarcoma Is Not Associated with Immunodeficiency. International Journal of Cancer, 36, 29-34.
https://doi.org/10.1002/ijc.2910360109

[5]   Foreman, K.E., Friborg, J., Kong, W.-P., et al. (1997) Propagation of a Human Herpesvirus from AIDS-Associated Kaposis Sarcoma. New England Journal of Medicine, 336, 163-171.
https://doi.org/10.1056/NEJM199701163360302

[6]   Sitas, F., Parkin, M., Chirenje, Z., Stein, L., Mqoqi, N. and Wabinga, H. (2013) Cancers. In: Disease and Mortality in Sub-Saharan Africa, 2nd Edition, The International Bank for Reconstruction and Development/The World Bank, Washington DC, Chapter 20.

[7]   Ahmed, R., et al. (1994) CD4+ T Cells Are Required to Sustain CD8+ Cytotoxic T-Cell Responses during Chronic Viral Infection. Journal of Virology, 68, 8056-8063.
https://doi.org/10.1128/JVI.68.12.8056-8063.1994

[8]   Szomolay, B. and Lungu, E.M. (2014) A Mathematical Model for the Treatment of AIDS-Related Kaposi’s Sarcoma. Journal of Biological Systems, 22, 495-522.
https://doi.org/10.1142/S0218339014500247

[9]   Gumel, A.B., Zhang, X., Shivakumar, P.N., Garba, M.L. and Sahai, B.M. (2002) A New Mathematical Model for Assessing Therapeutic Strategies for HIV Infection. Journal of Theoretical Medicine, 4, 147-155.
https://doi.org/10.1080/1027366021000003289

[10]   Huynh, G.T. and Adler, F.R. (2012) Mathematical Modelling the Age Dependence of Epstein-Barr Virus Associated Infectious Mononucleosis. Mathematical Medicine and Biology, 29, 245-261.
https://doi.org/10.1093/imammb/dqr007

[11]   Flore, O., Raffi, S., Ely, S., O’Leary, J.J., Hyjek, E.M. and Cesarman, E. (1998) Transformation of Primary Human Endothelial Cells by Kaposi’s Sarcoma-Associated Herpesvirus. Nature, 394, 588-592.
https://doi.org/10.1038/29093

[12]   Shapiro, M., Duca, K., Lee, K., Delgado-Eckert, E., Hawlins, J., Jarrah, A., Laubenbacher, R., Polys, N., Hadinoto, V. and Thorley-Lawson, D. (2008) A Virtual Look at Epstein-Barr Virus Infection: Simulation Mechanism. Journal of Theoretical Biology, 252, 633-648.
https://doi.org/10.1016/j.jtbi.2008.01.032

[13]   Castillo-Chavez, C., Feng, Z. and Huang, W. (2002) On the Computation of R0 and Its Role on Global Stability. In: Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, Vol. 125 of IMA, Springer, New York, 229-250.
https://doi.org/10.1007/978-1-4757-3667-0_13

[14]   LaSalle, J.P. (1976) The Stability of Dynamical Systems. Regional Conference Series in Applied Mathematics, SIAM, Philadelphia.

[15]   Willems, J.L. (1970) Stability Theory of Dynamical Systems. Nelson, New York.

 
 
Top