Modeling the viral dynamics of SARS-CoV-2 infection.

Coronavirus disease 2019 (COVID-19), an infectious disease caused by the infection of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is spreading and causing the global coronavirus pandemic. The viral dynamics of SARS-CoV-2 infection have not been quantitatively investigated. In this paper, we use mathematical models to study the pathogenic features of SARS-CoV-2 infection by examining the interaction between the virus, cells and immune responses. Models are fit to the data of SARS-CoV-2 infection in patients and non-human primates. Data fitting and numerical simulation show that viral dynamics of SARS-CoV-2 infection have a few distinct stages. In the initial stage, viral load increases rapidly and reaches the peak, followed by a plateau phase possibly generated by lymphocytes as a secondary target of infection. In the last stage, viral load declines due to the emergence of adaptive immune responses. When the initiation of seroconversion is late or slow, the model predicts viral rebound and prolonged viral persistence, consistent with the observation in non-human primates. Using the model we also evaluate the effect of several potential therapeutic interventions for SARS-CoV-2 infection. Model simulation shows that anti-inflammatory treatments or antiviral drugs combined with interferon are effective in reducing the duration of the viral plateau phase and diminishing the time to recovery. These results provide insights for understanding the infection dynamics and might help develop treatment strategies against COVID-19.

Introduction

In late December 2019, highly contagious pneumonia of unknown etiology was first reported in Wuhan, China [1], [2], [3], [4], [5], [6], [7], [8], [9]. A novel strain of coronavirus was isolated from patients and named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) by the World Health Organization [10]. Since then, SARS-CoV-2 infection has been detected in about 190 countries and caused more than 600,000 deaths worldwide [9], [11], [12], [13], [14], [15]. Infected individuals develop acute respiratory distress syndrome (ARDS) characterized by lymphocytopenia and hyperactive inflammatory response, which becomes the major cause of lung damage and consequent death [16]. High-throughput sequencing demonstrates that SARS-CoV-2 highly resembles the severe acute respiratory syndrome (SARS-CoV) and the Middle East respiratory syndrome (MERS-CoV), both classified as beta coronaviruses identified in bats [17], [18], [19], [20], [21]. Its genome sequence is 96.2% identical to the bat CoV RaTG13 [22]. It uses the angiotensin-converting enzyme 2 (ACE2), the same receptor used by SARS-CoV [23], [24], [25], [26], [27], to infect humans.

Lymphocytopenia is an indicator of the severity and mortality in patients with coronavirus disease of 2019 (COVI-19) [4], [28], [29], [30], [31], [32], [33]. Critically ill patients have significantly low lymphocyte count, possibly due to the direct infection by coronavirus, leading to cytoplasmic damage and destruction of cellular components [30], [31], [33], [34], [35]. In addition, COVID-19 patients have hyperinflammatory responses that may induce massive recruitment of lymphocytes to the lung and result in the consequent pulmonary injury [36]. These findings suggest that SARS-CoV-2 might rely on “cytokine storm” [37] to develop a vicious cycle to drive lymphocytes to lung for infection to maintain the viral replication and transmissibility.

Currently, there are no approved antiviral drugs for SARS-CoV-2 infection. Many patients were given symptomatic management such as empirical broad-spectrum antibiotics or anti-influenza therapies [38]. Mechanical ventilations were used in patients with respiratory failure. The mortality of COVID-19 patients is primarily caused by severe lung injury induced by hyperactive host inflammation. Corticosteroids were administrated to critically ill patients as immunosuppression agents [4], [29]. The adaptive immunity of a patient may be essential for the recovery from SARS-CoV-2 infection [39], [40]. Convalescent plasma containing neutralizing antibodies donated by recovered patients can improve the clinical status of severely ill patients [41]. However, a fraction of patients recovered from COVID-19 still tested positive or even had long-term viral shedding [29], [42]. The mechanisms underlying these pathogenic features are critically needed to mitigate the COVID-19 pandemic.

