A data driven analysis and forecast of an SEIARD epidemic model for COVID-19 in Mexico

1.   Introduction

The COVID-19 pandemic originated in Wuhan, China in December 2019. Since then, the number of cases has accelerated in China and subsequently all over the world. The causative agent is a new betacoronavirus related to the Middle East Respiratory Syndrome virus (MERS-CoV) and the Severe Acute Respiratory Syndrome virus (SARS-CoV).

On January 30, 2020, the World Health Organization (WHO) formally declared the outbreak of novel coronavirus a Global Public Health Emergency of International Concern. In ^[1], Cruz-Pacheco et al. estimated the arrival of the infectious outbreak to Mexico between March 20 and March 30, 2020. Other models for predicting the evolution of COVID-19 outbreak in Mexico have been proposed in ^[2,3,4,5].

Compartmental models have been used for studying the spread of the COVID-19 pandemic in several countries, such as China ^[6,7,8,9], Italy ^[8,9], India ^[10] and Brazil ^[11]. Since it is known that individuals infected with SARS-CoV-2 can spread the virus to other people even when they do not present any symptoms of the disease, many authors have developed models that consider compartments for both the symptomatic and asymptomatic infectious individuals, see, e.g., ^{[5,12,13,14,15,16]}. Moreover, the rate at which the virus is transmitted to the susceptible population by contact with an infected individual is different for the symptomatic and asymptomatic, due to the varying infectivity at different stages of the infectious period and the fact that persons with COVID-19 symptoms take preventive measures like staying in quarantine, which reduces the effective number of contacts. Hence, it is more realistic to consider different transmission rates for the population in the symptomatic and asymptomatic compartments.

2.   Materials and methods

In this paper, we propose a novel SEIARD compartmental model to investigate the current outbreak of coronavirus disease in Mexico. Our model includes compartments for the exposed subpopulation and for the symptomatic and asymptomatic individuals. We conduct a detailed mathematical analysis of this model and demonstrate its application using the publicly reported data for the outbreak.

2.1.   Data

We obtained the daily cumulative number of new infected cases for 2019-nCoV epidemics in Mexico provided by the Ministry of Health of said country, but we downloaded the data of the outbreak directly from the repository of Johns Hopkins University ^[17], which is updated daily at 7 pm. The data were collected in the period since the first case of COVID-19 in Mexico (February 28, 2020) to August 12, 2020.

2.2.   Model

There are many mathematical models used to explain the spread of infectious diseases. In our case, we will use a compartmental differential equation model for the spread of COVID-19 in Mexico ^[18,19]. Our model monitors the dynamics of six subpopulations, which are: susceptible ($S(t)$), exposed ($E(t)$), infected ($I(t)$), asymptomatic ($A(t)$), recovered ($R_I(t)$ and $R_A(t)$) and dead ($D(t)$).

The model simulations were carried out with the following assumptions:

(a) The susceptible and infected individuals are homogeneous in the population.

(b) At first, no interventions were applied to stop the spread of COVID-19.

(c) No births are allowed in the population, and we only take into account the fatalities associated to COVID-19.

A diagram of the flow through the compartmental subpopulations can be seen in Figure 1.

2.3.   Basic reproduction number

This number is of great importance in epidemiological modeling, because it represents the number of secondary infections derived from a single infected individual in a fully susceptible population. In other words, it determines the potential of the spread of a virus without any control measures. The basic reproduction number was obtained from the reported cases (infected, recovered and deaths) from the period of having roughly 12 infected individuals (March 12) until August 12, 2020. The estimation of the basic reproductive number from the SEIARD model was calculated using the next-generation matrix method ^[20].

2.4.   Forecast

We calibrated our mathematical model with the daily cumulative data for new infections, recovered individuals and deaths from COVID-19 in Mexico, obtained from the open source repository of Johns Hopkins University ^[17]. Even though the outbreak of the disease started in December 2019, the first reported case in Mexico was until February 28, 2020 (an imported case). We considered the data for the outbreak of COVID-19 in Mexico up to August 12, 2020.

The set of differential equations was solved using Matlab 2016b with the ode45 solver, which is based on an explicit Runge-Kutta (4, 5) formula.

The optimization of parameters to describe the outbreak of COVID-19 in Mexico was fitted in two steps: The first step was minimizing the Sum of Squared Errors (SSE), which is calculated as follows. For a given vector of parameters $\mathbf{x}$, we compute the numerical solutions for our model $(S(t), E(t), I(t), A(t), R_I(t), R_A(t), D(t))$ and the cumulative number of infected cases with symptoms $C(t) : = I(t) + R_I(t) + D(t)$. Then the Sum of Squared Errors is given by

where $C^{\text{exp}}_i$, $D^{\text{exp}}_i$ and $R^{\text{exp}}_i$ denote the experimental data for cumulative infections, deaths and recoveries, respectively, reported for day $t_i$ ($i = 1, \ldots, n$), while $k_1$, $k_2$ and $k_3$ are coefficients used to compensate the order of magnitude for the data. In our simulations, we used $k_1 = 20$, $k_2 = 10$ and $k_3 = 1$. We considered fixed values for the parameters $w$ and $p$, assuming a latent period of 4 days ^[21] and that the total number of infectious cases (symptomatic and asymptomatic) is around eight times larger than the reported cases ^[22]. The other parameter values were fitted by applying three searches to minimize the SSE function: A gradient-based method, a gradient-free algorithm, and finally, a gradient-based method. This method was necessary to obtain the global minimum. For this step, we adapted the code from ^[9] for our mathematical model.

For the second step, we used the set of parameters obtained by minimizing the SSE as the initial value of a Markov Chain Monte Carlo (MCMC) approach, where we set the iteration number to 8000 and the first 6000 iterations as burn-in periods. Then, we computed the solutions of the model for each set of parameters obtained after the burn-in period and computed the mean and standard deviation for these solutions. Finally, we used these results to plot the 95% confidence intervals for the time evolution of the solutions. The code is available at https://github.com/UgoAvila/COVID-19.Mexico.

3.   Mathematical analysis

Susceptible population $S(t)$: This subpopulation cannot increase, because we are not considering the natural recruitment or deaths that are not related to COVID-19. The susceptible population will only decrease after an infection, an acquired characteristic due to the interaction with an infected person or asymptomatic one. The transmission coefficients will be $\beta_II$ (for contact with symptomatic infectious individuals) and $\beta_AA$ (for contact with asymptomatic infectious individuals). The rate of change of the susceptible population is expressed in the following equation:

Exposed population $E(t)$: This subpopulation consists of individuals that are infected but cannot infect others. The population decreases at a rate $w$ to become infectious with or without symptoms. We remark that $1/w$ represents the average length of the latent period, which is the period between the infection with the virus and the moment when the host is able to transmit the pathogen to susceptible individuals. An exposed individual becomes symptomatically infectious with a probability $p$ or asymptomatically infectious with a probability $1-p$. Consequently,

Infected population $I(t)$: Infected (symptomatic) individuals are generated at a proportion $p$ from the exposed class. They recover at a rate $\gamma$ and die at a rate $\delta$. This is the only population that acknowledges death. Thus,

Asymptomatic population $A(t)$: This population is considered an infected population, but the individuals do not develop the common symptoms of COVID-19. Asymptomatic individuals are important to model because they have the ability to spread the virus without knowing; they are produced at a rate $1-p$ and recover at a rate $\gamma$. Consequently,

Recovered populations $R_I(t)$ and $R_A(t)$: All individuals infected with symptoms or not will recover at a rate $\gamma$. We subdivide the recovered population in two compartments: individuals who recover after having symptoms ($R_I$) and individuals who recover from asymptomatic infection ($R_A$). Hence

Dead population $D(t)$: Infected individuals with symptoms die at a rate $\delta$, that is,

Hence, the system of differential equations that will model the dynamics of coronavirus spread in Mexico is:

We also observe that $N: = S+E+I+A+R_I+R_A+D$ is constant, where $N$ is the size of the population modeled.

3.1.   Basic reproduction number

System (3.7) has a disease-free equilibrium, which is given by $S = N$, $E = I = A = R_I = R_A = D = 0$, and we will denote it by $x_0$. We calculate the basic reproduction number $R_0$ based on this steady state. We use the next-generation matrix method proposed by ^[20]. To find $R_0$, we must solve the equation $R_0 = \rho(FV^{-1})$, where $F$ and $V$ are the derivatives of the new infections matrix $\mathcal{F}$ and the transition matrix $\mathcal{V}$, respectively, evaluated at the disease-free equilibrium. Then

The derivative of $\mathcal{F}$ at $x_0$ is:

The transition matrix is

The derivative of $\mathcal{V}$ at $x_0$ is

The inverse of $V$ is

Then

We need to find the eigenvalues of $FV^{-1}$, which are $\lambda_1 = \frac{\beta_Ip}{\delta+\gamma} + \frac{\beta_A(1-p)}{\gamma}$, $\lambda_2 = 0$ and $\lambda_3 = 0$. Then, the basic reproduction number is given by the dominant eigenvalue, that is,

3.2.   Time-dependent parameters

To describe the evolution of the epidemic in Mexico taking into account the social distancing measures taken by the government, we will assume that the infection rate, recovery rate and death rate are time-dependent functions, similar to those used by ^[9].

To model the effect of epidemic control measures, which cause the number of contacts per person per unit time to decrease as the epidemic progresses, we describe the infection rates by the functions

where $\beta_{0I}+\beta_{1I}$ is the initial infection rate by contact with the symptomatically infectious (respectively, $\beta_{0A}+\beta_{1A}$, by contact with the asymptomatically infectious); this rate decreases exponentially to the value $\beta_{1I}$ (respectively, $\beta_{1A}$) with a characteristic time of decrease $\tau_\beta$.

The time of recovery for patients may also vary with time due to the medical staff improving their therapeutic procedures. Hence, we will assume that the recovery rate is modeled by the function

where $\gamma_0$ is the recovery rate at time zero, and $\gamma_0 + \gamma_1$ is the recovery rate at a later time, which is reached after $\tau_\gamma$ days of adaptation.

Lastly, the death rate may decrease with time due to the adaptation of the pathogen or the development of more advanced treatments. Hence, we can model this with the function

where $\delta_0+\delta_1$ is the initial death rate, which decreases to the value $\delta_1$ with a characteristic time $\tau_\delta$.

If we replace the constant parameters $\beta_I$, $\beta_A$, $\delta$ and $\gamma$ in Eq (3.8) with the aforementioned time-dependent functions, we can define

as the effective daily reproduction ratio, which measures the number of new infections produced by a single infected individual per day, taking into account the evolving public health interventions and available resources ^[23].

4.   Results

The parameters of the model were fitted with the experimental data provided by a daily update from the University of Johns Hopkins, which is consistent with the data handed over daily by the Mexican Ministry of Health. Adjusting the data from the period from March 12, 2020 to August 12, 2020, we simulated the daily new COVID-19 cases in Mexico until December 2020. In Figure 2 we simulated the predicted evolution of the outbreak of COVID-19 in Mexico. As we can see in Figure 2, the red dotted line represents the active cases of this disease: Mexico has succeeded in flattening the curve. The peak of the disease will be roughly between October and November of 2020. The asymptomatic infected individuals (blue dashed line) are a major concern in the Mexican population because this subpopulation presents an exponential behavior with respect to the active cases.

The values of the best fit parameters are given in Table 1. After performing the MCMC simulations, we obtained the 95% confidence intervals for the estimated parameters. Figure 3 shows the variation of the infection rate by contact with symptomatic individuals ($\beta_I(t)$) and the infection rate by contact with asymptomatic individuals ($\beta_A(t)$). Let us evaluate first the force of infection of the symptomatic subpopulation (Figure 3a), which varies over time. The first 50 days from March 12 were the most important, starting with a value around 0.6, due to the "Jornada de Sana Distancia'' it descended to a value of 0.2. Since then this rate has stayed at a value close to 0.1. Roughly speaking, the symptomatic individuals may not be the responsible for spreading the disease due to the fact that they are quarantined with mild symptomatology or hospitalized. The force of infection of the asymptomatic individuals (Figure 3b) behaves practically the same than the symptomatic subpopulation starting around 0.8. In the first 50 days, it descended to a value of 0.3 and by then it has continued to a value near 0.2. This being said, an asymptomatic individual has the capacity to infect more individuals than symptomatic individuals. It is important to mention that the transmission rate decreased at a characteristic time of $\tau_\beta\approx31.91$, which demonstrates the emphasis of the lockdown. In Mexico, it took roughly 50 to 60 days to substantially reduce the spread of the virus.

Figure 4 shows the variation of the recovery rate $\gamma(t)$ and death rate $\delta(t)$ with respect to time. The death rate is represented in Figure 4a. At the beginning of the pandemic, the death rate in Mexico was very high with a value between 0.08 and 0.12. The average time of death was roughly 10 days since the onset of symptoms. As time passed, the increase in the average time of death was for two main reasons: The addition of less severe cases and an improvement in the treatments used in the patients. Eventually, the death rate will stabilize at a long term lethality betweem 0.02 and 0.03. The recovery rate is depicted in Figure 4b: The first 10 to 15 days were the most difficult in Mexico showing a recovery rate lower than 0.1. This dynamic could be associated with a much higher death rate. At 40 days past March 12, the recovery of patients was higher compared to the beginning of the pandemic in Mexico.

