Impulsive Antibody Therapy for Hopf Bifurcation Control in Human SARS-CoV-2 Dynamics

1. Introduction

The outbreak of the novel coronavirus disease (COVID-19), caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has posed unprecedented challenges to global health systems [1]. The virus primarily enters human cells through the angiotensin-converting enzyme 2 (ACE2) receptor, initiating a cascade of interactions between susceptible cells, infected cells, viral particles, and the host immune response [2]. Understanding the within-host dynamics of SARS-CoV-2 infection is crucial for designing effective therapeutic strategies [3,4]. Mathematical modeling provides a powerful framework to capture these complex biological processes and to predict possible outcomes under different therapeutic interventions [5,6].

Recent studies have emphasized the importance of antibody mediated immune responses in controlling viral replication. Antibodies neutralize viral particles and reduce infection rates, while cytotoxic immune cells eliminate infected cells [7]. However, the immune response alone may be insufficient to eradicate the virus, particularly in cases of high viral load or impaired immunity [8]. In such scenarios, external therapeutic interventions, such as antibody-based drug therapy, play a critical role. Impulsive drug therapy, where antibodies are administered at discrete time intervals, reflects realistic medical practices and introduces periodic forcing into the system dynamics.

Mathematical modeling of the within-host dynamics of SARS-CoV-2 provides critical insights into viral replication, immune response, and treatment efficacy [9]. A widely used approach involves systems of ordinary differential equations (ODEs) that describe the interactions between uninfected epithelial cells, infected cells, and free virions. Nath et al. [10] analyzed a model incorporating innate and adaptive immune responses, deriving the basic reproduction number ${\mathcal{R}}_0$ and demonstrating that antiviral therapy can significantly reduce peak viral load and infection duration. Song et al. extended this framework by simulating the effects of different antiviral treatment timings, showing that early intervention leads to faster viral clearance and reduced tissue damage [11]. Another study by Kumar et al. [12] incorporated both immune responses and cytokine dynamics, revealing that delayed immune activation correlates with severe disease progression and prolonged viral shedding. These models not only help quantify key parameters such as viral clearance rate and immune efficacy but also support the design of optimized treatment strategies and vaccine protocols.

Hopf bifurcation plays a pivotal role in epidemiological modeling by elucidating the emergence of oscillatory dynamics in infectious diseases [13]. Classical epidemic models for diseases such as measles, influenza, and COVID-19 typically predict either disease extinction or stabilization at an endemic equilibrium [14,15]. A Hopf bifurcation arises when a pair of complex conjugate eigenvalues of the Jacobian matrix cross the imaginary axis as system parameters vary, thereby destabilizing the equilibrium and giving rise to periodic orbits [16]. In epidemiological contexts, parameters such as the infection transmission rate or viral replication rate can induce this bifurcation, leading to recurrent oscillatory outbreaks [17,18]. The occurrence of Hopf bifurcation in disease models provides a theoretical explanation for recurrent epidemics and identifies critical thresholds at which interventions may stabilize system dynamics. Consequently, Hopf bifurcation serves as a powerful analytical tool for interpreting and controlling epidemic behavior, enhancing our ability to predict and mitigate disease progression [19,20]. Despite its significance, this phenomenon has not been fully investigated in the context of within-host dynamics of SARS-CoV-2 infection. In this study, we address this gap by proposing and analyzing a mathematical model that captures such dynamics.

New antiviral drugs are broadening treatment options for COVID-19. The oral antiviral ensitrelvir, already approved in Japan, has shown promise in reducing infection risk after exposure [21,22,23]. It complements existing therapies such as Paxlovid (nirmatrelvir/ritonavir), Veklury (remdesivir), and Lagevrio (molnupiravir), which are most effective when started soon after symptoms appear [24,25]. Monoclonal antibody (mAb) treatments including casirivimab/imdevimab, sotrovimab, and Evusheld are widely used under Emergency Use Authorization [26,27], but their effectiveness has declined with the emergence of new variants like Omicron, leading to revoked or suspended authorizations for several agents [28,29]. This shifting landscape highlights the need to continually adapt treatment strategies, with guidance from regulatory bodies such as the FDA and IDSA ensuring therapies remain aligned with evolving variant susceptibility [30,31,32]. In our study, antibody therapy is considered within an impulsive treatment framework, reflecting its potential as a periodic intervention strategy.

The present study develops and analyzes a mathematical model of SARS-CoV-2 infection within the human host, explicitly incorporating ACE2 receptor dynamics, antibody response, and impulsive antibody drug therapy. Particular attention is given to the existence of equilibria, with emphasis on the conditions under which Hopf bifurcation occurs. The model is further extended to an impulsive framework that integrates drug administration dynamics. For this impulsive system, the stability of the disease-free periodic orbit is examined using Floquet theory. The analytical results demonstrate the potential of impulsive therapy to stabilize the disease-free state or suppress oscillatory viral persistence, thereby providing valuable insights for clinical strategies aimed at combating COVID-19.

In the following sections, the remaining parts of this article are presented. The Section 2 contains the formulation of the mathematical model under immune and antibody responses. The qualitative properties such as basic reproduction number, equilibrium points and their stability of the model have been provided in Section 3. Theoretical analysis namely the existence and stability of disease-free periodic orbit of the impulsive model are studied in Section 4. Dynamics of the impulsive system is analysed in Section 6. Numerical simulations of the analytical findings are presented in Section 6 outcomes. The discussion and concluding remarks are given in Section 7.

2. Mathematical Model Derivation

SARS-CoV-2 predominantly targets epithelial cells located in the upper respiratory tract and the initial segments of the bronchi. The severity of the disease increases when the virus spreads to the lower regions of the lungs. The infection process unfolds in two main stages. Initially, the virus binds to angiotensin-converting enzyme 2 (ACE2) receptors found on the surface of host cells [33,34]. Following the attachment, the viral genetic material penetrates the cell, initiating infection. These infected cells then begin producing new viral particles through a cycle known as the eclipse and burst phases. Zhou et al. [35] were successful in identifying the target tissue of SARS-CoV-2 infection in the human body as the epithelial tissue, with the S protein of the virus having a strong binding affinity to ACE2 receptors. They were of the view that the successful eradication of the virus from the body depends on optimal strategies being adopted to treat the disease.

To investigate the within-host dynamics of SARS-CoV-2 infection, a deterministic model is proposed that comprising six populations listed below:

(i)

Uninfected susceptible epithelial cells, $S\left(t\right)$, located in the respiratory tract (lungs, nasal, tracheal, and bronchial tissues),

(ii)

Infected virus-producing epithelial cells, $I\left(t\right)$,

(iii)

Free virus particles, $V\left(t\right)$,

(iv)

ACE2 receptors on epithelial cells, $R\left(t\right)$,

(v)

Cytotoxic T lymphocytes (CTLs) targeting infected cells, $C\left(t\right)$, and

(vi)

Antibodies produced by B cells in response to infection, $A\left(t\right)$.

• The model incorporates the following biological mechanisms:

(A1) Epithelial cells are recruited at a rate ${\lambda }_s$, while their natural death rate is given by ${\mu }_1$. A susceptible cell can be infected after interacting with viral particles and $ACE2$ receptors, which leads to the generation of infected cells $I\left(t\right)$ at a rate proportional to $\beta$.

(A2) Infected cells die from cytopathic effects, immune destruction, or a combination thereof at a rate proportional to ${\mu }_2$. Additionally, infected cells can be eliminated by cytotoxic T lymphocytes (CTLs) in a rate proportional to $pIC$.

(A3) Virus particles are generated by infected cells at a rate denoted by m. They are removed through natural clearance processes at rate ${\mu }_3$, and further neutralized at rate $qVA$. The ACE2 receptors are synthesized at rate ${\lambda }_r$, consumed during viral infection at rate $r\beta SRV$, and undergo natural turnover at rate ${\mu }_4$.

(A4) CTLs proliferate in response to infection at rate ${\displaystyle \frac{\alpha IC}{\theta +V}}$, are suppressed by viral interference at rate $\gamma pIC$, and decay at rate ${\mu }_5$. Antibodies are produced in proportion to CTL levels at rate $\eta C$ and decay at rate ${\mu }_6$.

The complete system of differential equations for the model is given by:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& {\lambda }_s-\beta SRV-{\mu }_{1}S,\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =& \beta SRV-{\mu }_{2}I-pIC,\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =& mI-{\mu }_{3}V-qVA,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =& {\lambda }_r-r\beta SRV-{\mu }_{4}R,\hfill \\ \hfill {\displaystyle \frac{dC}{dt}}& =& {\displaystyle \frac{\alpha IC}{\theta +V}}-\gamma pIC-{\mu }_{5}C,\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =& \eta C-{\mu }_{6}A,\hfill \end{array}$$

(1)

subjected to the initial cell concentrations:

$$\begin{array}{c}\hfill {\mathrm{E}}_S\left(0\right)\ge 0,{\mathrm{E}}_I\left(0\right)\ge 0,\mathrm{V}\left(0\right)\ge 0,\mathrm{R}\left(0\right)\ge 0,\mathrm{C}\left(0\right)\ge 0.\end{array}$$

(2)

3. Dynamics Study Using the Model

In this section, we conduct an analytical investigation of the dynamical behavior of the system (1). Specifically, we determine the equilibrium points and examine their local stability properties, including possible stability transitions. Furthermore, we derive the basic reproduction number ($R_0$) to characterize the stability of the disease-free equilibrium.

3.1. Equilibria

There are three equilibria of the system. They are given below:

(i) The disease-free equilibrium is $E_0(\overline{S},0,0,\overline{R},0,0)$ where, $\overline{R}={\displaystyle \frac{{\lambda }_r}{{\mu }_4}},\overline{S}={\displaystyle \frac{{\lambda }_s}{{\mu }_1}}.$

(ii) The immune response-free diseased equilibrium, $E_1(S_1,I_1,V_1,R_1,0,0)$, where

$$\begin{array}{cc}\hfill S_1& ={\displaystyle \frac{{\lambda }_s-X}{{\mu }_1}},R_1={\displaystyle \frac{{\lambda }_r-rX}{{\mu }_4}},I_1={\displaystyle \frac{X}{{\mu }_2}},V_1={\displaystyle \frac{mX}{{\mu }_2{\mu }_3}},\hfill \end{array}$$

with $X={\displaystyle \frac{(r{\lambda }_s+{\lambda }_r)\pm \sqrt{{(r{\lambda }_s+{\lambda }_r)}^2-4r({\lambda }_s{\lambda }_r-K)}}{2r}}$ and $K={\displaystyle \frac{{\mu }_1{\mu }_4{\mu }_2{\mu }_3}{\beta m}}.$

(iii) The diseased equilibrium is $E^{*}(S^{*},I^{*},V^{*},R^{*},C^{*},A^{*})$ with

$$\begin{array}{cc}\hfill A^{*}& ={\displaystyle \frac{\eta }{{\mu }_6}}{C}^{*},V^{*}={\displaystyle \frac{m{I}^{*}}{{\mu }_3+q{A}^{*}}},I^{*}={\displaystyle \frac{{\mu }_5(\theta +V^{*})}{\alpha -\gamma p(\theta +V^{*})}},\hfill \\ \hfill S^{*}& ={\displaystyle \frac{{\lambda }_s}{\beta R^{*}{V}^{*}+{\mu }_1}},R^{*}={\displaystyle \frac{{\lambda }_r}{r\beta S^{*}{V}^{*}+{\mu }_4}},\hfill \end{array}$$

where $C^{*}$ satisfies the cubic equation,

$$H\left(C\right)=a_3{C}^3+a_2{C}^2+a_{1}C+a_0=0,$$

(3)

with

$$\begin{array}{cc}\hfill a_3& ={(\gamma p\theta q\eta -\alpha q\eta )}^2,\hfill \\ \hfill a_2& =2(\gamma p\theta q\eta -\alpha q\eta )(\gamma p\theta {\mu }_6{\mu }_3+{\mu }_{5}m{\mu }_6-\alpha {\mu }_6{\mu }_3)-4\gamma pm{\mu }_6·{\mu }_5\theta q\eta ,\hfill \\ \hfill a_1& ={(\gamma p\theta {\mu }_6{\mu }_3+{\mu }_{5}m{\mu }_6-\alpha {\mu }_6{\mu }_3)}^2-8\gamma pm{\mu }_6{\mu }_5\theta {\mu }_6{\mu }_3,\hfill \\ \hfill a_0& =-4\gamma pm{\mu }_6{\mu }_5\theta {\mu }_6{\mu }_3.\hfill \end{array}$$

