Hongyan Kang^*, Yinji Huang
College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.
DOI: 10.4236/jamp.2026.148154   PDF    HTML   XML   10 Downloads   66 Views   Citations

Abstract

This paper investigates the Hopf bifurcation problem for a class of reaction-diffusion models of HTLV-I infection that incorporate CTL immune responses and immunochemotaxis. Using the chemotaxis coefficient as the bifurcation parameter, this paper derives the conditions for the existence of a Hopf bifurcation at a positive equilibrium point and the existence of periodic solutions. Furthermore, it determines the bifurcation direction and stability of the periodic solutions: if ${\xi }_{k_0}^{H^{″}}\left(0\right)>0$ , the bifurcation direction is supercritical, and the periodic solution is stable; if ${\xi }_{k_0}^{H^{″}}\left(0\right)<0$ , the bifurcation direction is subcritical, and the periodic solution is unstable.

Keywords

Immune Chemotaxis, Hopf Bifurcation, Periodic Solutions, Stability

Share and Cite:

1. Introduction

Communicable diseases are primarily caused by pathogens such as bacteria, viruses, fungi, and parasites, and are transmitted through the air, water, food, and other routes. Their spread affects not only human physical and mental health but also social and economic development. Human T-cell leukemia virus type I (HTLV-I) is a pathogenic retrovirus that is prevalent primarily in the Caribbean, Japan, Central Africa, and South America [1]. Following infection, HTLV-I can persist in the host for a long time, causing HTLV-I-associated myelopathy (HAM), also known as tropical spastic paraplegia (TSP), a chronic inflammatory disease of the central nervous system (CNS) [2] [3]. Additionally, HTLV-I infection can lead to malignant hematological diseases such as adult T-cell leukemia/lymphoma (ATL) [4] [5].

Individuals infected with HTLV-I typically display no obvious clinical manifestations; only 0.25% - 3.8% develop HAM/TSP, and 2% - 3% develop ATL. These patients have very high levels of CD8⁺ cytotoxic T lymphocytes (CTLs) in their peripheral blood [6]. Experiments have demonstrated that CTLs can reduce viral load, clear infected cells, and protect the body [7] [8]; however, when CTL levels are excessively high, they release large amounts of toxins, leading to the symptoms of HAM/TSP [9]. Therefore, incorporating the CTL immune response into the development of a dynamic model of HTLV-I infection is of great significance for gaining a deeper understanding of the pathogenesis of ATL and HAM/TSP and for formulating treatment strategies.

In 1999, Wodarz et al. first constructed an HTLV-I infection model that accounted for interactions between healthy and infected CD4⁺T cells [10]. Subsequently, Gómez-Acevedo and Li demonstrated the existence of the backward branch described in reference [10] through their analysis [11]. In 2002, Wodarz et al. established an HTLV-I infection model incorporating CTL immune responses to study the impact of CTL immunity on viral replication [12], using the bilinear function $cyz$ to describe the proliferation of CTL immune cells. Sun and Wei studied a class of HTLV-I infection models incorporating a CTL immune response and demonstrated that the global dynamical behavior of the system is determined by two threshold parameters, $R_0$ and $R_1$ [13].

In 2012, Li and Shu [14] established and studied a model of HTLV-I infection involving a CTL immune response; they introduced a time delay into the model and demonstrated that this delay can destabilize the HAM/TSP equilibrium, leading to a Hopf bifurcation and stable periodic oscillations. Jia and Xu studied a class of HTLV-I models featuring a Beddington-DeAngelis incidence rate of $\beta xy/\left(1+a_{1}y+a_{2}y\right)$ and infection delays, investigating the effects of single and double delays on the system dynamics [15]. Chen et al. studied a class of infection models featuring Logistic growth of healthy cells and CTL immune proliferation described by the function $f\left(y,z\right)=cyz/\left(1+qz\right)\left(1+\epsilon z\right)$ , and investigated the local and global stability of the model’s equilibrium points [16]. Li and Zhou constructed a mathematical model featuring Logistic growth of healthy cells and an infection rate given by $\beta xy/1+\omega x$ , and investigated the existence of backward branches in the model [17]. Zhang et al. studied the global dynamical behavior of a class of HTLV-I infection models featuring intracellular time delays and saturated CTL immune responses, proving the existence and local asymptotic stability of viable equilibria [18]. Xu and Yang studied the dynamics of an HTLV-I infection model featuring a saturated immune response and immune impairment. By calculating the unactivated and activated reproduction numbers, they discussed the stability of equilibrium points with (and without) time delays, as well as branching [19].

Building on the research in references [16], [17], [19] and based on the work in reference [14], we account for host resource limitations by incorporating a logistic growth model for the number of healthy cells; furthermore, since the immune response does not increase indefinitely with rising viral load, we adopt a saturated CTL immune response. We aim to study a model of HTLV-I infection characterized by logistic growth and a saturated CTL immune response

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(1-\frac{x}{K}\right)-\beta xy,\\ \frac{\text{d}y}{\text{d}t}=\beta xy-{\mu }_{1}y-pyz,\\ \frac{\text{d}z}{\text{d}t}=\frac{cyz}{1+qz}-{\mu }_{2}z,\end{array}$ (1.1)

In this model, $x$ , $y$ , and $z$ represent, respectively, the densities of healthy CD4⁺T cells, infected CD4⁺T cells, and HTLV-I-specific CD8⁺CTL cells at time $t$ ; $r$ represents the natural proliferation rate of healthy CD4⁺T cells; $K$ represents the carrying capacity of healthy CD4⁺T cells; $\beta$ represents the viral infection rate; $p$ represents the intensity of the CTL response, i.e., the rate at which the CTL immune response clears infected cells; $c$ represents the proliferation rate of CTL s; $q>0$ represents the CTL suppression rate, that is, CTL cells are stimulated by the virus and are generated at a rate of $\frac{cyz}{1+qz}$ ; ${\mu }_1$ and ${\mu }_2$ represent the natural death rates of infected CD4⁺T cells and CTL cells, respectively.

