Fractional optimal control in transmission dynamics of West Nile virus model with state and control time delay: a numerical approach

In this paper, an optimal control for a novel fractional West Nile virus model with time delay is presented. The proposed model is governed by a system of fractional delay differential equations, where the fractional derivative is defined in the Grünwald–Letnikov sense. Stability analysis of fixed points is studied. Corresponding fractional optimal control problem, with time delays in both state and control variables, is formulated and studied. Two simple numerical methods are used to study the nonlinear fractional delay optimal control problem. The methods are standard finite difference method and nonstandard finite difference method. Comparative studies are implemented, it is found that the nonstandard finite difference method is better than the standard finite difference method.


Introduction
Time delays play an important role in the modeling of real-life phenomena in different fields of applications, for example, engineering, chemical, and biological systems that are modeled by differential equations. Many methodologies have been proposed to deal with the control and optimal control problems of systems with time delay [1][2][3][4][5][6][7][8][9]. In [10], an introduction to time-delay control problems is given. In [11], Kharatishvili established a maximum principle for optimal control problems with a constant state delay. Similar results for control problems with pure control delays are presented in [12]. In [13], Halanay proved a maximum principle for optimal control problems with multiple constant delays in state and control variables. Soliman and Ray obtained the similar results in [14]. A simple method for obtaining necessary conditions for control problems with a constant delay in the state variable is introduced in [15]. A maximum principle for control systems with a time-dependent delay in the state variable is derived in [16]. The theory of necessary conditions for optimization problems in functional spaces is applied on a unified approach to control problems with delays in the state variable in [17]. Optimal control for a time-delay multi-strain tuberculosis fractional model is introduced in [18]. Recently, many interesting studies were presented in nonlinear and time-delay optimal control problems, see [19][20][21].
Fractional calculus has been the focus of many studies due to its appearance in various sciences [22][23][24][25][26][27][28][29][30]. The paths of the control variables that maximize or minimize a function of the state and control variables are described by fractional optimal control problems (FOCP). The theory of FOCPs has been under development. The approximation and numerical techniques are used to solve the FOCPs [25,31]. Since most of the fractional-order differential equations do not have exact analytic solutions, the approximation and numerical techniques must be used [26]. Delayed fractional differential equations (DEDEs) are also used to describe dynamical systems [32], for more details, see [33][34][35]. Recently, many papers have been devoted to DEDEs (see [36][37][38], and the references cited therein). Nonstandard finite difference method (NSFDM) was firstly developed in 1980s by Mickens and has been applied in many areas of science including biological and epidemic models [39][40][41]. It is important that the proposed numerical method can preserve the positivity of the solution [42].
In this paper, fractional time-delay optimal control in transmission dynamics of WNV model is presented as an extension of the WNV model given in [26]. Time delays in both state and control variables are formulated and studied. The aim of this paper is to study the fractional time-delay optimal control of this proposed model. We introduced NSFDM and standard finite difference method (SFDM) to study the behavior of this optimal control problem.
West Nile virus (WNV) is an arbovirus and single-standard RNA virus Flavivirus [43]. WNV is transmitted among mosquitoes, birds, humans, and other animals [44]. In 1937, the virus was first isolated in the West Nile district of Uganda [45]. It has spread in Africa, Europe, the Middle East, west and central Asia, Oceania (subtype Kunjin), and North America [46][47][48][49][50]. Owing to the absence of an effective diagnostic test therapeutic treatment against WNV, public health officials have focussed on the use of preventive measures in attempt to halt the spread of WNV in humans [43]. In 2001 [51], the first WNV model was presented by Thomas and Urena. In [52], a single season model with a system of differential equations for WNV transmission in the mosquito-bird population was presented by Wonham et al. They were using local stability results and simulations, showed that while mosquito control decreases WNV outbreak threshold, controlling birds increases it. In [53] Cruz-Pacheco et al. studied a mathematical model for the transmission of WNV infection between mosquito and avain populations, and by using experimental and field data as well as numerical simulations, they found the phenomenon of damped oscillations of the infected bird population. Also, the spatial spread of the virus was established by Lewis et al. in [54]. A comparative study of the discrete-time model in [51] and the continuoustime model in [55] was introduced by Lewis et al. in [56]. Kbenesh et al. determined the cost-effective strategies for combating the spread of WNV in a given population in [57]. Fractional optimal control of WNV model was introduced by using two numerical methods by Sweilam et al. in [26]. This paper is organized as follows. In Sect. 2, preliminaries and notations are presented. In Sect. 3, a mathematical model is formulated. Reproduction number R 0 and stability of the disease-free equilibrium and the endemic equilibrium are studied in Sect. 4 and Sect. 5, respectively. The optimal control problem is formulated in Sect. 6. In Sect. 7, NSFD method for fractional time-delay WNV model is presented. A numerical experiment is discussed in Sect. 8. Finally, in Sect. 9, conclusions are presented.

