Global dynamics of a multi-strain SEIR epidemic model with general incidence rates: application to COVID-19 pandemic.

This paper investigates the global stability analysis of two-strain epidemic model with two general incidence rates. The problem is modelled by a system of six nonlinear ordinary differential equations describing the evolution of susceptible, exposed, infected and removed individuals. The wellposedness of the suggested model is established in terms of existence, positivity and boundedness of solutions. Four equilibrium points are given, namely the disease-free equilibrium, the endemic equilibrium with respect to strain 1, the endemic equilibrium with respect to strain 2, and the last endemic equilibrium with respect to both strains. By constructing suitable Lyapunov functional, the global stability of the disease-free equilibrium is proved depending on the basic reproduction number R 0 . Furthermore, using other appropriate Lyapunov functionals, the global stability results of the endemic equilibria are established depending on the strain 1 reproduction number R 0 1 and the strain 2 reproduction number R 0 2 . Numerical simulations are performed in order to confirm the different theoretical results. It was observed that the model with a generalized incidence functions encompasses a large number of models with classical incidence functions and it gives a significant wide view about the equilibria stability. Numerical comparison between the model results and COVID-19 clinical data was conducted. Good fit of the model to the real clinical data was remarked. The impact of the quarantine strategy on controlling the infection spread is discussed. The generalization of the problem to a more complex compartmental model is illustrated at the end of this paper. © Springer Nature B.V. 2020.

Introduction

Nowadays, several infectious diseases are still targeting huge populations. They are considered amongst the principal causes of mortality, especially in many developing countries. Accordingly, mathematical modelling in epidemiology occupies more and more an increasingly preponderant place in public health research. This research discipline contributes indeed to well understand the studied epidemiological phenomenon and apprehend the different factors that can lead to a severe epidemic or even to a dangerous pandemic worldwide. The classical susceptible-infected-recovered (SIR) epidemic model was first introduced in [1]. Nevertheless, in many cases, the infection incubation period may take a long time interval. In this period, an incubated individual remains latent but not yet infectious. Therefore, another class of exposed individuals should be added to SIR and the new epidemic model will have SEIR abbreviation. Furthermore, epidemiological studies have revealed that the phenomenon of mutation causes more and more resistant viruses giving appearance of many new harmful epidemics or even new dangerous pandemics. Indeed, H1N1 flu virus is considered as mutation of the seasonal influenza [2, 3]. Also, the late coronavirus disease COVID-19 caused by the severe acute respiratory syndrome-related coronavirus SARS-Cov-2 is classified as a strain of SARS-CoV-1 [4]. Other processes of mutation were observed in many infections such as tuberculosis, human immunodeficiency virus and dengue fever [5-7]. For this reason, the multi-strain SEIR epidemic models present an important tool to study several infectious diseases that include a long incubation period and also various infection strains. The relevance of studying multi-strain models is to find out the different conditions permitting the coexistence of all acting strains. The global dynamics of one-strain SEIR model have been the subject of many investigations by considering either bilinear or nonlinear incidence rates [8-10]. The global stability of two-strain SEIR model have been tackled in [11], the authors include to the studied model two incidence functions, the first one is bilinear, while the second function is non-monotonic. Recently, the same problem was studied in [12] by assuming that the two incidence functions are non-monotonic. It is worthy to mention that the incidence rate gives more information of the disease transmission. Hence, the general incidence function has as goal to represent a large set of infection incidence rates. Accordingly, the purpose of this work is to generalize the previous models by taking into account a multi-strain SEIR model with two general incidence rates. Hence, our study will be carried out on the following two- strain generalized epidemic model:withwhere (S) is the number of susceptible individuals, and are, respectively, the numbers of each latent individuals class, and are, respectively, the numbers of each infectious individuals class and (R) is the number of removed individuals. The parameter is the recruitment rate, is the death rate of the population, and are, respectively, the latency rates of strain 1 and strain 2, and are, respectively, the two-strain transfer rates from infected to recovered. The general incidence functions and stand for the infection transmission rates for strain 1 and strain 2, respectively. The incidence functions and are assumed to be continuously differentiable in the interior of and satisfy the same properties as in [13, 14]:The properties , and , for the both functions f and g, are easily verified by several classical biological incidence rates such as the bilinear incidence function   [1, 15, 16], the saturated incidence function or [17, 18], Beddington–DeAngelis incidence function [19-21], Crowley–Martin incidence function [22-24], the specific nonlinear incidence function [25-29] and non-monotone incidence function [30-34]. The flowchart of the two-strain epidemiological SEIR model is illustrated in Fig. 1. Our main contribution centres around the global stability of multi-strain SEIR epidemic model with general incidence rates. In addition, a numerical comparison between our two-strain epidemic model results and COVID-19 clinical data will be conducted. It will be worthy to notice that the dynamics of SEIR COVID-19 epidemic model with two bilinear incidence functions was tackled in [35], and the authors introduce seasonality and stochasticity in order to describe the infection rate parameters. Taking into account non-monotone incidence function, the technique of sliding mode control was used to study an SEIR epidemic model describing COVID-19 disease [36]. The COVID-19 SEIR epidemiological model with Crowley–Martin incidence rate was studied [37], and the authors study the effect of different parameters on the disease spread. In our numerical comparison with COVID-19 clinical data, we will take into account the three latter incidence functions along with Beddington–DeAngelis incidence rate. A brief analysis to a more complex compartmental epidemic model will be established in Appendix of this paper.