For convenience in calculation, Equation (1.1) is nondimensionalized as follows:

$u=\frac{1}{K}x,v=\frac{1}{K}y,w=qz,\tau =rt.$

Define the dimensionless parameters:

$a=\frac{\beta K}{r},b=\frac{{\mu }_1}{r},d=\frac{p}{rq},f=\frac{cK}{r},e=\frac{{\mu }_2}{r}.$

Continuing to use $t$ to denote $\tau$ , Equation (1.1) can be written as

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=u\left(1-u\right)-auv,\\ \frac{\text{d}v}{\text{d}t}=auv-bv-dvw,\\ \frac{\text{d}w}{\text{d}t}=\frac{fvw}{1+w}-ew.\end{array}$ (1.2)

Individual cells often exhibit significant spatial heterogeneity; therefore, spatial heterogeneity plays an indispensable role in infectious disease modeling and helps elucidate the spatial complexity of infectious disease dynamics. References [20] [21] have investigated the role of convection in the dynamics of infectious disease models, while references [22]-[24] have discussed the impact of chemotaxis on infectious disease modeling. The introduction of chemotaxis makes the established infectious disease models more realistic and can significantly enhance the effectiveness of disease prevention. Therefore, considering the impact of spatial heterogeneity on the dynamics between immune cells and infected cells, this paper investigates an HTLV-I infection model incorporating immune chemotaxis:

$\{\begin{array}{ll}{u}_t=d_1\Delta u+u\left(1-u\right)-auv,\hfill & x\in \Omega ,t>0,\hfill \\ v_t=d_2\Delta v+auv-bv-dvw,\hfill & x\in \Omega ,t>0,\hfill \\ w_t=d_3\Delta w-\xi \nabla \cdot \left(w\nabla v\right)+\frac{fvw}{1+w}-ew,\hfill & x\in \Omega ,t>0,\hfill \\ \frac{\partial u}{\partial \nu }=\frac{\partial v}{\partial \nu }=\frac{\partial w}{\partial \nu }=0,\hfill & x\in \partial \Omega ,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),v\left(x,0\right)=v_0\left(x\right),w\left(x,0\right)=w_0\left(x\right),\hfill & x\in \Omega ,\hfill \end{array}$ (1.3)

Here, $-\xi \nabla \cdot \left(w\nabla v\right)$ describes the movement of CTL immune cells toward regions of high infected cell concentration. From a biological perspective, $\xi >0$ indicates that CTL immune cells move toward regions of high infected cell density, whereas $\xi <0$ implies that CTL immune cells move away from the habitats of infected cells to prevent large numbers of infected cells from aggregating and forming a defense.

This paper primarily investigates the existence of Hopf branches in model (1.3), as well as the branch direction and stability of periodic solutions on these branches.

2. Existence and Stability of Periodic Solutions

The positive equilibrium of model (1.2) is

$u^{\text{*}}=\frac{b+d{w}^{\text{*}}}{a},v^{\text{*}}=\frac{e\left(1+w^{\text{*}}\right)}{f},w^{\text{*}}=\frac{bf+a^{2}e}{df+a^{2}e}\left(\frac{af}{bf+a^{2}e}-1\right).$

When $\frac{af}{bf+a^{2}e}>1$ , the positive equilibrium of model (1.2) exists.

Let $U^{\text{*}}=\left(u^{\text{*}},v^{\text{*}},w^{\text{*}}\right)$ be a positive equilibrium point of model (1.3). At $U^{\text{*}}$ , the following conclusions hold:

Let $\Omega \subset {ℝ}^n$ $\left(n\ge 1\right)$ be a bounded domain with smooth boundary, and let the initial data $\left(u_0,v_0,w_0\right)\in {\left[C\left(\overline{\Omega }\right)\right]}^3$ be nonnegative and not identically zero. Under homogeneous Neumann boundary conditions, the theory of parabolic equations ensures that system (1.3) admits a unique local classical solution, and $u,v,w>0$ for all $x\in \Omega$ and $t\in \left(0,T_{\text{max}}\right)$ .

Theorem 2.1 Suppose that $\frac{af}{bf+a^{2}e}>1$ holds, then

1) When $\xi >0$ , $\left(u^{\text{*}},v^{\text{*}},w^{\text{*}}\right)$ is locally asymptotically stable;

2) When $0>\xi >{\xi }_0=\underset{k\in {\mathbb{N}}^{+}}{\text{max}}\left\{{\xi }_k^S,{\xi }_k^H\right\}$ , $\left(u^{\text{*}},v^{\text{*}},w^{\text{*}}\right)$ is locally asymptotically stable;

3) When $\xi <{\xi }_0=\underset{k\in {\mathbb{N}}^{+}}{\text{max}}\left\{{\xi }_k^S,{\xi }_k^H\right\}<0$ , then $\left(u^{\text{*}},v^{\text{*}},w^{\text{*}}\right)$ is unstable. Where

${\xi }_k^S=-\frac{d_1{d}_2{d}_3{\mu }_k^3+\left(u^{\text{*}}{d}_2{d}_3+m{d}_1{d}_2\right){\mu }_k^2+\left(u^{\text{*}}m{d}_2+dk{w}^{\text{*}}{d}_1+a^2{u}^{\text{*}}{v}^{\text{*}}{d}_3\right){\mu }_k+e_3}{d{v}^{\text{*}}{w}^{\text{*}}{\mu }_k\left(d_1{\mu }_k+u^{\text{*}}\right)},$

${\xi }_k^H=-\frac{{\beta }_1{\mu }_k^3+{\beta }_2{\mu }_k^2+{\beta }_3{\mu }_k+{\beta }_4}{{\mu }_{k}d{v}^{\text{*}}{w}^{\text{*}}\left(d_2{\mu }_k+d_3{\mu }_k+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)\right)},$

${\beta }_1=\left(d_1+d_2\right)\left(d_1+d_3\right)\left(d_2+d_3\right),$