Preliminaries and notations
Definition 2.1 ([58]) The left Caputo fractional derivative (LCFD) is defined as follows: the right Caputo fractional derivative (RCFD) is defined as follows: where Γ is the Euler gamma function.

Definition 2.2 ([58])
The Grünwald-Letnikov approximation of the fractional derivative of order α of the function V(t) is defined as follows: where [t] is the integer part of t, h is the step size, and the Grünwald-Letnikov coefficients are given by

Definition 2.3 ([59])
The numerical scheme is called NSFD discretization if at least one of the following conditions is satisfied: (i) Nonlocal approximation is used [39,59,60].
(ii) The discretization of derivative is not traditional and a nonnegative function , called a denominator function, is used [41]. The discretization is based on the approximations of temporal derivatives by a generalized forward scheme. If we have ϕ( t) = t, then the scheme is defined as finite difference or SFDM. Hence, if M(t) ∈ C 1 (R), the first derivative can be defined as follows:

Mathematical model
In this section, a new formulation for the transmission dynamics of WNV model given in [57] is presented. In the following, the transmission dynamics of WNV model is modified to be nonlinear time-delayed fractional differential equations. The WNV introduced in [57] consists of nine nonlinear ODEs which is described by the total mosquito population at time t, denoted by N M (t), split into the populations of susceptible M s (t) and infected M i (t) mosquitoes. Similarly, the total bird population at time t, denoted by N B (t), is subdivided into susceptible B s (t) and infected B i (t) birds. Finally, the total human population The population of infected mosquitoes N M (t) The total population of mosquitoes The population of susceptible birds The population of infected birds N B (t) The total population of birds The population of susceptible humans The population of exposed humans The population of infected humans The population of hospitalized humans The population of recovered humans N H (t) The total population of humans and  The probability of WNV transmission from an infected bird to a susceptible mosquito β α 2 The probability of WNV transmission from an infected mosquito to a susceptible bird β α 3 The probability of WNV transmission from mosquitoes to humans b α 1 The per capita bitting rate of mosquitoes on the primary host (birds) and b The per capita bitting rate of mosquitoes on the humans and b The average bitting rate of mosquitoes The rate of development of clinical symptoms of WNV r α The natural recovery rate τ α The treatment-induced recovery rate γ α The hospitalized rate of infectious humans Table 3 Parameter values used in system (6) Parameter Value Parameter Value with the following initial conditions: where c a D α t is the Caputo fractional-order derivative. Let us consider that . The variables in system (6) are defined in Table 1. The list of parameter values and their interpretation are introduced in Tables 2 and 3; for more details, see [57].

