Global Dynamics of SARS-CoV-2 Infection with Antibody Response and the Impact of Impulsive Drug Therapy

Mathematical modeling is crucial to investigating tthe ongoing coronavirus disease 2019 (COVID-19) pandemic. The primary target area of the SARS-CoV-2 virus is epithelial cells in the human lower respiratory tract. During this viral infection, infected cells can activate innate and adaptive immune responses to viral infection. Immune response in COVID-19 infection can lead to longer recovery time and more severe secondary complications. We formulate a micro-level mathematical model by incorporating a saturation term for SARS-CoV-2-infected epithelial cell loss reliant on infected cell levels. Forward and backward bifurcation between disease-free and endemic equilibrium points have been analyzed. Global stability of both disease-free and endemic equilibrium is provided. We have seen that the disease-free equilibrium is globally stable for R0<1, and endemic equilibrium exists and is globally stable for R0>1. Impulsive application of drug dosing has been applied for the treatment of COVID-19 patients. Additionally, the dynamics of the impulsive system are discussed when a patient takes drug holidays. Numerical simulations support the analytical findings and the dynamical regimes in the systems.


Introduction
COVID-19 is considered to be transmitted mainly between people who are close in contact with one another within about six feet, as well as through respiratory droplets created when an infected person coughs or sneezes. These droplets can enter the mouths or noses of people close by or possibly be inhaled into the lungs [1]. Virus spread depends on the possibility of touching virus-infected surfaces or objects and then touching one's own mouth, nose, or possibly eyes [2,3].
In a SARS-CoV-2-infected human, the innate and adaptive immune responses work together to neutralize the threat of SARS-CoV-2 infection [4][5][6]. When the virus enters the human body, the innate immune response starts immediately. Proteins of the natural immune system in a healthy cell also respond against the invading pathogens within the first minutes or hours of infection [7]. This response is of great importance in preventing new infections through the activation of the adaptive immune system [8,9]. Cytokines, which are small soluble proteins, are an essential component of the immune system [10]. They are secreted from different cells in the human body. They can be categorized into one of four families: (i) the hematopoietic family, (ii) the immunoglobin superfamily, (iii) the tumor necrosis factor family, and (iv) interferons (IFNs) [11]. Cytokines balance the innate and adaptive immune responses. Among cytokines, IFNs play a vital role in the innate immune response during viral infection. Thus, we consider the effect of adaptive immune responses in our mathematical model.
There are many research articles that include population modeling for the transmission dynamics of COVID-19 [12][13][14][15][16][17][18] These studies mainly focused on susceptible exposed populations and asymptomatic infected populations for a particular region. Population density is a major factor for disease transmission for this type of modeling. Some of these articles include vaccinations and optimal control [19][20][21]; additionally, some articles focus on population awareness through media [22,23].
However, for the SARS-CoV-2 dynamics at the micro-level (i.e., the dynamics of the disease within the human host), few model-based studies are available. The dynamics of SARS-CoV-2 infection can give insight to controlling the virus in a human host [11,24,25]. Many infectious disease dynamics are explored extensively by researchers with the help of mathematical modeling with real data at the cellular level [26][27][28][29]. Tang et al. [30] proposed a four-population host cell infection model for MERS-CoV mediated by DPP4 receptors. The infection processes of SARS-CoV-2, SARS-CoV, and MERS-CoV are similar. Researchers are still working on inter-host modeling for SARS-CoV-2 and target cell limitations under immune responses [31,32]. Hernandez et al. [31] proposed a model to examine cellular level dynamics and T cell responses against viral replication. Wang et al. [32] evaluated the effect of several potential interventions for SARS-CoV-2. Their study reveals that combining antiviral drugs with interferons effectively reduces the viral plateau phase and shortens the recovery time. Chatterjee and Bashir [29] formulated a mathematical model to examine the consequences of adaptive immune response to viral mutation to control disease transmission. They also studied the effect of antiviral drug therapy and its impact on model dynamics. Chatterjee et al. [11] proposed a set of fractional differential equation models considering uninfected epithelial cells, infected epithelial cells, SARS-CoV-2 virus, and CTL response cells, accounting for the lytic and non-lytic effects of immune responses [11]. They also studied the impact of a commonly used antiviral drug in COVID-19 treatment in an optimal control-theoretic approach. Wang et al. [28] studied the effect of antiviral drugs against SARS-CoV-2 viral dynamics during COVID-19 infection. In [33], within-host dynamics of SARS-CoV-2 infection were studied with potential treatments. Authors have used repurposed drugs (remdesivir) that inhibit the transcription of SARS-CoV-2.
Despite numerous therapeutic strategies, to date, there is no specific effective treatment for SARS-CoV-2 infection. Recently, all over the world, clinicians have been working on an effective therapy for coronavirus disease 2019 . Clinical observation suggests that cytokine levels enhance the hyperinflammatory response secondary to SARS-CoV-2 infection. This is the leading cause of multi-organ damage for COVID-19 patients [1]. For these reasons, numerous clinical trials are currently undergoing to explore the effectiveness of drugs such as interleukin-1 blockers and interleukin-6 inhibitors in COVID-19 [33][34][35].
The most useful method to study drug dynamics is the use of impulsive differential equations [28]. Perfect or imperfect drug adherence and drug holidays can make the development of resistance easy. Recently, the effects of perfect drug dosing on antiretroviral therapy have been studied by impulsive differential equations. Chatterjee and Basir [24] formulated a dynamic model of epithelial cells during SARS-CoV-2 infection and CTL responses. They established a new mathematical model considering epithelial cells and the role of the ACE2 receptor using impulsive differential equations, which describe the within-host dynamics of SARS-CoV-2 infection with treatments. The dosing period and threshold values of dosage can be obtained using this method.
In the present research, we study the dynamics of COVID-19 in humans using immunostimulant drugs. We explor the dynamics using the antibody-response model of SARS-CoV-2 infection by examining the interaction between viral replication. We consider uninfected target epithelial cells, infected epithelial cells, SARS-CoV-2 virus, and antibody responses in the modeling process, with an aim to reduce the infected epithelial cells and viral load using immunostimulant drugs. The local and global dynamics of the system without drugs have been provided. Forward transcritical bifurcation is also analyzed. Finally, we implement impulsive differential equations to observe the impact of drugs. The dosing rate and interval, and how many drug holidays a patient takes have been studied.
The drug data we have used here are those of monoclonal antibodies (mAbs). These have the capability to detect and prevent the disease propagation [36,37]. A SARS-CoV-2 patient who is at high risk of transmitting to another individual with SARS-CoV-2 for longer than 4 weeks, and who is unable to mount an adequate immune response to SARS-CoV-2 vaccination, can take an initial dose of 600 mg of casirivimab and 600 mg of imdevimab, then repeat doses of 300 mg of casirivimab and 300 mg of imdevimab once every 4 weeks [38].
The article is organized as follows. The next section (Section 2) contains the formulation of a mathematical model of immune responses. The qualitative properties of the model are provided in Section 3. Theoretical analysis of the impulsive model is carried out in Section 4. The numerical simulation is included based on the analytical findings in Section 5. Discussion and concluding remarks are given in Section 6.

Derivation of the Mathematical Model
The mathematical model helps us to understand the basic dynamics of viral infection. In general cases, modeling consist of a antibody response model with some variants [39]. Here, we consider the simplest version including three populations, namely: The SARS-CoV-2 dynamics with immune cells are proposed in [39] as the following The first equation of (1) shows the dynamics of uninfected epithelial cells (E S (t)); the second equation shows the dynamics of the infected epithelial cells (E I (t)). The replication of the SARS-CoV-2 virus (V(t)) in the third equation of (1) is considering, as SARS infection promotes endothelins on several organs as a direct consequence of viral involvement [31].
The growth rate of epithelial cell is denoted as Π. The virus infects the uninfected cells with a rate β (mL (RNA copies) −1 day −1 ). After a cell becomes infected, it behaves as a virus-producing cell and produces viruses at a rate p (day −1 ), and are virus particle cleared at a rate µ 3 (day −1 ). The uninfected susceptible cells are cleared at a rate µ 1 (day −1 ) due to their natural apoptosis, and the infected cells are removed from the system at a rate µ 2 (day −1 ) as a result of cytopathic viral effects and immune response [31].
Cytokine is vital in inhibiting viral replication and modulating downstream effects of the immune response. Specific cytokines activate natural killer cells (NK), which act against virally infected cells. In the case of SARS-CoV-2 infection, it is observed that viruses often target the JaK/STAT pathway (i.e., a chain of interactions between proteins in a cell) to decrease the production of IFNs. This immune suppressing mechanism observed in SARS-CoV-2 can be represented in the functional form of a decrease in the cytokine production rate, assumed to be αE I V+θ . Cytokines activate the adaptive immune system, mainly T cells and B lymphocytes, to produce an antibody that acts against the virus. B cells mainly secrete IgM and IgG antibodies that are released from blood and lymph fluid and neutralize the SARS-CoV-2 viral particles.
Considering the antibody responses A(t), we extend the antibody-response model to include the depletion of viruses modeled via the term rAV. The extended model of (1) reads as follows: with the initial condition The graphical representation of the above model is shown in Figure 1. Here, r is the rate at which the antibody neutralizes the viral particles, α is the antibody simulation rate constant, and θ is the strength of antibody suppression, in which the antibodies are lost at a rate of µ 4 .
We now impose impulsive drug dosing in the above model and analyze it in Section 4. Impulsive differential equations result from drug effects; the metabolites are assumed to decay with time in an exponential manner during each cycle and are assumed to change instantaneously at dosing times t j for different drug doses, which can result in either implicit or explicit models. In the presence of antibody-controlled therapy with perfect adherence, we consider Model (2). Before analyzing the system, we first discuss the one-dimensional impulse system as follows: D(t − n ) denotes the drug concentration immediately before the impulse drug dosing, D(t + n ) denotes the concentration after the impulses, and ω is the dose that is taken at each impulse time t n , n ∈ N. Here, ζ is the rate at which antibodies are produced due to the use of a drug.

Dynamics of the Model without Impulses
In this section, basic properties such as nonnegativity and boundedness, basic reproduction number, and equilibrium points and their stability properties are analyzed.

Non-Negativity and Boundedness
In this section we investigate the non-negativeness of the state variables of the To prove the non-negativity property, we establish the following theorem.
Proof. The first equation of Model (2) can be rewritten as where ξ 1 = βV(t) + µ 1 . Integrating (5), we obtain This implies that E S (t) is nonnegative for all t. For the second equation of Model (2), we have which gives The third equation of Model (2) can be written as where Integrating (8), we obtain This implies that V(t) is non-negative for all t. In a similar way, for the last equation, we can say that which gives The above results show that all the solution trajectories of Model (2) are non-negative for all t > 0.
To verify the boundedness of Model (2) with non-negative initial values, we use the following theorem.
Proof. From the positivity of the solution, we obtain This implies that Now, E(t) = E S (t) + E I (t); then, Hence, we can write lim t→∞ sup E(t) ≤ Π µ 1 . From the third equation of (2), we also have From the last equation of Model (2), we obtain Hence, Therefore, all the solution trajectories that start from R 4 + will enter the region U and never leave it.

Basic Reproduction Number
The basic reproductive number (R 0 ) is useful in estimating the ability of a new pathogen to be transmitted. It is defined as the average number of secondary transmissions from a single infected person. When R 0 is less than one, then the disease (epidemic) does not grow, but if it is greater than one, the disease grows. The basic reproductive number has important implications for disease control. It indicates the level of mitigation efforts needed to bring an epidemic under control [40].
The next-generation matrix method, introduced by Driessche, Pauline, and Watmough in [41], is used to determine the basic reproduction number. For this purpose, we consider the non-negative matrix G and non-negative M matrix H, which represents the production of the new infection and its transportation, respectively. The viral dynamical system of (2) is defined below: Now, the matrix G and H can be written as Additionally, G is non-negative and H is a non-singular M matrix; all eigenvalues of J 4 have positive real parts.
For our system, Therefore, the basic reproduction number, denoted by R 0 , is the spectral radius of the next generation matrix and is obtained as Remark 1. Notice that the basic reproduction number R 0 is proportional to the infection rate β and replication rate p, and inversely proportional to the death rate of the virus µ 3 . Therefore, the disease can be managed by reducing the infection rate and replication rate p or increasing the death rate of the virus. This can be done using antiviral drugs. We adopted impulsive periodic application of antiviral drugs.

Existence of Equilibrium Points
Model (2) has two equilibria: and V * satisfies the equation We have the following Theorem: Theorem 3. When R 0 > 1, one and only one endemic equilibrium P * exists. For R 0 < 1, there may exist two endemic points.
Using Descartes' rule of signs, we can say that there exist a unique endemic equilibrium if b 2 < 0 and two positive endemic equilibrium if b 2 > 0, b 1 < 0 and b 2 1 − 4b 0 b 2 > 0.

Stability of Equilibrium Points
In this section, we discuss the local and global stability of the equilibrium points. For the stability of disease-free equilibrium, we have the following theorem. , 0, 0, 0) is locally asymptotically stable for R 0 < 1; when R 0 > 1, the disease-free system becomes unstable.
Proof. To verify the local asymptotic stability atP, we compute the Jacobian matrix of (2) aroundP as given below The characteristic equation from det(JP − λI 4 ) = 0 is that is, with two eigenvalues λ 1 = −µ 1 < 0 and λ 2 = −µ 4 < 0; the rest of the spectrum is given by the roots of the transcendental equation which suggest that the two roots are negative real roots or have negative real parts. Hence, the disease-free equilibrium is locally asymptotically stable if R 0 < 1, and unstable if R 0 > 1.
We have already proven that when R 0 < 1, the disease-free equilibriumP is locally asymptotically stable. Now, we verity the global stability ofP. For this purpose, we construct the Lyapunov function following [42,43]. We have the following theorem for the global stability ofP.
Theorem 5. The disease-free equilibriumP is globally asymptotically stable if R 0 < 1 and it is a unique equilibrium. Otherwise,P is unstable and a unique endemic equilibrium P * exists.
The proof of Theorem 5 is provided in Appendix A. Now, we analyze the transcritical bifurcation between the disease-free and the endemic equilibrium points. We have the following theorem for this analysis. Theorem 6. The system exhibits forward transcritical bifurcation when R 0 = 1.

Proof.
To prove this theorem, we use the approach used by Castillo-Chavez and Song, 2004 [44] of applying the center manifold theory to analyze the dynamics of Model (2). The variables of Model (2) are transformed as (2) becomes: At R 0 = 1, we choose the bifurcation parameterβ such that Then, the Jacobian matrix of Equation (26) at the disease-free equilibriumP is given by To compute the right eigenvectors, w = (w 1 , w 2 , w 3 , w 4 ) T , we consider the system From Equation (29), we obtain Next, we compute the left eigenvector, v = (v 1 , v 2 , v 3 , v 4 ) from vJ = 0 and the system becomes From Equation (30), we obtain where v 3 is calculated to ensure that the eigenvectors satisfy the condition vẇ = 1. Since the first and fourth component of v are zero, we do not need the derivatives of f 1 and f 4 .
From the derivatives of f 2 and f 3 , the only ones that are nonzero are the following: All the other partial derivatives are zero. The direction of the bifurcation at R 0 = 1 is determined by the signs of the bifurcation coefficients a and b, obtained from the above partial derivatives, given respectively by: Therefore, Model (2) exhibits forward bifurcation at R 0 = 1.

Remark 3.
In case of reinfection, the global stability of the endemic equilibrium P * is not guaranteed when R 0 > 1; this is due to some external factors. In the next subsection, we deal with the global stability of Model (2) at the endemic equilibrium point P * when R 0 > 1. Now, we study the stability of P * . (2) is locally asymptotically stable at P * if R 0 > 1; otherwise, it is unstable.

Theorem 7. Model
Proof. At the endemic equilibrium P * , the Jacobian matrix of Model (2) is given by: The characteristics in λ at the endemic equilibrium P * are given by where Clearly, σ 1 > 0. Thus, using the Routh-Hurwitz criteria, we can say that the equilibrium P * of Model (2) is locally asymptotically stable if the following relations are true: We now analyze the global stability of Model (2) for the endemic equilibrium point P * when R 0 > 1. To show this, we use a Dulac function. We prove the global stability of endemic equilibrium in the following theorem. Theorem 8. When R 0 > 1, Model (2) is globally asymptotically stable at the endemic equilibrium point P * .
The proof of the above theorem is given in Appendix B.

Dynamics of the Drug
To analyze the dynamics of the drug dosing, we first analyze the following sub-system.
. Let τ = t k+1 − t n be the period of the drug dosing. The solution of Equation (35) is In the presence of impulsive dosing, we can obtain the recursion relation at the moments of impulse, as written below: Thus, the concentrations of the drug before and after the impulse are obtained respectively as and Thus, the limiting value of the drug concentration before and after one cycle are and

respectively.
We now require the following definitions and lemmas for this study [45,46]: ; then, we say that B belongs to class B 0 if the following conditions hold: (i) B is continuous on (t n , t k+1 ] × R 5 + , n ∈ N, and for all Λ ∈ R 5 , lim (t,µ)→(t + n ,Λ) B(t, µ) = B(t + n , Λ) exists; (ii) B is locally Lipschitzian in Λ.
The lemmas provided above give the following result: In the above section, we discussed perfect drug dosing. Now, we discuss imperfect drug dosing in the following subsection.

Impact of Imperfect Drug Dosing
Suppose a COVID-19 patient during the treatment stage takes a drug holiday after taking n = n 1 doses. Let us assume a positive quantity D 1 , which denotes the minimum difference between the concentration of the drug after n = n 1 doses and the normal concentration of drug D 2 .
Our consideration is based on the assumption that a patient usually misses his perfect drug dose when he achieves the almost cured stage. Thus, COVID-19 patients can take a drug holiday after taking n 1 doses when the difference between the drug concentrations after n = n 1 doses and normal thresholdÃ is less than a chosen small positive number, D 1 . Thus, we must have A patient may take a drug holiday once the drug concentration reaches a periodic orbit. Suppose h 1 doses are subsequently missed; then, we impose the condition that the drug concentration reach a high level and the patient can realize that the further treatment is highly needed. Now, we assume that the difference between the present drug concentration and its possible maximum response is less or equal to a small positive number ( ). Therefore, the inequality D(t − n 1 +h 1 ) ≥ D max − allows us to find the maximum number of doses a patient can miss.
Suppose a patient has missed h 1 doses. Now, in order to keep the drug concentration above D 2 after n 2 doses are taken, a new condition must be applied that forces the drug concentration level to D 3 away from the D 2 . The condition is that the following that must be satisfied: From (A4), n 2 can be determined. However, the calculations for determining n 2 are complicated. One can see [47] for detailed analysis. We determine n 2 from numerical simulations. (4) Using the result in Lemma 4, we establish the following theorem.

Dynamics of the Impulsive System
Theorem 9. The disease-free periodic orbit (Ẽ S , 0, 0,Ã,D) of Model (2) is locally asymptotically stable ifR Proof. Let the disease-free periodic solution of Model (4) be denoted byP(Ẽ S , 0, 0,Ã,D), whereD with the initial condition D(0 + ) as in Lemma 4. We now test the stability of the disease-free equilibrium point. The variational matrix M(t) at the disease-free periodic orbitP(Ẽ S , 0, 0,Ã,D) is calculated as The monodromy matrix P of the variational matrix M(t) is where I n is the identity matrix.

