Dynamics of SEIR epidemic model by optimal auxiliary functions method.

The aim of the present work is to establish an approximate analytical solution for the nonlinear Susceptible, Exposed, Infected, Recovered (SEIR) model applied to novel coronavirus COVID-19. The mathematical model depending of five nonlinear differential equations, is studied and approximate solutions are obtained using Optimal Auxiliary Functions Method (OAFM). Our technique ensures a fast convergence of the solutions after only one iteration. The nonstandard part of OAFM is described by the presence of so-called auxiliary functions and of the optimal convergence-control parameters. We have a great freedom to select the auxiliary functions and the number of optimal convergence-control parameters which are optimally determined. Our approach is independent of the presence of small or large parameters in the governing equations or in the initial/boundary conditions, is effective, simple and very efficient.

Introduction

After the appearance of the Severe Acute Respiratory Syndrome in 2012 (SARS-CoV) and of the Middle East Respiratory Syndrome coronavirus in 2012 (MERS-CoV) the new virus named Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), led to a considerable problems all around the world. Unfortunately, the novel coronavirus continues to proliferate in many regions of the world. The place of birth of the COVID-19 is Wuhan, Hubei Province of China where it has been reported the first case in December 2019. At this the moment there are over 71 millions confirmed cases and over 1.6 millions deaths (15 December 2020).

The origin of coronavirus is uncertain, but it is believed that the source of the epidemic is the seafood and live animals market [1]. Some investigations over the time confirmed that the coronavirus transmission occurs from person to person [2]. Mathematical models have become a crucial tool in studying the epidemiological dynamics of epidemics: mathematical models can give important informations to the public health authorities about mechanism of transmission for coronavirus, short and long term predictions, reductions of the number for infected people and deaths lockdown, minimizing the peak number of infected people, quarantine application, vaccination campaigns, different scenarios to reduce the epidemic risks and so on. Mathematical method for coronavirus has attracted the attentions of epidemiologists, pharmacists, biologists, chemists or mathematicians.

The first scientist who proposed a mathematical model for infectious diseases was Bernoulli [3]. Kermack and McKendrick proposed in [4] the mathematical epidemic model of Susceptible-Infected-Recovered (SIR) a model to describe the spread of diseases.

In the last years, a large number of mathematical models have been devoted to the outbreak and spread of diseases. Li and Cui [5] studied a SEIR model with vaccination strategy. Using La Salle's invariant set theorem and Lyapunov function has proved the global asymptotical stable results of the disease free equilibrium. It is showed that there exists a periodic solution if the system has three equilibrium points. Ibeas et al. [6] proposed a robust vaccination strategy capable of eradicating infectious diseases from a population group. For this aim, a control theoretical approach based on a sliding-mode control law is used. Harir et al. [2] applied the variational iteration method and differential transform method for solving the SEIR epidemic model. Widyaningsih et al. [7] considered the SEIR model with immigration, determining equilibrium points and stability. Dynamics of a fractional order mathematical model for COVID-19 is analyzed by Rezapour et al. [8] and by Zhang et al. [9]. The numerical simulations indicate that there is a good agreement between theoretical results and numerical ones.

The differential equations of SEIR epidemic model are solved numerical with a forward Euler scheme for different values of the reproduction ratio. The Markov Chain Monte Carlo method in simulation of SEIR model for 45 possible scenarios is used by Iwata and Miyakoshi [11]. Lin et al. [12] proposed conceptual models for COVID-19 outbreak in Wuhan with the consideration of individual behavioural reaction and governmental actions, e. g., holiday extension, travel restriction, hospitalization and quarantine. Piovella [13] presented an analytically study of SEIR model for the peak and asymptotic values and their characteristic times of the population affected by the COVID-19.

The SEIR model is improved by Godio et al. [14] by means of a stochastic solver which identifies a set of possible solutions predicting the epidemic evolution with the associated uncertainly assessment. Heng and Althaus [15] obtained an approximate semi-analytical solution of SEIR model depending on the population that is infections. Tsay et al. [16] introduced a novel optimization based decision-making framework for managing the COVID-19 outbreak in the US including modelling the dynamics of affected population and optimal control strategy for sequencing social distancing and testing events.

Kwuimy et al. [17] developed a genetic algorithm technique to identify key model parameters employing COVID-19. Stability, bifurcations, dynamic behaviour and influence of social and government behaviour on disease dynamics are analyzed. He et al. [18] introduced the seasonality and stochastic infection parameters and nonlinear dynamic including chaos are founded in the system. Also, the control strategies of the COVID-19 are based on the structure and parameters of this model. Yang and Wang [19] studied a new mathematical model of COVID-19 including multiple transmission pathways with both the environment-to-human and human-to-human routes. Pengpeng et al. [20] considered a new SEIR propagation dynamics model considering the weak transmission ability of the incubation period, the variation of the incubation period length and the government interaction measures to track and isolate comprehensively. Mwalilli et al. [21] discussed a modified SEIR model incorporating the environment and social distancing. The next generation matrix approach was used to obtain the reproduction ratio .

A well-mixed SEIR model was employed by Hou et al. [22] to describe the dynamic of the COVID-19 epidemic based on epidemiological characteristics of individuals, clinical progression and quarantine intervention measure of the authority. It is showed that by reducing the contact rate of latent individuals, intervention such as quarantine and isolation can reduce the potential peak number of infections. Annas et al. [23] analyzed the stability of SEIR model by considering vaccination and isolation factors as model parameters, to obtain the basic reproduction number.

The purpose of this research is to obtain explicit and accurate analytical approximate solutions for SEIR model, taking into consideration the death toll. The system of five nonlinear differential equations with five unknowns is numerically examined in [10] but is very hard to solve by means of classical analytical procedures, without the presence of small parameters.

In comparisons with other known techniques applied to find approximate solutions for nonlinear dynamical systems, our approach is based upon original construction of the solution using a moderate number of the convergence-control parameters which are components of the so-called optimal auxiliary functions. We remark that these parameters lead to a high precision by comparing our approximate solutions with numerical solutions. Any nonlinear differential equation is reduced only to two differential equations. New explicit and accurate analytical approximate solutions were finally obtained by using only one iteration, this higlights the power and efficiency of the proposed technique.

Governing equations

The total population is divided into the following four categories: Susceptible (S), Exposed (E), Infected (I) and Recovered (R), when no vaccine is available. The governing differential nonlinear equations are [7,10,14,19] and [23]:

The dead population is:where , the dot denotes differentiation with respect to time, and the parameters are defined as: (per-capita birth rate), (per-capita natural death rate), (virus induced average fatality rate), (probability of disease transmission per contact-dimensionless), (rate of progression from exposed to infections) and (recovery rate of infectious individuals).

The initial conditions for Eqs. (1)–(5) are:

Basics of the OAFM

We consider a general nonlinear differential equationwith the initial/boundary conditionsin which is a linear operator, is a nonlinear differential operator and is a known function, is independent variable, is an unknown function, is the domain of interest and is a boundary operator. The linear operator does not necessarily coincide in it's entirely with the linear part of the equation in study.

If will be the approximate solution of Eqs. (7) and (8), this is a sum of two terms [24], [25], [26], [27], [28], [29] where are unknown parameters, being an arbitrary positive integer numbers.

Inserting Eq. (9) into Eq. (7), one get:

is the initial approximation and can be determined from the linear equation

is the first approximation and is obtained from the nonlinear differential equationwith the initial/boundary conditions

The nonlinear term of Eq. (12) is developed in the formwhere and means the differentiation of order of the nonlinear operator . To avoid the difficulties that appear is solving of the nonlinear equation obtained from Eqs. (12) and (14):and to accelerate convergence of the first approximation and the implicitly of the approximate solution , instead to solve Eq. (15), we make the observations.

In general, the solution of linear differential Eq. (11) can be expressed as:where the coefficients , the functions and the positive integer are known.