$\begin{array}{c}{\beta }_2=e_1\left(d_1{d}_2+d_1{d}_3+d_2{d}_3\right)+e_1{d}_2\left(d_1+d_2+d_3\right)\\ +\left(u^{\text{*}}{d}_3+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_1\right)\left(d_1+d_3\right),\end{array}$

$\begin{array}{c}{\beta }_3=e_1\left(e_1{d}_2+u^{\text{*}}{d}_3+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_1\right)+u^{\text{*}}\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)\left(d_1+d_3\right)\\ +de{w}^{\text{*}}\left(d_2+d_3\right)+a^2{u}^{\text{*}}{v}^{\text{*}}\left(d_1+d_2\right),\end{array}$

$\begin{array}{c}{\beta }_4=u^{\text{*}2}\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)+a^2{u}^{\text{*}2}{v}^{\text{*}}+u^{\text{*}}{\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)}^2\\ +\frac{df{v}^{\text{*}}{w}^{\text{*}}}{1+w^{\text{*}}}\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right).\end{array}$

We now discuss the existence and stability of periodic solutions for model (1.3). For computational convenience, this paper primarily considers the one-dimensional case of (1.3), namely the model

$\{\begin{array}{ll}{u}_t=d_1{u}^{″}+u\left(1-u\right)-auv,\hfill & x\in \left(0,l\right),t>0,\hfill \\ v_t=d_2{v}^{″}+auv-bv-dvw,\hfill & x\in \left(0,l\right),t>0,\hfill \\ w_t={\left(d_3{w}^{\prime }-\xi w{v}^{\prime }\right)}^{\prime }+\frac{fvw}{1+w}-ew,\hfill & x\in \left(0,l\right),t>0,\hfill \\ u^{\prime }\left(x,t\right)=v^{\prime }\left(x,t\right)=w^{\prime }\left(x,t\right)=0,\hfill & x=0,l,t>0,\hfill \\ u\left(x,0\right)=u_0\left(x\right),v\left(x,0\right)=v_0\left(x\right),w\left(x,0\right)=w_0\left(x\right),\hfill & x=0,l,\hfill \end{array}$ (2.1)

The symbol $\text{'}$ denotes the derivative with respect to $x$ . Combining Theorem 2.1, we see that the periodic solution of model (2.1) branches off from the positive

equilibrium point at $\xi ={\xi }_k^H$ . When $\xi$ passes through ${\xi }_0=\underset{k\in {\mathbb{N}}^{+}}{\text{max}}\left\{{\xi }_k^S,{\xi }_k^H\right\}$ , the equilibrium point becomes unstable via a Hopf bifurcation. To apply bifurcation theory to model (2.1) at $\xi ={\xi }_k^H$ , this paper must verify that the real part of the characteristic root at $\xi ={\xi }_k^H$ lies on the imaginary axis. When $\xi ={\xi }_k^H$ and $Q\left(\xi ,k\right)>0$ , the linearization matrix of model (1.3) has the following purely imaginary roots

${\lambda }_1^H\left({\xi }_k^H,k\right)=-P\left(\xi ,k\right)<0,{\lambda }_{2,3}^H\left({\xi }_k^H,k\right)=\pm i\sqrt{-Q\left({\xi }_k^H,k\right)},$

where

$\begin{array}{l}P\left(\xi ,k\right)=\left(d_1+d_2+d_3\right){\mu }_k+e_1,\\ Q\left(\xi ,k\right)=\left(d_1{d}_2+d_1{d}_3+d_2{d}_3\right){\mu }_k^2+e_2+\left(e_1{d}_2+u^{\text{*}}{d}_3+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_1+d{v}^{\text{*}}{w}^{\text{*}}\xi \right).\end{array}$ (2.2)

Thus, model (1.3) may generate a Hopf branch at a positive equilibrium point. Below, we explain when $Q\left({\xi }_k^H,k\right)>0$ . Let ${\xi }_k^M$ be

the unique root of $Q\left({\xi }_k^H,k\right)=0$ , whose expression is

${\xi }_k^M=-\frac{\left(d_1{d}_2+d_1{d}_3+d_2{d}_3\right){\mu }_k^2+e_2+\left(e_1{d}_2+u^{\text{*}}{d}_3+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_1\right)}{d{v}^{\text{*}}{w}^{\text{*}}},$ (2.3)

where ${\mu }_k={\left(\frac{k\pi }{l}\right)}^2$ , $k=0,1,2,\cdots$ , Calculations show that the following order relations hold among ${\xi }_k^S$ , ${\xi }_k^M$ , and ${\xi }_k^H$ .

Lemma 2.1 For $\forall k\in {\mathbb{N}}^{+}$ , either ${\xi }_k^S<{\xi }_k^M<{\xi }_k^H$ or ${\xi }_k^H<{\xi }_k^M<{\xi }_k^S$ . Furthermore, if the former holds, then $Q\left({\xi }_k^H\right)>0>Q\left({\xi }_k^S\right)$ . If the latter holds, then $Q\left({\xi }_k^S\right)>0>Q\left({\xi }_k^H\right)$ .

By Lemma 2.1, the linearization matrix of model (1.3) has a pair of purely imaginary roots if and only if ${\xi }_k^S<{\xi }_k^H$ . This implies that a Hopf branch of model (2.1) at $\left(U^{\text{*}},{\xi }_k^H\right)$ can occur only if ${\xi }_k^S<{\xi }_k^H$ . Therefore, in the subsequent analysis of Hopf branches, this paper consistently assumes that ${\xi }_k^S<{\xi }_k^H$ .

2.1. Hopf Branch

This subsection primarily proves the existence of a Hopf branch for model (2.1) under the assumption that ${\xi }_k^S<{\xi }_k^H$ . To apply the branch theorem [25] at the point ${\xi }_k^H$ , it is necessary to prove that the transversality condition holds. Similar to Theorem 3.4 in [26], the existence of non-trivial periodic solutions for model (2.1) yields the following result.

Theorem 2.2 Assume that ${\xi }_k^S<{\xi }_k^H<0,k\in {\mathbb{N}}^{+}$ and ${\xi }_i^H\ne {\xi }_k^H,\forall i\ne k,i\in {\mathbb{N}}^{+}$ . Then, there exists a normal constant $\theta$ and a unique single-parameter family of non-trivial periodic orbits

${\omega }_k\left(ϵ\right)=\left(U_k\left(ϵ,x,t\right),T_k\left(ϵ\right),{\xi }_k\left(ϵ\right)\right):ϵ\in \left(-\theta ,\theta \right)\to C^3\left(ℝ,{\xi }^3\right)\times {ℝ}^{+}\times {ℝ}^{+},$

and satisfying $\left(U_k\left(0,x,t\right),T_k\left(0\right),{\xi }_k\left(0\right)\right)=\left(U^{\text{*}},\frac{2\pi }{{\nu }_0},{\xi }_k^H\right)$ and

$U_k\left(ϵ,x,t\right)=U^{\text{*}}+ϵ\left(V_k^{+}{\text{e}}^{i{\nu }_{0}t}+V_k^{-}{\text{e}}^{-i{\nu }_{0}t}\right)\cos \frac{k\pi x}{l}+o\left(ϵ\right).$

where $\left(U_k\left(ϵ,x,t\right),{\xi }_k\left(ϵ\right)\right)$ is a non-trivial solution to model (2.1), and $U_k\left(ϵ,x,t\right)$ is a periodic solution with respect to time $t$ , with period

$T_k\left(ϵ\right)\approx \frac{2\pi }{{\nu }_0},\text{where}{\nu }_0=\sqrt{Q\left({\xi }_k^H,k\right)},$

and $\left(V_k^{\pm },\pm i{\nu }_0\right)$ is a pair of eigenvalues and eigenvectors of the linearization matrix of model (1.3). Furthermore, for all ${ϵ}_1\ne {ϵ}_2\in \left(-\theta ,\theta \right)$ , we have ${\omega }_k\left({ϵ}_1\right)\ne {\omega }_k\left({ϵ}_2\right)$ , and all non-trivial periodic orbits of model (2.1) in the vicinity of $\left(U^{\text{*}},{\xi }_k^H\right)$ must lie on the orbits ${\omega }_k\left(ϵ\right)$ , $ϵ\in \left(-\theta ,\theta \right)$ . In other words, in the vicinity of ${\omega }_k\left(ϵ\right)$ , if for some $\xi \in {ℝ}^{+}$ and a small constant $\delta >0$ , the model (2.1) has a non-trivial periodic solution $\widehat{U}\left(x,t\right)$ of period $T$ such that $\left|\xi -{\xi }_k^H\right|<\delta$ , $\left|T-2\pi /{\nu }_0\right|<\delta$ , and ${\max }_{t\in {ℝ}^{+},x\in \overline{\Omega }}\left|\widehat{U}\left(x,t\right)-U^{\text{*}}\right|<\delta$. Then, there exist ${ϵ}_0\in \left(-\theta ,\theta \right)$ and some ${\sigma }_0\in \left[0,2\pi \right)$ , such that $\left(T,\xi \right)=\left(T_k\left({ϵ}_0\right),{\xi }_k^H\left({ϵ}_0\right)\right)$ , $\widehat{U}\left(x,t\right)=U_k\left({ϵ}_0,x,t+{\sigma }_0\right)$ .

Proof The method is inspired by Theorem 6.1 in Reference [26]. By Theorem 2.1, if $\xi ={\xi }_k^S<{\xi }_k^H$ , then $Q\left({\xi }_k^H,k\right)>0$ , in which case the linearization matrix of model (1.3) has a pair of purely imaginary roots ${\lambda }_{2,3}=\pm \text{i}\sqrt{Q\left({\xi }_k^H,k\right)}$ . Since ${\xi }_i^H\ne {\xi }_k^H$ for any $i\ne k$ , the linearized matrix of (1.3) has no eigenvalues of the form $im{\nu }_0$ , where $m\in {\mathbb{N}}^{+}\left\{\pm 1\right\}$ . Furthermore, since ${\xi }_k^H\ne {\xi }_k^S$ for all $k\in \mathbb{N}$ , 0 is also not an eigenvalue of the linearization matrix of model (1.3) when $\xi ={\xi }_k^H$ . For $\xi$ in the vicinity of ${\xi }_k^H$ , let the three roots of the linearized matrix of model (1.3) be $\kappa \left(\xi \right)$ , $\pm \text{i}\nu \left(\xi \right)$ , and $\beta \left(\xi \right)$ , respectively, such that $\kappa \left({\xi }_k^H\right)=0$ and $\nu \left({\xi }_k^H\right)={\nu }_0>0$ . Based on the relationship between the coefficients and roots of the characteristic equation of model (1.3), we obtain

$\{\begin{array}{l}-P\left(\xi \right)=2\kappa \left(\xi \right)+\beta \left(\xi \right),\\ Q\left(\xi \right)={\kappa }^2\left(\xi \right)+{\nu }^2\left(\xi \right)+2\kappa \left(\xi \right)\beta \left(\xi \right),\\ -R\left(\xi \right)=\left({\kappa }^2\left(\xi \right)+{\nu }^2\left(\xi \right)\right)\beta \left(\xi \right).\end{array}$ (2.4)

Here, $P\left(\xi \right)$ and $Q\left(\xi \right)$ are as shown in (2.2) above, and $R\left(\xi \right)$ is given by:

$\begin{array}{c}R\left(\xi \right)=d_1{d}_2{d}_3{\mu }_k^3+\left(u^{\text{*}}{d}_2{d}_3+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_1{d}_2+d{v}^{\text{*}}{w}^{\text{*}}{d}_1\xi \right){\mu }_k^2+e_3\\ +\left(u^{\text{*}}\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)d_2+de{w}^{\text{*}}{d}_1+a^2{u}^{\text{*}}{v}^{\text{*}}{d}_3+d{u}^{\text{*}}{v}^{\text{*}}{w}^{\text{*}}\xi \right){\mu }_k\end{array}$

Differentiating each term in Equation (2.4) with respect to $\xi$ yields

$2{\kappa }^{\prime }\left(\xi \right)+{\beta }^{\prime }\left(\xi \right)=0,$ (2.5)

$\begin{array}{l}\left(\begin{array}{cc}\kappa \left(\xi \right)-\beta \left(\xi \right)& 2\nu \left(\xi \right)\\ {\kappa }^2\left(\xi \right)+{\nu }^2\left(\xi \right)-\kappa \left(\xi \right)\beta \left(\xi \right)& 2\nu \left(\xi \right)\beta \left(\xi \right)\end{array}\right)\left(\begin{array}{c}{\beta }^{\prime }\left(\xi \right)\\ {\nu }^{\prime }\left(\xi \right)\end{array}\right)\\ =\left(\begin{array}{c}d{v}^{\text{*}}{w}^{\text{*}}{\left(\frac{k\pi }{l}\right)}^2\\ -d{v}^{\text{*}}{w}^{\text{*}}{\left(\frac{k\pi }{l}\right)}^2\left(d_1{\left(\frac{k\pi }{l}\right)}^2+u^{\text{*}}\right)\end{array}\right)\end{array}$ (2.6)

Given that $\kappa \left({\xi }_k^H\right)=0$ and $\beta \left({\xi }_k^H\right)=-P\left({\xi }_k^H\right)$ , solving Equation (2.6) at $\xi ={\xi }_k^H$ yields

${\beta }^{\prime }\left({\xi }_k^H\right)=-\frac{d{v}^{\text{*}}{w}^{\text{*}}{\left(\frac{k\pi }{l}\right)}^2\left(\left(d_2+d_3\right){\left(\frac{k\pi }{l}\right)}^2+\left(e-\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}\right)\right)}{{\nu }_0^2+P^2\left({\xi }_k^H\right)}<0$ (2.7)

and

${\kappa }^{\prime }\left({\xi }_k^H\right)=-\frac{1}{2}{\beta }^{\prime }\left({\xi }_k^H\right)>0.$ (2.8)

This establishes the transverse condition, and the theorem is thus proved.

From a dynamical perspective, under normal physiological conditions with $\xi >0$ , chemotactic factors secreted by infected cells guide CTL to migrate and accumulate at infection foci, establishing a stable negative feedback regulation between viral proliferation and immune clearance, and the system maintains a stable positive equilibrium. When $\xi <0$ , repulsive chemotaxis reduces the effective accumulation of CTL at infection sites and attenuates the negative feedback regulation of viral proliferation by immune clearance, thereby destabilizing the system and inducing periodic oscillations. This dynamical feature may correspond to the pathological scenario of chronic HTLV-I infection, in which an immunosuppressive microenvironment impairs the directed recruitment of CTLs, resulting in periodic fluctuations in viral load and immune response levels.

The subsequent numerical simulations are used to verify that model (2.1) produces spatially inhomogeneous periodic solutions.

Example 1 We choose the parameters $a=4.8986$ ; $b=0.6333$ ; $f=11.3466$ ; $k=0.9384$ ; $d=0.6333$ , $d_1=0.0111$ , $d_2=0.0527$ , $d_3=2.6002$ , we obtain $U^{\text{*}}\approx \left(0.4340,0.1155,0.3970\right)$ . Setting $\xi =-14<{\xi }_0\approx -13.12$ , the systems positive equilibrium point becomes unstable. On the one-dimensional interval $\Omega =\left(0,10\right)$ , the spatial domain is discretized using a uniform grid with 100 grid points, and the time interval is taken as $t\in \left[0,250\right]$ with 1501 output time levels. The initial condition is chosen as a small perturbation around the positive equilibrium $U^{\text{*}}$ with perturbation amplitude 0.02, and the spatial distribution is given by the eigenfunction corresponding to the most unstable mode of the linearized system. Plotting $\left(u,v,w\right)\left(x,t\right)$ yields the spatiotemporal pattern shown in Figure 1.

2.2. Stability of Branch-Periodic Solutions

We now investigate the stability of the periodic solutions $U_k\left(ϵ,x,t\right)$ established in Theorem 2.2. Suppose that ${\xi }_k^S<{\xi }_{k_0}^H=\underset{k\in {\mathbb{N}}^{+}}{\text{max}}{\xi }_k^H$ , and that all conditions in Theorem 2.2 hold; then the following stability conclusion holds.

Theorem 2.3 Assuming that all conditions of Theorem 2.2 are satisfied, the following conclusions hold for model (2.1):

1) If ${\xi }_{k_0}^H<\underset{k\in {\mathbb{N}}^{+}}{\text{max}}{\xi }_k^H$ , then when ${\xi }_{k_0}^{H^{″}}\left(0\right)>0$ (${\xi }_{k_0}^{H^{″}}\left(0\right)<0$ ), the branching direction is supercritical (subcritical), and the periodic solution is stable (unstable).

