Mathematical Models for Describing and Predicting the COVID-19 Pandemic Crisis

A simple mathematical model for disease epidemic can been built dividing the population in 3 groups: susceptible individuals (S), infected individuals (I) and recovered individuals (R). The model composed by these groups is called the SIR model. In this article, however, we consider also individuals who have died by the disease, denote by D. Following the same arguments of the SIR model, the SIRD model can be described by the set of four differential equations:

		$\displaystyle\frac{dS}{dt}=-\frac{\beta}{N}I(t)S(t)$		(1)
		$\displaystyle\frac{dI}{dt}=\frac{\beta}{N}I(t)S(t)-(\gamma+\mu)I(t)$		(2)
		$\displaystyle\frac{dR}{dt}=\gamma I(t)$		(3)
		$\displaystyle\frac{dD}{dt}=\mu I(t)$		(4)

Last equation is easily understood by thinking that the variation of the number of deaths may be proportional to the infected individuals, where the proportionality constant is denoted by $\mu$. The constants $\gamma$ and $\beta$ are, respectively, the recovery rate and the number of infected, where $\mu$ and $\gamma$ are given in terms of the infection fatality rate (IFR) or the case fatality rate (CFR); that is, the number of people who contracted the disease and died according to the total number of infections (IFR) or the registered number of infections (CFR), and the average time taken from symptoms onset to recovery, $\tau_{r}$, or death $\tau_{d}$, formally $\mu=P_{CFR}/\tau_{d}$ and $\gamma=(1-P_{CFR})/\tau_{r}$. The equations are simply a mathematical way to describe how individuals passes from one group to the other according to the following chain of events: A susceptible individual becomes infected by the virus, and from this point, it either dies or recovers (Figure 1).

Summing the four equations we get

$\displaystyle S(t)+I(t)+R(t)+D(t)=\text{const},$

(5)

where the constant may represent the total number of individuals, $N$. Before proceeding, we propose the initial condition that, when $t$ goes to zero, $I(t)=I_{0}$, $R(t)=D(t)=0$, and, therefore, $S(t)=S_{0}=N-I_{0}\approx N$. Such an assumption is based on the fact that the entire population is susceptible to the SARS-CoV-2 virus.

Since $R(t)$ and $D(t)$ are both data updated day by day in Germany and Korea, it would be helpful to write $I(t)$ as function of them so as to predict its behavior, obtaining, for example, the maximum number of infected individuals. For reasons that may be clear soon, we first write $I(t)$ in terms of $S(t)$. An intuitive step is to divide equation (2) by (1). Thus,

$\displaystyle\frac{dI/dt}{dS/dt}=-1+k\frac{1}{S},$

(6)

where $k=(\gamma+\mu)/\beta$. Eliminating the temporal dependence, we get a separable differential equation, that is,

$\displaystyle dI=-dS+k\frac{dS}{S},$

(7)

which the solution is easily verified to be

$\displaystyle I(t)=-S(t)+k\ln{S(t)}+\text{const}.$

(8)

Applying the initial condition, we obtain

	$\displaystyle I_{0}=-(N-I_{0})+k\ln(N-I_{0})+\text{const}$		(9)
	$\displaystyle\rightarrow\text{const}=N-k\ln(N-I_{0}).$

Hence, equation (8) may be written as

$$I(t)=N-S(t)+k\ln\frac{S(t)}{N-I_{0}}.$$

(10)

We can visualize here, that depending on the combination of $\gamma$ and $\mu$, $I$ reaches 0 before the entire population $S$ becomes infected (Figure 2).

Next, we may write $S(t)$ in terms of $R(t)$ and $D(t)$. For this purpose, we begin by dividing equation (1) by (3) and (1) by (4),

	$\displaystyle\frac{dS}{dR}=-\frac{\beta}{\gamma}S(t)\,\,\,\text{and}$		(11)
	$\displaystyle\frac{dS}{dD}=-\frac{\beta}{\mu}S(t).$		(12)

Adding these two equations and writing $S(R,D)$ as $S(R,D)=f(R)g(D)$, we get

$\displaystyle\frac{1}{f}\frac{df}{dR}+\frac{1}{g}\frac{dg}{dD}=-\beta\left(\frac{1}{\gamma}+\frac{1}{\mu}\right)$

(13)

Since (13) is a separable equation, the well-known solution is given by

		$\displaystyle f(R)=Ae^{aR}\,\,\,\text{and}$		(14)
		$\displaystyle g(D)=Be^{bD}.$		(15)