It can be noted that equation (3) is a cubic equation and its constant coefficient ($a_0<0$) and the coefficient of highest power ($a_3$) are of opposite sign implying the existence of a positive root. We only need $I^{*}>0$ i.e. $\alpha >\gamma p(\theta +V^{*})$ for the existence of an endemic equilibrium is always feasible for the system (1).

3.2. The Basic Reproduction Number

The basic reproduction number, ${\mathcal{R}}_0$, is determined using the next-generation approach [36]. Formally, ${\mathcal{R}}_0$ is given by the spectral radius (i.e., the dominant eigenvalue) of the matrix that represents the rate of new infections and the transitions out of infectious compartments. This matrix is evaluated at the disease-free equilibrium, denoted by $E_0=(\overline{S},0,0,\overline{R},0,0).$

We denote the infected subsystem as $x={(I,V)}^T$. Then the dynamics can be decomposed using the following subsystems

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& =\underset{F_I}{\underbrace{\beta SRV}}-\underset{V_I}{\underbrace{({\mu }_{2}I+pIC)}},\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =\underset{F_V}{\underbrace{0}}-\underset{V_V}{\underbrace{({\mu }_{3}V-mI+qVA)}}.\hfill \end{array}$$

The Jacobian matrices (i.e. the next generation matrices) are obtained at the DFE:

$$F=\left[\begin{array}{cc}{\displaystyle \frac{\partial F_I}{\partial I}}& {\displaystyle \frac{\partial F_I}{\partial V}}\\ {\displaystyle \frac{\partial V_I}{\partial I}}& {\displaystyle \frac{\partial V_I}{\partial V}}\end{array}\right]=\left[\begin{array}{cc}0& \beta \overline{S}\overline{R}\\ 0& 0\end{array}\right],$$

and

$$V=\left[\begin{array}{cc}{\displaystyle \frac{\partial F_V}{\partial I}}& {\displaystyle \frac{\partial F_V}{\partial V}}\\ {\displaystyle \frac{\partial V_V}{\partial I}}& {\displaystyle \frac{\partial V_V}{\partial V}}\end{array}\right]=\left[\begin{array}{cc}{\mu }_2& 0\\ -m& {\mu }_3\end{array}\right].$$

Now we derive

$$G=F{V}^{-1}={\displaystyle \frac{1}{{\mu }_2{\mu }_3}}\left[\begin{array}{cc}\beta m\overline{S}\overline{R}& \beta {\mu }_2\overline{S}\overline{R}\\ 0& 0\end{array}\right].$$

Here G is the next generation matrix. Since G is upper triangular, its eigenvalues are the diagonal entries. Therefore, the basic reproduction number is

$${\mathcal{R}}_0={\displaystyle \frac{\beta m\overline{S}\overline{R}}{{\mu }_2{\mu }_3}}={\displaystyle \frac{\beta m}{{\mu }_2{\mu }_3}}{\displaystyle \frac{{\lambda }_s}{{\mu }_1}}{\displaystyle \frac{{\lambda }_r}{{\mu }_4}}.$$

(4)

Remark 1.

This expression shows that ${\mathcal{R}}_0$ increases with the infection rate β, viral production rate m, and equilibrium levels of susceptible cells and receptors, while it decreases with the clearance rates of infected cells (${\mu }_2$) and free virus (${\mu }_3$).

3.3. Stability Analysis

Here, we provide the stability of the equilibrium points. We find the Jacobian matrix at the equilibrium and check the nature of the eigenvalues.

3.3.1. Stability Analysis of the Disease-Free Equilibrium

To analyze the local stability of the disease-free equilibrium (DFE), we evaluate the Jacobian matrix at

$$E=(\overline{S},\overline{I},\overline{V},\overline{R},\overline{C},\overline{A})=\left(\frac{{\lambda }_s}{{\mu }_1},0,0,\frac{{\lambda }_r}{{\mu }_4},0,0\right).$$

The Jacobian at the DFE is

$$J\left(E_0\right)=\left(\begin{array}{cccccc}-{\mu }_1& 0& -\beta \overline{S}\overline{R}& 0& 0& 0\\ 0& -{\mu }_2& \beta \overline{S}\overline{R}& 0& 0& 0\\ 0& m& -{\mu }_3& 0& 0& 0\\ 0& 0& -r\beta \overline{S}\overline{R}& -{\mu }_4& 0& 0\\ 0& 0& 0& 0& -{\mu }_5& 0\\ 0& 0& 0& 0& \eta & -{\mu }_6\end{array}\right).$$

It is evident that the four eigenvalues $-{\mu }_1,-{\mu }_4,-{\mu }_5,-{\mu }_6$, are strictly negative under biologically reasonable assumptions (${\mu }_i>0$). Therefore, the stability of the DFE depends solely on the remain eigenvalues which satisfy

$${\lambda }^2+({\mu }_2+{\mu }_3)\lambda +\left({\mu }_2{\mu }_3-m\beta \overline{S}\overline{R}\right)=0.$$

Since ${\mu }_2+{\mu }_3>0$, the sign of the constant term determines stability. Using (4), we obtain the following the results.

(i)

If ${\mathcal{R}}_0<1$, then ${\mu }_2{\mu }_3-m\beta \overline{S}\overline{R}>0$, and all eigenvalues have negative real parts. The DFE is locally asymptotically stable.

(ii)

If ${\mathcal{R}}_0>1$, then ${\mu }_2{\mu }_3-m\beta \overline{S}\overline{R}<0$, and one eigenvalue becomes positive. The DFE is unstable.

Theorem 1.

The DFE is locally asymptotically stable when ${\mathcal{R}}_0<1$. If ${\mathcal{R}}_0>1$, then the DFE is unstable. At the threshold ${\mathcal{R}}_0=1$, the system undergoes a bifurcation (transcritical forward type).

3.3.2. Stability of $E^{*}$

We derive the Jacobian Matrix at the endemic steady state $E^{*}$ as

$$J\left(E^{*}\right)=\left(\begin{array}{cccccc}-\beta RV-{\mu }_1& 0& -\beta SR& -\beta SV& 0& 0\\ \\ \beta RV& -{\mu }_2-pC& \beta SR& \beta SV& -pI& 0\\ \\ 0& m& -{\mu }_3-qA& 0& 0& qV\\ \\ -r\beta RV& 0& -r\beta SR& -r\beta SV-{\mu }_4& 0& 0\\ \\ 0& {\displaystyle \frac{\alpha C}{\theta +V}}& -{\displaystyle \frac{\alpha IC}{{(\theta +V)}^2}}& 0& -{\mu }_5& -{\gamma }_{1}pI\\ \\ 0& 0& -{\gamma }_{2}qA& 0& \eta & -{\mu }_6-{\gamma }_{2}qA\\ \end{array}\right).$$

The characteristic equation in x is obtained as,

$$\begin{array}{c}\hfill h\left(\rho \right)={\rho }^6+J_1{\rho }^5+J_2{\rho }^4+J_3{\rho }^3+J_4{\rho }^2+J_5\rho +J_6=0,\end{array}$$

(5)

Let us denote the matrix J as

$$J\left(E^{*}\right)={\left[a_{ij}\right]}_{1\le i,j\le 6},$$

then

$$\begin{array}{cc}& a_{11}=-\beta RV-{\mu }_1,a_{12}=0,a_{13}=-\beta SR,a_{14}=-\beta SV,a_{15}=0,a_{16}=0,\hfill \\ & a_{21}=\beta RV,a_{22}=-{\mu }_2-pC,a_{23}=\beta SR,a_{24}=\beta SV,a_{25}=-pI,a_{26}=0,\hfill \\ & a_{31}=0,a_{32}=m,a_{33}=-{\mu }_3-qA,a_{34}=0,a_{35}=0,a_{36}=qV,\hfill \\ & a_{41}=-r\beta RV,a_{42}=0,a_{43}=-r\beta SR,a_{44}=-r\beta SV-{\mu }_4,a_{45}=0,a_{46}=0,\hfill \\ & a_{51}=0,a_{52}={\displaystyle \frac{\alpha C}{\theta +V}},a_{53}=-{\displaystyle \frac{\alpha IC}{{(\theta +V)}^2}},a_{54}=0,a_{55}=-{\mu }_5,a_{56}=-{\gamma }_{1}pI,\hfill \\ & a_{61}=0,a_{62}=0,a_{63}=-{\gamma }_{2}qA,a_{64}=0,a_{65}=\eta ,a_{66}=-{\mu }_6-{\gamma }_{2}qA.\hfill \end{array}$$

then the coefficients of (5) can be written as

$$\begin{array}{cc}\hfill J_1& =-\sum _{i=1}^6{a}_{ii},\hfill \\ \hfill J_2& =\sum _{1\le i<j\le 6}\left(a_{ii}{a}_{jj}-a_{ij}{a}_{ji}\right),\hfill \\ \hfill J_3& =-\sum _{1\le i<j<k\le 6}\left|\begin{array}{ccc}{a}_{ii}& a_{ij}& a_{ik}\\ a_{ji}& a_{jj}& a_{jk}\\ a_{ki}& a_{kj}& a_{kk}\end{array}\right|,\hfill \\ \hfill J_4& =\sum _{1\le i<j<k<\ell \le 6}\left|\begin{array}{cccc}{a}_{ii}& a_{ij}& a_{ik}& a_{i\ell }\\ a_{ji}& a_{jj}& a_{jk}& a_{j\ell }\\ a_{ki}& a_{kj}& a_{kk}& a_{k\ell }\\ a_{\ell i}& a_{\ell j}& a_{\ell k}& a_{\ell \ell }\end{array}\right|,\hfill \\ \hfill J_5& =-\sum _{1\le i_1<i_2<i_3<i_4<i_5\le 6}det\left(A_{\{i_1,i_2,i_3,i_4,i_5\}}\right),\hfill \\ \hfill J_6& =det\left(J\right)=det\left(\left[a_{ij}\right]\right).\hfill \end{array}$$

According to Routh-Hurwitz conditions, characteristic equation (5) has roots with negative real parts if the following conditions are satisfied:

$$\begin{array}{ccc}\hfill & i)& J_i>0,i=1,2,...,6\hfill \\ & ii)& J_1{J}_2-J_3>0,\hfill \\ & iii)& J_3(J_1{J}_2-J_3)-J_1(J_1{J}_4-J_5)>0,\hfill \\ & iv)& J_1{J}_2{J}_3{J}_4-J_3^2{J}_4-J_1^2{J}_4^2-J_1{J}_2^2{J}_5+J_2{J}_3{J}_5+2J_1{J}_4{J}_5\hfill \\ & & -J_5^2+J_1^2{J}_2{J}_6-J_1{J}_3{J}_6>0,\hfill \\ & v)& J_1{J}_2{J}_3{J}_4{J}_5-J_3^2{J}_4{J}_5-J_1^2{J}_4^2{J}_5-J_1{J}_2^2{J}_5^2+J_2{J}_3{J}_5^2+2J_1{J}_4{J}_5^2-J_5^3\hfill \\ & & -J_1{J}_2{J}_3^2{J}_6+J_3^3{J}_6+J_1^2{J}_3{J}_4{J}_6+2J_1^2{J}_2{J}_3{J}_6-3J_1{J}_3{J}_5{J}_6-J_1^3{J}_6^2>0.\hfill \end{array}$$

(6)

We now investigate whether the system undergoes a change in stability. A Hopf bifurcation occurs when a fixed point of a dynamical system becomes unstable as a parameter varies, giving rise to a periodic orbit (limit cycle). For a six-dimensional dynamical system, the occurrence of a Hopf bifurcation can be detected by analyzing the eigenvalues of the Jacobian matrix evaluated at the equilibrium point. In particular, the criterion is established through the application of the theorem described in [37].

Theorem 2.

The endemic equilibrium $E^{*}$ undergoes a Hopf bifurcation corresponding to the generic bifurcation parameter ζ at ${\zeta }^{*}$, if the following conditions are satisfied:

i)

$J_i\left({\zeta }^{*}\right)>0,\forall i=1,2,\dots ,6$

ii)

${[J_1{J}_2-J_3]}_{\zeta ={\zeta }^{*}}=0,$

iii)

${[J_1{J}_2{J}_3-J_1^2{J}_4-J_3^2-J_5]}_{\zeta ={\zeta }^{*}}=0,$

iv)

${\left[(J_1{J}_2-J_3)(J_3^2-J_1^2{J}_4+J_5)+(J_1{J}_2-J_3)(J_1{J}_6-J_1^2{J}_5)\right]}_{\zeta ={\zeta }^{*}}=0,$

v)

${\left[H_5\right]}_{\zeta ={\zeta }^{*}}={\left[J_1{J}_2{J}_3{J}_4-(J_1{J}_2-J_3)\left(J_4{J}_5{J}_6+J_3{J}_6(J_2{J}_5-J_1{J}_6)+J_1{J}_2{J}_6^2-2J_1{J}_2{J}_6^2{J}_1{J}_6^2-J_5(J_1{J}_4-J_2{J}_3)\right)\right]}_{\zeta ={\zeta }^{*}}=0,$

vi)

Transversality condition given below should be satisfied:

$${\displaystyle \frac{1}{M_6}}{\left[{\displaystyle \frac{d{H}_5}{d\zeta }}\right]}_{\zeta ={\zeta }^{*}}\ne 0$$

where,

$$M_6=2J_1(J_1{J}_4-J_5)-J_3(J_1{J}_2-J_3)·\left|\begin{array}{ccccc}-J_1& 1& 1& 0& 0\\ J_2-2\varphi & -J_1& 0& 1& 0\\ J_1\varphi -J_3& J_2-\varphi & 0& 1& 1\\ {\varphi }^2-J_2\varphi +J_4& 0& \varphi & 0& \varphi \\ 0& {\varphi }^2-J_2\varphi +J_4& 0& 0& \varphi \end{array}\right|,$$

with

$$\varphi ={\displaystyle \frac{J_3{J}_5+J_1^2{J}_6-J_1{J}_2{J}_5}{J_1^2{J}_4-J_1{J}_3-J_1{J}_2{J}_3+J_3^2}}.$$

4. The System with Impulsive Control

Now we impose the periodic application of drug dosing in the model (1). Thus the model becomes a impulsive differential equations result if drug effect is assumed to change instantaneously at dosing times, $t_j$ for different drug doses. In presence of antibody controlled therapy through perfect adherence, we consider the model system (1) to obtain the following model:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& ={\lambda }_s-\beta SRV-{\mu }_{1}S,\hfill \\ \hfill {\displaystyle \frac{dI}{dt}}& =\beta SRV-{\mu }_{2}I-pIC,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dV}{dt}}& =mI-{\mu }_{3}V-qVA,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& ={\lambda }_r-r\beta SRV-{\mu }_{4}R,\hfill \\ \hfill {\displaystyle \frac{dC}{dt}}& ={\displaystyle \frac{\alpha IC}{\theta +V}}-\gamma pIC-{\mu }_{5}C,\hfill \\ \hfill {\displaystyle \frac{dA}{dt}}& =\eta C-{\mu }_{6}A,t\ne t_n\hfill \\ \hfill A\left(t_n^{+}\right)& =\omega +A\left(t_n^{-}\right),t=t_n,\hfill \end{array}$$

(8)

with a periodic impulse sequence $t_n=t_0+nT$, $T>0$. $A\left(t_n^{-}\right)$ denotes the antibody responses immediately before the impulse drug dosing, $A\left(t_n^{+}\right)$ denotes the concentration after the impulse and $\omega$ is the dose that is taken at each impulse time $t_n$, $k\in \mathbb{N}$.

5. Dynamics of the Impulsive System

For a system with impulsive differential equations, there exists periodic orbits like equilibrium points in system (1). We focus our study on the disease-free periodic orbit and its stability conditions.

5.1. The Disease-Free Periodic Orbit

We can consider $I=0$, $V=0$, $C=0$ at the disease-free equilibrium. Thus the remaining population can be written as below:

$$\dot{S}={\lambda }_s-{\mu }_{1}S,\dot{R}={\lambda }_r-{\mu }_{4}R,\dot{A}=-{\mu }_{6}A.$$

Hence, at the disease-free equilibrium,

$$\tilde{S}={\displaystyle \frac{{\lambda }_s}{{\mu }_1}},\tilde{R}={\displaystyle \frac{{\lambda }_r}{{\mu }_4}},\tilde{C}\left(t\right)=0,\tilde{I}\left(t\right)=0,\tilde{V}\left(t\right)=0.$$

The Poincaré (impulse-to-impulse) map for $A\left(t\right)$ is

$$A_{n+1}^{+}=e^{-{\mu }_{6}T}{A}_n^{+}+\omega ,$$

(9)

which has a unique fixed point

$$A_{+}^{*}={\displaystyle \frac{\omega }{1-e^{-{\mu }_{6}T}}},A_{-}^{*}=A_{+}^{*}{e}^{-{\mu }_{6}T}.$$

(10)

Thus, between impulses $t\in [t_n,t_{n+1})$, we have

$$\tilde{A}\left(t\right)=A_{+}^{*}{e}^{-{\mu }_6(t-t_n)}.$$

Therefore the unique disease-free periodic orbit can be given by

$${\tilde{E}}_0\left(\tilde{S},\tilde{I},\tilde{V},\tilde{R},\tilde{C},\tilde{A}\left(t\right)\right)=\left({\displaystyle \frac{{\lambda }_s}{{\mu }_1}},0,0,{\displaystyle \frac{{\lambda }_r}{{\mu }_4}},0,A_{+}^{*}{e}^{-{\mu }_6(t-t_n)}\right),t\in [t_n,t_{n+1}).$$

5.1.1. Impulsive Reproduction Number ${\mathcal{R}}_0^{\mathrm{imp}}$