Fig. 1

Flowchart of two-strain SEIR model

The rest of this paper is outlined as follows. In the next section, we will establish the wellposedness of the suggested model by proving the existence, positivity and boundedness of solutions. In Sect. 3, we present an analysis of the model, we calculate the basic reproduction number of our epidemic model and we prove the global stability of the equilibria. Numerical simulations are given in Sect. 4 by using different specific incidence functions, and the comparison between the numerical results and COVID-19 clinical data is conducted in the same section. The generalization of the problem to a more complex compartmental model as well as concluding remarks is given at the end of this paper.

The problem wellposedness and steady states

Positivity and boundedness of solutions

For the problems dealing with population dynamics, all the variables must be positive and bounded. We will assume first that all the model parameters are positive.

Proposition 1

For all non-negative initial data, the solutions of the problem (1.1) exist, remain bounded and non-negative.

Moreover, we have .

Proof

By the fundamental theory of differential equations functional framework (see for instance [38] and the references therein), we confirm that there exists a unique local solution to the problem (1.1).

In order to prove the positivity result, we will show that any solution starting from non-negative orthant remains there forever.

First, letLet us now prove that . Suppose that T is finite; by continuity of solutions, we haveIf before the other variables , , , , R, become zero. Therefore,From the first equation of system (1.1), we havethen,However, from we havewhich presents a contradiction.

If before S, , , and R. Then,Again, from the second equation of the system (1.1) with the fact , we will haveHowever, from and , is positive, then we will haveAlso, if before S, , , , R become zero thenBut from the fourth equation of the system (1.1) with , we will haveSince , we haveThis leads to contradiction.

Similar proofs for , and R(t).

We conclude that T could not be finite, which implies that for all positive times. This proves the positivity of solutions.

About boundedness, let the total populationFrom the system (1.1), we haveand therefore,at , we havethenconsequentlysince , for all , we conclude thatThis implies that N(t) is bounded, and so are Thus, the local solution can be prolonged to any positive time, which means that the unique solution exists globally.

The steady states

In this subsection, we show that there exist a disease-free equilibrium and three endemic equilibria. First, since the first five equations of the system (1.1) are independent of R and knowing that the number of the total population verifies the Eq. (2.14), so we can omit the sixth equation of the system (1.1). Therefore, the problem can be reduced to:withAs usual, the basic reproduction number can be defined as the average number of new cases of an infection caused by one infected individual when all the population individuals are susceptibles. We will use the next generation matrix to calculate the basic reproduction number . The formula that gives us the basic reproduction number is: where stands for the spectral radius, F is the non-negative matrix of new infection cases and V is the matrix of the transition of infections associated with the model (2.19). We haveSo,The basic reproduction number is the spectral radius of the matrix . This fact implies thatwithandWe denotethenandWe call the strain 1 reproduction number and the strain 2 reproduction number.

Theorem 1

The problem (2.19) have the disease-free equilibrium and three endemic equilibria , and . Moreover, we have

The strain 1 endemic equilibrium exists when .

The strain 2 endemic equilibrium exists when .

The third endemic equilibrium exists when and .

In order to find the steady states of the system (2.19), we solve the following equationsFrom where, we obtain

When , we find the disease-free equilibrium

When , we find the strain 1 endemic equilibrium defined as follows Define now a function as follows We have using the conditions and , we deduce that However, Therefore, for , we he have Hence, there exists a unique strain 1 endemic equilibrium with , , and .

When , we find the strain 2 endemic equilibrium defined as follows Define also a function as follows We have using the conditions and , we conclude that However, . So, for , we have Hence, there exists a unique strain 2 endemic equilibrium with , , and .

When , we find the third endemic equilibrium defined as follows where with , and .

Global stability of equilibria

Global stability of disease-free equilibrium

Theorem 2

If , then the disease-free equilibrium is globally asymptotically stable.

First, we consider the following Lyapunov function in :The time derivative is given byWe will discuss two cases:By the above discussion, we deduce that, if and (which means that ), thenThus, the disease-free equilibrium point is globally asymptotically stable when .

If , using we will have , then, Since , we obtain Otherwise, therefore

If , using we will have and then, Since , we obtain From , we have

Global stability of strain 1 endemic equilibrium

For the global stability of , we assume that the function f satisfies the following condition:

Theorem 3

The strain 1 endemic equilibrium is globally asymptotically stable when .

First, we consider the Lyapunov function defined by:The time derivative is given byWe haveTherefore,Then,From and , we have then,From (3.12), we haveby the relation between arithmetic and geometric means, we haveWe discuss two cases:We conclude that the steady state is globally asymptotically stable when and .

If , from , we will have since , we obtain, for the following

If from , we will have since we obtain which implies, By the above discussion, we deduce that if we will have

Global stability of strain 2 endemic equilibrium

For the global stability of , we assume that the function g satisfies the following condition:

Theorem 4

The strain 2 endemic equilibrium point is globally asymptotically stable when .

First, we consider the Lyapunov function defined by:The time derivative is given byWe haveTherefore,thenFrom and , we have then,From (3.22), we haveBy the relation between arithmetic and geometric means, we haveWe discuss two cases:We conclude that the steady state is globally asymptotically stable when and .

If then from , we will have since , from we will have

If then from we will obtain from we will get which implies, By the above discussion, we deduce that if then

Global stability of the third endemic equilibrium

For the global stability of , we assume that the functions f and g satisfy the following condition:

Theorem 5

The endemic equilibrium is globally asymptotically stable when and .

We consider the Lyapunov function defined by:thenIt is easy to verify thatAs a result,Using the following trivial inequalitiesthenby the relation between arithmetic and geometric means, we have from (3.12) we havealso from (3.32) we havewith

We conclude that the steady state is globally asymptotically stable when and .

Numerical simulations

In this section, we will perform some numerical simulations in order to examine numerically the SEIR infection dynamics under various incidence functions and also validating our theoretical results. Indeed, we will restrict ourselves to four cases, the first one is to consider the model (1.1) along with the simplest two bilinear incidence functions and , the second case deals with two Beddington–DeAngelis incidence functions, and , while the third case is devoted to check the impact of choosing both incidence functions under Crowley–Martin incidence form, and . The last case consists of incorporating into the model two non-monotonic incidence functions and .

The stability of the disease-free equilibrium

The disease-free equilibrium is characterized by the extinction of the infection, and the disease cannot invade the population. In this subsection, attention will be focused to the numerical stability of such disease-free equilibrium. Indeed, Theorem 2 points out that we can expect the stability of disease-free equilibrium when the basic reproduction number is less than unity. This permits us to look for the right model parameters in order to check numerically the stability of the first steady sate.

Figure 2 shows the behaviour of the infection for the following parameter values: , , , , , , , , , , , , , , , , and . We clearly see that the solutions of the system, under the various suggested incidence functions (bilinear, Beddington–DeAngelis, Crowley–Martin and non-monotonic functions), converge towards the same disease-free equilibrium point . In this situation, the nature of the incidence function has no effect on the steady-state value. Consequently, it was remarked that the disease-free equilibrium first nonzero component depends only on the birth and death rates of the susceptible individuals and does not depend, in any way, on the incidence functions parameters. Here, the disease dies out, the susceptible reaches their maximum value, and the other variables vanish. With the chosen parameters, we can easily calculate the two- strain basic reproduction numbers; in our case, we will have and for the bilinear incidence functions case; we will have also and for the Beddington–DeAngelis incidences case; and are calculated for Crowley–Martin incidence functions case, and finally, we have and for the non-monotonic incidence functions case. In all these cases, we have the basic reproduction number is less than unity which is consistent with the theoretical result concerning the stability of the disease-free equilibrium .

Fig. 2

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the disease-free equilibrium . The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

The stability of the endemic equilibria

