Open AccessArticle

Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response

Yuan Xue

Jinli Xu

and

Yuting Ding

Department of Mathematics, Northeast Forestry University, Harbin 150040, China

Author to whom correspondence should be addressed.

Axioms 2023, 12(12), 1085; https://doi.org/10.3390/axioms12121085

Submission received: 23 October 2023 / Revised: 23 November 2023 / Accepted: 24 November 2023 / Published: 27 November 2023

(This article belongs to the Special Issue Recent Advances in Applied Mathematics and Artificial Intelligence)

Download

Browse Figures

Versions Notes

Abstract

In this paper, we introduce the Crowley–Martin functional response and nonlocal competition into a reaction–diffusion immunosuppressive infection model. First, we analyze the existence and stability of the positive constant steady states of the systems with nonlocal competition and local competition, respectively. Second, we deduce the conditions for the occurrence of Turing, Hopf, and Turing–Hopf bifurcations of the system with nonlocal competition, as well as the conditions for the occurrence of Hopf bifurcations of the system with local competition. Furthermore, we employ the multiple time scales method to derive the normal forms of the Hopf bifurcations reduced on the center manifold for both systems. Finally, we conduct numerical simulations for both systems under the same parameter settings, compare the impact of nonlocal competition, and find that the nonlocal term can induce spatially inhomogeneous stable periodic solutions. We also provide corresponding biological explanations for the simulation results.

Keywords:

immunosuppressive infection model; Crowley–Martin functional response; nonlocal competition; Hopf bifurcation; delay

MSC:

37G05; 37N25; 35R10

1. Introduction

Continued infection with certain human pathogens can lead to the development of chronic diseases, such as Hepatitis B Virus (HBV), Hepatitis C Virus (HCV), and Human Immunodeficiency Virus (HIV) [1]. Although specific antiviral treatments are essential for managing these infections, the complete eradication of these pathogens remains an elusive goal due to the high genetic variability of these viruses and their ability to evade the host immune system. In recent years, a few HIV-infected individuals, known as elite controllers, have been discovered [2,3]. They can maintain low or undetectable viral loads without treatment, possibly due to their robust immune responses [4,5,6]. Enhancing virus-specific immune responses is now a key research focus.

Mathematical models have been used to model the dynamic behavior of immune suppression infections between viruses and the immune system [7,8,9]. Komarova et al. [1] proposed a basic immunosuppressive infection model in which the proliferation term of immune cells is considered to be influenced by both viruses and immune cells, expressed as

Φ (u, v) = \frac{c u v}{1 + η u}

, where c is the immune intensity,

η

stands for the inhibitory effects of the virus on the proliferation of immune cells, and u and v represent the density of viruses and immune cells at time t, respectively.

Following this, many scholars continuously modified the immunosuppressive infection model [10,11,12]. Due to the immune system requiring a time lag to generate new immune cells through a series of events [13], Shu et al. [10] proposed an assumption that the activation rate of immune cells at time t is dependent on both the virus load and the number of immune cells at time

t - τ

. They investigated the immune expansion function described as

Φ (u, v) = \frac{c u (t - τ) v (t - τ)}{1 + η u (t - τ)}

. In this way, the time delay

τ

can be seen as the period during which the immune cells undergo activation or maturation. In other words, these immune cells are produced at time

t - τ

but do not work until t. Subsequently, Li et al. [12] made further modification to the immune expansion function; they proposed that the expansion of immune response is achieved through the stimulation of the virus on the immune cells. This process involves the processing and presentation of viral antigens by antigen-presenting cells for recognition and activation of effector cells [14]. This implies that it takes a certain amount of time for the virus to effectively stimulate the immune cells, which is known as the viral stimulation delay, denoted by

τ

. In this way, the immune cells at time t are activated by the virus that stimulated them at time

t - τ

. They investigated the immune expansion function described as

Φ (u, v) = \frac{c u (t - τ) v}{1 + η u (t - τ)}

. The latest advancements in this model involve the work of Chen et al. [15], where they replaced the immune expansion function with the Beddington–DeAngelis (B-D) functional response given by:

Φ (u, v) = \frac{c u (t - τ) v}{1 + η u (t - τ) + s v},

where s represents the competition intensity among immune responses. Additionally, Xue et al. [16] introduced diffusion into the model proposed by Chen et al. [15].

In this paper, we consider replacing the Beddington–DeAngelis functional response with the Crowley–Martin (C-M) functional response:

Φ (u, v) = \frac{c u (t - τ) v}{(1 + η u (t - τ)) (1 + s v)} .

The Crowley–Martin functional response was first introduced by Crowley and Martin [17]. Like the B-D functional response, it was initially used to describe interactions and competitive relationships in predator–prey models [18,19,20]. Skalski et al. [21] analyzed an important distinction between the B-D functional response and the C-M functional response in predator–prey models through mathematical formulas; they suggested using the B-D functional response when the predator’s consumption rate is no longer correlated with predator density at high prey densities, whereas they advised employing the C-M functional response when the predator’s consumption rate continues to be influenced by higher predator densities at high prey densities.

Now, revisiting the immunosuppressive infection model, when the viral quantity is sufficiently high, the B-D functional response can be expressed mathematically as:

lim_{u \to \infty} Φ (u, v) = lim_{u \to \infty} \frac{c v}{\frac{1}{u} + η + \frac{s v}{u}} = \frac{c v}{η},

while the C-M functional response can be expressed as:

lim_{u \to \infty} Φ (u, v) = lim_{u \to \infty} \frac{c v}{\frac{1}{u} + η + \frac{s v}{u} + η s v} = \frac{c v}{η (1 + s v)} .

From this, it can be observed that as the density of viruses (denoted as u) increases, in the B-D functional response, the competition among immune cells (represented by s) has almost no effect on their proliferation rate. However, in the C-M functional response, competition among immune cells still affects their proliferation rate. In reality, as viral quantity rises, it can harm the immune system, leading to immune exhaustion. This exhaustion makes it challenging for immune cells to efficiently clear the virus, resulting in a weakened immune response (as mentioned in [22]). Even when only a small fraction of immune cells are actively responding to the viruses, some level of competition among them can influence their proliferation rate, functionality, and lifespan, thus affecting the overall immune response (as mentioned in [23]). This implies that even though the viral quantity reaches a certain threshold, competition s among immune cells still influences the proliferation term of immune cells, aligning more closely with the description of the C-M functional response. Therefore, it is reasonable to replace the B-D functional response with the C-M functional response. In fact, many scholars have incorporated the C-M functional response into virus infection models [24,25,26,27], primarily to simulate the infection rate. To the best of our knowledge, there have been relatively few instances of utilizing the C-M functional response to simulate the proliferation term of immune cells in these models.

In this paper, we will consider the following system:

