A new fractional mathematical modelling of COVID-19 with the availability of vaccine.

The most dangerous disease of this decade novel coronavirus or COVID-19 is yet not over. The whole world is facing this threat and trying to stand together to defeat this pandemic. Many countries have defeated this virus by their strong control strategies and many are still trying to do so. To date, some countries have prepared a vaccine against this virus but not in an enough amount. In this research article, we proposed a new SEIRS dynamical model by including the vaccine rate. First we formulate the model with integer order and after that we generalize it in Atangana-Baleanu derivative sense. The high motivation to apply Atangana-Baleanu fractional derivative on our model is to explore the dynamics of the model more clearly. We provide the analysis of the existence of solution for the given fractional SEIRS model. We use the famous Predictor-Corrector algorithm to derive the solution of the model. Also, the analysis for the stability of the given algorithm is established. We simulate number of graphs to see the role of vaccine on the dynamics of the population. For practical simulations, we use the parameter values which are based on real data of Spain. The main motivation or aim of this research study is to justify the role of vaccine in this tough time of COVID-19. A clear role of vaccine at this crucial time can be realized by this study.

Introduction

COVID-19 disease caused by the novel SARS-CoV-2 coronavirus has become a widespread epidemic without precedents all over the world affecting millions of individuals. The first registered case appeared in December 2019 in Wuhan, China, and it rapidly spread all around China and then all over the world [1]. The World Health Organization (WHO) proclaimed the COVID-19 as a international public health emergency on 31 January and declared it as a pandemic on 11 March. As stated by the WHO report, all around the world, as of February 6th, 2021, there have been 104,370,550 confirmed cases of COVID-19 and 2,271,180 deaths.

Compared to other viral infectious disease such as Influenza whose dynamics can be described using SIR models, COVID-19 has the following characteristics [2]: it has a long incubation period; the transmission of this malady is difficult to control due to the fact that it can be propagated by asymptomatic patients [3] and it can be detected through PCR (Polymerase Chain Reaction) tests. This virus can cause the death, specially in elderly individuals or those with previous affections such as Cancer, Diabetes or arterial hypertension. Most common symptoms of the COVID-19 are fever, fatigue, and a dry cough, pain, stuffy and runny nose, sore throat, and diarrhoea. Those patients who suffer the most severe form of the disease can develop pneumonia [4]. This pandemic has tested public health systems from all countries and different measures such as quarantines, lockdowns and social distancing measures have been taken in order to control and contain the advance of the disease. In order to combact this disease, mathematical models and simulations are playing a very important role for taking preventive measures and anticipating the advance of the disease since these models help to foresee the possible future scenarios, see [5], [6] and the references therein. An investigation about the influence of nonpharmaceutical measures on the dynamics of COVID-19 was considered in [7]. In [8] a model analysing the effects of the super-spreader class in the propagation of the disease was studied.

COVID-19 has caused havoc in many countries. Spain, whose data of infection is going to be considered in this paper, has been one of the most damaged by this epidemic presenting one of the highest rate of infected individuals and deaths per million inhabitants. The first case in Spain was reported at the end of January 2020. At the middle of March, a national lockdown was forced in order to stabilize the epidemic. Since then, more than 2.91 millions of infections have been reported and more than 60,000 people have dead due to Covid-19. Spanish government has imposed strong measures in order to avoid the expansion of the virus among the population. These actions include among others social distancing, compulsory use of nose mask and temporal lockdown of commercial establishments [9]. Since January, mass-vaccination campaigns have already started in several countries being Spain one of them. However, a worldwide distribution of COVID-19 vaccines is a difficult and long procedure and then the effects of such starting vaccination campaigns are not reflected in the number of infections reported daily yet.

