Bistability and Robustness for Virus Infection Models with Nonmonotonic Immune Responses in Viral Infection Systems

: Recently, bistable viral infection systems have attracted increased attention. In this paper, we study bistability and robustness for virus infection models with nonmonotonic immune responses in viral infection systems. The results show that the existing transcritical bifurcation undergoes backward or forward bifurcation in viral infection models with nonmonotonic immune responses. Our investigation demonstrates that the backward bifurcation threshold is the elite control threshold. When the immune intensity is greater than the elite control threshold, the virus will be under elite control; when the immune intensity is less than the elite control threshold, the virus may rebound. We also give a new deﬁnition of robustness to characterize bistable systems.


Introduction
Viral infection has always been a major threat to human health, and the human immune system plays a key role in fighting viruses. When a virus first begins to infect the human body, host cells rapidly activate natural killer cells, macrophages and so on, so as to participate in nonspecific immunity. Then, after the virus infects the human body for a period of time, the immune system will activate cytotoxic T lymphocyte (CTL) cells and antibody cells to produce a specific immune response.
In many viral infections, toxicity T lymphocyte (CTL) cells and antibody cells will attack infected cells to destroy the virus. These two kinds of cells play extremely important roles in antiviral responses. There are many mathematical models formulated to study the within-host virus dynamics with and without an immune response in recent studies [1][2][3][4][5][6]. Furthermore, there are many studies considering age-structure [7,8] or reaction-diffusion [9,10] in virus dynamics.
Many studies focus on nonmonotonic functional responses, such as in population models [11,12], virus models with immune response [2][3][4] and epidemic models [13,14]. The general form of nonmonotonic functions named Monod-Haldane functions was given by Andrews [15] in 1968. Recently, Wang et al. [1] used the Monod-Haldane function to model the nonmonotonic immune response in virus dynamics. They investigated and studied three-dimensional and two-dimensional viral infection systems with nonmonotonic immune responses as follows

Bifurcations Analysis
In this section, we prove the existence of backward and forward bifurcation in system (1) and system (2). In fact, backward/forward bifurcation was well studied in [23] for the first time and was also studied in [24]. In the following, we use the method in [23,24] to prove the existence of backward/forward bifurcation.
We first denote B = γb − c and c 2 = γb + 2b √ α. Then, for system (1) and (2), we give some thresholds, and the basic reproduction number R 0 gives the average number of infections transmitted by a single infected individual among fully susceptible individuals.
To find the R (1) 0 of system (1), we follow the next-generation matrix method proposed by van den Driessche and Watmough [23]. Thus, for system (1), let X = (x, y, z) T and rewrite system (1) as dX dt = F − V, where F is the rate at which new infections occur, and V is all other traffic inside and outside of each compartments. Thus, we have The form of the next generation matrix is Now, according to Theorem 2 in [23], the spectral radius ρ of the matrix FV −1 is the maximum eigenvalue of FV −1 , which gives the basic reproduction number R (1) 0 of the system (1). Thus, we obtain R (1) . In addition, we also denote the other thresholds as follows .
As for system (2), Using the same method in [23], we can find and the next generation matrix of system (2) is Thus, the basic reproduction number of system (2) is R Then, we denote the other thresholds and the equilibria E and These thresholds were introduced in [1], and the stability of the equilibria in system (1) and (2) has also been analyzed, which can be concluded from Tables 1 and 2. Where GAS represents global asymptotic stability, LAS represents local asymptotic stability, US represents instability, and-represents that the equilibrium does not exist. Table 1. The stabilities of the equilibria and the behaviors of system (1) in the case that R Table 2. The stabilities of the equilibria and the behaviors of system (2) in the case that R

Bifurcations in Three-Dimensional Model
In this section, we prove the existence of backward and forward bifurcations in system (1) using the theorem developed by Castillo-Chavez and Song [24]. From the results in [1], letẼ be any arbitrary equilibrium of system (1). The Jacobian matrix associated with system (1) is The characteristic equation of the linearized system of (1) atẼ is given by λI − J 1 = 0.
c > 1 and c = c * * 1 , then system (1) undergoes a backward bifurcation at E corresponding to the zero eigenvalue, respectively. Then, we can find a right eigenvector associated with zero eigenvalue, given by The left eigenvector is given by w = (0, 0, 1) from Theorem 3.10 in [1]. Let (x, y, z) = (x 1 , x 2 , x 3 ) and G = (g 1 , g 2 , g 3 ) T . From the proof of Theorem 3.10 [1], we have We know that, when R c and c = c * * 1 , we have y [1]. Then, b > 0 and α − (y (1) 1 ) 2 < 0, implying that a > 0. Therefore, according to Corollary 4.1 in [24], both a and b are positive, implying that system (1) undergoes a backward bifurcation at c = c * * 1 . Here, c is the bifurcation parameter. Table 3 in [24], we can see that, when the vertical axis represents z (Immune cells) and x (Activated CD4 + T cells), the bifurcation diagrams are drawn in Figure 1a,c. Considering Remark 1 in [24] then we can find a < 0 and b > 0 by calculation. Thus, when the vertical axis represents y (Infected CD4 + T cells), the bifurcation diagram is drawn in Figure 1b.  Table 3. The solid line represents stable equilibria, and the dotted line represents unstable equilibria. The post-treatment control threshold is c 2 ≈ 0.3, the elite control threshold is c * * ≈ 0.3837, and the bistable interval is (0.3, 0.3837).
c , there also exists forward bifurcation in system (1), which will be proven in the following.
c and c = c * * 1 , then system (1) undergoes a forward bifurcation at E (1) 1 . Here, c is the bifurcation parameter.
Proof. We know that, when 1 < R Therefore, according to Theorem 4.1 in [24], we know that, when a < 0 and b > 0, system (1) undergoes a forward bifurcation at c = c * * 1 . Here, c is the bifurcation parameter. Table 3 in [24], we can see that the vertical axis represents z (Immune cells) and x (Activated CD4 + T cells) in the bifurcation diagrams plotted in Figure 2a then we can find a > 0 and b > 0. Thus, when the vertical axis represents y (Infected CD4 + T cells), the bifurcation diagram is plotted in Figure 2b.    Table 3 and α = 10. The solid line represents stable equilibria, and the dotted line represents unstable equilibria. The elite control threshold is c * * 1 ≈ 0.7549 .

Bifurcations in Two-Dimensional Model
In this section, we use the same method to prove the existence of backward bifurcation in system (2). From the results in [1], letẼ be any arbitrary equilibrium of system (2). The Jacobian matrix associated with the system (2) is The characteristic equation of the linearized system of (2) atẼ is given by λI − J 2 = 0.
c > 1 and c = c * * 2 , then system (2) undergoes a backward bifurcation at E corresponding to the zero eigenvalue, respectively. Then, the right eigenvector associated with 0 eigenvalue is ξ = (− pK r , 1) T , and the left eigenvector θ is θ = (0, 1). Let F = (P 2 , Q 2 ) T . From the proof of Theorem 5.11 in [1], we can find We know that, when R . Thus, it is clear that b > 0 and α − (y (2) 1 ) 2 < 0 implying that a > 0. Therefore, according to Corollary 4.1 in [24], both a and b are positive, implying that system (2) undergoes a backward bifurcation at c = c * * 1 . Here, c is the bifurcation parameter.
c . From Table 3 in [24], we can see that, when the vertical axis represents z (Immune cells), the bifurcation diagram is plotted in Figure 3b. Considering Remark 1 in [24], if we choose ξ = (1, − r pK ) T and θ = (0, − pK r ), then we can find a < 0 and b > 0 by calculation. Thus, the vertical axis represents y (Infected CD4 + T cells) in the bifurcation diagram plotted in Figure 3a.   Table 4 and K = 6 cells/µL. The solid line represents stable equilibria, and the dotted line represents unstable equilibria. The post-treatment control threshold is c 2 ≈ 2.5, the elite control threshold c * * 2 ≈ 3.8333, and the bistable interval is (2.5, 3.8333).
c and c = c * * 2 , then system (2) undergoes a forward bifurcation at E (2) 1 . Here, c is the bifurcation parameter.
Proof. We use the same method as Theorem 2 and the results in Theorem 3. When 1 < R (2) 0 < R (2) c and c = c * * 2 , we have y Therefore, according to Theorem 4.1 in [24], we know that, when a < 0 and b > 0, system (2) undergoes a forward bifurcation at c = c * * 2 . Here, c is the bifurcation parameter.
c . From Table 3 in [24], we can see that the vertical axis represents z (Immune cells), and the bifurcation diagram is plotted in Figure 4b. Considering Remark 1 in [24], if we choose ξ = (1, − r pK ) T and θ = (0, − pK r ), then we can find a > 0 and b > 0 by calculation. Thus, the vertical axis represents y (Infected CD4 + T cells), and the bifurcation diagram is plotted in Figure 4a.    Table 4 and K = 1.4 cells/µL. The solid line represents stable equilibria, and the dotted line represents unstable equilibria. The elite control threshold c * * 2 ≈ 2.6286.

Robustness for Bistable System
When a system with backward bifurcation is bistable, it has two basins of attraction. Robustness is one of the fundamental characteristics of biological systems. A general definition of the robustness for biological systems can be used: robustness, the ability to maintain performance in the face of perturbations and uncertainty, is a long-recognized key property of living systems [25].
According to the general definition, we propose the definition of robustness for bistable systems, which is measured by the area of the basin of attraction, namely the definite integral of positive steady solution curve on the bistable interval. The robustness corresponding to x, y and z for bistable system (1) can be calculated by As shown in Figure 1, the area covered by cyan denotes R z (c), the robustness for bistable system (1) corresponding to x, y and z, respectively. Here, we use the default parameter in Table 3, and we can calculate R (1) Using the same method, we can define the robustness corresponding to y and z for bistable system (2) as follows Then, the area covered by cyan in Figure 3 denotes R (2) y (a) and R (2) z (b), the robustness for bistable system (2) corresponding to y and z, respectively. Here, we use the default parameter in Table 4, and we can calculate R (2) y ≈ 2.9014 and R (2) z ≈ 1.0986. When the system has some parameter and initial value perturbations, if the parameter and initial values are always in one basin of attraction, then the system will maintain its performance. However, if the parameter and initial values change the basin of attraction to another due to the perturbations, then the system will appear to have a regime shift and change the original performance. Thus, the larger the attraction is, the more difficult to find a regime shift. Considering systems (1) and (2), we can think of robustness R (1) x , R   y (R (2) y ) is, the better robustness for bistable system, and then the tumor cells have greater difficulty in rebounding. However, the greater R (1) x and R (1) z (R (2) z ) are, the easier the system to be in the steady state of high virus load-that is, the virus can rebound more easily.

Sensitive Analysis
Sensitive analysis is an important method to test the influence of parameter variations on system dynamics [26]. In the following, we conduct sensitive analysis with the aim of revealing the relationship between the elite control threshold and the basic immune reproduction number and system parameters in our model (1) and model (2). Here, we use Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCs) [27,28] to test the dependence of c * * 1 , c * * 2 and R 1− * , R 2− * . As a statistical sampling method, LHS provides an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [27].
PRCC, on the other hand, shows the level of significance for each parameter. The PRCC is obtained using the rank transformed LHS matrix and output matrix [28]. We performed 4000 simulations per run and used a uniform distribution function to test for the significance of PRCCs for all parameters with wide ranges. When |PRCC| > 0.4, there was a significant correlation between the input parameters and output variables. For |PRCC| ∈ (0.2, 0.4], the correlations were moderate. When |PRCC| ∈ [0, 0.2], we had weak correlations. The PRCC results in Figure 5 illustrate the dependence of c * * 1 and R 1− * on each parameter in system (1). Furthermore, Figure 6 illustrates the dependence of c * * 2 and R 2− * on each system parameter in system (2).  Table 3.  Figure 6. Partial rank correlation coefficients illustrating the dependence of c * * 2 (a) and R 2− * (b) for the system (2) on each parameter. The parameter values are listed in Table 4.
From Figure 5 for system model (1), we notice that the generated rate of CD4 + T cells s, the death rate of CD4 + T cells d , the infected rate β, the treatment rate , the death rate of infected cells a and the death rate of immune cells b have significant influences on the exile control threshold c * * 1 ; and the generated rate of CD4 + T cells s, the death rate of CD4 + T cells d , the infected rate β, the treatment rate and the death rate of infected cells a have significant influences on the basic immune reproduction number R 1− * . From Figure 6 for system model (2), the growth rate of infected cells r, the environmental capacity K, the death rate of immune cells b and the death rate of infected cells a have significant influences on the exile control threshold c * * 2 ; and the growth rate of infected cells r, the death rate of immune cells b, the death rate of infected cells a and the immune intensity c have significant influences on the basic immune reproduction number R 2− * .

Numerical Simulation
In this section, we perform numerical simulations to verify our analytical results. We perform numerical simulations for system (1) using the default parameters in Table 3.
1 is stable (see Figure 9). When we choose c = 0.7 < c * * 1 and different initial values, only the positive equilibrium E 1− * is stable (see Figure 10).   Table 3. Only E 1− * is stable. Other parameter values are listed in Table 3. Only E 1− * is stable. Other parameter values are listed in Table 3 and α = 10.
We perform numerical simulations for system (2). We use the default parameters in Table 4.
0 , c 2 = 2.5, and c * * 2 ≈ 3.8333. In this case, when we choose c 2 < c = 3 < c * * 2 and different initial values, both the immune-free equilibrium E (2) 1 and the positive equilibrium E 2− * are stable (see Figure 11a). When we choose c = 4 > c * * 1 and different initial values, only the positive equilibrium E 2− * is stable (see Figure 11b). If we choose K = 1.4 cells/µL and different initial values, then R 6286. In this case, when we choose c = 2.5 < c * * 2 and different initial values, only the immune-free equilibrium E (2) 1 is stable (see Figure 12a). When we choose c = 2.7 > c * * 1 and different initial values, only the positive equilibrium E 2− * is stable (see Figure 12b). Figures 11 and 12 verify the results in Figures 3 and 4.

Discussion
In this paper, we found that the nonmonotonic immune responses in the viral infection systems in [1] undergo backward and forward bifurcation at the elite control threshold if we select the immune intensity as the bifurcation parameter under different conditions. When the immune intensity is less than the elite control threshold, the viral load will be under control or rebound. Only when the immune intensity is greater than the elite control threshold will the viral load be under elite control. Such concepts can be applied to the dynamics of infection and immunity of EBV, HIV and SARS-CoV2.
The results show that a combination of antiviral therapy and immunotherapy may contribute to the functional cure of HIV, that is, the immune intensity is greater than the elite control threshold through immunotherapy. Then, with long term control of viral replication and low viral load, a functional cure of HIV can be realized. The HIV patients under elite control still need to be cautious and continue to receive immunotherapy. If the immune intensity is less than the post-treatment control threshold, the virus can rebound again.
Author Contributions: Conceptualization, T.W. and S.W.; methodology, S.W. and F.X.; formal analysis, T.W. and S.W.; writing-original draft preparation, T.W. and S.W.; writing-review and editing, S.W. an F.X.; supervision, S.W. All authors have read and agreed to the published version of the manuscript.