\{\begin{matrix} \frac{\partial u (x, t)}{\partial t} = d_{1} Δ u (x, t) + r u (x, t) (1 - \frac{u (x, t)}{K}) - a u (x, t) - p u (x, t) v (x, t), \\ \frac{\partial v (x, t)}{\partial t} = d_{2} Δ v (x, t) + \frac{c u (x, t - τ) v (x, t)}{(1 + η u (x, t - τ)) (1 + s v (x, t))} - b v (x, t) - q u (x, t) v (x, t), \\ \partial_{ν} u (x, t) = \partial_{ν} v (x, t) = 0, x \in \partial Ω, t > 0, \\ u (x, t) = ψ_{1} (x, t), v (x, t) = ψ_{2} (x, t), x \in Ω, t \in [- τ, 0], \end{matrix}

(1)

where

Ω = (0, π)

. The variables

u (x, t)

and

v (x, t)

represent the density of viruses and immune cells at the spatial location x and time t, respectively. The parameters a and b correspond to the natural mortality rates of the viruses and immune cells, respectively. Meanwhile, r and K describe the viral replication rate and viral load, respectively. Additionally, we use the parameters c, s, and

η

to represent the immune intensity, the competition intensity among immune responses, and the inhibitory effects of the viruses on the proliferation of immune cells. The parameter

τ

characterizes the viral stimulation delay, which means the immune cells at time t are activated by the virus that stimulated them at time

t - τ

. It is postulated that immune cells eliminate viruses at a rate of

p u v

, while viruses suppress immune cells at a rate of

q u v

. The diffusion rates for the viruses and immune cells are denoted as

d_{1}

and

d_{2}

, respectively. It is crucial to emphasize that all parameters in system (1) must remain nonnegative.

Nonlocal competition was first proposed by Furter and Grinfeld in their research in 1989 [28]. They introduced nonlocal interactions into a population dynamics model of a single species to describe how species interact when competing for a shared resource. In recent years, many reaction–diffusion equations with nonlocal competition have been used in ecological system modeling [29,30,31,32]. Research indicates that introducing nonlocal effects may lead to the emergence of multi-peaked population density distributions [32,33,34]. In system (1), the viral replication rate follows logistic growth,

r u (1 - \frac{u}{K})

. That implies that there is an implicit one-to-one correspondence between viral and infected cell, while Bessonov et al. [35] and Banerjee et al. [36] mentioned that, in biomedical applications, viral cells can also compete for resources, with the possibility of different types of viral cells infecting the same cell. In this case, the viral reproduction rate is described by the expression

r u (x, t) (1 - \frac{1}{K} \int_{Ω} G (x, y) u (y, t) d y)

, where

G (x, y)

is some reasonable kernel which characterizes the efficacy of host cell infection. The integral

\hat{u} (x, t) = \int_{Ω} G (x, y) u (y, t) d y

corresponds to nonlocal resource consumption at a particular spatial point x [31]. In this paper, we choose

G (x, y) = 1 / π

. In order to investigate the influence of nonlocal competition on the dynamic behavior of system (1), we also consider the following system with nonlocal competition:

\{\begin{matrix} \frac{\partial u (x, t)}{\partial t} = d_{1} Δ u (x, t) + r u (x, t) (1 - \frac{\int_{0}^{π} u (x, t) d x}{π K}) - a u (x, t) - p u (x, t) v (x, t), \\ \frac{\partial v (x, t)}{\partial t} = d_{2} Δ v (x, t) + \frac{c u (x, t - τ) v (x, t)}{(1 + η u (x, t - τ)) (1 + s v (x, t))} - b v (x, t) - q u (x, t) v (x, t), \\ \partial_{ν} u (x, t) = \partial_{ν} v (x, t) = 0, x \in \partial Ω, t > 0, \\ u (x, t) = ψ_{1} (x, t), v (x, t) = ψ_{2} (x, t), x \in Ω, t \in [- τ, 0], \end{matrix}

(2)

where all the parameters are consistent with those in system (1).

The remaining sections of the paper are structured as follows. In Section 2, we separately discuss the stability of positive constant steady states and the existence of possible bifurcations for systems with and without nonlocal competition. Section 3 presents the normal forms of Hopf bifurcations in both systems and analyzes the stability of bifurcating periodic solutions, along with their respective bifurcation directions. In Section 4, we conduct numerical simulations for the two systems to compare the distinct spatiotemporal dynamics arising from the introduction of nonlocal competition. Finally, conclusions are given in Section 5.

2. Stability and Bifurcation

2.1. Stability and Bifurcation of the System with Nonlocal Competition

Considering the biological aspect, we focus on the positive constant steady state

E_{*} = (u_{*}, v_{*})

of system (2). Then we have:

\begin{matrix} r u_{*} (1 - \frac{u_{*}}{K}) - a u_{*} - p u_{*} v_{*} = 0, \\ \frac{c u_{*} v_{*}}{(1 + η v_{*}) (1 + s v_{*})} - b v_{*} - q u_{*} v_{*} = 0 . \end{matrix}

(3)

Suppose

(H 0) a + p v_{*} < r .

Then by direct calculation, we can deduce the form of

u_{*}

as follows:

u_{*} = K (1 - \frac{a + p v_{*}}{r}) > 0 .

(4)

Substituting (4) into (3), then

v_{*}

must satisfy the following equations:

h (v) = A v^{3} + B v^{2} + C v + D = 0,

where

\begin{matrix} A & = - q η K^{2} p^{2} s, \\ B & = - 2 p K [η [q ((a - r) s + \frac{p}{2}) K - \frac{b r s}{2}] - \frac{q r s}{2}], \\ C & = - q η (a - r) (s a - s r + 2 p) K^{2} + r [- s (b η + q) r + b (s a + p) η + (s a + p) q - c p] K - b r^{2} s, \\ D & = - q η {(a - r)}^{2} K^{2} + r (b η - c + q) (a - r) K - b r^{2} . \end{matrix}

(5)

Denote

\begin{matrix} Δ_{1} = 4 B^{2} - 12 A C, v_{1}^{+} = \frac{- 2 B + \sqrt{Δ_{1}}}{6 A}, v_{1}^{-} = \frac{- 2 B - \sqrt{Δ_{1}}}{6 A} . \end{matrix}

(6)

For the sake of discussing the existence of positive constant steady states, we make the following assumptions:

\begin{matrix} (H 1) C < 0, B > 0, Δ_{1} > 0, D > 0, h (v_{1}^{-}) < 0, h (v_{1}^{+}) > 0, \\ (H 2) C \geq 0, D < 0, h (v_{1}^{+}) > 0, \\ (H 3) C < 0, B > 0, Δ_{1} > 0, D > 0, h (v_{1}^{-}) = 0 or h (v_{1}^{+}) = 0, \\ (H 4) C < 0, B > 0, Δ_{1} > 0, D \leq 0, h (v_{1}^{+}) > 0, \\ (H 5) C \geq 0, D \geq 0 or h (v_{1}^{+}) = 0, \\ (H 6) C < 0, B > 0, Δ_{1} > 0, D > 0, h (v_{1}^{-}) > 0 or h (v_{1}^{+}) < 0, \\ (H 7) C < 0, B > 0, Δ_{1} > 0, D \leq 0, h (v_{1}^{+}) = 0, \\ (H 8) C < 0, B > 0, Δ_{1} \leq 0, D > 0, \\ (H 9) C < 0, B \leq 0, D > 0 . \end{matrix}

Theorem 1.

Assume that (H0) holds; then for the positive constant steady state

E_{*} = (u_{*}, v_{*})

of system (2), the following statements hold true.

(i): If (H1) holds, then system (2) has three positive constant steady states;
(ii): If any one of the conditions (H2)–(H4) holds, then system (2) has two positive constant steady states;
(iii): If any one of the conditions (H5)–(H9) holds, then system (2) has one positive constant steady state.

Proof.

Note that

A < 0

and

h^{'} (v) = 3 A v^{2} + 2 B v + C

is a quadratic equation. Then we will analyze the monotonicity of

h (v)

based on the direction of

h^{'} (v)

and the presence of the positive roots of

h^{'} (v)

to determine the number of positive roots of

h (v)

C < 0, B > 0 a n d Δ_{1} > 0

, then the parabola of

h^{'} (v)

opens downward and

h^{'} (v) = 0

has two positive roots

v_{1}^{+}

and

v_{1}^{-}

. That means

h (v)

is monotonically increasing when

0 < v_{1}^{-} < v < v_{1}^{+}

, and monotonically decreasing when

v \in [0, v_{1}^{-}) \cup (v^{+}, \infty)

. At this time, we have the following results.

(a): If $D > 0, h (v_{1}^{-}) < 0$ , $h (v_{1}^{+}) > 0$ are satisfied, then $h (v)$ has three positive roots (see Figure 1).
(b): If $D > 0, h (v_{1}^{-}) = 0$ or $D > 0, h (v_{1}^{-}) < 0, h (v_{1}^{+}) = 0$ or $D \leq 0$ , $h (v_{1}^{+}) > 0$ hold, then $h (v)$ has two positive roots (see Figure 2a–c).

C \geq 0

, no matter the sign of B, the parabola of

h^{'} (v)

always opens downward, and

h^{'} (v) = 0

has one positive root

v_{1}^{+}

and one negative root

v_{1}^{-}

. That means

h (v)

is monotonically increasing when

v_{1}^{-} < 0 < v < v_{1}^{+}

, and monotonically decreasing when

v \in (- \infty, v_{1}^{-}) \cup (v^{+}, \infty)

. Combining this with

D < 0, h (v_{1}^{+}) > 0

, we can conclude that

h (v)

has two positive roots (see Figure 2d).

This proves (i) and (ii). Similarly, we conclude (iii), the remaining proof process is omitted, and we only provide graphs for one positive root (see Figure 3). □

Denote the positive constant steady states of system (2) as

E_{*}

. Therefore, when the parameters of system (2) satisfy any one of conditions (i)–(iii) of Theorem 1, the positive constant steady state

E_{*} = (u_{*}, v_{*})

exists.

Set

U (x, t) = {(u (x, t), v (x, t))}^{T}

; then the linearization of system (2) at

E_{*}

\frac{\partial U (x, t)}{\partial t} = D_{*} Δ U (x, t) + J U (x, t) + A_{*} \hat{U} (x, t) + B_{*} U (x, t - τ),

(7)

where

\begin{matrix} D_{*} = (\begin{matrix} d_{1} & 0 \\ 0 & d_{2} \end{matrix}), J = (\begin{matrix} 0 & - a_{12} \\ - a_{21} & - a_{22} \end{matrix}), A_{*} = (\begin{matrix} - a_{11} & 0 \\ 0 & 0 \end{matrix}), & B_{*} = (\begin{matrix} 0 & 0 \\ b_{21} & 0 \end{matrix}), \end{matrix}

with

a_{11} = \frac{r u_{*}}{K}, a_{12} = p u_{*}, a_{21} = q v_{*}, a_{22} = \frac{c s u_{*} v_{*}}{(1 + η u *) {(1 + s v_{*})}^{2}}, b_{21} = \frac{c v_{*}}{{(1 + η u_{*})}^{2} (1 + s v_{*})} .

(8)

Clearly,

a_{11}, a_{12}, a_{21}, a_{22}

, and

b_{21}

are all positive numbers.

Therefore, the associated characteristic equations at

E_{*}

are

Γ_{n} : \{\begin{matrix} λ^{2} + T_{0} λ + N_{0} + b_{0} e^{- λ τ} = 0, n = 0, \\ λ^{2} + T_{n} λ + N_{n} + b_{0} e^{- λ τ} = 0, n = 1, 2, \dots, \end{matrix}

(9)

with

T_{0} = a_{11} + a_{22}, N_{0} = a_{11} a_{22} - a_{12} a_{21}, b_{0} = a_{12} b_{21},

and

T_{n} = a_{22} + d_{2} n^{2} + d_{1} n^{2}, N_{n} = a_{22} d_{1} n^{2} + d_{1} d_{2} n^{4} - a_{12} a_{21} .

Case 1

τ = 0

When

τ = 0

, Equation (9) becomes

\{\begin{matrix} λ^{2} + T_{0} λ + N_{0} + b_{0} = 0, n = 0, \\ λ^{2} + T_{n} λ + N_{n} + b_{0} = 0, n = 1, 2, \dots . \end{matrix}

(10)

Suppose that there exists

K_{0} > 0

such that

N_{0} + b_{0} = 0

for

K = K_{0}

; then the characteristic Equation (10) has a zero root for

n = 0

. Futhermore, note that

N_{n} + b_{0} (n = 1, 2, 3, \dots)

is monotonically increasing with respect to n, so

N_{n} + b_{0} > 0 (n = 1, 2, 3, \dots)

if and only if

N_{1} + b_{0} > 0

. Take K as a bifurcation parameter, and assume that the transversality condition for the occurrence of 0-mode Turing bifurcation is satisfied. Then we have the following results.

Theorem 2.

N_{1} + b_{0} > 0

holds, and

{(Re \frac{d λ}{d K})}^{- 1} |_{K = K_{0}} \neq 0

, then we have the following statements.

(i): If $N_{0} + b_{0} > 0$ holds, then all the roots of characteristic Equation (10) have negative real parts; thus, the positive constant steady state $E_{*}$ of system (2) is locally asymptotically stable.
(ii): If $N_{0} + b_{0} < 0$ holds, the characteristic Equation (10) will exhibit characteristic roots with positive real parts for n = 0, then the positive constant steady state $E_{*}$ of system (2) is unstable.
(iii): If $N_{0} + b_{0} = 0$ holds, the characteristic Equation (10) has a zero root for $n = 0$ , then system (2) undergoes 0-mode Turing bifurcation at $E_{*}$ .

Case 2

τ > 0

Set

λ = i ω_{n} (ω_{n} > 0)

as a root of Equation (9); then we have

\{\begin{matrix} b_{0} sin (ω_{n} τ) = T_{n} ω_{n}, \\ b_{0} cos (ω_{n} τ) = ω_{n}^{2} - N_{n}, \end{matrix} for n = 0, 1, 2, 3, \dots .

(11)

Futhermore, we define

\begin{matrix} R_{n} \overset{▵}{=} sin (ω_{n} τ) = \frac{T_{n} ω_{n}}{b_{0}}, \\ Q_{n} \overset{▵}{=} cos (ω_{n} τ) = \frac{ω_{n}^{2} - N_{n}}{b_{0}} . \end{matrix}

(12)

By squaring both equations in Equation (11) and adding them together, then denoting

g \overset{▵}{=} ω_{n}^{2}

, we can obtain

φ (g) \overset{▵}{=} g^{2} + (T_{n} - 2 N_{n}) g + N_{n}^{2} - b_{0}^{2} = 0 .

(13)

Lemma 1.

When

τ > 0

, if

N_{n}^{2} - b_{0}^{2} < 0 (n = 0, 1, 2, 3, \dots)

holds, then the characteristic Equation (9) of the linearization of system (2) has a unique pair of purely imaginary roots

\pm i ω_{n}

τ_{n, j}

, where

\begin{matrix} ω_{n} = \sqrt{\frac{2 N_{n} - T_{n}^{2} + \sqrt{Δ_{2}}}{2}}, \\ τ_{n, j} = \frac{1}{ω_{n}} (arccos (Q_{n}) + 2 j π), for n \in I_{1}, j \in N_{0} = 0, 1, 2, 3, \dots, \end{matrix}

(14)

with

Δ_{2} = {(T_{n}^{2} - 2 N_{n})}^{2} - 4 (N_{n}^{2} - b_{0}^{2}), I_{1} = \{n \in N_{0} | N_{n}^{2} - b_{0}^{2} < 0\} .

Proof.

By simple calculation, we can obtain

\{\begin{matrix} T_{0}^{2} - 4 N_{0} = {(a_{11} - a_{22})}^{2} + 4 a_{12} a_{21} > 0, n = 0, \\ T_{n}^{2} - 4 N_{n} = {(a_{22} + d_{2} n^{2} - d_{1} n^{2})}^{2} + 4 a_{12} a_{21} > 0, n = 1, 2, 3, \dots, \end{matrix}

and

\{\begin{matrix} 2 N_{0} - T_{0}^{2} = - a_{11}^{2} - a_{22}^{2} - 2 a_{12} a_{21} < 0, n = 0, \\ 2 N_{n} - T_{n}^{2} = - {(a_{22} + b_{22})}^{2} - d_{1}^{2} n^{4} - 2 a_{12} a_{21} < 0, n = 1, 2, 3, \dots . \end{matrix}

Then we conclude that, for

n = 0, 1, 2, 3, \dots

, there always exist

Δ_{2} = {(T_{n}^{2} - 2 D_{n})}^{2} - 4 (N_{n}^{2} - b_{0}^{2}) = T_{n}^{2} (T_{n}^{2} - 4 N_{n}) + 4 b_{0}^{2} > 0

and

g_{1} + g_{2} = 2 N_{n} - T_{n}^{2} < 0

. Futhermore, if

N_{n}^{2} - b_{0}^{2} < 0 (n = 0, 1, 2, 3, \dots)

holds, we can obtain

g_{1} \cdot g_{2} = N_{n}^{2} - b_{0}^{2} < 0 (n = 0, 1, 2, 3, \dots)

The results mentioned above suggest that there is a unique positive root

g = \frac{2 N_{n} - T_{n}^{2} + \sqrt{Δ_{2}}}{2}

satisfying Equation (13), then Equation (9) has a unique pair of purely imaginary roots

λ = \pm i ω_{n}

, where

ω_{n} = \sqrt{g}

. Furthermore, due to

b_{0} > 0, T_{n} > 0 a n d ω_{n} > 0

, we have

R_{n} > 0

. From Equation (12), we can derive the form of

τ_{n, j}

as shown in Equation (14). This completes the proof. □

Remark 1.

Since we can deduce that

Δ_{2} > 0

and

g_{1} + g_{2} < 0

must hold, it implies that there are no instances of two positive roots for Equation (13); thus, the positive constant steady state of system (2) will not exhibit stability switches. If

N_{n}^{2} - b_{0}^{2} \geq 0 (n = 0, 1, 2, 3, \dots)

holds, then Equation (13) has no positive root.

We can deduce from Equation (14) that for a fixed value of

n_{*} \in I_{1}

τ_{n, j}

attains its minimum at

j = 0

. Define this smallest critical value as:

τ_{c} = τ_{n_{*}, 0} = min \{τ_{n, 0} | n \in I_{1}\} .

(15)

Lemma 2.

For system (2), the transversality condition for the occurrence of

n_{*}

-mode Hopf bifurcation is satisfied, that is,

{(Re \frac{d λ}{d τ})}^{- 1} |_{τ = τ_{n, j}} > 0

, for

n \in I_{1}, j \in N_{0} .

Proof.

Differentiating both sides of Equation (9) with respect to the time delay

τ

results in

{(\frac{d λ}{d τ})}^{- 1} = \frac{- 2 λ - T_{n}}{λ [λ^{2} + T_{n} λ + N_{n}]} - \frac{τ}{λ} .

From Equation (9) we can obtain that, for all

n = 0, 1, 2, 3, \dots

, there always exists

b_{0} e^{- λ τ} = - (λ^{2} + T_{n} λ + N_{n})

. Then substituting it into the above formula, we can obtain

{(Re \frac{d λ}{d τ})}^{- 1} |_{τ = τ_{n, j}} = \frac{(T_{n}^{2} - 2 N_{n}) + 2 ω^{2}}{T_{n}^{2} ω^{2} + {(ω^{2} - N_{n})}^{2}} > 0, n \in I_{1}, j \in N_{0} .

This completes the proof. □

Theorem 3.

Assuming that

N_{n}^{2} - b_{0}^{2} < 0 (n = 1, 2, 3, \dots)

and

N_{1} + b_{0} > 0

are satisfied, for system (2), we have the following statements.

(i): When $N_{0} + b_{0} > 0$ and $N_{0}^{2} - b_{0}^{2} < 0$ , the positive constant steady state $E_{*}$ is locally asymptotically stable for $τ \in [0, τ_{c})$ and unstable for $τ \in (τ_{c}, \infty)$ . Additionally, system (2) undergoes $n_{*}$ -mode Hopf bifurcation at $E_{*}$ when $τ = τ_{c}$ ;
(ii): When $N_{0} + b_{0} < 0$ , the positive constant steady state $E_{*}$ is unstable for $τ > 0$ ;
(iii): When $N_{0} + b_{0} = 0$ , then system (2) undergoes $(0, n_{*})$ -mode Turing–Hopf bifurcation at $(K, τ) = (K_{0}, τ_{c})$ ,

where

τ_{c}

is defined in Equation (15), and

K_{0} > 0

satisfies

N_{0} + b_{0} = 0

2.2. Stability and Bifurcation of the System without Nonlocal Competition

In fact, the positive constant steady states of system (2) with nonlocal competition are identical to those of system (1) without nonlocal competition; they all align with the description provided in Theorem 1. Then, the characteristic equation of the linearization of system (1) at

E_{*}

λ^{2} + T_{n}^{'} λ + N_{n}^{'} + b_{0} e^{- λ τ} = 0, n = 0, 1, 2, \dots,

(16)

where

\begin{matrix} T_{n}^{'} = a_{22} + a_{11} + d_{1} n^{2} + d_{2} n^{2}, \\ N_{n}^{'} = a_{22} d_{1} n^{2} + d_{1} d_{2} n^{4} - a_{12} a_{21} + a_{11} d_{2} n^{2} + a_{11} a_{22}, \\ b_{0} = a_{12} b_{21}, \end{matrix}

where

a_{11}, a_{12}, a_{21}, a_{22}, b_{21}

are defined as in Equation (8).

By the same approach as in Section 2.1, we have the following results:

Theorem 4.

{(N_{1}^{'})}^{2} + b_{0}^{2} > 0

and

{(N_{n}^{'})}^{2} - b_{0}^{2} < 0

are satisfied, then the positive constant steady state of system (1) is locally asymptotically stable when

τ \in [0, τ_{c}^{'})

, and unstable when

τ \in (τ_{c}^{'}, \infty)

. Moreover, system (1) undergoes the

n_{* *}

-mode Hopf bifurcation at

E_{*}

when

τ = τ_{c}^{'}

, where

τ_{c}^{'} = τ_{n_{* *}, 0}^{'} = min \{τ_{n, 0}^{'} | n \in I_{1}^{'}\},

(17)

with

\begin{matrix} τ_{n, j}^{'} = \frac{1}{ω_{n}^{'}} (arccos (Q_{n}^{'}) + 2 j π), for n \in I_{1}^{'}, j \in N_{0}, \end{matrix}

and

\begin{matrix} R_{n}^{'} \overset{▵}{=} \frac{T_{n}^{'} ω_{n}^{'}}{b_{0}}, Q_{n}^{'} \overset{▵}{=} \frac{{(ω_{n}^{'})}^{2} - N_{n}^{'}}{b_{0}}, ω_{n}^{'} = \sqrt{\frac{2 N_{n}^{'} - {(T_{n}^{2})}^{'} + \sqrt{{(Δ_{2})}^{'}}}{2}}, \\ I_{1}^{'} = \{n \in N_{0} | {(N_{n}^{'})}^{2} - b_{0}^{2} < 0\}, Δ_{2}^{'} = {({(T_{n}^{2})}^{'} - 2 N_{n}^{'})}^{2} - 4 ({(N_{n}^{2})}^{'} - b_{0}^{2}) . \end{matrix}

The specific proof process is omitted.

3. Stability and Direction of Hopf Bifurcation

In this section, we will use the method of multiple time scales to analyze the stability and direction of Hopf bifurcation for the two systems.

Denote

f (u (x, t - τ), v (x, t)) \overset{▵}{=} \frac{c u (x, t - τ) v (x, t)}{(1 + η u (x, t - τ)) (1 + s v (x, t))},

and

\begin{matrix} β = 1 + s v_{*}, \\ α = 1 + η u_{*} . \end{matrix}

(18)

Expand

f (u (x, t - τ), v (x, t))

by using a Taylor series around

E_{*}

, leading to

\begin{matrix} f (u (x, t - τ), v (x, t)) = & \frac{c u_{*} v_{*}}{α β} + \frac{c v_{*}}{α^{2} β} (u (x, t - τ) - u_{*}) + \frac{c u_{*}}{α β^{2}} (v (x, t) - v_{*}) \\ - \frac{c η v_{*}}{α^{3} β} {(u (x, t - τ) - u_{*})}^{2} + \frac{c}{α^{2} β^{2}} (u (x, t - τ) - u_{*}) (v (x, t) - v_{*}) \\ - \frac{c s u_{*}}{α β^{3}} {(v (x, t) - v_{*})}^{2} + \frac{c η^{2} v_{*}}{α^{4} β} {(u (x, t - τ) - u_{*})}^{3} \\ - \frac{c η}{α^{3} β^{2}} {(u (x, t - τ) - u_{*})}^{2} (v (x, t) - v_{*}) \\ - \frac{c s}{α^{2} β^{3}} (u (x, t - τ) - u_{*}) {(v (x, t) - v_{*})}^{2} + \frac{c s^{2} u_{*}}{α β^{4}} {(v (x, t) - v_{*})}^{3} + \dots . \end{matrix}

(19)

3.1. Stability and Direction of Hopf Bifurcation of the System with Nonlocal Competition

We have obtained that when

τ = τ_{c}

, the characteristic equation (9) of the linearization of system (2) has a unique pair of purely imaginary roots

\pm i ω_{n_{*}}

, where

τ_{c}

is defined in (15). At this time, system (2) undergoes

n_{*}

-mode Hopf bifurcation at

E_{*}

Choose the delay

τ

as a bifurcation parameter, by setting

t \to t / τ

and translating the equilibrium points of system (2) to the origin; we can express the truncation up to the third-order form of system (2) as follows:

\{\begin{matrix} \begin{matrix} \frac{\partial u (x, t)}{\partial t} = τ [d_{1} Δ u (x, t) + r (u (x, t) + u_{*}) (1 - \frac{\int_{0}^{π} (u (x, t) + u_{*}) d x}{π K}) - a (u (x, t) + u_{*}) \\ - p (u (x, t) + u_{*}) (v (x, t) + v_{*})], \\ \frac{\partial v (x, t)}{\partial t} = τ [d_{2} Δ v (x, t) + \frac{c u_{*} v_{*}}{α β} + \frac{c v_{*}}{α^{2} β} u (x, t - 1) + \frac{c u_{*}}{α β^{2}} v (x, t) - \frac{c η v_{*}}{α^{3} β} u^{2} (x, t - 1) \\ + \frac{c}{α^{2} β^{2}} u (x, t - 1) v (x, t) - \frac{c s u_{*}}{α β^{3}} v^{2} (x, t) + \frac{c η^{2} v_{*}}{α^{4} β} u^{3} (x, t - 1) \\ - \frac{c η}{α^{3} β^{2}} u^{2} (x, t - 1) v (x, t) - \frac{c s}{α^{2} β^{3}} u (x, t - 1) v^{2} (x, t) + \frac{c s^{2} u_{*}}{α β^{4}} v^{3} (x, t) \\ - b (v (x, t) + v_{*}) - q (u (x, t) + u_{*}) (v (x, t) + v_{*})], \\ \partial_{ν} u (x, t) = \partial_{ν} v (x, t) = 0, x \in \partial Ω, t > 0, \\ u (x, t) = ψ_{1} (x, t), v (x, t) = ψ_{2} (x, t), x \in Ω = (0, l π), t \in [- τ, 0], \end{matrix} \end{matrix}

(20)

where

α

and

β

are defined as Equation (18).

Denote the eigenvector of the characteristic matrix of system (20) at

E_{*}

corresponding to the eigenvalue

i ω_{n_{*}} τ

h = {(h_{11}, h_{12})}^{T}

, and the eigenvector of the adjoint matrix corresponding to the eigenvalue

- i ω_{n_{*}} τ

h^{*} = d {(h_{13}, h_{14})}^{T}

; they satisfy the inner product

〈h^{*}, h〉 = {\bar{h}}^{* T} \cdot h = 1 .

Then we have

\begin{matrix} h_{11} = 1, h_{12} = - \frac{1}{p u_{*}} (i ω_{n_{*}} + n_{*}^{2} d_{1}), \\ h_{14} = 1, h_{13} = - \frac{1}{p u_{*}} [- i ω_{n_{*}} + \frac{c s u_{*} v_{*}}{α β^{2}} + n_{*}^{2} d_{2}], d = {(h_{13} {\bar{h}}_{11} + h_{14} {\bar{h}}_{12})}^{- 1} . \end{matrix}

(21)

Set

τ = τ_{c} + ε μ

, where

τ_{c}

represents the critical point for the Hopf bifurcation, as defined in Equation (15). Here,

μ

serves as a small perturbation parameter, and

ε

is the characteristic scale used for non-dimensionalization. Next, we only present the main results regarding the computation of the normal form for Hopf bifurcation. The specific computational process is detailed in Appendix A.

By the method of multiple time scales, we deduce the normal form of the

n_{*}

-mode Hopf bifurcation for system (2) reduced on the center manifold as follows:

\dot{G} = M μ G + χ G^{2} \bar{G},

(22)

where M and

χ

are given in Equation (A11) and Equation (A17) in Appendix A.

Now introduce

G = r e^{i θ}

and substitute it into the normal form Equation (22). This yields:

\{\begin{matrix} \dot{r} = Re (M) μ r + Re (χ) r^{3}, \\ \dot{θ} = Im (M) μ + Im (χ) r^{2} . \end{matrix}

(23)

Regarding the stability of the bifurcating periodic solutions of system (2), the following result can be established:

Theorem 5.

When

\frac{Re (M) μ}{Re (χ)} < 0

, system (2) has periodic solutions near the positive constant steady state

E_{*}

. Moreover:

(i): If $Re (M) μ < 0$ , the bifurcating periodic solutions reduced on the center manifold are unstable, and when $μ > 0$ ( $μ < 0$ ), the direction of bifurcation is forward (backward);
(ii): If $Re (M) μ > 0$ , the bifurcating periodic solutions reduced on the center manifold are stable, and when $μ > 0$ ( $μ < 0$ ), the direction of bifurcation is forward (backward).

3.2. Stability and Bifurcation of the System without Nonlocal Competition

Similar to the analysis in reference [16], we can also obtain the normal form of the

n_{* *}

-mode Hopf bifurcation for system (1) reduced on the center manifold:

\{\begin{matrix} \dot{r} = Re (M^{'}) μ r + Re (χ^{'}) r^{3}, \\ \dot{θ} = Im (M^{'}) μ + Im (χ^{'}) r^{2}, \end{matrix}

(24)

where

\begin{matrix} M^{'} = \frac{{\bar{h^{'}}}_{13} M_{1}^{'} + {\bar{h}}_{14} M_{2}^{'}}{{\bar{h^{'}}}_{13} + {\bar{h}}_{14} (h_{12}^{'} + \frac{τ_{c}^{'} c v_{*}}{α^{2} β} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}})}, χ^{'} = \frac{R^{'}}{V^{'}}, \end{matrix}

(25)

with

\begin{matrix} M_{1}^{'} = - a + r - n_{* *}^{2} d_{1} - p v_{*} - p u_{*} h_{12}^{'} - \frac{2 r u_{*}}{K}, \\ M_{2}^{'} = - n_{* *}^{2} d_{2} h_{12} + \frac{c z_{*} e^{- i ω_{n_{* *}}^{'} τ_{c}}}{α^{2} β} + (\frac{c u_{*}}{α^{2} β} - b - q u_{*}) h_{12} - q z_{*}, \\ V^{'} = [{\bar{h^{'}}}_{13} + {\bar{h}}_{14} (h_{12}^{'} + \frac{τ_{c}^{'} c v_{*} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{2} β})] \int_{0}^{π} {cos}^{2} (n_{* *} x) d x, \\ R^{'} = {\bar{h}}_{13}^{'} [- \frac{2 τ_{c} r \sum_{k = 0}^{+ \infty} (η_{0, k}^{'} + η_{1, k}^{'})}{K} - τ_{c}^{'} p \sum_{k = 0}^{+ \infty} (ζ_{0, k}^{'} + ζ_{1, k}^{'} + h_{12}^{'} η_{0, k}^{'} + {\bar{h}}_{12}^{'} η_{1, k}^{'})] \cdot E_{k}^{'} \\ + {\bar{h}}_{14} [\frac{3 τ_{c}^{'} c η^{2} v_{*} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{4} β} - \frac{τ_{c}^{'} c η e^{- 2 i ω_{n_{* *}}^{'} τ_{c}^{'}} {\bar{h}}_{12}^{'}}{α^{3} β^{2}} - \frac{τ_{c}^{'} c s {(h_{12}^{'})}^{2} e^{i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{2} β^{3}} \\ + \frac{3 τ_{c} c s^{2} u_{*} {(h_{12}^{'})}^{2} {\bar{h}}_{12}^{'}}{α β^{4}}] \cdot \int_{0}^{π} {cos}^{4} (n_{* *} x) d x + {\bar{h}}_{14} [\frac{- 2 τ_{c}^{'} c η v_{*} \sum_{k = 0}^{+ \infty} (η_{0, k}^{'} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}} + η_{1, k}^{'} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}})}{α^{3} β} \\ + \frac{τ_{c}^{'} c \sum_{k = 0}^{+ \infty} (ζ_{0, k}^{'} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}} + ζ_{1, k}^{'} e^{i ω_{n_{* *}}^{'} τ_{c}^{'}} + h_{12}^{'} η_{0, k}^{'} + {\bar{h}}_{12}^{'} η_{1, k}^{'} e^{- 2 i ω_{n_{* *}}^{'} τ_{c}^{'}})}{α^{2} β^{2}} \\ - \frac{- 2 τ_{c}^{'} c u_{*} s \sum_{k = 0}^{+ \infty} (h_{12}^{'} ζ_{0, k}^{'} + {\bar{h}}_{12}^{'} ζ_{1, k}^{'})}{α β^{3}} - τ_{c}^{'} q \sum_{k = 0}^{+ \infty} (ζ_{0, k}^{'} + ζ_{1, k}^{'} + h_{12}^{'} η_{0, k}^{'} + {\bar{h}}_{12}^{'} η_{1, k})] \cdot E_{k}^{'}, \end{matrix}

and

\begin{matrix} h_{12}^{'} = - \frac{1}{p u_{*}} [i ω_{n}^{'} + n_{* *}^{2} d_{1} + \frac{r u_{*}}{K}], \\ h_{13}^{'} = - \frac{1}{p u_{*}} [- i ω_{n}^{'} + \frac{c s u_{*} v_{*}}{α β^{2}} + n_{* *}^{2} d_{2}], \\ η_{1, k}^{'} = \frac{(J_{2} + 2 i ω_{n_{* *}}^{'}) l_{1}^{'} - W^{'}}{\frac{(J_{2} + 2 i ω_{n_{* *}}^{'}) (J_{1}^{'} + 2 i ω_{n_{* *}}^{'})}{p u_{*}} - (- \frac{c v_{*} e^{- 2 i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{2} β} + q v_{*})} \cdot \frac{E_{k}^{'}}{C_{k}}, \\ η_{0, k}^{'} = \frac{J_{2} l_{2}^{'} - S^{'}}{\frac{J_{1}^{'} J_{2}}{p u_{*}} - (- \frac{c v_{*}}{α^{2} β} + q v_{*})} \cdot \frac{E_{k}^{'}}{C_{k}}, \\ ζ_{1, k}^{'} = l_{1}^{'} \cdot \frac{E_{k}}{C_{k}} - \frac{J_{1}^{'} + 2 i ω_{n_{* *}}^{'}}{p u_{*}} \cdot η_{1, k}^{'}, \\ ζ_{0, k}^{'} = l_{2}^{'} \cdot \frac{E_{k}}{C_{k}} - \frac{J_{1}^{'}}{p u_{*}} \cdot η_{0, k}^{'}, \\ W^{'} = \frac{- c η v_{*} e^{- 2 i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{3} β} + \frac{c h_{12}^{'} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}}}{α^{2} β^{2}} - \frac{c u_{*} s {(h_{12}^{'})}^{2}}{α β^{3}} - q h_{12}^{'}, \\ S^{'} = \frac{- 2 c η v_{*}}{α^{3} β} + \frac{c (h_{12}^{'} e^{i ω_{n_{* *}}^{'} τ_{c}^{'}} + {\bar{h^{'}}}_{12} e^{- i ω_{n_{* *}}^{'} τ_{c}^{'}})}{α^{2} β^{2}} - \frac{2 c u_{*} s h_{12}^{'} {\bar{h^{'}}}_{12}}{α β^{3}} - q (h_{12}^{'} + {\bar{h^{'}}}_{12}), \\ J_{1}^{'} = p v_{*} + a - r + \frac{2 r u_{*}}{K} + k^{2} d_{1}, \\ l_{1}^{'} = \frac{1}{p u_{*}} \cdot [\frac{- r}{K} - p h_{12}^{'}], \\ l_{2}^{'} = \frac{1}{p u_{*}} \cdot [- p (h_{12}^{'} + {\bar{h^{'}}}_{12}) - \frac{2 r}{K}], \\ E_{k}^{'} = \int_{0}^{π} {cos}^{2} (n_{* *} x) cos (k x) d x, \\ F_{k}^{'} = \int_{0}^{π} cos (n_{* *} x) cos (k x) d x \cdot \int_{0}^{π} cos (k x) d x . \end{matrix}

The remaining parameters

h_{14}, J_{2}, α, β, C_{k}, ω_{n_{* *}}^{'}

can be found in Appendix A and Equation (17).

4. Numerical Simulations

In this section, we will conduct numerical simulations to demonstrate the theoretical findings we have obtained in the preceding sections. For both system (2) and system (1), we select the same set of parameters as follows:

\begin{matrix} a = 3.0 {day}^{- 1}, p = 1.0 {mm}^{3} {cells}^{- 1} {day}^{- 1}, q = 1.0 {mm}^{3} {virus}^{- 1} {day}^{- 1}, \\ d_{1} = 0.4 {mm}^{- 3} {day}^{- 1}, K = 10 virus {mm}^{- 3}, b = 2 {day}^{- 1}, \\ d_{2} = 0.1 {mm}^{- 3} {day}^{- 1}, c = 4 {mm}^{3} {virus}^{- 1} {day}^{- 1}, r = 7.0 {day}^{- 1}, \\ η = 0.35 {mm}^{- 1} virus, s = 0.03 {mm}^{3} immune {cells}^{- 1} . \end{matrix}

(26)

4.1. Numerical Simulations of the System with Nonlocal Competition

With the parameters given in (26), we can obtain a positive constant steady state for system (2):

E_{*} = (1.337684786, 3.063620649) .

According to Lemma 1 and (15), we can obtain that

n_{*} = 1

ω_{1} = 1.7003

I_{1} = \{n = 0, 1, 2, 3\}

. And the critical point for the Hopf bifurcation is

τ_{c} = τ_{1, 0} = 0.1128

. Furthermore, according to (A11) and (A17), we can calculate that

M = 0.3695 + 1.6296 i, χ = - 0.0039 - 0.0186 i

. Hence, when

τ \in [0, τ_{c}) = [0, 0.1128)

, we can deduce from Theorem 3 that

E_{*}

is locally asymptotically stable, as shown in Figure 4; when

τ \in (0.1128, + \infty)

, we can infer from Theorems 3 and 5 that

E_{*}

is unstable and system (2) has stable forward bifurcating periodic solutions near

τ = τ_{1, 0} = 0.1128

, as shown in Figure 5. When

τ

significantly exceeds the critical value

τ_{c}

, system (2) no longer exhibits stable periodic solutions and will undergo substantial oscillations, as shown in Figure 6.

Remark 2.

From Figure 4, it can be observed that when the time delay τ required for immune cells to generate an immune response is less than the critical value

τ_{c}

, system (2) undergoes a period of oscillations before eventually converging to a stable state. This implies that when the immune system is stimulated by the virus, it takes some time to mount an immune response, but eventually manages to control the number of viruses within a certain range. In this scenario, the disease is brought under control, and the quantities of viruses and immune cells no longer exhibit significant fluctuations.

Remark 3.

From Figure 5, it can be observed that when the delay τ required for immune cells to generate an immune response is slightly greater than the critical value

τ_{c}

, system (2) loses its original stable state and generates spatially inhomogeneous periodic solutions. This occurrence encapsulates the intricate dynamics between viruses and the immune system. When viruses trigger an immune response, immune cells swiftly localize and eliminate local viruses, leading to a decline in the local virus population. However, the nonlocal effect implies that the immune system has a certain sensing range and can perceive the density of viruses in distant areas. This results in the immune system generating a stronger immune response when it senses areas with a high density of viruses in other locations, to combat viruses at a distance. Consequently, the dynamism in viruses and immune cells across the system becomes non-uniform, exhibiting periodic spatial fluctuations. This results in varying dominance; in some regions, the immune system prevails, effectively clearing viruses, while in others, viruses regain dominance. These spatially driven periodic oscillations stem from the intricate interplay between the immune system and viruses, underscoring the complex dynamics of their interaction. The detailed mathematical analysis emphasizes the crucial role played by nonlocal competition in the emergence of spatial inhomogeneity, highlighting its impact on the system’s dynamics.

We have also observed a stable number of cells in the central space, while the cell count on both sides exhibits periodic fluctuations. This phenomenon may be related to the boundary conditions. When immune cells localize to areas with higher virus density, they congregate at the boundary due to the inability to traverse it. As the virus density decreases at the boundary, the nonlocal competition suggests that the immune cells might hypothetically move in the opposite direction, resulting in periodic variations in the number of immune cells on both sides of the boundary. With nearly equivalent counts of immune cells in both forward and backward diffusion trends, the number of immune cells in the central space tends to stabilize.

Remark 4.

Once τ significantly exceeds the threshold, taking

τ = 0.12 > 0.1128

as an example, system (2) no longer possesses stable periodic solutions, and the virus proliferates rapidly (see Figure 6). This is the worst case scenario, where the virus could spread rapidly, posing significant risks to the host.

4.2. Numerical Simulations of the System without Nonlocal Competition

Based on the parameters in (26), the positive constant steady states in system (1) equal the positive constant steady state

E_{*}

in system (2), indicating that

E_{*} = (1.337684786, 3.063620649) .

Similarly, by calculation, we find

n_{* *} = 0, ω_{0}^{'} = 1.68, I_{1}^{'} = \{n = 0, 1, 2, 3\}, τ_{c}^{'} = τ_{0, 0}^{'} = 0.1774,

M^{'} = 0.5345 + 1.5215 i, χ^{'} = - 0.0008 - 0.1513 i

. Therefore,

E_{*}

is locally asymptotically stable when

τ \in [0, τ_{c}^{'}) = [0, 0.1774)

, as shown in Figure 7, and unstable when

τ \in (0.1774, + \infty)

. Furthermore, system (1) also has stable bifurcating periodic solutions and the direction of the Hopf bifurcation is also forward, as shown in Figure 8.

Remark 5.

Under the identical parameter settings, a noticeable difference is observed between the reaction–diffusion system (1) without nonlocal competition and system (2) incorporating nonlocal competition terms. Specifically, in the former case, the critical value

τ_{c}^{'}

is comparatively larger, whereas, in contrast, system (2), with nonlocal competition, exhibits a smaller critical value. Furthermore, Figure 8 demonstrates that system (1) also exhibits periodic and steady state solutions. However, a pivotal disparity emerges regarding the nature of these solutions.

In the context of system (1), the solutions appear spatially homogeneous, indicating an independence of the interaction between the immune system and viruses from spatial variables. This suggests that in the absence of nonlocal competition, the dynamics of the immune response to viruses remain consistent throughout space.

However, a more realistic scenario warrants consideration, suggesting that the interaction between the immune system and viruses should be influenced by spatial variables. In systems where nonlocal competition terms are present, the smaller critical value implies that spatial influences significantly impact the dynamics of the immune system–virus interaction. This signifies a spatial dependency, wherein the behavior of the immune system in response to viruses may vary across different spatial regions, emphasizing the importance of spatial considerations in understanding the dynamics of virus–immune system interactions.

5. Conclusions

This paper introduced the Crowley–Martin functional response into an immunosuppressive infection model and conducted a comparative analysis of the dynamic behavior between a delayed reaction–diffusion system with nonlocal viral competition and a delayed reaction–diffusion system without nonlocal viral competition. The nonlocal term represents the space-time weighted average of the viral population. We derived the normal forms of the Hopf bifurcations reduced on the center manifold for both systems, one with nonlocal competition and the other without, and performed numerical simulations. The results indicated that the introduction of nonlocal competition causes a reduction in the stability region, i.e.,

τ_{c} < τ_{c}^{'}

, and the system with nonlocal competition generated spatially inhomogeneous periodic solutions, whereas the reaction–diffusion system without nonlocal competition terms yielded spatially homogeneous solutions. This suggested that systems with nonlocal effects exhibited more diverse dynamic behaviors, highlighting the special biological significance of spatiotemporal nonlocality.

Author Contributions

Writing—original draft preparation, Y.X.; funding acquisition, J.X.; methodology, J.X.; supervision, Y.D. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by Fundamental Research Funds for the Central Universities of China (Grant No. 2572022DJ07).

Data Availability Statement

No new data were created in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

We assume that the solution of system (20) follows the following form:

U (x, t) = U (x, T_{0}, T_{1}, T_{2}, \dots) = \sum_{k = 1}^{+ \infty} ε^{k} U_{k} (x, T_{0}, T_{1}, T_{2}, \dots),

(A1)

with

\begin{matrix} U (x, T_{0}, T_{1}, T_{2}, \dots) = {(u (x, T_{0}, T_{1}, T_{2}, \dots), v (x, T_{0}, T_{1}, T_{2}, \dots))}^{T}, \\ U_{k} (x, T_{0}, T_{1}, T_{2}, \dots) = {(u_{k} (x, T_{0}, T_{1}, T_{2}, \dots), v_{k} (x, T_{0}, T_{1}, T_{2}, \dots))}^{T} . \end{matrix}

(A2)

The time derivative is now transformed into

\frac{\partial}{\partial t} = \frac{\partial}{\partial T_{0}} + ε \frac{\partial}{\partial T_{1}} + ε^{2} \frac{\partial}{\partial T_{2}} + \dots = D_{0} + ε D_{1} + ε^{2} D_{2} + \dots,

(A3)

with

D_{i} = \frac{\partial}{\partial T_{i}}, i \in N_{0} = 0, 1, 2, 3, \dots .

Set

\begin{matrix} u_{j} = u_{j} (x, T_{0}, T_{1}, T_{2}, \dots), u_{j, 1} = u_{j} (x, T_{0} - 1, T_{1}, T_{2}, \dots), \\ v_{j} = v_{j} (x, T_{0}, T_{1}, T_{2}, \dots), v_{j, 1} = v_{j} (x, T_{0} - 1, T_{1}, T_{2}, \dots), \end{matrix}

(A4)

where

j \in N = 1, 2, 3, \dots

. Then we have

\frac{\partial U}{\partial t} = ε D_{0} U_{1} + ε^{2} D_{1} U_{1} + ε^{3} D_{2} U_{1} + ε^{2} D_{0} U_{2} + ε^{3} D_{1} U_{2} + ε^{3} D_{0} U_{3} + \dots .

(A5)

To handle the delayed terms, we can approximate

u (x, t - 1)

u (x, T_{0} - 1, T_{1}, T_{2}, \dots)

; thus, we obtain:

u (x, t - 1) = ε u_{11} + ε^{2} u_{21} + ε^{3} u_{31} - ε^{2} D_{1} u_{11} - ε^{3} D_{1} u_{21} - ε^{3} D_{2} u_{11} + \dots,

(A6)

with

u_{j 1} = u_{j} (x, T_{0} - 1, T_{1}, T_{2}, \dots), j \in N

Substituting Equations (A1)–(A6) into system (20) results in a set of perturbation equations. Extract the coefficients of the

ε

terms:

\{\begin{matrix} D_{0} u_{1} - τ_{c} [d_{1} Δ u_{1} - p u_{1} v_{*} + r u_{1} - p u_{*} v_{1} - a u_{1} - \frac{r u_{1} u_{*}}{K} - \frac{r u_{*} \int_{0}^{π} u_{1} d x}{π K}] = 0, \\ D_{0} v_{1} - τ_{c} [d_{2} Δ v_{1} - b v_{1} - q u_{*} v_{1} - q u_{1} v_{*} + \frac{c v_{*} u_{11}}{α^{2} β} + \frac{c u_{*} v_{1}}{α β^{2}}] = 0, \end{matrix}

(A7)

where

α, β

are given by Equation (18). Hence, the solution to Equation (A7) can be represented in the following form:

\{\begin{matrix} u_{1} = G (T_{1}, T_{2}, \dots) e^{i ω_{n_{*}} τ_{c} T_{0}} h_{11} cos (n_{*} x) + c . c ., \\ v_{1} = G (T_{1}, T_{2}, \dots) e^{i ω_{n_{*}} τ_{c} T_{0}} h_{12} cos (n_{*} x) + c . c ., \end{matrix}

(A8)

where

h_{11}, h_{12}

are defined by Equation (21), and

c . c .

denotes the complex conjugate of the preceding term. Next, collecting the coefficient of

ε^{2}

, we can obtain

\{\begin{matrix} D_{0} u_{2} + & τ_{c} (p v_{*} + a - r + \frac{r u_{*}}{K}) u_{2} + τ_{c} p u_{*} v_{2} - τ_{c} d_{1} Δ u_{2} + \frac{τ_{c} r u_{*} \int_{0}^{π} u_{2} d x}{π K} \\ = & - D_{1} u_{1} - μ d_{1} Δ u_{1} - τ_{c} p u_{1} v_{1} - μ p u_{*} v_{1} + (- μ p v_{*} - μ a + μ r - \frac{μ r y_{*}}{K}) u_{1} \\ - \frac{τ_{c} r u_{1} \int_{0}^{π} u_{1} d x}{π K} - \frac{μ r u_{*} \int_{0}^{π} u_{1} d x}{π K}, \\ D_{0} v_{2} + & τ_{c} [- d_{2} Δ v_{2} + b v_{2} + q u_{*} v_{2} + q u_{2} v_{*} - \frac{c v_{*} u_{21}}{α^{2} β} - \frac{c u_{*} v_{2}}{α β^{2}}] \\ = & τ_{c} [- q u_{1} v_{1} - \frac{c v_{*} D_{1} u_{11}}{α^{2} β} - \frac{c η v_{*} u_{11}^{2}}{α^{3} β} + \frac{c u_{11} v_{1}}{α^{2} β^{2}} - \frac{c u_{*} s v_{1}^{2}}{α β^{3}}] \\ + μ [d_{2} Δ v_{1} - b v_{1} - q u_{*} v_{1} - q u_{1} v_{*} + \frac{c v_{*} u_{11}}{α^{2} β} + \frac{c u_{*} v_{1}}{α β^{2}}] - D_{1} v_{1}, \end{matrix}

(A9)

where

α

and

β

are determined by Equation (18). After substituting solution (A8) into the right side of Equation (A9) and simplifying, we define the coefficient vector of

e^{i ω_{n_{*}} τ_{c} T_{0}}

m_{1}

. It is important to note that

m_{1}

must satisfy the solvability condition:

〈h^{*}, (m_{1}, cos (n_{*} x))〉 = 0 .

Then one has

\frac{\partial G}{\partial T_{1}} = M μ G,

(A10)

where

M = \frac{{\bar{h}}_{13} M_{1} + {\bar{h}}_{14} M_{2} + {\bar{h}}_{13} M_{3}}{{\bar{h}}_{13} + {\bar{h}}_{14} (h_{12} + \frac{τ_{c} c v_{*}}{α^{2} β} e^{- i ω_{n_{*}} τ_{c}})},

(A11)

with

\begin{matrix} M_{1} = - n_{*}^{2} d_{1} - p v_{*} - a + r - p u_{*} h_{12} - \frac{r u_{*}}{K}, \\ M_{2} = - n_{*}^{2} d_{2} h_{12} + \frac{c z_{*} e^{- i ω_{n_{*}} τ_{c}}}{α^{2} β} + (\frac{c u_{*}}{α^{2} β} - b - q u_{*}) h_{12} - q z_{*}, \\ M_{3} = (- \frac{r u_{*}}{π K}) \cdot \frac{\int_{0}^{π} c o s (n_{*} x) d x \cdot \int_{0}^{π} c o s (n_{*} x) d x}{\int_{0}^{π} c o s^{2} (n_{*} x) d x} . \end{matrix}

Assume that the solution of Equation (A9) takes the form

\{\begin{matrix} u_{2} = \sum_{k = 0}^{+ \infty} (η_{0, k} G \bar{G} + η_{1, k} e^{2 i ω_{n_{*}} τ_{c} T_{0}} G^{2} + {\bar{η}}_{1, k} e^{- 2 i ω_{n_{*}} τ_{c} T_{0}} {\bar{G}}^{2}) cos (k x), \\ v_{2} = \sum_{k = 0}^{+ \infty} (ζ_{0, k} G \bar{G} + ζ_{1, k} e^{2 i ω_{n_{*}} τ_{c} T_{0}} G^{2} + {\bar{ζ}}_{1, k} e^{- 2 i ω_{n_{*}} τ_{c} T_{0}} {\bar{G}}^{2}) cos (k x) . \end{matrix}

(A12)

To facilitate the determination of the coefficients for the second-order solution (A12), we denote

\begin{matrix} C_{k} = \int_{0}^{π} cos (k x) cos (k x) d x, G_{k} = \int_{0}^{π} cos (k x) d x \cdot \int_{0}^{π} cos (k x) d x, \\ E_{k} = \int_{0}^{π} {cos}^{2} (n_{*} x) cos (k x) d x, F_{k} = \int_{0}^{π} cos (n_{*} x) cos (k x) d x \cdot \int_{0}^{π} cos (k x) d x . \end{matrix}

(A13)

After substituting Equations (A8) and (A12) into Equation (A9) and simplifying, we take the inner product with

cos (k x)

on both sides of the resulting equation, where

k \in N_{0}

. By summing up the obtained results, we can then compare the coefficients in front of the terms

G^{2}

and

G \bar{G}

, leading to

\begin{matrix} η_{1, k} = \frac{(J_{2} + 2 i ω_{n_{*}}) l_{1} - W \cdot \frac{E_{k}}{C_{k}}}{\frac{(J_{2} + 2 i ω_{n_{*}}) (J_{1} + 2 i ω_{n_{*}})}{p u_{*}} - (- \frac{c v_{*} e^{- 2 i ω_{n_{*}} τ_{c}}}{α^{2} β} + q v_{*})}, ζ_{1, k} = l_{1} - \frac{J_{1} + 2 i ω_{n_{*}}}{p u_{*}} \cdot η_{1, k}, \\ η_{0, k} = \frac{J_{2} l_{2} - S \cdot \frac{E_{k}}{C_{k}}}{\frac{J_{1} J_{2}}{p u_{*}} - (- \frac{c v_{*}}{α^{2} β} + q v_{*})}, ζ_{0, k} = l_{2} - \frac{J_{1}}{p u_{*}} \cdot η_{0, k}, \end{matrix}

(A14)

where

\begin{matrix} W = & \frac{- c η v_{*} e^{- 2 i ω_{n_{*}} τ_{c}}}{α^{3} β} + \frac{c h_{12} e^{- i ω_{n_{*}} τ_{c}}}{α^{2} β^{2}} - \frac{c u_{*} s h_{12}^{2}}{α β^{3}} - q h_{12}, \\ S = & \frac{- 2 c η v_{*}}{α^{3} β} + \frac{c (h_{12} e^{i ω_{n_{*}} τ_{c}} + {\bar{h}}_{12} e^{- i ω_{n_{*}} τ_{c}})}{α^{2} β^{2}} - \frac{2 c u_{*} s h_{12} {\bar{h}}_{12}}{α β^{3}} - q (h_{12} + {\bar{h}}_{12}), \\ J_{1} = & p v_{*} + a - r + \frac{r u_{*}}{K} + k^{2} d_{1} + \frac{r u_{*}}{π K} \cdot \frac{G_{k}}{C_{k}}, \\ J_{2} = & k^{2} d_{2} - \frac{c u_{*}}{α β^{2}} + b + q u_{*}, \\ l_{1} = & \frac{1}{p u_{*}} \cdot [- p h_{12} \cdot \frac{E_{k}}{C_{k}} - \frac{r}{π K} \cdot \frac{F_{k}}{C_{k}}], \\ l_{2} = & \frac{1}{p u_{*}} \cdot [- p (h_{12} + {\bar{h}}_{12}) \cdot \frac{E_{k}}{C_{k}} - \frac{r (h_{11} + {\bar{h}}_{11})}{π K} \cdot \frac{F_{k}}{C_{k}}], \\ β = & 1 + s v_{*}, \\ α = & 1 + η u_{*} . \end{matrix}

By calculating, we discover that for

n_{*} \in N_{0}

η_{1, k}, η_{0, k}, ζ_{1, k}, ζ_{0, k} \neq 0

if and only if

k = 0

. Therefore, we can deduce the specific form of solution (A12):

\{\begin{matrix} u_{2} = η_{0, 0} G \bar{G} + η_{1, 0} e^{2 i ω_{n_{*}} τ_{c} T_{0}} G^{2} + {\bar{η}}_{1, 0} e^{- 2 i ω_{n_{*}} τ_{c} T_{0}} {\bar{G}}^{2}, \\ v_{2} = ζ_{0, 0} G \bar{G} + ζ_{1, 0} e^{2 i ω_{n_{*}} τ_{c} T_{0}} G^{2} + {\bar{ζ}}_{1, 0} e^{- 2 i ω_{n_{*}} τ_{c} T_{0}} {\bar{G}}^{2} . \end{matrix}

(A15)

Collecting the coefficient of

ε^{3}

, we can obtain

\{\begin{matrix} D_{0} u_{3} & - τ_{c} d_{1} Δ u_{3} + τ_{c} (a - r + \frac{r u_{*}}{K} + p z_{*}) u_{3} + τ_{c} p u_{*} v_{3} + \frac{τ_{c} r u_{*} \int_{0}^{π} u_{3} d x}{π K} \\ = & - D_{2} u_{1} - D_{1} u_{2} + μ d_{1} Δ u_{2} - μ (a - r + \frac{r u_{*}}{K} + p v_{*}) u_{2} - μ p u_{1} v_{1} - μ p u_{*} v_{2} - τ_{c} p u_{1} v_{2} \\ - τ_{c} p u_{2} v_{1} - \frac{μ r u_{1} \int_{0}^{π} u_{1} d x}{π K} - \frac{τ_{c} r u_{2} \int_{0}^{π} u_{1} d x}{π K} - \frac{μ r y_{*} \int_{0}^{π} u_{2} d x}{π K} - \frac{τ_{c} r u_{1} \int_{0}^{π} u_{2} d x}{π K}, \\ D_{0} v_{3} - & τ_{c} d_{2} Δ v_{3} - \frac{τ_{c} c z_{*} u_{31}}{α^{2} β} + τ_{c} (\frac{- c v_{*}}{α β^{2}} + b + q u_{*}) v_{3} + τ_{c} q v_{*} u_{3} \\ = & - D_{2} v_{1} - D_{1} v_{2} + τ_{c} [\frac{c v_{*} (- D_{1} u_{21} - D_{2} u_{11})}{α^{2} β} - \frac{2 c η v_{*} u_{11} (- D_{1} u_{11} + u_{21})}{α^{3} β} + \frac{c u_{11} v_{2}}{α^{2} β^{2}} \\ + \frac{c (- D_{1} u_{11} + u_{21}) v_{1}}{α^{2} β^{2}} - \frac{2 c u_{*} s v_{1} v_{2}}{α β^{3}} + \frac{c η^{2} v_{*} u_{11}^{3}}{α^{4} β} - \frac{c η y_{11}^{2} v_{1}}{α^{3} β^{2}} - \frac{c s u_{11} v_{1}^{2}}{α^{2} β^{3}} + \frac{c s^{2} u_{*} v_{1}^{3}}{α β^{4}} \\ - q u_{1} v_{2} - q u_{2} v_{1}] + μ [d_{2} Δ v_{2} + \frac{c z_{*} (- D_{1} u_{11} + u_{21})}{α^{2} β} + (\frac{c y_{*}}{α β^{2}} - b - q u_{*}) v_{2} \\ - \frac{c η v_{*} u_{11}^{2}}{α^{3} β} + \frac{c u_{11} v_{1}}{α^{2} β^{2}} - \frac{c y_{*} s v_{1}^{2}}{α β^{3}} - q u_{1} v_{1} - q u_{2} v_{*}], \end{matrix}

(A16)

where

α

and

β

are determined by Equation (18). Substitute solutions (A8) and (A12) into the right side of Equation (A16) and simplify the equation. Next, we define the coefficient vector associated with

e^{i ω_{n_{*}} τ_{c} T_{0}}

m_{2}

. It is worth noting that

m_{2}

satisfies the solvability condition presented as

〈h^{*}, (m_{2}, cos (n_{*} x))〉 = 0 .

Given that

μ

is regarded as the disturbance parameter, and

μ^{2}

has a negligible impact on the outcome, we can ignore the terms containing

μ^{2} G

. Then we obtain

\frac{\partial G}{\partial T_{2}} = χ G^{2} \bar{G},

where

χ = \frac{R}{V},

(A17)

with

\begin{matrix} β = 1 + s z_{*}, \\ α = 1 + η y_{*}, \\ V = [{\bar{h}}_{13} + {\bar{h}}_{14} (h_{12} + \frac{τ_{c} c v_{*} e^{- i ω_{n_{*}} τ_{c}}}{α^{2} β})] \int_{0}^{π} {cos}^{2} (n_{*} x) d x, \\ R = {\bar{h}}_{13} [- τ_{c} p \sum_{k = 0}^{+ \infty} (ζ_{0, k} + ζ_{1, k} + h_{12} η_{0, k} + {\bar{h}}_{12} η_{1, k})] \cdot E_{k} + {\bar{h}}_{13} [- \frac{τ_{c} r \sum_{k = 0}^{+ \infty} (η_{0, k} + η_{1, k})}{π K}] \cdot F_{k} \\ + {\bar{h}}_{13} [- \frac{τ_{c} r \sum_{k = 0}^{+ \infty} (η_{0, k} + η_{1, k})}{π K}] \cdot \int_{0}^{π} {cos}^{2} (n_{*} x) d x \cdot \int_{0}^{π} cos (k x) d x \\ + {\bar{h}}_{14} [\frac{3 τ_{c} c η^{2} v_{*} e^{- i ω_{n_{*}} τ_{c}}}{α^{4} β} - \frac{τ_{c} c η e^{- 2 i ω_{n_{*}} τ_{c}} {\bar{h}}_{12}}{α^{3} β^{2}} - \frac{τ_{c} c s h_{12}^{2} e^{i ω_{n_{*}} τ_{c}}}{α^{2} β^{3}} + \frac{3 τ_{c} c s^{2} u_{*} h_{12}^{2} {\bar{h}}_{12}}{α β^{4}}] \cdot \int_{0}^{π} {cos}^{4} (n_{*} x) d x, \\ + {\bar{h}}_{14} [\frac{- 2 τ_{c} c η v_{*} \sum_{k = 0}^{+ \infty} (η_{0, k} e^{- i ω_{n_{*}} τ_{c}} + η_{1, k} e^{- i ω_{n_{*}} τ_{c}})}{α^{3} β} - \frac{- 2 τ_{c} c u_{*} s \sum_{k = 0}^{+ \infty} (h_{12} ζ_{0, k} + {\bar{h}}_{12} ζ_{1, k})}{α β^{3}} \\ + \frac{τ_{c} c \sum_{k = 0}^{+ \infty} (ζ_{0, k} e^{- i ω_{n_{*}} τ_{c}} + ζ_{1, k} e^{i ω_{n_{*}} τ_{c}} + h_{12} η_{0, k} + {\bar{h}}_{12} η_{1, k} e^{- 2 i ω_{n_{*}} τ_{c}})}{α^{2} β^{2}} \\ - τ_{c} q \sum_{k = 0}^{+ \infty} (ζ_{0, k} + ζ_{1, k} + h_{12} η_{0, k} + {\bar{h}}_{12} η_{1, k})] \cdot E_{k} . \end{matrix}

Thus we deduce the normal form of the

n_{*}

-mode Hopf bifurcation for system (2) reduced on the center manifold:

\dot{G} = M μ G + χ G^{2} \bar{G} .

References

Komarova, N.L.; Barnes, E.; Klenerman, P.; Wodarz, D. Boosting immunity by antiviral drug therapy: A simple relationship among timing, efficacy, and success. Proc. Natl. Acad. Sci. USA 2003, 100, 1855–1860. [Google Scholar] [CrossRef]
Sheppard, H.W.; Celum, C.; Michael, N.L.; O’Brien, S.; Dean, M.; Carrington, M.; Dondero, D.; Buchbinder, S.P. HIV-1 infection in individuals with the CCR5-Δ32/32/Δ32 genotype: Acquisition of syncytium-inducing virus at seroconversion. J. Acquir. Immune Defic. Syndr. 2002, 29, 307–313. [Google Scholar] [CrossRef]
Walker, B.D. Elite control of HIV Infection: Implications for vaccines and treatment. Top. HIV Med. 2007, 15, 134–136. [Google Scholar]
Hersperger, A.R.; Martin, J.N.; Shin, L.Y.; Sheth, P.M.; Kovacs, C.M.; Cosma, G.L.; Betts, M.R. Increased HIV-specific CD8+ T-cell cytotoxic potential in HIV elite controllers is associated with T-bet expression. Blood 2011, 117, 3799–3808. [Google Scholar] [CrossRef]
Saag, M.; Deeks, S.G. How do HIV elite controllers do what they do? Clin. Infect. Dis. 2010, 51, 239–241. [Google Scholar] [CrossRef]
Deeks, S.G.; Walker, B.D. Human immunodeficiency virus controllers: Mechanisms of durable virus control in the absence of antiretroviral therapy. Immunity 2007, 27, 406–416. [Google Scholar] [CrossRef]
Funk, G.A.; Jansen, V.A.; Bonhoeffer, S.; Killingback, T. Spatial models of virus-immune dynamics. J. Theor. Biol. 2005, 233, 221–236. [Google Scholar] [CrossRef]
Browne, C.J.; Smith, H.L. Dynamics of virus and immune response in multi-epitope network. J. Math. Biol. 2018, 77, 1833–1870. [Google Scholar] [CrossRef]
Wodarz, D.; Nowak, M.A. Mathematical models of HIV pathogenesis and treatment. Bioessays 2002, 24, 1178–1187. [Google Scholar] [CrossRef]
Shu, H.; Wang, L.; Watmough, J. Sustained and transient oscillations and chaos induced by delayed antiviral immune response in an immunosuppressive infection model. J. Math. Biol. 2014, 68, 477–503. [Google Scholar] [CrossRef]
Tian, C.; Gan, W.; Zhu, P. Stability analysis in a diffusional immunosuppressive infection model with delayed antiviral immune response. Math. Methods Appl. Sci. 2017, 40, 4001–4013. [Google Scholar] [CrossRef]
Li, J.; Ma, X.; Li, J.; Zhang, D. Dynamics of a chronic virus infection model with viral stimulation delay. Appl. Math. Lett. 2021, 122, 107547. [Google Scholar] [CrossRef]
Fenton, A.; Lello, J.; Bonsall, M.B. Pathogen responses to host immunity: The impact of time delays and memory on the evolution of virulence. Proc. R. Soc. B 2006, 273, 2083–2090. [Google Scholar] [CrossRef]
Johansen, P.; Storni, T.; Rettig, L.; Manolova, V.; Lang, K.S.; Qiu, Z. Antigen kinetics determines immune reactivity. Swiss Med. Wkly. 2007, 137, 23S. [Google Scholar] [CrossRef]
Chen, Y.; Wang, L.; Liu, Z.; Wang, Y. Complex dynamics for an immunosuppressive infection model with virus stimulation delay and nonlinear immune expansion. Qual. Theory Dyn. Syst. 2023, 22, 118. [Google Scholar] [CrossRef]
Xue, Y.; Xu, J.; Ding, Y. Dynamics analysis of a diffusional immunosuppressive infection model with Beddington-DeAngelis functional response. Electron. Res. Arch. 2023, 31, 6071–6088. [Google Scholar] [CrossRef]
Crowley, P.H.; Martin, E.K. Functional responses and interference within and between year classes of a dragonfly population. J. N. Am. Benthol. Soc. 1989, 8, 211–221. [Google Scholar] [CrossRef]
Hossain, S.; Haque, M.M.; Kabir, M.H.; Gani, M.O.; Sarwardi, S. Complex spatiotemporal dynamics of a harvested prey-predator model with Crowley-Martin response function. Results Control Optim. 2021, 5, 100059. [Google Scholar] [CrossRef]
Chen, S.; Wei, J.; Yu, J. Stationary patterns of a diffusive predator-prey model with Crowley-Martin functional response. Nonlinear Anal. Real World Appl. 2018, 39, 33–57. [Google Scholar] [CrossRef]
Upadhyay, R.K.; Raw, S.N.; Rai, V. Dynamical complexities in a tri-trophic hybrid food chain model with Holling type II and Crowley-Martin functional responses. Nonlinear Anal. Model. Control. 2010, 15, 361–375. [Google Scholar] [CrossRef]
Skalski, G.T.; Gilliam, J.F. Functional responses with predator interference: Viable alternatives to the holling type II model. Ecology 2001, 82, 3083–3092. [Google Scholar] [CrossRef]
Bocharov, G.A. Modelling the dynamics of LCMV infection in mice: Conventional and exhaustive CTL responses. J. Theor. Biol. 1998, 192, 283–308. [Google Scholar] [CrossRef]
Keşmir, C.; De Boer, R.J. Clonal exhaustion as a result of immune deviation. Bull. Math. Biol. 2003, 65, 359–374. [Google Scholar] [CrossRef]
Naik, P.A.; Zu, J.; Ghoreishi, M. Stability analysis and approximate solution of SIR epidemic model with Crowley-Martin type functional response and Holling type-II treatment rate by using homotopy analysis method. J. Appl. Anal. Comput. 2020, 10, 1482–1515. [Google Scholar] [CrossRef]
Jan, M.N.; Ali, N.; Zaman, G.; Ahmad, I.; Shah, Z.; Kumam, P. HIV-1 infection dynamics and optimal control with Crowley-Martin function response. Comput. Methods Programs Biomed. 2020, 193, 105503. [Google Scholar] [CrossRef]
Roomi, V.; Afshari, H.; Gharahasanlou, T.K. Global Stability of an HIV Dynamical Model with Crowley-Martin Functional Response. Lett. Nonlinear Anal. Appl. 2023, 1, 39–46. [Google Scholar]
Li, X.; Fu, S. Global stability of a virus dynamics model with intracellular delay and CTL immune response. Math. Methods Appl. Sci. 2015, 38, 420–430. [Google Scholar] [CrossRef]
Furter, J.; Grinfeld, M. Local vs. non-local interactions in population dynamics. J. Math. Biol. 1989, 27, 65–80. [Google Scholar] [CrossRef]
Geng, D.; Jiang, W.; Lou, Y.; Wang, H. Spatiotemporal patterns in a diffusive predator-prey system with nonlocal intraspecific prey competition. Stud. Appl. Math. 2022, 148, 396–432. [Google Scholar] [CrossRef]
Wu, S.; Song, Y. Spatiotemporal dynamics of a diffusive predator-prey model with nonlocal effect and delay. Commun. Nonlinear Sci. Numer. Simul. 2020, 89, 105310. [Google Scholar] [CrossRef]
Liu, Y.; Duan, D.; Niu, B. Spatiotemporal dynamics in a diffusive predator-prey model with group defense and nonlocal competition. Appl. Math. Lett. 2020, 103, 106175. [Google Scholar] [CrossRef]
Doebeli, M.; Killingback, T. Metapopulation dynamics with quasi-local competition. Theor. Popul. Biol. 2003, 64, 397–416. [Google Scholar] [CrossRef] [PubMed]
Yang, R.; Nie, C.; Jin, D. Spatiotemporal dynamics induced by nonlocal competition in a diffusive predator-prey system with habitat complexity. Nonlinear Dyn. 2022, 110, 879–900. [Google Scholar] [CrossRef]
Yang, R.; Wang, F.; Jin, D. Spatially inhomogeneous bifurcating periodic solutions induced by nonlocal competition in a predator–prey system with additional food. Math. Meth. Appl. Sci. 2022, 45, 9967–9978. [Google Scholar] [CrossRef]
Bessonov, N.; Bocharov, G.; Meyerhans, A.; Popov, V.; Volpert, V. Nonlocal reaction-diffusion model of viral evolution: Emergence of virus strains. Math 2020, 8, 117. [Google Scholar] [CrossRef]
Banerjee, M.; Kuznetsov, M.; Udovenko, O.; Volpert, V. Nonlocal Reaction-Diffusion Equations in Biomedical Applications. Acta Biotheor. 2022, 70, 12. [Google Scholar] [CrossRef]

Figure 1. If

(H 1)

holds, then system (2) has three positive constant steady states.

Figure 1. If

(H 1)

holds, then system (2) has three positive constant steady states.

Figure 2. System (2) has two positive constant steady states under various conditions. (a,b) Correspond to when

(H 3)

holds. (c) Corresponds to when

(H 4)

holds. (d) Corresponds to when

(H 2)

holds.

Figure 2. System (2) has two positive constant steady states under various conditions. (a,b) Correspond to when

(H 3)

holds. (c) Corresponds to when

(H 4)

holds. (d) Corresponds to when

(H 2)

holds.

Figure 3. System (2) has one positive constant steady state under various conditions. (a,b) Correspond to when

(H 5)

holds. (c,d) Correspond to when

(H 6)

holds. (e) Corresponds to when

(H 7)

holds. (f) Left corresponds to when

(H 8)

holds, and right corresponds to when

(H 9)

holds.

Figure 3. System (2) has one positive constant steady state under various conditions. (a,b) Correspond to when

(H 5)

holds. (c,d) Correspond to when

(H 6)

holds. (e) Corresponds to when

(H 7)

holds. (f) Left corresponds to when

(H 8)

holds, and right corresponds to when

(H 9)

holds.

Figure 4. Choose

τ = 0.1 < 0.1128

and the initial function is (

1.337684786 + 0.01 cos x, 3.063620649 + 0.025 cos x

). Then the positive constant steady state

E_{*} = (1.337684786, 3.063620649)

of system (2) is locally asymptotically stable.

Figure 4. Choose

τ = 0.1 < 0.1128

and the initial function is (

1.337684786 + 0.01 cos x, 3.063620649 + 0.025 cos x

). Then the positive constant steady state

E_{*} = (1.337684786, 3.063620649)

of system (2) is locally asymptotically stable.

Figure 5. Choose

τ = 0.1165 > 0.1128

and consider the initial function as (

1.337684786 + 0.01 cos x

3.063620649 + 0.025 cos x

). Then system (2) exhibits a spatially inhomogeneous stable periodic solution.

Figure 5. Choose

τ = 0.1165 > 0.1128

and consider the initial function as (

1.337684786 + 0.01 cos x

3.063620649 + 0.025 cos x

). Then system (2) exhibits a spatially inhomogeneous stable periodic solution.

Figure 6. Choose

τ = 0.12 > 0.1128

and consider the initial function as (

1.337684786 + 0.01 cos x

3.063620649 + 0.025 cos x

). The constant steady state solution

E_{*}

of system (2) is unstable.

Figure 6. Choose

τ = 0.12 > 0.1128

and consider the initial function as (

1.337684786 + 0.01 cos x

3.063620649 + 0.025 cos x

). The constant steady state solution

E_{*}

of system (2) is unstable.

Figure 7. Choose

τ = 0.165 < 0.1774

and the initial function is (

1.337684786 + 0.003 + 0.005 cos x

3.063620649 + 0.001 cos x

). Then the constant steady state

E_{*} = (1.337684786, 3.063620649)

of system (1) is locally asymptotically stable.

Figure 7. Choose

τ = 0.165 < 0.1774

and the initial function is (

1.337684786 + 0.003 + 0.005 cos x

3.063620649 + 0.001 cos x

). Then the constant steady state

E_{*} = (1.337684786, 3.063620649)

of system (1) is locally asymptotically stable.

Figure 8. Choose

τ = 0.178 > 0.1774

, and consider the initial function as (

1.337684786 + 0.003 + 0.005 cos x, 3.063620649 + 0.001 cos x

), then system (1) exhibits a spatially homogeneous stable periodic solution.

Figure 8. Choose

τ = 0.178 > 0.1774

, and consider the initial function as (

1.337684786 + 0.003 + 0.005 cos x, 3.063620649 + 0.001 cos x

), then system (1) exhibits a spatially homogeneous stable periodic solution.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Xue, Y.; Xu, J.; Ding, Y. Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response. Axioms 2023, 12, 1085. https://doi.org/10.3390/axioms12121085

AMA Style

Xue Y, Xu J, Ding Y. Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response. Axioms. 2023; 12(12):1085. https://doi.org/10.3390/axioms12121085

Chicago/Turabian Style

Xue, Yuan, Jinli Xu, and Yuting Ding. 2023. "Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response" Axioms 12, no. 12: 1085. https://doi.org/10.3390/axioms12121085

APA Style

Xue, Y., Xu, J., & Ding, Y. (2023). Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response. Axioms, 12(12), 1085. https://doi.org/10.3390/axioms12121085

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Spatiotemporal Dynamics of a Diffusive Immunosuppressive Infection Model with Nonlocal Competition and Crowley–Martin Functional Response

Abstract

1. Introduction

2. Stability and Bifurcation

2.1. Stability and Bifurcation of the System with Nonlocal Competition

2.2. Stability and Bifurcation of the System without Nonlocal Competition

3. Stability and Direction of Hopf Bifurcation

3.1. Stability and Direction of Hopf Bifurcation of the System with Nonlocal Competition

3.2. Stability and Bifurcation of the System without Nonlocal Competition

4. Numerical Simulations

4.1. Numerical Simulations of the System with Nonlocal Competition

4.2. Numerical Simulations of the System without Nonlocal Competition

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI