Adv Syst Sci Appl 2021; 02:83–103 Published online at https://ijassa.ipu.ru. Stability Analysis and Optimal Measure for Controlling Eco-Epidemiological Dynamics of Prey-Predator Model Jamiu A. Ademosu, Samson Olaniyi*, Sunday O. Adewale Department of Pure and Applied Mathematics, Ladoke Akintola University of Technology, Ogbomoso, Nigeria Abstract: An eco-epidemiological model representing the interactions between prey and predator populations affected by a disease in an ecosystem is presented. The model is governed by a five-dimensional nonlinear system of ordinary differential equations coupling both ecological and epidemiological features of interacting populations. The well-posedness of the model is established with respect to positivity and boundedness of solutions. Conditions for asymptotic stability of different equilibrium points are extensively investigated to determine the existence and coexistence of prey and predator species using local linearization and Lyapunov functions techniques. Additionally, the analysis of the model is extended to assess the effects of three time- dependent control functions, such as disease prevention, treatment and alternative resource for predator, on the population dynamics of the prey-predator coexistence in the system. Keywords: ecology, epidemiology, prey-predator model, stability, optimal control problem, simulations 1. INTRODUCTION Interaction between organisms of different species is a common phenomenon in an ecosystem. It is an intrinsic feature which cannot be downplayed in the study of ecological system. This interaction, especially between prey and predator species, often affects or changes the population sizes of organisms in the system [11]. Since the advent of the classical Lotka and Volterra prey-predator models [20, 31], many ecosystem models, which can be defined as mathematical representations of ecological system of prey and predator interactions, have been studied in literature to better understand the real system (see, e.g. [5, 11, 13, 17, 27]). In [11], analysis of modified Lotka-Volterra model was carried out by taking into consideration an environmental case containing two related populations of prey and predator species having influence on the size of each other. Kumar and Kumari [17] studied the impact of fear on a chaotic model describing the interactions of the species in the food chain of one prey and two predators. Prasad et al. [27] investigated the role of mutual interference between predators on the dynamics of additional food provided predator-prey system. In recent decades, a number of studies have considered ecosystems of interacting populations among which diseases spread. The resulting ecological systems are described by eco-epidemiological models [7–9, 21, 29]. Das [9] studied the predator-prey model with parasitic infection transmitting among the predator population only. Nandi et al. [21] developed a predator-prey system where prey species only are infected with a viral disease resulting into susceptible and two-stage infected classes. Bera et al. [8] presented theoretical analysis supported by simulation of the dynamical behaviors of a prey-predator system where both populations are affected by diseases. ∗Corresponding author: solaniyi@lautech.edu.ng 84 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE Interest has been growing recently in the application of optimal control theory to the dynamics of eco-epidemiological models for prey-predator system. AL-Nassir [6] studied a two-dimensional continuous prey-predator system incorporating an optimal control variable with the aim of reducing the number of the predator density in the population. In [10], the authors extended a prey-predator mathematical model with infection and harvesting on prey to include one time-dependent control for preventing infection in the prey population only. In another development, Hugo et al. [16] applied optimal control theory in minimizing the spread of Newcastle disease among chickens (prey) and humans (predator) populations using three control variables, such as vaccination of prey, education campaign and treatment of infected humans. In this work, a five-dimensional nonlinear prey-predator model is studied to gain further insights into the dynamics of eco-epidemiological interactions between prey and predator populations. The analysis of the model is extended to investigate the impact of three optimal control variables, namely disease prevention in both prey and predator populations, treatment control for prey, and provision of alternative resources control for predator population. The rest of the work is organized as follows: The model is formulated and its well-posedness is established in Section 2. In Section 3, existence and stability of possible equilibrium points are investigated. In Section 4, the model is extended to explore optimal control dynamics with numerical simulations. Section 5 is devoted to concluding remarks. 2. MODEL FORMULATION The dynamic eco-epidemiological prey-predator system is formulated by sub-dividing the prey population at time, t, into susceptible Xs(t), infected Xi(t) and recovered Xr(t). On the other hand, the predator population is sub-divided into two, namely susceptible Ys(t) and infected Yi(t). It is assumed that the susceptible prey population increases logistically with intrinsic growth rate r and carrying capacity k. Following effective contact with the infected prey, the susceptible prey becomes infected at rate β1. To avoid extinction of the prey population due to disease, the infected prey is allowed to recover at per capita rate γ. The natural mortality rate for the prey population is given by µ1. Further, it is assumed that the susceptible predator, being healthy, has the capacity to attack and consume both susceptible and infected prey, while the recovered prey is assumed to be under cover from predation. It is also assumed that predators depend on prey for sustenance. The susceptible predator contracts disease from infected predator at rate β2, while the infected predator, being unhealthy, has the capacity to attack and consume infected prey only as assumed in [8]. It has been assumed that only infected prey population recovers due to their accessibility to treatment. Hence, recovered class for predators is not considered. This assumption is based on some eco-epidemiological cases, for examples: carnivorous animals (predator) and humans (prey), and hawks (predator) and chicks (prey). The predation and conversion rates for the predator population are given, respectively, by Pi and Ci, for i = 1, 2, 3. The natural mortality rate for predator population is given by µ2. Concise description of the variables and the parameters of the prey-predator system is presented in Table 2.1 and the nonlinear ordinary differential equations describing the interacting prey and predator populations is given by Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 85 dXs dt =rXs ( 1− Xs k ) − β1XsXi − P1XsYs − µ1Xs dXi dt =β1XsXi − P2XiYs − P3XiYi − (µ1 + γ)Xi dXr dt =γXi − µ1Xr dYs dt =C1P1XsYs + C2P2XiYs − β2YsYi − µ2Ys dYi dt =β2YsYi + C3P3XiYi − µ2Yi. (2.1) Table 2.1. The variables and parameters of the malaria model (2.1) Variable Description Xs(t) Susceptible (healthy) prey Xi(t) Infected (unhealthy) prey Xr(t) Recovered prey Ys(t) Susceptible (healthy) predator Yi(t) Infected (unhealthy) predator Parameter Description r Intrinsic growth rate k Carrying capacity β1 Disease transmission rate in prey β2 Disease transmission rate in predator Pi, i = 1, 2, 3 Predation coefficients Ci, i = 1, 2, 3 Conversion rates γ Recovery rate of infected prey µ1 Natural mortality rate of prey µ2 Natural mortality rate of predator 2.1. Well-posedness of the model The mathematical and eco-epidemiological relevance of the prey-predator system (2.1) depends on the well-posedness of the model. Keeping in mind that all the parameters of the model are non-negative. Here, the boundedness and positivity of solutions of the model are investigated to establish the well-posedness of the model. 2.1.1. Boundedness of solutions The following result is required to establish the boundedness of solutions of the model (2.1). Lemma 2.1: The susceptible prey population Xs(t) is bounded. Proof It is clear from the first equation of the model (2.1) that dXs dt ≤ rXs − X2 s k . (2.2) Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 86 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE Solving the nonlinear first order differential inequality of Bernoulli type (2.2) yields Xs(t) ≤ kXs(0) ke−rt +Xs(0)(1− e−rt) . It follows that the lim supXs(t) < k as t→∞. Theorem 2.1: The solutions of the system (2.1) are uniformly bounded. Proof Let the total prey and predator populations be represented, respectively, by X and Y , so that the total population of both species be T = X + Y . Then it follows that dT dt ≤ rXs − µ1X − µ2Y. (2.3) Using Lemma 2.1 and let µ = min{µ1, µ2}, then (2.3) becomes dT dt + µT ≤ rXs ≤ rk, which on integration gives T (t) ≤ rk µ (1− e−µt) + T (0)e−µt. Consequently, lim supT (t) < (rk/µ) + ε, ∀ ε > 0 as t→∞. This ends the proof. 2.1.2. Positivity of solutions Theorem 2.2: The solution set {Xs, Xi, Xr, Ys, Yi} of the eco-epidemiological model (2.1) with non- negative initial conditions, Xs(0), Xi(0), Xr(0), Ys(0), Yi(0) in Ω, remain non-negative for all time t > 0. Proof SinceXs(t) < k as shown in Lemma 2.1, then (1−Xs/k) ≥ 0 with equality if the size of the susceptible prey is exactly its carrying capacity. Hence, the first equation of the model (2.1) implies that dXs(t) dt + (β1Xi + P1Ys + µ1)Xs(t) ≥ 0, (2.4) which on using integrating factor gives d dt [ Xs(t) exp (∫ t 0 (β1Xi(φ) + P1Ys(φ))dφ+ µ1t )] ≥ 0. (2.5) Further integration of (2.5) gives Xs(t) ≥ Xs(0) exp [ − (∫ t 0 (β1Xi(φ) + P1Ys(φ))dφ+ µ1t )] > 0, ∀t > 0. (2.6) All other variables can be proved to be positive in a similar approach. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 87 Hence, it is sufficient to consider the dynamics of the flow generated by the system (2.1) in a feasible region defined by Ω = { (Xs, Xi, Xr, Ys, Yi) ∈ R5 + : T < rk µ + ε } . In the region Ω, the prey-predator model (2.1) can be said to be eco-epidemiologically well- posed. 3. ANALYSIS OF THE MODEL In this section, the prey-predator model is analysed around the possible equilibrium points. 3.1. Existence of equilibrium points The prey-predator model (2.1) has the following possible equilibrium points: 3.1.1. Trivial equilibrium point (E0) This is a steady state in the absence of prey and predator populations. Hence, no interactions exist. The equilibrium point is given by E0 = (0, 0, 0, 0, 0). (3.7) 3.1.2. Axial equilibrium point (E1) This is a steady state where only the susceptible prey is present. In other words, it is the disease-free cum predator-free equilibrium point. It exists provided that the intrinsic growth rate r exceeds the natural mortality rate µ1. The equilibrium point is given by E1 = ( k [ 1− µ1 r ] , 0, 0, 0, 0 ) . (3.8) 3.1.3. Disease-free equilibrium point (E2) This is a steady state where there is no disease in both prey and predator populations. It exists provided the inequality r ( 1− µ2 kC1P1 ) > µ1 holds. Hence, it is obtained as E2 = ( µ2 C1P1 , 0, 0, 1 P1 [ r ( 1− µ2 kC1P1 ) − µ1 ] , 0 ) . (3.9) 3.1.4. Predator-free equilibrium point (E3) This is a steady state where there is no predator in the system. It exists provided the inequality r ( 1− ( µ1+γ kβ1 )) > µ1 holds. The equilibrium is given as E3 = ( µ1 + γ β1 , 1 β1 [ r ( 1− ( µ1 + γ kβ1 )) − µ1 ] , γ β1µ1 [ r ( 1− ( µ1 + γ kβ1 )) − µ1 ] , 0, 0 ) . (3.10) 3.1.5. Infected prey-free equilibrium (E4) This is a steady state where the infected prey does not exist. The equilibrium will exist if k > k r ( P1µ2 β2 + µ1 ) and C1P1 β2 [ k − k r ( P1µ2 β2 + µ1 )] > µ2. Thus, it is given by E4 = ( k − k r ( P1µ2 β2 + µ1 ) , 0, 0, µ2 β2 , C1P1 β2 [ k − k r ( P1µ2 β2 + µ1 )] − µ2 ) . (3.11) Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 88 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE 3.1.6. Infected predator-free equilibrium (E5) This is a steady state where the infected predator is absent in the ecosystem. The equilibrium is given by E5 = (X∗ s , X ∗ i , X ∗ r , Y ∗ s , 0), (3.12) such that X∗ s > 0, X∗ i > 0, X∗ r > 0, Y ∗ s > 0, and where X∗ s = C2[P2(r − µ1)− P1(µ1 + γ)]− β1µ2 C2P2 ( 1 k + P1β1 P2 ) − C1P1β1 , X∗ i = µ2 C2P2 − C1P1 C2P2 C2[P2(r − µ1)− P1(µ1 + γ)]− β1µ2 C2P2 ( 1 k + P1β1 P2 ) − C1P1β1  , X∗ r = γX∗ i µ1 , Y ∗ s = β1X ∗ s − (µ1 + γ) P2 . 3.1.7. Interior equilibrium point (E6) This is a steady state where all the populations for prey and predator species are present in the system. The equilibrium is given by E6 = (X∗∗ s , X ∗∗ i , X ∗∗ r , Y ∗∗ s , Y ∗∗ i ), (3.13) such that X∗∗ s > 0, X∗∗ i > 0, X∗∗ r > 0, Y ∗∗ s > 0, Y ∗∗ i > 0, and where X∗∗ s = 1 k [ (r − µ1)−X∗∗ i ( β1 − C3P1C3 β2 )] , X∗∗ i = β2[(P2 − P3)µ2 + β2(µ1 + γ)− (r − µ1)(β1β2 − C1P1P3)] β2(P2P3C3 + C2P2P3) + (β1β2 − C1P1P3)(β1β2 − C3P1P3) , X∗∗ r = γX∗∗ i µ1 , Y ∗∗ s = µ2 − C3P3X ∗∗ i β2 , Y ∗∗ i = β1X ∗∗ s − P2Y ∗∗ s − (µ1 + γ) P3 . 3.2. Stability analysis The Jacobian matrix of the prey-predator system (2.1) is given by J =  J1 −β1Xs 0 −P1Xs 0 β1Xi J2 0 −P2Xi −P3Xi 0 γ −µ1 0 0 C1P1Ys C2P2Ys 0 J3 −β2Ys 0 0 0 β2Yi β2Ys + C3P3Xi − µ2  , (3.14) Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 89 where J1 = r − 2rXs k − β1Xi − P1Ys − µ1, J2 = β1Xs − P2Ys − P3Yi − (µ1 + γ), J3 = C1P1Xs + C2P2Xi − β2Yi − µ2. The stability of each of the equilibrium points is analysed by finding the eigenvalues of the Jacobian matrix (3.14) evaluated at each point. 3.2.1. Stability of E0 The Jacobian matrix (3.14) is evaluated at the trivial equilibrium point (3.7). Hence, solving the corresponding characteristic equation |J(E0)− λI5)| = 0, where λ is the eigenvalue and I5 is the identity matrix of order five, the following eigenvalues are obtained: λ1 = r − µ1, λ2 = −µ1, λ3 = −µ2, λ4 = −(µ1 + γ), λ5 = −µ2. It can be seen that all the eigenvalues are unconditionally negative except λ1. Thus, E0 will be stable if r < µ1 and unstable if r > µ1. This result is theorized as follows. Theorem 3.1: The trivial equilibrium point E0 of the prey-predator model (2.1), given by (3.7), is stable if r < µ1 and unstable if r > µ1. The implication of Theorem 3.1 is that the coexistence of prey-predator population will not occur if the intrinsic growth rate r does not exceed the natural mortality rate µ1. In other words, stable equilibrium point E0 will result into extinction of prey and predator species in the ecosystem. 3.2.2. Stability of E1 The characteristic equation of the matrix (3.14) evaluated at the axial equilibrium point E1 is given by |J(E1)− λI5)| = 0, which gives the following eigenvalues: λ1 = µ1 − r, λ2 = −µ1, λ3 = β1k ( 1− µ1 r ) − (µ1 + γ), λ4 = C1P1k ( 1− µ1 r ) − µ2, λ5 = −µ2. The stability result is theorized as follows. Theorem 3.2: The axial equilibrium point E1 of the prey-predator model (2.1), given by (3.8), is stable if the inequalities: µ1 < r, β1k ( 1− µ1 r ) < (µ1 + γ) and C1P1k ( 1− µ1 r ) < µ2 hold. The implication of Theorem 3.2 is that the stable E1 will support the existence of prey population only in the ecosystem provided these threshold conditions are satisfied: R01 < 1,R02 < 1 andR03 < 1, where R01 = µ1 r , R02 = β1k ( 1− µ1 r ) µ1 + γ , and R03 = C1P1k ( 1− µ1 r ) µ2 . Now, suppose R0 = max{R01,R02,R03}, and considering the prey sub-model only such that the axial equilibrium point E1 becomes Ẽ1 = (X0 s , 0, 0), where X0 s = k [ 1− µ1 r ] . The following global stability result is proved. Theorem 3.3: The equilibrium point Ẽ1 is globally asymptotically stable ifR0 < 1. Proof Since P1 = 0 in the absence of predation, thenR0 = max{R01,R02}. Consider the Lyapunov Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 90 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE function (combination of Goh-Volterra [2, 3, 13, 15, 22] and linear types) given by L = Xs −X0 s −X0 s ln Xs X0 s +Xi. (3.15) The time derivative of Lyapunov function (3.15) along the trajectory of the prey only sub- model is given by L̇ = ( 1− X0 s Xs ) dXs dt + dXi dt . (3.16) Using dXs dt = rXs ( 1− Xs k ) − β1XsXi − µ1Xs and dXi dt = β1XsXi − (µ1 + γ)Xi in (3.16) with Xs ≤ X0 s , it follows that L̇ = (Xs −X0 s ) ( r [ 1− Xs k ] − β1Xi − µ1 ) + β1XsXi − (µ1 + γ)Xi ≤ r [ 1− X0 s k ] − β1Xi − µ1 + β1X 0 sXi − (µ1 + γ)Xi. (3.17) Since r [ 1− X0 s k ] = µ1 at the disease-free steady state of the prey only sub-model, then it follows from (3.17), by further simplifications, that L̇ ≤ −β1Xi − (µ1 + γ) ( 1− β1X 0 s µ1 + γ ) Xi. (3.18) SinceR02 = β1X 0 s µ1 + γ , consequently (3.18) becomes L̇ ≤ −β1Xi − (µ1 + γ) (1−R02)Xi. Therefore, L̇ ≤ 0 when R02 < 1 with equality, L̇ = 0, if and only if Xi = 0. It follows from LaSalle’s Invariance Principle [19] that the largest invariant set contained in {(Xs, Xi, Xr) ∈ R3 + : L̇ = 0} is the singleton set {Ẽ1}. Hence, the equilibrium point Ẽ1 = (k[1− (µ1/r)], 0, 0) is globally asymptotically stable. The implication of Theorem 3.3 is that the presence of disease in the prey population can be eliminated if R02 < 1, irrespective of the initial sizes of the infected prey in the ecosystem. Thus, every solution of the prey only sub-model will always converge to the point Ẽ1 wheneverR02 is less than unity. This theoretical result is demonstrated in Figure 3.1. 3.2.3. Stability of E2 Considering that there are only two disease classes in the prey-predator system (2.1), then using the notations in the next generation operator method [30], it follows that F = [ β1µ2 C1P1 0 0 β2 P1 m1 ] and V = [ P2 P1 m1 + µ1 + γ 0 0 µ2 ] , where m1 = r ( 1− µ2 kC1P1 ) − µ1. Hence, the basic reproduction number of the full system (2.1), denoted by R∗ 0, is given by R∗ 0 = ρ(FV −1) = max  β1µ2 C1P1 ( P2 P1 m1 + µ1 + γ ) , β2m1 P1µ2  . (3.19) Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 91 0 5 10 15 0 10 20 30 40 50 60 70 Time In fe ct ed P re y Fig. 3.1. Convergence of solution trajectories of prey sub-model to the equilibrium point Ẽ1 at different initial sizes. The parameter values used are: r = 11.2, β1 = 0.04, µ1 = 0.85, γ = 0.75, k = 30, so thatR01 = 0.0759 andR02 = 0.6931, implying thatR0 = R02 < 1. Therefore, by Theorem 2 of [30], the following stability result is claimed. Theorem 3.4: The disease-free equilibrium point E2 of the prey-predator model (2.1), given by (3.9), is stable provided that R∗ 0 < 1. The implication of Theorem 3.4 is that the presence of few species around the equilibrium point E2 will support the coexistence of both prey and predator species in a disease-free ecosystem whenever the basic reproduction number, R∗ 0, is less than unity. The global dynamics of the prey-predator around the disease-free equilibrium point E2 is explored next. Theorem 3.5: The disease-free equilibrium point of the prey-predator model (2.1), given by (3.9), is globally asymptotically stable if R∗ 0 > 1. Proof The proof is based on using the Lyapunov function defined by F = P1 P2m1 + P1(µ1 + γ) Xi + 1 µ2 Yi, (3.20) which when differentiated with respect to time, t, gives dF dt = P1 P2m1 + P1(µ1 + γ) [β1XsXi − P2XiYs − P3XiYi − (µ1 + γ)Xi] + 1 µ2 (β2YsYi + C3P3XiYi − µ2Yi). Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 92 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE Since Xs ≤ β1µ2 C1P1 and Ys ≤ 1 P1 [ r ( 1− µ2 kC1P1 ) − µ1 ] , noting that m1 = r ( 1− µ2 kC1P1 ) − µ1. Therefore, it follows, after simplification, that dF dt ≤  β1µ2 C1P1 ( P2 P1 m1 + µ1 + γ ) − 1 Xi − P1P3 P2m1 + P1(µ1 + γ) XiYi + [ β2m1 P1µ2 − 1 ] Yi + 1 µ2 (C3P3XiYi) ≤ max  β1µ2 C1P1 ( P2 P1 m1 + µ1 + γ ) , β2m1 P1µ2 − 1  (Xi + Yi) − [ P1 P2m1 + P1(µ1 + γ) − C3 µ2 ] P3XiYi = [R∗ 0 − 1](Xi + Yi)− [ P1 P2m1 + P1(µ1 + γ) − C3 µ2 ] P3XiYi. Consequently, dF/dt ≤ 0 if R∗ 0 < 1 and Xi = 0 or Yi = 0. Thus, Xs → µ2/(C1P1), Xr → 0 and Ys → m1/P1 as both Xi → 0 and Yi → 0. Therefore, by LaSalle’s invariance principle [19], the disease-free equilibrium point E2 is globally asymptotically stable. The implication of Theorem 3.5 is that elimination of disease in the ecosystem is independent of initial sizes of infected prey and infected predator whenever R∗ 0 < 1. This result is demonstrated in Figure 3.2(a) for infected prey and Figure 3.2(b) for infected predator population, where every solution trajectory converges to the disease-free equilibrium. 0 5 10 15 20 25 30 35 40 0 1 2 3 4 5 6 7 8 9 10 Time In fe c te d P re y (a) 0 5 10 15 0 5 10 15 20 25 Time In fe c te d P re d a to r (b) Fig. 3.2. Convergence of solution trajectories of prey-predator model to the disease-free equilibrium point E2 at different initial sizes. The parameter values used are: mu1 = 0.8, r = 3, β1 = 0.4, P1 = 1.0, γ = 0.35, P2 = 0.2, P3 = 1.0, µ2 = 0.85, k = 30, β2 = 0.3, C1 = 0.25, C2 = 0.25, C3 = 0.15, so that R∗ 0 = 0.6565 < 1. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 93 3.2.4. Stability of E3 Evaluating the Jacobian matrix (2.1) at predator-free equilibriumE3 leads to the characteristic equation |J(E3)− λI5| = 0. Thus, the following eigenvalues are obtained: λ1 = −µ1 λ2 = −µ2, λ3 = C1P1 β1 (µ1 + γ)− C2P2 β1 (m2 + µ2), λ4 = −r 2 ( µ1 + γ kβ1 ) + 1 2 √ r2 ( µ1 + γ kβ1 )2 + 4m2(µ1 + γ), λ5 = −r 2 ( µ1 + γ kβ1 ) − 1 2 √ r2 ( µ1 + γ kβ1 )2 + 4m2(µ1 + γ), where m2 = r ( 1− µ1+γ kβ1 ) − µ1. The stability of E3 is theorized as follows Theorem 3.6: The predator-free equilibrium point E3 of the prey-predator model (2.1), given by (3.10), is stable if the following conditions hold C1P1(µ1 + γ) < C2P2(m2 + µ2) and √ r2 ( µ1 + γ kβ1 )2 + 4m2(µ1 + γ) < r ( µ1 + γ kβ1 ) . This result implies that stable E3 will not be able to support the coexistence of both prey and predator populations. The presence of few species around the predator-free equilibrium point E3 will make predator population to go extinct in the ecosystem. 3.2.5. Stability of E4 The Jacobian matrix (3.14) is evaluated at the infected prey-free equilibrium point (3.11). Hence, solving the corresponding characteristic equation |J(E4)− λI5)| = 0 gives one eigenvalue, λ1 = −µ1, and the remaining four eigenvalues are obtainable from λ4 +W1λ 3 +W2λ 2 +W3λ+W4 = 0, (3.21) where W1 = r − µ1 − P1 µ2 β2 , W2 = d1d2 − 2µ2(d1 + d2) + µ2 + C1P1d3 − C1P1 µ2 β2 with d1 = r − 2r d3 k − P1 µ2 β2 − µ1, d2 = β1d3 − P2 µ2 β2 − C1P1P3 d3 β2 − (µ1 + µ2 + γ) and d3 = k − k r ( P1µ2 β2 + µ2 ) . W3 = C1P1 [ (d2 − µ2) µ2 β2 − d3(d2µ2 + d1) ] − d2µ2(2d1 + µ2) and W4 = d2µ 2 2 ( d1 + P1 β2 ) . Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 94 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE Applying Routh-Hurwitz’s criterion, the roots of (3.21) will have negative real parts if W1 > 0, W1W2 > W3, W1W2W3 > W 2 1W4 +W 2 3 . The stability result is hereby theorized as follows. Theorem 3.7: The infected prey-free equilibrium of the prey-predator model (2.1), given by (3.11), is stable if W1 > 0, W1W2 > W3 and W1W2W3 > W 2 1W4 +W 2 3 . The implication of Theorem 3.7 is that the presence of few species around the stable infected prey-free equilibrium E4 will support the coexistence of both prey and predator populations in the ecosystem. 3.2.6. Stability of E5 Two of the eigenvalues of the Jacobian matrix (3.14) evaluated at the infected predator-free equilibrium point E5 are given by λ1 = −µ1 and λ2 = −µ2, while the remaining can be obtained from the characteristic equation λ3 + e1λ 2 + e2λ+ e3 = 0, (3.22) where e1 = −(f1 + f2 + f3), with f1 = 2C1P1f̃1 − µ2(1 + C2P2) such that f̃1 = C2P2(r − µ1)− P1C2(µ1 + γ)− βµ2 C2P2 ( 1 k + P1β1 P2 ) − C1P1β1 , f2 = r − µ2 − [ 2r k + C2P2 + kP1β1(C2 − C1) k(C2P2 + kP1C2 − C1P1β1) ] f̃1 − β1 ( µ2 − C1P1 C2P2 f̃1 ) (1− P1) f3 = β1f̃1 − P2 [( C2P2 + kP1β1(C2 − C1) k(C2P2 + kP1C2 − C1P1β1) ) f̃1 − β1 ( µ2 − C1P1 C2P2 f̃1 )] − (µ1 + γ). e2 = f1f2 + f3(f1 + f2) + β2 1 f̃1 ( µ2 − C1P1 C2P2 f̃1 ) + P 2 1 C1 P2 f̃1(f3 + (µ1 + γ)− β1f̃1), e3 = 1 P2 (f3 + (µ1 + γ)− β1f̃1) [( µ2 − C1P1 C2P2 f̃1 ) (P2(C2[P1f̃1 − P2]− C1P1))− C1P 2 1 f̃1(f3 + f1f2) ] . By Routh-Hurwitz’s criterion, the roots of the characteristic equation (3.22) will have negative real parts if e1 > 0 and e1e2 > e3. This result is stated as follows. Theorem 3.8: The infected predator-free equilibrium point of the prey-predator system (2.1), given by (3.12), is stable if e1 > 0 and e1e2 > e3. The result in Theorem 3.8 implies that coexistence of prey and predator species in the ecosystem is possible in the absence of infected predator, provided that the given parametric conditions are valid. 3.2.7. Stability of E6 Recall from (3.19) that the basic reproduction number of the full system (2.1) is given by R∗ 0 = max  β1µ2 C1P1 ( P2 P1 m1 + µ1 + γ ) , β2m1 P1µ2  . The following global stability result is explored. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 95 Theorem 3.9: The interior equilibrium point E6 of the prey-predator model (2.1) is globally asymptotically stable if R∗ 0 > 1. Proof Let R∗ 0 > 1, such that the interior equilibrium point E6 exists. Noting from (3.13) that E6 = (X∗∗ s , X ∗∗ i , X ∗∗ r , Y ∗∗ s , Y ∗∗ i ), then the following Lyapunov function of quadratic type [28] is defined: V = 1 2 [(Xs −X∗∗ s ) + (Xi −X∗∗ i ) + (Xr −X∗∗ r ) + (Ys − Y ∗∗ s ) + (Yi − Y ∗∗ i )]2 (3.23) The time derivative of V in (3.23) along the solution path of the full system (2.1) is given by V̇ = [(Xs −X∗∗ s ) + (Xi −X∗∗ i ) + (Xr −X∗∗ r ) + (Ys − Y ∗∗ s ) + (Yi − Y ∗∗ i )] dT dt , (3.24) where T (t) = Xs +Xi +Xr + Ys + Yi. Recall from (2.3) of the boundedness result (Theorem 2.1) that dT dt ≤ rk − µT, and since X∗∗ s +X∗∗ i +X∗∗ r + Y ∗∗ s + Y ∗∗ i = rk µ . It follows from (3.24) that V̇ ≤ [ T (t)− rk µ ] [rk − µT (t)] = µ [ T (t)− rk µ ] [ rk µ − T (t) ] = −µ [ T (t)− rk µ ]2 = −µ[(Xs −X∗∗ s ) + (Xi −X∗∗ i ) + (Xr −X∗∗ r ) + (Ys − Y ∗∗ s ) + (Yi − Y ∗∗ i )]2. It can be seen that V̇ ≤ 0 with equality if and only if Xs = X∗∗ s , Xi = X∗∗ i , Xr = X∗∗ r , Ys = Y ∗∗ s and Yi = Y ∗∗ i . Hence, by LaSalle’s invariance principle [19], the largest invariant set in {(Xs, Xi, Xr, Ys, Yi) ∈ R5 + : V̇ = 0} is the only set {E6}. Thus, the interior equilibrium point E6 = (X∗∗ s , X ∗∗ i , X ∗∗ r , Y ∗∗ s , Y ∗∗ i ) is globally asymptotically stable. This completes the proof. The implication of Theorem 3.9 is that coexistence of prey and predator species in the ecosystem is possible while disease persists in the population. The presence of disease in the ecosystem, such that R∗ 0 > 1, is independent of initial population sizes of both infected prey and predator species. This global stability result is depicted in Figure 3.3(a) for infected prey population and Figure 3.3(b) for infected predator population. The phase portrait which illustrates the stable coexistence of prey and predator populations is presented in Figure 3.4. 4. OPTIMAL CONTROL MODEL In order to enhance the coexistence of both prey and predator species in the ecosystem, the prey-predator model (2.1) is extended by incorporating the following three time-dependent optimal control functions: Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 96 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE 0 5 10 15 20 25 30 35 40 0 2 4 6 8 10 12 14 16 18 20 Time In fe c te d P re y (a) 0 5 10 15 0 5 10 15 20 25 Time In fe c te d P re d a to r (b) Fig. 3.3. Convergence of solution trajectories of prey-predator model to the interior equilibrium point E6 at different initial sizes. The parameter values used are:mu1 = 0.5, r = 11.2, β1 = 0.4, P1 = 1.0, γ = 0.35, P2 = 0.2, P3 = 1.0, µ2 = 0.65, k = 30, β2 = 0.65, C1 = 0.25, C2 = 0.25, C3 = 0.15, so that R∗ 0 = 9.7293 > 1. 0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 Prey P re da to r Fig. 3.4. Phase portrait of predator against prey showing stable coexistence of both populations around the interior (positive) equilibrium point. 1. The control function u1(t) ∈ [0, 1], which represents disease prevention measure for both prey and predator populations. 2. The function u2(t) ∈ [0, 1], which represents treatment control for prey population. 3. The control function u3(t) ∈ [0, 1] representing alternative source of food for predator survival. Therefore, the prey-predator model (2.1) becomes a non-autonomous system of ordinary differential equations of the form: Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 97 dXs dt = rXs ( 1− Xs k ) − (1− u1(t))β1XsXi − P1XsYs − µ1Xs dXi dt = β1XsXi − P2XiYs − P3XiYi − (µ1 + γ)Xi − u2(t)r0Xi dXr dt = γXi + u2(t)r0Xi − µ1Xr dYs dt = C1P1XsYs + C2P2XiYs − (1− u1(t))β2YsYi − µ2Ys + u3(t)c0Ys dYi dt = (1− u1(t))β2YsYi + C3P3XiYi − µ2Yi, (4.25) where r0 and c0 are the rate constants for the optimal controls u2(t) and u3(t), respectively. Since the goal is to reduce the risk of both prey and predator species going extinct, such that the populations of infected prey and predator are minimized at minimum costs possible, then the following objective functional is formed J = ∫ tf 0 ( A1Xi + A2Yi + 1 2 [B1u 2 1 +B2u 2 2 +B3u 2 3] ) dt, (4.26) where A1, A2, B1, B2 and B3 are the positive weight constants for balancing the terms in the objective functional over the time interval [0, tf ]. Further, the cost function associated with the disease prevention measure for prey and predator populations is represented by the term 1/2B1u 2 1, while the treatment measure for prey population and the provision of alternative food source for predator are, respectively, represented by the terms 1/2B2u 2 2 and 1/2B3u 2 3. The objective functional is made of nonlinear cost on controls of quadratic type in line with the standards for optimal control problems (see, e.g. [1, 4, 14, 16, 24, 25]). Therefore, the optimal control problem for the prey-predator system (4.25) is to seek a control triple u∗ = (u∗1, u ∗ 2, u ∗ 3, ), such that J(u∗) = min {J(u1, u2, u3) : u1, u2, u3 ∈ U} , (4.27) where U = {ui(t) : 0 ≤ ui(t) ≤ 1,Lebesgue measurable, t ∈ [0, tf ], for i = 1, 2, 3} is a non- empty set of all admissible controls. 4.1. Analysis of the control model The necessary conditions that must be satisfied by the optimal control triple u∗ for the minimization problem are derivable from Pontryagin’s Maximum Principle [26]. This principle converts the problem in (4.27) with (4.26) subject to the state system (4.25) into a problem of minimizing pointwise a Hamiltonian H, with respect to u1(t), u2(t) and u3(t). Thus, the Hamiltonian for the control problem is given by Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 98 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE H = A1Xi + A2Yi + 1 2 [B1u 2 1 +B2u 2 2 +B3u 2 3] +η1 [ rXs ( 1− Xs k ) − (1− u2(t))β1XsXi − P1XsYs − µ1Xs ] +η2 [β1XsXi − P2XiYs − P3XiYi − (µ1 + γ)Xi − u2(t)r0Xi] +η3[γXi + u2(t)r0Xi − µ1Xr] +η4[C1P1XsYs + C2P2XiYs − (1− u1(t))β2YsYi − µ2Ys + u3(t)c0Ys] +η5[(1− u1(t))β2YsYi + C3P3XiYi − µ2Yi], (4.28) where η1, η2, η3, η4 and η5 are the adjoint or co-state variables corresponding to the state variables Xs, Xi, Xr, YS and Yi, respectively. Applying the existence results as given in [2, 12], the following result is obtained. Theorem 4.1: For the optimal control triple u∗ that minimizes the objective functional (4.26) over U subject to (4.25), there exist adjoint variables η1, η2, η3, η4 and η5 satisfying the system governing the adjoint variables dη1 dt = 2rXs k η1 + (1− u1)β1Xi(η1 − η2) + P1Ys(η1 − C1η4) + (µ1 − r)η1 dη2 dt = (1− u1)β1Xs(η1 − η2) + P2Ys(η2 − C2η4) + P3Yi(η2 − C3η5) +(γ + u2r0)(η2 − η3) + µ1η2 − A1 dη3 dt = µ1η3 dη4 dt = P1Xs(η1 − C1η4) + P2Xi(η2 − C2η4) + (1− u1)β2Yi(η4 − η5) +(µ2 − u3c0)η4 dη5 dt = P3Xi(η2 − C3η5) + (1− µ1)β2Ys(η4 − η5) + µ2η5 − A2 (4.29) and with the transversality conditions η1(tf ) = η2(tf ) = η3(tf ) = η4(tf ) = η5(tf ) = 0. (4.30) Moreover, the optimal control triple (u∗1, u ∗ 2, u ∗ 3) is characterized as follows u∗1 = max { 0,min { 1, 1 B1 [β1XsXi(η2 − η1) + β2YsYi(η5 − η4)] }} , u∗2 = max { 0,min { 1, 1 B2 [r0Xi(η2 − η3)] }} , u∗3 = max { 0,min { 1, −c0Ysη4 B3 }} . (4.31) Proof The system (4.29) governing the adjoint variables is obtained by finding the partial derivative of the Hamiltonian H (4.28) as follows: dη1 dt = − ∂H ∂Xs , dη2 dt = − ∂H ∂Xi , dη3 dt = − ∂H ∂Xr , dη4 dt = − ∂H ∂Ys , dη5 dt = − ∂H ∂Yi Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 99 with the transversality conditions: η1(tf ) = η2(tf ) = η3(tf ) = η4(tf ) = η5(tf ) = 0. Finally, differentiating the Hamiltonian H partially with respect to each of the three controls u1, u2 and u3, so that ∂H ∂u1 ≡ B1u1 + β1XsXi(η1 − η2) + β2YsYi(η4 − η5) = 0, ∂H ∂u2 ≡ B2u2 + r0Xi(η3 − η2) = 0, ∂H ∂u3 ≡ B3u3 + c0Ysη4 = 0, which when solved on the interior of the control set U gives the required characterization (4.31). Hence, the proof. 4.2. Simulations of the control model Here, the numerical simulations of the coupled system of state equations (4.25) with adjoint equations (4.29) are conducted. The standard iterative fourth order forward-backward Runge- Kutta method is used to solve the optimality system due to difference in time orientations of the state and adjoint equations [18, 23]. The optimal control problem is simulated over the time interval [0, 15] in years using these parameter values: µ1 = 0.5, r = 11.2, β1 = 0.4, P1 = 1.0, γ = 0.35, P2 = 0.2, P3 = 1.0, µ2 = 0.65, k = 30, β2 = 0.65, C1 = 0.25, C2 = 0.25, C3 = 0.15. In addition, the balancing weight constants are given as A1 = 1, A2 = 0.1, B1 = 0.01, B2 = 0.0001 and B3 = 0.008. The rate constants for the controls are chosen as r0 = 0.4 and c0 = 0.15. Thus, to minimize the objective functional (4.26) with disease prevention control only at low cost of implementation, the control u1 should be maintained at maximum (100%) effort for 8 years before being dropped to zero in final time as shown in the Figure 4.5(a). The control profile shown in Figure 4.5(b) indicates that single implementation of treatment control u2 requires maximum (100%) effort throughout the period of the control intervention. It is observed in Figure 4.5(c) that alternative food source control u3 should be at the upper bound for 14 years before coming to the lower bound. The combined implementation of the three controls for minimizing the objective functional is depicted in Figure 4.5(d), where it is can be observed that prevention control u1 should be maintained at the upper bound for 8 years before declining to the lower bound, treatment control u2 and provision of alternative food source u3 should be maintained at maximum (100%) for about 2 years and almost 1 year, respectively. The effects of combining all the three optimal controls on the dynamical behaviors of prey and predator populations are illustrated in Figure 4.6. The size of susceptible prey population with control reduces due to predation when compared with the size without control, but the population does not go extinct due to the presence of the combined optimal control as shown in Figure 4.6(a). This suggests the possibility of stable coexistence of predators with low density of prey population as indicated in [27]. The infected prey population with optimal control (Figure 4.6(b)) decreases sharply when compared with the case without control. Further, in Figure 4.6(c), it can be seen that the size of susceptible predator population with combined optimal control is higher than the size without control. Whereas, the presence of the combined optimal control reduces the population size of the infected predator in the ecosystem. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 100 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 O p tim a l c o n tr o l u 1 Time (a) 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 O p tim a l c o n tr o l u 2 Time (b) 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 O p tim a l c o n tr o l u 3 Time (c) 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time O p ti m a l c o n tr o l p ro fi le s u 1 u 2 u 3 (d) Fig. 4.5. Control profiles for single and combined implementations of u1, u2 and u3. 5. CONCLUSION A nonlinear mathematical model of prey-predator interacting populations in an ecological system with disease spread is studied. The five-dimensional eco-epidemiological model is analyzed with a view to providing insights into the behavior of ecosystem in the presence of disease. Seven possible equilibrium points, namely trivial, axial, disease-free, predator- free, infected prey-free, infected predator-free and interior equilibrium points are analytically determined. Asymptotic stability of the eco-epidemiological system is analyzed to investigate the behavior of the system around each of the obtained equilibrium points. Some conditions that guarantee the existence of each prey and predator species, as well as coexistence of both prey and predator populations are given using local linearization and Lyapunov functions techniques. In addition, the eco-epidemiological model is extended to include three time-dependent optimal controls, such as disease prevention in both species, treatment strategy for prey and provision of alternative food source for predator survival. Hence, the optimal control prey- predator model is analyzed using Pontryagin’s maximum principle in order to enhance the coexistence of both prey and predator species by reducing the risks of their extinction in the ecosystem, such that the populations of infected prey and predator species are minimized at minimum costs possible. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 101 0 5 10 15 0 5 10 15 20 25 30 Time S u s c e p ti b le P re y u 1 ≠ 0, u 2 ≠ 0, u 3 ≠ 0 u 1 =u 2 =u 3 =0 (a) 0 5 10 15 0 0.5 1 1.5 2 2.5 Time In fe c te d P re y u 1 ≠ 0, u 2 ≠ 0, u 3 ≠ 0 u 1 =u 2 =u 3 =0 (b) 0 5 10 15 0 2 4 6 8 10 12 14 Time S u s c e p ti b le P re d a to r u 1 ≠ 0, u 2 ≠ 0, u 3 ≠0 u 1 =u 2 =u 3 =0 (c) 0 5 10 15 0 2 4 6 8 10 12 Time In fe c te d P re d a to r u 1 ≠ 0, u 2 ≠ 0, u 3 ≠0 u 1 =u 2 =u 3 =0 (d) Fig. 4.6. Combined effects of optimal controls u1, u2 and u3 on the dynamics of the prey-predator model. The analysis of the eco-epidemiological model presented in this study can be applied to a wide spectrum of interactions between different kinds of living beings in the environment. ACKNOWLEDGEMENT The authors are thankful for the constructive comments and suggestions of the editor and the anonymous reviewers. Initial discussion on the original manuscript with Professor Mini Ghosh of Vellore Institute of Technology, Chennai, India is also acknowledged with thanks. REFERENCES 1. Abidemi, A. and Aziz, N.A.B. (2020). Optimal control strategies for dengue fever spread in Johor, Malaysia, Comput. Meth. Prog. Bio., 196, 105585. https://doi.org/10.1016/j.cmpb.2020.105585 2. Abimbade, S.F., Olaniyi, S., Ajala, O.A. and Ibrahim, M.O. (2020). Optimal control analysis of a tuberculosis model with exogenous re-infection and incomplete treatment, Optim. Control Appl. Meth., (41), 2349–2368. https://doi.org/10.1002/oca.2658 Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) 102 J.A. ADEMOSU, S. OLANIYI, S.O. ADEWALE 3. Adewale, S.O., Podder, C.N. and Gumel, A.B. (2009). Mathematical analysis of a TB transmission model with DOTS, Can. Appl. Math. Q., 17, 1–36. 4. Akanni, J.O., Akinpelu, F.O., Olaniyi, S., Oladipo, A.T., and Ogunsola, A.W. (2020). Modelling financial crime population dynamics: optimal control and cost-effectiveness analysis, Int. J. Dynam. Control, 8, 531–544. 5. AL-Darabsah, I., Tang, X. and Yaun, Y. (2016). A prey-predator model with migrations and delays, Discrete Contin Dyn Syst Ser B., 21(3), 737–761. 6. AL-Nassir, S. (2017). The dynamics and optimal control of a prey-predator system, Global J. Pure Appl. Math., 13, 5287–5298. 7. Bakare, E.A., Adekunle, Y.A. and Nwagwo, A. (2012). Mathematical analysis of the control of the spread of infectious disease in a prey-predator ecosystem, Int. J. Comp. Organ. Trends, 2(1), 27–32. 8. Bera, S.P., Maiti, A. and Samanta, G.P. (2015). A prey-predator model with infection in both prey and predator, Filomat, 29(8), 1753–1767. 9. Das, K.P. (2011). A mathematical study of a predator-prey dynamics with disease in predator, ISRN Appl. Math., Article ID 807486, 1–16. 10. Diva, A.R., Fatmawati, Windarto and Didik, K.A. (2018). Optimal control of predator- prey mathematical model with infection and harvesting on prey, Journal of Physics: Conf. Series, 974, 012050. 11. Doust, M.H.R. and Gholizade, S. (2014). An analysis of the modified Lotka-Volterra predator-prey model, Gen. Math. Notes, 25(2), 1–5. 12. Fleming, W.H. and Richel, R.W. (1975). Deterministic and Stochastic Optimal Control. New York: Springer. 13. Ghosh, M. (2010). Modeling prey-predator type fishery with reserve area, Int. J. Biomath., 3, 351–365. 14. Ghosh, M., Olaniyi, S. and Obabiyi, O.S. (2020). Mathematical analysis of reinfection and relapse in malaria dynamics, Appl. Math. Comput., 373, 125044. https://doi.org/10.1016/j.amc.2020.125044 15. Goswami, N.K. and Shanmukha, B. (2021). A mathematical analysis of Zika virus transmission with optimal control strategies, Comput. Meth. Diff. Equ., 9(1), 117–145. 16. Hugo, A., Makinde, O.D., Kumar, S. and Chibwana, F.F. (2017). Optimal control and cost effectiveness analysis for Newcastle disease eco-epidemiological model in Tanzania, J. Biol. Dyn., 11, 190–209. 17. Kumar, V. and Kumari, N. (2019). Controlling chaos in three species food chain model with fear effect, AIMS Math., 5(2), 828–842. 18. Lenhart, S. and Workman, J.T. (2007). Optimal Control Applied to Biological Models. London: Chapman & Hall. 19. LaSalle, J.P. (1976). The Stability of Dynamical Dystems, Regional Conference Series in Applied Mathematics. Philadelphia PA: SIAM. 20. Lotka, A.J. (1925). Elements of Physical Biology. Baltimore: Williams and Wilkins. 21. Nandi, S.K., Mondal, P.K., Jana, S., Haldar, P. and Kar, T.K. (2015). Prey-predator model with two-stage infection in prey: concerning pest control, J. Nonlinear Dynam., Article ID 948728 (2015), 1–13. 22. Obabiyi, O.S. and Olaniyi, S. (2019). Global stability analysis of malaria transmission dynamics with vigilant compartment, Electron. J. Differential Equations., 2019(09), 1– 10. 23. Okyere, E., Olaniyi, S. and Bonyah, E. (2020). Analysis of Zika virus dynamics with sexual transmission route using multiple optimal controls. Scientific African, 9, e00532. https://doi.org/10.1016/j.sciaf.2020.e00532 24. Olaniyi, S., Okosun, K.O., Adesanya, S.O. and Lebelo, R.S. (2020). Modelling malaria dynamics with partial immunity and protected travelers: optimal control and cost- effectiveness analysis. J. Biol. Dyn., 14, 90–115. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) STABILITY ANALYSIS AND OPTIMAL MEASURE 103 25. Olaniyi, S., Obabiyi, O.S., Okosun, K.O., Oladipo, A.T. and Adewale, S.O. (2020). Mathematical modelling and optimal cost-effective control of COVID-19 transmission dynamics, Eur. Phys. J. Plus, 135, 938. https://doi.org/10.1140/epjp/s13360-020-00954- z 26. Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V. and Mishchenko, E.F. (1962). The Mathematical Theory of Optimal Processes, New York: Wiley. 27. Prasad, B.S.R.V., Banerjee, M. and Srinivasu, P.D.N. (2013). Dynamics of additional food provided predator-prey system with mutually interfering predators. Math. Biosci., 246, 176–190. 28. Safi, M.A. and DarAssi, M.H. (2018). Mathematical analysis of a model for ectoparasite-borne diseases, Math. Meth. Appl. Sci., 41(17), 8248–8257. 29. Sani, A., Cahyono, E., Mukhsar and Rahman, G.A. (2014). Dynamics of disease spread in a predator-prey system, Adv. Stud. Biol. 6, 169–179. 30. van den Driessche, P. and Watmough, J. (2002). Reproduction numbers and sub- threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180, 29–48. 31. Volterra, V. (1926). Variazioni e Fluttaazioni of Numero di Individul in Specie Animali Conviventi, Mem. Accd. Linc., 2, 31–113. Copyright © 2021 ASSA. Adv Syst Sci Appl (2021) Introduction Model formulation Well-posedness of the model Boundedness of solutions Positivity of solutions Analysis of the model Existence of equilibrium points Trivial equilibrium point (E0) Axial equilibrium point (E1) Disease-free equilibrium point (E2) Predator-free equilibrium point (E3) Infected prey-free equilibrium (E4) Infected predator-free equilibrium (E5) Interior equilibrium point (E6) Stability analysis Stability of E0 Stability of E1 Stability of E2 Stability of E3 Stability of E4 Stability of E5 Stability of E6 Optimal control model Analysis of the control model Simulations of the control model Conclusion