Therefore, $S(t)$ can be written as

$\displaystyle S(t)=Ce^{aR(t)\,+\,bD(t)},$

(16)

where we absorbed both constants $A$ and $B$ into $C$. By the initial condition, we find that $C=N-I_{0}$. To find $a$ and $b$, we must derive (16) in time under the condition that it may return to equation (1). In this way, we see that

$\displaystyle a\gamma+b\mu=-\beta.$

(17)

By the other hand, substituting equations (14) and (15) in (13), we get

$\displaystyle a+b=-\beta\left(\frac{1}{\gamma}+\frac{1}{\mu}\right).$

(18)

Solving this system,

	$\displaystyle a=-\frac{\beta}{\gamma(1-\gamma/\mu)}$		(19)
	$\displaystyle b=-\frac{\beta}{\mu(1-\mu/\gamma)}$		(20)

Hence, $I(t)$ can be finally written as

$\displaystyle I(t)=N-(N-I_{0})e^{aR(t)+bD(t)}+\gamma\mu\frac{\gamma+\mu}{\gamma-\mu}\left(\frac{R(t)}{\gamma^{2}}-\frac{M(t)}{\mu^{2}}\right)$

(21)

With this equation, see that as $t\rightarrow\infty$, $I$ does not approaches $N$ necessarily, depending on the recovery and death rates, $I$ does not reach $N$.

The last important quantity extracted from this model is the basic reproduction number $\mathcal{R}_{0}$, given by [14]:

$\displaystyle\mathcal{R}_{0}=\frac{\beta}{\gamma+\mu}.$

(22)

This quantity, is of vital importance of the study of a disease outbreak.

Another deterministic mathematical model possible is the SEIRD model, in which we consider the population $N$ of a given region as divided in 5 groups. At time $t$, there are those who are susceptible to get infected $S(t)$, the ones who have already been exposed the virus but does not present symptoms yet $E(t)$, people who are already infected and present the symptoms $I(t)$, the ones that have already recovered from the disease $R(t)$ and those who are dead due to the infection $D(t)$. This model is a good approximation to a short epidemic, so the population of a region is roughly constant throughout the epidemic period. Also, since this is a deterministic model, we assume $N$ to be a big number compared to the number of people associated with the infection of a single person. The final consideration is that we also assume that people that are recovered from the disease acquire immensity and does not become susceptible to become infected again.

The rate of infection $\lambda$ is proportional to the number of people infected, $\lambda(t)=\beta I(t)$, where the constant $\beta$ represents the effectiveness of the infection, the rate of cure $\gamma=P_{:)}\tau_{r}^{-1}$, where $P_{:)}$ is the probability of recovery and $\tau_{r}$ is the average time taken for an infected person to recover. Similarly the rate of death is $\mu=P_{CFR}\tau_{d}^{-1}$, where $P_{CFR}=1-P_{:)}$ is the probability of death, given by the CFR and $\tau_{d}$ is the average time taken for an infected person to die. Figure 3 carries an visual representation of the SEIRD model.

The differential equations representing the evolution of the populations are given by

		$\displaystyle\frac{dS}{dt}=-\frac{(1-P_{exp})\beta}{N}I(t)S(t)-\frac{P_{exp}\beta}{N}E(t)S(t)$		(23)
		$\displaystyle\frac{dE}{dt}=\frac{(1-P_{exp})\beta}{N}I(t)S(t)+\frac{P_{exp}\beta}{N}E(t)S(t)-cE(t)$		(24)
		$\displaystyle\frac{dI}{dt}=cE(t)-\gamma I(t)-\mu I(t)$		(25)
		$\displaystyle\frac{dR}{dt}=\gamma I(t)$		(26)
		$\displaystyle\frac{dD}{dt}=\mu I(t)$		(27)

We first turn our attention to the construction of an appropriate formula for calculating $\mathcal{R}_{0}$ with this model. For that we follow the method derived on [20]. The study develops a mathematical generalization for writing $\mathcal{R}_{0}$ depending on the type of epidemiological model. $\mathcal{R}_{0}$ is defined as

$\displaystyle\mathcal{R}_{0}=\rho(FV^{-1})$

(28)

where $\rho(X)$ means the spectral radius of the matrix X, that is, the largest absolute eigenvalue. Both $F$ and $V$ are the matrices of the derivatives of the functions defining the behavior of the disease population, with respect to each population compartment.

