Modeling the transmission dynamics and the impact of the control interventions for the COVID-19 epidemic outbreak

1.   Introduction

In late December 2019, the World Health Organization (WHO) received notification of up to 27 possible cases of pneumonia of unknown etiology, including 7 severe cases, in the Chinese city of Wuhan. Within a few days, the novel coronavirus, SARS-CoV-2, provisionally named 2019-nCoV, was identified as the causative agent. Since the first report in Wuhan, China, many countries have now reported cases of infection, affecting people of all ages from different origins. Most people with coronavirus disease 2019 (COVID-19), will experience mild to moderate respiratory illness and recover without requiring special treatment. The most common symptoms at the onset of COVID-19 illness are fever, cough, and fatigue, while other symptoms include sputum production, headache, diarrhea, dyspnoea, and lymphopenia, see ^[1] and the references therein. Older people and those with underlying medical problems are more likely to develop serious illness.

On January 30, 2020, the WHO declared COVID-19 as an international emergency and on March 11, 2020, the WHO declared the global COVID-19 outbreak a pandemic, pointing to the over 118,000 cases of the coronavirus illness in over 110 countries around the world and the constant risk of further spread ^[2]. The COVID-19 pandemic was confirmed to have reached Mexico in February 2020. On February 28, Mexico confirmed its first three cases. According to the WHO, Mexico entered Phase 2 of the coronavirus pandemic on March 23, 2020, with 367 confirmed cases. Phase 2 includes cases where the sick individuals did not have direct contact with someone who had recently been in another country. As of April 18, there had been 7497 confirmed cases of COVID-19 in Mexico and 650 reported deaths.

According to the Center for Disease Control (CDC), the main transmission route for COVID-19 is from person-to-person, either among people in close proximity or through respiratory droplets produced when an infected person coughs or sneezes. Although it is not precisely known the importance of infections caused by contact with contaminated surfaces, the environment-to-human transmission route is also possible, so a person can get COVID-19 by touching a surface or object that has the virus on it and then touching their mouth, nose or eyes ^[3]. Extensive measures to reduce both person-to-person and environment-to-human transmission of COVID-19 are essential to control the current outbreak. Several countries, including China and the US, have implemented major control interventions, including travel bans and airport screening. However, the impact of such interventions is probably minor on COVID-19 containment given the potentially large number of asymptomatic individuals and the possibility of transmission before the onset of symptoms ^[4].

Analysis of epidemiological changes in COVID-19 infection is of paramount importance to boost awareness and public health efforts to control the COVID-19 outbreak. In recent years, mathematical modeling has become a valuable tool for the analysis of dynamics of infectious disease and for the support of control strategies development ^[5]. Mathematical and statistical models are especially useful to estimate key epidemiological parameters such as the basic reproduction number, $\mathcal{R}_{0}$, which is an indicator of the potential severity of an epidemic and provides a powerful tool to estimate the control effort needed to eradicate the disease. Several models, most of them using extensions of the Susceptible-Exposed-Infected-Recovered (SEIR) structure, have been proposed to investigate the spread of COVID-19 in different regions ^[4,6,7,8,9]. In ^[10], the authors review current estimates for the basic reproduction number of COVID-19 from 1 January 2020 to 7 February 2020. They found that the estimates range from $1.4$ to $6.49$, with a mean of $3.28$, a median of $2.79$ and an interquartile range of $1.16$.

In this study, we use a mathematical model to investigate the dynamics of the on-going epidemic outbreak of COVID-19. The rest of the paper is structured as follows. In the next section, we formulate our model and develop the analysis to compute the basic reproduction number. In Section 3, we calibrate our model using a Bayesian approach and officially published data by the Secretariat of Health, Mexico, corresponding to the daily cumulative cases of infected individuals. In Section 4, we use extensive numerical simulations to investigate the theoretical impact of several control interventions against the spread of COVID-19 and compute the effective reproduction number. The last section contains a discussion of the obtained results.

2.   Model formulation without control

Based on the clinical progression of the disease, we propose a deterministic compartmental epidemic model under the SEIR structure. One important aspect in our model is that, in addition to human-to-human transmission, we consider the indirect infections caused by contact with a contaminated environment.

For our model formulation, we divide the total human population (denoted $N$) into five compartments: susceptible individuals (denoted $S$), exposed/latent individuals (denoted $E$), infectious asymptomatic individuals (denoted $A$), infectious with symptoms (denoted $I$), and recovered (denoted $R$). Finally, we consider a compartment for the free-living SARS-CoV-2 in the environment (denoted $V$).

For our model formulation, we consider a short time horizon in which the total human population is relatively fixed. Therefore, demographic dynamics are not considered in the model. The susceptible population $S$ can acquire the infection when they come in contact with asymptomatic $A$ and symptomatic $I$ infectious individuals at rates $\beta_{A}$ and $\beta_{I}$, respectively. They also can be infected by contact with contaminated surfaces with coronavirus at a rate $\beta_{V}$. A proportion $p$ of the exposed individuals $E$ will transition to the symptomatic infectious class $I$ at a rate $\sigma$, while the other proportion $1-p$ will enter the asymptomatic infectious class $A$. The recovery rates for individuals in the classes $A$, $I$ are $\gamma_{A}$, $\gamma_{I}$, respectively. These individuals gain permanent immunity and move to the recovered class $R$. However, individuals in the symptomatic infectious class $I$ can die due to the disease at a rate $\mu$. Asymptomatic and symptomatic infected individuals release virus into the environment with shedding rates $c_{1}$ and $c_{2}$, respectively. Hence, the free-living virus in the environment grows with a factor $c_{1}A+c_{2}I$. The parameter $\mu_{V}$ represents the mortality rate of the free-living virus in the environment.

These assumptions lead to the following system of differential equations:

where $\lambda = \beta_{A}A + \beta_{I}I + \beta_{V}V$ is the force of the infection.

According to the WHO, the SARS-CoV-2 is primarily transmitted between people through respiratory droplets and contact routes. Droplet transmission occurs when a person is in close contact (within 1 m) with an infectious individual and is therefore at risk of having his mouth, nose, or eyes exposed to potentially infectious respiratory droplets ^[3]. The parameters $\beta_k$ ($k = A, I$) model this direct person-to-person transmission and are of the form $\beta_{k} = b_{k}\phi_{k}$, where $b_{k}$ is the average number of contacts per person per unit of time and $\phi_{k}$ is the probability of successful infection given a contact. For example, $\beta_{A} = b_{A}\phi_{A}$, where $b_{A}$ is the average number of close contacts in which a susceptible is exposed to respiratory droplets produced when an asymptomatic infected person coughs, sneezes or talks, and $\phi_{A}$ is the probability of successful infection given this contact. Since the virus can survive on inanimate surfaces ^[11], transmission may also occur through contaminated fomites. The parameter $\beta_{V}$ models these indirect infections caused by touching an object or surface contaminated (due to an infected person) with the virus. Therefore, $\beta_{V} = b_{V}\phi_{V}$, where $b_{V}$ is the average number of times a susceptible person touches a surface contaminated with SARS-CoV-2, and $\phi_{V}$ is the probability of infection given this contact. The parameters $c_{i}\geq 0$ ($i = 1, 2$) measure the number of virus particles released through respiratory droplets produced per infected individual (during coughs or sneezes) that remain alive and infectious on surfaces or objects per unit of time.

We remark that for the starting model (2.1), we are not including the current intervention measures against COVID-19. This will allow us to focus first on the predictions of the model without control. In Section 4, we incorporate control interventions into our model and investigate the extent of the influence of the controls to prevent SARS-CoV-2 spread comparing with the case without control.

2.1.   Disease-free equilibrium and the basic reproduction number $\mathcal{R}_{0}$

The biologically feasible region for model (2.1) is

Let $X(t)$ be the solution of system (2.1) for a well-defined initial condition $X(0)\in \Omega$. Since $X_{i} = 0$, implies $\dot{X}_{i}\geq 0$ for any state variable, then $X(t)\in \Omega$ for all $t > 0$. Thus, solutions trajectories satisfy the usual positiveness and continuity properties and the model is both epidemiologically and mathematically well posed ^[12].

To compute the coordinates of the disease-free equilibrium, we set the rate of change of all state variables equal to zero. Solving the system of algebraic equations we find a unique disease-free equilibrium with the following coordinates:

where $N_{0}$ is the value of the total population at equilibrium.

To compute the basic reproduction number $\mathcal{R}_{0}$, we use the next-generation operator introduced by Diekmann et al. ^[13]. Under this approach, it is necessary to study the subsystem that describes the production of new infections and changes among infected individuals. The Jacobian matrix ${\bf{J}}$ of this subsystem at the disease-free equilibrium is decomposed as ${\bf{J}} = {\bf{F}}-{\bf{V}}$, where ${\bf{F}}$ is the transmission part and ${\bf{V}}$ describe changes in the infection status. The next-generation matrix is defined as ${\bf{K}} = {\bf{F}}{\bf{V}}^{-1}$, and $\mathcal{R}_{0} = \rho({\bf{K}})$, where $\rho(\cdot)$ denotes spectral radius.

For system (2.1), we obtain

Therefore, the basic reproduction number is given by

To interpret the biological meaning of the basic reproduction number (2.4), we need the following components. During his infection period, $1/\gamma_{A}$, an asymptomatic infectious individual produces on average $\beta_{A}S_{0}$ infections and $c_{1}$ virus particles into the environment. Since the coronavirus survives in the environment a mean time of $1/\mu_{V}$, the average number of infections caused by the virus is $\beta_{V}/\mu_{V}$. Hence,

measure the contribution of asymptomatic infectious individuals to the production of new infections taking into account the environment-to-human transmission route for virus released by asymptomatic individuals. Analogously,

is the contribution of symptomatic infectious individuals to the production of new infections. Therefore, the basic reproduction number (2.4) is the weighted sum of the terms $T_{A}$ and $T_{I}$, that is,

As a consequence of Theorem 2 in ^[14], we establish the following result regarding the local stability of the disease-free equilibrium.

Corollary 1. The disease-free equilibrium of system (2.1) is locally asymptotically stable for $\mathcal{R}_{0} < 1$ and unstable for $\mathcal{R}_{0} > 1$.

In this study, we are interested in the early dynamics of the infection process. Therefore, we did not consider demographic dynamics and the study of the asymptotic behavior for endemic equilibria.

3.   Parameter estimates

The compartmental epidemic model (2.1) for the transmission dynamics of SARS-CoV-2 has $11$ parameters. First, we gather some parameter values from the literature. Next, we estimate those parameters that are not found in the literature or that depend on the population under study. We assume the time unit is days and estimate the parameters as follows.

(ⅰ) Recovery rates. The estimated mean value for the recovery rates $\gamma_{A}$, $\gamma_{I}$, for asymptomatic and symptomatic infectious individuals, respectively, have been estimated to be $\gamma_{A} = 0.13978$ and $\gamma_{I} = 0.33029$ ^[8].

(ⅱ) Mean incubation period. The mean incubation period ($1/\sigma$) for coronavirus infection has been estimated to be $6.4$ days, ranging from 2.1 to 11.1 days ^[15]. Therefore, we assume $\sigma = 1/6.4$.

(ⅲ) Fraction of individuals which develop symptoms. The probability of having symptoms after the infection has been estimated to be $p = 0.868343$ ^[8].

(ⅳ) Mortality rate of coronavirus in the environment. Some studies have estimated that coronaviruses can remain infectious on inanimate surfaces at room temperature from a few hours up to 9 days ^[11]. Here, we assume an average survival rate of $1$ day which implies $\mu_{V} = 1$.

(ⅴ) Disease induced death rate. The estimated mean value for the disease induced death rate is $\mu = 1.7826\times 10^{-5}$ ^[8].

The rest of the parameters, that is, the transmission rates $\beta_{A}$, $\beta_{I}$, and $\beta_{V}$, in addition to the shedding rates $c_1$ and $c_2$, will be estimated using Bayesian inference. We focus on this set of parameters for the estimation because transmission parameters depend highly on population-level characteristics and it can be unreliable to take estimations from different data.

We consider data corresponding to the daily cumulative cases of infected individuals in Mexico. The data were obtained from the daily report of the Mexican Secretaria of Health from March 11, 2020, to March 25, 2020 ^[16]. It is important to remark that this data corresponds to the confirmed cases; therefore, it is highly possible that the real epidemic curve is higher than the total infected cases presented in the data. In other words, since in México there is not massive testing, the initial data on the confirmed cases from the pandemic corresponds to symptomatic infections. Therefore, as an attempt to avoid estimates biased down, we fit the data using only the individuals in the symptomatic infectious class, $I$, without considering the asymptomatic infectious class $A$.

The following values were taken as initial conditions: the initial total population was taken as the approximate Mexican population at the year 2020, i.e., $N(0) = 128,000,000$; the initially symptomatic infectious individuals as $I(0) = 4$, which is equal to the initial number of confirmed cases in the data. No recovered individuals are considered at the initial time, thus $R(0) = 0$. Finally, we assumed $E(0) = 4$, $A(0) = 1$, $V(0) = 10$, and $S(0) = N(0)-E(0)-A(0)-I(0)$.

For the parameter inference, we use a Bayesian approach. We run a Markov Chain Monte Carlo (MCMC) using twalk, introduced in ^[17] (see Appendix A for details). We consider the time in days and $t_0 = 0$ for the first data on March 11. The resulting total infected cases for the maximum a posteriori (MAP) and the posterior mean estimates are shown in Figure 1(a). The corresponding values for the parameter estimates are presented in Table 1. The results show that the transmission rate for symptomatic infectious, $\beta_I$, is greater than the transmission rate for asymptomatic infectious individuals, $\beta_A$. This result can be counter-intuitive considering that people may get more contacts with asymptomatic infectious than symptomatic infectious. Nevertheless, this only implies $b_A > b_I$, but the probability of infection is way higher for symptomatic infectious because they have higher infectiousness, thus $\phi_I > \phi_A$ which explains why $\beta_I$ is greater than $\beta_A$. Note that both estimates values for $\beta_V$ almost match. The posterior distribution on $\beta_V$ ensures the identifiability of this parameter (see Appendix A, Figure 9(c)) which provides support for the modeling approach around the state variable $V$ and the environment-to-human transmission route.

In Figure 1(b), we show the values for $\mathcal{R}_0$ corresponding to the parameter estimates and the elements of the chain. In particular, the value of the basic reproduction number for the MAP estimate is $\mathcal{R}_{0}^{MAP} = 2.5$, and for the posterior mean estimate is $\mathcal{R}_{0}^{CM} = 2.7$; hence, $\mathcal{R}^{MAP}_0 < \mathcal{R}^{CM}_0$. These values are in the range of the current $\mathcal{R}_{0}$ estimates ^[10]. Moreover, please observe the heavy tail to the right of these values; this heavy tail implies that there exist possible scenarios with higher $\mathcal{R}_0$ values (see Figure 2(b)) for the uncertainty region on the fitted data. Note that even when the estimated values of $\beta_A$ and $c_1$ have different orders of magnitude for the MAP and the posterior mean, the corresponding values of $\mathcal{R}_{0}$ do not change significantly. This is due to the low value of the assigned weight $(1-p)$ of the asymptomatic infectious class contribution $T_{A}$ to the production of new infections. Moreover, using the MCMC and the analytic expression for $R_0$, we estimate that, on average, the environment-to-human transmission route contributes with $36\%$ of the value of $\mathcal{R}_{0}$. Finally, in Figure 2(a) we present the curves of the symptomatic infected class corresponding to the MAP and posterior mean estimates, $I(t)$, for a time horizon of 200 days.

The posterior distribution obtained allow us to compute posterior predictive marginals for future data (after March 25). The probability of a future observation ${\bf{z}}$ given the data $y$ is $p({\bf{z}}|y)$ and can be computed as follows

where ${\bf{x}}$ denotes our vector of parameters. Figure 3(a) shows the predictive posterior marginal for the total cumulative infections on March 31 (See Figure 10 for more predictive posterior marginals on other dates). Figure 3(b) shows a comparison between the predicted values for the cumulative number of infections and the officially published data from March 26 to April 9. The dashed lines represent the interval with 98 percent of the mass for the predictive marginal.

According to our estimations, the value of the basic reproduction number $\mathcal{R}_{0}$ in the absence of control is above unity. Under a non-intervention scenario, we expect that the number of individuals in the infectious class $I(t)$ to have a high peak that can produce a collapse in the health care system. Therefore, it is of paramount importance the application of effective control measures to limit the spread of SARS-CoV-2 and flatten the epidemic curve.