Figure 1. Space-time diagram of system (2.1) in one dimension.

2) For all $k\ne k_0$ , the periodic solution is always unstable.

Proof Define $U_k\left(ϵ,x,t\right)=\left(u_k\left(ϵ,x,t\right),v_k\left(ϵ,x,t\right),w_k\left(ϵ,x,t\right)\right)$ , and let $\left(U_k\left(ϵ,x,t\right),T_k\left(ϵ\right),{\xi }_k^H\left(ϵ\right)\right)$ be the periodic solution on the branch ${\omega }_k\left(ϵ\right)$ obtained from Theorem 2.2. Rewrite model (2.1) in the following abstract form

$\frac{\text{d}{U}_k}{\text{d}t}=\mathcal{G}\left(U_k,{\xi }_k^H\left(ϵ\right)\right),$

where

$\mathcal{G}\left(U_k,{\xi }_k^H\left(ϵ\right)\right)=\left(\begin{array}{c}{d}_1{u^{″}}_k+u_k\left(1-u_k\right)-a{u}_k{v}_k\\ d_2{v^{″}}_k+a{u}_k{v}_k-b{v}_k-d{v}_k{w}_k\\ d_3{w^{\prime }}_k-{\xi }_k^H{\left(w_k{v}_k\right)}^{\prime }+\frac{f{v}_k{w}_k}{1+w_k}-e{w}_k\end{array}\right).$