The nonlinear operator calculated for may be written in the formwhere the coefficients , the functions and the positive integer are known and depend on the initial approximation and also on the nonlinear operator .

In the following, because Eq. (15) is very difficult to solve, we do not solve this nonlinear equation, but from the theory of differential equations and taking into consideration especially the method of variation of parameters, Cauchy method, the method of influence functions, the operator method and so on [30], it is more convenient to consider the unknown first approximation as depends on and . More precisely, can be determined from the linear equation:where , are unknown parameters and are so-called auxiliary functions depending on the functions and defined through Eqs. (16) and (17) respectively. These functions and are source for the auxiliary functions . We have a great freedom to choose the values for integer positive and of the auxiliary functions . We note that the initial/boundary conditions could be fulfilled by (18) so that finally Eq. (9) responds to all initial/boundary conditions given by Eq. (8). The approximate analytical solutions of Eqs. (7) and (8) is obtained from Eqs. (10) and (18).

Finally, the unknown parameters , can be optimally identified via rigorous mathematical procedures such as: the least square method, Ritz method, Galerkin method, collocation method and Kantorovich method [24-29]. The preferred approach would be the minimize the square residual error, by computing the functionalwhere is domain interest and the residual is given byand is obtained from Eq. (9).

The unknown parameters , can be identified from the conditions:

In this way, the parameters (namely optimal convergence-control parameters or convergence-control parameters) and the auxiliary functions are known.

With these parameters known, the approximate solution is well determined.

In general, for the classical mathematical methods known in the literature, the equilibrium and stability points are known. In sub cases 4.1 and 4.2 we will specify the existence of equilibrium points from which only one is asymptotically stable. Our procedure applied to the stable point which will lead to the validation and full understanding of the presented procedure. The method is inapplicable to points of instability, because Eq. (21) does not admit solutions in this case.

Equilibrium points

Stability analysis of governing equations is carried out to the disease free equilibrium and endemic equilibrium points. A system has an equilibrium points if there no change in Eqs. (1)–(5) at all time. This means that:or

The disease free equilibrium can be determined from the condition and therefore the equilibrium point of disease free for coronavirus is

Endemic equilibrium point can be obtained from conditions , , , and in used to indicate appearance of the disease spread.

The endemic point becomes:where

Stability of equilibrium points

To study the stability of the equilibrium points and we find the eigenvalues of the Jacobian matrix in points and , we denote -unit matrix of order four and - Jacobian matrix. For the equilibrium of disease-free , we obtain

If , then Eq. (30), has eigenvalue and therefore point is unstable. If then point is asymptotic stable.

For the endemic point we obtainwhere

We remark that all coefficients in Eq. (31) are positive. This is a necessary condition for point to be asymptomatic stabile but it is an insufficient condition. The point in asymptomatic stabile if or

Where is given by Eq. (28). We remark that the condition of stability (33) is independent of , and .

The point is unstable if the condition (33) is not fulfilled.

Reproduction ratio

For SEIR model the reproduction ratio is associated with the reproductive power of the disease and is defined as [10,23]:

If , the disease dies out and if an epidemic occurs, and therefore the reproduction ratio show the threshold for the stability of the disease-free equilibrium point. can be used to estimate the growth of the virus outbreak.

Application of OAFM for SEIR model (1)-(6)

Taking into account Eqs. (1)–(5), the function for Eq. (7), becomes:

The approximate analytic solution is: where , , , and , are unknown parameters at this moment. It is obviously that:

It is should be emphasize that the linear operator for every variable , , , and is not unique. For example the first variable , the linear operator and the function g(t) can be choosen in the forms: where is an unknown real parameter.

Considering only Eq. (42), the initial approximation is obtained from Eqs. (11) and (42) where is known from Eq. (6), and follows that

The nonlinear operator for Eq. (1) is

For Eqs. (2)–(5), the linear operators and the function are respectively where and are unknown parameters. Into expression (43) and (48)–(51) we took into consideration Eq. (41).