Near the disease-free orbit, the $(I,V)$ subsystem governs the invasion dynamics (the C-equation is exponentially stable and decoupled at first order):

$${\displaystyle \frac{d}{dt}}\left(\begin{array}{c}I\\ V\end{array}\right)=\left(\begin{array}{cc}-{\mu }_2& \beta \tilde{S}\tilde{R}\\ m& -{\mu }_3-q\tilde{A}\left(t\right)\end{array}\right)\left(\begin{array}{c}I\\ V\end{array}\right).$$

Let $\Phi \left(t\right)$ be the fundamental matrix with $\Phi \left(0\right)=I_2$. Over one impulse period T, the monodromy matrix is $M=\Phi \left(T\right)$.

Let $L\left(t\right)$ denote the time-periodic Jacobian of the $(I,V)$ subsystem along the disease-free periodic orbit:

$$L\left(t\right)=\left(\begin{array}{cc}-{\mu }_2& \beta \tilde{S}\tilde{R}\\ m& -{\mu }_3-q\tilde{A}\left(t\right)\end{array}\right),t\in [0,T],$$

and let $\Phi \left(t\right)$ be the fundamental matrix solution of $\dot{X}=L\left(t\right)X$ with $\Phi \left(0\right)=I_2$. The monodromy matrix over one impulse period is

$$M=\Phi \left(T\right).$$

The spectral radius of M is

$$\rho \left(M\right)=\underset{i=1,2}{max}\left|{\lambda }_i\left(M\right)\right|,$$

where ${\lambda }_i\left(M\right)$ are the (complex) eigenvalues of M. Equivalently, using trace and determinant for a $2\times 2$ matrix,

$${\lambda }_{1,2}\left(M\right)={\displaystyle \frac{tr\left(M\right)\pm \sqrt{tr{\left(M\right)}^2-4det\left(M\right)}}{2}},\rho \left(M\right)=max\left\{\left|{\lambda }_1\left(M\right)\right|,\left|{\lambda }_2\left(M\right)\right|\right\}.$$

5.1.2. Floquet Stability Criterion for ${\tilde{E}}_0$

The disease-free periodic orbit is locally asymptotically stable iff all Floquet multipliers lie inside the unit circle, i.e.,

$$\rho \left(M\right)<1,$$

where $\rho \left(M\right)$ denotes the spectral radius of the monodromy matrix. Equivalently, instability (invasion) occurs when the principal Floquet multiplier exceeds unity.

We introduced the impulse reproduction number, denoted ${\mathcal{R}}_0^{\mathrm{imp}}$, which is the leading Floquet multiplier that governs how infections evolve across one full cycle. In simpler terms, it measures the growth potential of the infection when the system is periodically driven by impulses. To make this precise, we look at the linear flow that describes how infections generate virus over the time interval $[0,T]$:

$$\dot{I}=-{\mu }_{2}I+\beta \tilde{S}\tilde{R}V,\dot{V}=mI-\left({\mu }_3+q\tilde{A}\left(t\right)\right)V.$$

The mapping that takes a small initial condition $\left(I\right(0),V(0\left)\right)$ to $\left(I\right(T),V(T\left)\right)$ is linear, with matrix M. Then,

$${\mathcal{R}}_0^{\mathrm{imp}}:=\rho \left(M\right).$$

where,

$$M=\left(\begin{array}{cc}-{\mu }_2& \beta \tilde{S}\tilde{R}\\ m& -{\mu }_3-q\tilde{A}\end{array}\right).$$

Hence, we have the following theorem.

Theorem 3.

The disease-free periodic orbit is locally asymptotically stable iff ${\mathcal{R}}_0^{\mathrm{imp}}<1,$ and unstable ${\mathcal{R}}_0^{\mathrm{imp}}>1$.

5.1.3. Antibody Periodic Level and Stability Analysis

To quantify the effect of impulsive antibody boosts on viral dynamics, we define the impulse-period average of the antibody concentration along the disease-free orbit. Specifically,

$${\tilde{A}}_m={\displaystyle \frac{1}{T}}{\int }_{t_n}^{t_{n+1}}\tilde{A}\left(t\right)dt={\displaystyle \frac{1}{T}}{\int }_0^T{A}_{+}^{*}{e}^{-{\mu }_{6}s}ds={\displaystyle \frac{A_{+}^{*}}{{\mu }_{6}T}}\left(1-e^{-{\mu }_{6}T}\right)={\displaystyle \frac{\omega }{{\mu }_{6}T}}.$$

This average provides a tractable measure of the effective antibody level sustained between impulses, and will serve as the basis for constructing a conservative constant-coefficient proxy for the time-periodic stability threshold.

Replacing the time-varying antibody concentration $\tilde{A}\left(t\right)$ with its average ${\tilde{A}}_m$ yields a simplified linear system with constant coefficients. The resulting Jacobian matrix is

$$\left(\begin{array}{cc}-{\mu }_2& \beta \tilde{S}\tilde{R}\\ m& -{\mu }_3-q{\tilde{A}}_m\end{array}\right).$$

For such constant systems, the basic reproduction number reduces to the ratio of infection and production terms to decay terms:

$${\widehat{\mathcal{R}}}_0={\displaystyle \frac{\beta \tilde{S}\tilde{R}m}{{\mu }_2\left({\mu }_3+q{\tilde{A}}_m\right)}}={\displaystyle \frac{\beta ({\lambda }_s/{\mu }_1)({\lambda }_r/{\mu }_4)m}{{\mu }_2\left({\mu }_3+q\omega /\left({\mu }_{6}T\right)\right)}}.$$

Now, the monodromy matrix governing stability over one impulse period is given by

$$M_T=Iexp\left({\int }_0^{T}A\left(t\right)dt\right),$$

where I denotes the $2\times 2$ identity matrix. If $\tilde{A}\left(t\right)$ is approximated by its constant-coefficient proxy $tilde{A}_m={\displaystyle \frac{\omega }{{\mu }_{6}T}},$ then

$${\widehat{A}}_{2\times 2}=\left(\begin{array}{cc}-{\mu }_2& \beta \tilde{S}\tilde{R}\\ m& -{\mu }_3-q{\tilde{A}}_m\end{array}\right),$$

and

