In-host Mathematical Modelling of COVID-19 in Humans.

COVID-19 pandemic has underlined the impact of emergent pathogens as a major threat to human health. The development of quantitative approaches to advance comprehension of the current outbreak is urgently needed to tackle this severe disease. Considering different starting times of infection, mathematical models are proposed to represent SARS-CoV-2 dynamics in infected patients. Based on the target cell limited model, the within-host reproductive number for SARS-CoV-2 is consistent with the broad values of human influenza infection. The best model to fit the data was including immune cell response, which suggests a slow immune response peaking between 5 to 10 days post-onset of symptoms. The model with the eclipse phase, time in a latent phase before becoming productively infected cells, was not supported. Interestingly, model simulations predict that SARS-CoV-2 may replicate very slowly in the first days after infection, and viral load could be below detection levels during the first 4 days post infection. A quantitative comprehension of SARS-CoV-2 dynamics and the estimation of standard parameters of viral infections is the key contribution of this pioneering work. These models can serve for future evaluation of control theoretical approaches to tailor new drugs against COVID-19.

INTRODUCTION

Epidemics by infectious pathogens are a major threat to humankind. The year 2020 has uncovered one of the biggest pandemics in history, the novel coronavirus SARS-CoV-2 that was first reported in Wuhan, Hubei Province, China in December 2019 (CDC, 2020). While China did a large effort to shrink the outbreak, COVID-19 developed into a pandemic in more than 210 countries moving the epicentre from China to Europe and consequently to America(CDC, 2020). Several countries are planning to relax the strict social distancing regulations. Nevertheless, epidemic rebound risks are latent(Lopez, Rodo, 2020, Ricardo-Azanza, Vargas-Hernandez, 2020).

Coronaviruses are found in different species of animals (e.g. bats and camels) and can evolve to infect humans by droplets from coughing or sneezing. In February 2003, the Severe Acute Respiratory Syndrome (SARS-CoV) was reported in Asia resulting in 8422 cases with a case-fatality rate of 11% (CDC, 2020). Later, in 2012, the Middle East respiratory syndrome (MERS-CoV) was identified in Saudi Arabia with about 2506 cases, killing 862 between 2012 and 2020 (CDC, 2020). Metagenomics studies previous to the COVID-19 outbreak envisaged the possibility of future threats due to the identification of several sequences closely related SARS-like viruses circulating in the Chinese bat populations (He, Zhang, Xu, Yang, Yang, Feng, Xia, Zhou, Zhen, Feng, Guo, Zhang, Tu, 2014, Menachery, Yount, Debbink, Agnihothram, Gralinski, Plante, Graham, Scobey, Ge, Donaldson, Randell, Lanzavecchia, Marasco, Shi, Baric, 2015).

So far, no vaccine or antiviral drug is likely to be available soon. Either monoclonal antibody or vaccine approaches have failed to neutralize and protect from coronavirus infections (Menachery et al., 2015). Therefore, individual behaviour (e.g early self-isolation and social distancing), as well as preventive measures such as hand washing and covering when coughing are critical to control the spread of COVID-19 (Anderson, Heesterbeek, Klinkenberg, Hollingsworth, 2020, Hernandez-Mejia, Hernandez-Vargas, 2020, Prather, Wang, Schooley, 2020). Additionally to these measures, several travel restrictions and quarantines have taken place in many countries around the globe.

Epidemiological models have highlighted that social distancing interventions to mitigate the epidemic is a key aspect (Anderson, Heesterbeek, Klinkenberg, Hollingsworth, 2020, Lopez, Rodo, 2020, Hernandez-Mejia, Hernandez-Vargas, 2020, Ricardo-Azanza, Vargas-Hernandez, 2020). Nevertheless, there are many epidemiological unknowns with the COVID-19 pandemic (Anderson et al., 2020). The case fatality rate for COVID-19 is about 0.3-1% (CDC, 2020). However, adjusted estimations by Baud et al. (2020) indicates that the COVID-19 mortality rate could be as high as 20% in Wuhan. In its early stages, the epidemic has doubled in size every 7.4 days (Li et al., 2020). Moreover, the basic reproductive number was estimated to be 2.2 (95% CI, 1.4 to 3.9) (Li et al., 2020). Based on the relatively long incubation period for COVID-19, about 5–6 days, (Anderson, Heesterbeek, Klinkenberg, Hollingsworth, 2020, CDC) suggested a considerable pre-symptomatic infectiousness.

