We consider a pursuit-evasion differential game problem in which countably many pursuers chase one evader in the Hilbert space ℓ2 and for a fixed period of time. Dynamic of each of the pursuer is governed by first order differential equations and that of the evader by a second order differential equation. The control function for each of the player satisfies an integral constraint. The distance between the evader and the closest pursuer at the stoppage time of the game is the payoff of the game. The goal of the pursuers is to minimize the distance to the evader and that of the evader is the opposite. We constructed optimal strategies of the players and find value of the game.
1.
Introduction
Epidemics have posed a significant threat to global public health over the years. The emergence of COVID-19 in 2019 has had a profound impact on human health, the global economy, and social behavior. Nevertheless, the effective addressing of disease transmission remains a challenge. Mathematical modeling has emerged as a crucial tool in tackling this challenge. Numerous disease models have been developed in existing literature to study and control the spread of epidemics. It is important to note that mathematical models based on ordinary differential equations (i.e., classical derivatives) have their limitations and may not accurately capture biological phenomena. On the other hand, fractional models can offer a relatively more accurate understanding of disease outbreaks. Therefore, they are increasingly being used to simulate disease transmission with higher accuracy. Relevant literature on fractional models can be found in [1,2,3,4,5]. In addition, as another epidemic that endangers human health, the acute viral infection caused by the human immunodeficiency virus (HIV) studied in this paper is also a hot issue in society. Many researchers aim to capture the dynamics between viral and antiviral immune responses through mathematical modeling.
As we all know, HIV, which causes AIDS, can directly infect the immune system (mainly regulating CD4+ T cells). The consequences of this impact are a continuous decrease in the number of CD4+ T cells, ultimately leading to the death of infected individuals due to immune system collapse. The highly active antiretroviral therapy (HAART) currently in widespread use has been shown to improve the survival probability of HIV patients and reduce the incidence rate [6,7]. This therapy effectively suppresses the plasma virus to levels below the standard detection for extended periods and even halts viral evolution [8,9]. However, many complex problems arise after long-term use, such as obvious drug resistance, and, due to the side effects of antiviral drugs, some AIDS patients have poor compliance with antiviral therapy [10,11,12].
In fact, numerous mathematical models have been put forth to describe the dynamics of HIV and elucidate various phenomena. The effect of antiviral therapy has been investigated by some researchers [13,14,15,16]. In their work, Xiao et al.[13] analyzed the free terminal time optimal tracking control problem to determine the optimal multidrug therapy for HIV, considering both the optimal time frame and therapeutic strategies. The literature [17,18,19,20,21,22,23,24] incorporated the expansion delay of immune cells to discuss the local and global stability of equilibrium solutions. In particular, [17] indicates that such an unstable equilibrium will exhibit oscillatory solutions of increasing amplitude. In recent years, realizing gradually the multiple effects of spatial heterogeneity and mobility, many scholars utilized the reaction-diffusion equation to investigate the spatial effect of viral infection [23,24,25].
Several recent clinical studies have exhibited that structured treatment interruptions (STIs) can be used for early treatment of HIV infection to achieve sustained specific immunity. For some chronically infected individuals who may require lifelong medication, this may be a beneficial option as it can help patients rebuild their immune system during periods of non-medication[26]. While numerous mathematical models have been employed for simulating continuous therapy [27,28,29], there is scant investigation on modeling structured interruptions in treatment.
To investigate strategies for STIs, Tang et al.[30] suggested a piecewise model for delineating CD4 cell-guided STIs. The system provides an explanation for some controversial clinical research results. In 2017, Tang et al.[31] proposed a mathematical model to describe the dynamics of the interplay between the virus and the immune system. This model takes into consideration the structured treatment guided by effector cells while also incorporating the use of combined antiretroviral therapy and interleukin (IL)-2 treatment. However, they posit a linear growth pattern for the HIV virus [31], which does not accurately reflect the true dynamics of the virus. Some clinical facts show that the growth of HIV may have a saturation effect[32].
To better illustrate the non-linear evolutionary characteristics of the interplay between virus and immune response, the immunosuppressive infection model was devised by Komarova et al.[33]. The model is given as follows:
The assumptions in model (1.1) are as follows:
❑ y and z represent the population sizes of the virus and immune cells, respectively. The virus population is assumed to grow logistically: The replication rate at low viral loads, denoted as r, is expected to decrease linearly with an increase in the viral load until it becomes zero at a viral load K.
❑ c represents the immune intensity, while the proliferation term is denoted as cyz/(1+hy). Thus, the assessment of immune cell proliferation depends on both the immune cells and the virus. The inhibitory effect of the virus on the proliferation of immune cells is represented by the variable η.
❑ The viral elimination rate, denoted as a, is a result of natural decay and antiretroviral therapy. Immune cells, which have the ability to kill the virus at a rate pyz, also have a death rate b. Furthermore, these immune cells can be inhibited by the virus at a rate qyz.
The study conducted in [33] aimed to investigate the optimal timing and duration of antiviral treatment. The research elucidates the presence of bistability dynamics, wherein a stable state without immunity coexists alongside a stable state with immunity. Meanwhile, this model demonstrates the attainment of sustained immunity following the interruption of therapy. Additionally, Wang et al.[34] expanded on this model and uncovered that bistability arises within the range delineated by the post-treatment control threshold and the elite control threshold.
Following the pioneer works above[30,31,32,33,34], in this thesis, we extend model (1.1) by proposing a Filippov immunosuppressive infection model with viral logistic growth and effector cell-guided therapy. We have proposed the following model:
We assume that the sole course of action is to administer antiretroviral therapy to the patient if the number of effector cells exceeds the critical value ET. Conversely, when the count falls below the ET threshold, a combination of antiretroviral therapy and immune therapy is concurrently implemented. In this context, ε represents the rate at which effector cells grow as a result of immune therapy, like the treatment of interleukin (IL)-2. As the interpretation of other parameters is consistent with a model (1.1), all parameters in (1.2) remain nonnegative.
This paper presents a switching model with viral load logistic growth to analyze effector cell-guided treatment and assess the threshold strategy's effectiveness. The following section provides an overview of the model, defining the switching system, and summarizing the dynamic behavior of the subsystems. Additionally, in Section 3, there is a discussion on sliding mode and dynamics, exploring the presence of a sliding domain and pseudo-equilibrium. The global dynamics of the proposed model are examined in Section 4, while Section 5 focuses on the boundary equilibrium bifurcation of the system. Finally, the paper concludes with discussions and biological implications.
2.
Filippov model and preliminaries
2.1. Model formulation
By rearranging the system (1.2), we can obtain a generic planar system in the form of Filippov given by
with
Systems (2.1) and (2.2) describe a Filippov immunosuppressive infection model where (2.1) is considered a free system when ϕ=0(i.e.,z>ET), indicating that the patient receives antiretroviral therapy. On the other hand, (2.1) as a control system when ϕ=1(i.e.,z<ET) reflects the simultaneous utilization of antiretroviral therapy and immune therapy.
Let R2+={X=(y,z)T|y⩾0,z⩾0}, S1={X∈R2+|H(X)>0}, and S2={X∈R2+|H(X)<0} with H(X) being a smooth scale function. For convenience, we further denote
We can rewrite model (1.2) to represent the Filippov system as follows:
The discontinuous boundary Σ that separates the two areas can be represented as:
It is evident that R2+=S1∪Σ∪S2. Henceforth, we shall designate the Filippov system (2.4) as subsystem S1 when it is defined within region S1, and as subsystem S2 when defined within region S2.
Let
where ⟨⋅,⋅⟩ represents the standard scalar product and HX(X) denotes the gradient of H(X) that remains nonvanishing on Σ. FSiH(X)=FSi⋅gradH(X) is the Lie derivative of H with respect to the vector field FSi (i=1,2) at X. To analyze the direction of the vector field [FS1(X),FS2(X)], through a specific point X∈Σ, we categorize the areas on Σ based on whether the vector field points towards it:
(a) Crossing region:
(b) Sliding region:
(c) Escaping region:
Throughout the paper, it is crucial to have a clear understanding of the definitions regarding all types of equilibria in Filippov systems [35,36].
Definition 2.1. If FS1(X∗)=0,H(X∗)>0, or FS2(X∗)=0, H(X∗)<0, then X∗ is defined as a real equilibrium of the Filippov system (2.4). Analogously, if FS1(X∗)=0, H(X∗)<0, or FS2(X∗)=0,H(X∗)>0, then X∗ is a virtual equilibrium. Both the real and virtual equilibriums are named as regular equilibria.
Definition 2.2. If X∗ is an equilibrium of the sliding mode of system (2.4), and satisfies (1−λ)FS1(X∗)+λFS2(X∗)=0,H(X∗)=0 with 0<λ<1, then X∗ is a pseudo-equilibrium, where
Definition 2.3. If FS1(X∗)=0,H(X∗)=0, or FS2(X∗)=0,H(X∗)=0, then X∗ is defined a boundary equilibrium of Filippov system (2.4).
Definition 2.4. If FS1H(X∗)=0 but F2S1H(X∗)>0(F2S1H(X∗)<0), then X∗ is a visible (invisible) Σ-fold point of FS1. The same definition applies to FS2.
Definition 2.5. If X∗∈Σs and FS1H(X∗)=0 or FS2H(X∗)=0, then X∗ is defined a tangent point of Filippov system (2.4).
2.2. Qualitative analysis of subsystems
The model equation for the free system S1 is as follows:
We can easily define a threshold R0=ra. If R0<1, subsystem S1 has only one uninfected equilibrium E10=(0,0); if R0>1, then subsystem S1 also has an immune-free equilibrium E11=(y1,0)=(K(1−ar),0).
We certainly get an equation for y
It can be confirmed that S1(y)=0 has a sole solution when c=q+bη±2√bqη. Denote c1=q+bη−2√bqη and c2=q+bη+2√bqη. Thus, we have two possible positive roots
when c>c2, where B=c−q−bη. Substituting y11 or y12 into the first equation of (2.11), we get
Let R1i=ra(1−y1iK)(i=1,2). Hence, the subsystem S1 has two immune equilibriums E11=(y11,z11) and E12=(y12,z12) when R1i>1 and c>c2 are satisfied.
In fact, subsystem S1 has been extensively examined in a previous study. Therefore, we will provide an overview of the main findings without delving into specific calculations. For more details on the discussion of the stability analysis of this system, please consult [34]. Here, we define three thresholds by referring [34], i.e., c∗=q+bη+2qηK(1−ar), c∗∗=q+bη+bK(1−ar)+qηK(1−ar), and threshold Rc1=1+r√bqηaqηK. Moreover, we have the following
Lemma 2.1. If R0<1, the infection-free equilibrium E10 is globally asymptotically stable (GAS). If R0>1, we have E11 as locally asymptotically stable (LAS) when 0<c<c2 or c2<c<c∗∗, and E11 is unstable when c>c∗∗. The immune equilibrium E11 is LAS when Rc1>R0>1 and c>c∗∗ or R0>Rc1 and c2<c. Suppose R0>Rc1>1 and c2<c<c∗, then positive equilibrium E11 is LAS and E12 is a unstable saddle.
Remark 2.1. For R0>Rc1>1 and c2<c<c∗∗, subsystem S1 has bistable behavior, i.e., E11 and E11 are bistable but the equilibrium E12 is an unstable saddle. In other cases, subsystem S1 does not exhibit bistable behavior. Note that the post-treatment control threshold is represented by c2, while the elite control threshold is denoted as c∗∗. The range between c2 and c∗∗ is referred to as the bistable interval.
Dynamical analysis of the subsystem S1 is presented in Table 1.
Control system S2 gives
The basic reproduction number is also R0=ra for subsystem S2. Similarly, we can get uninfected equilibrium E20=(0,0) and the immune-free equilibrium E21=(K(1−ar),0). Thus, we use E10=E20=E0 and E21=E11=E1 in the following.
Using the same method as subsystem S1, we can derive the following quadratic equation for y
Let c3=q+(b−ε)η−2√qη(b−ε) and c4=q+(b−ε)η+2√qη(b−ε), then there are two possible positive roots
when c>c4, where D=c−q−bη+εη. It follows that
In fact, we can obtain y21<y11<y12<y22 by doing a simple calculation. Now, we define R2i=ra(1−y2iK)(i=1,2). For subsystem S2, we have two positive equilibriums E21=(y21,z21) and E22=(y22,z22) when c>c4, b>ε and R2i>1(i=1,2) hold. Furthermore, following a similar approach to subsystem S1, we can also define thresholds c†=q+(b−ε)η+2qηK(1−ar), c††=q+(b−ε)η+b−εK(1−ar)+qηK(1−ar), and Rc2=1+r√qη(b−ε)aqηK. Meanwhile, we have the following results about the behaviors of subsystem S2 by using the consistent method with S1.
Lemma 2.2. Suppose R0<1, the equilibrium E20 is GAS. Suppose R0>1, then E21 is LAS when 0<c<c4 or c4<c<c††, and E21 is unstable when c>c††. If Rc2>R0>1 and c>c†† or R0>Rc2 and c4<c, the subsystem S2 has a locally asymptotically stable immune equilibrium E21; if R0>Rc2>1 and c4<c<c†, then positive equilibrium E21 is LAS and E22 is an unstable saddle.
Dynamical analysis of the subsystem S2 is shown in Table 2.
3.
Sliding mode dynamics
In this section, we will provide a brief overview of the definitions pertaining to the sliding segment and crossing segment discussed in Section 2. We have σ(X)=⟨HX(X),FS1(X)⟩⋅⟨HX(X),FS2(X)⟩. Here, HX(X)=(∂H∂y,∂H∂z) is the non-vanishing gradient on the discontinuity boundary Σ, where H=z−ET. Therefore, we denote
and calculating the inequality σ(X)<0 obtains y21<y<y11 and y12<y<y22. Naturally, we can verify that there are ⟨FS1,HX(X)⟩=cyz1+ηy−bz−qyz<0 and ⟨FS2,HX(X)⟩=cyz1+ηy−bz−qyz+εz>0 for y21<y<y11 or y12<y<y22. Therefore, the Filippov system (2.4) always comprises two sliding segments, which can be obtained as
Naturally, the crossing region we can get is Σc={(y,z)∈R2+|0<y<y21, or y11<y<y12, or y>y22,z=ET}. Notably, every trajectory within the segment {(y,z)∈R2+|0<y<y21 or y>y22,z=ET} will intersect the z=ET line, moving from region S1 to S2. Similarly, the trajectory within the segment {(y,z)∈R2+|y11<y<y12,z=ET} will cross the z=ET line, transitioning from region S2 to S1.
The Filippov convex method is employed in this study to analyze the sliding domain and sliding mode dynamics of the switching system (2.4). According to Definition 2.2, we have
By a straightforward calculation, one can get
It follows that
Therefore, the dynamic equation of the switching system (2.4) on the sliding mode domain is
There exists one positive equilibrium Ec=(yc,ET), where yc=K(1−a+pETr). We can easily obtain that r>a+pET. According to Definition 2.2, the equilibrium Ec is referred to as a pseudo-equilibrium.
When y21<yc<y11, we have
and
Clearly, y21<yc<y11 is equivalent to
(i.e.,z11<ET<z21). Thus, the equilibrium Ec is located on the sliding segment Σ1s when z11<ET<z21. It can be readily confirmed that cyz/(1+hy)<0 holds on the segment Σ2s. Similariy, if y12<yc<y22, we conclude that
(i.e.,z22<ET<z12). Under this circumstance, the pseudo-equilibrium Ec is located on the sliding segment Σ2s={(y,z)|y12<y<y22,z=ET} for z22<ET<z12, meanwhile, we can get cyz/(1+hy)>0 holds on the segment Σ1s.
Theorem 3.1. If z11<ET<z21, system (2.4) has one pseudo-equilibrium Ec=(K(1−a+pETr),ET) on the sliding segment Σ1s, which is always LAS when it exists; analogously, if z22<ET<z12, then pseudo-equilibrium Ec is LAS on the sliding segment Σ2s.
Proof. Without loss of generality, we only consider the proof for the first case. Let y′=ry(1−yK)−ay−pyz.=g(y). Substituting H=z−ET=0 into the g(y), we will get g(y)=−rKy2+(r−a−pET)y. Based on the function g(y), we know g(y21)>0,g(y11)<0. Further, we note that y′=g(y), then g′(y)=−2rKy+(r−a−pET). Naturally g′(yc)=−2rK⋅K(1−a+pETr)+(r−a−pET)=−r+a+pET<0, which implies that Ec is LAS. Thus, the equilibrium Ec is locally stable provided it is feasible. Likewise, we can use a similar method to prove Ec is LAS on the Σ2s.
4.
Global dynamics of system (2.3)
This part focuses on examining the global dynamics of a switching system. To certify the global stability of the equilibrium of the system (2.4), it is necessary to rule out the presence of limit cycles. Firstly, we let the Dulac function be V(y,z)=1/yz. For subsystem Si(i=1,2), this gives ∂∂y(V(y,z)dydt)+∂∂z(V(y,z)dzdt)=−rzK<0. Based on the Dulac-Bendixson criterion, it can be concluded that there are no limit cycles present. Consequently, we can derive the following Lemma 4.1.
Lemma 4.1. There exists no limit cycle that is entirely situated within the region Si(i=1,2).
Next, we will exclude limit cycles that intersect with the sliding segment or surrounding the whole sliding segment. Note that this exclusion is necessary for us to follow up with better explanations.
Lemma 4.2. There does not exist a limit cycle that includes a portion of the sliding segment.
Proof. We need to establish the proofs of Lemma 4.2 for the cases ET>z21, z11<ET<z21, z12<ET<z11, z22<ET<z12, and 0<ET<z22.
If ET>z21, the absence of pseudo-equilibrium is deduced from the demonstration of Theorem 3.1, which indicates the presence of dy/dt<0 on the sliding segments Σis(i=1,2) under such circumstances. Therefore, any trajectory reaching Σ1s or Σ2s will approach the boundary points (y21,ET), and E21 is a real stable node. Note that (y21,ET) is visible Σ-fold points of subsystem S2 (see Definition 2.4 in Section 2). Thus, the trajectory initiating at (y21,ET) tends to either approach the stable state E21 directly or in a spiral manner [shown in Figure 1(a)], without touching the switching line again. Therefore, there are no closed orbits that include any part of the sliding segment.
If z11<ET<z21, we know that the pseudo-equilibrium Ec is LAS on the segment Σ1s (see Theorem 3.1 in Section 3) under this scenario, which means the nonexistence of a limit cycle that contains part of the sliding segment Σ1s. Beyond that, there is dy/dt<0 on the segment {(y,z)|y12<y<y22,z=ET} when z11<ET<z21. This implies that the orbits reaching sliding segment Σ2s firstly slide towards the boundary point (y12,ET), which is a visible Σ-fold point as well, enter into the S1, and then tend to Ec. Therefore, there exists no limit cycle incorporating any portion of the sliding segment.
If z12<ET<z11, there exists one real stable equilibrium E11 in region S1. Thus, we can employ a way similar to that applied in the first case to demonstrate the conclusion [shown in Figure 1(b)].
If z22<ET<z12, we have that pseudo-equilibrium Ec is LAS on the sliding segment Σ2s. We get dy/dt>0 on the segment Σ1s in this case. Clearly, the proof process is similar to the second case, so we have omitted here.
If 0<ET<z22, we know that there is no pseudo-equilibrium, and dy/dt>0 holds on the sliding segments Σis(i=1,2), since there exists a stable equilibrium E1, and (y22,ET) is visible Σ-fold point. Thus, the orbits starting from segments Σis(i=1,2) move from the left to the right along the sliding line to the boundary point (y22,ET). Then, they proceed into region S1 and eventually converge to E1, without experiencing hitting the switching line z=ET again. Thus, the proof for Lemma 4.2 is thereby completed.
Significantly, the following lemma is similar to previous studies [37,38], which is obtained by constructing a cycle around the sliding segment and then using the method of counter-evidence to draw a contradiction.
Lemma 4.3. There exists no closed orbit containing the whole sliding segment.
Combined with analysis in Section 2, we obtain the two sets of key parameters, i.e., R0,Rc1,Rc2, and c2,c4,c∗,c∗∗,c†,c††. They are the crucial elements that determine the dynamic behavior of the system (2.4). Next, we will discuss the global dynamics of the proposed piecewise system based on the relationships between the parameters mentioned earlier. Considering that there exists only one uninfected equilibrium E0 for both subsystem S1 and S2 when R0<1, we will not delve further into it. Here, we will focus on the following situations:
Case(a): R0>Rc1>Rc2>1,
Case(b): Rc1>Rc2>R0>1,
Case(c): Rc1>R0>Rc2>1.
4.1. The global dynamics for Case (a)
The global dynamics for Case (a) are the main focus of this section. Note that the analysis methods for the three cases are analogous, thus we have omitted detailed proofs for the other cases. According to the results in Section 2, it can be seen that the relationship between immune intensity c and parameters c2,c4,c∗,c∗∗,c†,c†† have a significant impact on the dynamics of the system (2.4). Through direct analysis, we will consider the under situations:
To understand the dynamics of the system (2.4) more comprehensively, we present the stabilities of various equilibriums completely in Table 3, considering that the dynamics of the system also depend on the relationship between threshold ET and z11,z12,z21, and z22. Thus, we have the following threshold levels:
4.1.1. Case(a1):ET>z21
In such a case, if equilibriums E11,E12,E21, and E22 exist, it can be noted that E11 and E12 are virtual, whereas E21 and E22 are real, with E21 being LAS. We will examine the following six cases in the light of the connections between c and c4,c2,c††,c∗∗,c†,c∗.
Subcase (ⅰ): Assume c4<c<c2, here we know E1 is LAS. Note that equilibria E11 and E12 do not exist under this scenario. All the orbits in S1 will reach the switching line z=ET within a finite amount of time, then firstly enter S2 and ultimately approach either E21 or E1 [shown in Figure 2(a)], contingent upon their initial positions. Thus, the E21 and E1 are bistable in this particular scenario.
Subcase (ⅱ): Assume c2<c<c††, we have immune-free equilibrium E1 is LAS in subsystem S2. According to Lemma 4.2, there exists dy/dt<0 on the Σis(i=1,2) when ET exceeds z21. Therefore, all the orbits of subsystem S1 will slide from right to left to the (y21,ET) or (y12,ET) when they reach the sliding segments, and finally tend to E21. Hence, E21 and E1 are bistable [illustrated in Figure 2(b)].
Subcase (ⅲ): Assume c††<c<c∗∗, here we have equilibrium E1, real equilibrium E21 are LAS. In this subcase, system (2.4) has the same bistable behavior as the previous case [shown in Figure 2(c)].
Subcase (ⅳ): Assume c∗∗<c<c†, we get that E1 is US, and there is only one stable equilibrium E21. Thus, any orbit starting from region S1 firstly crosses the switching line z=ET and enters S2, following the dynamics of S2 tends to E21. Additionally, Lemmas 4.1–4.3 confirms the absence of a limit cycle, implying that all trajectories ultimately converge to E21. Therefore, it can be concluded that the equilibrium E21 is GAS [shown in Figure 2(d)].
Subcase (ⅴ): Assume c†<c<c∗. In such a subcase, we know equilibrium E22 does not exist. In consideration of the nonexistence of the limit cycle, E21 becomes a glocally stable equilibrium [shown in Figure 2(e)].
Subcase (ⅵ): Assume c>c∗. Under this scenario, we know equilibria E12 and E22 do not exist. That is, there exists not any other stable equilibrium other than E21. Considering that the existence of the limit cycle is excluded, endemic equilibrium E21 is GAS [shown in Figure 2(f)]. The dynamics of the system (2.4) can be summarized below when ET>z21.
Theorem 4.1. Suppose ET>z21, it can be deduced that both E21 and E1 are bistable for c4<c<c2,c2<c<c††,c††<c<c∗∗; positive equilibrium E21 is GAS for c∗∗<c<c†, c†<c<c∗ and c>c∗.
4.1.2. Case(a2): z11<ET<z21
In such a case, if equilibriums E11,E12,E21, and E22 exist, we have E11, E12, E21 are virtual, while E22 is real but US. Meanwhile, our analysis reveals the occurrence of sliding mode and the emergence of pseudo-equilibrium within the sliding segment Σ1s={(y,z)|y21<y<y11,z=ET}. By Theorem 3.1, we know that the pseudo-equilibrium Ec is LAS when it exists. Likewise, we will analyze the following six scenarios.
Subcase (ⅰ): Assume c4<c<c2. We know equilibria E11 and E12 do not exist, then we have omitted the description for this case.
Subcase (ⅱ): Assume c2<c<c††. In such a subcase, E1 is LAS in the subsystem S2 and the pseudo-equilibrium will appear on the sliding segment Σ1s. Thus, the partial orbits starting from region S1 will follow the dynamics of S1 to the sliding segment {(y,z)|y21<y<yc,z=ET}, whereas partial trajectories starting from subsystem S2 will arrive at the segment {(y,z)|yc<y<y11,z=ET} along the S2, depending on the initial point, and both types of orbits will eventually converge toward pseudo-equilibrium Ec. Therefore, we conclude that Ec and E1 are bistable [shown in Figure 3(a)].
Subcase (ⅲ): Assume c††<c<c∗∗. In such a subcase, we get that both pseudo-equilibrium Ec on the sliding segment Σ1s and equilibrium E1 are LAS. Thus, we can obtain the coincident conclusion by using similar ways of discussion. That is, the orbits will either go to Ec or tend to E1 along the subsystem S1 and S2, respectively [shown in Figure 3(b)].
Subcase (ⅳ): Assume c∗∗<c<c†. In such a subcase, Ec is the only stable equilibrium for the system (2.4). Our findings reveal that all trajectories intersecting with the line z=ET and following the sliding segment Σ1s reach the pseudo-equilibrium Ec. Considering that a limit cycle does not exist for the entire system, we can conclude that Ec is GAS [illustrated in Figure 3(c)].
Subcase (ⅴ): Assume c†<c<c∗. In such a subcase, we know E22 does not exist and the pseudo-equilibrium Ec is LAS. Here, we have immune-free equilibrium E1 is US. Equally, system (2.4) does not exist any limit cycle. Then the equilibrium Ec is GAS [shown in Figure 3(d)].
Subcase (ⅵ): Assume c>c∗. Under this condition, equilibria E12 and E22 do not exist. There is only one stable equilibrium Ec in the switching system (2.4). Considering the exclusion of closed orbit, then Ec is GAS [illustrated in Figure 3(e)]. Hence, we summarize the aforementioned conclusion of system (2.4) to the following when z11<ET<z21.
Theorem 4.2. Suppose z11<ET<z21, we can conclude that system (2.4) has bistable behavior, i.e., immune-free equilibrium E1 and pseudo-equilibrium Ec are LAS for c2<c<c††,c††<c<c∗∗; equilibrium Ec is GAS for c∗∗<c<c†, c†<c<c∗ and c>c∗.
4.1.3. Case(a3): z12<ET<z11
In such a case, if equilibriums E11,E12,E21,E22 exist, we get E12 and E21 are virtual, but E11 and E22 are real, where E11 is LAS in the S1. We prove the existence of sliding mode but there is no pseudo-equilibrium. Next, we will analyze the stability of equilibria based on the connections between c and c4,c2,c††,c∗∗,c†,c∗.
Subcase (ⅰ): Assume c4<c<c2. Under this scenario, it follows from Section 2 that the two equilibria E11 and E12 are not feasible. So we can omit the description for this case.
Subcase (ⅱ): Assume c2<c<c††. In such a subcase, E1 is LAS and E11 is a real and stable equilibrium, since subsystem S1 has only one stable endemic state. Thus, trajectories initiating from S1 will intend to approach the equilibrium E11. In this scenario, the endemic equilibrium E11 can coexist with the immune-free equilibrium E1. That is, system (2.4) has bistable behavior [shown in Figure 4(a)].
Subcase (ⅲ): Assume c††<c<c∗∗. Regarding the presence and stability of the equilibriums, the dynamics exhibited by subsystems S1 and S2 resemble those of the former scenario [shown in Figure 4(b)]. Thus, we get that E11 and E1 are bistable for the Filippov system (2.4).
Subcase (ⅳ): Assume c∗∗<c<c†. In this subcase, the equilibrium E21 is virtual and the system (2.4) only exhibits one stable endemic equilibrium E11. As there is no limit cycle, then E11 is GAS [shown in Figure 4(c)].
Subcase (ⅴ): Assume c†<c<c∗. Analogously, there exists only one stable equilibrium point E11 in subsystem S1. It should be noted that the endemic equilibrium E11 functions as an attractor. We have excluded the existence of limit cycles. Consequently, any orbit starting from region S1 or S2 will approach the equilibrium point E11 [shown in Figure 4(d)]. Thus, the equilibrium E11 is GAS.
Subcase (ⅵ): Assume c>c∗. In this scenario, since the equilibria E12 and E22 do not exist, we have omitted the description for this case and we get the following conclusion.
Theorem 4.3. Suppose z12<ET<z11, we can conclude that immune-free equilibrium E1 and endemic equilibrium E11 are bistable for c2<c<c†† and c††<c<c∗∗; equilibrium E11 is GAS for c∗∗<c<c† and c†<c<c∗.
4.1.4. Case(a4): z22<ET<z12
In such a case, if equilibriums E11,E12,E21, and E22 exist, we have E21 is virtual, while E11, E12, and E22 are real. Note that both equilibriums E12 and E22 are US. Furthermore, the sliding mode does exist, and the pseudo-equilibrium Ec is LAS on the Σ2s. A similar discussion works for z22<ET<z12.
Subcase (ⅰ): Assume c4<c<c2, we know equilibria E11 and E12 do not exist according to Section 2. Therefore, we can ignore the explanation of this situation.
Subcase (ⅱ): Assume c2<c<c††. In such a subcase, E1 is LAS and E11 is real and stable in S1. As mentioned above, the pseudo-equilibrium Ec is LAS. From the dynamics of subsystems S2 and S1, it can be deduced that the orbits will either directly reach the E11, E1 or firstly arrive at the line z=ET on the segment {(y,z)|y12<y<y22,z=ET}, then slide to the Ec along sliding segment Σ2s, depending on the initiating points. Hence, as is shown in Figure 5(a), the equilibriums E11, Ec, and E1 are tristable.
Subcase (ⅲ): Assume c††<c<c∗∗. In such a subcase, we get that the equilibrium E1 is LAS, and there exist two equilibria E11 and Ec, which are locally stable in their respective regions. Thus, we can conclude that the dynamic behaviors in this scenario are consistent with the former subcase. That is, system (2.4) has tristable behavior in this case [shown in Figure 5(b)].
Subcase (ⅳ): Assume c∗∗<c<c†. Note that the locally stable equilibrium E21 is virtual, then cannot be attained. Part of the trajectories starting from subsystem S1 will approach the equilibrium E11, while certain trajectories will collide on the switching line at finite time, and the locally stable pseudo-equilibrium Ec appears ultimately on the segment Σ2s={(y,z)|y12<y<y22,z=ET}. Hence, they are bistable for the switching system (2.4) [shown in Figure 5(c)].
Subcase (ⅴ): Assume c†<c<c∗. We know equilibrium E22 does not exist, then we rule this out.
Subcase (ⅵ): Assume c>c∗. In such a subcase, we know equilibria E12 and E22 do not exist, so we have omitted the description in this case as well. To sum up, we can derive the conclusion as follows.
Theorem 4.4. Suppose z22<ET<z12, we can conclude that the real equilibrium E11, pseudo-equilibrium Ec, and the immune-free equilibrium E1 are tristable for c2<c<c†† and c††<c<c∗∗; immune-free equilibrium E11 and pseudo-equilibrium Ec are bistable for c∗∗<c<c†.
4.1.5. Case(a5): 0<ET<z22
In such a case, if equilibriums E11,E12,E21, and E22 exist, we have E21 and E22 are virtual, but E11, E12 are real, where E11 is LAS in the S1. Under this scenario, there is no pseudo-equilibrium for the Filippov system (2.3). Subsequently, we will analyze the following situations based on the relationships between immune intensity c and c4,c2,c††,c∗∗,c†,c∗.
Subcase (ⅰ): Assume c4<c<c2. In this condition, we get that equilibria E11 and E12 do not exist on the basis of the discussion of Section 2. Note that the locally stable equilibrium E21 is virtual. Thus, there is not any other stable equilibrium besides E1, as all the possible limit cycles have been excluded. Hence, we have E1 is GAS [shown in Figure 6(a)].
Subcase (ⅱ): Assume c2<c<c††. In such a subcase, E1 is LAS and E11 is a real and stable equilibrium. However, E21 is virtual in S1. Thus, positive equilibrium E11 and immune-free equilibrium E1 are bistable for c2<c<c†† being satisfied [shown in Figure 6(b)].
Subcase (ⅲ): Assume c††<c<c∗∗. We know that equilibria E11 and E1 are LAS as well. Thus, the dynamics of this case are similar to the former, i.e., the bistable behavior occurs [shown in Figure 6(c)].
Subcase (ⅳ): Assume c∗∗<c<c†. In such a subcase, note that locally stable equilibrium E21 is virtual and then cannot be attained. There exists only one stable equilibrium point E11 in subsystem S1. As there is no limit cycle, we obtain that the equilibrium E11 is GAS [illustrated in Figure 6(d)].
Subcase (ⅴ): Assume c†<c<c∗. In such a subcase, we know that equilibrium E22 does not exist, so we have ignored the description of this situation.
Subcase (ⅵ): Assume c>c∗. In such a subcase, both E12 and E22 do not exist, therefore, the explanation for this case is omitted. Consequently, we arrive at the subsequent conclusion.
Theorem 4.5. Suppose 0<ET<z22, we can conclude that immune-free equilibrium E1 is GAS for c4<c<c2; the real equilibrium E11 and immune-free equilibrium E1 are bistable for c2<c<c†† and c††<c<c∗∗; the equilibrium E11 is GAS for c∗∗<c<c†.
In fact, extensive discussions have been conducted on the stability of various equilibriums for the switching system (2.4) under Case(a). Furthermore, Table 4 provides a comprehensive summary of the global dynamics associated with these specific cases.
5.
Boundary equilibrium bifurcation
Notably, the collision of pseudo-equilibrium, tangent point, and regular equilibrium (or tangent point and regular equilibrium) in switching systems occurs when ET reaches a critical value, leading to boundary equilibrium bifurcations on the discontinuity surface [as shown in Figure 7]. The understanding and analysis of these boundary equilibrium bifurcations are crucial in studying the dynamical behavior of the Filippov system. To verify the boundary equilibrium bifurcation, we select ET as the bifurcation parameter while keeping all other parameters constant. Detailed explanations of the boundary equilibrium and the tangent point are shown in Definitions 2.3 and 2.5.
Boundary equilibrium of the Filippov system (2.4) satisfies
where ϕ=0 or ϕ=1. By solving the equations provided in (5.1), it is possible to obtain four potential boundary equilibria
Tangent points of the Filippov system (2.4) satisfy
Thus, the potential tangent points can be denoted as
which are the solutions of (5.4) corresponding to ϕ=0 and ϕ=1.
Figure 7 examines a series of boundary equilibrium bifurcations when c2<c<c††. In this case, both subsystems have two positive equilibria. The real and stable equilibrium E21 coexists simultaneously with the visible tangent point T21 when ET>z21 [shown in Figure 7(a)]. As ET decreases from ET>z21 to z21, E21 collides with T21 [shown in Figure 7(b)]. With the threshold ET decreasing further to z11<ET<z21, a stable pseudo-equilibrium Ec emerges and T21 transforms into an invisible tangent point [as depicted in Figure 7(c)]. This bifurcation exhibits the progress of the formation of Ec. Furthermore, boundary bifurcation takes place again when ET through the critical value z11. This case results in the collision of the tangent point T11, the equilibrium point E11 with the pseudo-equilibrium point Ec [as depicted in Figure 7(d)]. Subsequently, the Ec vanishes, and the stable point E11 transforms into the local attractor [as illustrated in Figure 7(e)]. When ET drops consistently until z12, the third boundary bifurcation takes place, leading to the collision of the visible tangent point T12 with the equilibriums E12 [as depicted in Figure 7(f)].
Provided ET continues to decrease until z22<ET<z12, a locally stable pseudo-equilibrium Ec appears [as shown in Figure 7(g)] and a tristable phenomenon (Ec, E11, and the immune-free equilibrium E1) occurs. When ET passes E22, the fourth boundary equilibrium bifurcation occurs. In this scenario, the pseudo-equilibrium Ec will collide with the equilibrium point E22 and tangent point T22 if ET=z22 [illustrated in Figure 7(h)]. However, as the threshold ET continues to decrease, the equilibrium Ec disappears, and the tangent point T22 becomes invisible [illustrated in Figure 7(i)]. Here, the equilibriums E1 and E11 exhibit bistability under this scenario.
6.
The influence of the pivotal parameters of the system (2.4)
In order to stabilize the HIV viral loads and effector cell counts within the required predetermined level, it is crucial to implement a control strategy for the switching system (2.4) by setting an appropriate threshold ET. The dynamics of the system are influenced by two key parameters: The immune intensity c and the inhibition of the virus on the proliferation of immune cells η. Therefore, it is important to study the impact of these parameters on the system's dynamics. Note that the parameters in this study are based on the findings of Komarova[33] and colleagues. A direct calculation reveals that the bistable interval is (2.7666, 3.2333) [shown in Figure 8(a)]. From this figure, it is evident that system (2.4) has the potential to have either one or two nontrivial LAS equilibria depending on the value of c. In detail, two stable equilibriums E21 and E21 coexist when c4<c<c††. As established in Section 2, there is a unique LAS equilibrium E21 if c>c††. It is important to note that a saddle-node bifurcation occurs at c=c4. Generally, a longer bistable interval implies a wider range of variation for the proliferation coefficient c of immune cells. In this context, a lower viral inhibition intensity η [shown in Figure 8(b)] is more advantageous for immune control. This indicates the necessity of developing medications aimed at diminishing the viral inhibitory effect on immune cells. The phase portrait of this system reveals that both variables y and z will gradually tend to a stable value over time t changed [shown in Figures 8(c) and (d)]. In subsystem S2, the combination of antiretroviral therapy and immunotherapy not only controlled the number of viruses better but also stabilized the immune cell count at a more reliable level than the free system, as can be seen in Figures 8(e) and (f).
7.
Conclusions
Given the immunosuppressive HIV infection model[32,33,34], it is suggested that the threshold policy control (TPC) method for treating infected patients should be based on the number of immune cells. Thus, we propose a piecewise immunosuppressive infection system. We call model (2.4) a discontinuous right-hand side dynamical system under TPC, which consists of two subsystems. Specifically, it is recommended to initiate a combination of antiretroviral therapy and immune therapy when the value drops below a certain threshold ET, which includes the usage of antiretroviral medications and interleukin (IL)-2. This leads to the emergence of a nonsmooth system. Unlike [31], we consider the logical growth of HIV rather than linear growth in this paper. Some clinical facts indicate that the saturated growth of HIV is reasonable.
Initially, we provide a concise overview of the dynamics exhibited by the two subsystems. Through the subsystems Si(i=1,2), we derive the thresholds R0, Rc1, and Rc2. It becomes evident that R0 plays a pivotal part in determining the eradication of the virus. We also obtain the post-treatment control thresholds c2(c4) and the elite control threshold c∗∗(c††) for subsystems S1 (S2). According to [34], there exists a bistable behavior between these two threshold intervals. The sliding dynamics and sliding domain of the system (2.4) are studied in the subsequent analysis. Our purpose is to demonstrate the existence of two sliding segments Σis(i=1,2). By employing the Filippov convex approach, we investigate the possibility and local asymptotic stability of the pseudo-equilibrium Ec on the sliding segment Σ1s under the z11<ET<z21, or on the sliding segment Σ2s under the z22<ET<z12. Significantly, we have primarily focused on Case (a) and discussed the global dynamics of the system in this paper. To investigate the global dynamic behavior of the system, we have excluded the existence of three types of limit cycles. It is important to understand the relationship not only among R0,Rc1,Rc2, and 1 but also among immune intensity c and c4,c2,c††,c∗∗,c†,c∗. Subsequently, the bifurcation theories were utilized to address the dynamics of sliding mode and local sliding bifurcations.
The analysis reveals that the system can demonstrate diverse and complex dynamic behaviors: (ⅰ) One of the equilibria in the system is GAS, which can manifest as the immune-free equilibrium E1, pseudo-equilibrium Ec, or even as the positive equilibrium E11 or E21 within subsystems S1 or S2; (ⅱ) There are two possible equilibria in this system that exhibit bistability, namely the immune-free equilibrium E1 and equilibrium E21 (or Ec or E11), or the positive equilibrium E11, which is bistable alongside the pseudo-equilibrium Ec; (ⅳ) Three equilibria are tristable, i.e., immune-free equilibrium E1, positive equilibrium E11, and the pseudo-equilibrium Ec are stable for z22<ET<z12 and c2<c<c∗∗. Our work demonstrates that the utilization of effector cell-guided therapy leads to an expansion of the controllable area of initial values for patients, generating a more complex Filippov dynamics system when compared with [35]. Interestingly, we find that there exists an optimal threshold interval for immune intensity that can maximize the controllable area of initial values. This highlights the importance of considering the effects of effector cell-guided therapy and immune intensity when studying the dynamics of the switching system. It suggests that maximizing the controllable area of initial values can potentially improve the effectiveness of treatment strategies for patients.
In this paper, the existence of three types of equilibria including pseudo-equilibrium is explored. These equilibria can exhibit bistability or tristability, meaning that the HIV viral loads and effector cell counts can be stabilized at a preset level. Achieving these stable states depends on factors such as the threshold level, immune intensity, and the initial values of the system. Consequently, determining the optimal strategy for immune intensity and the threshold conditions should still take into account the individual characteristics of the patients. In the case of boundary H(X), the selection of parameter values is crucial in stabilizing different equilibria within the system. From a biological standpoint, employing rational control intensity and intervention is highly effective in ensuring the control and management of diseases.
Although our research has an impact on HIV disease control, it is still insufficient. We only consider the relationship between the number of effector cells and the threshold level to construct switching conditions. The actual disease control strategy should also take into account the change rate of effector cell count, which will be our next work. By considering these factors, we aim to further provide insights into the effectiveness and impact of this treatment approach on the virus-immune system dynamics. This mathematical model could potentially contribute to the improvement of treatment strategies for viral infections.
Use of AI tools declaration
The authors declare that no Artificial Intelligence (AI) tools were used in the creation of this article.
Acknowledgments
This work was supported by The National Natural Science Foundation of China (Grant No. 12261033).
The authors would like to thank the editor and anonymous referees for their valuable comments and suggestions, which have led to an improvement of the paper.
Conflict of interest
The authors declare no conflict of interest.