Some modeling studies have shed light on the transmission dynamics of COVID-19 at the population level [8], [43], [44], [45], [46], [47], [48], [49]. However, the within-host viral dynamics of SARS-CoV-2 infection have not been quantitatively investigated. Here we use mathematical models to study the pathogenic characteristics of SARS-CoV-2 infection by examining the interaction between viral replication and the host immune responses. We fit models to available data of SARS-CoV-2 infected patients and non-human primates in different countries [14], [50], [51], [52]. We also use model simulations to evaluate some potential therapies against COVID-19. Combination of antiviral drugs and type I interferon (IFN) may control the SARS-CoV-2 replication. This is consistent with the observation that SARS-CoV-2 is sensitive to IFN [53]. These findings may provide new insights into the pathogenesis of SARS-CoV-2 infection and development of new treatment strategies.

Methods

Patient and experiment data

The viral load data in COVID-19 patients with and without treatment we analyzed in this paper were from Germany, South Korea and China. We also studied the viral load data from a SARS-CoV-2 infected adult rhesus macaque in the US. The viral load data of eight patients from Germany were from a larger cohort of COVID-19 cases in Munich [54]. Viral RNA levels of all patients were measured every day after onset of symptoms via reverse transcription polymerase chain reaction (RT-PCR) from specimens of throat swab, sputum, and stool. Seroconversion was detected by IgM and IgG immunofluorescence and virus neutralization assay within cells expressing the spike proteins of SARS-CoV-2 [14]. The data from China were collected by RT-PCR from serial samples (throat swabs, sputum, urine and stool) of two patients after they were hospitalized in Beijing [51]. We fit the models to one patient due to the lack of enough measurements in another. The viral load data from South Korea were from a 35-year-old Chinese woman traveling from Wuhan, China [50]. She was the first case reported in South Korea. The viral load was also measured with RT-PCR every day after diagnosis from samples of nasopharyngeal and oropharyngeal swabs and sputum. This patient was treated with lopinavir/ritonavir on day 5 since symptom’s onset. The rectal viral load data from a rhesus macaque having prolonged viral shedding were collected by RT-PCR every day from its rectal swab samples after inoculation [52].

A basic viral dynamic model

We use three models, with increasing complexity, to study the viral dynamics of SARS-CoV-2 infection. The first is a very basic model, which only includes target cells, infected cells and the virus. Variations of the model have been used to study the within-host dynamics of many virus infections, such as HIV, hepatitis, and influenza virus [55]. The second model includes an eclipse phase in the basic model and the third one includes a secondary target cell population of SARS-CoV-2 infection. The basic model is given by the following equations.

The exact target of SARS-CoV-2 infection is not fully clear. Based on the findings from SARS-CoV, MERS-CoV [18], [23], [24], [25], [26], [27], [56] and results from SARS-CoV-2 infection [36], [57], [58], SARS-CoV-2 can infect type II pneumocytes [59]. Target cells (T) are assumed to be infected by SARS-CoV-2 (V) to generate infected cells (I) at a mass action rate VT. Infected cells die with a death rate . They produce coronaviruses at a rate p per infected cell. The virus is cleared at a constant rate c.

SARS-CoV-2 infection can lead to the diffuse alveolar damage (DAD) [60]. DAD usually has three phases: exudative phase, proliferative phase and fibrotic phase [61]. In the exudative phase, immune cells are recruited to the lung, causing accumulated fluid in the airspace. This prevents alveoli from remaining open for gas exchange. The proliferative phase is characterized by regeneration of alveolar type II epithelial cells or pneumocytes  [62]. When replicating, type II epithelial cells transport fluid out of the airspace. SARS-CoV-2 infects pneumocytes and induces hyperactive inflammation that greatly enhances the exudative reaction. Thus, pneumocytes may experience a significant delay in proliferation [63]. In addition, about 30% of COVID-19 patients are asymptomatic [13], [64], [65]. They do not develop DAD so the proliferation of pneumocytes could be at a minimum level. Therefore, when studying the viral dynamics during the early stage of infection, we do not include the proliferation and death of pneumocytes. This also keeps the model in the simplest form for data fitting.

The model with an eclipse phase

The model with an eclipse phase includes four variables: Target cells (T), infected cells in the eclipse phase (), productively infected cells () and viruses (V). Target cells are infected and enter the eclipse phase, in which they do not produce viruses. They progress to productively infected cells at a transition rate k. The other parameters are the same as those in the basic model.

The model with a secondary target of infection

Previous studies [30], [31], [33] detected a large number of viral particles within lymphocytes of SARS-CoV and MERS-CoV patients. A new study also found the infection of T lymphocytes by SARS-CoV-2 [35]. We include lymphocytes as a secondary target cell population in the basic model.

