A REACTION DIFFUSION SYSTEM MODELING VIRUS DYNAMICS AND CTLS RESPONSE WITH CHEMOTAXIS

In this paper, we study the effect of chemotactic movement of CTLs on HIV-1 infection dynamics by a reaction diffusion system with chemotaxis. Choosing a typical chemosensitive function, we find that chemoattractive movement of CTLs due to HIV infection does not change stability of a positive steady state of the model, meaning that the stability/instability of the model remains the same as in the model without spatial effect. However, chemorepulsion movement of CTLs can destabilize the positive steady state as the strength of the chemotactic sensitivity increases. In this case, Turing instability occurs, which may result in Hopf bifurcation or steady state bifurcation, and spatial inhomogeneous pattern forms.

1. Introduction.Some living organisms or cells, such as somatic cells and lymphocytes, have the ability to detect certain chemicals in their environment and adapt their movement accordingly, moving either toward or away from the chemical stimulus.This phenomenon is called chemotaxis or generally chemosensitive movement.In the mathematical literature, the term chemotaxis is used broadly to describe general chemosensitive movement responses, including chemoattraction (positive chemotaxis) and chemorepulsion (negative chemotaxis).However, in the experimental community, for example, in leukocytes trafficking mechanism, the term chemotaxis is defined only as chemoattraction, that is, an active movement of leukocytes toward chemokinetic agents, while chemorepulsion is referred to as fugetaxis, describing the active movement of leukocytes away from chemokinetic agents.In this paper, we use chemotaxis to mean the general chemosensitive movement, either chemoattraction (positive chemotaxis) or chemorepulsion (positive chemotaxis).
Cytotoxic T lymphocytes (CTLs), or effector CD8 + T cells, play a critical role in host defense against human immunodeficiency virus type 1 (HIV-1) infection.Normally, effector T cells leave lymph nodes and traffic to peripheral sites of infection.However, in HIV-1 infection, the majority of HIV-1 replication occurs in lymphoid tissues.To implement their antiviral activity, CTLs must migrate reversely back into infected lymphoid tissues, and remain within them .Thus the recruitment of CTLs is very important for the clearance of HIV-1.The movement of lymphocytes between the circulatory system and specific tissues is coordinated by Chemokines and their receptors.For example, inflammatory chemokines guide effector T cells to exit lymphoid tissues and home to peripheral sites of infection.HIV-1 infection and replication in the lymphoid tissues changes the chemotactic and cellular environments.As HIV-1 disease progresses, the homing ability of CTLs to infected lymph nodes may be disrupted, due either to reduced lymph node chemokine levels or reduced CTLs chemokine receptor expression, and thus affect the cytotoxic effect of CTLs in advanced HIV-1 infection [5].
Many viruses encode chemotactically active proteins.For instance, the envelope protein gp120 of HIV-1 has been shown to act as a T-cell chemoattractant via binding to the chemokine receptor and HIV-1 coreceptor CXCR4.However, some studies [4,18] showed that HIV-specific CTLs move toward or away from the CXCR4-binding HIV-1 gp120 in a concentration-dependent manner.The high concentration of CXCR4-binding HIV-1 gp120 repels HIV-specific CTLs, while low concentration of gp120 attracts CTLs with specific interaction with CXCR4.The repellant activity of HIV-1 gp120 on CTLs causes the active movement of HIV-1-specific CTLs away from the site of infection, which allows the virus to evade immune recognition and invade immune system.
In this paper, we study the effect of chemotactic movement of CTLs during HIV-1 infection by mathematical modeling.We denote by T (x, t), I(x, t) and E(x, t) the population densities of uninfected CD4 + T cells, infected CD4 + T cells and CTLs at location x at time t respectively.Assuming that the virus population is at a quasi-steady state [15], we consider the following model.
Here, we assume that uninfected CD4 In this model, we assume that uninfected CD4 + T cells and infected CD4 + T cells move randomly, with the diffusion coefficients D T and D I respectively.In contrast, the diffusion of CTLs consists of two parts, the random diffusion and the chemotactic movement.The random diffusion coefficient is assumed to be D E , that is, the diffusion flux of CTLs is proportional to their density gradient J D = −D E ∇E.As mentioned above, the HIV-1 viral protein gp120, binding with the coreceptor CXCR4, acts as a chemoattractant or chemorepellant for CTLs.Thus, the chemotaxis flux J C of CTLs depends on the their own density, the density of HIV-1 viral protein gp120, and the density gradient of this protein.Here we assume that the density of viral protein gp120 binding to CXCR4 is proportional to the density of infected CD4 + T cells, and the chemotaxis flux of CTLs is J C = −EΨ(E, I)∇I.The derivation of this chemotactic term can be referred from the works of Painter [12] and Hillen [8].The function Ψ(E, I) represents the chemotactic response, which denotes chemoattraction (chemorepulsion) if it is negative (positive).Notice that movement of CD4 + T cells may be very slow comparing with CTLs, that is, D T D E , or they even do not diffuse at all in the lymphoid tissue.Here we assume that D T > 0 but very small compared with D E for the sake of mathematical consideration.
For the PDE model system (1), we consider the no-flux boundary conditions: The rest of the paper is organized as follows.In section 2 we discuss the wellposedness of the model.Linear stability of the steady states are shown in section 3. The conditions for Turing instability and pattern formation are derived in section 4. Numerical simulation about the stability of positive steady state, steady state bifurcation, Hopf bifurcation and pattern formation are shown in section 5. Finally, we present conclusions and discussions in section 6.
We consider (1)-( 2) with the following initial conditions Existence and uniqueness, as well as nonnegativity of solution of (1)-( 2) are obtained in the following two theorems.
The following theorem confirms that the existence of solution is indeed global if Theorem 2.2.Assume that D T = D I .Suppose that (T, I, E) is the solution obtained in Theorem 1, then it is a global solution to system (1).
Proof.In order to prove the global existence of the solution, by (ii) and (iii) of Theorem 1, we only need to prove that (T, I, E) is bounded from above.
From the T and I equations of (1), we see that where d m = min{d T , d T * }.By Lemma 1 in Lou and Zhao [9], h/d m is globally attractive in C( Ω, R) for the scalar parabolic equation The parabolic comparison theorem (see, e.g., Theorem 7.3.4 in [13]) implies that T + I is bounded on [0, τ ).This together with the non-negativity of T (t, x) and I(t, x) further implies that T (t, x) and I(t, x) are both bounded.We assume 0 ≤ T (t, x) ≤ T M , 0 ≤ I(t, x) ≤ I M .
Let E M = c ηd E I M .Define the differential operator P by For any solution of system (1), we have PE = 0.However, for E M = c ηd E I M sufficiently large, by (H), we have On the boundary ∂Ω, ∂E M ∂ν = 0. Therefore, E = E M is a upper solution of the E equation in system (1).By the comparison principle, we then obtain that E(t, x) ≤ E M .Therefore, the solution (T, I, E) is bounded from above, completing the proof.
Notice that we assumed D T = D I in the proof of Theorem 2, and this will be assumed in what follows.
3. Linear stability analysis.System (1) has three possible spatially homogeneous steady states: (i) the infection-free steady state S 0 = (h/d T , 0, 0) always exists; (ii) if R 0 := hβ/d I d T > 1, there exists a virus-established steady state S 1 = (T 1 , I 1 , 0), where there exists a positive steady state S * = (T * , I * , E * ), where Notice that if there is no saturation of CTL, that is η → 0, the positive steady state tends to where It can be verified that under the assumption cI 1 > d E , the following holds It is straightforward to verify that cI 1 > d E is equivalent to The formulas of the positive steady state (T * , I * , E * ) are very complex.Before discussing the stability of the steady states, we investigate the dependence of the positive steady state on parameters.One can easily verify that T * is an increasing function of p, c, d I , and a decreasing function of β, η, d T , d E ; I * is an increasing function of β, η, d E , and a decreasing function of p, c, d T , d I ; E * is an increasing function of c, and a decreasing function of η, d E .Here, we take the parameter c as an examples to demonstrate that as c increases, T * and E * increases while Furthermore, Therefore, In a similar way, we can show the dependence of T * , I * and E * on other parameters.
In what follows, we discuss the linear stability of the steady states.In the absence of spatial effect, we know that S 0 is locally asymptotically stable if R 0 < 1; S 1 is locally asymptotically stable if R 0 > 1 and cI 1 < d E ; S * is locally asymptotically stable if cI 1 > d E .In fact, these results can be obtained from the discussion below as a special case when the spatial effect disappears.
Let Ŝ = ( T , Î, Ê) be a steady state of the system (1).Then the linearization of system (1) at Ŝ is given by where and u = (u 1 , u 2 , u 3 , u 4 ) T , Ψ( Ŝ) := ÊΨ( Ê, Î).Notice that Ψ(S 0 ) = 0, Ψ(S 1 ) = 0 and The corresponding characteristic polynomial of the linearized system ( 5) is where k ≥ 0, called the wavenumbers or the wave modes, are the eigenvalues of Laplace operator on a finite domain with no-flux boundary conditions.For instance, in one-dimensional domain where n and m are integers.λ is the eigenvalue which determines temporal growth.The steady state Ŝ is linearly stable if Reλ < 0, for every eigenvalue λ of (6) (see, e.g., [11]).
Theorem 3.1.The infection-free steady state S 0 = (h/d T , 0, 0) is linearly stable if R 0 < 1, and unstable if R 0 > 1 Proof.For the infection-free steady state S 0 = (h/d T , 0, 0), we have and the characteristic equation of the linearized system at S 0 is Note that λ 1 < 0, λ 3 < 0, and then λ 2 > 0 for some small modes k, including k = 0, which means S 0 is unstable, completing the proof.
Theorem 3.2.The virus-persistence steady state Proof.For the steady state S 1 = (T 1 , I 1 , 0), we have Noticing that βT 1 = d I and d T + βI 1 = hβ/d I , the characteristic polynomial of the linearized system at S 1 is given by Other eigenvalues are determined by where Thus the roots of ( 7) have negative real parts for all k.Therefore, if R 0 > 1 and cI 1 < d E , the steady state 0 for some small k, which implies ( 7) has at least one eigenvalue with positive real part, completing the proof.9), (10) and (11) respectively; otherwise, it is unstable.
Proof.For the positive steady state S * , the characteristic polynomial is given by Here where where where By Routh-Hurwitz Criteria, Re(λ) < 0 for every eigenvalue λ if and only if b Notice that, in the absence of diffusion, that is, in the spatial homogeneous case In contrast, under the same condition, the homogeneous steady state S * can be unstable to small spatial perturbations when diffusion is present, ) < 0 for some modes k.This diffusion-driven instability is called Turing instability [11].
4. Turing instability and pattern formation.Among the two conditions b 3 (k 2 ) > 0 and b 1 (k 2 )b 2 (k 2 ) − b 3 (k 2 ) for k ≥ 0 in Theorem 3.3 that ensure the stability of S * , violation of each will lead to different consequences: when b 3 (k 2 ) changes from positive to negative, typically there will be steady state bifurcation resulting in occurrence of a spatially heterogeneous steady state; whereas when b 1 (k 2 )b 2 (k 2 ) − b 3 (k 2 ) changes from positive to negative, Hopf bifurcation will occur causing temporally periodic solution oscillating about the constant steady state S * (see, e.g., Yu [19]).Thus, it is worthwhile to explore the possible sign changes for those two quantities.
We expect spatial pattern formation when the homogeneous positive steady state S * loses its stability, that is when b To investigate conditions for Turing instability, we denote φ d (s) =: b 3 (s) and φ h (s) =: where s is assumed to be s ∈ R. We see that φ d (s) has a negative root, since and there is no bifurcation.We only need to worry the case of either d 2 < 0 or d 3 < 0 for which φ d (s) may have two (up to multiplicity) or no positive roots, by the Descarte's rule of signs, and in the latter case, b 3 (k 2 ) = φ d (s) remains positive for all s ≥ 0. A steady state bifurcation occurs if φ d (s) has two distinct positive roots.The conditions for the existence of two distinct positive roots of φ d (s) can be determined by the sign of φ d (s d + ), where s d + is one of the roots of φ d (s) = 3d 1 s 2 + 2d 2 s + d 3 , say According the the signs of d Then, φ h (s) can become negative only when it has two positive real roots, and this is implied by one of the following two conditions: Turing instability may occur if one of the conditions (C1)-(C4) is satisfied, resulting formation of patterns.Notice that the formula of φ d (s d + ) can be simplified to the following form applying the fact that φ d (s d + ) = 3d 1 (s d + ) 2 + 2d 2 s d + + d 3 = 0.A challenging problem for chemotaxis models is the blow-up of solutions in finite time.Blow-up of solutions implies the whole population concentrate in a single point in a finite time.To avoid a blow-up, various mechanisms are introduced [7,8,14,17].One mechanism is to consider the volume-filling effect [12,17], which supposes only a finite number of cells can be accommodated at any site.We assume the chemotactic response function has the form Ψ(E, I) = χ(I)q(E), where the chemotactic sensitivity function χ(I) is chosen to be a constant χ.The volumefilling effect is described by q(E) satisfying q(E M ) = 0, and q(E) ≥ 0 for all 0 ≤ E < E M , where E M is the upper bound of the CTLs population E in the model, which shows the maximum number of cells those can be accommodated at any site.The function q(E) is considered as the probability of the cells finding a space at its neighboring location depends upon the availability of space.In this paper, we choose a simple form q(E) = 1 − E E M for the volume-filling effect of CTLs to prevent the blow-up of solutions of the model.Thus we adopt the following form for the chemosensitivity function Ψ(E, I): , for a chemorepulsion system (13) and , for a chemoattraction system.( 14) Furthermore, we have , for a chemorepulsion system (15) and For a chemoattraction system, Ψ * < 0 implies that d 2 , d 4 , h 2 and h 4 are all positive, meaning that S * is always linearly stable and thus, Turing instability will not occur.However, for a chemorepulsion system, Turing instability may occur.To demonstrate this we choose χ, the strength of the chemotactic sensitivity, as the bifurcation parameter and numerically explore the occurrence of Turing instability.Naturally, one would expect that there exists a critical value χ c such that there is no pattern formation if χ < χ c , while pattern will be formed if χ > χ c .Notice that d 1 , d 4 , h 1 and h 4 are independent of χ, while d 2 (χ), d 3 (χ), h 2 (χ) and h 3 (χ) are linear strictly decreasing functions of χ, with d j (0) > 0, h j (0) > 0, d j (+∞) = −∞, and h j (+∞) = −∞, (j = 2, 3).Thus, d 2 2 (χ) − 3d 1 d 3 (χ) and h 2 2 (χ) − 3h 1 h 3 (χ) are quadratic functions, which tend to positive infinity as χ → +∞.Furthermore, Therefore, as χ increases, at least one of the conditions (C1)-(C4) holds.The Turing instability threshold values can be determined according to the conditions (C1)-(C4).
For example, to determine the threshold for steady state bifurcation, let χ d 2 and χ d 3 be the roots of d 2 (χ) = 0 and d 3 (χ) = 0 respectively, that is, ).Notice that when we derive the bifurcation thresholds, we assumed s (i.e k 2 ) to be continuum.However, with finite spatial domains, there is a discrete set of possible modes k as mentioned above.Therefore, the threshold values χ S c and χ H c obtained here may be not the exact bifurcation values.χ S c and χ H c give the lower bound of the bifurcation values.The exact bifurcation values may be somewhat greater than χ S c and χ H c , depending on the size of the domain and the shapes of φ d (s) and φ h (s).
In a special case, if the uninfected cells and infected cells cannot diffuse (D T = 0), we have It is easy to see that the Turing instability conditions change to (H1) If (H1) holds, a steady state bifurcation occurs, while there is a Hopf bifurcation if any of (H2) and (H3) is satisfied.The threshold value for steady bifurcation is χ S c = χ d 3 , the root of d 3 (χ) = 0. From the formulas of h 2 and h 3 , we know that the root χ h 2 of h 2 (χ) = 0 is always smaller than the root χ h 3 of h 3 (χ) = 0. Thus, the threshold value for Hopf bifurcation is χ H c = χ h 2 , according to the conditions (H2) and (H3).
5. Numerical simulation.In this section, we use the framework of section 4 to perform some numerical simulations to explore the aforementioned bifurcations.To this end, we choose χ as the bifurcation parameter and other baseline parameters as h = 10, d T = 0.1, d I = 0.2, d E = 0.1, β = 0.1, p = 0.5, c = 0.2, η = 0.01 and E M = 2000, then R 0 = 50 > 1, I 1 − d E c = 48.5 > 0, and S * = (64.0197,0.5620, 12.4039).The space is assumed to be one-dimensional and Ω = [0, lπ], then k 2 = n 2 /l 2 , where n is an integer and l = 3.In what follows, we suppose that k assumes continuous values (i.e.s is continuous), in order to derive the threshold values of χ for bifurcation.We choose Ψ(E, I) as ( 13), and χ as the bifurcation parameter.Fig. 5 shows that the positive steady state S * is stable as χ = 0.02 < χ H c .As χ increases so that χ > χ H c but χ < χ S c , the positive steady state S * becomes unstable, and Hopf bifurcation occurs.Fig. 6 shows the temporal and spatial evolution of T (x, t), I(x, t) and E(x, t).From the end part of the these time evolution figures (i.e., the figures in second row of Fig. 6), we see the form of spatial patterns and temporal periodicity of the solutions.Fig. 7 shows the periodic solutions at space locations x = 1.5π and x = 3π.The amplitudes of the periodic solutions vary in different space locations as shown in Fig. 7.We see that the amplitudes at x = 1.5π are greater than those at x = 3π for T (x, t), I(x, t) and E(x, t) respectively.When χ = 0.033 which exceeds the threshold value χ S c , the positive steady state also loses its stability and spatial heterogeneous pattern forms, which may be not periodic, as shown in Fig. 8.     I(x, t) and E(x, t).The homogenous positive steady state is unstable and spatial patterns form when χ = 0.06 (> χ S c ).
6. Conclusion and discussion.In this paper, we have studied the effect of chemotactic movement, including chemoattaction and chemorepulsion, of CTLs on the HIV-1 infection dynamics by a reaction diffusion model with chemotaxis.In the absence of spatial effect, that is, without random diffusion and chemotactic movement, the homogeneous positive steady state is locally stable.Choosing the typical chemotactic sensitivity function ( 14), we found that chemoattraction movement of CTLs also cannot destabilize the homogeneous positive steady state, and there is no heterogeneous pattern formation.In contrast, chemorepulsion movement of CTLs can lead to instability of the homogeneous positive steady state and formation of spatially heterogeneous patterns.Using the typical chemotactic sensitivity function (13) and choosing the strength of chemotactic sensitivity χ as the bifurcation parameter, we found that Turing instability occurs when χ exceeds some threshold.The bifurcation may be steady state bifurcation or Hopf bifurcation depending on χ and other parameters, such as the diffusion coefficients D T and D E .We identified four conditions (C1)-(C4) from which bifurcation thresholds may be derived.With a finite domain in one dimensional space, we demonstrated the stability of positive steady state, steady state bifurcation and Hopf bifurcation for different diffusion coefficients D T and D E and chemotactic sensitivities χ.
The chemotactic movement of CTLs leads to the instability of the homogeneous positive steady state and pattern formation or periodic solutions of the system for some range of parameter values.With stable patterns or periodic solutions, the concentration of infected cells and virus load cannot stabilize at a constant level, but show heterogeneous distributions or oscillations.This is important for experimental or clinic estimation of virus load.Due to the heterogeneous distributions or periodic oscillations, lower (or higher) virus load detected at a site and at a moment does not indicate the same lower (or higher) load at different sites for a long time.The oscillations of viral load levels in the plasma are also plausible under the effects of immune responses or delays in the virus infection dynamics (e.g., see [6]).
In order for CTLs to successfully control HIV-1 infection, they must home efficiently to infected tissue sites and migrate within the infected tissue to the virusinfected cells.In this paper, we only considered the chemotactic sensitivity functions ( 13) and ( 14) for chemorepulsion and chemoattraction models respectively.For the chemoattraction model with (14), the negativity of Ψ * drives the positive steady state to be stable.If other chemotactic sensitivity functions are chosen so that Ψ * is not always negative, then it is possible that Turing bifurcations may also occur for chemoattraction model.For example, for the Macroscopic form of Lapidus and Schiller (with receptor response law) [10], we have which is positive for some K.For small K, there may be Turing bifurcation and pattern formation.More intensive investigation is needed for this, which is left as possible future work.
In our chemorepulsion model, the chemotactic sensitivity is assumed to be constant and positive.However, some experiments (e.g., [4]) demonstrate that gp120 elicits bidirectional T-cell movement in a receptor-mediated, concentration-dependent manner, attracting CD8 + lymphocytes and HIV-specific CTLs maximally at low concentration of gp120 and repelling the same cells at a higher concentration of gp120.Therefore, the chemotactic sensitivity function should be negative (chemoattraction) for low concentration of gp120 and positive (chemorepulsion) for high concentration of gp120.Relating to our simplified model, the chemotactic sensitivity function should be negative (chemoattraction) for low level of I(x, t) and positive for high level of I(x, t).Considering this in our model, the analysis would be more difficult mathematically, and is thus also left for our future work.
∂T ∂t = D T ∆T + h − d T T − βT I, ∂I ∂t = D I ∆I + βT I − d I I − pIE,

2 and d 3 ,
we have the following three cases.(i) If d 3 < 0, then d 2 2 − 3d 1 d 3 > |d 2 |, s d + > 0 and s d − < 0. In this case, φ d (s) has two distinct positive roots, if and only if φ d (s d + ) < 0. (ii) If d 3 > 0 and d 2 < 0, then for different cases of d 2 2 − 3d 1 d 3 , we have that s) > 0 for s ∈ R. In the case (ii)-1, φ d (s) has two positive roots, if and only if φ d (s d + ) < 0, and in the cases (ii)-2 and (ii)-3, φ d (s) does not have any positive root and φ d (s) > 0 for s ≥ 0. (iii) If d 3 = 0 and d 2 < 0, then s d + > 0, s d − = 0, and φ d (s) has two distinct positive roots if and only if φ d (s d + ) < 0. From the above, we conclude that φ d (s) has two positive roots if one of the following conditions is satisfied.

Figure 10 .
Figure 10.Temporal and spatial evolution of T (x, t), I(x, t) and E(x, t).The homogenous positive steady state is unstable and spatial patterns form when χ = 0.06 (> χ S c ).
+ T cells are recruited at a constant rate h, and infected at a rate βT I.The infected cells are cleared by CTLs at a rate pIE.The CTLs proliferate in response to antigenic simulation with a rate cEI, and the rate of CTLs expansion saturates as the number of CTLs grows to relatively high numbers.1/η represents the saturation level.Uninfected CD4 + T cells, infected CD4 + T cells and CTLs are lost at rates d T T , d I I and d E E respectively.
the threshold for steady state bifurcation.Similarly, the threshold value χ H 1 or χ H 2 for Hopf bifurcation can be derived according to the conditions (C3) and (C4), which is determined by the roots of [φ h (s h + )](χ) = 0. From the above analysis, we know that S * may lose its stability leading to bifurcation when φ d (s d + ) = 0 or φ h (s h + ) = 0.The bifurcation occurs at critical values χ = χ S c and χ = χ H c which are determined by φ d (s d + ) χ S c = 0 and φ h (s h + ) χ H c = 0.Among these two critical values, χ S c is a steady state bifurcation value, while χ H c is a Hopf bifurcation value.If χ S c < χ H c and χ S c < χ, φ d (s) has two positive solutions s d 1 and s d 2 which give a range of unstable wave numbers k 2 ∈ (s