Fractional derivatives have been extensively used for modelling due to their nonlocal character. This property makes fractional differential equations a really useful tool for describing phenomena in nature which have memory and hereditary properties. These derivatives present several applications in many fields of science such as chemistry, physics, or engineering, see for instance [10], [11], [12], [13] and the monographs [14], [15], [16]. In one of the research fields in which its use has stood out the most is epidemiology. For instance in [17] Mouaouine et al. introduced a Caputo fractional SIR epidemic model considering a nonlinear incidence rate. The same fractional derivative has been implemented in [18] in a SEIR model where the death rate of the population is density dependent. Recently, researchers in [19] have analysed a non-linear fractional tuberculosis epidemic model. Akgul et al. have proposed some new applications on the dynamics of Coronavirus in [20]. Kumar et al.  [21] have used the modified version of Caputo type fractional derivative for studying the outbreaks of COVID-19 in India by using the real raw data based parameter values. In [22], authors gave a prediction based study on the COVID-19 structure for the Brazilian data. Peter et al.  [23] organized a brief discussion on the structure of non-classical non-linear mathematical framework of Coronavirus in Nigeria by using the Atangana–Baleanu type non-classical derivative with Mittag-Leffler memory effects. Recently, Kumar et al.  [24] defined the structure of an ecological model via new generalized Caputo type non-integer order operator by using a recent version of Predictor–Corrector algorithm. In [25], authors solved a non-linear model of Malaria disease via Caputo–Fabrizio fractional derivative. Every day, number of novel analysis are being obtained by the mathematicians in the field of fractional calculus. Early, Researchers in [26] estimated some new novel results by the generalizations in the frame of weighted non-singular non-integer order integral operators.

The high interest in modelling dynamics of COVID-19 disease via fractional derivatives can be reflected in the literature. In [27] the authors considered a fractional model in terms of the fractal-fractional Atangana–Baleanu derivative for describing the spread of COVID-19 taking into account quarantine and isolation. Using the modified predictor–corrector scheme, a new generalized fractional model was given by Erturk and Kumar in [10] for describing the dynamics of COVID-19. In [28], this scheme is considered for solving a time delay fractional COVID-19 SEIR epidemic model via Caputo fractional derivatives. In [29], a SIR fractional epidemic model in terms of the Mittag-Leffler fractional derivative is considered. Another fractional compartmental mathematical model for the spread of the COVID-19 taking into account super-spreaders was given in [30] using the Caputo derivative and supported by the infection data from Spain and Portugal. The same Caputo–Fabrizio fractional derivative is examined in a model for the transmission of COVID-19 in Wuhan by employing the Adams–Bashforth numerical scheme in [31].

After getting the high motivation and brief ideas of modelling approaches, dynamical behaviours and applications of the fractional derivatives from the above mentioned works, the given research paper is structured as follows: In ‘Premilaries’, we provide some important definitions and lemmas. The brief discussion of a new model structure is given in ‘Model structure in classical sense’, where we first define the model in classical sense, in which we calculate the disease-free equilibrium, endemic equilibrium and the basic reproductive number . In ‘Model structure in fractional sense’, we formulate the model in Atangana–Baleanu fractional derivative sense. We also show the existence of a unique solution by using Picard–Lindelof method. The derivation of the solution for the proposed model is obtained by the application of Predictor–Corrector algorithm with its stability analysis. A brief discussion on the model dynamics by the graphical work is given in ‘Graphical simulations’. We simulate the graphs at different transmission rate and vaccine rate by calculating the different basic reproductive numbers. A conclusion section is included at the end of the paper.

Preliminaries

Some preliminary results are discussed here.

[32]

For any given function , where and , the non-integer order Atangana–Baleanu (AB) derivative is given as follows: where satisfying = 1 denotes the normalization function and is the one-parameter Mittag-Leffler function.

The non-integer order AB integral for normalization function is given as

[33]

If and is an integer (non-negative), then there exist the constants and only dependent on , s.t and

Let us assume & for & is a integer. Let for . If then where is a constant independent of and .

Model structure in classical sense

We divide the total population into four compartments: for susceptible humans, for infected human in latent stage i.e. who are infected but not yet become infectious and cannot pass the infection to other susceptible individuals, for infectious and symptomatic humans and who are capable to transmit the infection to any susceptible individuals, and for recovered humans who are temporarily immune from the infection. We add a last compartment namely to denote the total number of individuals who dying in the studied country. is the transmission rate of infection from an infectious human to a susceptible human, and is the total number of human at time . Exposed individuals become infectious at a rate . Infectious humans, after appropriate treatment move to the recovered class at a rate . Recovered individuals can lose their immunity and become susceptible at a rate . We denote the natural death rate of susceptible, exposed, infected and recovered individuals by .

The above assumptions lead to the following non-linear system of five ordinary differential equations given by The description of state variables and parameters of the model are consigned in Table 1, Table 3, respectively.

Table 1

Description of state variables of the SEIRS disease model (3).

