Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL)

John P. Elliott; John C. Morrison; James A. Bofill; John P. Elliott; John C. Morrison; James A. Bofill

doi:10.3934/publichealth.2016.2.348

AIMS Public Health

2016, Volume 3, Issue 2: 348-356. doi: 10.3934/publichealth.2016.2.348

Previous Article Next Article

Editorial

Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL)

1.
Banner University Medical Center, Phoenix, 1111 East McDowell Road, Phoenix, AZ 85005
2.
Department of Obstetrics and Gynecology, University of Mississippi Medical Center, 2500 North State Street, Jackson, Mississippi, 39216

Received: 04 December 2015 Accepted: 23 May 2016 Published: 30 May 2016

The U.S. Food and Drug Administration issued a drug safety communication on 05/30/2013 recommending “against prolonged use of magnesium sulfate to stop preterm labor (PTL) due to bone changes in exposed babies.” In September of 2013, The American Congress of Obstetrics and Gynecologists issued Committee Opinion No. 573 “ Magnesium Sulfate Use in Obstetrics” , which supports the short term use of MgSO₄ to prolong pregnancy (up to 48 hrs.) to allow for the administration of antenatal corticosteroids.” Are these pronouncements by respected organizations short sighted and will potentially result in more harm than good? The FDA safety communication focuses on bone demineralization (a few cases with fractures) with prolonged administration of MgSO₄ (beyond 5–7 days). It cites 18 case reports in the Adverse Event Reporting System with an average duration of magnesium exposure of 9.6 weeks (range 8–12 wks). Other epidemiologic studies showed transient changes in bone density which resolved in the short duration of follow up. Interestingly, the report fails to acknowledge the fact that these 18 fetuses were in danger of PTD and the pregnancy was prolonged by 9.6 weeks (e.g. extending 25 weeks to 34.6 wks), thus significantly reducing mortality and morbidity. Evidence does support the efficacy of MgSO₄ as a tocolytic medication. The decision to use magnesium, the dosage to administer, the duration of use, and alternative therapies are physician judgments. These decisions should be made based on a reasonable assessment of the risks of the clinical situation (PTL) and the treatments available versus the benefits of significantly prolonging pregnancy.

Keywords:

magnesium sulfate; tocolysis; preterm labor; FDA

Citation: John P. Elliott, John C. Morrison, James A. Bofill. Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL)[J]. AIMS Public Health, 2016, 3(2): 348-356. doi: 10.3934/publichealth.2016.2.348

Related Papers:

[1]	Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006
[2]	Yan Wang, Minmin Lu, Daqing Jiang . Viral dynamics of a latent HIV infection model with Beddington-DeAngelis incidence function, B-cell immune response and multiple delays. Mathematical Biosciences and Engineering, 2021, 18(1): 274-299. doi: 10.3934/mbe.2021014
[3]	Yu Yang, Gang Huang, Yueping Dong . Stability and Hopf bifurcation of an HIV infection model with two time delays. Mathematical Biosciences and Engineering, 2023, 20(2): 1938-1959. doi: 10.3934/mbe.2023089
[4]	Jiawei Deng, Ping Jiang, Hongying Shu . Viral infection dynamics with mitosis, intracellular delays and immune response. Mathematical Biosciences and Engineering, 2023, 20(2): 2937-2963. doi: 10.3934/mbe.2023139
[5]	Yuting Ding, Gaoyang Liu, Yong An . Stability and bifurcation analysis of a tumor-immune system with two delays and diffusion. Mathematical Biosciences and Engineering, 2022, 19(2): 1154-1173. doi: 10.3934/mbe.2022053
[6]	Kaushik Dehingia, Anusmita Das, Evren Hincal, Kamyar Hosseini, Sayed M. El Din . Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses. Mathematical Biosciences and Engineering, 2023, 20(11): 20025-20049. doi: 10.3934/mbe.2023887
[7]	A. M. Elaiw, A. S. Shflot, A. D. Hobiny . Stability analysis of general delayed HTLV-I dynamics model with mitosis and CTL immunity. Mathematical Biosciences and Engineering, 2022, 19(12): 12693-12729. doi: 10.3934/mbe.2022593
[8]	Juan Wang, Chunyang Qin, Yuming Chen, Xia Wang . Hopf bifurcation in a CTL-inclusive HIV-1 infection model with two time delays. Mathematical Biosciences and Engineering, 2019, 16(4): 2587-2612. doi: 10.3934/mbe.2019130
[9]	Hongying Shu, Wanxiao Xu, Zenghui Hao . Global dynamics of an immunosuppressive infection model with stage structure. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102. doi: 10.3934/mbe.2020111
[10]	Xuejuan Lu, Lulu Hui, Shengqiang Liu, Jia Li . A mathematical model of HTLV-I infection with two time delays. Mathematical Biosciences and Engineering, 2015, 12(3): 431-449. doi: 10.3934/mbe.2015.12.431

Abstract

1. Introduction

Recently, different kinds of viral infection such as HIV (human immunodeficiency virus), HCV (hepatitis C virus) and HBV (hepatitis B virus) have received many attention. Different mathematical models have been formulated to describe the dynamics of virus population in vivo ^{[1,2,3,4,5,6]}. The basic viral infection model of within-host can be described by the following three dimensional system ^[7]:

$\begin{equation} \begin{cases} \dfrac{dx}{dt} = s-\beta xv-d_{1}x, \\ \dfrac{dy}{dt} = \beta xv-d_{2}y,\\ \dfrac{dv}{dt} = ky-u v. \end{cases} \end{equation}$ $

(1.1)

Here $ x(t), y(t), $ and $ v(t) $ represent uninfected target cells, infected cells and free virus, respectively. Uninfected cells are produced at a constant rate $ s $, die at rate $ d_{1}x $, and become infected at rate $ \beta xv $. Infected cells are produced at rate $ \beta xv $ and die at rate $ by $. Free virus are produced from infected cells at rate $ ky $ and die at rate $ u v $. It was showed that, for system (1.1), there exist a critical threshold named as the basic reproduction number of virus $ R_{0} = \frac{\beta sk}{d_{1}d_{2}} $ to determine its global dynamical behavior ^[8].

During viral infections, the adaptive immune response is mediated by lymphocytes expressing antigen specific receptors, T and B lymphocytes, namely, humoral and cellular immunity. To study the population dynamics of immune response, Nowak et al. ^[9] introduced the CTL population into system (1.1) and obtain the following four dimensional system:

$\begin{equation} \begin{cases} \dfrac{dx}{dt} = s-\beta xv-d_{1}x, \\ \dfrac{dy}{dt} = \beta xv-d_{2}y-pyz,\\ \dfrac{dv}{dt} = ky-u v,\\ \dfrac{dz}{dt} = cyz-d_{3}z.\\ \end{cases} \end{equation}$ $

(1.2)

Here $ z(t) $ denotes CTLs, which is produced at rate $ cyz $ because of the stimulation of infected cells, and die at rate $ d_{3}z $. The infected cells are eliminated by CTLs at rate $ pyz $. Following ^[9], many authors present and develop mathematical models for the cell-mediated immune response ^{[10,11,12,13]} and the humoral immunity ^{[14,15,16,17]}.

Recently study indicates that the self-proliferation of immune cells can not be neglected besides the stimulation of infected cells. In order to mimic the spontaneous proliferation of CTLs, a logistic proliferation term for CTLs was incorporated in virus infection models in ^[18], where a rigorous mathematical analysis of the effect of self-proliferation of CTLs on the dynamics of viral infection is necessary. In order to understand the effect of self-proliferation and delayed activation of immune cells in a virus model, we propose the following virus infection model:

$\begin{eqnarray} \begin{cases} \begin{split} \frac{dx}{dt}& = s-\beta xv-d_{1}x ,\\ \frac{dy}{dt}& = \beta xv-d_{2}y-pyz,\\ \dfrac{dv}{dt}& = ky-u v,\\ \frac{dz}{dt}& = cy(t-\tau)z(t-\tau)+rz(1-\frac{z}{m})-d_{3}z. \end{split} \end{cases} \end{eqnarray}$ $

(1.3)

Here the logistic proliferation term $ rz(1-z/m) $ describes the self-proliferation of CTLs, in which parameter $ r $ denotes a per capita self-proliferation rate, and $ m $ means the capacity of CTLs population. Time delay term $ cy(t-\tau)z(t-\tau) $ represents a sequence of events such as antigenic activation and selection. It has been shown that time delays cannot be ignored in models for viral immune response including intracellular delay during viral infection ^{[19,20,21,22]} and time delay of CTLs stimulating proliferation ^[23,24].

Note that the turnover of virus is much faster than that of infection cells within-host ^[25,26]. A plausible quasi steady-state assumption is proposed to mimic the fast time-scale^[27,28]. In other words, $ v $ can be replaced by $ \frac{ky}{u} $ in (1.3). Let $ \hat{\beta} = \frac{\beta k}{u} $ and also note as $ \beta $ to simplify the parameter, system (1.3) can be written as:

$\begin{eqnarray} \begin{cases} \begin{split} \frac{dx}{dt}& = s-\beta xy-d_{1}x ,\\ \frac{dy}{dt}& = \beta xy-d_{2}y-pyz,\\ \frac{dz}{dt}& = cy(t-\tau)z(t-\tau)+rz(1-\frac{z}{m})-d_{3}z. \end{split} \end{cases} \end{eqnarray}$ $

(1.4)

This paper is organized as follows. In section 2, some preliminary results are obtained, including non-negativity and boundedness of the solution of system (1.4), the existence of equilibria under different self-proliferation intensities of CTLs. Section 3 investigate the local and global stability of the boundary equilibria. Section 4 study the effects of delayed activation of immune cells on the existence of Hopf bifurcation. The direction and stability of Hopf bifurcation is investigated in section 5. Some numerical simulations are given to quantify the impact of self-proliferation and delayed activation of CTLs in section 6. Finally, some conclusions and discusses are presented in section 7.

2. Preliminary results

In this section, we first discuss the non-negativity and boundedness of solutions of system (1.4). For $ \tau > 0 $, let $ \mathcal{C} = \mathcal{C}([-\tau, 0], \mathbb{R}^{3}_{+}) $ denote the Banach space of continuous function mapping the interval $ [-\tau, 0] $ into $ \mathbb{R}^{3}_{+} $ with the topology of uniform convergence. The initial conditions are given by

$\begin{equation} x(\xi) = \phi_{1}(\xi),y(\xi) = \phi_{2}(\xi),z(\xi) = \phi_{3}(\xi) \end{equation}$ $

(2.1)

with $ \phi_{i}(\xi)\geq 0 $, $ \xi\in[-\tau, 0] $ and $ \phi_{i}(0) > 0 $ $ (i = 1, 2, 3) $.

Theorem 2.1. The solutions of system (1.4) satisfying the initial conditions (2.1) are non-negative and ultimately bounded.

Proof. We first define the right-hand side function of system (1.4) as

$ G(t,K(t)) = \left(

$\begin{array}{c} s-\beta xy-d_{1}x \\ \beta xy-d_{2}y-pyz \\ cy(t-\tau)z(t-\tau)+rz(1-\dfrac{z}{m})-d_{3}z \\ \end{array}$ \right) , $

where $ K(t) = (K_{1}(t), K_{2}(t), K_{3}(t))^{T} $ and $ K_{1}(t) = x $, $ K_{2}(t) = y $, $ K_{3}(t) = z. $ It is obvious that the function $ G(t, K(t)) $ is locally Lipschitz and by the standard theory of functional differential equation, we know that there exists a unique solution for a given initial conditions. To prove the non-negativity of solutions of system (1.4), we first prove the non-negativity over the time interval $ [0, \tau] $. Considering the right-hand side functions of system (1.4) over the time interval $ [0, \tau] $, we have

$\begin{equation*} G(t,K(t)) = \left( \begin{array}{c} G_{1}(t,K(t)) \\ G_{2}(t,K(t)) \\ G_{3}(t,K(t)) \\ \end{array} \right)_{K_{i}(t) = 0} = \left( \begin{array}{c} s \\ 0 \\ cy(t-\tau)z(t-\tau) \\ \end{array} \right) = \left( \begin{array}{c} s \\ 0 \\ c\phi_{2}(\xi)\phi_{3}(\xi) \\ \end{array} \right), \end{equation*}$ $

where $ \xi\in[-\tau, 0] $. Note that the positivity of parameters and nonnegativity of initial functions. It can be seen from the above expression that each component $ G_{i}(t, K(t))\geq0 $. It follows that the solutions of system (1.4) remain nonnegative in $ [0, \tau] $. Similarly, we can repeat the process over $ [\tau, 2\tau] $ and so on by using the method of steps ^[29], it can be proved that for any finite interval $ [0, t] $, the solutions of system (1.4) remains non-negative.

Now we consider the boundeness of the solutions. Define a new variable $ X(t) = x(t)+y(t)+\dfrac{p}{c}z(t+\tau) $, and let $ d = \min\{1, d_{1}, d_{2}\} $. By the non-negativity of the solutions of system (1.4), we have

$\begin{align} \dfrac{dX}{dt}& = s-d_{1}x-d_{2}y+\frac{pr}{c}z(t+\tau)\Big(1-\frac{z(t+\tau)}{m}\Big)-d_{3}z(t+\tau){}\\ {}& = s-d_{1}x-d_{2}y-\frac{p}{c}z(t+\tau)+\frac{p}{c}z(t+\tau)\Big(1+r-\frac{rz(t+\tau)}{m}\Big)-d_{3}z(t+\tau){}\\ {}&\leq s-d_{1}x-d_{2}y-\frac{p}{c}z(t+\tau)+\frac{pm(1+r)^2}{4rc}{}\\ {}&\leq s+\frac{pm(1+r)^2}{4rc}-dX(t). \end{align}$ $

(2.2)

Taking $ M = s+\dfrac{pm(1+r)^2}{4rc} $, we know $ \limsup_{t\rightarrow\infty}X(t)\leq\dfrac{M}{d} $. So the solutions of system(1.4) are ultimately bounded. The equilibria of (1.4) are the solutions of the following algebraic equations:

$\begin{eqnarray} \begin{cases} \begin{split} & s-\beta xy-d_{1}x = 0 ,\\ & \beta xy-d_{2}y-pyz = 0,\\ & cyz+rz(1-\frac{z}{m})-d_{3}z = 0. \end{split} \end{cases} \end{eqnarray}$ $

(2.3)

It is easy to see that system (1.4) always has infection-free equilibrium $ E_{0}^{yz} = (x_{0}, 0, 0) $, where $ x_{0} = \dfrac{s}{d_{1}} $. According to the definition in ^[30], we obtain the basic reproduction number of virus $ R_{0} = \dfrac{\beta s}{d_{1}d_{2}} $. Now we define $ x_{1} = \dfrac{d_{2}}{\beta}, y_{1} = \dfrac{d_{1}(R_{0}-1)}{\beta}, x_{2} = \dfrac{s}{d_{1}}, z_{2} = \dfrac{m(r-d_{3})}{r} $. Based on (2.3), after some simple calculations, it is easy to get the following results.

Proposition 2.2. (i) Suppose that $ R_{0} > 1 $. The immunity-inactivated infection equilibrium $ E_{0}^{z} = (x_{1}, y_{1}, 0) $ always exists. Especially, besides $ E_{0}^{z} $, an infection-free but immunity-activated equilibrium $ E_{0}^{y} = (x_{2}, 0, z_{2}) $ will appear if $ r > d_{3} $. (ii) Suppose that $ 0\leq r\leq d_{3} $. If $ R_{0} > 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $, system (1.4) has a unique immunity-activated infection equilibrium $ E_{1}^{*} = (x_{1}^{*}, y_{1}^{*}, z_{1}^{*}) $, where

$ x_{1}^{*} = \dfrac{d_{2}+pz_{1}^{*}}{\beta},\; \; \; \; \; \; \; y_{1}^{*} = \dfrac{rz_{1}^{*}+m(d_{3}-r)}{mc}, $

and $ z_{1}^{*} $ is the unique positive root of the quadratic equation $ A_{1}z^{2}+B_{1}z+C_{1} = 0 $, in which

$ A_{1} = -pr\beta, B_{1} = -(d_{2}r\beta+pmcd_{1}+\beta pm(d_{3}-r)), C_{1} = mcd_{1}d_{2}(R_{0}-1)-\beta md_{2}(d_{3}-r). $

(iii) Suppose that $ r > d_{3} $. If $ R_{0} > 1+\dfrac{pm(r-d_{3})}{rd_{2}} $, system (1.4) has a unique immunity-activated infection equilibrium $ E_{2}^{*} = (x_{2}^{*}, y_{2}^{*}, z_{2}^{*}) $, where

$ x_{2}^{*} = \frac{s}{d_{1}+\beta y_{2}^{*}},\; \; \; \; \; \; \; z_{2}^{*} = \frac{m(cy_{2}^{*}+r-d_{3})}{r}, $

and $ y_{2}^{*} $ is the unique positive root of the quadratic equation $ A_{2}y^{2}+B_{2}y+C_{2} = 0 $, in which

$ A_{2} = pmc\beta, B_{2} = pm\beta (r-d_{3})+d_{2}r\beta+pmcd_{1}, C_{2} = -rd_{1}d_{2}(R_{0}-1)+pmd_{1}(r-d_{3}). $

3. Global stability analysis of boundary equilibria

In order to study the stability of equilibrium, we linearize system (1.4) about one equilibrium $ E^* = (x^{*}, y^{*}, z^{*}) $ and obtain the following linear system

$\begin{eqnarray} \begin{cases} \begin{split} \frac{dx}{dt}& = (-\beta y^{*}-d_{1})x-\beta x^{*}y,\\ \frac{dy}{dt}& = \beta y^{*}x+(\beta x^{*}-d_{2}-pz^{*})y-py^{*}z,\\ \frac{dz}{dt}& = cz^{*}y(t-\tau)+cy^{*}z(t-\tau)+(r-d_{3}-\dfrac{2rz^{*}}{m})z. \end{split} \end{cases} \end{eqnarray}$ $

(3.1)

When $ r = 0 $, it follows from the study of Michael Y. Li and Hongying Shu ^[23] that time delay $ \tau $ does not change the stability of the boundary equilibria.

When $ r > 0 $, by using the linear system (3.1) and constructing Lyapunov function, we can obtain the following results.

Proposition 3.1. Suppose that $ 0 < r < d_{3} $. Then we have

(i) The infection-free equilibrium $ E_{0}^{yz} $ is globally asymptotically stable if $ R_{0} < 1 $, and it is unstable when $ R_{0} > 1 $.

(ii) The immunity-inactivated infection equilibrium $ E_{0}^{z} $ is globally asymptotically stable if $ 1 < R_{0} < 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $, and it is unstable when $ R_{0} > 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $.

Proof. (ⅰ) By the linear system (3.1), we can obtain the characteristic equation of system(1.4) at $ E_{0}^{yz} $ as

$\begin{equation} (\lambda+d_{1})(\lambda-(r-d_{3}))\big(\lambda+d_{1}d_{2}(1-R_{0})\big) = 0. \end{equation}$ $

(3.2)

So the eigenvalues are

$ \lambda_{1} = -d_{1} \lt 0,\lambda_{2} = r-d_{3} \lt 0,\lambda_{3} = d_{1}d_{2}(R_{0}-1). $

As a result, if $ R_{0} < 1 $, the infection-free equilibrium $ E_{0}^{yz} $ is locally asymptotically stable, and $ E_{0}^{yz} $ is unstable if $ R_{0} > 1 $.

In order to obtain the global stability of equilibrium $ E_{0}^{yz} $, we consider a Lyapunov function given by

$ L_{1} = x-x_{0}-x_{0}\ln\frac{x}{x_{0}}+y+\frac{p}{c}z+p\int_{-\tau}^{0}y(t+\theta)z(t+\theta)d\theta. $

Taking the time derivative of $ L_{1} $ along the solution of system (1.4), we have

$\begin{align*} \frac{dL_{1}}{dt}& = (1-\frac{x_{0}}{x})(s-\beta xy-d_{1}x)+\beta xy-d_{2}y-pyz\\ &+\frac{p}{c}[cy(t-\tau)z(t-\tau)+rz(1-\dfrac{z}{m})-d_{3}z]+pyz-py(t-\tau)z(t-\tau)\\ & = -\frac{d_{1}(x-x_{0})^{2}}{x}+(\beta x_{0}-d_{2})y+\frac{p}{c}(r-d_{3})z-\dfrac{prz^{2}}{cm}\\ & = -\frac{d_{1}(x-x_{0})^{2}}{x}+d_{2}(R_{0}-1)y+\frac{p}{c}(r-d_{3})z-\dfrac{prz^{2}}{cm}. \end{align*}$ $

Notice that $ r < d_{3} $ and $ R_{0} < 1 $. We have $ L_{1}^{'}\leq0 $ for all $ x(t), y(t), z(t) > 0 $, and $ L_{1}^{'} = 0 $ only if $ x = x_{0} $, $ y = 0 $ and $ z = 0 $. It can be verified that the maximal compact invariant set in $ {L_{1}^{'} = 0} $ is the singleton $ E_{0}^{yz} $. By the LaSalle's invariance principle, we know the infection-free equilibrium $ E_{0}^{yz} $ is globally asymptotically stable.

(ⅱ) The characteristic equation of system(1.4) at $ E_{0}^{z} $ is

$\begin{equation} (\lambda^{2}+d_{1}R_{0}\lambda+\beta^{2}x_{1}y_{1})\big(\lambda-(r-d_{3}+cy_{1}e^{-\lambda\tau})\big) = 0. \end{equation}$ $

(3.3)

When $ \tau = 0 $ and $ R_{0} < 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $, the roots of (3.3) are negative. When $ \tau > 0 $, the local stability of equilibrium $ E_{0}^{z} $ is solely determined by

$\begin{equation} \lambda-(r-d_{3}+cy_{1}e^{-\lambda\tau}) = 0. \end{equation}$ $

(3.4)

Let $ \lambda = i\omega(\tau) (\omega(\tau) > 0) $ be the root of Eq (3.4). Separating real and imaginary parts yields

$\begin{equation} \left\{ \begin{aligned} & d_{3}-r = cy_{1}\cos(\omega\tau),\\ & \omega = cy_{1}\sin(\omega\tau). \end{aligned} \right. \end{equation}$ $

(3.5)

Squaring and adding Eq (3.5) gives

$ \omega^{2}-(cy_{1})^{2}+(d_{3}-r)^{2} = 0. $

Note that

$ (cy_{1})^{2}-(d_{3}-r)^{2} = \bigg(\dfrac{cd_{1}(R_{0}-1)}{\beta(d_{3}-r)}\bigg)^{2}-1 \lt 0 $

if $ R_{0} < 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $. Then we konw that (3.3) has no root that can across the imaginary axis, which indicate that the immunity-inactivated infection equilibrium $ E_{0}^{z} $ is locally asymptotically stable when $ R_{0} < 1+\dfrac{\beta(d_{3}-r)}{cd_{1}} $.

In order to obtain the global stability of $ E_{0}^{z} $, we construct the following Lyapunov functions:

$ L_{2} = x-x_{1}-x_{1}\ln\frac{x}{x_{1}}+y-y_{1}-y_{1}\ln\frac{y}{y_{1}}+\frac{p}{c}z+p\int_{-\tau}^{0}y(t+\theta)z(t+\theta)d\theta. $

Taking the time derivative of $ L_{2} $ along the solution of system(1.4), we get

$\begin{align*} \frac{dL_{2}}{dt}& = (1-\frac{x_{1}}{x})(s-\beta xy-d_{1}x)+(1-\frac{y_{1}}{y})(\beta xy-d_{2}y-pyz)\\ &+\frac{p}{c}\big(cy(t-\tau)z(t-\tau)+rz(1-\dfrac{z}{m})-d_{3}z\big)+pyz-py(t-\tau)z(t-\tau)\\ & = s-d_{1}x-s\frac{x_{1}}{x}+\beta x_{1}y+d_{1}x_{1}-d_{2}y-\beta xy_{1}+d_{2}y_{1}+py_{1}z+\frac{p}{c}(r-d_{3})z-\frac{prz^{2}}{cm}. \end{align*}$ $

Using $ s = \beta x_{1}y_{1}+d_{1}x_{1} $ and $ d_{2} = \beta x_{1} $, we have

$\begin{align*} \frac{dL_{2}}{dt}& = (d_{1}x_{1}+\beta x_{1}y_{1})(2-\frac{x}{x_{1}}-\frac{x_{1}}{x})+p(y_{1}-\frac{d_{3}-r}{c})z-\frac{prz^{2}}{cm}\\ & = (d_{1}x_{1}+\beta x_{1}y_{1})(2-\frac{x}{x_{1}}-\frac{x_{1}}{x})+\frac{p d_{1}}{\beta}(R_{0}-1-\frac{\beta(d_{3}-r)}{cd_{1}})z-\frac{prz^{2}}{cm}. \end{align*}$ $

Using the LaSalle's invariance principle, we can obtain that $ E_{0}^{z} $ is globally asymptotically stable.

Proposition 3.2. Suppose that $ r = d_{3} $. Then we have

(i) The infection-free equilibrium $ E_{0}^{yz} $ is globally asymptotically stable if $ R_{0} < 1 $, and it is unstable when $ R_{0} > 1 $.

(ii) The immunity-inactivated infection equilibrium $ E_{0}^{z} $ is unstable as long as it appears, i.e., $ R_{0} > 1 $.

Proof. (ⅰ) The characteristic equation of system (1.4) at $ E_{0}^{yz} $ is

$\begin{equation} \lambda(\lambda+d_{1})\big(\lambda+d_{1}d_{2}(1-R_{0})\big) = 0. \end{equation}$ $

(3.6)

So the eigenvalues are

$ \lambda_{1} = 0,\lambda_{2} = -d_{1} \lt 0,\lambda_{3} = d_{1}d_{2}(R_{0}-1). $

So if $ R_{0} < 1 $, the infection-free equilibrium $ E_{0}^{yz} $ is locally asymptotically stable and if $ R_{0} > 1 $, the equilibrium $ E_{0}^{yz} $ is unstable.

In order to obtain the global stability of $ E_{0}^{yz} $, we construct the following Lyapunov functions:

$ L_{3} = x-x_{0}-x_{0}\ln\frac{x}{x_{0}}+y+\frac{p}{c}z+p\int_{-\tau}^{0}y(t+\theta)z(t+\theta)d\theta. $

Taking the time derivative of $ L_{3} $ along the solution of system (1.4), we have

$\begin{align*} \frac{dL_{3}}{dt}& = (1-\frac{x_{0}}{x})(s-\beta xy-d_{1}x)+\beta xy-d_{2}y-pyz\\ &+\frac{p}{c}(cy(t-\tau)z(t-\tau)-\dfrac{rz^{2}}{m})+pyz-py(t-\tau)z(t-\tau)\\ & = -\frac{d_{1}(x-x_{0})^{2}}{x}+(\beta x_{0}-d_{2})y-\dfrac{prz^{2}}{cm}\\ & = -\frac{d_{1}(x-x_{0})^{2}}{x}+d_{2}(R_{0}-1)y-\dfrac{prz^{2}}{cm}. \end{align*}$ $

Using the LaSalle's invariance principle, we can obtain that $ E_{0}^{yz} $ is globally asymptotically stable.

(ⅱ)The characteristic equation of system (1.4) at $ E_{0}^{z} $ is

$\begin{equation} (\lambda^{2}+d_{1}R_{0}\lambda+\beta^{2}x_{1}y_{1})(\lambda-cy_{1}e^{-\lambda\tau}) = 0. \end{equation}$ $

(3.7)

It can be showed that there exist positive real root for the characteristic Eq (4.3), which indicates that the infection-free equilibrium $ E_{0}^{z} $ is unstable. This completes the proof.

Proposition 3.3. Suppose that $ r > d_{3} $. Then we have

(i) The infection-free equilibrium $ E_{0}^{yz} $ is unstable.

(ii) The immunity-inactivated infection equilibrium $ E_{0}^{z} $ is unstable.

(iii) If $ R_{0} < 1+\dfrac{pm(r-d_{3})}{rd_{2}} $, the infection-free equilibrium $ E_{0}^{y} $ is locally asymptotically stable, moreover if $ R_{0} < 1 $, the infection-free equilibrium $ E_{0}^{y} $ is globally asymptotically stable; if $ R_{0} > 1+\dfrac{pm(r-d_{3})}{rd_{2}} $, $ E_{0}^{y} $ is unstable.

Proof. (ⅰ) The characteristic equation of system (1.4) at $ E_{0}^{yz} $ is

$\begin{equation} (\lambda+d_{1})(\lambda-(r-d_{3}))\big(\lambda+d_{1}d_{2}(1-R_{0})\big) = 0. \end{equation}$ $

(3.8)

So the eigenvalues are

$ \lambda_{1} = -d_{1},\lambda_{2} = r-d_{3},\lambda_{3} = d_{1}d_{2}(R_{0}-1). $

It follows from $ \lambda_{2} = r-d_{3} > 0 $ that the infection-free equilibrium $ E_{0}^{yz} $ is always unstable. So we can easily get the infection-free equilibrium $ E_{0}^{yz} $ is always unstable when $ \tau > 0 $.

(ⅱ) The characteristic equation of system (1.4) at $ E_{0}^{z} $ is

$\begin{equation} (\lambda^{2}+d_{1}R_{0}\lambda+\beta^{2}x_{1}y_{1})\big(\lambda-(r-d_{3}+cy_{1}e^{-\lambda\tau})\big) = 0. \end{equation}$ $

(3.9)

Note that $ r > d_{3} $. It can be shown that there exist positive real root for the characteristic equation above, which indicates that the immunity-inactivated infection equilibrium $ E_{0}^{z} $ is unstable.

(ⅲ) The characteristic equation of system (1.4) at $ E_{0}^{y} $ is

$\begin{equation} (\lambda+d_{1})(\lambda+r-d_{3})(\lambda-\frac{s\beta}{d_{1}}+d_{2}+\dfrac{pm(r-d_{3})}{r}) = 0. \end{equation}$ $

(3.10)

The eigenvalues are

$ \lambda_{1} = -d_{1},\lambda_{2} = -(r-d_{3}),\lambda_{3} = \frac{s\beta}{d_{1}}-d_{2}-\dfrac{pm(r-d_{3})}{r}. $

Note that

$ \frac{s\beta}{d_{1}}-d_{2}-\dfrac{pm(r-d_{3})}{r} = d_{2}(R_{0}-1)-\dfrac{pm(r-d_{3})}{r} \lt 0\Leftrightarrow R_{0} \lt 1+\dfrac{pm(r-d_{3})}{rd_{2}}. $

Thus, the infection-free equilibrium $ E_{0}^{y} $ is locally asymptotically stable if $ R_{0} < 1+\dfrac{pm(r-d_{3})}{rd_{2}} $, and $ E_{0}^{y} $ is unstable if $ R_{0} > 1+\dfrac{pm(r-d_{3})}{rd_{2}} $.

In order to obtain the global stability of $ E_{0}^{y} $, we construct the following Lyapunov function:

$ L_{4} = x-x_{2}-x_{2}\ln\frac{x}{x_{2}}+y+\frac{p}{c}(z-z_{2}-z_{2}\ln\frac{z}{z_{2}})+\frac{p}{c}z+p\int_{-\tau}^{0}y(t+\theta)z(t+\theta)d\theta. $

Taking the time derivative of $ L_{4} $ along the solution of system (1.4) and using a computation process similar to that of Theorem 3.9, we have

$\begin{align*} \frac{dL_{4}}{dt} & = (1-\frac{x_{2}}{x})(s-\beta xy-d_{1}x)+\beta xy-d_{2}y-pyz \\ & \;\;\;\; +\frac{p}{c}(1-\frac{z_{2}}{z})\bigg(cy(t-\tau)z(t-\tau)+rz(1-\dfrac{z}{m})-d_{3}z\bigg)\\& = -\frac{d_{1}(x-x_{2})^{2}}{x}+d_{2}(R_{0}-1)y-\dfrac{pz_{2}y(t-\tau)z(t-\tau)}{z}-\frac{rp}{cm}(z-z_{2})^{2}. \end{align*}$ $

Thus, $ \dfrac{dL_{4}}{dt}\leq0 $ if $ R_{0} < 1 $, and $ \dfrac{dL_{4}}{dt} = 0 $ only if $ x = x_{2}, y = 0 $ and $ z = z_{2} $, i.e., the maximal invariant subset in $ \{(x, y, z): \dfrac{dL_{4}}{dt}\Big|_{(1.4)} = 0\} $ is the singleton $ \{E_{0}^{y}\} $. As a result, $ E_{0}^{y} $ is globally asymptotically stable based on the LaSalle's invariance principle. This complete the proof.

Remark 3.4. From Proposition 3.1 to Proposition 3.3, we can see that the delay $ \tau $ does not affect the stability of infection-free equilibrium $ E_{0}^{yz} $, $ E_{0}^{y} $ and immune-unactivated $ E_{0}^{z} $ equilibrium.

4. Stability of positive equilibrium and Hopf bifurcation

In this section, we take the discrete delay $ \tau $ as a bifurcation parameter and show that when the positive equilibrium $ E_{2}^{*} $ loses its stability and a Hopf bifurcation appears when the delay $ \tau $ passes through a critical value. We point out there also exists a Hopf bifurcation at positive equilibrium $ E_{1}^{*} $ as the delay $ \tau $ passes through a critical value and the proof is similar.

The characteristic equation of system (1.4) at $ E_{2}^{*} $ is

$\begin{equation} \lambda^{3}+B_{1}\lambda^{2}+B_{2}\lambda+B_{3}+(C_{1}\lambda^{2}+C_{2}\lambda+C_{3})e^{-\lambda\tau} = 0, \end{equation}$ $

(4.1)

where

$\begin{align*} B_{1}& = \frac{2rz_{2}^{*}}{m}-r+d_{3}+d_{1}+\beta y_{2}^{*},\; \; \; \; \; \; \; \; \; B_{2} = (d_{1}+\beta y_{2}^{*})(\frac{2rz_{2}^{*}}{m}-r+d_{3})+\beta^{2} x_{2}^{*}y_{2}^{*},\\ B_{3}& = (\frac{2rz_{2}^{*}}{m}-r+d_{3})\beta^{2} x_{2}^{*}y_{2}^{*},\; \; \; \; \; \; \; \; \; \; \; \; \; C_{1} = -cy_{2}^{*},\\ C_{2}& = -(d_{1}+\beta y_{2}^{*})cy_{2}^{*}+pcy_{2}^{*}z_{2}^{*},\; \; \; \; \; \; \; \; \; \; C_{3} = -cy_{2}^{*}\beta^{2} x_{2}^{*}y_{2}^{*}+(d_{1}+\beta y_{2}^{*})pcy_{2}^{*}z_{2}^{*}. \end{align*}$ $

When $ \tau = 0, $ the characteristic equation of system (1.4) at $ E_{2}^{*} $ is

$\begin{equation*} \label{4.3} \lambda^{3}+A_{1}\lambda^{2}+A_{2}\lambda+A_{3} = 0, \end{equation*}$ $

where

$\begin{align*} A_{1}&\equiv B_{1}+C_{1} = d_{1}+\beta y_{2}^{*}+cy_{2}^{*}+r-d_{3} \gt 0, \\A_{2}&\equiv B_{2}+C_{2} = (d_{1}+\beta y_{2}^{*})(cy_{2}^{*}+r-d_{3})+\beta^{2} x_{2}^{*}y_{2}^{*}+pcy_{2}^{*}z_{2}^{*} \gt 0,\\ A_{3}&\equiv B_{3}+C_{3} = pcy_{2}^{*}z_{2}^{*}(d_{1}+\beta y_{2}^{*})+(cy_{2}^{*}+r-d_{3})\beta^{2} x_{2}^{*}y_{2}^{*} \gt 0. \end{align*}$ $

After some calculations, we have

$\begin{align*} A_{1}A_{2}-A_{3}& = (cy_{2}^{*}+r-d_{3})\big(pcy_{2}^{*}z_{2}^{*}\\&+(d_{1}+\beta y_{2}^{*})^{2}\big) +(d_{1}+\beta y_{2}^{*})\big(\beta^{2} x_{2}^{*}y_{2}^{*}+(cy_{2}^{*}+r-d_{3})^{2}\big) \gt 0. \end{align*}$ $

It follows from the Routh-Hurwitz criterion that $ E_{2}^{*} $ is locally asymptotically stable when time delay is absent.

When $ \tau > 0, $ putting $ \lambda = i\omega $ into characteristic Eq (4.1) and separating real and imaginary parts, we have

$\begin{equation} B_{1}\omega^{2}-B_{3} = (C_{3}-C_{1}\omega^{2})\cos(\omega\tau)+C_{2}\omega\sin(\omega\tau), \end{equation}$ $

(4.2)

$\begin{equation} \omega^{3}-B_{2}\omega = -(C_{3}-C_{1}\omega^{2})\sin(\omega\tau)+C_{2}\omega\cos(\omega\tau). \end{equation}$ $

(4.3)

Since $ \sin^{2}(\omega\tau)+\cos^{2}(\omega\tau) = 1 $, squaring and adding the two Eqs (4.2) and (4.3), we have

$\begin{equation} \omega^{6}+p_{1}\omega^{4}+p_{2}\omega^{2}+p_{3} = 0, \end{equation}$ $

(4.4)

where

$\begin{align*} p_{1}& = B_{1}^{2}-2B_{2}-C_{1}^{2},\\ p_{2}& = B_{2}^{2}+2C_{1}C_{3}-C_{2}^{2}-2B_{1}B_{3},\\ p_{3}& = B_{3}^{2}-C_{3}^{2}.\\ \end{align*}$ $

Let $ u = \omega^{2} $, Eq (4.4) becomes

$\begin{equation} G(u) = u^{3}+p_{1}u^{2}+p_{2}u+p_{3} = 0. \end{equation}$ $

(4.5)

If Eq (4.5) has a positive real root $ u $, the characteristic Eq (4.1) has a purely imaginary root $ i\omega = i\sqrt{u} $; otherwise, Eq (4.1) has no purely imaginary root.

Note that

$\begin{equation} G^{'}(u) = 3u^{2}+2p_{1}u+p_{2}. \end{equation}$ $

(4.6)

Let

$ \Delta = p_{1}^{2}-3p_{2}. $

Then we know that $ G(u) $ is monotonically increasing if $ \Delta\leq0 $, which indicates that Eq (4.5) has no positive root when $ p_{3}\geq0 $ and $ \Delta\leq0 $.

If $ \Delta > 0 $, the function $ G(u) $ has two critical points

$ u^{*} = \frac{-p_{1}+\sqrt{\Delta}}{3},\; \; \; \; \; \; \; u^{**} = \frac{-p_{1}-\sqrt{\Delta}}{3}. $

Thus we know that Eq (4.5) has unique positive root $ u_{0} $ with $ G'(u_{0}) > 0 $ if one of the following two conditions hold:

$ (H_{1})\; $$ \Delta > 0 $, $ p_{3} < 0 $, $ u^{**} < 0 $ and $ u^{*} > 0 $;

$ (H_{2}) \; $$ \Delta > 0 $, $ p_{3} < 0 $, $ u^{*} > 0 $ and $ G(u^{**}) < 0 $.

Equation (4.5) has two positive roots $ u_{1} < u_{2} $ with $ G'(u_{1}) < 0 $ and $ G'(u_{2}) > 0 $ if the following assumption is satisfied:

$ (H_{3})\; $$ \Delta > 0 $, $ p_{3} > 0 $, $ u^{*} > 0 $ and $ G(u^{*}) < 0 $.

Let $ \omega_{k} = \sqrt{u_{k}}, k = 0, 1, 2 $. It follows from Eqs (4.2) and (4.3) that the value of $ \tau $ associated with the purely imaginary root $ i\omega_{k} $ should satisfy

$\begin{equation*} (B_{1}\omega_{k}^{2}-B_{3})(C_{3}-C_{1}\omega_{k}^{2})+ (\omega_{k}^{3}-B_{2}\omega_{k})C_{2}\omega_{k} = \big((C_{3}-C_{1}\omega_{k}^{2})^{2}+C_{2}^{2}\omega_{k}^{2}\big)\cos(\omega_{k}\tau). \end{equation*}$ $

For $ k = 0, 1, 2 $, we define

$\begin{equation} \tau^{k}_{n} = \frac{1}{\omega_{k}}\arccos\Big(\frac{(B_{1}\omega_{k}^{2}-B_{3})(C_{3}-C_{1}\omega_{k}^{2}) +(\omega_{k}^{3}-B_{2}\omega_{k})C_{2}\omega_{k}}{(C_{3}-C_{1}\omega_{k}^{2})^{2}+C_{2}^{2} \omega_{k}^{2}}\Big)+\frac{2n\pi}{\omega_{k}},n = 0,1,2\cdots. \end{equation}$ $

(4.7)

Then at increasing sequences of $ \tau $ values,

$\begin{equation*} \begin{aligned} &\tau_{0}^{0} \lt \tau_{1}^{0} \lt \tau_{2}^{0} \lt \cdot\cdot\cdot \lt \tau_{n}^{0}\cdot\cdot\cdot,\\ &\tau_{0}^{1} \lt \tau_{1}^{1} \lt \tau_{2}^{1} \lt \cdot\cdot\cdot \lt \tau_{n}^{1}\cdot\cdot\cdot,\\ &\tau_{0}^{2} \lt \tau_{1}^{2} \lt \tau_{2}^{2} \lt \cdot\cdot\cdot \lt \tau_{n}^{2}\cdot\cdot\cdot, \end{aligned} \end{equation*}$ $

Eq (4.1) has purely imaginary roots $ i\omega_{k}, k = 0, 1, 2 $.

Now we consider the transversality conditions associated with Hopf bifurcation. Substituting $ \lambda(\tau) $ into Eq (4.1) and differentiating the resulting equation in $ \tau $, we obtain

$\begin{align*} \bigg(\frac{d\lambda}{d\tau}\bigg)^{-1}& = \frac{3\lambda^{2}+2B_{1}\lambda+B_{2}+(2C_{1}\lambda+C_{2})e^{-\lambda\tau}+(C_{1}\lambda^{2}+C_{2}\lambda+C_{3})\cdot(-\tau)e^{-\lambda\tau}}{\lambda e^{-\lambda\tau}(C_{1}\lambda^{2}+C_{2}\lambda+C_{3})}\\ & = \frac{(3\lambda^{2}+2B_{1}\lambda+B_{2})e^{\lambda\tau}}{\lambda(C_{1}\lambda^{2}+C_{2}\lambda+C_{3})} +\frac{2C_{1}\lambda+C_{2}}{\lambda(C_{1}\lambda^{2}+C_{2}\lambda+C_{3})}-\frac{\tau}{\lambda}. \end{align*}$ $

which indicates that

$\begin{align*} \text{sign}\bigg\{\frac{d(\text{Re}\lambda)}{d\tau}\bigg\}\Big|_{\lambda = i\omega_{k}} & = \text{sign}\bigg\{\text{Re}\Big(\frac{d\lambda}{d\tau}\Big)^{-1}\bigg\}\Big|_{\lambda = i\omega_{k}}\\ & = \text{sign}\bigg\{\frac{3\omega_{k}^{4}+2(B_{1}^{2}-2B_{2}-C_{1}^{2})\omega_{k}^{2}+(B_{2}^{2}+2C_{1}C_{3}-C_{2}^{2}-2B_{1}B_{3})}{(C_{3} -C_{1}\omega_{k}^{2})^{2}+C_{2}^{2}\omega_{k}^{2}}\bigg\}\\ & = \text{sign}\bigg\{\frac{G^{'}(\omega_{k}^{2})}{(C_{3}-C_{1}\omega_{k}^{2})^{2}+C_{2}^{2}\omega_{k}^{2}}\bigg\}. \end{align*}$ $

Since $ G'(u_{0}) > 0 $, $ G'(u_{1}) < 0 $ and $ G'(u_{2}) > 0 $, then for $ n = 0, 1, 2, ..., $ we have

$ \text{sign}\Big\{\frac{d(\text{Re}\lambda)}{d\tau}\Big|_{\tau = \tau_{n}^{k}}\Big\} = \text{sign}\Big\{\text{Re}\Big(\frac{d\lambda}{d\tau}\Big)^{-1}\Big|_{\tau = \tau_{n}^{k}}\Big\} \gt 0, (k = 0,2), $

and

Then we know that at each $ \tau_{n}^{0} $ and $ \tau_{n}^{2} $, $ n = 0, 1, 2, \cdot\cdot\cdot $, a pair of characteristic roots of Eq (4.1) cross the imaginary axis to the right. At each $ \tau_{n}^{1} $, $ n = 0, 1, 2, \cdot\cdot\cdot $, a pair of characteristic roots of Eq (4.1) cross the imaginary axis to the left. Then the transversality condition required by the Hopf bifurcation theorem is satisfied.

Note that $ \omega_{k} = \sqrt{u_{k}} $ and $ u_{1} < u_{2} $. We have $ \omega_{1} < \omega_{2} $, which indicates that

$ \tau^{1}_{n}-\tau^{1}_{n-1} = \frac{2\pi}{\omega_{1}} \gt \frac{2\pi}{\omega_{2}} = \tau^{2}_{n}-\tau^{2}_{n-1}, n = 0,1,2.... $

Then there exists a positive integer $ k $ such that

$ \tau_{0}^{2} \lt \tau_{0}^{1} \lt \tau_{1}^{2} \lt \tau_{1}^{1}\cdot\cdot\cdot \lt \tau_{k}^{2} \lt \tau_{k+1}^{2} \lt \tau_{k}^{1}. $

We thus obtain the following result.

Theorem 4.1. For system (1.4), we have

(i) If $ \Delta\leq0 $ and $ p_{3}\geq0 $ holds, $ E_{2}^{*} $ is asymptotically stable for all $ \tau > 0 $.

(ii) If one of the conditions ($ H_{1} $) and ($ H_{2} $) holds, $ E_{2}^{*} $ is asymptotically stable when $ \tau\in[0, \tau^{0}_{0}) $ and unstable when $ \tau > \tau^{0}_{0} $, that is system (1.4) undergoes a Hopf bifurcation at $ E_{2}^{*} $ when $ \tau = \tau^{0}_{0} $.

(iii) If the condition ($ H_{3} $) holds, system (1.4) undergoes a Hopf bifurcation at $ E_{2}^{*} $ along two sequences of $ \tau $ values $ \tau_{n}^{1, 2} $, $ n = 0, 1, 2 $.... Furthermore there exists a positive integer $ k $ such that $ E_{2}^{*} $ is stable when

$ \tau\in[0,\tau_{0}^{2})\cup(\tau_{0}^{1},\tau_{1}^{2})\cup\cdot\cdot\cdot\cdot\cdot\cup(\tau_{k-1}^{1},\tau_{k}^{2}); $

$ E_{2}^{*} $ is unstable when

$ \tau\in(\tau_{0}^{2},\tau_{0}^{1})\cup(\tau_{1}^{2},\tau_{1}^{1})\cup\cdot\cdot\cdot\cdot\cdot\cup(\tau_{k-1}^{2},\tau_{k-1}^{1})\cup(\tau_{k}^{2},+\infty). $

5. Direction and stability of the Hopf bifurcation

In the previous section, we know that system (1.4) undergoes a Hopf bifurcation at $ E_{2}^{*} $ along some sequences of $ \tau $ values. Let $ \tau^{*} $ be one of the Hopf bifurcation points. As pointed out in Hassard et al. ^[31], it is interesting to determine the direction, stability and period of these periodic solutions.

Let $ \mu = \tau-\tau^{*} $, and then $ \mu $ is new bifurcation parameter of the system. Define

$\begin{equation*} \label{9} X(t) = (x-x_{2}^{*},y-y_{2}^{*},z-z_{2}^{*})^{\mathrm{T}}, X_{t}(\theta) = X(t+\theta), \theta\in[-\tau,0]. \end{equation*}$ $

System (1.4) may be written as:

$\begin{equation} X^{'}(t) = L_{\mu}X_{t}+f(X_{t}(.),\mu),\ \ \ L_{\mu}\phi = F_{1}\phi(0)+F_{2}\phi(-\tau), \end{equation}$ $

(5.1)

where

$f(\phi,\mu) = \left[ {

$\begin{array}{*{20}{cccc}} -\beta\phi_{1}(0)\phi_{2}(0)\\ \beta\phi_{1}(0)\phi_{2}(0)-p\phi_{2}(0)\phi_{3}(0) \\ c\phi_{2}(-\tau)\phi_{3}(-\tau)-\dfrac{r}{m}{\phi_{3}}^{2}(0) \\ \end{array}$ } \right],$

$F_{1} = \left[

$\begin{array}{cccc} -\beta y_{2}^{*}-d_{1} & -\beta x_{2}^{*} & 0 \\ \beta y_{2}^{*} & -d_{2}-pz_{2}^{*} & -py_{2}^{*} \\ 0 & 0 & r-d_{3}-\dfrac{2rz_{2}^{*}}{m} \\ \end{array}$ \right],$

and

$F_{2} = \left[ {

$\begin{array}{*{20}{cccc}} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & cz_{2}^{*} & cy_{2}^{*} \\ \end{array}$ } \right].$

By Riesz representation theorem, there exists a matrix $ \eta(\theta, \mu):[-\tau, 0]\rightarrow \mathbb{ R}^{3} $ whose components are bounded variation functions such that

$ L_{\mu}\phi = \int^{0}_{-\tau}d\eta(\theta,\mu)\phi(\theta), $

where

$ d\eta(\theta,\mu) = F_{1}\delta(\theta)d\theta+F_{2}\delta(\theta+\tau)d\theta $

and $ \delta(\theta) $ is the Dirac delta function.

For $ \phi\in C([-\tau, 0], R^{3}) $, we define

$\begin{equation*} A(\mu)\phi(\theta) = \left\{ \begin{aligned} & \frac{d\phi(\theta)}{d\theta}, &\theta &\in[-\tau,0), \\ & \int^{0}_{-\tau}d\eta(\xi,\mu)\phi(\xi), &\theta& = 0, \end{aligned} \right. \end{equation*}$ $

and

$\begin{equation*} R(\mu)\phi(\theta) = \left\{ \begin{aligned} & 0, &\theta&\in[-\tau,0), \\ & f(\phi,\mu), &\theta& = 0. \end{aligned} \right. \end{equation*}$ $

Then system (1.4) is equivalent to the following operator equation

$\begin{equation} X_{t}^{'}(\theta) = A(\mu)X_{t}(\theta)+R(\mu)X_{t}(\theta). \end{equation}$ $

(5.2)

For $ \varphi\in C([0, \tau], R^{3}) $, we define the adjoint operator of $ A(0) $ as $ A^{\ast}(0) $, where

$\begin{equation*} A^{\ast}(0)\varphi(s) = \left\{ \begin{aligned} & -\frac{d\varphi(s)}{ds}, &s&\in(0,\tau], \\ & \int^{0}_{-\tau}d\eta^{\mathrm{T}}(\xi,0)\phi(-\xi), &s& = 0. \end{aligned} \right. \end{equation*}$ $

For the convenience of research, we simply write $ A $ for $ A(0) $, $ A^{\ast} $ for $ A^{\ast}(0) $, $ R $ for $ R(0) $, $ \eta(\theta) $ for $ \eta(\theta, 0) $, for $ \varphi\in C([0, \tau], R^{3}) $ and $ \phi\in C([-\tau, 0], R^{3}) $.

Define a bilinear form as

$\begin{equation*} \label{12} \langle\varphi,\phi\rangle = \bar{\varphi}^{\mathrm{T}}(0)\phi(0)-\int^{0}_{\theta = -\tau}\int^{\theta}_{\xi = 0}\bar{\varphi}^{\top}(\xi-\theta)d\eta(\theta)\phi(\xi)d\xi. \end{equation*}$ $

Let $ q(\theta) $ and $ q^{\ast}(s) $ to be the eigenvectors of matric $ A $ and $ A^{\ast} $ corresponding to eigenvalue $ i\omega_{0} $ and $ -i\omega_{0} $, respectively. Then

$\begin{equation*} \label{13} Aq(\theta) = i\omega_{0}q(\theta),\ \ \ \ \ \ \ A^{\ast}q^{\ast}(s) = -i\omega_{0}q^{\ast}(s). \end{equation*}$ $

We can choose appropriate $ q(\theta) $ and $ q^{\ast}(s) $ such that $ < q(\theta), q^{\ast}(s) > = 1 $,

where

$ q(\theta) = (1,q_{2},q_{3})^{\mathrm{T}}e^{i\omega_{0}\theta},\ \ \ q^{\ast}(s) = D(1,q^{\ast}_{2},q^{\ast}_{3})^{\mathrm{T}}e^{i\omega_{0}s}, $

and

$\begin{align*} q_{2}& = \frac{\beta y_{2}^{*}+d_{1}+i\omega_{0}}{-\beta x_{2}^{*}},\; \; \; \; \; \; \; \; \; q_{3} = \frac{\beta^{2} y_{2}^{*}x_{2}^{*}+(\beta y_{2}^{*}+d_{1}+i\omega_{0})(d_{2}+pz_{2}^{*}+i\omega_{0})}{py_{2}^{*}\beta x_{2}^{*}},\\ q_{2}^{\ast}& = \frac{\beta y_{2}^{*}+d_{1}-i\omega_{0}}{\beta y_{2}^{*}},\; \; \; \; \; \; \; \; \; q_{3}^{\ast} = \frac{[\beta^{2} y_{2}^{*} x_{2}^{*}+(\beta y_{2}^{*}+d_{1}-i\omega_{0})(d_{2}+pz_{2}^{*}-i\omega_{0})]e^{i\omega_{0}\tau^{*}}}{c\beta y_{2}^{*}z_{2}^{*}}. \end{align*}$ $

Note that

$\begin{align*} \lt q,q^{\ast} \gt & = \bar{D}\bar{q}^{\mathrm{T}}(0)q^{\ast}(0)-\int^{0}_{\theta = -\tau^{*}}\int^{\theta}_{s = 0}\bar{D}\bar{q}^{\mathrm{T}}(s-\theta)d\eta(\theta)q^{\ast}(s)ds\\ & = \bar{D}(1+\bar{q}_{2}q_{2}^{\ast}+\bar{q}_{3}q_{3}^{\ast}-\int^{0}_{\theta = -\tau^{*}}\int^{\theta}_{s = 0}(1,\bar{q}_{2},\bar{q}_{3})e^{-(s-\theta)i\omega_{0}}(1,q_{2}^{\ast},q_{3}^{\ast})^{\mathrm{T}}e^{i\omega_{0}s}dsd\eta(\theta))\\ & = \bar{D}(1+\bar{q}_{2}q_{2}^{\ast}+\bar{q}_{3}q_{3}^{\ast}-\int_{\theta = -\tau^{*}}^{0}(c\bar{q}_{2}z^{*}+cq_{3}^{\ast}y^{*})\theta e^{i\omega_{0}\theta}d\eta(\theta))\\ & = \bar{D}(1+\bar{q}_{2}q_{2}^{\ast}+\bar{q}_{3}q_{3}^{\ast}+(c\bar{q}_{2}z^{*}+cq_{3}^{\ast}y^{*})q_{3}^{\ast}\tau^{*}e^{i\omega_{0}\tau^{*}}). \end{align*}$ $

Then we can choose $ \bar{D} $ as

$ \bar{D} = [1+\bar{q}_{2}q_{2}^{\ast}+\bar{q}_{3}q_{3}^{\ast}+(c\bar{q}_{2}z_{2}^{*} +cq_{3}^{\ast}y_{2}^{*})q_{3}^{\ast}\tau^{*}e^{i\omega_{0}\tau^{*}}]^{-1} $

such that $ < q(\theta), q^{\ast}(s) > = 1 $.

According to the notations in Hassard et al. ^[31], we need to compute the center manifold $ C_{0} $ at $ \mu = 0 $. Let $ u_{t} $ be the solution of Eq (5.1) when $ \mu = 0 $ and define

$\begin{equation} Z(t) = \lt q^{\ast},u_{t} \gt ,\ \ W(t,0) = u_{t}(0)-2Re\{Z(t)q(\theta)\}. \end{equation}$ $

(5.3)

On the center manifold $ C_{0} $, we have $ W(t, 0) = W(Z(t), \bar{Z}(t), \theta) $, where

$ W(Z,\bar{Z},\theta) = W_{20}(\theta)\frac{Z^{2}}{2}+W_{11}(\theta)Z\bar{Z}+W_{02}(\theta)\frac{\bar{Z}^{2}}{2}+\cdots, $

$ Z $ and $ \bar{Z} $ are local coordinates for center manifold $ C_{0} $ in the direction of $ q^{\ast} $ and $ \bar{q}^{\ast} $. Note that if $ X_{t} $ is real, $ W $ is real. We only consider the real solutions. For the solution of the Eq (5.1) $ X_{t}\in C_{0} $, we have

$\begin{equation*} \dot{Z}(t) = i\omega Z+\bar{q}^{\ast}(\theta)f(0,W(Z,\bar{Z},\theta))+2\text{Re}\{Zq(\theta)\}\overset{\triangle}{ = }i\omega Z+\bar{q}^{\ast}(0)f_{0}(Z,\bar{Z}). \end{equation*}$ $

Now we rewrite this equation as $ \dot{Z}(t) = i\omega Z+g(Z, \bar{Z}) $, where

$\begin{equation} g(Z,\bar{Z}) = \bar{q}^{\ast}(0)f_{0}(Z,\bar{Z}) = g_{20}\frac{Z^{2}}{2}+g_{11}Z\bar{Z}+g_{02}\frac{\bar{Z}^{2}}{2}+g_{21}\frac{Z^{2}\bar{Z}^{}}{2}+\cdots. \end{equation}$ $

(5.4)

Note that $ u_{t}(\theta) = W(t, \theta)+Zq(\theta)+\bar{Z}\bar{q}(\theta) $. We have

$\begin{align*} &g(Z,\bar{Z}) = \bar{q}^{\ast}(0)f_{0}(Z,\bar{Z})\\ & = \bar{D}\big\{-\beta u_{1t}(0)u_{2t}(0)+q_{2}^{\ast}\big(\beta u_{1t}(0)u_{2t}(0)-pu_{2t}(0)u_{3t}(0)\big)+q^{\ast}_{3}\big(cu_{2t}(-\tau)u_{3t}(-\tau)-\frac{r}{m}u_{3t}^{2}(0)\big)\big\}\\ & = \bar{D}\Big\{(q^{\ast}_{2}-\beta)(Z+\bar{Z}+W_{20}^{(1)}(0)\frac{Z^{2}}{2}+W_{11}^{(1)}(0)Z\bar{Z}+W_{02}^{(1)}(0)\frac{\bar{Z}^{2}}{2})(q_{2}Z+\bar{q_{2}}\bar{Z}+W_{20}^{(3)}(0)\frac{Z^{2}}{2}\\ &+W_{11}^{(3)}(0)Z\bar{Z}+W_{02}^{(3)}(0)\frac{\bar{Z}^{2}}{2})-pq^{\ast}_{2}(q_{2}Z+\bar{q_{2}}\bar{Z}+W_{20}^{(2)}(0)\frac{Z^{2}}{2}+W_{11}^{(2)}(0)Z\bar{Z}+W_{02}^{(2)}(0)\frac{\bar{Z}^{2}}{2})(q_{3}Z\\ &+\bar{q_{3}}\bar{Z}+W_{20}^{(3)}(0)\frac{Z^{2}}{2}+W_{11}^{(3)}(0)Z\bar{Z}+W_{02}^{(3)}(0)\frac{\bar{Z}^{2}}{2})+q^{\ast}_{3}c(q_{2}Ze^{-i\omega_{0}\tau}+\bar{q_{2}}\bar{Z}e^{-i\omega\tau}+W_{20}^{(2)}(-\tau)\frac{Z^{2}}{2}\\ &+W_{11}^{(2)}(-\tau)Z\bar{Z}+W_{02}^{(2)}(-\tau)\frac{\bar{Z}^{2}}{2})(q_{4}Ze^{-i\omega_{0}\tau}+\bar{q_{4}}\bar{Z}e^{-i\omega\tau}+W_{20}^{(3)}(-\tau)\frac{Z^{2}}{2}+W_{11}^{(3)}(-\tau)Z\bar{Z}\\ &+W_{02}^{(3)}(-\tau)\frac{\bar{Z}^{2}}{2})-\frac{r}{m}q^{\ast}_{3}\big(q_{3}Z+\bar{q_{3}}\bar{Z}+W_{20}^{(3)}(0)\frac{Z^{2}}{2} +W_{11}^{(3)}(0)Z\bar{Z}+W_{02}^{(3)}(0)\frac{\bar{Z}^{2}}{2}\big)^{2}\Big\}. \end{align*}$ $

Comparing the coefficients with Eq (5.4), we have

$\begin{equation} \begin{aligned} g_{20}& = 2\bar{D}\big\{(q^{\ast}_{2}-\beta)q_{2}+q^{\ast}_{3}q_{2}q_{3}ce^{-2i\omega\tau}-q^{\ast}_{2}q_{2}q_{3}p-\frac{r}{m}q^{\ast}_{3}q_{3}^{2}\big\},\\ g_{11}& = 2\bar{D}\big\{(q^{\ast}_{2}-\beta)\text{Re}\{q_{2}\}+q^{\ast}_{3}c\text{Re}\{\bar{q_{2}}q_{3}\} e^{-2i\omega\tau}-q^{\ast}_{2}p\text{Re}\{\bar{q_{2}}q_{3}\}-\frac{r}{m}q^{\ast}_{3}q_{3}\bar{q_{3}} \big\},\\ g_{02}& = 2\bar{D}\big\{(q^{\ast}_{2}-\beta)\bar{q_{2}}+q^{\ast}_{3}\bar{q_{2}}\bar{q_{3}}ce^{-2i\omega\tau} -q^{\ast}_{2}\bar{q_{2}}\bar{q_{3}}p-\frac{r}{m}q^{\ast}_{3}\bar{q_{3}}^{2}\big\},\\ g_{21}& = 2\bar{D}\Big\{(q^{\ast}_{2}-\beta)\big(W^{(2)}_{11}(0)+\frac{W^{(2)}_{20}(0)}{2}+\frac{\bar{q_{2}}W^{(1)}_{20}(0)}{2}+q_{2}W^{(1)}_{11}(0)\big)\\ &+q^{\ast}_{3}c\big(q_{2}e^{-i\omega\tau}W^{(3)}_{11}(-\tau)+q_{3}e^{-i\omega\tau}W^{(2)}_{11}(-\tau) +\frac{\bar{q_{2}}W^{(3)}_{20}(-\tau)}{2}e^{-i\omega\tau} +\frac{\bar{q_{3}}W^{(2)}_{20}(-\tau)}{2}e^{-i\omega\tau}\big)\\&-q_{2}^{*}p\big(q_{2}W^{(3)}_{11}(0) +\frac{\bar{q_{2}}W^{(3)}_{20}(0)}{2}+\frac{\bar{q_{3}}W^{(2)}_{20}(0)}{2} +q_{3}W_{11}^{(2)}(0)\big)\\&-\frac{r}{m}q^{\ast}_{3}\big(2q_{3}W^{(3)}_{11}(0)+\bar{q_{3}}W^{(3)}_{20}(0)\big)\Big\}. \end{aligned} \end{equation}$ $

(5.5)

It remains to compute $ W_{11}(0) $ and $ W_{20}(0) $ in $ g_{21} $. From Eqs (5.2) and (5.3), we have

$\begin{equation} \dot{W} = \dot{x_{t}}-\dot{Z}q-\overline{\dot{Z}q} = \left\{ \begin{aligned} & AW-2\text{Re}\{\bar{q}^{\ast}(0)f_{0}q(\theta)\}, &\theta&\in[-1,0), \\ & AW-2\text{Re}\{\bar{q}^{\ast}(0)f_{0}q(0)\}+f_{0}, &\theta& = 0. \end{aligned} \right. \end{equation}$ $

(5.6)

We rewrite the Eq (5.6) as

$\begin{equation} \dot{W}\overset{\triangle}{ = }AW+H(Z,\bar{Z},\theta) \end{equation}$ $

(5.7)

where

$\begin{equation*} \label{16} H(Z,\bar{Z},\theta) = H_{20}(\theta)\frac{Z^{2}}{2}+H_{11}(\theta)Z\bar{Z}+H_{02}(\theta)\frac{\bar{Z}^{2}}{2}+\cdots. \end{equation*}$ $

Thus

$\begin{equation} (A-2i\omega)W_{20}(\theta) = -H_{20}(\theta),\ \ \ \ AW_{11}(\theta) = -H_{11}(\theta). \end{equation}$ $

(5.8)

From Eq (5.7), we know that for $ \theta\in[-1, 0) $

$\begin{equation*} \label{18} H(Z,\bar{Z},\theta) = -\bar{q^{\ast}}(0)f_{0}q(\theta)-q^{\ast}(0)\bar{f_{0}}\bar{q}(\theta) = -gq(\theta)-\bar{g}\bar{q}(\theta). \end{equation*}$ $

Comparing the coefficients with Eq (5.8), we obtain

$\begin{equation} H_{20}(\theta) = -g_{20}q(\theta)-\bar{g}_{20}\bar{q}(\theta), \end{equation}$ $

(5.9)

$\begin{equation} H_{11}(\theta) = -g_{11}q(\theta)-\bar{g}_{11}\bar{q}(\theta). \end{equation}$ $

(5.10)

From Eqs (5.8) and (5.9) and the definition of $ A $, we have

$ \dot{W}_{20}(\theta) = 2i\omega W_{20}(\theta)+g_{20}q(\theta)+\bar{g}_{02}\bar{q}(\theta), $

where $ q(\theta) = (1, q_{2}, q_{3})^{\mathrm{T}}e^{i\theta\omega\tau} $. Hence

$\begin{equation*} W_{20}(\theta) = \frac{ig_{20}}{\omega}e^{i\omega\theta}q(0)+\frac{i\bar{g}_{02}}{3\omega}e^{-i\omega\theta}\bar{q}(0)+E_{20}e^{2i\omega\theta} \end{equation*}$ $

where

$ E_{20} = 2(2i\omega-F_{1}-F_{2}e^{-2i\omega\tau})^{-1}\left[

$\begin{array}{c} -\beta q_{2}\\ \beta q_{2}-pq_{2}q_{3}\\ cq_{2}q_{3}e^{-2i\omega\tau}-\dfrac{r}{m}q_{3}^{2} \\ \end{array}$ \right]. $

Similarly, from Eqs (5.8) and (5.10) we obtain

$\begin{equation*} \label{19} W_{11}(\theta) = \frac{ig_{11}}{\omega}e^{i\omega\theta}q(0)+\frac{i\bar{g}_{11}}{\omega}e^{-i\omega\theta}\bar{q}(0)+E_{11}, \end{equation*}$ $

where

$ E_{11} = 2(-F_{1}-F_{2})^{-1}\left[

$\begin{array}{c} -\beta Re\{q_{2}\}\\ \beta Re\{q_{2}\}-pRe\{q_{2}\bar{q}_{3}\}\\ cRe\{q_{2}\bar{q_{3}}\}-\dfrac{r}{m}\{q_{3}^{2}\} \\ \end{array}$ \right]. $

So far, we have calculated $ g_{20} $, $ g_{11} $, $ g_{02} $, $ g_{21} $ in Eq (5.5) and then we can obtain

$\begin{equation} \begin{aligned} c_{1}(0)& = \frac{i}{2\omega}\Big(g_{11}g_{20}-2\mid g_{11}\mid^{2}-\frac{\mid g_{02}\mid^{2}}{3}\Big)+\frac{g_{21}}{2}, \\ \nu_{2} & = -\frac{\text{Re}(c_{1}(0))}{\text{Re}(\lambda^{'}(\tau^{*}))}, \\ \beta_{2} & = 2\text{Re}(c_{1}(0)),\\ T_{2}& = \frac{-\big(\text{Im}\{c_{1}(0)\}+\nu_{2}\text{Im}\{\lambda^{'}(\tau^{*})\}\big)}{\omega}. \end{aligned} \end{equation}$ $

(5.11)

It is well known that $ \nu_{2} $ and $ \beta_{2} $ will determine the direction and stability of the Hopf bifurcation, and $ T_{2} $ determines the period of the bifurcated periodic solutions, respectively. In particular, the Hopf bifurcation is supercritical (subcritical) if $ \nu_{2} > 0 (\nu_{2} < 0) $, and the bifurcated periodic solutions exist for $ \tau > \tau_{0} (\tau < \tau_{0}) $. The bifurcated periodic solutions are stable (unstable) if $ \beta_{2} < 0 $ ($ \beta_{2} > 0 $) and the period will become longer (shorter) if $ T_{2} > 0 $ ($ T_{2} < 0 $).

6. Numerical simulations

In this section, we carry out some numerical simulations to display some qualitative behaviours of system (1.4). Table $ 1 $ lists the values or ranges of parameters for system (1.4) referring to ^[21].

We first study the effect of logistic growth on the dynamics of system (1.4). According to the ranges of parameters in Table 1, we take the following parameters

$\begin{equation} s = 10,\beta = 0.02,d_{1} = 0.2,d_{2} = 0.5,d_{3} = 0.2,p = 0.05,c = 0.4,m = 20. \end{equation}$ $

(6.1)

Table 1. Meanings and units of parameters.

Parameters	Descriptions	Ranges	Units
$ s $	Source rate of uninfected cell	$ 1 $–$ 10 $	$ \text{cells} $ $ \text{mL}^{-1} \text{day}^{-1} $
$ \beta $	Virus-to-cell infection rate	$ 0.00025 $–$ 0.5 $	$ \text{mL} $ $ \text{virion}^{-1} \text{day}^{-1} $
$ p $	Predation rate of infection cell by CTLs	$ 1 $–$ 4.048\times10^{-4} $	$ \text{mL} $ $ \text{cell}^{-1} \text{day}^{-1} $
$ r $	Breath rate of CTLs	$ 0.0051 $–$ 3.912 $	$ \text{mL} $ $ \text{cell}^{-1} \text{day}^{-1} $
$ m $	Capacity of immune cells	$ 6.25 $–$ 235999.9 $	$ \text{mL} $ $ \text{cell}^{-1} $
$ c $	Development rate of CTLs	$ 0.0051 $–$ 3.912 $	$ \text{mL} $ $ \text{cell}^{-1} \text{day}^{-1} $
$ d_{1} $	Death rate of uninfected cell	$ 0.007 $–$ 0.1 $	$ \text{day}^{-1} $
$ d_{2} $	Death rate of infected cell	$ 0.2 $–$ 0.3 $	$ \text{day}^{-1} $
$ d_{3} $	Death rate of immune cell	$ 0.001 $–$ 0.3 $	$ \text{day}^{-1} $

| Show Table

DownLoad: CSV

Figure 1 presents the bifurcation diagrams of the solutions for system (1.4) with respect to $ \tau $ and different logistic growth rate $ r $. One can see that the positive equilibrium $ E_{1}^{*} $ of system (1.4) is local asymptotically stable when $ \tau < \tau^{*} = 0.8 $ and the bifurcated periodic solutions occurs through Hopf bifurcations when $ \tau > \tau^{*} $. Similarly, the positive equilibrium $ E_{2}^{*} = (49.42, 0.12, 9.77) $ of system (1.4) is local asymptotically stable when $ \tau = 8 < \tau^{*} = 8.655 $ and the bifurcated periodic solutions occur through Hopf bifurcations when $ \tau = 9 > \tau^{*} $. At the same time, one can note that, except $ r = 0 $, the values of Hopf bifurcation points increase with the increase of $ r $. On the other hand, the amplitude of the bifurcated periodic solution increase with the increase of time delay $ \tau $ and decrease with the increase of $ r $.

Figure 1. Bifurcation diagrams of system (1.4) with respect to $ \tau $ and different $ r $. Parameter values are given by (6.1).

DownLoad: Full-Size Img PowerPoint

Figure 2 presents phase diagrams of the solutions for system (1.4) with $ r = 0.3 $ and different values of $ \tau $. One can see that the positive equilibrium $ E_{2}^{*} = (49.42, 0.12, 9.77) $ is local asymptotically stable when $ \tau = 8 < \tau^{*} = 8.655 $ and the bifurcated periodic solutions occur through Hopf bifurcations when $ \tau = 9 > \tau^{*} $. Furthermore, using the given parameter values (6.1), we can obtain $ c_{1}(0) = -19.651+26.769 i $ by some calculations and then we know that $ \mu_{2} > 0, $ $ \beta_{2} < 0, $ $ T_{2} < 0 $ by (5.11). According to the conclusion in ^[31], we know that the Hopf bifurcation at $ \tau^{*} = 8.655 $ is supercritical, and the bifurcating periodic solution is stable.

In order to investigate the stability switch of equilibrium, we take the following parameters

$\begin{equation} s = 10,\beta = 0.2,d_{1} = 0.2,d_{2} = 0.3,d_{3} = 0.2,p = 0.05,c = 0.5,r = 1.6,m = 50. \end{equation}$ $

(6.2)

Figure 2. Phase diagrams of the solutions for system (1.4) with $ r = 0.3 $ and different values of $ \tau $. The other parameters values are given by (6.1).

DownLoad: Full-Size Img PowerPoint

We obtain there exists positive equilibrium $ E_{2}^{*} = (18.8775, 1.6487, 69.5102) $ and the characteristic equation of system (1.4) have two pure virtual roots $ \omega_{1} = 0.4464 $ and $ \omega_{2} = 1.4748 $. Then by (4.7) we get

$ \tau^{1}_{0} = 0.8365, \tau^{1}_{1} = 8.1911, \tau^{1}_{2} = 12.4513; $

and

$ \tau^{2}_{0} = 3.9309, \tau^{2}_{1} = 14.9113. $

Noting that $ \tau^{2}_{1} > \tau^{1}_{2} $, we know that $ E_{2}^{*} $ is asymptotically stable when $ \tau\in[0, \tau_{0}^{1})\cup(\tau_{0}^{2}, \tau_{1}^{1}) $ and is unstable when $ \tau\in(\tau_{0}^{1}, \tau_{0}^{2})\cup(\tau_{1}^{1}, +\infty) $. Figure 3 presents time series diagrams of the solutions for system (1.4) with different values of $ \tau $. One can see that there exists stability switch of equilibrium $ E_{2}^{*} $ as $ \tau $ increases and some periodic solutions are bifurcated by Hopf bifurcation.

Figure 3. Time series diagrams of the solutions for system (1.4) showing stability switch with increase of $ \tau $. Parameter values are given by (6.2).

DownLoad: Full-Size Img PowerPoint

7. Conclusions and discussions

In this paper, a viral infection model with self-proliferation of cytotoxic T lymphocytes (CTLs) and activation time delay of immune cells is proposed. We mainly focus on two topics. First, we study the global dynamics of system (1.4) through constructing appropriate Lyapunov functions. Then we examine the impact of the activation time delay of immune cells on the existence of periodic solutions. The global dynamical properties of system (1.4) can be summarized in the following Table 2. Numerical simulations verify the theoretical analyses. In particular, we find that for different values of the delay in CTL response, the system can stabilize at the positive equilibrium when the delay is small, or stabilize at a stable periodic oscillation when the delay is large. The amplitudes of bifurcated periodic solutions increase with the increase of activation time delay of immune cells and decrease with the increase of $ r $ (Figure 1). It is also found that there exists stability switch phenomenon under some conditions (Figure 3). The results indicate that the self-proliferation intensity and activation time delay of immune cells can significantly affect the kinetics of viral infection.

Table 2. Global properties of system (1.4) with $ \tau > 0 $.

Case	Conditions	Equilibria and its stability
$ r=0 $	$ R_{0} < 1 $	$ E_{0}^{yz} $ is GAS
	$ 1 < R_{0} < 1+\frac{\beta d_{3}}{cd_{1}} $	$ E_{0}^{yz} $ is US, $ E_{0}^{z} $ is GAS
	$ R_{0} > 1+\frac{\beta d_{3}}{cd_{1}} $	$ E_{0}^{yz} $ and $ E_{0}^{z} $ are US, $ E^{*} $ (Hopf)
$ r < d_{3} $	$ R_{0} < 1 $	$ E_{0}^{yz} $ is GAS
	$ 1 < R_{0} < 1+\frac{\beta (d_{3}-r)}{cd_{1}} $	$ E_{0}^{yz} $ is US, $ E_{0}^{z} $ is GAS
	$ R_{0} > 1+\frac{\beta (d_{3}-r)}{cd_{1}} $	$ E_{0}^{yz} $ and $ E_{0}^{z} $ are US, $ E^{*} $ (Hopf)
$ r=d_{3} $	$ R_{0} < 1 $	$ E_{0}^{yz} $ is GAS
	$ R_{0} > 1 $	$ E_{0}^{yz} $ and $ E_{0}^{z} $ are US, $ E^{*} $ (Hopf)
$ r > d_{3} $	$ R_{0} < 1 $	$ E_{0}^{yz} $ is US, $ E_{0}^{y} $ is GAS
	$ 1 < R_{0} < 1+\frac{pm(r-d_{3})}{rd_{2}} $	$ E_{0}^{yz} $ and $ E_{0}^{z} $ are US, $ E_{0}^{y} $ is GAS
	$ R_{0} > 1+\frac{pm(r-d_{3})}{rd_{2}} $	$ E_{0}^{yz} $, $ E_{0}^{z} $ and $ E_{0}^{y} $ are US, $ E^{*} $ (Hopf)
Note: GAS means globally asymptotically stable; US means unstable.

| Show Table

DownLoad: CSV

Some aspects of the viral infection problem remain to be studied in the future. For instance, we plan to extend our analysis to global Hopf bifurcation analysis. Some other factors that influence the dynamics of viral infection including the heterogeneity of space and movement of cells may be investigated.

Acknowledgments

The authors would like to thank the reviewers and the editor for their careful reading, helpful comments and suggestions that greatly improved the paper. This work is supported by the National Natural Science Foundation of China (Grant Nos. 11701472, 11771448, 11871403).

Conflict of interest

The authors declare that they have no conflict of interest.

References

[1]	ACOG Practice Bulletin No. 127 Management of Preterm Labor. June 2012.
[2]	U.S. Food and Drug Administration Drug Safety Communication. FDA recommends against prolonged use of magnesium sulfate to stop pre-term labor due to bone changes in exposed babies. 05/30/2013.
[3]	ACOG Committee Opinion, No. 573. (2013) Magnesium Sulfate Use in Obstetrics. Obstet Gynecol 122: 727-8.
[4]	Vintzileos, AM. (2009) Evidence-Based Compared With Reality-Based Medicine in Obstetrics. Obstet Gynecol 113: 1335-40.
[5]	Evidence-Based Medicine Working Group. (1992) Evidence-Based Medicine: A New Approach to Teaching the Practice of Medicine. JAMA 268:2420-5. doi: 10.1001/jama.1992.03490170092032
[6]	Sackett, DL. (1997) Evidence-Based Medicine. Seminar Perinatol 21: 3-5. doi: 10.1016/S0146-0005(97)80013-4
[7]	Steer CM, Petrie, RH. (1977) A comparison of magnesium sulfate and alcohol for the prevention of premature labor. Am J Obstet Gynecol 29: 1-4.
[8]	Elliott JP. (1983) Magnesium sulfate as a tocolytic agent. Am J Obstet Gynecol 147: 277-84. doi: 10.1016/0002-9378(83)91111-0
[9]	Elliott JP, Lewis DF, Morrison JC. (2009) In Defense of Magnesium Sulfate. Obstet Gynecol 113: 1341-7.
[10]	Ma L. (1992) Magnesium sulfate in the prevention of preterm labor (translation). Chung Hua I Hsveh Chin Taipla 72: 158-61.
[11]	Fox MD, Allbert JR, McCaul JF, Martin RW, McLaughlin BN, Morrison JC. (1993) Neonatal morbidity between 34 and 37 weeks gestation. J Perinatol XIII: 349-53.
[12]	Cox SM, Sherman ML, Leveno KJ. (1990) Randomized investigation of magnesium sulfate for prevention of preterm birth. Am J Obstet Gynecol 163 (3): 568-72.
[13]	Elliott JP. (1990) Letter to the editor- Sub therapeutic doses of magnesium sulfate do not inhibit preterm labor. Am J Obstet Gynecol 163: 568.
[14]	Mercer BM and Merlino AA. (2009) Magnesium Sulfate Preterm Labor and Preterm Birth. Am J Obstet Gynecol 114: 650-66.
[15]	Hung J, Tsai M, Yang B. (2005) Maternal Osteoporosis After Prolonged Magnesium Sulfate Tocolysis therapy: A Case Report. Arch Phys Med Rehabil 86: 146-9. doi: 10.1016/j.apmr.2004.02.016
[16]	Conde-Agudelo A, Romero R. (2009) Antenatal magnesium sulfate for the prevention of cerebral palsy in preterm infants less than 34 weeks gestation. Am J Obstet Gynecol 200: 595-609. doi: 10.1016/j.ajog.2009.04.005
[17]	ACOG Committee Opinion No. 455. Magnesium sulfate before anticipated preterm birth for neuroprotection. March 2010.
[18]	Scott JR. (2009) Magnesium Sulfate for neuroprotection. What do we do now? Obstet Gynecol 114: 500-01.
[19]	Reeves SA, Gibbs RS, and Clark SL. (2011) Magnesium for fetal neuroprotection. Am J Obstet Gynecol 204: 202.e1-4.
[20]	Mittendorf R, Covert R, Bowman J, et al. (1997) Is tocolytic magnesium sulfate associated with increased total pediatric mortality. Lancet 350: 1517-8. doi: 10.1016/S0140-6736(97)24047-X
[21]	Crowther CA, Hiller JE, Doyle LW, Haslam RR. (2003) Effect of magnesium sulfate given for neuroprotection before preterm birth: a randomized controlled trial. JAMA 290: 2669-76. doi: 10.1001/jama.290.20.2669
[22]	Farkouh LJ,Thorp JA, Jones PG, Clark RH, Knox GE. (2001) Antenatal magnesium exposure and neonatal demise. Am J Obstet Gynecol 185: 869-72. doi: 10.1067/mob.2001.117362
[23]	Mittendorf R, Dambrosia J, Pryde PG, et al. (2002) Associationbetween the use of antenatal magnesium sulfate in preterm labor and adverse health outcomes in infants. Am J Obstet Gynecol 186: 1111-8. doi: 10.1067/mob.2002.123544
[24]	Elliott JP, Morrison JC. The Evidence Regarding Maintenance Tocolysis. Obstetrics and Gynecology International Hindawi Publishing Corporation. Vol 2013, Article ID 708023, 11 pages. Available from: http://dx.doi.org/10.1155/2013/708023.
[25]	Yokoyama K, Takahashi N, Yada Y. (2010) Prolonged maternal magnesium administration and bone metabolism in neonates. Early Hum Dev 86: 187-91. doi: 10.1016/j.earlhumdev.2010.02.007
[26]	Wedig KE, Kogan J, Schorry EK et al. (2006) Skeletal demineralization and fractures caused by fetal magnesium toxicity. J Perinatol 26: 371-4. doi: 10.1038/sj.jp.7211508
[27]	Malaeb SN, Rassi A, Haddad MC. (2001) Bone mineralization in newborns whose mothers received magnesium sulfate for tocolysis of premature labor. Pediatr Radiol 34: 384-6.
[28]	Kaplan W, Haymond MW, Mckay S, Karaviti LP. (2006) Osteopenic effects of magnesium sulfate in multiple pregnancies. J Pediatric Endocrinology and Metabolism 19:1225-30.
[29]	Nassar AH, Sakhel K, Maarouf H et al. (2006) Adverse maternal and neonatal outcome of prolonged course of magnesium sulfate tocolysis. Acta Obstet Gynecol Scan 19: 1225-30.
[30]	Schanier RJ, Smith LG, Burns PA. (1997) Effects of long-term maternal intravenous magnesium sulfate therapy on neonatal calcium metabolism and bone mineral content. Gynecol Obstet Invest 43: 236-41. doi: 10.1159/000291864
[31]	Doyle LE, Anderson PJ, Haslam R, et al. (2014) School-age outcomes of very preterm infants after antenatal treatment with magnesium sulfate vs. placebo. JAMA 312(11): 1105-13.
[32]	ACOG Practice Bulletin No. 33. (2002) Diagnosis and management of preeclampsia and eclampsia. American College of Obstetricians and Gynecologists. Obstet Gynecol 99: 159-67.
[33]	Keys S, Elliott JP. (1997) The therapeutic use of vaginal pessary as adjunctive treatment in patients with multiple gestations presenting with preterm labor and low fetal station. J Reprent Med 42: 751-5.

This article has been cited by:

1.	Abdul Qadeer Khan, Ayesha Yaqoob, Ateq Alsaadi, Neimark-Sacker bifurcation, chaos, and local stability of a discrete Hepatitis C virus model, 2024, 9, 2473-6988, 31985, 10.3934/math.20241537
2.	Yuequn Gao, Ning Li, Fractional order PD control of the Hopf bifurcation of HBV viral systems with multiple time delays, 2023, 83, 11100168, 1, 10.1016/j.aej.2023.10.032

Reader Comments

Your name:*

Email:*
© 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Public Health

2.7 5.3

Metrics

Article views(7161) PDF downloads(1275) Cited by(19)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

AIMS Public Health

Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL)

Related Papers:

Abstract

1. Introduction

2. Preliminary results

3. Global stability analysis of boundary equilibria

4. Stability of positive equilibrium and Hopf bifurcation

5. Direction and stability of the Hopf bifurcation

6. Numerical simulations

7. Conclusions and discussions

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Other Articles By Authors

Catalog

AIMS Public Health

Risks and Benefits of Magnesium Sulfate Tocolysis in Preterm Labor (PTL)

Related Papers:

Abstract

1. Introduction

2. Preliminary results

3. Global stability analysis of boundary equilibria

4. Stability of positive equilibrium and Hopf bifurcation

5. Direction and stability of the Hopf bifurcation

6. Numerical simulations

7. Conclusions and discussions

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog