key: cord-274209-n0aast22
authors: Yaro, David; Apeanti, Wilson Osafo; Akuamoah, Saviour Worlanyo; Lu, Dianchen
title: Analysis and Optimal Control of Fractional-Order Transmission of a Respiratory Epidemic Model
date: 2019-07-15
journal: Int J Appl Comput Math
DOI: 10.1007/s40819-019-0699-7
sha: 
doc_id: 274209
cord_uid: n0aast22

The World Health Organization is yet to realise the global aim of achieving future-free and eliminating the transmission of respiratory diseases such as H1N1, SARS and Ebola since the recent reemergence of Ebola in the Democratic Republic of Congo. In this paper, a Caputo fractional-order derivative is applied to a system of non-integer order differential equation to model the transmission dynamics of respiratory diseases. The nonnegative solutions of the system are obtained by using the Generalized Mean Value Theorem. The next generation matrix approach is used to obtain the basic reproduction number [Formula: see text] . We discuss the stability of the disease-free equilibrium when [Formula: see text] , and the necessary conditions for the stability of the endemic equilibrium when [Formula: see text] . A sensitivity analysis shows that [Formula: see text] is most sensitive to the probability of the disease transmission rate. The results from the numerical simulations of optimal control strategies disclose that the utmost way of controlling or probably eradicating the transmission of respiratory diseases should be quarantining the exposed individuals, monitoring and treating infected people for a substantial period.

The respiratory syncytial virus, influenza virus and parainfluenza virus are some viruses that cause respiratory diseases. Influenza [1] . Influenza viruses are group into Type A and Type B. The viruses are often transmitted from people and cause seasonal influenza epidemics each year. The evolving prevalence of infectious diseases is increasing every year. It occurs through the respiratory tract. The spread of the disease is rapid and widespread. Since 1980, the disease has appeared just like bird flu in Asia. The outbreak of the 2014 Ebola virus in West Africa, severe acute respiratory infectious (SARS) disease and the Middle East respiratory syndrome (MERS) caused severe impact on the health system. Most respiratory diseases have no vaccines, such as SARS and Ebola. These diseases spread quickly and may be re-infected. After the flu season, there are several different types of influenza (Types A and B) and subtypes (Type A) that circulate and cause disease. For example, the bird flu virus has several subtypes such as H5N1, H7N3, H7N7, H7N9, and H9N2, which can infect humans [1] .

Epidemiological mathematical models have proven to be a valuable tool for understanding, analyzing influenza virus infection dynamics, disseminating and recommending control strategies. Although [2] [3] [4] [5] [6] has completed a great deal of work on dynamic modeling of influenza, it is limited to ordinary differential equations. However, currently, it has been found that the use of fractional differential equations to model many different fields of phenomena has been very successful [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] . For instance, in mathematical epidemiology, Ebola virus epidemic has been modeled with fractional-order differential equations by [19] . They used the SEIR fractional-order model to analyze data published by the WHO to provide projections of outbreaks in three countries in West Africa. The model considered fits precisely with the real data. Their findings revealed that the outbreak will last about two years, with an estimated 9 million infected people. Although individuals heredity play an important role in mortality rate of the disease, data analysis shows that the predicted death toll is very high.

Goufo et al. [20] took into account the stability analysis of a non-linear spread Ebola hemorrhagic fever epidemic model. They used conventional time derivatives to express models that contain new parameters that happen to be fractional. They proved that the model is well defined, poses non-negative solutions and also established the conditions for boundedness. The Routh-Hurwitz criterion was used to show the existence and stability analysis of Ebola virus model equilibrium states and showed that they strongly depend on non-linear propagation. Also, they provide the conditions for the persistence of the Ebola virus in the system. In addition, numerical simulations of non-linear propagation are provided and the results obtained are significant for combating and preventing Ebola hemorrhagic fever, which has so far caused the deaths of hundreds of families and continues to infect many people in West Africa.

