Mathematical Modelling and optimal control of pneumonia disease in sheep and goats in Al-Baha region with cost-effective strategies

1.   Introduction

Pneumonia is a lung tissue infection with many causes. It is a major medical problem for sheep and goats of all ages. It is the most widespread problem associated with the lower respiratory tract of sheep that can be acute, chronic, or progressive and can be caused by bacteria, viruses, or parasites. It is an inflammatory response of the alveoli in lungs to infective agents, resulting in lung consolidation. It is a common disease of sheep in all major sheep-producing countries ^[1]. The disease is multi-factorial, resulting from dynamic interactions between host, infectious agent, and environmental factors. When the bacterial herds reach a certain threshold, host susceptibility, non-specific defense mechanisms, and lung defense mechanisms become compromised, allowing diseases to occur ^[2]. In younger animals, separate upper and lower respiratory tract bacteria, viruses, and parasites are frequently involved in developing pneumonia. These same disease-causing agents can cause pneumonia in adults. In sheep, a systemic virus known as Ovine Progressive Pneumonia virus (OPPV) can play an important role. In goats, pneumonia may be caused by a related systemic virus, the Caprine Arthritis, and Encephalitis Virus (CAEV). The term "systemic" means that OPPV and CAEV are viruses that, like the lungs, can affect multiple organs. These viruses can affect the brain, udder, and joints as well. Parasites (worms) can travel from the gastrointestinal tract into the lungs in certain climates, causing pneumonia. Sheep could convert and diversify several types of forages into valuable products for mankind, such as mutton, milk, and wool ^[3]. The importance of sheep, to the socio-economic well-being of animals in developing countries, cannot be emphasized hence it is necessary to study the diseases and syndromes that affect this species to enhance and sustain their productivity to meet the demand of the human herd's places upon them ^[4]. Sheep herds have incurred a steady increase worldwide because of the low initial investments and its dual purpose for (meat and wool) and milk; however, their production is often limited by various infectious diseases of which respiratory infections are important ^[5]. The severity of respiratory infections in sheep-like goats has been associated with several factors including environmental, nutritional, transportation, and management. Respiratory infection resulting in fatal outcomes usually occurs when a primary viral infection compromises host defenses and enhances the complications from a secondary bacterial infection ^[6]. Among infectious agents, bacterial pneumonia is responsible for outbreaks, sudden death, and high mortality in lambs ^[7]. The lamb losses relate to bacterial infections leading to pneumonia, diarrhea, and subsequent sepsis, which is a potential complication of pneumonia ^[8]. Over the years, pneumonia in small ruminants has caused enormous production losses with mortality between 71 and 90 among goats, but with less mortality in sheep ^[9]. Respiratory diseases in most sheep-producing countries result in lamb mortality, reduced growth rate with a substantial negative economic impact on animal husbandry because of the need to employ chemotherapeutic and vaccination programs ^[10]. Finally, it is the most common epidemic widespread among Sheep in Al-Baha region according to the Ministry of Environment, Water, and Agriculture data ^[11].

Nowadays, the mathematical modeling of the spread of diseases and optimal control problems for epidemics, models have attracted the attention of many researchers. It is playing a vital role in predicting the health and environmental impacts of epidemics and diseases for animals and humans. Little essential research has been done in the last decade on the dynamics of pneumonia. Some of them are proposed a model on pneumonia dynamics ^{[12,13,14,15,16,17,18,19,20,21]}. Most of these scientists have considered and estimated the basic reproductive number as a random variable by first developing and analyzing a deterministic model for transmission patterns of pneumonia. Tilahun et al. ^[21], proposed a co-infection model for pneumonia-typhoid and mathematically analyzed their characteristic relationship in case of cure and medical strategies. In 2018, Tilahun et al. ^[22] described a model of pneumonia-meningitis co-infection with the help of ordinary differential equations and some of the theorems. Tilahun et al. divided the population into compartments in accordance with the state of their health, such as susceptible $(S)$, vaccinated $(V)$, carrier $(C)$, infected $(I)$, and recovered $(R)$.

Fractional derivatives have been used efficiently to fine-tune many existing models of physical and natural phenomena ^{[23,24,25,26,27,28,29,30,31,32,33,34,35,36]}. In the field of Bio-Modeling, fractional derivatives were used to generalize the Hodgkin-Huxley model ^[37]. The model using fractional derivatives gave more accurate results than the model using classical derivatives. In the field of physical science, fractional derivatives were used to model viscoelastic substances with good agreement experimental results ^[38,39,40]. Also, Sherief et al. ^[41] have constructed a successful fractional model for thermoelasticity (see also ^[42,43]).

In this work we shall consider formulation of Tilahun et al. model using fractional derivatives in sense of Caputo. Which is used to improve a mathematical model for the transmission dynamics of pneumonia in the Al-Baha region of the Kingdom of Saudi Arabia. Al-Baha is a located in the southwestern part of Saudi Arabia. It is characterized by varying terrain that affects climate and vegetation. The geographic characteristics of the Al-Baha are described by distinct forests, high Juniperus procera woodland, and diverse vegetation ^[34]. The existence, uniqueness, non-negativity and boundedness of the solutions of the proposed model is proved. The basic reproduction number $R_{0}$ is determined, which helped us to determine the dynamical behavior of the system. Moreover, sufficient conditions of stability or instability of the model are obtained. The endemic equilibrium stable globally has been derived and analyzed under some conditions. Global dynamics have also been established to be completely determined. For the main parameters, the model bifurcation potential for the basic reproduction number is determined. Numerical simulations that we provided support the analytical results, which are obtained using generalized Euler method to handle the fractional derivatives. On the other hand, using Pontryagin's Maximum Principle ^[44], the optimal control problem is formulated with three control strategies: Disease prevention through education, treatment, and screening. We expect that the cost-effectiveness analysis of the adopted control strategies reveals that the combination of prevention and treatment is the most cost-effective intervention strategies to combat the pneumonia pandemic. We perform numerical simulation and present relevant results graphically. The results of this project will provide quantitative pneumonia knowledge that helps us to design disease monitoring tools.

The main advantage behind the use of fractional derivatives stems from its memory effect. This means that its definition as an integral over past times makes the effects of the stimuli retarded not instantaneous as in the classical models. This is in accord with what is found in nature. Other fractional operators such as the Caputo-Fabrizio and Atangana-Baleanu will be discussed in the future and compared with the Caputo one ^[45,46]. The important mathematical findings for the proposed model's dynamical behavior have also been numerically verified for an Al-Baha region.

2.   Model analysis

2.1.   Preliminaries

In this section, we give necessary definitions and some properties of fractional calculus. Denote by $\mathbb{R}_+$ the set of all non-negative real numbers. If $\nu\in \mathbb{R}_+$ is a non integer order, the fractional integral $J^{\nu}_mf(t)$ of the function $f: \mathbb{R}_+\longrightarrow \mathbb{R}$ is defined as

where $\Gamma(z) = \int_0^\infty e^{-t}t^{z-1}dt$ is the Euler Gamma function.

Definition 1 (^[47]). The Riemann-Liouville fractional derivative $D^{\nu}f(t)$ of order $\nu > 0$, $n-1 < \nu < n$, $n\in \mathbb N$ is defined as

where $n$ is a natural number satisfying $n = [\nu] + 1$ and $[\nu]$ denotes the integer part of $\nu$.

Definition 2 (^[47]). The Caputo fractional derivative $D^{\nu}f(t)$ of order $\nu > 0$, $n-1 < \nu < n$, $n\in \mathbb N$ is defined as

The Riemann-Liouville and Caputo derivatives for $0 < \nu < 1$ are related by the following formulas:

Consider the initial value problems for fractional differential equations in the form of

where the fractional derivative $D^{\nu} = {}_m^{c}D^{\nu}$ is in sense of Caputo's definition, and the function $f(X(t)) : \mathbb{R}\times \mathbb{R}^2 \longrightarrow\mathbb{R}^2$ is called a vector field. Particularly, $\mathbb{R}^2$ endowed a proper norm $\parallel.\parallel$ becomes a complete metric space. Denote by

Lemma 1 (^[48]). Assume that the function $f:\mathcal{D}\longrightarrow\mathbb{R}^2$ satisfies the following conditions:

(1) $f(X(t))$ is Lebesgue measurable with respect to $t$ on $\mathcal{J}$.

(2) $f(X(t))$ is continuous with respect to $X$ on $\hbar$.

(3) There exists a real-valued function $g(t)\in L^2(\mathcal{J})$ such that

for almost every $t\in\mathcal{J}$, $X\in\hbar$.

(4) There exists a real-valued function $\mu(t) \in L^4(\mathcal{J})$ such that

for almost every $t\in\mathcal{J}$ and all $X_1, X_2\in\hbar$.

Then there exists a unique solution of the initial value problem (2.2) on $[t_0-\varepsilon, t_0+\varepsilon]$ for some positive number $\varepsilon$.

Lemma 2 (^[49]). Assume that the function $f: \mathbb{R}_+\times \mathbb{R}^2\longrightarrow \mathbb{R}^2$ satisfies the following conditions:

(1) The function $f(X(t))$ is Lebesgue measurable with respect to $t$ on $\mathbb{R}$.

(2) The function $f(X(t))$ is continuous with respect to $X$ on $\mathbb{R}^2$.

(3) The function $\frac{\partial f(X(t))}{\partial X}$ is continuous with respect to $X$ on $\mathbb{R}^2$.

(4) There exists two positive constants $\omega$, $\lambda$ such that

for almost every $t\in\mathbb{R}$, and all $X\in\mathbb{R}^2$.

Then, there exists a unique function $X(t)$ on $(-\infty, +\infty)$ solving the initial value problem (2.2).

Lemma 3 (^[50], Lemma 3). Let $u(t)$ be a function such that $D^{\nu}u(t)$ exists for all $t$ and satisfying

where $0 < \nu < 1$, $(\lambda, \mu)\in\mathbb{R}^2$ and $\lambda\neq0$ and $t_0\geq0$ is the initial time. Then

where $E_{\nu, 1}(t) = \sum_{k = 0}^\infty\frac{t^k}{\Gamma(kq+1)} > 0$ is the Mittag Leffler function.

To prove the global stability of $E^*$, we need the following two lemmas.

Lemma 4 (^[51,52]). Let $u(t)\in \mathbb{R}^+$ be a continuous and derivable function. Thus, for any time instant $t \geq t_0$,

Lemma 5 (^[53]). (Fractional LaSalle's invariance principle). Let $\Omega$ be a bounded closed set. Every solution of $D^\nu X(t) = f(X(t))$ that starts from a point in $\Omega$ and remains in $\Omega$ for all time. If $\exists$ $L(X(t)): \Omega\longrightarrow \mathbb{R}$ with continuous first partial derivatives satisfies the following condition:

Let $E = \{X\in \Omega: D^\nu L(X(t)) = 0\}$ and $\Delta$ be the largest invariant set of $E$. Then every solution $X(t)$ originating in $\Omega$ tends to $\Delta$ as $t \longrightarrow \infty$. Particularly, when $\Delta = \{X^*\}$, then $X\longrightarrow X^*$, as $t\longrightarrow \infty$.

Proof. Let $\tilde{L}(X(t)) = J^{1-\nu}_0 (L(X(t))-L(X(0)))$. In (2.1), the relationship between Riemann-Liouville and Caputo fractional derivatives yields

Also, $E = \{X\in\Omega: \dot{\tilde{L}}(X(t)) = 0\}$ and $\Delta$ is the largest invariant set in $E$. Thus, using the classical LaSalle's invariance principle, Theorem 1 in ^[54], every solution starting in $\Omega$ tends to $\Delta$ as $t\longrightarrow \infty$.

2.2.   Bifurcation analysis

In the present work, we will investigate the nature of the bifurcation of the solution trajectories of the proposed model by using the method introduced by Carlos and Song ^[55], and the method is based on the use of center manifold theory.

Theorem 1 (^[31] Centre manifold theory). : Let $f: R^n\times R^n\longrightarrow R$ and $f\in C^2(R^n\times R^n)$. Then consider the following general system of ODEs with a parameter $\varphi$.

where $y = 0$ is an equilibrium point of the system (2.2) (that is, $f(0, \varphi) = 0, \forall \varphi$. Assume that:

A1: $A = \Big(\frac{\partial f_i}{\partial y_j}(0, 0)\Big)$ is the linearization matrix of (2.2) around the equilibrium point 0 with $f$ evaluated at 0. Zero is a simple eigenvalue of $A$ and other eigenvalues of $A$ have negative real parts;

A2: Matrix $A$ has a right eigenvector $w$ and a left vector $v$ (each corresponding to the zero eigenvalue).

Let $f_k$ be the k-th component of $f$ and

The local dynamics of the system (2.2) around 0 is totally determined by the signs of $\ell$ and $\Bbbk$:

ⅰ. $\ell > 0$, $\Bbbk > 0$. When $\varphi < 0$ with $|\varphi|\leqslant1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0 < \varphi < < 1$, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium,

ⅱ. $\ell < 0$, $\Bbbk < 0$. When $\varphi < 0$ with $|\varphi|\leqslant1$, 0 is unstable; when $0 < \varphi < < 1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium,

ⅲ. $\ell > 0$, $\Bbbk < 0$. When $\varphi < 0$ with $|\varphi|\leqslant1$, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0 < \varphi < < 1$, 0 is stable, and a positive unstable equilibrium appears,

ⅳ. $\ell < 0$, $\Bbbk > 0$. When $\varphi$ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

2.3.   Model description and formulation

In this work, we introduce a fractional-order for the model developed by Tilahun GT et al. in ^[21], which describes the dynamics of transmission of pneumonia disease. Here, we study the dynamics of pneumonia disease in sheep and goats in Al-Baha region between five compartments based on the disease status, that is: Susceptible $(S)$, vaccinated $(V)$, carrier $(C)$, infected $(I)$, and recovered $(R)$, namely SVCIR model. More precisely, a fractional order SVCIR model is considered in sense of Caputo derivative $D^\nu$, $0 < \nu \leq 1$, and is given by

with initial data $S(0) = S_0$, $V(0) = V_0$, $C(0) = C_0$, $I(0) = I_0$, and $R(0) = R_0$.

The model assumes that before the disease outbreak, a fraction of the herds was vaccinated at a rate of $(p)$, with a fraction of herds susceptible at $(1-p)$. The Susceptible class is increased from the vaccinated class, which contains animals who have been vaccinated but did not respond to vaccination due to a waning rate of $\phi$, and from recovered class in which those individuals who loss their temporary immunity by $\delta$ rate. Animals from the susceptible class, on the other hand, move to the vaccinated class at a rate of $\vartheta$. The infection force is defined as $\lambda = \xi (\frac{I(t)+\varUpsilon C(t)}{N})$, where $N = N(t)$ represents the total population of sheep or goats at the moment $t\geq 0$. Thus, $N = S(t)+E(t)+C(t)+I(t)+R(t)$ is a non-constant time-dependent quantity (except if $\alpha = 0$, which is not the case since $\alpha = 0.057$ as appears in Table 1). The susceptible class is infected with a force of infection $\lambda$ by either carriers or symptomatically infected animals, where $\xi = k\tau$, $k$ is contact rate, $\tau$ is the probability that a contact is effective to cause infection and $\varUpsilon$ is transmission coefficient for the carrier. If $\varUpsilon > 1$, the carriers are more likely to infect susceptible than infective. If $\varUpsilon = 1$, both carriers and infective have an equal chance of spreading the susceptible, but if $\varUpsilon < 1$, infective have a greater chance than carriers of contacting the susceptible. Because the model assumes that vaccination isn't 100 percent effective, vaccinated classes (V) have a small chance of being infectious or carriers, and the force of infection for the vaccinated class is $\lambda_\nu = \varepsilon\lambda$, where $0\leqslant \varepsilon < 1$ and $\varepsilon$ is the proportion of the serotype not covered by the vaccine. By the force of infection, newly infected animals become either carriers with a probability of $\rho$ to join the carrier class $C$ or infected with a probability of $1-\rho$ to join the infected class $I$. The carrier class can develop disease symptoms, or they can screen and join the infected class at a rate of $\chi$, or they can recover by gaining natural immunity at a rate of $\beta$. Animals in the infected class move to the recovered compartment at a per capita rate of $\eta$ by treatment, with treatment efficacy of $q$ proportion of animals join the recovered class or join the carrier class at a rate of $(1-q)$ by adapting the treatment, or die from the disease at a rate of $\alpha$. All of the parameters are positive, and $\mu$ is the natural death rate of the animals in all compartments.

2.4.   Existence and uniqueness

This section deals with the existence, uniqueness, non-negativity and boundedness of the solutions of the system (2.4). Denote by $\Bbb{R}_+$ the set of all non-negative real numbers and assuming that $\Omega = \{(S, V, C, I, R)\in\Bbb{R}^5_+: S\geq0, \, V\geq0, \, C\geq0, \, I\geq0, \, R\geq0, \max (|S|, |V|, |C|, |I|, |R|)\leq N\}$. Following ^[50], the existence and uniqueness of the solutions of the fractional-order model (2.4) is proved in the region $\Omega \times (0, T]$ (see Figure 1).

Theorem 2. If $X_0 = (S(0), V(0), C(0), I(0), R(0))\in \Omega$ is an initial condition, then the solution $X = (S(t), V(t), C(t), I(t), R(t))\in \Omega$ to the fractional-order model (2.4) is unique, for $t \geq 0$.

Proof. Let

Hence, for $X, \overline{X}\in \Omega$, one obtains

where

Therefore $H(X)$ satisfies the Lipschitz condition, which imply that the solution of fractional-order model (2.4) exists and unique.

2.5.   Non-negativity and boundedness of the solution

Theorem 3. If $p < 1$, $q < 1$, $\rho < 1$, the fractional-order model (2.4) has non-negative solutions.

Proof. One has

Thus, by using Lemmas 5 and 6 in ^[56], the solutions of fractional-order model (2.4) are non-negative.

Theorem 4. The fractional-order model (2.4) has uniformly bounded solutions start in

Proof. Denote by

the total herds at time $t$. By direct sum to the equations in fractional-order model (2.4), one obtains

Thus

Following to ^[57], one obtains

where $M_\nu$ is the Mittag-Leffler function. According to ^[56], one obtains:

Thus, the solutions of fractional-order model (2.4), starting in $\Omega^+$, are uniformly bounded in the region $\Omega^+$.

Remark 1. Obviously, the set $\Omega^+$ is positively invariant with respect to the fractional-order model (2.4).

3.   Stability analysis

The basic reproduction number $R_{0}$ has been computed in this section. Moreover, the local and global asymptotic stability of the infectious-free equilibrium point is studied. In addition, by building appropriate Lyapunov functions, the local and global asymptotic stability of the endemic equilibrium point is being studied.

3.1.   The infectious-free equilibrium point

For $I = C = 0$, one obtains

Therefore

where $\Phi = \frac{\mu+\phi-p\mu}{\mu+\phi+\vartheta}$ and $\Theta = \frac{\vartheta+p\mu}{\mu+\phi+\vartheta}$. We can easily obtain the infection-free equilibrium $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$.

3.2.   The endemic-equilibrium point

For $I^\ast > 0$, the endemic equilibrium is denoted by $E^\ast = (S^\ast, V^\ast, C^\ast, I^\ast, R^\ast)$ and defined as a steady state solutions for the model (2.4). This can occur when there is a persistence of the disease. It can be obtained by equating the system of Eq (2.4) to zero, that is,

We used the third and fourth equations of model (2.4) to get the following result:

where

By substituting in the fifth equation of model (2.4), one obtains

where

The second equation gives the following:

where

The first equation of model (2.4) gives the following:

where

Thus

Substituting in the second equation, one obtains

Thus, one obtains

where $I^\ast$ is given by (3.1).

3.3.   The basic reproductive number

The basic reproductive number in the classical sense is the expected number of secondary infections produced by one infectious individual in a totally susceptible population during the individual's period of infectiousness. This is a static quantity (independent of time) obtained by the standard "next-generation matrix approach" (see Driessche et al. ^[58]).

Obviously, the state variables $(S(t), V(t), C(t), I(t), R(t))$ remain in the biologically meaningful region $\Omega$ is a positively invariant for model (2.4). The model equations are re-written starting with newly infective classes:

Then by the principle of next generation matrix we can obtained, $x = (S, V, R)^{\mathcal T}$ and $y = (C, I)^{\mathcal T}$, then we have

where

The next is obtaining the Jacobian matrix of $\aleph(y)$ and $\hbar(y)$ with respect to $C$ and $I$ at the disease free equilibrium $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$. To simplify our work, we assigned $\aleph(y)$ and $\hbar(y)$ in the following way,

where

The Jacobian matrices of the two matrices $\aleph(y)$ and $\hbar(y)$, at the infection-free equilibrium point $E_0$, are, respectively. The Jacobian matrix of $\aleph(y)$ and $\hbar(y)$ is obtained, at the infection-free equilibrium point $E_0$, by $F$ and $\Gamma$ respectively,

Since $\lambda = a C+ b I$, where $a = \frac{k\tau \varUpsilon}{N}$ is the transmission coefficient for the carrier compartment and $b = \frac{k\tau}{N}$ is also the transmission co-efficient for the infective compartment. Then

Similarly the entry members of $\Gamma$ is also obtained

The entry members of $F$ is obtained in the following way:

Thus

The matrix $F.\Gamma^{-1}$ is given by

where

The eigenvalue of $F.\Gamma^{-1}$ can be obtained,

Then the eigenvalues are

Thus the basic reproductive number is given by

By substitution of $d_1$, $d_2$, $d_3$, and $d_4$, the basic reproduction number becomes

3.4.   Local stability

Lemma 6. The infection-free equilibrium $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$ is locally asymptotically stable in $\Omega$ if $R_{0} < 1$ and unstable if $R_{0} > 1$.

Proof. For the fractional-order model (2.4) at $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$, the Jacobian matrix $J(E_0)$ is given by

Thus, the eigenvalue is given by

Then

Therefore

or

From Eq (3.2), one obtains

or

Then, by Routh-Hurwitz criteria, Eq (3.4) have two strictly negative root $\sigma_2, \sigma_3$.

On the other hand, the determinant of Eq (3.3) can be obtained from the equation

where

Thus, according to the Routh-Hurwitz criterion, for $R_{0} < 1$, the free equilibrium $E_0$ is locally asymptotically stable if $h_1 > 0$ and $h_2 > 0$. So, we need to prove that $h_1 > 0$ and $h_2 > 0$ when it holds that $R_{0} < 1$.

If $R_{0} < 1$, one obtains

which implies

That is

Thus, one obtains

If $R_{0} < 1$, one obtains

That is

If $\Gamma = min\{\mu+\alpha+\eta, \mu+\beta+\chi\}$, one obtains

Thus, one obtains

Then

Thus, the disease free equilibrium is locally asymptotically stable if $R_{0} < 1$. Otherwise, when $R_{0} > 1$ the system (2.4) is unstable.

Lemma 7. For $R_{0} > 1$, a unique endemic equilibrium point $E^\ast$ exist and no endemic equilibrium otherwise.

Proof. For the disease to be endemic, we must have $D^{\nu}C > 0$ and $D^{\nu}I > 0$, that is:

From the first inequality of (3.5),

Then

Using the fact $\frac{(S(t)+\varepsilon V(t))}{N} \leqslant1$, one obtains

From the second inequality of (3.5),

Then

Using the fact $\frac{(S(t)+\varepsilon V(t))}{N} \leqslant1$, one obtains

Substituting (3.6) into (3.7), one obtains

Then, by rearranging and cancelling of $I(t)$ in both sides, one obtains

Thus a unique endemic equilibrium exist when $R_0 > 1$.

Now, we investigate the local stability of the endemic-equilibrium point $E^\ast$. We use the center manifold theory to purpose stability from endemic equilibrium if $R_{0} > 1$. The solution of equilibrium point on the model assumed that occur around a point bifurcation, when $R_{0} = 1$. $\beta = \beta^*$ as bifurcation parameter of the equation of $R_0$ with $\beta = \beta^*$, obtained. This is done by renaming the variables as follows: Assume that

further by introducing the vector notation

Thus, the model can be written in the form of

where

as follows:

where

Then, the Jacobian system at the disease free is

Suppose that $\xi = \xi^*$ is a bifurcation parameter; the system (3.8) is linearized at the disease free equilibrium point when $\xi = \xi^*$ with $R_{0} = 1$, solving for $\xi^*$ for $R_{0} = 1$ from:

Thus

The system (3.8) with $\xi = \xi^*$ has a simple zero eigenvalues, hence the centre manifold theory will be used to analyse the dynamics of the system near $\xi = \xi^*$. The Jacobean matrix near $\xi = \xi^*$ has a right eigenvector associated with the zero eigenvalue $w = (w_1, w_2, w_3, w_4, w_5)^T$ is given by

The system of equation become

By solving (3.9), one obtains

The left eigenvectors of $J_E$ associated with the zero eigenvalue at $\xi = \xi^*$ is given by

System (3.9) gives

Thus

By solving (3.10), one obtains

To compute $\ell$ and $\Bbbk$, we use a formula explained in (2.3) as follows

where

Taking into account system (3.11) and considering only the non-zero components of the left eigenvectors $v_3$ and $v_4$, then one obtains

where

Since the coefficient $\Bbbk$ is always positive, it is the sign of the coefficient a and consequently the sign of the quantity $\ell_0$ which determines the local dynamics of the disease around disease-free equilibrium. Therefore, $\ell > 0$ depending on whether $\ell_0$ is greater or less than 0. Thus we have established the following result.

Lemma 8. If $\ell_0 > 0$, $\ell > 0$ then model system (2.4) has a backward bifurcation at $R_0 = 1$, otherwise, $\ell < 0$ and a unique endemic equilibrium point $E^\ast = (S^\ast, V^\ast, C^\ast, I^\ast, R^\ast)$ locally asymptotically stable for $R_{0} > 1$, but close to 1.

According to Manifold Center theorem ^[31], if $\ell < 0$ and $\Bbbk > 0$ so locally asymptotically stable. The system of equations will experience bifurcation transcritical at $R_{0} = 1$, if $\ell < 0$ and $\Bbbk > 0$. So, because $R_{0} > 1$ of the endemic equilibrium $E^\ast = (S^\ast, V^\ast, C^\ast, I^\ast, R^\ast)$ locally asymptotically stable.

3.5.   Global stability

Theorem 5. The infection-free equilibrium point $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$ of the fractional-order model (2.4) is globally asymptotically stable.

Proof. Consider the following positive definite Lyapunov function:

From Lemma 4, one obtains

At $E_0 = \left(\frac{\pi\Phi}{\mu}, \frac{\pi\Theta}{\mu}, 0, 0, 0\right)$, one obtains

Thus, $D^\nu L_1(S, V, C, I, R) < 0$ for all $(S, V, C, I, R)\in \Omega$. Furthermore $D^\nu L_1(S, V, C, I, R) = 0$ implies that $S = S_0$, $V = V_0$, $C = C_0$, $I = I_0$, and $R = R_0$. Therefore, the singleton $\{E_0\}$ is the only invariant set such that $D^\nu L_1(S, V, C, I, R) = 0$. Therefore, by fractional La-Salle Invariance Principle, Lemma 5 gives conclusion that $E_0$ is globally asymptotically stable on $\Omega$.

Theorem 6. The positive endemic equilibrium point $E^\ast = (S^\ast, V^\ast, C^\ast, I^\ast, R^\ast)$ of the fractional-order model (2.4) is globally asymptotically stable.

Proof. Consider the following positive definite Lyapunov function:

From Lemma 4, one obtains

Thus, $D^\nu L_2(S, V, C, I, R) < 0$ for all $(S, V, C, I, R)\in \Omega$. Furthermore $D^\nu L_2(S, V, C, I, R) = 0$ implies that $S = S^*$, $V = V^*$, $C = C^*$, $I = I^*$, and $R = R^*$. Therefore, the singleton $\{E^*\}$ is the only invariant set such that $D^\nu L_2(S, V, C, I, R) = 0$. Therefore, by fractional La-Salle Invariance Principle, Lemma 5 gives conclusion that $E^*$ is globally asymptotically stable on $\Omega$.

4.   Extension the model to an optimal control for $\nu=1$

In this section, as in ^[21], we use optimal control strategies on the model (2.4) for $\nu = 1$. This helped us in determining the most effective intervention options for eradicating the disease in the period allowed. The best control model is an extension pneumonia model that includes the three controls listed below:

ⅰ. $u_1$ a prevention effort, that protect susceptible from contacting the disease.

ⅱ. $u_2$ a treatment effort, to minimize infection by treating infectious.

ⅲ. $u_3$ a screening effort, to help carriers to screen themselves.

We get the following optimal control model of pneumonia after including $u_1$, $u_2$ and $u_3$ in pneumonia model (2.4):

To study the optimal levels of the controls, the control set as

is Lebesgue measurable.

Our aim is to obtain a control $U$ and $S$, $V$, $C$, $I$ and $R$ that minimize the proposed objective function

where $b_1$, $b_2$ and $w_i$ are positive (see ^[59]). The expression $\frac{1}{2}w_iu_i^2$ represents cost which is associated with the controls $u_i$. We seek to find an optimal triple controls $(u^*_1, u^*_2, u^*_3)$ such that:

4.1.   Existence of an optimal control

The boundedness of solutions of system (4.1) for a finite time interval is used to prove the existence of optimal control (see ^[60]). To prove that, we check that the following properties are satisfied.

(1) The set of controls and corresponding state variables is non-empty.

(2) The control set $U$ is convex and closed.

(3) The right hand side of the state system is bounded by a linear function in the state and control.

(4) The integrand of the objective functional is concave on $U$.

(5) The function is bounded below by $a_2-a_1(|u_1^2|+|u_2^2|+|u_3^2|)^{\frac{\alpha}{2}}$ where $a_1, a_2 > 0$ and $\alpha > 1$.

An existence result in (^[61], Theorem 9.2.1, p. 182) for the state system (4.1) with bounded coefficients is used to give condition (1). The control set $U$ is convex and closed by definition. The right hand side of the state system (4.1) satisfies condition (3) as the state solutions are a priori bounded. The integrand in the objective functional, $b_1S+b_2C+b_3I+\frac{1}{2}\sum_{i = 1}^{3}w_iu_i^2$ is clearly concave on $U$. Moreover, there are $a_1, a_2 > 0$ and $\alpha > 1$ satisfying

because, the state variables are bounded. Finally under assumption (5), there exists an optimal control $(u_1, u_2, u_3)$ that minimizes the objective functional $J(u_1, u_2, u_3)$.

4.2.   The Hamiltonian and optimal system

By using the principle of Pontryagin et al. ^[44], we got the necessary conditions which is satisfied by optimal pair. By this principle, a Hamiltonian $H$ is given by

where

and $\lambda_i$, $i = 1, 2, 3, 4, 5$ are the adjoint variable functions that be determined suitably by applying Pontryagin's maximal principle ^[44], and Fleming and Rishel ^[60,62] for existence of the optimal control pairs.

Theorem 7 (^[21]). For an optimal control set $u_1$, $u_2$, $u_3$ that minimizes $J$ over $U$, there is an adjoint variables, $\lambda_i$, $i = 1, 2, 3, 4, 5$, such that

With transversality conditions, $\lambda_i(t_f) = 0$, $i = 1, ..., 5$. Furthermore, one obtains the control set $(u^*_1, u^*_2, u^*_3)$ characterized by

where

4.3.   Uniqueness of the optimality system

As in ^[21], we can obtain the uniqueness of solutions of the optimal system (4.1) for small time interval.

Lemma 9 (^[21]). The function $u^*(s) = \min(\max(s, a), b)$ is Lipschitz continuous in $s$, where $a < b$ are fixed positive constants.

Theorem 8 (^[21]). For $t \in [0, t_f]$, the bounded solutions to the optimal system (4.1) are unique.

Using different combinations of the controls, like one control only at a time, two controls at a time and also all controls at a time, that we analyze and compare numerical results from simulations. We used $b_1 = 300$, $b_2 = 150$, $w_1 = 2$, $w_2 = 2$ and $w_3 = 6$ for simulation of Pneumonia model with optimal control and also for cost-effectiveness analysis. Additionally we used $S(0) = 900$, $V(0) = 200$, $C(0) = 100$, $I(0) = 100$, $R(0) = 10$ as initial values.

4.4.   Control with prevention only

Using Prevention effort $(u_1)$ of susceptible without treatment ($u_2$ = 0) and with no screening ($u_3 = 0$). Only Preventive intervention is used to replicate the model.

Figure 2 depicts the reduction in infectious and carrier herds as a result of preventative efforts. This can be attributed to the fact that prevention reduces the level of animals joining infective and carrier compartments. This means that better prevention reduces the burden of pneumonia infection.

4.5.   Control with treatment only

Using treatment effort $(u_2)$ without prevention $(u_1 = 0)$ and with no screening $(u_3 = 0)$.

Figure 3 shows a decrease in infectious herds for the first four months, after which it begins to increase. Those who had the disease earlier are being treated, which is why the number of infective animals is decreasing over the first four months. Then, due to a lack of protection, newly infected animals began to join both the infective and carrier classes. As a result, the number of infective begins to rise after four months of decrease, and the number of carriers begins to rise after five months.

4.6.   Control with screening only

Using screening ($u_3$) but without prevention ($u_1 = 0$) and no treatment of infectious $u_2 = 0$. Carriers can move into the infective classes and begin treatment with the use of screening.

Figure 4(b) shows a decrease in the carrier herds for the first five months, then an increase due to a lack of protection. Susceptible become infected and join both the carrier and infective classes. As a result, screening alone may not be sufficient to reduce the burden of pneumonia infection.

4.7.   Control with prevention and treatment

Using prevention ($u_1$) and treatment ($u_2$) and without screening ($u_3 = 0$). As an intervention technique, we utilized prevention and therapy.

Figure 5 shows that the number of infective and also carriers decreases over time. As a result, this strategy is effective in eliminating the disease from the herds within a set time frame.

4.8.   Control with prevention and screening

Using prevention $u_1$ and screening ($u_3$) and without treatment ($u_2 = 0$). We used prevention and screening in this technique.

In Figure 6, the curve for optimal control is above the curve for no control in the first. Because there is no therapy, but animals from carrier groups are joining the infective compartment through screening, and there are a number of infected persons in the compartment prior to prevention who are not receiving treatment, the curve is expected to rise for the time being. After a period of time, the number of infectious animals decreases because new infections are no longer occurring as a result of prevention methods, and because there is no treatment, the number of infective animals decreases owing to disease-causing mortality and natural death rates.

4.9.   Control with treatment and screening

Using treatment ($u_2$) and screening $u_3$ and without prevention ($u_1$). As an intervention, we used treatment and screening controls.

Figure 7 shows how optimal control of the combination of treatment and screening helps to reduce the infectious as well as the carrier herds, resulting in the disease being eliminated in the community.

4.10.   Control with prevention, treatment and screening

Using all the three controls, prevention $u_1$ treatment of infective ($u_2$) and screening of carriers $u_3$. We implement all three controls interventions, which helps in the reduction of the objective function.

Figure 8 shows that as a result of the intervention measures, the number of infectious and carrier herdss decreases at the prescribed period. As a result, implementing this technique helps in the elimination of pneumonia within a set time frame.

5.   Numerical results

5.1.   Application of generalized Euler method

In sense of Odibat and Momani ^[63], which is a generalization of Euler's classical method, we consider the following problem

The interval $[0, b]$ over which we want to find the problem's solution (5.1) can be subdivided into $k$ subintervals $[t_j, t_{j+1}]$ of equivalent width $h = b/k$, by using the nodes $t_j = jh$ for $j = 0, 1, ..., k$. Expand $x$ about $t = t_0 = 0$ and using the generalized Taylor's formula ^[64]. Suppose that $x$, $\mathcal{D}^\nu x$, and $\mathcal{D}^{2\nu}x$ are continuous on $[0, b]$. Thus, there is a value $a_1$ for each value $t$, implying that

When $(\mathcal{D}^\nu\, y(t))(t_0) = f(t_0, x(t_0))$ and $h = t_1$ are substituted into Eq (5.2), the result is an expression for $y(t_1)$:

If the step size $h$ is small enough, the second-order term (involving $h^{2\nu}$) can be neglected and the result obtained.

Repeat the process until a series of points approximates the solution $x(t)$. When $t_{j+1} = t_j + h$, the general formula for GEM is

for $j = 0, 1, ..., k-1$. It is clear that if $\nu = 1$, then the GEM (5.3) reduces to the classical Euler's method. Thus

In this way, one obtains

where $0 < m = \frac{h^\nu}{\Gamma(\nu+1)} < 1$ is the step size and subject to initial conditions, $S(0) = S_0$, $V(0) = V_0,$ $C(0) = C_0$, $I(0) = I_0,$ $R(0) = R_0$. We have $S_{n+1} = S_{n} = S$, $V_{n+1} = V_{n} = V$, $C_{n+1} = C_{n} = C$, $I_{n+1} = I_{n} = I$, and $R_{n+1} = R_{n} = R$ at a fixed point.

5.2.   Numerical simulations

Using a set of parameter values that existed in Table 1, the numerical SVCIR model proposed for different fractional-order values $\nu$ is simulated by using the generalized Euler method to support analytical results and to evaluate the control strategy effectiveness. We constructed an algorithm to simulate the results by using the Mathematica program. The numerical results have been presented for the initial 10 months $\approx 300$ days.

Figures 9 and 10 show the global stability of the model's infection-free equilibrium (2.4) using potentially different conditions when $\mathcal{R}_0 = 0.0083 < 1$. The solutions of (2.4) converge to the unique disease-free equilibrium $E_0 = (0.1192, 0.2961, 0, 0, 0)$, as predicted. $E_0$ is asymptotically stable globally. We have shown the dynamical behavior of each state variable from the proposed pneumonia disease model in Figures 9 and 10 for varying values of the fractional-order parameter $\nu$. Figure 10(a–e) shows the phase plot ($S-V-C-I-R$) for different values of $\nu = 0.7, 0.75, 0.85, 0.95, 1$, respectively. Although the number of infected animals has decreased significantly for smaller fractional orders, the number of susceptible animals has increased, as shown in Figure 9(a), (d).

Figure 9(a) shows decreasing behavior of the susceptible of certain groups of animals for decreasing values of $\nu$ with integer-order derivatives operator predicting more susceptibility.

Figure 9(b) depicts the significance of the fractional-order parameter $\nu$ for the vaccinated individuals of the model. This plot shows a substantial decrease in the number of vaccinated individuals for decreasing values of $\nu$ under integer-order predicting 50 more cases in comparison with Caputo operator assuming $\nu = 0.70$.

As shown in Figure 9(c), the Caputo operator for $\nu = 0.7$ predicts 350 carrier cases, whereas integer-order suggests 150 cases.

Figure 9(d), in the first 70 days, shows decreasing behavior of the infectious of certain groups of animals for increasing values of $\nu$. While, in the interval from 70–300 days, it shows decreasing behavior of the infectious of certain groups of animals for decreasing values of $\nu$. Thus, the Caputo operator gives better agreement with the real situation in this regard.

Further 9(e), the integer-order derivatives operator suggests that 20 recovered cases. While we observed that the Caputo operator with $\nu = 0.7$ generates 60 additional recovered cases. Thus, the Caputo operator gives better agreement with the real situation in this regard. Thereby recommending the use of the integer-order derivative operator over the Caputo.

6.   Discussion of results

A fractional-order mathematical model was intended for describing the pneumonia disease in sheep and goats in the Al-Baha region of the Kingdom of Saudi Arabia. We have divided the total herds into five compartments, namely susceptible, vaccinated, carrier, infective, recovered herds, and investigated the dynamical behavior of this model. As it was explained in the introduction, over the years, pneumonia in small herds has caused enormous production losses with mortality between 71 and 90 among goats, but with less mortality in sheep. In the analysis of the dynamic of the virus, serious attention is needed. As a result, we provided insights into the complex behavior of the Pneumonia infection model in a fractional sense using the Caputo operator. In light of this, we employ the initial conditions as $S(0) = 900, V(0) = 200, C(0) = 100, I(0) = 100$ and $R(0) = 10$, and the values of the parameters of the model as in Table 1.

One of the great advantages of the Caputo fractional derivative is that it minimizes the number of susceptible and infected while maximizing the number of recovered population from Pneumonia infection. Also, we note that when the $S(t)$ decreases, $C(t)$, as well as $I(t)$, get increases. However, at some point, both the $I(t)$ start to decrease. For the fractional variant the curves continue to decrease. And this gets $R(t)$ increased as seen in Figures 9(e). This is one of the advantages of fractional order over integer order.

As is well known, epidemiological models which incorporate the control strategies can be useful to both controls the spread of disease and minimize the intervention costs. For our model, it is natural to consider the vaccination rate coefficient as a control to reduce the disease burden. Thus, it is important and interesting to prove the existence of optimal control, characterize the optimal control, prove the uniqueness of optimal control, compute the optimal control numerically and investigate how the optimal control depends on various parameters in the models. Furthermore, one of the most important and effective parameters called the force of infection $\lambda$ which is the rate at which susceptible individuals acquire an infectious disease has been varied for increasing and decreasing values in order to see further the dynamical characteristics of the model. The results depicted in Figures 2–8, with the decreasing increasing values of $\lambda$, can clearly be observed how effective is the force of infection in the behavior of the model.

7.   Conclusions

The proposed fractional-order model of pneumonia infection in sheep and goats has been applied in the Al-Baha region of the Kingdom of Saudi Arabia. We have obtained real data related to the pneumonia disease in sheep and goats from the Ministry of Environment, Agriculture, and Water in the Al-Baha region of the Kingdom of Saudi Arabia. And this ministry is mainly responsible for collecting this data annually and there is no authentic source that has the ability to provide such data on daily basis. It should be noted that the proposed fractional-order model understands the disease's actions more correctly than the integer-order version. The reason to introduce such the fractional order in the model and the dynamical analysis of the proposed model has illustrated in Section 5.2.

By using the Routh-Hurwitz criteria and constructing Lyapunov functions, the local and the global stability of free equilibrium and endemic equilibrium points are obtained. Our main aim is to minimize the infection and death rate due to pneumonia disease. For our proposed model, we founded that if $R_{0} < 1$, the disease-free equilibrium is locally asymptotically stable. Most importantly, the goal of this model is to estimate in the future the number of total infected, active cases, deaths as well as recoveries from pneumonia disease. This helps to decide the future decision so that we control or minimize the above cases.

By means of an effective numerical scheme, various numerical simulations were carried out in order to shed more light on the model's characteristics. Moreover, we examined the optimality system using the following strategies:

● Using just one control, preventative, treatment, or screening method.

● Using two types of control: Preventive and treatment.

● Using two types of controls: Prevention and screening.

● Using two types of controls: Treatment and screening.

● By implementing all preventative, therapeutic, and screening measures.

Furthermore, we studied cost-effectiveness to determine the least and most expensive solutions. Based on the findings of this study's pairwise analysis, we conclude that combining prevention and treatment is the most cost-effective strategy in terms of both cost and health benefits.The proposed fractional order model of pneumonia infection in sheep and goats has been applied in Al-Baha region of the Kingdom of Saudi Arabia. We have obtained real data whic related to the pneumonia disease in sheep and goats from Ministry of Environment, Agriculture and Water in Al-Baha region of the Kingdom of Saudi Arabia. And this ministry is mainly responsible for collecting this data annually and there is no authentic source has the ability to provide such data on daily basis. It should be noted that the proposed fractional order model understands the disease's actions more correctly than the integer-order version.

8.   Future research

The work presented in this paper can be extended in a lot of respects including the following:

● We note that the models in this paper are not exhaustive.

● Incorporating drug resistance compartment could be studied for Pneumonia diseases.

● An extension of the models in this paper to take care of the period between contacting disease-causing organisms and the development of clinical symptoms will be interesting. In this case, delay differential equation models will be developed.

● Analyzing the models based upon immune level dynamics of the diseases.

Acknowledgments

This research is a part of a project entitled "Mathematical Modelling and optimal control of pneumonia disease in sheep and goats in Al-Baha region with cost-effective strategies". This project was funded by the Deanship of Scientific Research, Al-Baha University, KSA (Grant No. 1442/7). The assistance of the deanship is gratefully acknowledged. Also, we would like to thank the reviewers for their careful reading of the manuscript and for the detailed comments and corrections they suggested.

Conflict of interest

The authors declare no conflicts of interest.