Mathematical Modeling and Analysis of the Dynamics of RNA Viruses in Presence of Immunity and Treatment: A Case Study of SARS-CoV-2

The emergence of novel RNA viruses like SARS-CoV-2 poses a greater threat to human health. Thus, the main objective of this article is to develop a new mathematical model with a view to better understand the evolutionary behavior of such viruses inside the human body and to determine control strategies to deal with this type of threat. The developed model takes into account two modes of transmission and both classes of infected cells that are latently infected cells and actively infected cells that produce virus particles. The cure of infected cells in latent period as well as the lytic and non-lytic immune response are considered into the model. We first show that the developed model is well-posed from the biological point of view by proving the non-negativity and boundedness of model’s solutions. Our analytical results show that the dynamical behavior of the model is fully determined by two threshold parameters one for viral infection and the other for humoral immunity. The effect of antiviral treatment is also investigated. Furthermore, numerical simulations are presented in order to illustrate our analytical results.


Introduction
RNA viruses are one of the main causes of human infectious diseases and continue to be a major public health problem worldwide, especially after the emergence of new types of such viruses. For instance, coronavirus disease 2019 (COVID- 19) is an infectious disease caused by a novel coronavirus named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) by the International Committee on Taxonomy of Viruses (ICTV) [1]. Since the first case emerged in Wuhan, China at the end of 2019, the disease has spread rapidly from country to country, causing enormous economic damage and many deaths around the world. According to World Health Organization (WHO) statistics as 18 September 2022 [2], SARS-CoV-2 has infected more than 609 million people worldwide and caused more than 6.5 million deaths.
Despite compliance with preventive measures and the administration of vaccines in several countries, SARS-CoV-2 keeps on spreading. This has led researchers to develop treatments or vaccines to prevent or stop the progression of this disease. Antiviral treatment for COVID-19 remains a challenge. Some drugs have shown promising antiviral activity. This is the case of Paxlovid which has been shown to be highly effective against SARS-CoV-2 with an effectiveness rate of 89% in high-risk patients [3].

Mathematical Formulation
We use a new mathematical model to describe the dynamics of RNA viruses such as SARS-CoV-2 that incorporates two modes of transmission (virus-to-cell and cell-to-cell), two classes of infected cells, humoral immunity and antiviral treatment. The model is formulated by the following nonlinear system of ordinary differential equations (ODEs): where the uninfected cells (S) are generated at rate σ, die at rate µ 1 S and become infected either by free virus particles at rate β 1 SV or by direct contact with infected cells at rate β 2 SI. The two modes of transmission are inhibited by non-lytic humoral immune response at rate 1 + q 1 W and 1 + q 2 W, respectively. The latently infected cells (L) die at rate µ 2 L and become productively infected cells rate δL. Also, the latently infected cells are assumed to be cured at rate ρL, resulting from the clearance of virus through the non-cytolytic process as for HCV infection in [20] and HIV in [21,22]. The cure of infected epithelial cells was also considered in a recent work of SARS-CoV-2 [23]. The productively infected cells (I) die at rate µ 3 I. Free viruses (V) are produced by infected cells at rate kI, cleared at rate µ 4 V and neutralized by antibodies at rate pVW. Antibodies develop in response to free virus at rate rVW and decay at rate µ 5 W. Here, the parameter represents the effectiveness of the antiviral treatment which blocks the production of viral particles. The flow diagram of the model is shown in Figure 1. Most viruses can spread via two modes: by virus-to-cell infection and through direct cell-cell contact [24][25][26]. A recent study provided evidence that SARS-CoV-2 spreads through cell-cell contact in cultures, mediated by the spike glycoprotein [27]. Furthermore, it is known that antibodies neutralize free virus particles and inhibit the infection of susceptible cells [28]. They also contribute significantly to non-lytic antiviral activity [29]. For this reason, both modes of transmission with the lytic and non-lytic immune response are considered into the model.
On the other hand, it is very important to note that the SARS-CoV-2 model presented by system (1) includes many mathematical models for viral infection existing in the literature. For instance, we get the model of Rong et al. [21] when q 1 = 0, β 2 = 0 and both treatment and humoral immunity are ignored. The global stability of the model [21] was investigated in [30]. In addition, the model of Baccam et al. [31] is a special case of system (1), it suffices to neglect immunity and take σ = 0, µ 1 = µ 2 = ρ = 0, β 2 = 0 and = 0. The last model presented in [31] was recently used by Rodriguez and Dobrovolny [32] to determine viral kinetics parameters for young and aged SARS-CoV-2 infected macaques. In the case where latently infected cells not revert back to susceptible and when antibodies do not reduce cell-to-cell transmission, we have ρ = 0, q 2 = 0 and system (1) reduces to the following model:

Equilibria and Threshold Parameters
The equilibria of the model are the stationary solutions of system (1). Then they satisfied the following equations: From Equation (7), we have W = 0 or V = µ 5 r . So, we discuss two cases. (ii) When I = 0, we have S = µ 3 µ 4 (µ 2 + δ + ρ) δk(1 − )β 1 + δµ 4 β 2 . By summing (3) and (4), we The threshold parameter R 0 is called the basic reproduction number. Biologically, this threshold parameter represents the average number of secondary infections produced by one productively infected cell at the beginning of infection. It can be rewritten as R 01 + R 02 , where R 01 = kδσβ 1 (1 − ) is the basic reproduction number of the virus-to-cell transmission mode and R 02 = σδβ 2 µ 1 µ 3 (µ 2 + δ + ρ) is the basic reproduction number of the cell-to-cell transmission mode. When R 0 > 1, model (1) has another biological equilibrium called the infection equilibrium without humoral immunity of the form E 1 = (S 1 , L 1 , I 1 , V 1 , 0), where . On the other hand, we have .
When the humoral immune response has not been established, we have rV 1 − µ 5 ≤ 0. Hence, we define another threshold parameter called the reproduction number for humoral immunity as follows where 1 µ 5 is the average life span of antibodies and V 1 is the quantity of viruses at the steady state E 1 . So, the number R W 1 can biologically determine the average number of antibodies activated by viral particles.
This establishes the uniqueness of S 2 and therefore model (1) has an unique infection equilibrium point with humoral immunity E 2 = (S 2 , L 2 , Summarizing the cases discussed above in the following result.

Analytical Results
This section is devoted to the analytical results including the well-posedness of the model and the stability analysis of equilibria.

Well-Posedness
The model presented by system (1) describes the dynamics of cells and virus which are nonnegative and bounded biological quantities. Since the second member of the system (1) is continuously differentiable, we deduce that the solution of (1) exists and it is unique. Next, we prove the non-negativity and boundedness of the solutions of (1) by establishing the following theorem.
According to Proposition B.7 of [33], we deduce that the solutions S(t), L(t), I(t), V(t) and For the boundedness of solutions, we consider the following function Then This ensures the boundedness of solutions.

Stability Analysis
In this subsection, we analyze the dynamical behaviors of our model. First, we investigate the stability of the infection-free steady state E 0 . Theorem 3. The infection-free steady state E 0 is globally asymptotically stable if R 0 ≤ 1 and becomes unstable if R 0 > 1.
Proof. Let u = (S, L, I, V, W) be a solution of (1) and construct a Lyapunov functional as follows Hence, dH 0 dt ≤ 0 for R 0 ≤ 1 and with equality if and only if S = S 0 , L = 0, I = 0, V = 0 and W = 0. It follows from LaSalle's invariance principle [34] that E 0 is globally asymptotically stable if R 0 ≤ 1.
When R 0 > 1, the characteristic equation of model (1) at E 0 is given by where We have lim Therefore, the characteristic Equation (10) has at least one positive eigenvalue if R 0 > 1, which implies that E 0 is unstable when R 0 > 1.
Next, we establish the asymptotic stability of the two infection steady states E 1 and E 2 by supposing that R 0 > 1 and the following further hypothesis Proof. The characteristic equation at E 1 is given by where Hence, λ 1 = rV 1 − µ 5 is a root of the characteristic equation (11).
Theorem 5. Assume that (H) is satisfied for E 2 . If R W 1 > 1 and µ 1 S 2 − ρL 2 ≥ 0, then the infection steady state with humoral immunity E 2 is globally asymptotically stable. 5 and using the same technique of computation as in H 2 , we obtain

Proof. Consider the following Lyapunov functional
It follows from (12) and the condition µ 1 S 2 − ρL 2 ≥ 0 that dH 2 dt ≤ 0 with equality if and only if S = S 2 , L = L 2 , I = I 2 , V = V 2 and W = W 2 . Hence, E 2 is globally asymptotically stable.
The rate to become productively infected cells, δ, was assumed to be 4 day −1 in [36,37], between 1 and 5 day −1 in [38]. Further, the parameter δ was estimated to be 7.88 day −1 in [39]. Then δ can be estimated between 1 and 7.88 day −1 , which implies that the latently infected cells started producing virus about 3 to 24 h after they were infected.

Sensitivity Analysis
Sensitivity analysis helps understand how changes in model parameters affect the dynamics of SARS-CoV-2 infection. Since the value of the basic reproduction number R 0 determines the eradication or the persistence of the virus in the human body, we so compute the sensitivity index of R 0 against model parameters. By using the explicit formula of R 0 given in (8), such index with respect a parameter q is calculated as follows Note that, by applying above Equation (13) and from the data in Table 1, we conclude that the most sensitive parameters to the basic reproduction number R 0 of the SARS-CoV-2 model are σ, β 1 , β 2 and κ. Clearly, an increase of the value of any one of these parameters will increase the basic reproduction number. While, an increase of the value of other parameters ρ, µ 1 , µ 2 and µ 3 will decrease R 0 ( Table 2). Then, we get the following results and illustrate it in Figures 2 and 3.
Finally, we study the impact of antiviral treatment on of the dynamics of SARS-CoV-2 infection. According to the explicit formula of the basic reproduction number R 0 presented in (8), we notice that R 0 is a decreasing function with respect to . From the numerical results (see, Figure 7), we notice the value of R 0 becomes less than 1 when the effectiveness of the antiviral treatment exceeds 70%. This biologically implies that SARS-CoV-2 infection will disappear from the human body when R 0 < 1.

Conclusions and Discussion
In this study, we have developed a new mathematical model that describes the dynamics of RNA viruses such as SARS-CoV-2 inside the human body. The developed model takes into account the two modes of transmission and both classes of infected cells that are latently infected cells and actively infected cells that produce virus particles. The cure of infected cells in latent period, both the lytic and non-lytic immune response are considered in our model. We first presented the result of the proposed model's well-posedness in terms of non-negativity and boundedness of the solutions. We found two threshold parameters that completely determine the dynamics of this SARS-CoV-2 infection model while looking for biologically feasible equilibria. The first is the basic reproduction number, denoted by R 0 and the second is the reproduction number for humoral immunity, denoted by R W 1 . Our results indicate that the infection-free equilibrium E 0 is globally asymptotically stable when R 0 ≤ 1. This means that SARS-CoV-2 has been completely eradicated from human lungs. However, when R 0 > 1, E 0 becomes unstable and the virus persists in human lungs. Also, two scenarios appear depending on the value of the reproduction number for humoral immunity R W 1 . More precisely, the infection equilibrium without humoral immunity E 1 is globally asymptotically stable if R W 1 ≤ 1. Whereas, when R W 1 > 1, E 1 becomes unstable, and the infection-immunity equilibrium E 2 is globally asymptotically stable. Moreover, the sensitivity analysis showed that the evolution of the infection is influenced by the different parameters of the model. Finally, from the numerical simulations it is deduced that the use of the drugs makes the basic reproduction number R 0 of the patient less than or equal to 1, which automatically leads to the eradication of SARS-CoV-2 from the patient's lungs. This shows the impact of treatment on the dynamics of infection by RNA viruses, in particular SARS-CoV-2. As well as, there are more results that we got it from sensitivity analysis as some parameters influence on the reproduction number. And then we get an accurate understanding of the behavior of the SARS-CoV-2 for example: σ, β 1 , β 2 and k. Clearly, an increase of the value of any one of these parameters will increase R 0 . While, an increase of the value of other parameters ρ, µ 1 , µ 2 and µ 3 will decrease R 0 (see Figures 2 and 3).
Immunologic memory is a major feature of adaptive immunity. It consists in the ability of the immune system to recognize antigens, to which it was formerly exposed, and triggers a more robust and more efficient response than the first exposure [41]. For that reason, it will be more interesting to investigate the effect of immunological memory on the dynamics of the developed model as presented in [8,42,43]. This will be the our future scope.