The dynamics behaviour concerning the stability of the strain 1 endemic equilibrium is shown in Fig. 3. Indeed, we have chosen the following parameter values: , , , , , , , , , , , , , , , , and . For the bilinear incidence functions case, the solution converges towards the strain 1 endemic equilibrium (1.8, 1.0667, 0, 0.7111, 0, 1.4222) and with the chosen parameters we have and . For the model with Beddington–DeAngelis, we observe the convergence towards and we have and ; for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium (2.3664, 0.8779, 0, 0.5852, 0, 1.1705) and we have and . For the last non-monotonic incidence functions case, we see in the figure the convergence towards the steady state (2.7223, 0.7592, 0, 0.5062, 0, 1.0123) and we have and . We remark that, for the all four cases, the strain 1 basic reproduction number is greater than unity while the other strain 2 basic reproduction number is less than one which confirms our theoretical findings concerning the stability of the strain 1 endemic equilibrium . This endemic equilibrium is characterized by the vanishing the strain 2 latent and infectious individuals.

Fig. 3

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the strain 1 endemic equilibrium . The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

The strain 2 endemic equilibrium stability is illustrated in Fig. 4. In this figure, the following parameters are chosen: , , , , , , , , , , , , , , , , and . We observe the convergence of the solution towards the steady state (2.9167, 0, 0.8333, 0, 0.3571, 0.8928) for the bilinear incidence functions case and with the adopted parameters we have and . For the model with Beddington–DeAngelis, we observe the convergence towards and we have and ; for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium (3.6876, 0, 0.5250, 0, 0.2250, 0.5624) and we have and . For the last non-monotonic incidence functions case, we see in the figure the convergence towards the steady state (3.2918, 0, 0.6833, 0, 0.2928, 0.7321) and we have and . We remark that, for all the four cases, the strain 1 basic reproduction number is less than unity while the other strain 2 basic reproduction number is greater than one which is in good agreement with our theoretical findings concerning the stability of the strain 2 endemic equilibrium . From the components of this endemic equilibrium, we remark the clearance of the strain 1 latent and infectious individuals.

Fig. 4

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the strain 2 endemic equilibrium . The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

The last endemic equilibrium stability behaviour is depicted in Fig. 5 for , , , , , , , , , , , , , , , , and . As the previous figures, the dynamics for the different suggested incidence functions are illustrated. For the first one, the bilinear incidence function, we observe the convergence towards (0.8167, 0.5196, 0.6757, 0.7422, 0.9652, 1.2806) and with the adopted parameters we have and . For Beddington–DeAngelis case, we observe the convergence towards and we have and , for the third case, Crowley–Martin incidence type, we observe the convergence towards the equilibrium and we have and . About the last non-monotonic incidence functions case, we see the convergence towards the steady state and we have and . We remark that both the two-strain basic reproduction numbers and are greater than unity which confirms our theoretical findings concerning the stability of the last endemic equilibrium . This last endemic equilibrium is characterized by the persistence of all strains latent and infectious individuals. The numerical simulations confirm that the model with general incidence functions encompasses a large number of classical well-known incidence functions. Therefore, the generalized mathematical model can give a wide view about the stability of the different problem equilibria.

Fig. 5

Time evolution of susceptible (top left), the strain 1 latent individuals (top middle), the strain 2 latent individuals (top right), the recovered (bottom left), the strain 1 infectious individuals (bottom middle) and the strain 2 infectious individuals (bottom right) illustrating the stability of the endemic equilibrium . The bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). (Color figure online)

Comparison with COVID-19 clinical data

As we have mentioned in the introduction, the recent pandemic COVID-19 is a multi-strain infection. Therefore, the main interest of this subsection is to compare the numerical simulations resulting from our multi-strain epidemic model with COVID-19 clinical data. We have chosen to make our comparison the Moroccan clinical data during the year 2020 in the period between March 31 and June 20 [39, 40].

Figure 6 shows the time evolution of infected cases for the following parameter values: , , , , , , , , , , , , , , , , and . We observe that a significant relationship exists between the curve representing the COVID-19 clinical data and the numerical simulations resulting from our mathematical model for the different incidence functions. Indeed, a good fit between the infected cases given by the mathematical model and the clinical data is observed. Each numerical result with the incidence function (bilinear, Beddington–DeAngelis, Crowley–Martin or non-monotonic function) fits well the clinical data, especially for a certain period of observation. Hence, the multi-strain mathematical model with generalized incidence functions is more appropriate to represent the studied disease.

