Modeling Dynamics of RNA Virus Transmission

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): $$\left\{\begin{array}{ccc}\frac{dS}{dt}\hfill & =\hfill & \sigma -{\mu }_{1}S-{\displaystyle \frac{{\beta }_{1}SV}{1+q_{1}W}}-{\displaystyle \frac{{\beta }_{2}SI}{1+q_{2}W}}+\rho L,\hfill \\ \frac{dL}{dt}\hfill & =\hfill & {\displaystyle \frac{{\beta }_{1}SV}{1+q_{1}W}}+{\displaystyle \frac{{\beta }_{2}SI}{1+q_{2}W}}-({\mu }_2+\delta +\rho )L,\hfill \\ \frac{dI}{dt}\hfill & =\hfill & \delta L-{\mu }_{3}I,\hfill \\ \frac{dV}{dt}\hfill & =\hfill & k(1-ϵ)I-{\mu }_{4}V-pVW,\hfill \\ \frac{dW}{dt}\hfill & =\hfill & rVW-{\mu }_{5}W,\hfill \end{array}\right.$$
where the uninfected cells $\left(S\right)$ are generated at rate $\sigma$ , die at rate ${\mu }_{1}S$ and become infected either by free virus particles at rate ${\beta }_{1}SV$ or by direct contact with infected cells at rate ${\beta }_{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 $\left(L\right)$ die at rate ${\mu }_{2}L$ and become productively infected cells rate $\delta L$ . Also, the latently infected cells are assumed to be cured at rate $\rho L$ , resulting from the clearance of virus through the non-cytolytic process as for HCV infection in [20 (link)] and HIV in [21 (link),22 (link)]. The cure of infected epithelial cells was also considered in a recent work of SARS-CoV-2 [23 (link)]. The productively infected cells $\left(I\right)$ die at rate ${\mu }_{3}I$ . Free viruses $\left(V\right)$ are produced by infected cells at rate $kI$ , cleared at rate ${\mu }_{4}V$ and neutralized by antibodies at rate $pVW$ . Antibodies develop in response to free virus at rate $rVW$ and decay at rate ${\mu }_{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 (link),25 (link),26 (link)]. A recent study provided evidence that SARS-CoV-2 spreads through cell-cell contact in cultures, mediated by the spike glycoprotein [27 (link)]. Furthermore, it is known that antibodies neutralize free virus particles and inhibit the infection of susceptible cells [28 (link)]. 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 (link)] when $q_1=0$ , ${\beta }_2=0$ and both treatment and humoral immunity are ignored. The global stability of the model [21 (link)] was investigated in [30 (link)]. In addition, the model of Baccam et al. [31 (link)] is a special case of system (1), it suffices to neglect immunity and take $\sigma =0$ , ${\mu }_1={\mu }_2=\rho =0$ , ${\beta }_2=0$ and $ϵ=0$ . The last model presented in [31 (link)] was recently used by Rodriguez and Dobrovolny [32 (link)] 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 $\rho =0$ , $q_2=0$ and system (1) reduces to the following model: $$\left\{\begin{array}{ccc}\frac{dS}{dt}\hfill & =\hfill & \sigma -{\mu }_{1}S-{\displaystyle \frac{{\beta }_{1}SV}{1+q_{1}W}}-{\beta }_{2}SI,\hfill \\ \frac{dL}{dt}\hfill & =\hfill & {\displaystyle \frac{{\beta }_{1}SV}{1+q_{1}W}}+{\beta }_{2}SI-({\mu }_2+\delta )L,\hfill \\ \frac{dI}{dt}\hfill & =\hfill & \delta L-{\mu }_{3}I,\hfill \\ \frac{dV}{dt}\hfill & =\hfill & k(1-ϵ)I-{\mu }_{4}V-pVW,\hfill \\ \frac{dW}{dt}\hfill & =\hfill & rVW-{\mu }_{5}W,\hfill \end{array}\right.$$