Reproduction number R 0 and stability of the disease-free equilibrium
In this section, we compute the reproduction number R 0 . The basic reproduction number is the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [32]. Method in [32] is applied to derive R 0 . The order of infected variables is The new infections into the infected classes ς are defined by the vector F, and the other flows within and out of the infected classes ς are defined by the vector V, where these vectors are given by Let us differentiate the above matrices F and V with respect to the infected variables ς and evaluate at the disease-free equilibrium then the basic reproduction number R 0 for system (6) with spectral radius ρ is defined by
Proof To analyze the stability of the disease-free equilibrium (DFE) point, let us have that The characteristic equation can be calculated at the DFE point. It can be written as follows: . Thus the characteristic equation corresponding to the above matrix is [62] |J ε eq -λI| = 0, If we have d = 0, then the eigenvalues of the Jacobian matrix J ε 0 at the DFE point ε 0 are , and λ 9 = -(d α I + γ α + μ α H + r α ). According to the Routh-Hurwitz theorem [35], these roots are negative or have negative real parts and the eigenvalues satisfy Matignon's conditions [63] defined by | arg λ i | > απ 2 . Thus, the DFE point ε eq is asymptotically stable. If we consider d > 0, then (λ + d α H + μ α H + τ α e -λd ) has no pure imaginary roots for any value of the delay d. Hence, all roots of the characteristic equation have negative real parts, and the DFE point ε eq is locally asymptotically stable according to Matignon's conditions [63] defined by | arg λ i | > απ 2 .

Stability analysis of the endemic equilibrium
In this section, the stability analysis of the endemic equilibrium point is introduced. If at least one of the infected variables is not zero, then system (6) Table 3.

The optimal control problem
In this section, we present the fractional optimal control in transmission dynamics of WNV model (6) with a time delay in the state variables M i (t), S(t), and H(t). Two control functions are considered u 1 and u 2 . The control u 1 represents the level of larvacide and adulticide used for mosquito control administered at mosquito breeding sites. The control u 2 measures the level of successful prevention (personal protection) efforts. Hence, we introduce a time delay d u in the control u 2 that represents the time required for the WNV prevention. The resulting model of WNV is given by the following system of nonlinear ordinary delay fractional-order differential equations: We consider that d u = d. Let us denote the state and control variables of the control system (16) by X = (M s , M i , B s , B i , S, E, I, H, R) ∈ R 9 and U = (u 1 , u 2 ) ∈ R 2 . We consider the following objective function which is quadratic in the control variable: The optimal control problem is defined by determining a control function U = (u 1 , u 2 ) ∈ R 2 that minimizes the cost functional J * (X, U) subject to the constraints (16). To get the necessary optimality conditions of the optimal control problem described by time delay in state and control variables, see [64]. To get the Hamiltonian for the delayed fractional optimal control problem, let us introduce the delayed state variables The Hamiltonian for the objective functional (17) and the control system (16) is given by where the adjoint variable is λ * = (λ * 1 , λ * 2 , λ * 3 , λ * 4 , λ * 5 , λ * 6 , λ * 7 , λ * 8 , λ * 9 ) ∈ R 9 . The adjoint equations are given by where 0, otherwise.
The equations for λ * 2 (t), λ * 5 (t), and λ * 8 (t) contain the advanced time t + d. Therefore, the Lagrange multipliers can be written as follows: Thus, the terminal state is free, the transversality conditions are and the characterization of the control variable u 1 is To characterize the control variable u 2 , we introduce the following switching function: Thus, the maximum condition for the control variable u 2 (t) is equivalent to the maximization ð 2 (t)u 2 (t) = max 0≤u 2 ≤1 ð 2 (t)u 2 (t).

NSFD for fractional time-delay West Nile virus model
In this section, the discretization of fractional delay WNV model is presented. We apply the NSFD method with Grünwald-Letnikov definition, so the discretization can be written as follows:

Figure 2
The approximate solutions of the fractional delay optimal control problem by using SFDM. When 0 ≤ u 1 ≤ 0.8, 0 ≤ u 2 ≤ 0.5, with d = 10 and α = 1, and t = 5 the discretized system (24) can be written as follows: Figure 4 The approximate solutions of the fractional delay WNV model without control (i.e., u 1 = u 2 = 0) by using NSFDM. With d = 1 at different α, and t = 1 N n H = S n + E n + I n + H n + R n .  (17), (i.e., the minimization of the number of exposed and infected humans is given more importance than the reduction of the total number of mosquitoes). Respectively, we use the upper bound of 0.8 and 0.5 on u 1 and u 2 . All parameters are collected and given in Table 3. In Table 4  The main goal in this study is to reduce the infected and exposed human populations from the increase due to delay in treatment. Through the results obtained in Table 4, it can be noted that, if the time delay d is increased, the infected and exposed human populations are increased at different values of α. For example, in case of non-delay optimal control problem at d = 0, it can be seen that the infected and exposed human populations are increased if the value of α is increased, and the better results occur at α = 0.70. While in the case of delay optimal control problem, it can be noted that at d = 1, the lowest value of the infected and exposed human populations occurs at α = 0.70, similarly for d = 3, 5, 7. In Table 5 Thus, the objective function J * (X, U) obtained using NSFDM is better than that obtained using SFDM. It can be concluded that the numerical treatments for the effec- Figure 6 The approximate solutions of the fractional delay optimal control problem by using NSFDM. When 0 ≤ u 1 ≤ 0.8, 0 ≤ u 2 ≤ 0.5, with d = 5 and different α, and t = 0.1 tiveness of nonstandard finite difference technique for the fractional time-delay optimal control problem are better than the numerical treatments for the time-delay optimal control.

Numerical experiment
Also, from Table 5 we see that NSFDM helps us to achieve the goal of this research and to obtain the lowest value of the infected and exposed human populations better than SFDM, where the lowest value of these populations occurs at α = 0.90 and d = 3, 5, 7.
The numerical simulations of the optimal control problem (16) and (17) using two numerical methods, NSFDM and SFDM, are shown in Figs. 1 and 2. Figure 1 presents the comparison between these proposed methods at d = 10 and α = 1, and t = 1. It can be observed that the infected human population I(t) is about 280 and 320 for NSFDM and  Fig. 4, the approximate solutions of the fractional delay WNV model (6) without control (i.e., u 1 = u 2 = 0) with d = 1 at different α, and t = 1 by using NSFDM are shown. It can be seen that the infected human population I(t) is about 750, 650, and 400 for α = 1, α = 0.98, and α = 0.90, respectively. Numerical comparisons between the approximate solutions of the fractional delay optimal control problem (16) and (17) with and without control (i.e., u 1 = u 2 = 0) with d = 1 at α = 0.90, and t = 1 by using NSFDM, are introduced in Fig. 5. It can be shown that the infected human population I(t) is about 250 and 450 in two cases with and without control, respectively. The behavior of the approximate solutions of the fractional delay optimal control problem (16) and (17) at different time delay d and different α at t = 0.1 is included in Fig. 6. We notice that the delayed infected human population I(t -d) is about 589.0570, 507.9876, and 278.5177 for α = 1, α = 0.98, and α = 0.90, respectively, when d = 5. Figure 7 shows numerical comparisons between the approximate solutions of the optimal control problem (16) and (17) in two cases, the fractional delay optimal control problem at d = 5, and the fraction non-delay optimal control problem at d = 0 when α = 0.98. The infected human population I(t) is about 650 and 400 in two cases with and without delay, respectively. Figure 8 shows the relationship between state  H(t -d) and H(t), and the control variable u 2 (t -d) and u 2 (t) when d = 1 for α = 0.98. Figure 9 describes the approximate solutions of the fractional delay optimal control problem at the same time delay d = 1, d = 5, and d = 10 in case of fractional order at α = 0.98. The infected human population I(t) is about 400, 500, and 650 for d = 1, d = 5, and d = 10, respectively.

Conclusions
In this paper, fractional time-delay optimal control in transmission dynamics of WNV model is presented. Time delays in both state and control variables are formulated and studied. NSFDM and SFDM are introduced to study the behavior of the proposed model. According to the results implemented in this paper, it can be observed that the numerical Table 4 Comparisons between the value of objective function J * (X, U) with and without control cases and E(T f ) + I(T f ) when 0 ≤ u 1 ≤ 0.8, 0 ≤ u 2 ≤ 0.5, at different time delay d and different α with t = 0.1, T f = 100 by using NSFDM d α J * (X, U) with control E(T f ) + I(T f ) J * (X, U) without control (i.e., u 1 = u 2 = 0) 0   Tables 4 and 5, it is found that the values of the objective functional obtained by NSFDM are better than the results obtained by SFDM, it can be applied simply and effectively to solve optimal control of fractional delay.