State variable	Description
S	Number of susceptible individuals
E	Number of exposed individuals
I	Number of infectious individuals
R	Number of recovered individuals
D	Number of death individuals
N	Total Population

Table 3

Description and values of model parameters.

Parameter	Description	Values	Source
Λ	Rate of birth	0	Fitted
ν	Newborns vaccination proportion	0	Fitted
ν_r	Rate of vaccination of susceptible individuals	0.05,0.1,0.5,0.8	Fitted
ϖ	Rate of transmission from recovered to susceptible class	0.2	[34]
ψ	Rate of transmission from susceptible to exposed	0.542,0.315,0.271,0.192	[34]
ϱ	Rate of transmission from exposed to infected	0.1924	[34]
κ	Recovery rate of infected individuals	0.095	[34]
d	Natural death rate or deaths by any other disease	0	Fitted
l	Death rate due to COVID-19	0.11	[34]
S(0)	Initial population of S	47,431,256	[34]
E(0)	Initial population of E	100	[34]
I(0)	Initial population of I	20	[34]
R(0)	Initial population of R	2	[34]
D(0)	Initial population of D	1	[34]

The total population size is defined by . Then i.e using the linearity of the ordinary derivatives, we can define the feasible region of the system as follows: So the population is given as: with time-varying population .

So the following set is positively invariant for the model system (3).

For further simulations in an easy description, let , , , , then, .

Using this notation, the system (3) becomes Using the fact that , the above system can be written as For readability, let us set , , , and . Thus, (5) becomes

The disease-free equilibrium (DFE)

The local stability is analysed to determine the disease-free equilibrium (DFE):

Through substitution of (7) into of the system in (6), the DFE in (7) is computed to generate or .

At , the Jacobian matrix is given as: with eigenvalues which are solutions of the following equation The characteristic polynomial of the above Jacobian matrix can be written as: Let , and .

Based on the Routh–Hurwitz criteria for a cubic polynomial , the three conditions must be satisfied by the DFE in (7) to be locally stable. So, and and From , we can write the relation as After dividing by , (11) is simplified as: where As a result of the Routh–Hurwitz criteria, all the eigenvalues in the polynomial (9) have negative real part. Then we can conclude that the DFE in (7) is locally stable with (10). Thus we claim the following result.

The disease-free-equilibrium of the system (6) is locally asymptotically stable whenever the basic reproduction number is less than unity.

Endemic equilibrium (EE)

To find the endemic equilibrium points, we assume that the system is time independent such that using the properties of proposed derivative we have: From the last two equations, we obtain Now replacing these values into the first equation of the system, we obtain Substituting this value (18) in (16) we get which only makes physical significance if Number of researchers have described the use of basic reproduction number in analysing the stability of the endemic equilibrium states. According to this, when then the disease free equilibrium is unstable and locally asymptotically stable when .

So if (13) then infection disease dies out eventually. Conversely if (21) then the endemic equilibrium is stable and disease free equilibrium is unstable.

The Jacobian matrix for EE can be written as: and Similarly we calculated the for the above equation and as a consequence of the Routh–Hurwitz criteria, all the eigenvalues in the cubic polynomial have negative real parts and we can conclude that the in (20) is locally stable with (21).

Model structure in fractional sense

Number of research articles have been proposed by researchers to study the dynamics of COVID-19 via fractional derivatives. In [35], solution of a non-linear COVID-19 model via Caputo fractional derivative by the application of q-homotopy analysis transform method has been derived. In [10], authors solved a fractional COVID-19 model by a recent Predictor–Corrector algorithm. Kumar et al.  [28], defined a time-delay non-linear model of COVID-19 in the sense of Caputo type non-integer order operator. In [36], the authors studied the dynamics of COVID-19 in Cameroon by a realistic data with the proposal of a new non-linear model. Also in [37], researchers derived the solution of the model with optimal control analysis in Atangana–Baleanu fractional derivative sense. So after getting a strong motivation from these works, the re-formulation of the proposed integer order model (3) in the Atangana–Baleanu fractional derivative sense is defined as follows: Here, is the Atangana–Baleanu non-integer order derivative operator of order .

Existence of solution

After defining any model in fractional sense, the first major step always becomes to prove the existence of the solution. In this section, we prove the existence and uniqueness of solution for the proposed non-linear model (24). Here first we write the above fractional model in a short form by following the generality. The compact form of the model (24) is given as follows: where and are the contraction for and respectively.

Now, we are following the well-known Picard–Lindelof technique to fulfil the existence of unique solution requirements. Now, by applying the AB integral operator, the given system (25), becomes the following Volterra type integral system of order . Now, we get the subsequent iterative algorithm Here we assume that we can get the exact solution by taking the limit as tends to infinity.

For ABC sense we discuss the existence and uniqueness analysis of the given model by using the Picard–Lindelof scheme.

Let , where We consider the Picard operator , defined as follows: where and   . Next we assume the solution of the given fractional non-linear model being bounded within a time period, , where we demand that . Now by the application of fixed point theorem pertaining to Banach space along with the metric, we get .

Now we have, with . Since is a contraction, we have , hence the given operator is also a contraction. Therefore, the given fractional ABC model (24) has a unique set of solution. □

Derivation of solution

Since last few decades, many computational methods have been investigated for solving different dynamical systems or real world phenomena. In these numerical methods, Predictor–Corrector is one of the very famous and effective methods. This method has been derived for mostly all fractional order derivatives. Here we apply this method in Atangana–Baleanu sense for finding the solution of the given model. The complete description of this method form is given in [38]. Here we derive the solution of the proposed fractional model by assuming . Let us consider

The corresponding fractional Volterra integral equation is given by: So following the method mentioned in [38] for , and setting and , for , the Corrector formula of the projected problem is where and The predictor formula is derived as where Thus the corrector formulae for the proposed model (24) are where

Stability analysis

The numerical method (36) – (37) is conditionally stable.

Let and be perturbations of , respectively. Then, the following perturbation equations are obtained by using Eqs. (36), (37) Using the Lipschitz condition, we obtain where . Also, from Eq.(3.18) in [33] we write where . Substituting from Eq. (41) into Eq. (40) results where . is a positive constant which only depends on (see Lemma 1) and is assumed to be small enough. Applying Lemma 2 we conclude , which completes the proof.

Graphical simulations

In this section, we perform the graphical simulation by using the numerical values based on real data from Spain [34]. We have also fitted some parameter values for studying the role of vaccine clearly. Because to the date, vaccines have just come to the world and they are not present in the enough amount yet. In that case, we fixed the population size by taking birth rate and natural death rate equal to zero. If birth rate is zero then the value of newborns vaccination proportion is automatically zero. For justifying the vaccine availability in the given fixed population size, we used four different values of vaccine rate . Also, we have taken four different values of transmission rate which are based on real COVID-19 cases in Spain from February to September. We calculated the value of basic reproductive number for the all given sets of parameters specified in Table 2.

Table 2

Different sets of parameter values.

Set	ψ	ν_r	R_0
1	0.542	0.05	1.3457
		0.1	1.1214
		0.5	0.4806
		0.8	0.3364
2	0.315	0.05	0.7821
		0.1	0.6518
		0.5	0.2793
		0.8	0.1955
3	0.271	0.05	0.6728
		0.1	0.5607
		0.5	0.2403
		0.8	0.1682
4	0.192	0.05	0.4768
		0.1	0.3973
		0.5	0.1703
		0.8	0.1192

From Fig. 1, we show the effect of vaccine rate on the given model classes at different fractional order values . In that case, we choose the transmission rate . We exemplified that when the order of proposed fractional derivative decreases then the population of exposed and infected individuals decrease. Deaths are also plotted in sub- Fig. 1(e). In this case, the value of basic reproductive number is which is the case of endemic equilibrium. As we have already mentioned, our target is to show the role of availability of vaccine at the given transmission rate of COVID-19. In order to complete this target, we change the value of vaccine rate at the fixed value of transmission rate . In Fig. 2, we used and performed the graphs at various fractional order values . We found that in this case, endemic equilibrium condition again exists but when we set in the company of Fig. 3, then the condition of endemic equilibrium changed to disease free equilibrium with the numerical value of basic reproductive number . When we set then the value of which is the best optimal control case for the population at the transmission rate . This case is exemplified graphically in Fig. 4.

Fig. 1

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 2

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 3

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 4

Behaviour of the given fractional model for vaccine rate and transmission rate .