$$M\approx e^{\widehat{A}T},\rho \left(M\right)\approx \underset{i}{max}\left|e^{{\lambda }_i\left(\widehat{A}\right)T}\right|=exp\left(T·\underset{i}{max}\Re \left({\lambda }_i\left(\widehat{A}\right)\right)\right).$$

This formulation gives the conservative proxy threshold value, ${\widehat{\mathcal{R}}}_0$ as

$${\widehat{\mathcal{R}}}_0={\displaystyle \frac{\beta \tilde{S}\tilde{R}m}{{\mu }_2\left({\mu }_3+q{\tilde{A}}_m\right)}}$$

(11)

This provides an explicit stability condition, while the exact stability is determined by the spectral radius $\rho \left(M\right)$ computed from the full time-periodic $A\left(t\right)$ i.e. using $\tilde{A}\left(t\right)$. The condition for stability is ${\widehat{\mathcal{R}}}_0<1$. This bound highlights how stronger or more frequent impulses (larger $\omega$, smaller T) and faster antibody decay (larger ${\mu }_6$) reduce viral persistence. The net effect enters through the averaged antibody level ${\tilde{A}}_m=\omega /\left({\mu }_{6}T\right)$.

Remark 2.(i) In general, ${\mathcal{R}}_0^{\mathrm{imp}}$ is computed by numerically integrating the time-periodic linear system over one period to obtain M and then evaluating $\rho \left(M\right)$.

(ii) For the full system, when ${\mathcal{R}}_0^{\mathrm{imp}}<1$ the disease-free periodic orbit is global asymptotic stabile; whereas ${\mathcal{R}}_0^{\mathrm{imp}}>1$ leads to uniform persistence and existence of a positive periodic orbit ${\tilde{E}}^{*}$, driven by the periodic antibody drug applications.

6. Numerical Simulation

In this section, we performed the numerical simulations to study the dynamics of the system by visualization in figures. We begin by presenting the sensitivity analysis of the model parameters. Subsequently, the system dynamics are examined in the absence of impulsive drug administration. Finally, numerical simulations are carried out for the impulsive system (7).

6.1. Sensitivity Analysis Using ${\mathcal{R}}_0$

The basic reproduction number ${\mathcal{R}}_0$ serves as a threshold parameter that determines whether SARS-CoV-2 infection can successfully establish within the human host. For the model (1), ${\mathcal{R}}_0$ is calculated as

$${\mathcal{R}}_0={\displaystyle \frac{\beta m{\lambda }_s{\lambda }_r}{{\mu }_1{\mu }_4{\mu }_2{\mu }_3}}.$$

To assess the relative importance of each parameter, sensitivity analysis is performed using ${\mathcal{R}}_0$. Sensitivity measures the proportional change in ${\mathcal{R}}_0$ resulting from a proportional change in a given parameter.

The results shown in Figure 1 revealed that all parameters contribute equally in magnitude, with sensitivities of $\pm 1$. Specifically, parameters $\beta$, m, ${\lambda }_s$, and ${\lambda }_r$ have positive sensitivities, indicating that increases in these values directly increase ${\mathcal{R}}_0$. Conversely, parameters ${\mu }_1$, ${\mu }_4$, ${\mu }_2$, and ${\mu }_3$ have negative sensitivities, meaning that increases in these clearance or death rates reduce ${\mathcal{R}}_0$. Thus, a 10% increase in viral production or infection rate results in a 10% increase in ${\mathcal{R}}_0$, while a 10% increase in viral clearance or infected cell death leads to a 10% reduction. Overall, the sensitivity analysis demonstrates that both viral and host parameters exert proportional influence on ${\mathcal{R}}_0$.

Table 1. Biological interpretation of the variables and parameters applied in the model [38,39,40].

Variables/Parameters	Explanation	Assigned Values
S	Count of susceptible (uninfected target) cells	—
I	Population of SARS-CoV-2 infected (virus-producing) cells	—
V	Quantity of free SARS-CoV-2 viral particles	—
R	Number of ACE2 receptors present on epithelial cells	—
C	Population of cytotoxic T lymphocytes (CTLs)	—
${\lambda }_r$	Rate of ACE2 receptor production	12
$\beta$	Infection transmission rate	${10}^{-7}-{10}^{-5}$
${\lambda }_s$	Production rate of epithelial cells	15
${\mu }_1$	Natural death rate of uninfected epithelial cells	0.1
${\mu }_2$	Death rate of infected epithelial cells	0.1
${\mu }_3$	Clearance rate of viral particles	1.67
${\mu }_4$	Clearance rate of ACE2 receptors	0.02
${\mu }_5$	Elimination rate of immune cells	0.02
p	Intracellular replication rate of virus	${10}^{-5}$
q	Rate of virus neutralization by antibodies	${10}^{-3}$
m	Number of virions generated per cell	$1-10$
$\alpha$	Antibody response rate from immune cells	0.82
${\mu }_6$	Clearance rate of antibodies	0.02

Table 1. Biological interpretation of the variables and parameters applied in the model [38,39,40].

Variables/Parameters	Explanation	Assigned Values
S	Count of susceptible (uninfected target) cells	—
I	Population of SARS-CoV-2 infected (virus-producing) cells	—
V	Quantity of free SARS-CoV-2 viral particles	—
R	Number of ACE2 receptors present on epithelial cells	—
C	Population of cytotoxic T lymphocytes (CTLs)	—
${\lambda }_r$	Rate of ACE2 receptor production	12
$\beta$	Infection transmission rate	${10}^{-7}-{10}^{-5}$
${\lambda }_s$	Production rate of epithelial cells	15
${\mu }_1$	Natural death rate of uninfected epithelial cells	0.1
${\mu }_2$	Death rate of infected epithelial cells	0.1
${\mu }_3$	Clearance rate of viral particles	1.67
${\mu }_4$	Clearance rate of ACE2 receptors	0.02
${\mu }_5$	Elimination rate of immune cells	0.02
p	Intracellular replication rate of virus	${10}^{-5}$
q	Rate of virus neutralization by antibodies	${10}^{-3}$
m	Number of virions generated per cell	$1-10$
$\alpha$	Antibody response rate from immune cells	0.82
${\mu }_6$	Clearance rate of antibodies	0.02

6.2. Dynamics Without Impulses

Figure 2 contains the forward bifurcation diagram. Disease-free equilibrium (DFE) is stable when ${\mathcal{R}}_0$ is less than unity and unstable otherwise. Consequently a forward bifurcation occurs at ${\mathcal{R}}_0=1$, marking the transition between disease extinction and persistence (Theorem 1).

Regions of stability of the disease-free equilibrium is shown in Figure 3 in different parameter planes. When the product of the infection rate and viral replication rates crosses a threshold value, then disease-free equilibrium becomes unstable via forward transcritical bifurcation. Results of this figures is important for disease management. If we can control the viral load or infection rate, the disease can be eradicated.

Numerical simulations are presented in Figure 4, where the influence of the infection rate on the system dynamics is examined. For relatively low values of the infection rate, the endemic equilibrium remains stable. However, when the infection rate increases to $\beta =2\times {10}^{-6}$, the system exhibits periodic oscillations. This behavior indicates a transition in stability, specifically a change from a stable focus to a periodic orbit, which arises through a Hopf bifurcation.

Furthermore, the occurrence of Hopf bifurcation is demonstrated by considering both the infection rate and the viral replication rate as bifurcation parameters (Figure 5 and Figure 6 respectively). As these parameters exceed their respective critical thresholds, the system undergoes a qualitative change in dynamics, leading to the emergence of periodic solutions characteristic of Hopf bifurcation (Theorem 2).

The regions of stability corresponding to the equilibrium points are illustrated in Figure 7. When the infection rate ($\beta$), viral replication rate (m), and the production rate of the ACE2 receptor $\left({\lambda }_r\right)$ are relatively low, the disease fails to persist. As these parameters increase, the system transitions to an endemic state, and for sufficiently high values, the equilibrium point $E_1$ becomes stable. Furthermore, when the infection rate is elevated and the immune response is diminished, stability of $E_1$ is again observed.

These results imply that minor changes in infection or replication rates can significantly impact the course of a disease. Chronic infections that maintain a balanced viral load and immune response are known as stable equilibria. On the other hand, as seen in some chronic viral infections, periodic oscillations could indicate recurrent flare-ups or cyclical patterns of viral activity. The discovery of Hopf bifurcation provides insight into how therapeutic interventions that lower infection or replication rates could stop oscillatory, potentially harmful dynamics by highlighting critical thresholds beyond which the host system loses stability.

6.3. Dynamics with Impulses

Figure 8 displays the behaviour of the impulsive system. The computed value of ${\widehat{\mathcal{R}}}_0$ for the impulse-incorporating system is $0.429272$, which is less than unity. The disease-free periodic orbit is asymptotically stable, according to this result. Consequently, the endemic system can be driven toward a disease-free state through the application of impulsive drug dosing (Theorem 3).

Moreover, different dosing strategies may be employed to stabilize the disease-free orbit, as illustrated in Figure 9. The syability to the disease-free orbit can be achieved either by adopting a longer impulsive interval with a reduced drug dosage or, alternatively, by implementing a shorter interval with an increased dosage.

These results emphasize how important treatment scheduling is to infection control. Even in situations where the system would otherwise maintain an endemic equilibrium, timely interventions may be able to eradicate infection due to the stability of the disease-free orbit under impulsive dosing. Proper drug administration i.e. frequency (T) and dosage ($\omega$) can minimize viral persistence, lessen treatment burden, and improve patient outcomes, all of which have practical implications for therapeutic design. Crucially, the findings highlight that timing and dose magnitude function as forward bifurcation parameters, dictating whether the host system successfully eliminates the infection or stays infected.

7. Discussion and Conclusions

Mathematical modeling serves as a rigorous and effective framework for analyzing the dynamics of infectious diseases and for designing strategies to control their spread. Hopf bifurcation analysis has emerged as a fundamental tool in mathematical biology for investigating the onset of oscillatory behavior in nonlinear systems. Within the context of viral infections, Hopf bifurcation reveals the conditions under which the system transitions from a stable equilibrium to sustained oscillations in viral load and immune response. Such oscillatory dynamics may correspond to recurrent infection patterns, fluctuating antibody levels, or periodic viral persistence, thereby offering valuable insight into the mechanisms underlying chronic or relapsing infections.

We have studied the model dynamics using analytical and numerical methods. We calculate the equilibrium points of the system without drug dosing and found three equilibrium points. We derive the basic reproduction number of the system without and with impulses to find the disease-free situation. We perform the sensitivity analysis using the basic reproduction number.

From a biological perspective, the basic reproduction number and the sensitivity analysis highlights the therapeutic significance of targeting viral replication and infection processes. For instance, entry inhibitors that reduce the infection rate ($\beta$) or polymerase/protease inhibitors that suppress viral replication (m) are effective strategies for lowering the basic reproduction number ${\mathcal{R}}_0$. Similarly, therapies that enhance viral clearance (${\mu }_3$) or promote infected cell death (${\mu }_2$) can prevent the establishment of infection. Therapeutic strategies should prioritize reducing viral replication while enhancing clearance mechanisms to stabilize the disease-free equilibrium.

In addition to continuous-time modeling, this research incorporates impulsive differential equations to capture the effects of discrete therapeutic interventions. Using Floquet theory, the existence and stability of the disease-free periodic orbit are established. Numerical simulations confirm that impulsive drug dosing can drive the system from an endemic state to a disease-free state. Importantly, the results show that stabilization can be achieved either through longer impulsive intervals with smaller doses or shorter intervals with larger doses, thereby offering flexibility in treatment design.

In conclusion, this study not only advances the theoretical understanding of within-host SARS-CoV-2 dynamics but also provides practical insights into the timing and efficacy of antibody-based interventions. By integrating impulsive antibody therapy into the model, the system exhibits hybrid dynamics that combine continuous viral-host interactions with discrete therapeutic effects. This hybrid structure enriches the bifurcation landscape and emphasizes the critical role of treatment scheduling in infection control. Collectively, the findings suggest that mathematical models can guide the optimization of therapeutic strategies, ensuring effective suppression of viral persistence and promoting long-term disease-free outcomes.