The variable represents the concentration of pneumocytes and represents the concentration of lymphocytes. There is a relatively constant level of lymphocytes in the body of a healthy person. The turnover of uninfected lymphocytes is also slow [66], [67]. Thus, we use a constant rate to represent the recruitment of lymphocytes to the infection site due to the lung inflammatory response. This was also assumed in the viral dynamic models that studied HIV primary infection and influenza virus infection [55], [68]. To minimize the number of parameters for data fitting, we also assume that lymphocytes and pneumocytes are infected by the virus at the same rate . After infection, they are cleared at the same rate (which will be described later) and produce virus at the same rate p. With these assumptions we can use one variable (I) to represent the combined population of infected lymphocytes and pneumocytes. The schematic diagram of model (3) is shown in Fig. 1.

Fig. 1

Schematic diagram of model (3). The secondary target cells () are recruited to the foci at a constant rate . The two target cell populations are infected by SARS-CoV-2 with the same infection rate . Infected cells can be removed by the innate response with a rate . There is also the death rate of infected cells, (t), due to the adaptive immune response activated on day . The death rate (t) increases from the base rate , following an exponential function , where determines how fast the death rate increases.

Lymphocytes primarily consist of T cells, B cells and natural killer (NK) cells. T cells regulate both innate and adaptive immune responses. B cells produce antibodies and NK cells can recognize and kill stressed cells in the absence of antibodies [70], [71], [72], [73]. Because of the lack of data on the innate and adaptive immune responses for model fitting, we assume that the concentration of lymphocytes is proportional to the innate immune response and use the term T2 I to represent the clearance of infected cells by the innate immune response. The adaptive immune response is assumed to be activated on day and the death rate of infected cells increases from the base rate , following an exponential function , where determines how fast the death rate increases. Thus, is when and it is when . Although we can use a saturating function to set the maximum death rate, as used by Pawelek et al. to study the control of influenza virus infection by the adaptive immune response [68], it will introduce an extra parameter for estimation from data fitting. Because we aim to model the viral elimination in the end of infection, the exponentially increasing function is reasonable for this purpose.

Fig. 2

Best fits of model (3) to the LRT viral load data in patients from Germany [14]. The patients #1, #2, #3, #4, #7, #8, #10 and #14 are from the cases in [14], [54]. The detection limit is the blue dotted line according to [14]. We fixed according to the emergence time of seroconversion.

Model parameters

We fix the base death rate of infected cells () to be 2 cells day−1 [68] and the recruitment rate of lymphocyte () to be 104 cells ml−1 day−1 [69]. The total number of alveoli cells within a human is about 6108 cells [74], [75]. Type II pneumocytes constitute 60% of the total alveolar cells [76] and the volume of lung is about 6000 ml. Thus, the initial value of pneumocytes is chosen to be T(0)  (6108 60%)/6000  6104 cells/ml. The initial value of lymphocytes at the infection site is 0. There are no infected cells at the initial time. We fit models to the data of viral load extracted from sputum/saliva. Because much of the early infection occurs in the lung and only a small proportion of the viruses transport to the fluid immediately after infection, we assume V(0) to be an effective initial viral concentration and choose the value of V(0) from 10−6 to 10−3 RNA copies/ml, as used in [77], that can generate the best fits. We estimate the rest of the parameters by fitting models to the data of viral load. All the variables and parameters are summarized in Table 1.

Table 1

Parameters and variables.

Variables/parameters	Description	Value	Ref
T	Target cells	–	–
I_1	Nonproductively infected cells	–	–
I_2	Productively infected cells	–	–
T_1	Primary target cells/pneumocytes	–	–
T_2	Secondary target cells/lymphocytes	–	–
V	Virus	–	–
β	Viral infection rate	Fitted
δ_I	Base death rate of productively infected cells	2 day−1	[68]
k	Transition rate from exposed to productive stage	Fitted
c	Clearance rate of virus	Fitted
p	Viral production rate	Fitted
ω	Killing rate by innate immune response	Fitted
λ	Recruitment of lymphocytes to the infection site	104 day−1	[69]
σ	Parameter in the killing by adaptive immune response	Fitted
μ	Time of emergence of adaptive immune response	Fixed

Model comparison