While there are many mathematical models developed at epidemiological level for COVID-19 to discuss the transmission of SARS-CoV-2 and de-confinement strategies (Anderson, Heesterbeek, Klinkenberg, Hollingsworth, 2020, Ferretti, Wymant, Kendall, Zhao, Nurtay, Abeler-Dörner, Parker, Bonsall, Fraser, 2020, Lopez, Rodo, 2020, Hernandez-Mejia, Hernandez-Vargas, 2020, Peng, Yang, Zhang, Zhuge, & Hong, Ricardo-Azanza, Vargas-Hernandez, 2020), there are too few models at within-host level to understand SARS-CoV-2 replication cycle, interactions with the immune system, and drug effects (Du, Yuan, 2020, Ejima, Kim, Ito, Iwanami, Ohashi, Koizumi, Watashi, Bento, Aihara, Iwami, 2020, Gonçalves, Bertrand, Ke, Comets, de Lamballerie, Malvy, Pizzorno, Terrier, Calatrava, Mentré, Smith, Perelson, Guedj, 2020, Goyal, Cardozo-Ojeda, Schiffer, 2020, Su, Ejima, Ito, Iwanami, Ohashi, 2020, Wang, Heiland, Craig, Davis, Ford Versypt, Jenner, Ozik, Collier, Cockrell, Becker, An, Glazier, Narayanan, Smith, Macklin, 2020). Among different model structures to represent viral dynamics, the target cell limited model has served to represent several diseases such as HIV (Hernandez-Vargas and Middleton (2013); Perelson and Ribeiro (2013); Pinkevych et al. (2016); Rong and à (2009)), Hepatitis (Graw and Perelson (2015); Reluga, Dahari, and Perelson (2009)), Ebola (Nguyen, Binder, Boianelli, Meyer-Hermann, and Hernandez-Vargas (2015); Nguyen and Hernandez-Vargas (2017)), influenza (Baccam, Beauchemin, Macken, Hayden, and Perelson (2006); Handel, Longini, and Antia (2007); Hernandez-Vargas et al. (2014b); Pawelek, Dor, Salmeron, and Handel (2016)), among many others. A detailed reference for viral modelling can be found in Hernandez-Vargas (2019). Very recent data from infected patients with COVID-19 has enlightened the within-host viral dynamics. Zou et al. (2020) presented the viral load in nasal and throat swabs of symptomatic patients. Interestingly, SARS-CoV-2 replication cycles may last longer than flu, about 10 days or more after the incubation period (Anderson, Heesterbeek, Klinkenberg, Hollingsworth, 2020, Zou, Ruan, Huang, Liang, Huang, Hong, Yu, Kang, Song, Xia, Guo, Song, He, Yen, Peiris, Wu, 2020). Here, we contribute to the mathematical study of SARS-CoV-2 dynamics at the within-host level based on data presented by Wölfel et al. (2020).

MODELLING SARS-CoV-2 DYNAMICS IN THE HOST

Using ordinary differential equations (ODEs), different mathematical models are presented to adjust the viral kinetics reported by Wölfel et al. (2020) in patients with COVID-19. ODEs are solved using the MATLAB library ode45, which is considered for solving non-stiff differential equations (Mathworks, 2020).

Note that the viral load (Wölfel et al., 2020) was sampled from throat swab cultures and measured in copies/ml, g Swab, at Log10 scale. The clinical data set of 9 individuals is from Wölfel et al. (2020). Due to close contact with index cases and an initial diagnostic test before admission, patients were hospitalized in Munich (Wölfel et al. (2020)). Viral load kinetics were reported in copies/ml per the whole swab for 9 individual cases. All samples were taken about 2 to 4 days post symptoms. Further details can be found in Wölfel et al. (2020).

Parameter identification aims to estimate model parameters based on minimizing the error between model prediction and experimental data. The viral load is measured in Log10 scales, thus, parameter fitting is performed using the following cost function:

This is the difference on Log10 scales between the model output (), and the experimental measurement (y). This writes as follows where n is the number of measurements.

