Analytical solution of SEIR model describing the free spread of the COVID-19 pandemic.

We analytically study the SEIR (Susceptible Exposed Infectious Removed) epidemic model. The aim is to provide simple analytical expressions for the peak and asymptotic values and their characteristic times of the populations affected by the COVID-19 pandemic.

Introduction

The COVID-19 outbreak has motivated a large number of numerical studies using epidemiology models [1], [2]. A commonly used model is the Susceptible-Exposed-Infected-Removed (SEIR) model [3]. This model is formulated as a system of nonlinear ordinary differential equations, for which no exact analytic solution has yet been found. For this reason, most of the recent works focus on the numerical analysis of statistical ensembles of initial data for these equations. However, due to the uncertainty and often unreliability of the clinical data, the prediction about the real evolution of the epidemic is rather difficult, if not impossible [4], [5]. On the other hand, the SEIR epidemic model provides a deterministic evolution for some given initial state. Therefore, the aim of this work is to provide simple expressions of the main characteristics of the population of individuals that have been in contact with the disease, as of instance the peak of the infected population and the time after which it occurs, the final number of individuals who have contracted the disease and the temporal shape of the infectious population’s curves. These analytical expressions can become useful through their application to the COVID-19, to obtain fundamental parameters as the reproduction number r and the epidemic starting time.

The paper is organized as follow, In Section 2 we recall the SEIR model; in Section 3 we study the linear regime with the exponential growing and decaying evolution, depending on the reproduction number r; in Section 4 we investigate the nonlinear regime in the free spread evolution with r > 1. We approximate the exact model of equations by a reduced model where the decaying mode is adiabatically eliminated. This reduced model allows to obtain analytical results which have been seen to be in good agreement with the exact numerical solution. Section 5 summarizes the results and draws the conclusions.

The SEIR model

We used the susceptible-exposed-infected-removed (SEIR) compartment model [3], [6], [7], [8] to characterize the early spreading of COVID-19, where each individual could be in one of the following states: susceptible (S), exposed (E, being infected but without infectiousness), infected (I, with infectiousness), recovered (R) and dead (D). At later times a susceptible individual in the state S would turn to be an individual in the exposed state E with a rate r/τ, where r is the reproduction number (i.e. the average number of infected people generated by each infected person during the desease) and is the average time in the infected state I. An exposed individual in the state E becomes infected, i.e. in the state I in an average time . Then the infected individual is removed from the total population with the rate γ 2 either by recovering (R) or dying (D) with a mean case fatality proportion p. The dynamical process of SEIR is described by the following set of equations: Here S(t), E(t), I(t), R(t) and D(t) respectively represent the number of individuals in the susceptible, exposed, infectious, recovered and death states at time t and N is the total number of individuals in the system such that . Finally, the cumulative population C isequal to the total population of individuals who have contracted the infection.

Linear regime

If E(t), I(t), R(t) ≪ N(t), then the susceptible population S can be approximated by the total population N (i.e. S ~ N) and the equations for the exposed and infected population are linear:

General solution of the linear equations

Introducing the Laplace transformswith Reλ > 0, Eqs. (7) and (8) becomeswhere E(0) and I(0) are the initial conditions. The eigenvalues λ are solution ofgivingwith solutionswhereSince Δ > 0 the eigenvalues are real. Depending on r, we distinguish three cases:

If r > 1 then so that and . The solution grows exponentially (explosive regime);

If r < 1 then so that both and . The solution decays exponentially (relaxation regime);

If then so that and . The solution remains partially constant (marginally stable regime).

For the cases (a) and (b) the solution is whereas in the case (c) () the solution is

Analysis

The only parameter which can be controlled by confinement measures is the reproduction number r. In the following we assume that for COVID-19 the characteristic times are days and days [9]. We consider the time evolution of the population E and I for r > 1, and r < 1, corresponding to the explosive, marginally stable and relaxation regimes, respectively.

Explosive regime

For r > 1 and with .

Marginally stable regime

When in the asymptotic limit t ≫ τ and I are constant,and the death population grows linearly in timewhere D(0), E(0) and I(0) are the values taken at time when r starts to be .

Relaxation regime

When r < 1, is negative and E and I tend to zero, whereas D tends to the following constant value,where D(0), E(0) and I(0) are the values taken at time when r starts to be r < 1. Fig. 1 shows a typical temporal evolution of I(t) and D(t) starting with r > 1, then subsequently changed to and later on to a value r < 1. The regime is linear (i.e. with E, I ≪ N), the initial values are and and . The red dashed line is for (explosive regime). The green dashed-dotted line is for r changed from to at (marginally stable regime) and the blue solid line is for until then between and and finally for t > 30 (relaxation regime). Notice the asymmetry of the curve of I(t) due to the different growing and decaying rates.

Fig. 1

Evolution of I(t) and D(t) with (red dashed line), after days (green dashed-dotted line) and after days (blue solid line). Initial conditions: . ; .

Nonlinear regime

In the following, we investigate the nonlinear regime with a constant reproduction number r > 1. This corresponds to a free spread of the infection, with an initial exponential growth of the exposed population E, and so also of I and D. The exponential growth stops when susceptible population S becomes sensibly less then the total number N of the individuals. This regime is similar to the saturation in a single-mode laser, where steady-state is reached when the gain of emitted photons equals the losses by the cavity [10]. Notice thatis a constant of motion and . However, if p ≪ 1 we always have D ≪ N 0, so that with a good approximation we can approximate N by N 0. Introducing the removed population we can eliminate using the constant of motion and obtain We normalize the variables by N 0 defining and where is the cumulative population, i.e. the total number of individuals who have contracted the infection. Then the equations become

These equations have a single steady-state solution (i.e. ) with (end of the epidemics) and with 0 < z 0 < 1. This solution is stable if . We see that the stability condition impliesIn Fig. 2 we plot E/N, I/N and C/N for days, days and initial conditions . We observe that C/N tends to a steady-state value of about 0.6, whereas the peak of I/N is about 0.03: it means that for these parameters the 60% of the total population has contracted the infection and the peak the infected population is about 3% of the total population. Note that these results are independent on p and depend only on τ and r.

Fig. 2

Simulation with and . Initial conditions: . I/N (a) and C/N (b) vs. time from the numerical solution (solid black line) and from the analytic expressions, Eqs. (63) and (62) (dashed blue line). The time t is in units of days and days, days.

Reduced model

In this section we find an approximated analytic solution of Eqs. (27)–(29) in the free spread evolution with r > 1. The idea is to adiabatically eliminate the decaying mode with negative eigenvalue . To this aim, it is convenient to write Eqs. (27)–(29) in the basis of the eigenvalues λ ±. Writing again the linear Eqs. (7) and (8) in the formthe normalized eigenvectors associated to the eigenvalues λ ± of Eq. (12) arewhereHence, in the new basisand the inverse isIn the new basis Eqs. (27)–(29) take the form: Notice that as expected in the linear regime the dynamics of and are uncoupled. Now we consider the free spread regime with r > 1 such that is positive and is negative. If is small, then and we can adiabatically eliminate the ’slave’ variable . Neglecting in (37) we obtainwhich when inserted in Eqs. (36) and (38) yields Since and where and . Finally, the original variables are

Analytical solution

Eqs. (42) and (43) may provide some analytical result. Rescaling the time asand defining Eqs. (42) and (43) take the form: In the limit β → 0 they have the form of Lotka-Volterra equations [11]. From them, dividing member by member, it resultswhich when integrated yieldswhere we assumed s → 0 when z → 0. On the other hand, s → 0 when z → z ∞ (see Fig. 3 ), where z ∞ is the solution of the transcendental equationThe same transcendental Eq. (53) for z ∞ has been obtained for the SIR compartmental model [12], [13]. Here we have demonstrated its validity also for the SEIR model.

Fig. 3

Plot of s vs. z for from Eq. (52).

We see from Fig. 3 that for and . The maximum value of s occurs when so thatThese simple equations provide two analytic expressions for the asymptotic value of C/N and for the peak of I/N.