Numerical Simulation
In this section, we study the mathematical models (2) through numerical simulation. The values of the model parameters used in the numerical simulations are taken from Table 1; some are varied for to study different dynamical regimes.
For the dynamical simulation of the model without drugs (i.e, Model (2)), we take the initial conditions as: H(0) = 4 × 10 5 cells mL −1 , I(0) = 5 × 10 −4 cells mL −1 , V(0) = 300 RNA copies mL −1 , and A = 0 molecules mL −1 . Figure 2 describes the forward bifurcation of Model (2) using Theorem 6 and the parameters value in Table 1. When R 0 > 1, the endemic equilibrium exists. The disease-free state loses its stability when R 0 > 1. The disease-free equilibrium is stable when R 0 < 1, while the endemic equilibrium starts to rise with R 0 > 1 (Theorem 4). Figure 3 represents a region of stability of the equilibrium points in parametric planes. In Figure 3a, as both β and p increase from lower to higher values, the disease equilibrium becomes unstable in the region where R 0 > 1, and the endemic state becomes feasible. In Figure 3b, we observe that as the clarence rate of the viral load is increased, the area of stability of the disease-free state is increased. This can be accomplished using drug dosing. Color represents the value of R 0 . DFE is stable for R 0 < 1 and unstable for R 0 > 1. EE exists and is stable for R 0 > 1. Virus clearance rate (0-1) day −1 [24] α Rate of antibody response from immune cells 0-1 day −1 [48] r Viral particles' rate of neutralization by antibodies 0-1 mL (molecules) −1 day −1 [48] θ Half-maximal simulation threshold 0.5 (RNA copies) ml −1 [48] µ 4 Antibody clearance rate 0.07 day −1 [48] ζ Antibody production rate by drug 6 molecules day −1 gm −1 Assumed µ 5 Decay rate of drug 0.1 mg day −1 Assumed Figure 4 shows the solution trajectories with time. For the set of values used, R 0 is greater than one. That means the system is in an endemic state. Now, as the infection rate is increased, susceptible cells decrease and both the infected cell population and the viral load increase accordingly ( Figure 5).   Table 1.  The numerical solution of Model (4) is plotted in Figure 7. The virus population is increasing, whereas susceptible cell density decreases due to the infection. Without drug application, the antibody response is low (Figure 7b) and, thus, infected cells or viruses are increased.
Numerical solutions of the impulsive system for dosing rates (ω = 100 mg) are plotted in Figure 7 for a fixed impulse interval τ = 7 day. We may conclude that for quick recovery, higher doses should be taken. Figure 8 plots two different intervals of impulses. We can see that the lower interval (τ = 7 day) is capable of achieving disease-free periodic orbit sooner than the higher interval, in τ = 14 days. Additionally, we found that with a higher interval of impulse (τ = 14 days) with higher dosing (200 mg), the system does not converge to the disease-free periodic orbit (Figure 8).  In Figure 9, we can observe the dynamics of the drug for imperfect drug dosing corroborated with the Section 4.1.1. In order to show the impact of taking drug holidays, we find the time at which a SARS-CoV-2 patient can take the required number of doses and then miss the maximum number of doses. From this figure, we obtain the number of possible maximum drug holidays during the treatment period for a fixed drug dose and dosing interval. Taking the drug dose ω = 100 mg and a dosing interval τ = 7 day, the maximum number of holidays is fourteen days; i.e., two doses can be missed after 9 doses (n 1 = 9). This figure also shows that after two consecutive drug holidays (h 1 = 2), if the patient again takes five doses (i.e., n 2 = 5), then the antibody response will gain its previous equilibrium position (periodic orbit).  Figure 9. Dynamics of the drug with and without drug holidays, taking an impulse interval of (τ = 7 days) and a fixed drug dosing of ω = 100 mg with two consecutive drug holidays, h 1 = 2.

Discussion and Conclusions
In this study, we used classical susceptible or uninfected cells, infected cells, and virus population in the presence of adaptive immune responses as a functional form. More attention is given to antibody-response modeling and the role of immune responses against the invading SARS-CoV-2 virus in our respiratory system, which is the primary target area.
We computed the basic reproduction number R 0 for our model. We observed that the model system has two equilibrium points; one is disease-free equilibrium (P) and the other is endemic equilibrium (P * ). The disease-free equilibrium is stable asymptotically when the basic reproduction number R 0 is below the unity. When R 0 is greater than unity, the disease-free equilibrium becomes unstable and endemic equilibrium becomes feasible. Here, R 0 = 1 is the forward transcritical bifurcation point at which the system switches its stability from disease-free to endemic equilibrium.
Finally, we studied the effects of taking the drug in impulsive mode and with holidays during treatment. The numerical simulation of an impulse dosing interval of τ = 7 days and ω = 100 mg dosing rate can achieve a disease-free state in a short time. This study also shows that for a short treatment period, instead of taking the drug at every one-day interval for the entire length of the induction period, it would be better if the patient takes two one-day drug holidays after taking the first twenty-two doses.
In a nutshell, the proposed impulsive mathematical model is functional. It successfully describes the dynamics of SARS-CoV-2 within humans. The results obtained from this study can further guide development of a cost-effective drug regimen for SARS-CoV-2 patients with fewer side effects.  We use (26) for the endemic equilibrium point P * . Then, Then,