Now we choose the value of transmission rate . From Fig. 5, we show the effect of vaccine rate on the given model classes at different fractional order values . We exemplified that when the order of proposed fractional derivative decreases then the population of exposed and infected individuals went to decrease at a later time period respectively. Deaths are also plotted in sub- Fig. 5(e). In this case, the value of basic reproductive number is which is the case of disease-free equilibrium. Also at this transmission rate, the disease-free equilibrium case exists for every values of vaccine rate. In the group of Figs. 6, 7 to 8, we used several values of and performed the graphs at various fractional order values . We found that in all the cases, infection is under control and population can easily recover from this virus at the fixed transmission rate.

Fig. 5

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 6

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 7

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 8

Behaviour of the given fractional model for vaccine rate and transmission rate .

Similarly, when we choose the value of transmission rate (see Figs. 9, 10, 11, 12) then the disease-free conditions again exist for each values of vaccine rate . From Fig. 9, we show the effects of vaccine rate on the given model classes at different fractional order values . Again we evaluated that when the order of proposed fractional derivative decreases then the population of exposed and infected individuals went to decrease at the later time period respectively. Deaths are also plotted in sub- Fig. 9, Fig. 10, Fig. 11, Fig. 12. In this case, the several values of basic reproductive number are given in Table 2. Again we analysed that in these all cases of set-3, infection is under control and population can easily recover from this virus at the fixed transmission rate.

Fig. 9

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 10

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 11

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 12

Behaviour of the given fractional model for vaccine rate and transmission rate .

Graphical simulations for set-4 (see Table 2) are established in the range of Fig. 13, Fig. 14, Fig. 15, Fig. 16 where the role of vaccine in disease-free equilibrium case can be easily analysed. In the group of sub-figures, we established the graphs of every classes at different fractional order values . It is cleared that when the values of the fractional order derivative vary then the nature of the given classes , and also went to differ versus time variable .

Fig. 13

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 14

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 15

Behaviour of the given fractional model for vaccine rate and transmission rate .

Fig. 16

Behaviour of the given fractional model for vaccine rate and transmission rate .

A clear role of vaccine on all classes of the given model can be seen from Fig. 17, Fig. 18 at a fixed value of fractional order . From the sub- Figs. 17(b)–17(d), we observed that when the vaccine rate increases then the infection decreases and the recovered population increases with respect to time . In that case, the transmission rate is fixed at . Similarly the role of vaccine at the transmission rate can be observed from the family of Fig. 18 for fractional order value .

Fig. 17

Behaviour of the given classes for various vaccine rate when fractional order and transmission rate .

Fig. 18

Behaviour of the given classes for various vaccine rate when fractional order and transmission rate .

From all above graphical simulations, we exemplified that the given model works well to show the role of vaccine on the dynamics of COVID-19 in the population. We observed that in both cases (disease-free equilibrium and endemic equilibrium) vaccine is effective and when the vaccine rate increases then the population goes into a comfortable environment. We used Mathematica software to simulate the above graphs.

Conclusion

The main target of this paper which was to justify the role of vaccine in this tough time of COVID-19 pandemic has been over. In this research article, we have investigated a fresh framed SEIRS dynamical model by including the vaccine rate along with all other necessary parameters. First we have derived the model in integer order sense and after that we modified it in Atangana–Baleanu fractional derivative sense. The main reason to recall non-integer order AB derivative for exploring the dynamics of the given model has been justified more clearly by using plots. We have provided the analysis of the existence of solution for the proposed fractional SEIRS model. We have used the famous Predictor–Corrector algorithm to derive the solution and we have also analysed its stability. We have simulated number of graphs to see the role of vaccine on the population. For graphical investigations, we have used the parameter values which are based on real data of Spain. In the graphs, we have justified how the model will behave when vaccine rate and the fractional order values change. The comparative plots between the given classes are performed at the sufficient time-range. A clear role of vaccine at this crucial time can be observed by this study. In future, this study can be extended to study the role of vaccine with a real rate of vaccine. Also in the future works, present model can be investigated by considering some other fractional derivatives with the description of optimal control problem. This research work may become useful for the medical organizations to predict the future impact of Coronavirus on the world when the vaccine is available.

CRediT authorship contribution statement

Pushpendra Kumar: Investigation, Conceptualization, Data curation, Formal analysis, Methodology, Visualization, Project administration, Resources, Writing - original draft. Vedat Suat Erturk: Conceptualization, Investigation, Supervision, Software, Visualization, Validation, Writing - review & editing. Marina Murillo-Arcila: Formal analysis, Funding acquisition, Writing - review & editing.

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.