From Eqs. (48)–(51) we can determine the initial approximations:

The nonlinear operators corresponding to Eqs. (2)–(5) are respectively: where , , and are unknown parameters at this moment.

Inserting the initial approximations (46) and (52)–(55) into Eqs. (47) and (56)–(59), it follow that

In the expression (60)–(64), we use the identity:such that , , , , are functions depending of , .

The auxiliary function from Eq. (18) can be considered as combination of exponential functions (from Eqs. (46), (52)–(55)) and polynomial functions (from Eqs. (60)–(65)). In this way, Eq. (18) for the first approximations and for every variable , , , , can be written as: where , , , and are unknown parameters.

From Eqs. (66)–(70), we obtain:

The approximate analytical solutions of Eqs. (1)–(6) are of the form

Numerical example

In order to prove Fig. 1, Fig. 2, Fig. 3, Fig. 4 and 5 the accuracy of our technique, we consider that the data for system (1 – 6) are for Italian region of Lombardy with the reproduction ratio , corresponding to epidemic situation. The parameters are [10]:

Fig. 1

Comparison between numerical solution [10] and approximate solution (82)

Numerical solution [10] ________; Approx. solution (82) _ _ _ _ _ _ _ _

Fig. 2

Comparison between numerical solution [10] and approximate solution (83)

Numerical solution [10] ________; Approx. solution (83) _ _ _ _ _ _ _ _

Fig. 3

Comparison between numerical solution [10] and approximate solution (84)

Numerical solution [10] ________; Approx. solution (84) _ _ _ _ _ _ _ _

Fig. 4

Comparison between numerical solution [10] and approximate solution (85)

Numerical solution [10] ________; Approx. solution (85) _ _ _ _ _ _ _ _

Fig. 5

Comparison between numerical solution [10] and approximate solution (86)

Numerical solution [10] ________; Approx. solution (86) _ _ _ _ _ _ _ _

In Fig. 1, Fig. 2, Fig. 3, Fig. 4 and 5 are graphically presented the approximate solution (82)–(86) of the Eqs. (1)–(6) in comparisons with corresponding numerical integration results obtained in [10]:

For every variable, the error where is given by Eq. (35) is , , , and . It is clear that the OAFM is very efficient.

Conclusions

In this work, we have developed the SEIR model applied to the death toll of the coronavirus epidemic that incorporates key feature of this pandemic. The mathematical modelling and dynamics, corresponds to reported data in the Italian region of Lombardy with a reproduction rate of . We assumed that the parameters from the governing equations do not change during the period of three months. We mention that the susceptible population decreases after 20 days with 33 percent and then after 35 days the susceptible curve starts to bend downwards. The peak of the exposed population appears at day 25 and the peak of the infected population corresponds to 30 days. The recovered population increases abruptly after 15 days and then is stabilizing after 50 days. The dead population increases abruptly after 18 days but it is stabilizing after two months.

This mathematical model may be the beginning of the epidemiological dynamics in any other regions or countries. The SEIR model of COVID-19 given by Eqs. (1)–(6) was investigated by means of Optimal Auxiliary Functions Method (OAFM). The approximate analytical results are in very good agreement with the numerical results presented in [10]. Our approach is based upon original construction of the solutions using some auxiliary functions and optimal convergence-control parameters which lead to an excellent precision if comparing with the results obtained by OAFM with numerical solutions.

Any nonlinear differential equations is reduced to only two linear differential equation. The construction of the equations from which it can be determine initial approximation and first approximation it is made in original manner. We have a great freedom to choose the number of the auxiliary functions and of the optimal convergence-control parameters. The values of the these parameters are determined by means of rigorous mathematical tools. Our technique leads to a very accurate results using only one approximations and allows us to control and adjust the convergence of the solutions. We remark the construction of the linear operator and the auxiliary functions. Our procedure is effective, explicit and can be applied to any nonlinear dynamical systems.

Authors' contribution

All the authors have equal contribution in this work.

Funding

No source exists of funding this work.

Declaration of Competing Interest

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