Turn out that the minimization of a cost function implies a nonlinear optimization problem, complex, with several variables, and multiple minima. This complexity can be tackled using evolutionary optimization algorithms such as the Differential evolution (DE) algorithm (Storn & Price, 1997). Note that several optimization solvers were considered in previous modeling work (Hernandez-Vargas et al., 2014b), including both deterministic (fmincon Matlab routine) and stochastic (e.g Genetic and Annealing algorithm) methods. Simulation results revealed that the DE global optimization algorithm is robust to initial guesses of parameters than other mentioned methods (Torres-Cerna, Alanis, Poblete-Castro, Bermejo-Jambrina, & Hernandez-vargas, 2016).

Note that multiple models can provide the same fit with observed experimental data. Thus, it becomes necessary to choose between different models. The standard approach to model selection is first estimate all model parameters from the data, then select the model with the best-fit error and some penalties on model complexity. A very used model selection criteria is defined next.

The corrected (AIC) writes as follows: where N is the number of data points, M is the number of unknown parameters and RSS is the residual sum of squares obtained from the fitting routine. AIC is used here to compare the goodness-of-fit for models that evaluate different hypotheses (Burnham & Anderson, 2002). A lower AIC value means that a given model describes the data better than other models with higher AIC values. Small differences in AIC scores (e.g. < 2) are not significant (Burnham & Anderson, 2002).

Redundant parametrization in models provides difficulties to estimate the parameters uniquely. Ambiguous parameters θ ⊂ θ may be varied without changing the output y resulting in constant values for the cost function to minimize e.g residuals. This is particularly important for biological systems, where a large variability is presented from one host to another host, limiting the prediction value of mathematical models and estimated parameters (Xia & Moog, 2003).

Xia & Moog, 2003

A mathematical model is identifiable if θ can be uniquely determined from the measurable output y(t); otherwise, the system is unidentifiable.

Note that a mathematical model is algebraically identifiable may still be practically non-identifiable if the amount and quality of the data are insufficient and the data manifest large deviations. The computational approach by Raue et al. (2009) exploits the profile likelihood to determine identifiability and is considered here. This methodology can detect both structurally and practically non-identifiable parameters (Hernandez-Vargas, 2019).

Target cell limited model for SARS-CoV-2

The mathematical model used here to represent SARS-CoV-2 dynamics is based on the target cell-limited model (Ciupe, Heffernan, 2017, Hernandez-Vargas, 2019, Perelson, 2002), which writes as follows:

Equation (3) represents the dynamics of susceptible cells (U) while the equation (4) represents the dynamics of infected cells (I). Note that SARS-CoV-2 can replicate in a variety of cell types, including epithelial cells. In fact, SARS-CoV-2 infection may promote the induction of endotheliitis in several organs as a direct consequence of viral involvement (Varga et al., 2020).

Viral dynamics are represented by (5). Viral particles (V) infect susceptible cells with a rate β ((copies/ml)1 day1). Once cells are productively infected, they release virus at a rate p (copies/ml day1 cell1) and virus particles are cleared with a rate c (day1). Infected cells are cleared at rate δ (day1) as a consequence of cytopathic viral effects and immune responses.

It is still debatable in the literature which compartments SARS-CoV-2 can infect, however, there is a common agreement that the infection mainly takes place in respiratory epithelial cells (Tyrrell & Myint, 1996). Previous mathematical modelling work for human influenza infection has computed the number of target cells in an adult, U(0), is about 4 × 108 cells (Baccam et al., 2006). Initial values for infected cells (I(0)) are taken as zero. Note that V(0) cannot be measured as it is below detectable levels (about 100 copies/ml) (Wölfel et al., 2020). Previous modelling work has suggested using half of the detection levels (less than 50 copies/ml) (Hernandez-Vargas, Wilk, Canini, Toapanta, Binder, Uvarovskii, Ross, Guzman, Perelson, Meyer-Hermann, 2014b, Thiebaut, Guedj, Jacqmin-Gadda, Chene, Trimoulet, Neau, Commenges, 2006). Here, using a regression model, the initial viral concentration V(0) was estimated to be about 0.31 copies/ml.

SARS-CoV-2 kinetics are measured after the onset of symptoms (Wölfel et al., 2020). However, it is unknown when the initial infection took place. Patients infected with MERS-CoV in Oh et al. (2016) showed that the virus peaked during the second week of illness, which indicated that the median incubation period was 7 days (range, 2 to 14) (Oh et al., 2016). The incubation period for SARS-CoV-2 has a median time of 2-5 days from exposure to symptoms onset (Lauer et al., 2020). Therefore, for simulation purposes, we explore different scenarios of initial infection day (t), that is, -7, -3 days before the onset of symptoms.

Infectivity. This is the ability of a pathogen to establish an infection (Diekmann, Heesterbeek, & Metz, 1990). To quantify infectivity, the within-host reproductive number (R 0) is employed to compute the expected number of secondary infections produced by an infected cell (Heffernan, Smith, & Wahl, 2005). When R 0 < 1, one infected cell would infect less than one cell. Thus, the infection would be cleared from the population. Otherwise, if R 0 > 1, the pathogen would be able to invade the target cell population. For the model (3)-(5), the reproductive number is

Using only the viral titters for parameter fitting procedures in the target cell limited model would result in identifiability problems (Nguyen, Klawonn, Mikolajczyk, & Hernandez-Vargas, 2016). That is, parameter values are rescaled and consequently the validity of parameters would be doubtful to extract biological meaning. As the viral clearance is attributed to a process not directly related to the immune system or viral particle per se, the viral clearance parameter (c) is fixed here with previous estimates in humans, e.g approximately 2.4 for influenza and HIV (Baccam, Beauchemin, Macken, Hayden, Perelson, 2006, Hernandez-Vargas, Middleton, 2013).

Assuming the day of infection closer to the post symptom onset (pso), day 0, numerical results show high reproductive numbers (R 0) and higher infection rates (β) as presented in Table 1 . Considering the Remark 1, it is assumed the initial day of infection is -7 or -3 pso, therefore, the rate of infection (β) would be slow but associated with a high replication rate (p). Note that individual parameter values should not be interpreted because of identifiability problems. Fig. 1 reveals that viral replication is below detectable levels from 3 to 4 days post-infection. Independently of the starting infection time (t), numerical results at the Table 1 reveal high reproductive numbers (more than 4.91), implying that SARS-CoV-2 would invade most of the susceptible target cells.

Table 1

Estimations for the target cell model (3)-(5). The “Mean” arrow represents the average of parameter fitting results to the data sets of 9 patients from Wölfel et al. (2020). The 95% confidence intervals (95% CI) of the parameters is computed based on the fitting results for the 9 patients. The parameter c is fixed in 2.4.

ti		β	δ	p	R0	AIC
-3	Mean	4.71×10^−8	1.07	3.07	22.53	11.64
	95% CI	[0.075-21.3]×10^−8	[0.71-1.91]	[0.2-360]	[9.13-70]	[6.81-20.91]
-7	Mean	1.58×10^−8	1.04	5.36	13.51	7.59
	95% CI	[0.03-13.5]×10^−8	[0.61-2.01]	[0.2-362]	[4.5-45.12]	[4.91-20.8]

Fig. 1

Target cell model for SARS-CoV-2. Continuous line are simulation based on the taget cell model (3)-(5). Blue circles represents the data from Wölfel et al. (2020). The most complete data sets to represent the exponential viral growth in Wölfel et al. (2020) were for the patient A in panel (A) and the patient B in panel (B), respectively. Infection time was assumed at -7 and -3 days post symptoms onset.

An additional state known as the “Eclipse Phase” has been proposed by Beauchemin et al. (2008); Holder et al. (2011) to represent the time frame of the infection more adequately. Newly infected cells spend time in a latent phase (E) before becoming productively infected cells (I), this can be written as follows:

Cells in the eclipse phase (E) can become productively infected at rate k. For SARS-CoV-2, we found that the eclipse phase model (AIC ≈ 34) does not improve the fitting respect to the target cell model. This can be attributed to identifiability problems as we only have data for the viral titter (Nguyen et al., 2016).

Model for SARS-CoV-2 and its immune response

Previous modelling studies have acknowledged the relevance of the immune T-cell response to clear influenza virus (Hancioglu, Swigon, Clermont, 2007, Hernandez-Vargas, Wilk, Canini, Toapanta, Binder, Uvarovskii, Ross, Guzman, Perelson, Meyer-Hermann, 2014b, Lee, Topham, Park, Hollenbaugh, Treanor, Mosmann, Jin, Ward, Miao, Holden-Wiltse, Others, 2009, Miao, Hollenbaugh, Zand, Holden-Wiltse, Mosmann, Perelson, Wu, Topham, 2010, Pawelek, Huynh, Quinlivan, Cullinane, Rong, Perelson, 2012, Saenz, Quinlivan, Elton, MacRae, Blunden, Mumford, Daly, Digard, Cullinane, Grenfell, McCauley, Wood, Gog, 2010). Here, we adapted a minimalistic model derived by Almocera, Nguyen, and Hernandez-Vargas (2018); Boianelli et al. (2015) to represent the interaction between influenza and immune response dynamics. The model assumes that the virus (V) level induces the proliferation of T cells (T) as follows:

Equation (11) refers to SARS-CoV-2 dynamic. Viral replication is modelled with a logistic function with a maximum carrying capacity K and a replication rate p. K is the maximum viral load for each of the patients in Wölfel et al. (2020). The initial viral concentration V(0) is 0.31 copies/ml. The virus is cleared at a rate c, which is considered as in the Remark 2. The term c represents the rate of killing infected cells by the immune response.

Equation (12) represents the T cell response against SARS-CoV-2. T cell homoeostasis is represented by where T(0) is the initial number of T cells and δ is the half life of T cells. It is assumed T(0) is about 106 cells. The steady state condition must be satisfied to guarantee the T cell homeostatic value in the absence of viral infection. The half life of T cells is approximately 4-34 days (McDonagh & Bell, 1995), therefore we take . T cells can proliferate at a rate r. It is assumed that the activation of T cell proliferation by the virus follows a log-sigmoidal form with half saturation constant k. The coefficient m relates to the width of the sigmoidal function. While different values of m were tested, rendered the best fit.

The summary of fitting procedures at dpso is presented in Table 2 , which highlights that the mean of AICs values of the model (11)-(12) are smaller than those presented in the target cell limited model. Results of parameter fitting are portrayed in Fig. 2 . Independently of the starting infection day, the immune response by T cells peaks between 5 to 10 dpso as shown in Fig. 3 . Note that this model represents viral clearance only for patients D, H and I. In the other patients the viral clearance was not reported in the data (Wölfel et al., 2020).

Table 2

Estimations for model with immune system (11)-(12) using experimental data from Wölfel et al. (2020) assuming and infection time -3 dpso.

Patient	r	cT	p	kT	AIC
A	0.794	1.58×10^−6	6.31	7.94 × 107	-1.33
B	0.126	7.94×10^−6	12.58	1.99 × 106	11.43
C	0.020	1.58×10^−5	19.95	1.58 × 103	14.28
D	0.251	1.58×10^−6	6.31	3.16 × 104	7.16
E	0.316	1.00×10^−6	5.01	5.01 × 105	16.04
F	0.398	1.26×10^−6	6.31	1.00 × 107	7.07
G	0.158	1.99×10^−6	6.31	7.94 × 104	24.59
H	0.050	2.52×10^−5	31.62	2.52 × 103	10.23
I	0.794	1.58×10^−8	3.98	1.00 × 103	-42.49
Mean	0.194	1.89×10^−6	8.57	1.26 × 105	5.22
95% CI	[0.05-0.79]	[1-15.8] ×  10^−6	[5.01-12.58 ]	[1.58-2000] ×  103	[-1.33-24.59]

Fig. 2

SARS-CoV-2 dynamics with immune responses. Continuous lines are the simulations based on (11)-(12). Blue circles represent the data from Wölfel et al. (2020). Infection time was assumed at -3 days post symptom onset.

Fig. 3

T cell immune response dynamics. Continuous lines are the simulations based on (11)-(12). Infection time was assumed at -3 days post symptom onset.

CONTROL INSIGHTS INTO COVID-19 TREATMENT

The novel coronavirus SARS-CoV-2 first reported in Wuhan in December 2019 has paralysed our societies, leading to self-isolation and quarantine for several days. Indeed, COVID-19 is a major threat to humans, with alarming levels of spread and death tolls, in particular on the elderly. COVID-19 is the first pandemic after the H1N1 ”swine flu” in 2009 (CDC, 2020). While many mathematical models have concentrated on the epidemiological level predicting how SARS-CoV-2 would spread, this paper aims to model SARS-CoV-2 dynamics at the within-host level to quantify SARS-CoV-2 infection kinetics in humans.