Fractional diffusion mimics the human mobility network by simulating disease outbreaks [21] . Human mobility networks can smooth the spread of infectious diseases from the sidewalk to the flight route. The effort of control and removal depends on describing these networks in terms of individual connections and flux rates between the contact nodes. In some cases, transportation can be parametrized by a gravity-type model or approximated by diffuse random walks. Alternatively, they separated domestic commercial air traffic into a case study of the utility of non-diffusion heavy-tailed transport models. A new stochastic simulation of typical influenza-like infections was adopted, targeting the dense, highly connected America air travel network. They show that the mobility on the network can be defined mainly by a power law, which is consistent with the previous study. They observed that the global evolution of an outbreak on the network is precisely reproduced by a twoparameter space-fractional diffusion equation, which is determined by the air travel network. The dynamic changes of cytotoxic T-lymphocyte (CTL) responses in Ebola virus infection in vivo are described by a fractional-order model with time delay [22] . They introduced a time-delay during the CTL response period to represent the time required to simulate the immune system. Stability and Hopf bifurcation were obtained from the model through fractional Laplace transform. Moreover, the stability conditions show that the dynamics of the model can be improved by the fractional-order delay.

In [23] SIR epidemic model with non-integer order under the conditions of external noise is studied. The behavior of the system changes with the introduction of seasonality and noise force. It revealed that different non-integer order and parameters improve the dynamical behavior of the system. Multi-scale fuzzy entropy is used to study the complexity of stochastic models. They designed a hard limiter control system and simulation results which showed that with effective medical and health measures the proportion of infected people can be controlled to significantly small numbers.

Compared with integer order models, non-integer order models are more effective in modeling biological systems which have long-range and temporal memories [24] . In [25] , González-Parra et al. proposed a nonlinear order system to discuss the outbreaks of H1N1 influenza. The fractional model does not only rely on the present stage but also on all its history which makes fractional models more general than the integer-order models. They also determined that the nonlinear fractional epidemiological model can be well matched to provide numerical results that are in good agreement with the actual data for H1N1 influenza. The model suggested gives valuable facts for understanding, predicting and controlling the spread of different epidemics across the globe.

For these reasons, the fact that fractional (non-integer) order models of many phenomena have recently been shown to be more realistic than the integer order, and the fact that the fractional order also has long-range memories motivated this study. In addition, it is also obvious that modeling by fractional order derivative can capture any natural phenomena or a rich variety of dynamics observed in the system. Many of the past models focused only on the analysis of the influenza outbreaks in human populations. But, a model in which the order is an integer to describe the transmission of a respiratory epidemic has recently been proposed in [26] . However, modeling the transmission of a respiratory epidemic by fractional order model might provide a feasible alternative in controlling or probably eradicating the disease.

The present study introduces Caputo fractional-order into SEIR epidemic type of model proposed by [26] . The proposed model is analyzed without control measures to investigate its stability conditions. The sensitivity of the basic reproductive number R 0 is analyzed to determine the most sensitive parameters. The interpretation of the sensitivity led us to two control measures. The optimal control theory is then used to investigate the efficiency of incorporating the control measures, namely, quarantine of exposed population and monitoring and treating of the infected population.

The remaining part of the paper is structured as follows: We provide some basic and necessary definitions about the fractional calculus in section "Fractional Calculus". The description of the model is discussed in section "The Model". The nonnegative solution and the stability analysis of the fractional-order differential equations are discussed in Sections "Non-negative Solutions" and "Model Equilibrium States and Stability" respectively. The optimal control of the epidemic is discussed in section "Optimal Control". A numerical solution of the fractional model using Atanackovic and Stankovic method is discussed in section "Numerical Method". We support our theoretical analysis with numerical simulations in section "Numerical Simulation and Discussion". Finally, the paper concludes in section "Conclusion".

In this section, we define fractional integral and fractional derivative of Riemann-Liouville and Caputo respectively which are applied in this work.

Definition 1 [27] . The order σ > 0 of fractional integral operator function h : R + −→ R defined by

is called the Riemann-Liouville fractional integral.

Here and elsewhere Γ is known as the Euler gamma function and is defined as