Applying differentiation with respect to $t$ to the above abstract system, and letting ${\dot{U}}_k=\frac{\text{d}{U}_k}{\text{d}t}$, we obtain

$\frac{\text{d}{\dot{U}}_k}{\text{d}t}=\mathcal{G}\left(U_k,{\xi }_k\left(ϵ\right)\right){\dot{U}}_k.$

It can be observed that 0 is the Floquet index, and 1 is the Floquet multiplier of $U_k$ .

Substituting the perturbation solution $U_k+w{\text{e}}^{-\kappa t}$ and linearizing the periodic approximation on the branch ${\omega }_k\left(ϵ\right)$ where $w$ is a sufficiently small $T$ -periodic function and $\kappa =\kappa \left(ϵ\right)$ is a continuously differentiable function of $ϵ$ we obtain

$\frac{\text{d}w\left(ϵ,t\right)}{\text{d}t}={\mathcal{G}}_U\left(U_k,{\xi }_k^H\left(ϵ\right)\right)w\left(ϵ,t\right)+\kappa \left(ϵ\right)w\left(ϵ,t\right),$ (2.9)

where ${\mathcal{G}}_U$ is defined by

$\begin{array}{c}{\mathcal{G}}_U\left(U_k,{\xi }_k^H\left(ϵ\right)\right)=D_{\left(u,v,w\right)}\mathcal{G}\left(u_k,v_k,w_k,{\xi }_k^H\left(ϵ\right)\right)\left(u,v,w\right)\\ =\left(\begin{array}{c}{d}_1{u}^{″}-u_{k}u-a{u}_{k}v\\ d_2{v}^{″}+a{v}_{k}u-d{v}_{k}w\\ d_3{w}^{\prime }-{\xi }_k^H{\left(w_{k}v+{v^{\prime }}_{k}w\right)}^{\prime }+\frac{f{w}_k}{1+w_k}v+\left(\frac{f{v}_k}{{\left(1+w_k\right)}^2}-e\right)w\end{array}\right)\end{array}$