Data from Oh et al. (2016) showed that MERS-CoV levels peak during the second week with a median value of 7.21 (log10 copies/mL) in the severe patient group, and about 5.54 (log10 copies/mL) in the mild group. For SARS, the virus peaked at 5.7 (log10 copies/mL) between 7 to 10 days after onset (Peiris et al., 2003). For COVID-19, the viral peak was approximately 8.85 (log10 copies/mL) before 5 dpso (Wölfel et al., 2020). Liu et al. (2020) found that patients with severe disease reported a mean viral load on admission 60 times higher than mild disease cases. Additionally, high viral levels persisted in severe patients for 12 days after onset (Liu et al., 2020).

The reproductive number for human influenza ranges from 3.5 - 75 (Baccam, Beauchemin, Macken, Hayden, Perelson, 2006, Hernandez-Vargas, Colaneri, Middleton, 2014a), which is consistent with the values reported here for COVID-19. Interestingly, both of our models when fitted to the data set of patient A predict that the virus could replicate below detection levels for the first 4 dpi. This may explain why infected patients with SARS-CoV-2 would take from 2-14 dpi to exhibit symptoms.

The model with immune system (Fig. 3) highlights that the T cell response is slowly mounted against SARS-CoV-2 (Anderson et al., 2020). Thus, the slow T cell response may promote low inflammation levels during the first days post infection (Hernandez-Vargas et al., 2014a), which might be a reason to the observations during COVID-19 pandemic of the detrimental outcome on French patients that used non-steroidal anti-inflammatory drugs (NADs) such as ibuprofen. However, so far, there is not any conclusive clinical evidence on the adverse effects by NADs on SARS-CoV-2 infected patients.

Accelerating therapeutic and prophylactic medication against SARS-CoV-2 is unprecedented. Among the different medical strategies, antivirals are central in inhibiting one or several parts of the viral cycle, see Fig. 4 . Network proximity analyses of drug targets by Zhou et al. (2020) prioritized 16 potential anti-SARS-CoV-2 drugs among them melatonin, mercaptopurine, and sirolimus. Additionally, Zhou et al. (2020) identified three potential drug combinations (e.g sirolimus plus dactinomycin, mercaptopurine plus melatonin, and toremifene plus emodin). Based on a rhesus macaque model of MERS-CoV infection, (de Wit et al., 2020) revealed that early remdesivir treatment may provide a clear clinical benefit with a decrease virus replication in the lungs.

Fig. 4

SARS-CoV-2 Replication Cycle. After the binding to receptors of the host cell, the virus RNA is uncoated in the cytoplasm. Then, transcription/translation processes take place to generate new viral RNA material and proteins. Virus assembly occurs within vesicles followed by virus release. Once the virus is released can infect other cells.

Several stages would take place to optimize protocols for the different drugs pointed out against SARS-CoV-2. Mathematical models can help to evaluate in silico the potential of hypothetical drugs. For instance, mathematical terms to represent antiviral effects in the model (3)-(5) writes as follow where u and u would represent the effect of inhibiting the replication cycle and the release of virus, respectively. In this direction, (Gonçalves et al., 2020) found that in order to reduce the peak viral load by more than 2 logs, drug efficacy needs to be greater than 80% if treatment is administered after symptom onset. In a similar direction, the modelling work by Goyal et al. (2020) predicted that to lower viral area under the curve therapies need to be given before the SARS-CoV-2 viral load peaks.

Additionally, mathematical models have served to investigate the effect of immune therapies such as in vivo neutralization of pro-inflammatory cytokines during secondary streptococcus pneumoniae infection post influenza infection (Sharma-Chawla et al., 2019). In the context of COVID-19, our model with immune response (11)-(12) could integrate a term to represent antiviral effects (u) as well as immune modulation (u) to promote the proliferation of T cells, this would be represented as follows:

In other viral infections such as Ebola Nguyen and Hernandez-Vargas (2017), monoclonal antibodies (mAbs) played an important role to inhibit viral replication, e.g the term while increasing viral clearance e.g cV. Another way to suppress SARS-CoV-2 replication is when in co-infections with other viruses because of SARS-CoV-2 has a slower growth rate than the other viruses examined in Pinky and Dobrovolny (2020). Our parameters fitted with the data sets from Wölfel et al. (2020) support also the idea of a slow replication by SARS-CoV-2.