To get to these matrices, we first note that the set of equations regarding the dynamics of the SEIRD model can be expressed as follow: Consider $\vec{x}$ the vector of populations, that is $\vec{x}=(x_{1},x_{2},x_{3},x_{4},x_{5})$ where $x_{1}=E$, $x_{2}=I$, $x_{3}=S$, $x_{4}=R$ and $x_{5}=D$. Analogously, $d\vec{x}/dt$ is the vector of the first derivatives. Then, we can write the dynamics of the populations as

$\displaystyle\frac{d\vec{x}}{dt}=\mathcal{F}-\mathcal{V}$

(29)

where $\mathcal{F}$ is the vector that relates the appearance of new infections on the disease populations due to contamination, and $\mathcal{V}$ is the input and output of members in all populations due to all other causes, such as recovery from the disease, development of symptoms after an incubation period, etc. In our case, since all newly infected members go to the $E$ population

$\displaystyle\mathcal{F}=\begin{pmatrix}(1-P_{exp})\beta IS+P_{exp}\beta ES\\ 0\\ 0\\ 0\\ 0\end{pmatrix}$

(30)

while

$\displaystyle\mathcal{V}=\begin{pmatrix}\beta cE\\ -cE+\gamma I+\mu I\\ (1-P_{exp})\beta IS+P_{exp}\beta ES\\ -\gamma I\\ -\mu I\end{pmatrix}.$

(31)

Now, we know that the situation of a disease free equilibrium (DFE), meaning no disease is happening, is achived by the vector $\vec{x}_{0}=(0,0,S_{0},0,0)$, where $S_{0}=N$. According now to [20] we can calculate $F$ and $V$ as

	$\displaystyle F=\left.\left(\frac{\partial\mathcal{F}_{i}}{\partial x_{j}}\right)\right|_{x=x_{0}}\hskip 28.45274pt1\leq i$		(32)
	$\displaystyle V=\left.\left(\frac{\partial\mathcal{V}_{i}}{\partial x_{j}}\right)\right|_{x=x_{0}}\hskip 28.45274ptj\leq m$		(33)

being $x_{j}$ the vector components of $\vec{x}$ related to the populations with the disease, in our case $E$ and $I$, and $m$ is the number of populations related to infectious beings. Here $m=2$. Performing the derivatives, we conclude

	$\displaystyle F=\begin{pmatrix}P_{exp}\beta&(1-P_{exp})\beta\\ 0&0\end{pmatrix}$		(35)
	$\displaystyle V=\begin{pmatrix}c&0\\ -c&\gamma+\mu\end{pmatrix}$		(36)

The next step is to find the inverse matrix of $V$, fortunately $V$ is a 2x2 matrix and the formula for it’s inverse is straightforward

$\displaystyle V^{-1}=\frac{1}{c(\gamma+\mu)}\begin{pmatrix}\gamma+\mu&0\\ c&c\end{pmatrix},$

(37)

and we proceed to the last step of combining $FV^{-1}$ in order to retrieve $\rho(FV^{-1})$ and find $\mathcal{R}_{0}$.

	$\displaystyle FV^{-1}=\frac{1}{c(\gamma+\mu)}\times$		(38)
	$\displaystyle\times\begin{pmatrix}P_{exp}\beta(\gamma+\mu)+(1-P_{exp})\beta c&(1-P_{exp})\beta c\\ 0&0\end{pmatrix},$

therefore, by computing the eingenvalues of $FV^{-1}$ we find

$\displaystyle\mathcal{R}_{0}=\frac{P_{exp}\beta\left[\gamma+(1-P_{exp})\beta\right]+(1-P_{exp})\beta c}{c(\gamma+\mu)}$

(39)

Having $\mathcal{R}_{0}$ in our hands, we continue to the study of some behaviors of this model.

The set of equations describing the model is subjected to the initial conditions. When $t\rightarrow 0$, $I(t)-\rightarrow I_{0}$, $S(t)\rightarrow S=N-I_{0}$, $R(t)\rightarrow 0$, $D(t)\rightarrow 0$ and $E(t)\rightarrow E_{0}$, where $I_{0}$ is the initial number of infected, $E_{0}$ is the initial number of exposed in the population and no deaths or recoveries are assumed at $t=0$.

Without vaccines or efficient medicine against the disease, non-pharmaceutical interventions are the only effective way to prevent further increase of the pandemic [13]. These interventions take different approaches such as social distancing, social isolation and lockdown of the population. Despite the differences, they all carry the same objective, decreasing the infection rate $\beta$. It is convenient to implement the effect of these interventions on the model, when making predictions. Here, we model this effect by a logistic function, where $\beta$ starts at a initial value $\beta_{i}$ and at some critical time $t_{c}$ a intervention is imposed and beta decreases to $\beta_{f}=P_{dec}\beta_{i}$, where $P_{dec}$ is the fraction of $\beta_{i}$ decreased by the intervention. In France, studies estimate that the intervention decreased $\beta_{i}$ by 77% [21], therefore $P_{dec}=0.77$ in France.

$\displaystyle\beta(t)=\frac{(1-P_{dec})\beta_{i}}{1+\tau e^{t-t_{c}}}+P_{dec}\beta_{i}$

(40)

where $\tau$ is a constant related to the time taken for the intervention to have the effect desired. Such model reconstruct the general behavior of interventions against the spread of the disease (Figure 4)

Since the case fatality rate (CFR) of COVID-19 is different among age groups [6, 7, 22], we propose here a modification on both models, including the age distribution of the population and the social aspects of close contact between members of the population. The modification is describe as follow: Each compartment is divided into $M$ age groups, where each $i$-th group has a $P_{CFR_{i}}$ associated to it, that is, the probability of death associated to the $i$-th age group. The $\beta$ parameter is now described as the average number of daily contacts between a member belonging to the $i$-th age group to the $j$-th age group, multiplied by the infection probability $P_{infc}$

	$\displaystyle\beta_{i_{I}}=\sum_{j=1}^{N}C_{ij}I_{j}P_{infc},$		(41)
	$\displaystyle\beta_{i_{E}}=\sum_{j=1}^{N}C_{ij}E_{j}P_{infc},$		(42)

where $C_{ij}$ is called the social contact matrix and we included $I_{j}$ and $E_{j}$ inside $\beta$ now to place everything on the same sum. The age distribution among the population is retrieved from the UN prospects [23] and the social contact matrix for those countries was measured on previous studies [24, 25]. The specific contact matrix for the Republic of Korea was not found, however, [26] finds evidences of cultural clusters in the world, where countries belonging to the same cluster share cultural similarities; thus, we use this fact to justify the use of Hong Kong’s social contact matrix to describe the Republic of Korea. That way, we include cultural and population aspects for each of those countries, increasing the odds of a successful prediction. This type of model was used recently to describe the coronavirus outbreak on large cities in Brazil [27].

For Germany, the fitting data corresponds to the cases and deaths until 17th March. However, until the peak is reached, both models find presents very large margin of error for $N$, to avoid this problem, we varied $N$ manually, from 0 to 10% of the local population; which was taken from a united nation prospect for the year of 2020 [23], at steps of 0.05%. At each step, we fit the initial data and reject the fitting if the $\chi^{2}$ value is lower than 0.995, this $\chi^{2}$ method for validating the goodness of a data adjustment had already been used for epidemiological models [41]. We than plotted the maximum and minimum acceptable fits to generate the margin of prediction, comparing it with the complete dataset. We also decided to use $\tau_{r}$ and $\tau_{d}$ according to clinical studies when performing the prediction, instead of leave them as free parameters for the fitting, $I_{0}$ was also determined a priori according to the first registered number on 15/02/2020. The resulting free parameters for fitting the training set are $P_{infc}$ and $E_{0}$.

With the SEIRD model, we found the maximum and minimum values of $N$ to be $N_{min}=0.25\%$ of the German population and $N_{max}=0.40\%$ of the German population. The $P_{infc}$ parameter varied from 15.5% to 16.2%

Using the SIRD, leaving $\beta$ and $I_{0}$ as free parameters. The limit values of $N$ were $N_{min}=0.15\%$ and $N_{max}=0.2\%$, while $\beta$ varied from $0.247$ to $0.279$ and $I_{0}$ from 11 to 30. Figures 8 and 8 show the result of both simulations, with the maximum $N_{max}$ and minimum $N_{min}$ curves. The shaded region is the region between $N_{max}$ and $N_{min}$.

The training set consisted of 20 days, corresponding to the infections from 15/02 to 06/03. The SEIRD model found $N_{min}=0.03\%$ and $N_{max}=0.04\%$ of the total Korean population, while $P_{infc}$ varied between 80 to 85%.

The simple SIRD model found $N_{min}=0.02\%$ and $N_{max}=0.03\%$, $\beta$ went from 0.345 to 0.436, and $I_{0}$ was between 23 to 65. Figures 10 and 10 present the result for prediction of both models.