Definition 2 [12] . The derivative function of fractional order σ > 0 (where σ lie in the half-open interval (0,1]), defined by

is called the Caputo fractional derivative, where a is the starting point.

We use the Caputo derivative definition in this work. The initial conditions for Caputo defined fractional differential equation is the same as ordinary differential equations which is a core advantage.

The model we studied in this work is proposed by [26] . The model considers four variables, namely, the population which is susceptible ( S(ξ )), the population which is exposed (E(ξ )), the infected population (I (ξ )), and the recovered population (R(ξ )). According to [26] , b(ξ ) represents the rate of new susceptible people entering the population at time ξ , β represents the probability of disease transmission, ν represents the seroconversion rate, α represents the recovery rate, μ N represents the natural mortality rate, μ D represents the disease induced death rate, and κ represents the rate at which the recovered return to the susceptible population (due to the loss of immunity). These assumptions lead to the following integer-order differential equation presented by [26] ,

The system (4) in fractional order is given by

with initial conditions

where D σ * a = d σ dξ σ is the Caputo fractional derivative of order σ . It is important to notice that when the fractional-order σ −→ 1 , the system (5) becomes the integer-order system (4).

We need the following Lemma in [28] to proof non-negative solutions of system (5).

* a h(m) ∈ C(d, e] for 0 < σ ≤ 1 and d, e ∈ R then h(m) = h(d) + 1 Γ (σ ) D σ * a h( )(m − d) σ , with d ≤ ≤ 1, ∀m ∈ (d, e]. Remark 4.1 Assume h(m) ∈ C[0, e] and D σ * a h(m) ∈ C(0, e] for 0 < σ ≤ 1. It follows from Lemma 4.1 that if D σ * a h(m) ≥ 0, ∀m ∈ (0, e) then h(m) is non-decreasing ∀m ∈ [0, e], and if D σ * a h(m) ≤ 0, ∀m ∈ (0, e) then h(m) is non-increasing ∀m ∈ [0, e].

Proof Using Theorem 3.1 together with Remark 3.2 in [29] , we can obtain the existence and uniqueness of the solution of (5)- (6) in (0, ∞). We prove that the domain R 4 + for the model is positively invariant, since

around non-negative domain or neighborhood on each hyperplane, the vector field points to R 4 + .

For the equilibrium states of the system of fractional order model (5), let

Then the equilibrium states are

The diseases free equilibrium state K 0 , is where the infectives equal to zero (I = 0) and the endemic equilibrium state K 1 , is where the infectives is nonzero (I = 0).

The basic reproduction number of the system (5) using the next generation matrix approach, given by

where F is non-negative and V is a non-singular M-matrix. Applying this method on the system (5), where

Now, finding F and V , at disease free equilibrium, K 0 , and using K = F V −1 we have

The Jacobian matrix J (K 0 ) for the system (5) evaluated at K 0 is given by

The disease free equilibrium state K 0 of system (5) is locally asymptotically

Proof The disease free equilibrium state K 0 is asymptotically stable locally given that the eigenvalues Λ l , l = 1, 2, 3, 4, of J (K 0 ) satisfy the conditions [30, 31] :

We can evaluate these eigenvalues by solving the following characteristic equation

this leads to the equation

The characteristic equation gives the roots

Obviously, P + Q > 0, and if P Q > W , then all the eigenvalues Λ l , l = 1, 2, 3, 4, satisfy the condition given by (9) .

The basic reproduction number denoted by R 0 is

value and is defined as the number of cases occurring in a population which is completely susceptible. Biologically, if R 0 is less than one, then the infection will disappear, but if it is more than one, the infection still exists. For the discussion of the asymptotic stability of the persistence of the disease of system (5), we need the following definition and Lemma 5.1.

Definition 3 [32] .The discriminant D(h) of a polynomial h(y) = y m + c 1 y m−1 + c 2 y m−2 + · · · + c m (10)