4.   Control interventions

In this section, we extend the compartmental epidemic model for COVID-19 transmission dynamics (2.1) including appropriate compartments to take into account some of the current intervention measures for COVID-19 control. In particular, we consider the following control interventions:

(ⅰ) Social distancing and home quarantine.

(ⅱ) Isolation of infected individuals.

(ⅲ) Environmental cleaning and disinfection.

We want to differentiate susceptible individuals depending on the probability of being infected. We assume that part of the population changes behavior due to different reasons, including external measures. That change can be because they are following rules of social distancing or because they are afraid and naturally keep a greater distance from others, etc. In any case, we consider two types of behavior and obtain two categories of susceptible individuals. Hence, we assume susceptible individuals $S$ change their behavior and become cautious susceptible individuals (denoted $S_{c}$) at a rate $\alpha$. The parameter $\alpha$ is the rate of behavioral change. This may be increased through mass communication (TV, social networks, etc.). Cautious susceptible individuals will reduce their probability of infection by a factor $1-\theta\in(0, 1)$ taking appropriate measures such as self-imposed home quarantine, social distancing, hand washing, and mask-wearing.

To model isolation, we assume that symptomatic infected individuals, $I$, are screened at a rate $d_{2}$ and moved to a diagnosed compartment $D$. Likewise, individuals in the exposed and asymptomatic classes are diagnosed at a rate $d_{1}$. It should be easier to identify infected people with strong symptoms in comparison to asymptomatic individuals, therefore, $d_{2} > d_{1}\geq 0$. We assume individuals in the $D$ class are being isolated and treated. Finally, we consider cleaning of visibly dirty surfaces followed by disinfection which is an important practice measure for the prevention of COVID-19. We model this by considering an additional mortality rate $m$ for the free virus $V$.

From the above considerations, the control model for the transmission dynamics of COVID-19 is governed by the following equations:

where $\lambda = \beta_{A}A + \beta_{I}I + \beta_{V}V$ is the force of the infection.

It is important to remark that if factors such as lockdown or some external measures affect social behavior of susceptible individuals, one can easily think that this may also happens in other compartments. Hence, regardless of their infection status, individuals in all classes may reduce infecting others. Therefore, from a modeling perspective, it would appear that not only in the susceptible class $S$, but in all other classes, individuals should be allowed to become cautious. In particular, if individuals in the infectious classes become cautious, their contact rates indeed change, but this is reflected mainly in the average contacts of the susceptible class. Hence, since susceptible individuals are the ones who catch the disease and continue the spread of the infection, we only consider a cautious sub-class for the susceptible group.

The disease-free equilibrium for system (4.1) is of the form

Defining the vector of constant controls $u = (\alpha, 1-\theta, d_{1}, d_{2}, m)$ and using the next-generation matrix, we obtain the following expression for the effective reproduction number $\mathcal{R}_{e}$:

where

is the contribution of asymptomatic infectious individuals $A$ to the production of new infections, and

is the contribution of the symptomatic infectious individuals $I$ to the incidence. Please note that $\tilde{T}_{j}(0) = T_{j}$ ($j = \left\lbrace A, I\right\rbrace$), therefore, in the absence of controls, the effective reproduction number is equal to the basic reproduction number, $\mathcal{R}_{e}(0) = \mathcal{R}_{0}$.

4.1.   The impact of social distancing and home quarantine

Here, we investigate the impact of cautious behavior of susceptible individuals, which results in self-imposed prevention measures such as social distancing and home quarantine. In mathematical terms, we explore how our model dynamics depends on the control parameters $\alpha$ and $\theta$. The rest of the control parameters are not considered here.

We assume susceptible individuals in the susceptible cautious class reduce their probability of infection by a factor $1-\theta$ for different values of the parameter $\theta$, and we set the rate of behavioral change as $\alpha = 0.01$. It should be pointed out that the values of these parameters used in the simulations are theoretical as they were chosen with the purpose of highlight the possible impact of the control measures proposed in this study. This will allow us to focus on the investigation of the role played by the initial times for the application of the intervention to flatten the prevalence curve. In particular, we explore three possible initial times for the application of the intervention: (ⅰ) since day one ($t = 1$), (ⅱ) since two weeks after the first confirmed cases ($t = 15$), and (ⅲ) since a month after the first cases ($t = 30$). For the sake of simplicity, we assume that after the initial time of application the intervention is maintained for the whole time horizon. Since some of the posterior distributions have heavy tails (see Appendix A), for all the subsequent numerical simulations, the non-control parameters are fixed using the posterior mean estimates.

The results are shown in Figure 4. Dashed lines represent hypothetical health-care system capacity. Please observe that social distancing and home quarantine as control measures have the potential to reduce the maximum number of infected individuals at the peak of the outbreak and also delay the time of peak. Hence, this intervention has the potential to flatten the epidemic curve. From Figure 4(a), one can also notice that a delayed introduction of control measures increases a lot the size of the peak. For our parameters, a fifteen days delay in the use of the control causes, roughly, the number of cases at the peak to double. From Figure 4(b), it can be observed that small variations in the parameter $\theta$ have a big effect on the epidemic curve. Hence, increasing the effectiveness of social distancing and home quarantine will produce a huge benefit to reduce the prevalence of the infection.

4.2.   The impact of isolation of infected individuals

In this section, we analyze the effect of isolation of infected individuals. Therefore, we focus on the screening/diagnosed rates $d_{1}$ and $d_{2}$. In Mexico, according to some early reports, only $10$ percent of mild suspected cases are tested for COVID-19 ^[18]. On the other hand, for severe cases, $100$ percent of patients are tested. Hence, we assume $10 d_{1} = d_{2}$. In particular, we take $d_{2} = 0.2$, and $d_{1} = 0.02$ and explore how the initial time of control application influence the possible prevalence of the infection. For our simulations, the parameters which are already described by the model without control (2.1) are fixed with their posterior mean estimate. In addition, we take $\gamma_{D} = 0.1162$ from ^[8].

The results are presented in Figure 5. We can see that the diagnosis and isolation of infected individuals will reduce the maximum number of infected individuals in comparison with the no control case. It is noteworthy to mention (see Figure 5(a)) that for this intervention, the initial application of the control does not influence significantly neither the size of the peak number of diagnoses nor the peak's timing. In Figure 5(a), it can be observed that the number of individuals at the peak is less for the intervention that start at $t = 30$. However, this result is due to the fact that, for this intervention. there are less diagnoses than for the other interventions, but the number of individuals in the class $I$ is actually the same for the three interventions. Moreover, the value of the control parameters has the potential to reduce the size of the peak number of diagnoses and delay its occurrence (Figure 5(b)).

4.3.   The impact of environmental cleaning and disinfection

It has been documented that SARS-CoV-2 can be deposited onto everyday surfaces in a household or hospital setting by an infected person through coughing or touching objects and that the virus is transmissible through relatively casual contact with contaminated surfaces ^[11]. Hence, we analyze the effect of environmental cleaning and disinfection as a measure to prevent COVID-19 spread. Hence, we study the effect of the parameter $m$ related to an increase in the death rate of the virus that remains in contaminated surfaces.

The simulations in Figure 6 show that analogously to the case of the isolation measure, for this intervention the initial application of the control does not influence the size of peak and only moves the peak's time (see Figure 6(a)). An increase in the value of $m$ has the potential to flatten the epidemic curve, however, big increments in $m$ are needed to reduce substantially the prevalence of the infection below the theoretical health-care system capacity (see Figure 6(b)). Therefore, under this strategy alone will be difficult to successfully prevent further spread of COVID-19.

4.4.   Combination of control strategies

In this section, we investigate the extent of the impact of using our three control interventions simultaneously. As China and South Korea have demonstrated ^[19], social distancing is an effective measure to slow the spread of the virus and limit how many people are infected at one time. However, there is a lot of uncertainty about how long social distancing will have to last to reduce the spread of COVID-19 to near zero. Therefore, we focus on exploring the effect of different quarantine's duration on the reduction of the prevalence.

In particular, in the numerical simulations, three possible social distancing and home quarantine's duration are analyzed: one month, two months and three months. The illustrative simulations of these scenarios are presented in Figure 7(a). The values of the control parameters during active social distancing are $\alpha = 0.01$, $\theta = 0.4$, $d_{1} = 0.015$, $d_{2} = 0.15$, $m = 5$. After the application of home quarantine and social distancing the parameters $d_1$, $d_2$, $m$ keep the same value but the parameters $\alpha$ and $\theta$ are turned off. The simulations (see Figure 7(a)) show the unexpected result that extending quarantine duration does not reduce the size of the peak number of diagnoses and only moves the peak's timing. Considering the results in Figure 4, one can deduce that the most influencing factor for the efficacy of social distancing, home quarantine, and other lockdown measures is the timing. Hence, these measures must be put in use as soon as possible by health authorities.

Another important scenario that is of interest for public health officials is if there will be a need for several rounds of social distancing and home quarantine. As an illustrative example, we simulate the periodic application of the control interventions for four days, one week, two weeks and a month. That is, for example, the interventions are used for one month, then turned off for the next month and then turned on for the next month periodically. When the control is on, the values for the parameters are $\alpha = 0.01$, $\theta = 0.3$, $d_{1} = 0.02$, $d_{2} = 0.2$, $m = 5$, the results are shown in Figure 7(b). Please observe that this brings oscillations in the prevalence of the infection. The oscillations increase their altitude and amplitude with a larger time frame for the control interventions. The results in the simulations imply that, for longer periods of intermittent home quarantine, the number of infected individuals at the peak is much higher. This outcome coincides with the results obtained in ^[20].

5.   Discussion

In this study, we have proposed a compartmental epidemic model to model the transmission dynamics of the COVID-19 epidemic. Our model formulation is based on the SEIR structure augmented with appropriate compartments to take into account the current intervention measures against the spread of SARS-CoV-2. Moreover, in addition to human-to-human transmission, our model considers indirect infections caused by contact with contaminated surfaces using an extra compartment for the free-living coronavirus in the environment. We used a Bayesian approach and officially published data to calibrate the model and estimate the basic reproduction number $\mathcal{R}_{0}$.

The results of our Bayesian inference show that the value of the basic reproduction number for the MAP estimate is $\mathcal{R}_{0}^{MAP} = 2.5$, and for the posterior mean estimate is $\mathcal{R}_{0}^{CM} = 2.7$. Moreover, under a non-intervention scenario, the model outcome shows that the maximum number of infected individuals at the peak of the outbreak will be very high and can produce a collapse in the health care system. Therefore the importance of prompt implementation of effective interventions to prevent the further spread of COVID-19.

After our model calibration, we incorporated some of the current control interventions against COVID-19 into our model: (ⅰ) social distancing and home quarantine, (ⅱ) isolation of infected individuals, and (ⅲ) environmental cleaning and disinfection. We present illustrative numerical simulations as a tool to evaluate the theoretical impact of our control interventions for plausible scenarios related to the effectiveness and duration of the control application. In particular, we first study the effect of each of our interventions alone and the role played by the initial times of the application of the control to flatten the epidemic curve.

The results of our numerical simulations suggest that social distancing and home quarantine as control measures have the potential to reduce the amplitude and delay the appearance of the peak of maximum number of infected individual of the outbreak. Hence, this intervention alone has the potential to flatten the epidemic curve. However, this intervention should be implemented as soon as possible because a delayed introduction increases a lot the size of the peak of the infected. In particular, a fifteen days delay in the use of this intervention causes, roughly, the number of cases at the peak to double. The simulations also show that the diagnosis and isolation of infected individuals will reduce the size of the peak number of diagnoses and delay its occurrence. Nevertheless, to successfully control the infection more effort is needed under this intervention in comparison with social distancing and quarantine. The impact of environmental cleaning and disinfection to reduce the prevalence is low, so this strategy alone will be very difficult to achieve disease eradication. We also investigated scenarios related to an intermittent administration of quarantine and social distancing, including weekly and monthly intermittent interventions. This analysis indicates that for longer periods of intermittent quarantine and social distancing, the maximum number of infected individuals at the peak of the outbreak is higher. Hence, short periods of intermittent quarantine are more likely to prevent the collapse of the health care system. These results may be relevant when suggesting exit strategies after a long time quarantine.