Let’s now find an approximated solution of z as a function of the scaled time τ. Using Eq. (52) in Eq. (50) we obtain a differential equation for z:From the numerical analysis and assuming βz ≪ 1, we find that z(τ) is well approximated by the following function:where τ depends on the initial conditions. From (49) it follows for βz ≪ 1This equation can be integrated to giveSince kτ ≫ 1 and, from Eqs. (52) and (56), we can write Eq. (58) in the following form:The time at which s(τ) is maximum can be evaluated from the condition which, using Eq. (56), yieldswhere . For instance, for and we obtain and .

Results and conclusions

We have obtained analytical expressions for the asymptotic value of the cumulative population fraction C/N and the peak of the infectious population fraction I/N in the case of free spread evolution of COVID-19. Furthermore, we have obtained approximated expressions of these quantities as a function of time and the times at which the peak and the end of the epidemics is expected. We summarize here below these results:

The asymptotic value of the cumulative population fraction is where z ∞ is the solution of the transcendental Eq. (53). A comparison between the exact solution obtained by integrating Eqs. (1)–(5) and the solution of Eq. (53) is shown in Fig. 4 (a). Notice that this value depends only on the reproduction number r.

Fig. 4

(a): Plot of C∞/N vs. r, from the numerical solution of Eqs. (1)-(5) (dashed line) and from the analytical result of Eq. (53) (continuous line). The dotted red line is the threshold value . (b) Peak value of I/N vs. r from the numerical solution of Eqs. (1)-(5) (dashed line) and from Eq. (61) (continuous line).

The peak value of the infectious population fraction is, from Eqs. (45), (48) and (54),The agreement of this expression with the exact result shown in Fig. 4(b) is better for values of r closer to the threshold .

We have obtained an approximated temporal profile of C(t)/N,where and where x 0 and y 0 are the initial values of x and y. From this expression we have obtained the expression of I(t)/N as a function of time:where andThe good agreement of Eqs. (62) and (63) with the exact numerical solution of Eqs. (1)–(5) is shown in Fig. 2.

The time at which the peak of I/N is reached is Fig. 5 shows t peak (in units of days) as a function of r for an initial value of and .

Fig. 5

Plot of peak time tpeak (in units of days) vs. r for initial values of and and days, days.

These analytic expressions can be useful for deriving the uncertainty in the estimates of COVID-19 caused by the fluctuations of the values of the control parameters, as for instance the reproduction number r. In fact, the results of Ref. [4] suggest that uncertainties in both parameters and initial conditions rapidly propagate in the model and can result in different outcomes of the epidemics. For instance, Fig. 4a and b show the dependence of the fraction of the final cumulative fraction, C ∞/N, and the daily infections peak, I peak/N, as a function of r. We observe that the sensitivity of C ∞/N on r variations is larger when r is close to unity (with approximately ) whereas it decreases for increasing values of r. On the other hand, I peak/N grows almost linearly with r (approximately as ), so that its sensitivity to r variations is almost constant. Finally, the uncertainty of the peak time t peak (see Fig. 5) on r variations is very large for r close to unity and it reduces strongly at larger r.

In conclusions, we have obtained analytical expressions for the peak and asymptotic values of COVID-19 pandemic curves in the free spread as a function of the reproduction number and the two average times in the exposed and infected states. The results have been obtained by reducing the exact nonlinear model by adiabatically eliminating the decaying mode of the linear regime. This allows to reduce the SEIR model of a set of two equations similar to the Lotka-Volterra equations, from which exact and approximated solutions can be obtained. The analytical results have been compared with the exact numerical solution, showing good agreement. Particular interesting is the asymptotic fraction of the removed (recovered+deaths) population fraction, which depends only on the reproduction number r. Finally, the infected population curve is an almost symmetric function described by an hyperbolic secant function.

CRediT authorship contribution statement

Nicola Piovella: Conceptualization, Formal analysis, Writing - original draft.

Declaration of Competing Interest

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