The humoral response against SARS-CoV-2 is urgently needed to evaluate the protection to reinfections. A longitudinal study in rhesus monkeys by Bao et al. (2020) uncovered that infected monkeys presented viral replication at 7 days post-infection (dpi). A significant increase of specific IgG was detected at 14, 21 or 28 dpi. Infected monkeys were re-challenged after specific antibody tested positively and symptoms vanished. Monkeys with re-exposure presented no recurrence of COVID-19, highlighting that protection can be presented to subsequent exposures. Regarding antiviral drugs, Remdesivir treatment has shown a good prophylactic effect during the first 24 hours post MERS-CoV infection in a non-human primate model (de Wit et al., 2020). Furthermore, several benefits were reported for treatment if provided during 12 hours post MERS-CoV infection (de Wit et al., 2020). Our study here mainly addressed T cell responses, therefore, future modelling attempts should be directed to establish a more detailed model of antibody production and cross-reaction (Hernandez-Mejia & Hernandez-Vargas, 2020) as well as in silico testing of different antivirals (Hernandez-Mejia, Alanis, Hernandez-Gonzalez, Findeisen, & Hernandez-vargas, 2019).

There are technical limitations in this study that need to be highlighted. The data for SARS-CoV-2 kinetics in Wölfel et al. (2020) is at the onset of symptoms. This is a key aspect that can bias parameter estimation as mathematical models initiate on the day of infection. In fact, we could miss viral dynamics at the onset of symptoms as well as the SARS-CoV-2 viral load peak. For example, from throat samples in Rhesus macaques infected with SARS-CoV-2, two peaks were reported on most animals at 1 and 5 dpi (Shan, 2020).

In a more technical aspect using only viral load on the target cell limited model to estimate parameters may lead to identifiability problems (Miao, Xia, Perelson, Wu, 2011, Nguyen, Hernandez-Vargas, 2015, Nguyen, Klawonn, Mikolajczyk, Hernandez-Vargas, 2016, Xia, 2003). Thus, our parameter values should be taken with caution when parameter quantifications are interpreted to address within-host mechanisms. For the model with the immune system, there is not data confrontation with immune response predictions, thus, new measurements on cytokines and T cell responses would uncover new information.

The race to develop the first vaccine to tackle COVID-19 has started with the first clinical trial just 60 days after the genetic sequence of the virus. Modelling work developed in this paper paves the way for future mathematical models of COVID-19 to reveal prophylactic and therapeutic interventions at multi-scale levels (Almocera, Hernandez-Vargas, 2019, Feng, Cen, Zhao, Velasco-Hernandez, 2015, Feng, Velasco-Hernandez, Tapia-Santos, Leite, 2012, Handel, Rohani, 2015, Nguyen, Mikolajczyk, Hernandez-Vargas, 2018, Parra-Rojas, Nguyen, Hernandez-Mejia, Hernandez-Vargas, 2018). Further insights into immunology and pathogenesis of SARS-CoV-2 will help to improve the outcome of this and future pandemics.

Ultimately, previous modelling efforts here and from others (Gonçalves, Bertrand, Ke, Comets, de Lamballerie, Malvy, Pizzorno, Terrier, Calatrava, Mentré, Smith, Perelson, Guedj, 2020, Goyal, Cardozo-Ojeda, Schiffer, 2020, Pinky, Dobrovolny, 2020) to represent SARS-CoV-2 could be further extended with control theoretical approaches such as optimal control and model predictive control to schedule drug candidate as well as immune modulators. Control theory in cooperation with Pharmacokinetic (PK)/Pharmacodynamic (PD) modelling have served for optimizing therapies in HIV and influenza infection (Chang, Astolfi, 2008, Hernandez-Mejia, Alanis, Hernandez-Gonzalez, Findeisen, Hernandez-vargas, 2019, Rivadeneira, Caicedo, Ferramosca, Gonzalez, 2017, Rivadeneira, Moog, Stan, Brunet, Raffi, Ferré, Costanza, Mhawej, Biafore, Ouattara, Ernst, Fonteneau, Xia, 2014).

CRediT authorship contribution statement

Esteban A. Hernandez-Vargas: Conceptualization, Formal analysis, Writing - original draft. Jorge X. Velasco-Hernandez: Writing - original draft.

Declaration of Competing Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.