Codimension-one bifurcation and stability analysis in an immunosuppressive infection model

One of the important medical problems is infectious diseases such as HIV and hepatitis which annually causes the death of many people. So it is important to study infectious diseases parametric models. In this paper, we investigate differential equations system of HIV and hepatitis (with delay and without delay) from the stability and codimension-one bifurcation point of view. We show that their dynamical behaviour will change when the parameters vary. We prove that this model has a saddle-node bifurcation and transcritical bifurcation when the delay parameter is absent. Also by using the center manifold theory, we show that the delay model has a saddle-node bifurcation.

The bistability in this model leads to sustained immunity when the treatment is stopped, because a solution from the basis of the attraction of the virus dominant equilibrium can be lifted to that of the immune control equilibrium via a single phase of therapy.
After that, Shu et al. (2014) incorporated the time lag during the immune response process into Komarova et al. 's model and studied the dynamics between an immunosuppressive infection and antiviral immune response.
To formulate their model, they followed the line in Komarova et al. (2003) and Fenton et al. (2006). They considered the following model where y and z denote the virus population size and population size of immune cells, respectively. The virus population is assumed to grow logistically: r is the viral replication rate and a is clearance rate. In addition, they assumed virus is killed by immune cells at a rate pyz and immune cells are assumed to be inhibited by the virus at a rate qyz and died at a rate b. The activation rate of immune cells at time t is assumed to depend on the virus load and the number of immune cells at time t − τ. Here, τ is the time lag accounting for the time needed for the immune system to trigger a sequence of events such as antigenic activation, selection and proliferation of immune cells to produce new immune cells. In model 1, it is important to note that f(y), function of immune expansion by virus load, is considered as follows (Shu et al. 2014) Note that if the time lag is ignored, τ = 0, model 2 reduces to the following model: They studied the local and global stability of the most of equilibria. By using bifurcation theory, they only found Hopf bifurcation in the model when τ = τ bif .
In this paper, we follow the line in Shu et al. (2014). It should be noted that, we detect another equilibrium point which is not considered in Shu et al. (2014). Furthermore, we choose another parameters, r and c, as bifurcation parameters. The parameter r is the viral replication rate and the parameter c is a coefficient in the function of immune expansion by virus load. We consider r and c as bifurcation parameters and obtain the following result: (1) (i) if r = r bif , then the transcritical bifurcation occurs in system 4, (2) (ii) if c = c bif , then the saddle-node bifurcation occurs in system 4, (3) (iii) if c = c bif , then the saddle-node bifurcation occurs in system 2. (2) As we mentioned, Shu et al. (2014) only investigated Hopf bifurcation by considering τ as bifurcation parameter. But we find new equilibrium in their model and obtain new dynamical behaviours in the model. Furthermore, we find other important parameters in studying dynamics of this model. To the best of our knowledge, this is the first time that these results are obtained in this immunosuppressive infection model.
The rest of paper is organized as follows. In the next section, we obtain the necessary condition of existence of equilibria in immunosuppressive infection model. In "Dynamics of the model without delay (system 4)" section, we will consider the dynamics of model 4. The dynamical behaviour of model 2 is investigated in "Dynamics of the model with delay (system 2)" section. In "Numerical simulation" section, the validity of the main results is illustrated by numerical simulations. Finally, we state some main conclusions.

Existence of equilibrium points
For any τ > 0, let C := {φ : [−τ , 0] → R is continuous} be Banach space of continuous function on [−τ , 0] with the norm is defined as �φ� = sup −τ ≤θ ≤0 φ(θ). We denote the nonnegative cone of C by C + . Shu et al. (2014) showed that system 2 with any initial condition (φ, ψ) ∈ C + × C + admits an unique solution and the solution (y(t), z(t)) remains nonnegative for t ≥ 0 and is bounded in C + × C + . Furthermore, they showed that the bounded region where µ = min{a, b} > 0, is positively invariant with respect to system 2 and the system is well posed (Shu et al. 2014). Now we find the equilibria of system 2. We then investigate their stability. As we said in "Background" section, we obtain an equilibrium point that it is not considered in Shu et al. (2014).
Cleary E 0 = (0, 0) is a trivial equilibrium of system 2, this equilibrium means that any virus cell and immune response do not exist in the body. There exists an equilibrium E 1 = (ȳ, 0) = ( k(r−a) r , 0) provided r > a. At equilibrium E 1 does not exist any immune response, also viruses are with positive size. Therefore, we call the equilibrium E 1 the virus dominante equilibrium (VDE). Assume that E * = (y * , z * ) is another equilibrium point of system 2 with y * > 0 and z * > 0 which means immune response and virus cells are present at the same time. Therefore, the virus cells can be controlled. Now, we consider the following equations The first equation of 6 follows that