Fig. 6

Time evolution of infected cases with the bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). The clinical infected cases are illustrated by magenta circles. (Color figure online)

Figure 7 (left-hand side) shows the time evolution of infected cases for the following parameter values: , , , , , , , , , , , , , , , , and . We see that the solution of our model, under the various suggested incidence functions, converges towards the same disease-free equilibrium point . In this situation, the disease dies out, the susceptible reaches their maximum value and the other variables vanish. Within the chosen parameters, we can easily calculate the basic reproduction number; in our case we will have for the bilinear incidence functions case; we will have also for the Beddington–DeAngelis incidences case; are calculated for Crowley–Martin incidence functions case, and finally, we have for the non-monotonic incidence functions case.

Fig. 7

Time evolution of infected cases with the bilinear incidence functions (dotted red), Beddington–DeAngelis incidence functions (yellow), Crowley–Martin incidence functions (green) and the non-monotone incidence functions (blue). The clinical infected cases are illustrated by magenta circles. (Color figure online)

The COVID-19 disease persistence case is illustrated in Fig. 7 (right-hand side). In this figure, the following parameters are chosen: , , , , , , , , , , , , , , , , and . We remark the convergence of the solution towards the endemic equilibrium for all the taken incidence functions. Indeed, for the case with bilinear incidence functions, we have the basic reproduction number is greater than unity . For other cases, the basic reproduction number is also greater than one; indeed for Beddington–DeAngelis case, we have ; for the third case, Crowley–Martin incidence type, we have . For the last non-monotonic incidence functions case, we have which means that the disease persists. We can conclude that the numerical simulations fit well COVID-19 clinical data. Our numerical simulations reveal two scenarios of evolution for this pandemic, and within the first scenario the disease will die out. The other scenario happens when the basic reproduction number is greater than unity; in this case the disease will persist. In this situation, it will be important to eventually undertake some strategies like quarantine, isolation, wearing of masks, disinfection and if it will be possible vaccination.

The effect of quarantine strategy

In this subsection, we will study the effect of quarantine strategy in controlling the infection spread. To illustrate this effect and for simplicity, we will restrict ourselves to the case of the bilinear incidence function. More precisely, we will take the following incidence forms:where stands for the efficiency of the quarantine strategy concerning the first strain infection rate, while represents the efficacy of the quarantine measures concerning the second-strain infection rate. Our numerical simulations will be performed in order to check the impact of the quarantine strategy for the case of the last endemic equilibrium.

Figure 8 shows the time evolution of the susceptible, the two- strain latent individuals, two-strain infected individuals and the recovered for the following parameter values: , , , , , , , and for different values of the strategy controls and . When no strategy is applied, we find the same result as in Fig. 5. By increasing the efficiency of the quarantine measures, we observe an interesting result. Indeed, the number of the susceptible individuals increases when and increase. For higher values of these two control parameters, the number of two-strain latent and infected individuals is reduced considerably, which means that the quarantine strategy can reduce the infection in efficient manner.

Fig. 8

Effect of quarantine strategy on the SEIR model dynamics

Conclusion

In this paper, we have studied the global stability of two-strain epidemic model with two general incidence functions. The model included six compartments, namely the susceptible, two categories of the exposed, two categories of the infected and the removed individuals, this kind of model takes the abbreviation SEIR. We have established the existence, positivity and boundedness of solutions results which guarantee the wellposedness of our SEIR model. The disease-free equilibrium, the endemic equilibrium with respect to strain 1, the endemic equilibrium with respect to strain 2 and the endemic equilibrium with respect to both strains are given. By using an appropriate Lyapunov functionals, the global stability of the equilibria is established depending on the basic reproduction number , the strain 1 reproduction number and the strain 2 reproduction number . Numerical simulations are performed in order to confirm our different theoretical results. It was observed that the model with a generalized incidence function encompasses a large number of classical incidence functions and it can give more clear view about the equilibria stability. In addition, a numerical comparison between our model results and COVID-19 clinical data is conducted. We have observed a good fit between our numerical simulations and the clinical data which indicates that our multi-strain mathematical model can fit and predict the evolution of the infection. We have observed that strict quarantine measures can reduce significantly the infection spread. It is worthy to notice that a generalization form of this problem to a more complex compartmental model is proposed in Appendix of this paper. For this complex model, we have given its basic reproduction number and each strain reproduction number. Furthermore, the forms of the Lyapunov functionals for this complex compartmental model are also formulated.