The given Frchet derivatives with respect to $U$ . By constructing eigenvalues, we can determine the stability of the branch solution near the branch point ${\xi }_k^H$ . When $ϵ=0$ , Equation (2.8) corresponds to an eigenvalue problem, where

${\mathcal{G}}_0\left(k\right)w=-\kappa \left(0\right)w$ (2.10)

${\mathcal{G}}_0\left(k\right)={\mathcal{G}}_U\left(U^{\text{*}},{\xi }_k^H\right)=\left(\begin{array}{ccc}{d}_1\frac{{\text{d}}^2}{\text{d}{x}^2}-u^{\text{*}}& -a{u}^{\text{*}}& 0\\ a{v}^{\text{*}}& d_2\frac{{\text{d}}^2}{\text{d}{x}^2}& -d{v}^{\text{*}}\\ 0& -{\xi }_k^H{w}^{\text{*}}\frac{{\text{d}}^2}{\text{d}{x}^2}& d_3\frac{{\text{d}}^2}{\text{d}{x}^2}+\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}-e\end{array}\right).$

It is immediately apparent that the spectrum of ${\mathcal{G}}_0$ is infinite-dimensional; in particular, ${\mathcal{G}}_0$ corresponds to the linearization matrix of the chemotaxis model (1.3) at $U^{\text{*}}$ :

$A_k\left({\xi }_k^H\right)=\left(\begin{array}{ccc}-d_1{\left(\frac{k\pi }{l}\right)}^2-u^{\text{*}}& -a{u}^{\text{*}}& 0\\ a{v}^{\text{*}}& -d_2{\left(\frac{k\pi }{l}\right)}^2& -d{v}^{\text{*}}\\ 0& {\xi }_k^H{w}^{\text{*}}{\left(\frac{k\pi }{l}\right)}^2& -d_3{\left(\frac{k\pi }{l}\right)}^2+\frac{f{v}^{\text{*}}}{{\left(1+w^{\text{*}}\right)}^2}-e\end{array}\right).$ (2.11)

Suppose that for some $k_0\in {\mathbb{N}}^{+}$ , $\underset{k\in {\mathbb{N}}^{+}}{\text{max}}\left\{{\xi }_k^H,{\xi }_k^S\right\}={\xi }_{k_0}^H$ . First, we prove that for any $k\ne k_0$ , the branch curve ${\omega }_k\left(ϵ\right)$ in the vicinity of ${\xi }_k^H$ is unstable. In fact, let the eigenvalues of $A_k\left({\xi }_k^H\right)$ be ${\lambda }_1^H\left({\xi }_k^H,k\right)$ ,${\lambda }_2^H\left({\xi }_k^H,k\right)$ , and ${\lambda }_3^H\left({\xi }_k^H,k\right)$ . By Theorem 2.1, if $\xi <{\xi }_0$ , then $A_k\left({\xi }_k^H\right)$ has at least one distinct eigenvalue with a positive real part. Therefore, for any positive integer $k\ne k_0$ , ${\mathcal{G}}_0\left(k\right)$ has at least one eigenvalue with a positive real part; that is, if $k\ne k_0$ , then $\kappa \left(0\right)<0$ . By fundamental perturbation theory for eigenvalues, if $k\ne k_0$ , then for sufficiently small $ϵ$ , we have $\kappa \left(ϵ\right)<0$ . Therefore, when $k\ne k_0$ , all branch curves ${\omega }_k\left(ϵ\right)$ in the vicinity of $U^{\text{*}}$ are unstable. This means that for a periodic solution to be stable, it necessarily lies on the $k_0$ -th branch, where ${\xi }_k^H$ attains its minimum for all ${\mathbb{N}}^{+}$ (i.e., it is the leftmost branch), and all branches to its right are always unstable.

Further discussion follows regarding the stability of the branch ${\omega }_{k_0}\left(ϵ\right)$ in the vicinity of $\left(U^{\text{*}},{\xi }_{k_0}^H\right)$ . According to [27], Lemma 2.10, the eigenvalue $\kappa \left(ϵ\right)$ is a continuous real-valued function of $ϵ$ in the vicinity of the origin. For $ϵ$ in the vicinity of ${\xi }_{k_0}^H$ , the eigenvalues corresponding to matrix (2.11) are ${\lambda }_1\left(\xi ,k\right)=\kappa \left(\xi ,k\right)\pm \text{i}\nu \left(\xi ,k\right)$ . From [27], Theorem 2.13, it follows that $-\kappa \left(ϵ\right)$ and $ϵ{\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)$ have the same zeros in a small neighborhood of $ϵ=0$ , $\kappa \left(ϵ\right)$ and $ϵ{\kappa }^{\prime }\left({\xi }_{k_0}^H\right){\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)$ have the same sign if $ϵ{\kappa }^{\prime }\left({\xi }_{k_0}^H\right){\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)\ne 0$ ,$\kappa \left(ϵ\right)\ne 0$ , and $\left|-\kappa \left(ϵ\right)+ϵ{\kappa }^{\prime }\left({\xi }_{k_0}^H\right){\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)\right|\le \left|{\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)\right|$ . From [27], we know that if $\kappa \left(ϵ\right)>0$ , the periodic solution is stable; if $\kappa \left(ϵ\right)<0$ , the periodic solution is unstable. Furthermore, since it has already been proven in Theorem 2.2 that ${\kappa }^{\prime }\left({\xi }_{k_0}^H\right)>0$ , $\kappa \left(ϵ\right)$ and ${\xi }_{k_0}^{H^{\prime }}\left(ϵ\right)$ have the same sign. Therefore, we need to perform a perturbation analysis in the neighborhood of the critical value ${\xi }_{k_0}^{H^{\prime }}$ , expanding ${\xi }_k^H\left(ϵ\right)$ and the periodic solution $U\left(ϵ,x,t\right)$ into Taylor series in $ϵ$ and substituting them into model (2.1). Combining this with [28], it can be shown through calculation that ${\xi }_{k_0}^H\left(0\right)={\xi }_{k_0}^H$ and ${\xi }_{k_0}^{H^{\prime }}\left(0\right)=0$ . Therefore, the branching direction and stability of the branching periodic solution are determined by the sign of ${\xi }_{k_0}^{H^{″}}\left(0\right)$ . If ${\xi }_{k_0}^{H^{″}}\left(0\right)>0$ , the branch direction is supercritical and the periodic solution is stable; if ${\xi }_{k_0}^{H^{″}}\left(0\right)<0$ , the branch direction is subcritical and the periodic solution is unstable.

Note It bears mentioning that ${\xi }_{k_0}^{H^{″}}\left(0\right)$ can be computed using the canonical form and the central manifold theorem in reference [28]. While straightforward, this calculation is rather tedious; therefore, the specific computational process is omitted here.

The scenario described by Theorem (2.3(1)) is illustrated in Figure 2 using a simple diagram.

Figure 2. Branching diagram for theorem 2.3.

In the illustration, the solid red line indicates that the branch curve is stable, while the dashed blue line indicates that the branch curve is unstable.

3. Conclusions

This manuscript investigates the Hopf bifurcation problem in a reaction-diffusion model of HTLV-I infection that incorporates CTL immune responses and chemotactic effects, and systematically analyzes the existence of Hopf bifurcations at positive equilibrium points as well as the stability of bifurcation periodic solutions.

Mathematically, this paper derives the critical conditions for the occurrence of a Hopf bifurcation, proves the existence of periodic solutions near the bifurcation point, and determines the bifurcation direction and stability of the periodic solutions: if ${\xi }_{k_0}^{H^{″}}\left(0\right)>0$ , the bifurcation direction is supercritical, the periodic solutions are stable; if ${\xi }_{k_0}^{H^{″}}\left(0\right)<0$ , the bifurcation direction is subcritical, the periodic solutions are unstable.

From a biological perspective, this result implies that the system undergoes a transition from a stable equilibrium to sustained periodic oscillations as the CTL chemotaxis intensity varies.

Author Contributions

First and corresponding author, Hongyan Kang: Conceptualization, mathematical derivation for bifurcation analysis, numerical simulation, original draft writing, interpretation of results, manuscript editing and review. Co‑author, Yinji Huang: Model nondimensionalization, equilibrium point analysis.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Gessain, A. and Cassar, O. (2012) Epidemiological Aspects and World Distribution of HTLV-1 Infection. Frontiers in Microbiology, 3, Article 388.[CrossRef] [PubMed]
[2]	Román, G. and Osame, M. (1988) Identity of HTLV-I-Associated Tropical Spastic Paraparesis and Htlv-I-Associated Myelopathy. The Lancet, 331, 651.[CrossRef] [PubMed]
[3]	Osame, M., Janssen, R., Kubota, H., Nishitani, H., Igata, A., Nagataki, S., et al. (1990) Nationwide Survey of HTLV-I—Associated Myelopathy in Japan: Association with Blood Transfusion. Annals of Neurology, 28, 50-56.[CrossRef] [PubMed]
[4]	Gallo, R.C. (2005) History of the Discoveries of the First Human Retroviruses: HTLV-1 and HTLV-2. Oncogene, 24, 5926-5930.[CrossRef] [PubMed]
[5]	Shichijo, T. and Yasunaga, J. (2025) Stratagems of HTLV-1 for Persistent Infection and the Resultant Oncogenesis: Immune Evasion and Clonal Expansion. Leukemia Research, 152, Article ID: 107680.[CrossRef] [PubMed]
[6]	Proietti, F.A., Carneiro-Proietti, A.B.F., Catalan-Soares, B.C. and Murphy, E.L. (2005) Global Epidemiology of HTLV-I Infection and Associated Diseases. Oncogene, 24, 6058-6068.[CrossRef] [PubMed]
[7]	Bangham, C.R. (2000) The Immune Response to HTLV-I. Current Opinion in Immunology, 12, 397-402.[CrossRef] [PubMed]
[8]	Bangham, C.R.M. (2003) The Immune Control and Cell-to-Cell Spread of Human T-Lymphotropic Virus Type 1. Journal of General Virology, 84, 3177-3189.[CrossRef] [PubMed]
[9]	Greten, T.F., Slansky, J.E., Kubota, R., Soldan, S.S., Jaffee, E.M., Leist, T.P., et al. (1998) Direct Visualization of Antigen-Specific T Cells: HTLV-1 Tax11-19-Specific CD8+ T Cells Are Activated in Peripheral Blood and Accumulate in Cerebrospinal Fluid from HAM/TSP Patients. Proceedings of the National Academy of Sciences of the United States of America, 95, 7568-7573.[CrossRef] [PubMed]
[10]	Wodarz, D., Nowak, M.A. and Bangham, C.R.M. (1999) The Dynamics of HTLV-I and the CTL Response. Immunology Today, 20, 220-227.[CrossRef] [PubMed]
[11]	Gomezacevedo, H. and Li, M. (2005) Backward Bifurcation in a Model for HTLV-I Infection of CD4 T Cells. Bulletin of Mathematical Biology, 67, 101-114.[CrossRef] [PubMed]
[12]	Wodarz, D., Christensen, J.P. and Thomsen, A.R. (2002) The Importance of Lytic and Nonlytic Immune Responses in Viral Infections. Trends in Immunology, 23, 194-200.[CrossRef] [PubMed]
[13]	Sun, X.G. and Wei, J.J. (2013) Global Dynamics of a HTLV-I Infection Model with CTL Response. Electronic Journal of Qualitative Theory of Differential Equations, No. 40, 1-15.[CrossRef]
[14]	Li, M.Y. and Shu, H.Y. (2012) Global Dynamics of a Mathematical Model for HTLV-I Infection of CD4+ T Cells with Delayed CTL Response. Nonlinear Analysis: Real World Applications, 13, 1080-1092.[CrossRef]
[15]	Jia, X.J. and Xu, R. (2022) Global Dynamics of a Delayed HTLV-I Infection Model with Beddington-Deangelis Incidence and Immune Impairment. Chaos, Solitons & Fractals, 155, Article ID: 111733.[CrossRef]
[16]	Chen, S.Y., Liu, Z.J., Wang, L.W. and Zhang, X. (2023) Global Dynamics Analysis for a Nonlinear HTLV-I Model with Logistic Proliferation and CTL Response. International Journal of Biomathematics, 17, Article ID: 2350023.[CrossRef]
[17]	Zhou, Y. and Li, S. (2016) Backward Bifurcation of an HTLV-I Model with Immune Response. Discrete and Continuous Dynamical Systems—Series B, 21, 863-881.[CrossRef]
[18]	Zhang, L.R. and Xu, R. (2021) Global Dynamical Behavior of a Model of HTLV-I Infection Featuring Intracellular Delay and Saturated CTL Immune Responses. Journal of Applied Mathematics in Higher Education, 36, 300-308. (In Chinese)
[19]	Xu, R. and Yang, Y. (2022) Dynamics of an HTLV-I Infection Model with Delayed and Saturated CTL Immune Response and Immune Impairment. Acta Mathematica Scientia, Series A, 42, 1836-1848.
[20]	Cui, R., Lam, K. and Lou, Y. (2017) Dynamics and Asymptotic Profiles of Steady States of an Epidemic Model in Advective Environments. Journal of Differential Equations, 263, 2343-2373.[CrossRef]
[21]	Cui, R. and Lou, Y. (2016) A Spatial SIS Model in Advective Heterogeneous Environments. Journal of Differential Equations, 261, 3305-3343.[CrossRef]
[22]	Li, H., Peng, R. and Xiang, T. (2018) Dynamics and Asymptotic Profiles of Endemic Equilibrium for Two Frequency-Dependent SIS Epidemic Models with Cross-Diffusion. European Journal of Applied Mathematics, 31, 26-56.[CrossRef]
[23]	Tao, Y. and Winkler, M. (2023) Analysis of a Chemotaxis-Sis Epidemic Model with Unbounded Infection Force. Nonlinear Analysis: Real World Applications, 71, Article ID: 103820.[CrossRef]
[24]	Tao, Y. and Winkler, M. (2023) Global Smooth Solutions in a Three-Dimensional Cross-Diffusive SIS Epidemic Model with Saturated Taxis at Large Densities. Evolution Equations and Control Theory, 12, 1676-1687.[CrossRef]
[25]	Amann, H. (1991) Hopf Bifurcation in Quasilinear Reaction-Diffusion Systems. In: Busenberg, S. and Martelli, M., Eds., Delay Differential Equations and Dynamical Systems, Springer, 53-63.[CrossRef]
[26]	Liu, P., Shi, J.P. and Wang, Z.A. (2013) Pattern Formation of the Attraction-Repulsion Keller-Segel System. Discrete and Continuous Dynamical Systems—B, 18, 2597-2625.[CrossRef]
[27]	Crandall, M.G. and Rabinowitz, P.H. (1977) The Hopf Bifurcation Theorem in Infinite Dimensions. Archive for Rational Mechanics and Analysis, 67, 53-72.[CrossRef]
[28]	Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Applications of Hopf Bifurcation. Cambridge University Press.