Comparing the three strategies presented here, we observe that social distancing and quick isolation of infected individuals are better strategies. Environmental cleaning can also be relevant, but its application level and efficacy required to bring the maximum of the outbreak under control indicate that it might be too high to achieve in real life scenarios. Although other control measures are more effective, the high percentage in $\mathcal{R}_0$ coming from the environment-to-human transmission route makes us believe that the number of infections that can occur through this route should not be underestimated. We highly recommend a periodic home-environmental cleaning. It is noteworthy to mention that the initial application of the control does not influence the maximum number of infected individuals at the peak of the outbreak for the isolation and environmental cleaning strategies, so the winning strategy, besides being applied as soon as possible, seems to be social distancing and home quarantine. In separate work, we compare the effect of similar percentage changes for each of the three parameters and their relative effect on the outcome of the epidemics. We expect social distancing and quarantine to also be the best strategies.

Acknowledgments

Fernando Saldaña acknowledge the support by The National Council of Science and Technology, CONACyT, Mexico under grant number 331194. Ariel Camacho acknowledge SEP-PRODEP-UABC for Post-doctoral Fellowship support. Hugo Flores Arguedas acknowledge the support by The National Council of Science and Technology, CONACyT, Mexico under grant number 334204. The authors would like to thank three anonymous reviewers for improving this paper with their comments and suggestions.

Conflict of interest

The authors declare no conflicts of interest in this paper.

Source code

The source code used for this work may be found in: https://github.com/arielcam27/Saldana-et-al_2020_COVID.

Appendix.   Bayesian Inference

For the Bayesian inference, we assume the following model for the data $y_i$

where $C(t_i; {\bf{x}})$ denote the cumulative cases built from the solution $I$ of the ODE's system at time $t_i$ and ${\bf{x}} = (\beta_{A}, \beta_{I}, \beta_{V}, c_1, c_2)$ is the vector of parameters to estimate. We assume independence in the realizations of ${\bf{x}}$ and $\eta$ and $\eta_i \sim \mathcal{N}(0, \sigma^2)$. Our cumulative infected cases satisfies $C(0; {\bf{x}}) = y_0$ and

where $I(t_k; {\bf{x}})$ denote the infected cases at time $t_k$ given by the solution for the $I$ class in our model given the parameter ${\bf{x}}$. Moreover, we define by $\pi_0({\bf{x}})$ the prior distribution for ${\bf{x}}$. We assume independence of the parameters, hence

Due to the uncertainty in the values of the vector of parameters ${\bf{x}}$, we perform a sensitivity analysis of the quantity $C(t_i; {\bf{x}})$ with respect to ${\bf{x}}$ for several time values $t_i$. In particular, we use Sobol method which is a variance-based global sensitivity analysis that decomposes the variance of the output of the model into fractions which can be attributed to sets of inputs using a sensitivity index We perform the Sobol sensitivity analysis (see Figure 8) on the python library Salib ^[21] and found that the results do not depend on the times $t_{i}$. The simulations show that the Sobol indices of the parameters $c_i$ ($i = 1, 2$) are null. From our modeling, we expect $c_1 < c_2$. Actually, we expect a low value for $c_1$ since corresponds to respiratory droplets of asymptomatic individuals. We propose a Gamma distribution for $c_1$ with mean close to 0. For the transmission parameters $\beta_{j}$ ($j = A, I, V$), the lowest first order Sobol index is the index of $\beta_A$. Since we are using mass action incidence, we can expect to have values of $\beta_{j}$ ($j = A, I, V$) around the negative order of the initial susceptible population. Hence, we propose a Gamma distribution for $\beta_A$ with mean $10^{-8}$. Recall that the gamma distribution is denoted by $\Gamma (\alpha, \beta)$ with $\alpha$ the shape parameter and $\beta$ the inverse scale parameter. If $Z \sim \Gamma (\alpha, \beta)$ then $\mathbb{E}[Z] = \alpha / \beta$ and $\mathbb{V}ar[Z] = \alpha /{\beta^2}$. We propose

where $U(a, b)$ denote the uniform distribution in the interval $(a, b)$. We run a MCMC using twalk for 2000000 samples with 1000000 of burnin. The posterior distribution for each parameter are shown in Figure 9.