Dynamics of the model without delay (system 4)
In this section, we provide a complete description about dynamics of system 4. To this end, we begin with the following result on local stability of system 4.
Proof Suppose that (ỹ,z) is an equilibrium of system 4. The associated characteristic equation is given by Then characteristic equation 15 at E 0 has two roots, ξ 1 = −b < 0 and ξ 2 = −(a − r) . If 12 holds, then E 0 is stable and if 13 or 14 holds then E 0 is saddle point.
Also, characteristic equation 15 at E 1 = (ȳ, 0) has two roots, Note to the graph of g(y), it is obvious that g(ȳ) > 0. Now, if 13 or 14 holds, then the equilibrium E 1 is asymptotically stable. We suppose that 15 holds, by substituting E * at Eq. 15, we have Suppose that Condition H 1 follows that ŷ = y * , then g 1 (y * ) = 0 . Therefore, the roots of Eq. 15 are ξ 1 = 0, ξ 2 = − r k y * , or other words E * is locally stable. When r ≥ a, the infection can not spread in body of patient, so there is no virus cell and immune response. In this case, system 4 converges to E 0 . We know viral cells infect the host without immune response as r increases from a to r t . In this case, system 4 converges to E 1 and the equilibrium point E 1 is locally asymptotically stable. By increasing r from r t , immune response increases and controls viral cells. In this case, E * and E 1 exist. Therefore, to obtain the better conditions and control of virus cells, we should converge the system to the equilibrium point E * .

Lemma 3
Assume that H 1 is satisfied, therefore system 4 has a saddle node bifurcation at equilibrium E * when the parameter c varies.
Proof By Lemma 2, characteristic equation 15 at E * has two simple roots ξ 1 = 0 and where the two conditions (a) and (b) are opposed zero. By Sotomayor Theorem (Guckenhiemer and Holmes 1993;Perko 1991), system 4 has a saddle-node bifurcation at E * whene c = c bif .

Lemma 4
If H 1 is satisfied, then system 4 has a trancscritical bifurcation at equilibrium E 0 when r = r bif = a.
Proof By Lemma 2, characteristic equation 15 has two roots ξ 1 = −b and ξ 2 = −(a − r) . Therefore (E 0 , r bif ) is a bifurcation point where r bif = a and ξ 2 = 0 is a simple zero of 15. Now, we assume A = Df (E 0 , r bif ) then the eigenvectors of A and A T at zero eigenvalue are Therefore, we have the following quantities By Sotomayor Theorem (Guckenhiemer and Holmes 1993; Perko 1991), system 4 has a trancscritical bifurcation at E 0 when r bif = a.
According to Lemma 4, we know that system 4 has a transcritical bifurcation at E 0 , when r = r bif . For r ≤ r bif , only equilibrium point E 0 is stable. In this case, the patient , s body does not have virus cells and immune response. Also, with increasing r (r > r bif ), the equilibrium E 1 occurs; in this case the system has a branch of stable equilibrium E 1 and a branch of the unstable equilibrum E 0 that express the transcritical bifurcation. In the branch of the stable equilibirum E 1 , the patient has a viral cells without any immune response. Therefore as shown if the viral replication rate r is greater than the threshold r t , then the two equilibrium points E 1 and E * at the same time are stable and the bistability phenomenon occurs. Also, we know that for c < c bif , there is no equilibrium E * and according to assumption H 1 at c = c bif , the equilibrium E * will be found. After passing through c bif (c > c bif ); according to Shu et al. (2014), the system has two equilibrium E * 1 and E * 2 . This means that there is a saddle-node bifurcation. With finding quantity of bifurcation parameter and rising it, we should try the patient's condition set in the stable branch of saddle-node bifurcation. In this case virus cells are controlled and patient is in the path of recuperation.

Dynamics of the model with delay (system 2)
In this section, we would like to investigate dynamics of system 2 with τ > 0.

Stability of equilibria
The first, we study the equilibrium E 0 in following theorem.
Theorem 1 if r ≤ a, then E 0 is locally stable; while if r > a then E 0 is unstable.
Proof By computing the characteristic equation of system 2 at E 0 , we have It complets the proof. Now, we consider the characteristic equation associated with the linearization of system 2 at E 1 Note that E 1 exists only if r > a, thus one root is ξ 1 = a − r < 0. Therefore, the dynamic of E 1 is depend on distribution of roots of the following equation