We compare the data fits with different models by calculating the Akaike information criterion (AIC), given by the following equation [78], [79] where m denotes the number of parameters in the model used to fit data and n is the number of observations (data points). The root mean square (RMS) between the prediction of model and data is calculated using the following formula where V() represents the viral load level at time predicted by the model, is the corresponding data at . Parameter estimates are based on the best fitting that achieves the minimum RMS. Data fitting is performed using the R programming language.

Results

Viral dynamics of SARS-CoV-2 infection

We fit each of the three models Eqs. ((1)–(3)) to the viral load data of SARS-CoV-2 infection in patients from Germany and China [14], [51]. The best fits to the data from lower respiratory tract (LRT) and upper respiratory tract (URT) using model (3) are shown in Fig. 2, Fig. 3, Fig. 4. All the parameters of the best fits are listed in Table 2, Table 3, Table 4. The best fits using models (1) and (2) are shown in the supplemental materials. We find that model (3) significantly improves the fit. The major problem in the fitting with models (1) and (2) is that these models show a single-phase viral decline after the peak. This is not consistent with the viral plateau or even a second peak observed in some patients [14]. Models (1) and (2) also use a constant rate to describe the clearance of infected cells. The killing by the immune responses is complicated. Using a constant rate fails to explain the viral plateau and subsequent decline of the viral load observed in patients [14], [51].

Fig. 3

Best fits of model (3) to the URT viral load data in patients from Germany [14]. The other information is the same as in Fig. 2.

Fig. 4

Best fits of model (3) to the URT (A) and LRT (B) viral load data in a patient from China [51]. The detection limit is the blue dotted line.

Table 2

Best fits of model (3) to the LRT viral load data in patients from Germany.

Patient	β (ml/virus/day)	p (day−1)	c (day−1)	ω (ml/cell/day)	σ (day−1)	V(0) RNA/ml	μ	AIC value
#1	4.8×10−8	7.9×105	17	2.1×10−4	10−3	10−3	10	−4.042036
#2	5.9×10−8	1.3×105	39	1.6×10−4	0.1	10−4	8	−1.708228
#3	3.6×10−6	1.2×103	4.4	10−3	0.1	10−3	9	−10.92767
#4	7.3×10−7	5.7×103	4.4	1.1×10−3	0.5	10−4	6	5.980355
#7	1×10−6	1.1×105	209	4.5×10−4	0.11	10−3	9	17.5142
#8	9.1×10−8	1.2×105	108	4×10−5	0.1	10−4	6	6.927418
#10	4.6×10−7	1.66×102	0.89	10−9	0.93	10−4	15	17.56648
#14	4×10−5	6.5×103	20	1.4×10−2	1.8	10−5	6	0.105858

Table 3

Best fits of model (3) to the URT viral load data in patients from Germany.

Patient	β (ml/virus/day)	p (day−1)	c (day−1)	ω (ml/cell/day)	σ	V(0) RNA/ml	μ	AIC value
#1	1.1×10−5	9.4×102	50	1.9×10−4	0.13	10−4	10	15.41636
#2	2.2×10−6	4.8×103	35	3×10−4	0.01	10−4	8	1.839616
#3	9.8×10−6	1.3×103	58	2×10−4	0.04	10−4	9	−7.151683
#4	4.5×10−5	8.8×10	5	10−3	10−3	10−4	6	8.759159
#7	9.9×10−7	1.08×104	48	2.1×10−4	10−4	10−4	9	7.385558
#8	7×10−6	3.1×103	118	1.5×10−4	0.07	10−4	6	−11.03853
#10	1.8×10−6	4.9×103	3.5	3×10−3	0	10−4	15	15.05884
#14	5.7×10−6	8×102	10	5.3×10−4	0	10−4	6	−0.552598

Table 4

Best fits of model (3) to the URT and LRT data of a patient from China.

	β (ml/virus/day)	p (day−1)	c (day−1)	ω (ml/cell/day)	σ	V(0) RNA/ml	μ	AIC value
URT	7.5×10−7	5×103	12	5.6×10−4	1.5	10−3	6	5.673903
LRT	1.7×10−7	2.9×104	19	2.9×10−4	0.06	10−4	6	2.669148

Numerical simulations of model (3) show that the viral dynamics of SARS-CoV-2 infection can exhibit a few distinct phases (Fig. 5A). During the first stage, the primary target cells are infected and consumed rapidly. The viral load experiences a substantial increase to the peak level. This agrees with the clinical data that onsets of ARDS are observed around this time [14]. Following the viral peak, the viral load declines slightly and enters a plateau phase in which the viral load remains approximately unchanged or declines slowly. In our model, the source of this plateau phase is the infection of lymphocytes as a secondary target. Without the infection of lymphocytes, SARS-CoV-2 would continue to infect and deplete pneumocytes. As a consequence, viral load would decline quickly and the infected individual would lose the transmissibility (Fig. 5A). After viral load persists for a period of time (e.g. a week in some patients), the next phase follows as the adaptive immune response emerges. Seroconversion was also detected about a week after the onset of ARDS in COVID-19 patients [14]. During this stage, viral load declines rapidly to an undetectable level. This also agrees with observation in patients [14].

Fig. 5

Predicted viral dynamics by model (3) under different scenarios. We use the fitting to the LRT viral load data of patient #2 in [14], [54] as an example to study the model prediction under different assumptions. (A) Model prediction with and without lymphocytes (i.e. the secondary target of infection) and seroconversion. The adaptive immune response, if included, is assumed to occur on day 8. (B) The adaptive immune response is assumed to be activated on day 5, 8, 10 and 12, respectively.

Due to the lack of approved antiviral drugs, many patients with severe symptoms are given oxygen treatment, mechanical ventilation, corticosteroids and antibiotics to maintain the respiratory function and ease the inflammatory response. This provides time for the adaptive immune response to emerge to control the virus [4], [29], [37], [38], [39], [40]. In Fig. 5B, we simulate the dynamics of the viral load with adaptive immune response activated at different time. The later the adaptive immune response takes place, the longer the viral load persists. Many senior patients die from COVID-19 because they could not sustain prolonged life support treatment due to multiple significant underlying diseases [80].

The best fits with the model show that the viral infection rate is significantly higher in URT, but the viral production rate and the immune response (described by and ) in LRT are much larger (>10 folds) than that in URT (Table 2, Table 3, Table 4). This is consistent with the route of human-to-human spread of SARS-CoV-2. Virus is transmitted from one host to another via exhaled air, respiratory droplets or fine-particle aerosols [81]. The primary infection of SARS-CoV-2 occurs in the URT while the systematic replication takes place within lung. This suggests that wearing face mask is an effective way to prevent transmission of SARS-CoV-2 from infected individuals [82].

Viral rebound and persistence

Some patients still tested positive after recovering from COVID-19 [42], [83]. The biological process underlying viral rebound is unknown. Our simulation in Fig. 6A shows that if the value of is small, i.e. the clearance rate of infected cells increases slowly as the adaptive immune response emerges, then a second viral peak is predicted because the adaptive immune response is not sufficient to control the viral replication. Thus, viral rebound that lasts for a period of time can be observed. As (t) increases to above a certain level, the viral load will decline again and eventually be eliminated (Fig. 6A).

Fig. 6

Predicted viral rebound and persistence under different scenarios. We use the fitting to the URT viral load data of patient #3 in [14], [54] as an example to study viral rebound and persistence with different assumptions. (A) Predicted viral dynamics with different values of , the parameter determining how fast the death rate of infected cells increases as the adaptive immune response emerges. (B) Predicted viral rebound and persistence assuming the adaptive immune response is activated on day 9, 15, 20 and 25, respectively. The parameter is fixed to be 0.03 day−1.

The timing of the emergence of adaptive immune response also affects viral rebound and persistence. If the adaptive immune system is activated later, then the viral load can persist at a higher level for a longer time (Fig. 6B). This is consistent with a report that a SARS-CoV-2 infected non-human primate with the lowest and latest antibody response had prolonged viral shedding from the intestinal tract [52]. We also fit our model to the rectal viral load data of this animal and the fit is excellent (Fig. 7A). The parameter estimates are shown in Table 5. The finding suggests that for infected individuals with compromised or delayed immune responses, the viral load will rebound following the first peak and even persist for a long time [84]. If the virus persists at a low level around the detection limit, the patient may be asymptomatic or only develop very mild symptoms. However, the patient can still test positive and spread virus. This has been reported in a study that suggested the transmission from an asymptomatic contact in Germany [13].

Fig. 7

Best fit of model (3) to the data of a non-human primate and a patient in South Korea treated with lopinavir/ritonavir. (A) Best fit to long-term viral shedding in an adult rhesus macaque [52]; (B) Best fit to the LRT viral load data of the Korean patient treated with lopinavir/ritonavir [50]; (C) Best fit to the URT data of the same patient in (B).

Table 5

Best fits of model (3) to the rectal viral load data of a non-human primate.

	β (ml/virus/day)	p (day−1)	c (day−1)	ω (ml/cell/day)	σ	V(0) RNA/ml	μ	AIC value
Rectal	4.9×10−6	350	11	5.8×10−5	0.17	10−2	7	−5.835853

Potential treatment

Lopinavir/ritonavir and remdesivir were suggested to be remedy that might inhibit SARS-CoV-2 viral replication [85], [86]. Lopinavir/ritonavir is a combination of protease inhibitors used to treat HIV-1 infection [87], while remdesivir was developed to treat Ebola and Marburg virus infection, acting as an analog of adenosine nucleotide to confuse viral RNA polymerase and evade proofreading by viral exoribonuclease (ExoN) [88]. We assume that they reduce the production of SARS-CoV-2 by a factor (1-) with a constant efficacy . We compare the prediction of model (3) with the LRT and URT viral load data of a patient in South Korea who received lopinavir/ritonavir on day 5 after the onset of symptoms (Fig. 7B, C). The best fits show that the efficacy of lopinavir/ritonavir is 60% effective in reducing viral replication in the LRT while it does not have effect in the URT (Table 6). We also simulate the dynamics with different efficacies (60% in Fig. 8A vs. 90% in Fig. 8B) on 10th, 15th and 20th day, respectively. We find that a higher efficacy of antiviral drugs only makes viral load decline slightly faster. However, early treatment can significantly diminish the time to recovery.

Table 6

Best fits of model (3) to the viral data of a patient treated with lopinavir/ritonavir.

	β (ml/virus/day)	p (day−1)	c (day−1)	ω (ml/cell/day)	σ	V(0) RNA/ml	μ	ε	AIC value
URT	9.6×10−8	7.7×104	4.4	1.8×10−3	1.2	10−4	11	0	−6.860464
LRT	9×10−9	8.5×104	4.3	4×10−5	0.69	10−4	11	0.6	5.335016

Fig. 8

Predicted viral dynamics under treatment. (A) The efficacy of antiviral drug is assumed to be and (B) . We performed simulation using the fit to the LRT viral load data of patient #2 in [14], [54] as an example. The treatment is assumed to be administrated on day 10, 15 and 20, respectively; (C) For comparison with (A) and (B), the value of is increased by 5 folds to represent IFN intensification. The antiviral drug efficacy is the same as that in (A), i.e. ; (D) Predicted dynamics under anti-inflammation therapy. The value of is fixed to be 104, 103, 102 and 10 cells ml−1 day−1, respectively.

SARS-CoV-2 is shown to be more sensitive to type I interferon (IFN) than SARS-CoV due to the absence of two important IFN antagonists [53]. We examine the impact of IFN as treatment intensification with the above antiviral therapy by increasing the killing due to the innate immune response (i.e. the value of in the model) by 5 folds. We find that increasing the killing by the innate immune response can largely suppress the viral replication (compare Fig. 8A,C). Under IFN treatment, the viral load could be significantly reduced with 60% drug efficacy (lopinavir/ritonavir or remdesivir). The result suggests that early treatment with a combination of IFN and other antiviral drugs might be promising against novel coronavirus infection.

The severe acute respiratory syndrome induced by SARS-CoV-2 infection is characterized by over-exuberant inflammatory response. Baricitinib, fedratinib and ruxolitinib that are known to inhibit JAK-STAT signaling pathway might reduce the levels of cytokines [3], [16]. Very recently, the National Institute of Health has started a clinical trial to test the combination of remdesivir and baricitinib in COVID-19 patients (ClinicalTrials.gov Identifier: NCT04401579). To describe the relieved inflammatory response that may be induced from the anti-inflammation treatment, we decrease the recruitment rate of lymphocytes () in the model. We find that the duration of the plateau phase will be significantly reduced as the recruitment rate decreases (Fig. 8D). The result implies that anti-inflammation drugs can accelerate viral decline because of the limited availability of the secondary target of infection.

Discussion

Within-host dynamics of SARS-CoV-2 infection have not been quantitatively investigated. We developed mathematical models to study the viral dynamics and compared modeling prediction with the viral load data in COVID-19 patients and non-human primates from different countries. Our study provides explanations to several pathogenic features of SARS-CoV-2 infection. Firstly, SARS-CoV-2’s very high transmission capacity is likely to come from a viral plateau phase after the peak, which provides a sustainable source of infectious viruses. The duration of the plateau stage is determined by how fast the seroconversion occurs and how strong the immune responses are on the basis of our data fitting and model simulations. Secondly, the model shows that the current best possible control of SARS-CoV-2 infection might be the adaptive immune response of the infected individual. This supports the legitimacy of current treatment with oxygen, mechanical ventilation, corticosteroids and antibiotics that let severely ill patient sustain the life till the activation of adaptive immunity. Thirdly, we explained the scenario of prolonged viral shedding and viral rebound. They originate from the oscillation caused by interplays between viral replication and the adaptive immune response. If seroconversion emerges late, then prolonged viral persistence is generated. As the adaptive immune response is above a threshold, the virus will be eradicated. These results may improve our understanding of the pathogenic features of SARS-CoV-2 infection.

Comparison of the prediction by the three models and the viral load data suggests that SARS-CoV-2 infection might have a secondary target. Without the second target cell population, the model shows a single-phase viral load decline after the peak, which does not agree with the data in many patients. The secondary target of SARS-CoV-2 infection might be lymphocytes. Lymphocytopenia is a predominant feature of COVID-19 patients in which a large number of viral RNAs within lymphocytes were detected [31], [33]. A recent study also found that SARS-CoV-2 can invade T lymphocytes via spike protein mediated membrane fusion [35]. The molecular mechanism underlying the entry of SARS-CoV-2 into lymphocytes is not fully understood. Lymphocytes do not express ACE2, but SARS-CoV-2 has a polybasic furin-type cleavage site at the S1–S2 junction in the spike protein, which can intensify the capacity of virus to penetrate tissues with low ACE2 expression [14], [89]. Lymphocytes may also be infected via the assistance by other co-receptors. The tropism of SARS-CoV could be altered towards human immune cells enhanced by Fc gamma II receptor [90]. The initiation of seroconversion in COVID-19 patients is slightly earlier than in SARS-CoV infected patients [14]. This might facilitate the entry of SARS-CoV-2 into lymphocytes in humans. In addition to lymphocytes, SARS-CoV-2 may also infect some other cells. The virus replicates to a level in the lung where it hits a threshold for which virus spills over into the blood, causing a viremia. This viremia will allow the virus to travel to the gut where it can further replicate, as gut enterocytes have the ACE2 receptor and have been reported to be infected in human COVID-19 infection [91].

Type I IFN was used to treat MERS-CoV, SARS-CoV, and SARS-CoV-2 infection [53], [92], [93]. The early treatment of MERS-CoV infection protected mice from lethal infection but late induction could result in serious pneumonia [92]. SARS-CoV-2 was shown to be more sensitive to type-I IFN treatment than SARS-CoV [53], [93]. Our simulation shows that the antiviral drug efficacy needs to be close to 90% to suppress viral replication effectively. However, intensification with interferon can more effectively reduce the duration of viral plateau with a 60% efficacy of the antiviral drug. Comparison of the genomic sequences of SARS-CoV-2 and SARS-CoV shows that some changes, such as the absence of open reading frame (ORF) 3b and a short truncation of ORF6 in SARS-CoV-2, can reduce the virus’ capacity to interfere with type I IFN [53], [93]. ORF3b is a 154 amino acid (AA) long protein that blocks IRF3 phosphorylation [94]. SARS-CoV-2 ORF3b has a premature terminal codon that produces a truncated protein of 22 AA. The N-terminal domain of ORF6 can disrupt karyopherin transportation of transcriptions factors like STAT1 [53], [95]. SARS-CoV-2 has a two amino acid truncation in its ORF6 [53]. These alterations might make SARS-CoV-2 infection more sensitive to IFN. IFN was also administrated with other drugs to treat COVID-19 patients but the outcomes have not been available [96]. Because IFN might activate immunity to attract more immune cells, further studies are needed to determine if IFN can exacerbate the cytokine storm and evaluate the effectiveness of using IFN with and without other antiviral drugs to treat SARS-CoV-2 infection in humans.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.