Using these values for the parameters, we can calculate the effective daily reproductive ratio $R_d(t)$ for each day (Figure 5). At the beginning of March, $R_d(t)$ was very high. After the month of March passed and the "Jornada de Sana Distancia'' was implemented, the value descended rapidly 60% percent with respect to the initial days of the pandemic in Mexico. By June, 2020, the value was near 1, and it will continue very close to that value as time passes.

In Figures 6–10, we carried out the simulation with the best fit parameters. We show a comparison of the cumulative number of infections in Figure 6). The death toll (Figure 7) in Mexico caused by COVID-19 may be around 250 000 if the outbreak of the disease maintains its course modeled by our differential equations, with a lower bound of 140 000 given by the 95% confidence interval. In Figure 8 we plotted the number of active symptomatic infections, in this graph we can acknowledge the peak of the disease will be approximately between August and Novermber 2020. At the moment we are writing this article in mid-August, Mexico is in the plateau of the exponential growth of the number of active cases. Asymptomatic individuals maintain the same behavior than the active cases (Figure 9). Finally, the recovered subpopulation until December 2020 will be between 800 000 and 2 million people (Figure 10).

In our simulations, we assumed that the population of asymptomatic individuals represents roughly 88% of all possible infections. The behavior of the asymptomatic population is important because they have the capacity to spread the virus without developing any symptoms, and public health interventions should focus on these individuals. By declaring self-quarantine to all individuals, we could separate the healthy population from the asymptomatic population and prevent the spread of the virus more rapidly. This action will only gain time in the number of infections and help avoiding the saturation of hospitals. If we control the population of infected individuals, we will prevent the death of many individuals, and the hospitals will have sufficient supplies to mitigate the severity of symptoms in patients with any type of chronic degenerative disease.

5.   Discussions

This work makes a forecast of the spread of COVID-19 in Mexico starting on August 12, 2020. The epidemic of COVID-19 in Mexico is far from over, we estimate that the force of infection has already passed and it is on the plateau of the exponential growth of the disease.

The Mexican government went ahead in establishing quarantine for all the population, a preventive measure specific to the Phase 2 of an epidemic. This measure may explain why in Mexico the number of infected individuals is different than in other countries ^[6,7,8,9]. The asymptomatic individuals are of great concern in Mexico and the world. In our case, their force of infection is greater than for symptomatic individuals because the initial value of the infection kinetics is $\beta_{0A}+\beta_{1A} = 0.83$ compared with $\beta_{0I}+\beta_{1I} = 0.66$ for the best fit parameters. This being said, Mexico should focus much more their actions on discovering rapidly this subpopulation, because they are spreading the virus without knowing and this is why the virus has been a success in infecting the entire world, but by advancing the quarantine theoretically, this subpopulation is isolated, ergo reducing considerably its kinetic infection at infinite time.

On June 1, Mexico ended the "Jornada de Sana Distancia'' and developed a system very similar to a traffic light system. This new normality was developed for the reactivation of the Mexican economy. It was determined based on 4 colors (red, orange, yellow and green), where red is the highest risk level, while green is the lowest. Due to the heterogeneity of the outbreak in Mexico, the first half of June the entire country remained at the highest risk level, which meant 15 days more of isolation.

For the 90% of those infected with COVID-19 to decrease, Mexico must have been quarantined for 93 days. Instead, they only stayed for 86 days, which resulted in the orderly reopening of quarantined individuals when the infected population had only decreased by 83%, which is a reasonable value for the reopening of non-essential economies in Mexico. Our mathematical model can be useful for real-time preparedness, anticipating the number of hospital beds and artificial respirators that will be needed once we reach the peak of infections or once Mexico enters Phase 3. The number of beds could be obtained from the number of infected symptomatic individuals, because this population can die from different illnesses caused by 2019-nCoV. Future directions of this work would be extending the SEIARD mathematical model, where we could incorporate how the isolation of infected symptomatic individuals and the quarantine of asymptomatic individuals are sufficient measures to flatten the curve, as Mexico did prior to entering Phase 2.

Acknowledgments

Ugo Avila Ponce de León is a doctoral student from Programa de Doctorado en Ciencias Biológicas of the Universidad Nacional Autónoma de México (UNAM). This paper was developed in the period of his PhD studies. Ugo Avila Ponce de León also received a fellowship (CVU: 774988) from Consejo Nacional de Ciencia y Tecnología (CONACYT).

This article was supported in part by Mexican SNI under grant number 15284.

Conflict of interest

All authors declare no conflict of interest in this paper.