Theorem 2 The equilibrium point E 1 is locally asymptotically stable.
Proof By Lemma 2, the conclusion is true for τ = 0. We have to prove that all roots of g 2 (ξ ) have only negative real parts. Suppose that ξ = α + iω is a zero of g 2 (ξ ). After substituting in g 2 (ξ ), we obtain Therefore Note that α � = 0, by the above discussion, we assume α > 0. The right hand side convergent to zero but left hand side is perfectly elder of zero. Therefore, we have a contradiction or α < 0 and the proof is complete.
Theorems 1 and 2 and Lemma 1 show that if 14 holds, then E 0 is unstable and E 1 is stable, and E * exists. We now study the stability of E * . The characteristic equation at E * is By Lemma 2, when τ = 0, E * is asymptotically stable, i.e, all roots of the characteristic equation 33 have negative real parts. We want to prove E * is locally stable. With inverse process, we suppose that iω (ω > 0) is the root of G(ξ ), then we have which yields where Since g 1 (ŷ) = g 1 (y * ) = 0, thus c ′ 0 = 0, and Therefore, Eq. 33 has non purely imaginary root. On the other hand, we know g 1 (y * ) = 0 or c (1+dy * ) 2 = q. Therefore a 0 + b 0 = 0, and ξ = 0 is a simple zero of G(ξ ). Now we can state the following theorem.
Note that α � = 0. By the above discussion, we assume α > 0. By conditions of (1) and (2) , the right hand side is convergent to zero but left hand side is perfectly elder of zero. Hence, we have a contradiction or α < 0. This completes the proof.

Saddle-node bifurcation of system 2
In this subsection, we want to study codimension-one bifurcations of system 2. For this aim, we consider c as bifurcation parameter. By Remark 1, we know that E * exists if c = ( √ q + √ bd) 2 . Also, we know that E * is locally stable by Theorem 3, and codimension-one bifurcation can occur in system 2 at E * . Define c bif = ( √ q + √ bd) 2 . Now, we assume µ = c − c bif as bifurcation parameter and rewrite system 2 as follows Below we state the important theorem about existence saddle-node bifurcation of system 36. For this aim, we use center manifold theory of DDE, see "Appendix".

Theorem 4 System 36 has a saddle-node bifurcation at
Proof We consider the linearization of system 36 at E * new where The characteristic equation associated with system 37 is where G(ξ ) is defined by 33.
By Theorem 3, G(ξ ) has ξ = 0 as a root. Thus, G 2 (ξ ) has double zero roots. We want to obtain the center manifold associated with 37. To this end, we compute basis of a center subspace associated with 37 and adjoint system as follows Dadi and Alizade SpringerPlus (2016) 5:106 By using inner multiplication, we have With normalization ψ relation to φ, we obtain where then Now, suppose that local coordinates at center manifold is U = (u, µ) T . The terms of nonlinear system 36 and matrix B are Therefore, we have the following system by using the center manifold Define and (41) µ 1 := µ, µ 2 := Bµ + D.

Then
Also, by assumption u new = u − µ 2 A 2 , we have where Thus, studying dynamics of system 36 is equivalent to studying the following system Hence, by assumption of theorem, system 53 or other words system 36 has a saddlenode bifurcation.

Numerical simulation
By considering the following parameters: we have the saddle-node bifurcation in system 36, see Fig. 3.

Conclusion
An immunosuppressive infection model with discrete delays and without delay is considered. We have analyzed this model without delay in this paper and showed that the model has transcritical and saddle-node bifurcation at different parameters. We obtained a new equilibrium in our model with delay. Then, we have shown that this model undergoes saddle node bifurcation at this equilibrium. We then compute its normal form. Finally, the presented numerical simulations have demonstrated the correctness of the theoretical analysis.
(54) y * = 1.003, z * = 1.003, A = 1.012, B = 0.25 D = −1, µ = 0.249, µ 2 = −1 Fig. 3 Saddle-node bifurcation diagram Therefore, for every ϕ = (ϕ 1 , ϕ 2 ) and φ = (ϕ, ϕ 3 ) T ∈ C, we have and The stability of the trivial solution of the Eq. (55) can be studied by the DDE of the following form Substituting Y (t) = Ce t in the system (59), gives the following characteristic equation Obviously, the Eq. (60) always has one eigenvalue on the imaginary axis. We assume that this characteristic equation has m + 1 eigenvalues (counting multiplicity) on the imaginary axis and all other eigenvalues have negative real parts. Therefore, the space C can be split as C = P ⊕ Q where Q ⊂ C is infinite-dimensional stable subspace and P ⊂ C is an (m + 1)-dimensional center subspace tangent to the center manifold. We will denote a basis for P by the 3 × (m + 1) matrix Φ; the columns of Φ are the basis vectors. Also, we will consider the transpose of the Eq. (59) with (m +1)-dimensional center subspace P ′ . We will denote a basis for P ′ by the (m + 1) × 3 matrix Ψ ′ . Also, we define a new basis Ψ by Ψ =< Ψ ′ , Φ > −1 Ψ ′ which implies < Ψ , Φ >= I. This bilinear form is defined where This kind of basis Ψ can help us to decompose the space C and also reduce the Eq. (57) on the local center manifold W c loc which is defined by where h(z, F ) ∈ Q for each z and is a C r−1 function with respect to z. Moreover, z satisfies the following ordinary differential equation (59) Ẋ (t) = A 0 (µ)X(t) + A 1 (µ)X(t − τ ) µ = 0.