is the determinant of the corresponding Sylvester (m + ξ) ⊗ (m + ξ) matrix, the Sylvester matrix is formed by filling the matrix beginning with the upper left corner with the coefficients of h(y) and then shifting down one row and one column to the right side. The process is then repeated for the coefficients of f (y).

Lemma 5.1 [32] . For the polynomial equation,

the conditions displayed below make all the roots of (11) satisfy (9):

1. for m = 1, the condition for (11) is b 1 > 0 2. for m = 2, the conditions for (11) are either Routh-Hurwitz conditions or b 1 < 0

3. (11) is satisfied for all σ ∈ [0, 1). 6. For general m, b m > 0 is a necessary condition for condition (11) to be satisfied.

The Jacobian matrix J (K 1 ) calculated at the disease persistence equilibrium is given as:

By using the characteristic equation det(J (K 1 ) − ΛI ) = 0, the linearized system above is in the form

where 

Lemma 5.2 From condition (6) in Lemma 5.1 the positive equilibrium point K 1 is locally asymptotically stable, since the polynomial H (Λ) as given in (12) has a coefficients d 1 , d 2 , d 3 , and d 4 positives.

Here, we investigate the response of R 0 to parameter changes and determine the effect of each parameter on R 0 and the potential for effective control and elimination of the disease. It is straightforward to calculate the partial derivatives of the value of R 0 using Eq. (8) with respect to the parameters β, ν σ , b σ , μ σ N , μ σ D and the recovery rate α σ . With all other parameters held constant, the elasticity E x (or the variable's normalized forward sensitivity index) approximates the fractional change in R 0 that results from a unit fractional change in parameter x, defined as

This index shows how sensitive R 0 is to changes of parameter x. Specifically, a positive (negative) index shows that an increase in the parameter value results in an increase (decrease) of R 0 [33] . The elasticities for the quantities of interest are Figure 1 indicate that, R 0 is most sensitive to β the disease transmission rate, followed by α σ the recovery rate. The seroconversion rate ν σ and the natural death rate μ σ N have the same sensitivity index. It can also be observed that, R 0 is least sensitive to μ σ D , the disease induced death rate. In detail, the sensitivity indexes for α σ , β, ν σ , μ σ D , and μ σ N , are found to be −0.9999988, 1, 0.0058, −0.00000122, and −0.0058 respectively, once all parameters are fixed at their baseline values (Fig. 1) . Thus, for instance, if the rate of recovery were to increase or (decrease) by 10%, then the value of R 0 would decrease or (increase) by 9.999988%. Likewise, a 10% increase or (decrease) of the disease transmission rate would correspond to a 10% decrease or (increase) of the R 0 , 10% increase or (decrease) of seroconversion rate would decrease or (increase) the R 0 by 0.058%, 10% increase or (decrease) of the disease induced death rate would correspond to decrease or (increase) the value of R 0 by 0.0000122% and 10% Fig. 1 The sensitivity analysis of the basic reproductive number increase or (decrease) of the natural death rate would correspond to a 0.058% decrease or (increase) in the value of R 0 . Therefore, the above interpretations recommend that control strategies that can efficiently decrease the probability of disease transmission β, natural death rate μ σ N , disease induced death rate μ σ D , should be used to control the disease transmission effectively. Additionally, increase the rate of recovery α σ will lead to a decrease in R 0 , thus all control strategies that can effectively help reduce the transmission of the respiratory diseases should be applied. The mathematical perspective for these strategies would be detailed in our subsequent model.

In this section, we extend our model in Eq. (6) by introducing two time-dependent control measures, namely u 1 (ξ ) (quarantine of exposed population groups) and u 2 (ξ ) (monitoring and treatment of infected populations ). It is assumed that the exposed population is reduced by the factor (1 − u 1 (ξ )) as they are quarantine. Furthermore, the infected population is reduced by a factor of (1 − u 2 (ξ )) as they are monitored and treated by health professionals. The model system (6) becomes

where E is the exposed population and I is the infected population. T is the final time and the coefficients c 1 , c 2 , c 3 , c 4 are positive weights. Our aim is to minimize the exposed and infected population while minimizing the cost of control u 1 , u 2 . Thus, we search for an optimal control u * 1 , u * 2 , such that

where the control set is

The terms c 1 E and c 2 I represent the cost of reducing the exposed and infected population respectively, while c 3 u 2 1 is the cost of quarantine and also, c 4 u 2 2 is the cost of monitoring and treatment. The necessary conditions that an optimal control must satisfy come from the Pontryagin's Minimum Principle [34] [35] [36] [37] . This principle converts Eqs. (6) and (18) into a problem of point-wise minimizing a Hamiltonian M with respect to (u 1 , u 2 ) stated as follows:

where λ S , λ E , λ I and λ R , adjoint variables or co-state variables [34] [35] [36] [37] 

The transversality conditions are λ S (T ) = λ E (T ) = λ I (T ) = λ R = 0. On the interior of the control set, where 0 < u i < 1, for i = 1, 2 we have

We obtain

The control parameters (u * 1 , u * 2 ) that minimizes J (u 1 , u 2 ) over U is given by

where λ S , λ E , λ I and λ R are the adjoint variables satisfying (6) and the following transversality conditions:

Atanackovic and Stankovic [38] numerical method FDE is discussed in this section. This method indicates that the fractional derivative of a function f (ξ ) with order σ may be stated as

where

with the following properties:

By using K terms with the sums appearing in (20) , we can approximate D σ g(ξ ) as

The above equation can be written as D σ * a g(ξ ) Ω(σ, ξ, K )g (1) 

where

,

We set

for k = 2, 3, . . . We can rewrite system (7) as

where

Now we can rewrite (23) and (25) as follows: 

The system (27) with (28) can be solved numerically by using the Runge-Kutta fourth order method. Fig. 4 The simulations displaying the result of u 1 (quarantine of the exposed population groups) and u 2 (monitoring and treatment of the infected people) on a exposed populations, b infected populations

In this section, we present numerical simulations to confirm the theoretical results obtained in the preceding section. By using the well known generalized Euler method (GEM) [39] with values in Table 1 , we simulate system (5) .

For the parameter values in Table 1 above and by calculation, we obtained R 0 = 1.5143 and the endemic equilibrium K 1 = (1679.128, 858.647, 4.883, 3 .712). We obtained R 0 = 1.5143, D(P) = 2.389 > 0 and the simulations in Fig. 2 show that the endemic equilibrium K 1 is positive and locally asymptotically stable for σ = 0.70, σ = 0.80, σ = 0.90 and σ = 1.0, satisfying condition 6 of Lemma 5.1. It can obviously be seen in Fig. 2 that, compared with the situation of order σ = 0.70 and σ = 0.80, the trajectory of the system with order σ = 0.90 is nearer to the trajectory of the system with the order σ = 1.0. Thus, the bigger of the trajectory difference, the distant from σ to 1.0.

It can be seen from Fig. 3 when σ = 1 system (5) is the classical integer-order system (4). The display of the trajectories indicates the behavior of the approximate solutions for system (5) obtained for the value of σ = 1. 

It can be seen from Fig. 4 that the optimal control u 1 and u 2 has a significant effect on the exposed and the infected populations respectively. Infection level is reduced rapidly but not eliminated. This suggests that monitoring and treatment strategies that can allow the immune response to rebuilding should also be well-thought-out (Figs. 5, 6 ). 

We introduced Caputo fractional-order into a classical integer-order model proposed by [26] to model the transmission dynamics of respiratory disease. The nonnegative solution of the model is provided by using the generalized mean value theorem. We obtained the basic reproductive number R 0 , which perform as a threshold parameter in the disease control. We established and investigated the stability analysis of the fractional-order model with respect to the values of R 0 . The disease-free equilibrium is locally asymptotically stable if R 0 < 1. For R 0 > 1, using Lemma 5.1 and condition 6, we investigated the local stability of the positive endemic equilibrium state. Sensitivity analysis shows that R 0 is most sensitive to the disease transmission rate β. This recommends that periodic monitoring by medical professionals and researchers should be done to control the transmission of the disease. Additionally, we investigated the optimal control problem by the application of the optimal control theory. We used the Pontryagin's Minimum Principle to provide the necessary conditions needed for the existence of the optimal solution to the optimal control problem. Also, we applied the Atanackovic and Stankovic method to provide a numerical solution to system (5) . Lastly, the theoretical results were verified by numerical simulations to measure the efficacy and impact of control on the transmission of the respiratory diseases. From the numerical simulation, the size of the infected population is significantly reduced under the controlled conditions. This proposes that if all two control measures u 1 (quarantine of exposed population groups) and u 2 (monitoring and treatment of infected populations) are employed for the same period of time and continue for a considerable period of time, the future of free from transmission of a respiratory disease could be achieved. In this manner, the fractional-order optimal control method can progress the value of the treatment (Fig. 7) . The major advantages of our proposed fractional order model which cannot be exhibited by classical order model are:

-It's highly effective and efficient which help us to obtain better results. -It's easy to implement.

-It provides improved precision of the process model by offering more flexibility in model identification. -By modeling system by fractional order we can model a higher system by low order model. -It has the effect of memory, which is essential factor in many biological processes.

Simple model for respiratory diseases

Modeling the initial transmission dynamics of influenza A H1N1 in Guangdong Province

A dynamic compartmental model for the Middle East respiratory syndrome outbreak in the Republic of Korea: a retrospective analysis on control interventions and superspreading events

Application of the CDC EbolaResponse modeling tool to disease predictions

Parameter estimation of Influenza epidemic model

Global dynamics of avian influenza epidemic models with psychological effect

A new fractional analysis on the interaction of HIV with CD4+ T-cells

New aspects of the poor nutrition in the life cycle within the fractional calculus

New aspects of the motion of a particle in a circular cavity

Suboptimal control of fractional-order dynamic systems with delay argument

On the nonlinear dynamical systems within the generalized fractional derivatives with Mittag-Leffler kernel

A fractional calculus based model for the simulation of an outbreak of dengue fever

Two-strain epidemic model involving fractional derivative with Mittag-Leffler kernel

Functionality of circuit via modern fractional differentiations

A fractional model of vertical transmission and cure of vector-borne diseases pertaining to the Atangana-Baleanu fractional derivatives

The effect of anti-viral drug treatment of human immunodeficiency virus type 1 (HIV-1) described by a fractional order model

A fractional order SEIR model with vertical transmission

Some properties of the Kermack-McKendrick epidemic model with fractional derivative and nonlinear incidence

On a fractional-order Ebola epidemic model

Stability analysis of epidemic models of Ebola hemorrhagic fever with non-linear transmission

Fractional diffusion emulates a human mobility network during a simulated disease outbreak. Institute for Disease Modeling, Intellectual Ventures

A fractional-order delay differential model for Ebola infection and CD8 + T-cells response: stability analysis and Hopf bifurcation

Epidemic outbreaks and its control using a fractional order model with seasonality and stochastic infection

A fractional-order model for Ebola virus infection with delayed immune response on heterogeneous complex networks

A fractional-order epidemic model for the simulation of outbreaks of influenza A(H1N1)

Applied Mathematics for the Analysis of Biomedical Data: Models, Methods, and MATLAB

Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation

Generalized Taylor's formula

Global existence theory and chaos control of fractional differential equations

Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models

Stability results for fractional differential equations with applications to control processing

On fractional-order differential equations model for nonlocal epidemics

Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model

An off-line NMPC strategy for continuous-time nonlinear systems using an extended modal series method

A novel feedforward-feedback suboptimal control of linear timedelay systems

An efficient finite difference method for the time-delay optimal control problems with time-varying delay

The Mathematical Theory of Optimal Processes

Stability analysis for a fractional HIV infection model with nonlinear incidence

An algorithm for the numerical solution of differential equations of fractionalorder

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: a prime for parameter uncertainty, identifiability, and forecasts

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations