# A Deterministic Model for Q Fever Transmission Dynamics within Dairy Cattle Herds: Using Sensitivity Analysis and Optimal Controls

Authors: Joshua Kiddy K. Asamoah; Zhen Jin; Gui-Quan Sun; Michael Y. Li
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

This paper presents a differential equation model which describes a possible transmission route for Q fever dynamics in cattle herds. The model seeks to ascertain epidemiological and theoretical inferences in understanding how to avert an outbreak of Q fever in dairy cattle herds (livestock). To prove the stability of the model’s equilibria, we use a matrix-theoretic method and a Lyapunov function which establishes the local and global asymptotic behaviour of the model. We introduce time-dependent vaccination, environmental hygiene, and culling and then solve for optimal strategies. The optimal control strategies are necessary management practices that may increase animal health in aCoxiella burnetii-induced environment and may also reduce the transmission of the disease from livestock into the human population. The sensitivity analysis presents the relative importance of the various generic parameters in the model. We hope that the description of the results and the optimality trajectories provides some guidelines for animal owners and veterinary officers on how to effectively minimize the bacteria and control cost before/during an outbreak.


---

## Body

## 1. Introduction

Q fever is a bacterial disease which is caused byCoxiella burnetii (C. burnetii) [1, 2]. Q fever is a potential biological warfare agent being very infectious and very durable in the environment as well as capable of windborne spread [3]. C. burnetii affects a wide range of animals and also causes illness in humans [1, 2]. The bacteria are mostly found in farm animals such as goats, sheep, and cattle, but can also be found in cats, dogs, rodents, birds, and other wildlife [4, 5]. Coxiella burnetii in ruminants result in abortion, stillbirth, mastitis, infertility, premature delivery, and weak offspring [6–8]. Other clinical signs in animals may include, fever, mild coughing, anorexia, and rhinitis [9]. The bacteria are shed from an infected animal into the environment through urine, faeces, milk, and vaginal fluids, but most commonly, the bacterium is in the amniotic fluids and the placenta discharge during abortion or parturition of an infected animal [5]. Q fever is noted as the second most commonly reported laboratory infection with several recorded outbreaks involving 15 or more persons [3]. The disease dynamics in both animals and humans start primarily through inhalation of contaminated dust, contact with placenta discharge during an abortion, parturition of an infected animal, drinking unpasteurised infected milk, ingestion meat containing Coxiella burnetii, or contact with contaminated wool [3, 5, 10, 11]. Outbreaks typically occur following birth or abortion where the environment becomes contaminated with birthing fluids of an infected animal [3]. The bacteria contaminate dust and then are spread by the wind for long distances [3].Q fever is a reemerging zoonosis in most parts of Europe, which has seen a sharp rise in recent years, especially in the Netherlands with a large number of human cases, attributed to livestock [12–14]. An outbreak of Q fever in the Netherlands between 2007 and 2010 resulted in about 4000 reported human cases, and efforts to end the epidemic resulted in the culling of more than 50,000 small ruminants and a temporary ban on breeding of animals [10]. A survey indicates that 82% of cows in some California dairies were seropositive, as well as 78% of coyotes, 55% of foxes, 53% of brush rabbits, and 22% of deer in Northern California [9]. Australia meat industry makes about $1 million losses annually as a result of Q fever [15]. Geographically, Q fever infection in animals has been detected in most countries, except, New Zealand, Iceland, French Polynesia, and Norway [3, 10].In recent time, Q fever has received an increase in publicity on the potential source of transmission and possible preventative measures in both livestock and humans (see, for example, [11, 16] and the references therein). The research by De Rooij et al. [17] suggested that during any future epidemics, attention should be placed on a rapid source of identification, quantification of emission, accurate data collection, and smooth data exchange amongst relevant actors to enable effective risk assessment and risk management. Mori and Roest [18] pointed out that in achieving an effective control and management, there should be collective expertise from other scientific disciplines and stakeholder so that the expectations of farmers and the wider community can be met on the spread of the disease. Courcoul et al. [19] used the Bayesian approach to obtain key epidemiological parameters from field data and suggested that the values of the parameters obtained can be used to provide information for calibrating simulation models to assess control strategies for C. burnetii. Courcoul et al. [20] further proposed a model to assess the effectiveness of vaccination in dairy cattle and suggested that their model should be modified to simulate various control strategies such as environmental and hygienic measures.In this sense, our goal is to use the parameter values presented in [19, 20] and some assumptions to study the impact of vaccination, environmental hygiene (cleaning/burying of placenta discharge after birth or abortion), and culling/isolation using ordinary differential equation model, sensitivity analysis, and optimal control theory. Thus, we systematically analysed the environmental transmission of the disease and then obtained mathematical and epidemiological properties, such as basic reproduction number, equilibrium points, local and global stability analysis, sensitivity analysis, and optimality conditions. The optimal control introduced is to obtain optimal trajectories that depict the most effective control measures and also account for the costs involved. Optimal control is a useful mathematical tool which has been recently used to determine optimal strategies for other infectious diseases, although not for Q fever.The rest of the paper is arranged as follows: Section2 presents a vaccination-induced model and mathematical analysis, which establishes the stability of the proposed model; in Section 3, we study the local and global sensitivity analysis of the model’s parameters; in Section 4, we set up an optimal control problem; finally, in Section 5, we draw conclusions from this paper and assess what new insights this work gives to the body of knowledge on Q fever transmission dynamics and the way forward.


## 2. Basic Model for Indirect Transmission

The model is divided into susceptibleSt, asymptomatic (expose) At, vaccinated Vt, and symptomatic It cattle, respectively, with an assumed bacteria- (C. burnetii-) induced environment Bt, where t represents time. The rate at which new cattle enter the susceptible population is denoted as Λ, μ is the constant rate of death, and thus 1/μ is the average lifetime. The asymptomatic cattle become symptomatic at a constant rate α, and so 1/α is the average asymptomatic period. The susceptible cattle get vaccinated at the rate ν, θ is the rate at which vaccinated cattle lose immunity and regain susceptibility. The symptomatic cattle are affected by an additional exit which may serve as a control measure at a constant rate c (culling/isolation rate). The rate of natural decay of the bacteria from the environment is denoted by δ and ε represents the rate of environmental hygiene, thus cleaning/burying of placenta discharge after birth or abortion. Asymptomatic cattle who develop temporal resistance to the bacteria become susceptible again at the rate ρ, if ρ=0 means no temporal resistance to the bacteria, and ρ=1 means 100% resistant to the bacteria (thus I=0), and hence throughout this paper, we assume 0≤ρ<1. We assume a mass action for the transmission of the disease, and thus susceptible cattle become asymptomatic at a proportional rate βS, where β is the effective rate of contracting the bacteria through indirect means (thus environmental transmission rate through inhaling contaminated dust/through grazing). The rate of inflow of the bacteria into the environment by both asymptomatic and symptomatic cattle is η (thus an assumed shedding rate). Here, the total cattle population for the vaccination-induced model is denoted as N=S+V+A+I. The following assumptions govern the vaccination-induced model:(i)
The parametersΛ,β,η,α,ν,θ,ρ,μ,c,ε,δ are all positive, and thus Λ,μ,η,α,δ,β>0 and c,ρ,ε,θ,ν≥0(ii)
Newly purchased or new birth goes into the susceptible class only(iii)
No permanent recovery(iv)
A successfully vaccinated animal stays immune for a period of timeFigure1 depicts the compartmental diagrams of the disease in a community with control measures; model (1) gives the dynamical nature of the individual compartments with respect to time.Figure 1
A flowchart of Q fever transmission dynamics within dairy cattle herds. The light blue box = susceptible cattleS; red box = asymptomatic cattle A; green box = symptomatic cattle I; white box = vaccinated cattle V; the presence of Coxiella burnetii in the environment B is encircled in blue.


### 2.1. Basic Model and Analysis

(1)
d
S
d
t
=
Λ
−
β
B
S
−
μ
+
ν
S
+
θ
V
+
ρ
A
,
d
V
d
t
=
ν
S
−
θ
+
μ
V
,
d
A
d
t
=
β
B
S
−
μ
+
α
+
ρ
A
,
d
I
d
t
=
α
A
−
c
+
μ
I
,
d
B
d
t
=
η
A
+
I
−
ε
+
δ
B
,with initial conditions(2)S0=S0>0,V0=V0>0,A0=A0≥0,I0=I0≥0,B0=B0≥0,and total differential population(3)dNdt=Λ−μN−cI,dBdt=ηA+I−ε+δB.


### 2.2. Positivity and Feasible Region

Lemma 1.
For any nonnegative initial data,(4)S0=S0>0,V0=V0>0,A0=A0>0,I0=I0>0,B0=B0>0,model (1) has a unique global solution for all t≥0.Lemma 2.
From the initial conditions and for epidemiological meaningfulness, we assume that all parameters are positive. Hence, the proposed model (1) will then be analysed in a feasible region given as(5)Ω=S,V,A,I,B∈ℛ+5Nt≤Λμ,Bt≤ηΛε+δμ.Now, we need to show that the regionΩ is positively invariant, to support the assertion that the model is epidemiologically useful and that all its state variables are nonnegative as time increases and stays in Ω.Proof.
LetX=S,V,A,I,BT where T denotes transposition, also if we replace βB with M, then model (1) can be written as dX/dt=DX+Y, where(6)D=−M+μ+νθρ00ν−θ+μ000M0−μ+α+ρ0000α−c+μ000ηη−ε+δ,Y=Λ0000.
Now, from the second equation in model (1), we have(7)dVdt=νS−θ+μV,applying the integration factor method and change of variables, we have(8)Vt=V0e−θ+μt+e−θ+μtν∫0tSpeθ+μpdp=e−θ+μtV0+ν∫0tSpeθ+μpdp.
From the third equation of model (1), we also obtain(9)At=e−μ+α+ρtA0+β∫0tBpSpeμ+α+ρpdp.
It follows from standard property thatdV/dtt=t0≥0, dA/dtt=t0≥0, and dS/dtt=t0≥0, provided St0=0, with similar properties being valid for It, Bt, which ensures that the state variable remains positive during the entire scope of the study (see [21] and Appendix A of [22]). Based on this, we see that matrix D has all its off-diagonal entries to be nonnegative and that matrix Y≥0, which proves the properties of Metzler matrix [23]. Therefore, it implies that model (1) is positively invariant in ℝ+5. Finally, using equation (3), it can be shown that all solution sets of model (1) remain in the feasible region Ω, with the notation that 0<A+I≤N and S+A+I+V=N. Hence, the proposed model is considered to be epidemiologically meaningful (mathematically well posed) and it is sufficient to say that all solutions of the Q fever model with initial conditions remain in the feasible region Ω.


### 2.3. Equilibrium Points and Basic Reproduction Number

The equilibrium points are the constant solutions to the differential equations in model (1), and these equilibrium points and the basic reproduction number play an important role in the long-term behaviour of the solution to the epidemiological model (1). The equilibrium points of model (1) are obtained by setting the right-hand side of model (1) to zero, which gives(10)Λ−βBS−μ+νS+θV+ρA=0,(11)νS−θ+μV=0,(12)βBS−μ+α+ρA=0,(13)αA−c+μI=0,(14)ηA+I−ε+δB=0.The two main equilibria considered here are the disease-free equilibrium and endemic equilibrium. We obtained the disease-free equilibrium by assuming that the disease-induced statesA=I=0 and a free C. burnetii environment and thus B=0. Therefore, by solving for the nonzero state variables in equations (10) and (11), a unique disease-free equilibrium is obtained as follows:(15)E0≔S0,V0,A0,I0,B0=Λμ+θμμ+ν+θ,Λνμμ+ν+θ,0,0,0.


#### 2.3.1. Basic Reproduction NumbersℛC and ℛ0

The basic reproduction numberℛ0 is obtained by using the disease-free equilibrium points in terms of the original parameters of the model; epidemiologically, the reproduction number informs us on the number of secondary infections produced when a single infectious individual is introduced into a completely susceptible population. Here, we define the basic reproduction number as the number of secondary infections produced in a completely susceptible cattle population when the organism is released into the environment. Also, the analytic nature of the basic reproduction number gives clues on parameters in the model that contributes to the spread of the disease. Hence, with the aid of the next-generation matrix FV−1 as introduced by [24, 25] and the disease-free equilibrium (15), the epidemiological reproduction numbers for model (1) are as follows:(16)ℛC=βηΛθ+μα+c+μμμ+α+ρθ+μ+νε+δc+μ,(17)ℛ0=βηΛα+μμ2μ+α+ρδ,where equation (16) is called the control reproduction number ℛC and equation (17) is the basic reproduction number without controls ℛ0 (thus c=0,θ=0,ν=0,ε=0). From the basic reproduction numbers, it is evidential that there exists a linear relationship between the rate of environmental transmission β and shedding rate η, and hence, this shows that an increase in bacterial load in the environment may have a corresponding increase in the transmission process. This assertion has been numerically varied using Figure 2(b) in Section 3.2.1.Figure 2
(a) The scatter diagrams for some selected parameters inℛC (scatter plots of β,Λ,η,ρ,α,θ,μ, and ν in ℛC). (b) Perfect rank relations of the parameters in ℛC (Tornado plot of parameters in ℛC), (c) The coupling effect of both vaccination and environmental hygiene on ℛC (3D plot of ε and v ℛC).
(a)
(b)
(c)


#### 2.3.2. Endemic Equilibrium Point

Now, takingS,V,A,I,B>0 in the feasible set Ω, and using equations (10)–(14) together with (16), we obtain the following endemic equilibrium point E∗≔S∗,V∗,A∗,I∗,B∗. Thus(18)S∗=μ+α+ρε+δc+μβηα+c+μ,V∗=S∗νθ+μ,A∗=ℛC−1βηθ+μμ+αα+c+μ,I∗=αℛC−1βηc+μθ+μμ+αα+c+μ,B∗=ℛC−1βμ+cε+δθ+μμ+α.The endemic equilibrium has a unique equilibrium point if the vaccination-induced reproduction number is greater than unityℛC>1. In the next subsection, we investigate the stability of the equilibrium points.


### 2.4. Stability Analysis

In this section, we focus on establishing the local and global stability of the model equilibria. To obtain the local stability of the disease-free equilibrium, we constructed a linearized Jacobian matrixJ1, evaluated at the disease-free equilibrium, which gives(19)J1=−μ+α+ρ0βΛθ+μμθ+μ+να−c+μ0ηη−ε+δ,where J1=F−V, J1 is irreducible and nonnegative. The simple eigenvalue of J1 is denoted as sJ1 which is obtained from the relation sJ1=maxReλ:Reλis the real part of eigenvalue of matrixJ1. The following equivalence holds:(20)ℛC>1⟺sJ1>0,ℛC<1⟺sJ1<0,see,Theorem2of24Section3.Using Theorem2 of van den Driessche and Watmough [24], we see that assumptions A1−A4 are satisfied, and also A5 is satisfied if all the eigenvalues of matrix JE0 have negative real parts, and thus(21)JE0=J10J3J4,where J4=−F, and J4=−μ+νθνθ+μ.The simple eigenvalue ofJ4 is sJ4=max−μ+ν+θ,−μ. Hence, model (1) fully satisfies the assumptions A1−A5 of van den Driessche and Watmough [24], and therefore, the disease-free equilibrium of the Q fever model is locally asymptotically stable if ℛC<1 and unstable if ℛC>1. The below theorem addresses the global stability of the disease-free equilibrium.Theorem 1.
IfℛC≤1, then the disease-free equilibrium of model (1) is globally asymptotically stable.Proof.
Applying Theorem 2.1 of the matrix-theoretic method in [26], the global stability of the disease-free equilibrium E0 can be established by constructing a Lyapunov function Φ=wTV−1x, since V−1F is reducible (the third column is the only nonzero column), wT is the left Perron eigenvector of the matrix V−1F, and x=A,I,BT is the disease-induced classes. Algebraic operation of V−1F gives(22)V−1F=00βS0μ+α+ρ00αβS0μ+α+ρc+μ00βηS0α+c+μμ+α+ρc+με+δ,where S0=Λμ+θ/μμ+ν+θ. Therefore, by Theorem 5.1 of [26], the left Perron eigenvector wT of V−1F for the largest eigenvalue is defined as(23)0,0,w3V−1F=ℛC0,0,w3.
This gives(24)ℛCw3=w3βηS0α+c+μμ+α+ρc+με+δ,which implies that w3=1; therefore, wT=0,0,1. Now simplifying Φ=wTV−1x yields(25)Φ=ηα+c+μA+ημ+α+ρI+μ+α+ρc+μBμ+α+ρc+με+δ,=ℛCβS0A+μ+α+ρα+c+μI+μ+α+ρc+μηα+c+μB.
Using [26], the derivative of Φ along the solution of model (1) is illustrated as follows:(26)dΦdt=ℛC−1wTx−wTV−1fx,S,where fx,S≔F−Vx−ℱx,S+Vx,S. Therefore, after some algebra, we have(27)fx,S=βBS0−S00,forS0≥S.
Now using equation (26) and some algebraic simplification leads to(28)dΦdt=ℛC−1B−ℛCS0BS0−S.
Therefore,dΦ/dt≤0 provided ℛC≤1 and S0≥S. Furthermore, dΦ/dt≤0 implies that ℛC=1 and B=0 or S=S0. It can be shown that the only invariant set where B=0 or S=S0 is a singleton of the set S0,0,0,0. Since the four state variables used in constructing the Lyapunov functional lead to global stability, it then means that Vt⟶V0 as t⟶∞. Hence, by [27] (LaSalle’s invariance principle), the disease-free equilibrium for the vaccination model is globally asymptotically stable. The ramification of this global stability in the presence of vaccination implies that no matter how many initial asymptomatic and symptomatic cattle in the population, the disease cannot be established (persist) and eventually dies out. It also indicates that no matter how many initial bacterial loads in the environment, the disease will die out from the population in the long round.


#### 2.4.1. Global Stability of Endemic Equilibrium

The nature of the endemic equilibriumE∗ indicates that for ℛC>1, then there exists a unique endemic equilibrium resulting in an unstable disease-free equilibrium. To establish the global stability of E∗, here we explore a special case (simpler case) of the model (1), where ρ=0. Let Ω0=S,V,A,I∈Ω:A=I=0 be the stable manifold of E0. Therefore, we claim the following result.Theorem 2.
The unique endemic equilibrium of the model (1) is globally asymptotically stable in Ω/Ω0 whenever ℛC>1 and ρ=0.Proof.
At the endemic equilibrium state of model (1) with ρ=0, the following equalities(29)Λ−βB∗S∗−μ+νS∗+θV∗=0,(30)νS∗−θ+μV∗=0,(31)βB∗S∗−μ+αA∗=0,(32)αA∗−c+μI∗=0,(33)ηA∗+I∗−ε+δB∗=0,hold, which further leads to(34)μ+ν=Λ−βB∗S∗+θV∗S∗,θ+μ=νS∗V∗,μ+α=βB∗S∗A∗,c+μ=αA∗I∗,ε+δ=ηA∗+I∗B∗.
Substituting equation (34) into model (1) result in model (35), in which its stability is the same as model (1):(35)dSdt=SΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1,dVdt=νS∗VV∗SV∗VS∗−1,dAdt=βB∗S∗AA∗BA∗SB∗S∗A−1,dIdt=αA∗II∗AI∗IA∗−1,dBdt=ηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Now, we define a Lyapunov function as(36)L1=S−S∗+S∗lnSS∗+V−V∗+V∗lnVV∗+A−A∗+A∗lnAA∗+k1I−I∗+I∗lnII∗+k2B−B∗+B∗lnBB∗,where k1 and k2 are positive constants to be determined, and thus k1>0,k2>0. The time derivative of the Lyapunov function along with the solution of model (35) yields the following equation:(37)dL1dt=1−S∗SSΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1+1−V∗VνS∗VV∗SV∗VS∗−1+1−A∗AβB∗S∗AA∗BA∗SB∗S∗A−1+k11−I∗IαA∗II∗AI∗IA∗−1+k21−B∗BηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Simplifying the above expression and rearranging lead to(38)dL1dt=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
We now let(39)L1=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗,(40)L2=−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
MakingΛ and ν the subject in equations (29) and (30), respectively, and substituting into (39) with series of simplification give(41)L1=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗V.
Further elaboration onL2 results in(42)L2=k2ηA∗1+AA∗−BB∗−AB∗BA∗+I∗1+II∗−BB∗−B∗IBI∗+αk1A∗AA∗−II∗−AI∗IA∗+1+βB∗S∗BB∗−AA∗+SS∗−BSA∗B∗S∗A.
To obtaink1,k2, we make the following assertion by equating the constant coefficients of the terms A/A∗,I/I∗ and B/B∗ in L2 to zero, thus(43)k2ηA∗+αk1A∗−βB∗S∗=0,(44)k2ηI∗−αk1A∗=0,(45)βB∗S∗−k2ηA∗+I∗=0.
Solving (44) and (45), we obtain k1=βB∗S∗I∗/αA∗A∗+I∗ and k2=βB∗S∗/ηA∗+I∗. Replacing k1,k2 into (42) and simplifying yield(46)L2=βB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Hence, adding equations (41) and (46), we have(47)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A,(48)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗D0+βB∗S∗A∗A∗+I∗D1+βBS∗D2,where(49)D0=2+AA∗−BB∗−B∗IBI∗−AI∗A∗I,D1=1+AA∗−BB∗−B∗ABA∗,D2=1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Thus, ifD0≤0,D1≤0,D2≤0, and from the fact that the arithmetic mean is greater than or equal to the geometric mean, then it follows that dL1/dt≤0, provided the coefficients μS∗,μV∗,θV∗,βB∗S∗I∗/A∗+I∗,βB∗S∗A∗/A∗+I∗,βBS∗ are nonnegative. With the claim that(50)wx=x−1−lnx,it leads to(51)D0=wAA∗−wBB∗−wB∗IBI∗−wAI∗A∗I,D1=wAA∗−wBB∗−wB∗ABA∗,D2=wB∗SBS∗−wB∗AA∗B−wSA∗S∗A.
To verify the above expression, the definition ofwx and properties of logarithms were used to simplify the expression on the right hand of Di,i=0,1,2. Then applying the proposition A1 of [28] gives(52)D0≤−wBB∗−wB∗IBI∗−wAI∗A∗I≤0,D1≤−wBB∗−wB∗ABA∗≤0,D2≤−wB∗AA∗B−wSA∗S∗A≤0.
It is evidential thatdL1S,V,A,I,B/dt=0 only holds when S=S∗,V=V∗,A=A∗,I=I∗,B=B∗ and that E∗ is the only equilibrium state of model (1) on its line. Therefore, by Lyapunov–LaSalle asymptotic stability theorem [27], the positive equilibrium state E∗ is globally asymptotically stable in the positive region ℝ+5 whenever ℛC>1 and ρ=0.


## 2.1. Basic Model and Analysis

(1)
d
S
d
t
=
Λ
−
β
B
S
−
μ
+
ν
S
+
θ
V
+
ρ
A
,
d
V
d
t
=
ν
S
−
θ
+
μ
V
,
d
A
d
t
=
β
B
S
−
μ
+
α
+
ρ
A
,
d
I
d
t
=
α
A
−
c
+
μ
I
,
d
B
d
t
=
η
A
+
I
−
ε
+
δ
B
,with initial conditions(2)S0=S0>0,V0=V0>0,A0=A0≥0,I0=I0≥0,B0=B0≥0,and total differential population(3)dNdt=Λ−μN−cI,dBdt=ηA+I−ε+δB.


## 2.2. Positivity and Feasible Region

Lemma 1.
For any nonnegative initial data,(4)S0=S0>0,V0=V0>0,A0=A0>0,I0=I0>0,B0=B0>0,model (1) has a unique global solution for all t≥0.Lemma 2.
From the initial conditions and for epidemiological meaningfulness, we assume that all parameters are positive. Hence, the proposed model (1) will then be analysed in a feasible region given as(5)Ω=S,V,A,I,B∈ℛ+5Nt≤Λμ,Bt≤ηΛε+δμ.Now, we need to show that the regionΩ is positively invariant, to support the assertion that the model is epidemiologically useful and that all its state variables are nonnegative as time increases and stays in Ω.Proof.
LetX=S,V,A,I,BT where T denotes transposition, also if we replace βB with M, then model (1) can be written as dX/dt=DX+Y, where(6)D=−M+μ+νθρ00ν−θ+μ000M0−μ+α+ρ0000α−c+μ000ηη−ε+δ,Y=Λ0000.
Now, from the second equation in model (1), we have(7)dVdt=νS−θ+μV,applying the integration factor method and change of variables, we have(8)Vt=V0e−θ+μt+e−θ+μtν∫0tSpeθ+μpdp=e−θ+μtV0+ν∫0tSpeθ+μpdp.
From the third equation of model (1), we also obtain(9)At=e−μ+α+ρtA0+β∫0tBpSpeμ+α+ρpdp.
It follows from standard property thatdV/dtt=t0≥0, dA/dtt=t0≥0, and dS/dtt=t0≥0, provided St0=0, with similar properties being valid for It, Bt, which ensures that the state variable remains positive during the entire scope of the study (see [21] and Appendix A of [22]). Based on this, we see that matrix D has all its off-diagonal entries to be nonnegative and that matrix Y≥0, which proves the properties of Metzler matrix [23]. Therefore, it implies that model (1) is positively invariant in ℝ+5. Finally, using equation (3), it can be shown that all solution sets of model (1) remain in the feasible region Ω, with the notation that 0<A+I≤N and S+A+I+V=N. Hence, the proposed model is considered to be epidemiologically meaningful (mathematically well posed) and it is sufficient to say that all solutions of the Q fever model with initial conditions remain in the feasible region Ω.


## 2.3. Equilibrium Points and Basic Reproduction Number

The equilibrium points are the constant solutions to the differential equations in model (1), and these equilibrium points and the basic reproduction number play an important role in the long-term behaviour of the solution to the epidemiological model (1). The equilibrium points of model (1) are obtained by setting the right-hand side of model (1) to zero, which gives(10)Λ−βBS−μ+νS+θV+ρA=0,(11)νS−θ+μV=0,(12)βBS−μ+α+ρA=0,(13)αA−c+μI=0,(14)ηA+I−ε+δB=0.The two main equilibria considered here are the disease-free equilibrium and endemic equilibrium. We obtained the disease-free equilibrium by assuming that the disease-induced statesA=I=0 and a free C. burnetii environment and thus B=0. Therefore, by solving for the nonzero state variables in equations (10) and (11), a unique disease-free equilibrium is obtained as follows:(15)E0≔S0,V0,A0,I0,B0=Λμ+θμμ+ν+θ,Λνμμ+ν+θ,0,0,0.


### 2.3.1. Basic Reproduction NumbersℛC and ℛ0

The basic reproduction numberℛ0 is obtained by using the disease-free equilibrium points in terms of the original parameters of the model; epidemiologically, the reproduction number informs us on the number of secondary infections produced when a single infectious individual is introduced into a completely susceptible population. Here, we define the basic reproduction number as the number of secondary infections produced in a completely susceptible cattle population when the organism is released into the environment. Also, the analytic nature of the basic reproduction number gives clues on parameters in the model that contributes to the spread of the disease. Hence, with the aid of the next-generation matrix FV−1 as introduced by [24, 25] and the disease-free equilibrium (15), the epidemiological reproduction numbers for model (1) are as follows:(16)ℛC=βηΛθ+μα+c+μμμ+α+ρθ+μ+νε+δc+μ,(17)ℛ0=βηΛα+μμ2μ+α+ρδ,where equation (16) is called the control reproduction number ℛC and equation (17) is the basic reproduction number without controls ℛ0 (thus c=0,θ=0,ν=0,ε=0). From the basic reproduction numbers, it is evidential that there exists a linear relationship between the rate of environmental transmission β and shedding rate η, and hence, this shows that an increase in bacterial load in the environment may have a corresponding increase in the transmission process. This assertion has been numerically varied using Figure 2(b) in Section 3.2.1.Figure 2
(a) The scatter diagrams for some selected parameters inℛC (scatter plots of β,Λ,η,ρ,α,θ,μ, and ν in ℛC). (b) Perfect rank relations of the parameters in ℛC (Tornado plot of parameters in ℛC), (c) The coupling effect of both vaccination and environmental hygiene on ℛC (3D plot of ε and v ℛC).
(a)
(b)
(c)


### 2.3.2. Endemic Equilibrium Point

Now, takingS,V,A,I,B>0 in the feasible set Ω, and using equations (10)–(14) together with (16), we obtain the following endemic equilibrium point E∗≔S∗,V∗,A∗,I∗,B∗. Thus(18)S∗=μ+α+ρε+δc+μβηα+c+μ,V∗=S∗νθ+μ,A∗=ℛC−1βηθ+μμ+αα+c+μ,I∗=αℛC−1βηc+μθ+μμ+αα+c+μ,B∗=ℛC−1βμ+cε+δθ+μμ+α.The endemic equilibrium has a unique equilibrium point if the vaccination-induced reproduction number is greater than unityℛC>1. In the next subsection, we investigate the stability of the equilibrium points.


## 2.3.1. Basic Reproduction NumbersℛC and ℛ0

The basic reproduction numberℛ0 is obtained by using the disease-free equilibrium points in terms of the original parameters of the model; epidemiologically, the reproduction number informs us on the number of secondary infections produced when a single infectious individual is introduced into a completely susceptible population. Here, we define the basic reproduction number as the number of secondary infections produced in a completely susceptible cattle population when the organism is released into the environment. Also, the analytic nature of the basic reproduction number gives clues on parameters in the model that contributes to the spread of the disease. Hence, with the aid of the next-generation matrix FV−1 as introduced by [24, 25] and the disease-free equilibrium (15), the epidemiological reproduction numbers for model (1) are as follows:(16)ℛC=βηΛθ+μα+c+μμμ+α+ρθ+μ+νε+δc+μ,(17)ℛ0=βηΛα+μμ2μ+α+ρδ,where equation (16) is called the control reproduction number ℛC and equation (17) is the basic reproduction number without controls ℛ0 (thus c=0,θ=0,ν=0,ε=0). From the basic reproduction numbers, it is evidential that there exists a linear relationship between the rate of environmental transmission β and shedding rate η, and hence, this shows that an increase in bacterial load in the environment may have a corresponding increase in the transmission process. This assertion has been numerically varied using Figure 2(b) in Section 3.2.1.Figure 2
(a) The scatter diagrams for some selected parameters inℛC (scatter plots of β,Λ,η,ρ,α,θ,μ, and ν in ℛC). (b) Perfect rank relations of the parameters in ℛC (Tornado plot of parameters in ℛC), (c) The coupling effect of both vaccination and environmental hygiene on ℛC (3D plot of ε and v ℛC).
(a)
(b)
(c)


## 2.3.2. Endemic Equilibrium Point

Now, takingS,V,A,I,B>0 in the feasible set Ω, and using equations (10)–(14) together with (16), we obtain the following endemic equilibrium point E∗≔S∗,V∗,A∗,I∗,B∗. Thus(18)S∗=μ+α+ρε+δc+μβηα+c+μ,V∗=S∗νθ+μ,A∗=ℛC−1βηθ+μμ+αα+c+μ,I∗=αℛC−1βηc+μθ+μμ+αα+c+μ,B∗=ℛC−1βμ+cε+δθ+μμ+α.The endemic equilibrium has a unique equilibrium point if the vaccination-induced reproduction number is greater than unityℛC>1. In the next subsection, we investigate the stability of the equilibrium points.


## 2.4. Stability Analysis

In this section, we focus on establishing the local and global stability of the model equilibria. To obtain the local stability of the disease-free equilibrium, we constructed a linearized Jacobian matrixJ1, evaluated at the disease-free equilibrium, which gives(19)J1=−μ+α+ρ0βΛθ+μμθ+μ+να−c+μ0ηη−ε+δ,where J1=F−V, J1 is irreducible and nonnegative. The simple eigenvalue of J1 is denoted as sJ1 which is obtained from the relation sJ1=maxReλ:Reλis the real part of eigenvalue of matrixJ1. The following equivalence holds:(20)ℛC>1⟺sJ1>0,ℛC<1⟺sJ1<0,see,Theorem2of24Section3.Using Theorem2 of van den Driessche and Watmough [24], we see that assumptions A1−A4 are satisfied, and also A5 is satisfied if all the eigenvalues of matrix JE0 have negative real parts, and thus(21)JE0=J10J3J4,where J4=−F, and J4=−μ+νθνθ+μ.The simple eigenvalue ofJ4 is sJ4=max−μ+ν+θ,−μ. Hence, model (1) fully satisfies the assumptions A1−A5 of van den Driessche and Watmough [24], and therefore, the disease-free equilibrium of the Q fever model is locally asymptotically stable if ℛC<1 and unstable if ℛC>1. The below theorem addresses the global stability of the disease-free equilibrium.Theorem 1.
IfℛC≤1, then the disease-free equilibrium of model (1) is globally asymptotically stable.Proof.
Applying Theorem 2.1 of the matrix-theoretic method in [26], the global stability of the disease-free equilibrium E0 can be established by constructing a Lyapunov function Φ=wTV−1x, since V−1F is reducible (the third column is the only nonzero column), wT is the left Perron eigenvector of the matrix V−1F, and x=A,I,BT is the disease-induced classes. Algebraic operation of V−1F gives(22)V−1F=00βS0μ+α+ρ00αβS0μ+α+ρc+μ00βηS0α+c+μμ+α+ρc+με+δ,where S0=Λμ+θ/μμ+ν+θ. Therefore, by Theorem 5.1 of [26], the left Perron eigenvector wT of V−1F for the largest eigenvalue is defined as(23)0,0,w3V−1F=ℛC0,0,w3.
This gives(24)ℛCw3=w3βηS0α+c+μμ+α+ρc+με+δ,which implies that w3=1; therefore, wT=0,0,1. Now simplifying Φ=wTV−1x yields(25)Φ=ηα+c+μA+ημ+α+ρI+μ+α+ρc+μBμ+α+ρc+με+δ,=ℛCβS0A+μ+α+ρα+c+μI+μ+α+ρc+μηα+c+μB.
Using [26], the derivative of Φ along the solution of model (1) is illustrated as follows:(26)dΦdt=ℛC−1wTx−wTV−1fx,S,where fx,S≔F−Vx−ℱx,S+Vx,S. Therefore, after some algebra, we have(27)fx,S=βBS0−S00,forS0≥S.
Now using equation (26) and some algebraic simplification leads to(28)dΦdt=ℛC−1B−ℛCS0BS0−S.
Therefore,dΦ/dt≤0 provided ℛC≤1 and S0≥S. Furthermore, dΦ/dt≤0 implies that ℛC=1 and B=0 or S=S0. It can be shown that the only invariant set where B=0 or S=S0 is a singleton of the set S0,0,0,0. Since the four state variables used in constructing the Lyapunov functional lead to global stability, it then means that Vt⟶V0 as t⟶∞. Hence, by [27] (LaSalle’s invariance principle), the disease-free equilibrium for the vaccination model is globally asymptotically stable. The ramification of this global stability in the presence of vaccination implies that no matter how many initial asymptomatic and symptomatic cattle in the population, the disease cannot be established (persist) and eventually dies out. It also indicates that no matter how many initial bacterial loads in the environment, the disease will die out from the population in the long round.


### 2.4.1. Global Stability of Endemic Equilibrium

The nature of the endemic equilibriumE∗ indicates that for ℛC>1, then there exists a unique endemic equilibrium resulting in an unstable disease-free equilibrium. To establish the global stability of E∗, here we explore a special case (simpler case) of the model (1), where ρ=0. Let Ω0=S,V,A,I∈Ω:A=I=0 be the stable manifold of E0. Therefore, we claim the following result.Theorem 2.
The unique endemic equilibrium of the model (1) is globally asymptotically stable in Ω/Ω0 whenever ℛC>1 and ρ=0.Proof.
At the endemic equilibrium state of model (1) with ρ=0, the following equalities(29)Λ−βB∗S∗−μ+νS∗+θV∗=0,(30)νS∗−θ+μV∗=0,(31)βB∗S∗−μ+αA∗=0,(32)αA∗−c+μI∗=0,(33)ηA∗+I∗−ε+δB∗=0,hold, which further leads to(34)μ+ν=Λ−βB∗S∗+θV∗S∗,θ+μ=νS∗V∗,μ+α=βB∗S∗A∗,c+μ=αA∗I∗,ε+δ=ηA∗+I∗B∗.
Substituting equation (34) into model (1) result in model (35), in which its stability is the same as model (1):(35)dSdt=SΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1,dVdt=νS∗VV∗SV∗VS∗−1,dAdt=βB∗S∗AA∗BA∗SB∗S∗A−1,dIdt=αA∗II∗AI∗IA∗−1,dBdt=ηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Now, we define a Lyapunov function as(36)L1=S−S∗+S∗lnSS∗+V−V∗+V∗lnVV∗+A−A∗+A∗lnAA∗+k1I−I∗+I∗lnII∗+k2B−B∗+B∗lnBB∗,where k1 and k2 are positive constants to be determined, and thus k1>0,k2>0. The time derivative of the Lyapunov function along with the solution of model (35) yields the following equation:(37)dL1dt=1−S∗SSΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1+1−V∗VνS∗VV∗SV∗VS∗−1+1−A∗AβB∗S∗AA∗BA∗SB∗S∗A−1+k11−I∗IαA∗II∗AI∗IA∗−1+k21−B∗BηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Simplifying the above expression and rearranging lead to(38)dL1dt=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
We now let(39)L1=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗,(40)L2=−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
MakingΛ and ν the subject in equations (29) and (30), respectively, and substituting into (39) with series of simplification give(41)L1=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗V.
Further elaboration onL2 results in(42)L2=k2ηA∗1+AA∗−BB∗−AB∗BA∗+I∗1+II∗−BB∗−B∗IBI∗+αk1A∗AA∗−II∗−AI∗IA∗+1+βB∗S∗BB∗−AA∗+SS∗−BSA∗B∗S∗A.
To obtaink1,k2, we make the following assertion by equating the constant coefficients of the terms A/A∗,I/I∗ and B/B∗ in L2 to zero, thus(43)k2ηA∗+αk1A∗−βB∗S∗=0,(44)k2ηI∗−αk1A∗=0,(45)βB∗S∗−k2ηA∗+I∗=0.
Solving (44) and (45), we obtain k1=βB∗S∗I∗/αA∗A∗+I∗ and k2=βB∗S∗/ηA∗+I∗. Replacing k1,k2 into (42) and simplifying yield(46)L2=βB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Hence, adding equations (41) and (46), we have(47)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A,(48)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗D0+βB∗S∗A∗A∗+I∗D1+βBS∗D2,where(49)D0=2+AA∗−BB∗−B∗IBI∗−AI∗A∗I,D1=1+AA∗−BB∗−B∗ABA∗,D2=1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Thus, ifD0≤0,D1≤0,D2≤0, and from the fact that the arithmetic mean is greater than or equal to the geometric mean, then it follows that dL1/dt≤0, provided the coefficients μS∗,μV∗,θV∗,βB∗S∗I∗/A∗+I∗,βB∗S∗A∗/A∗+I∗,βBS∗ are nonnegative. With the claim that(50)wx=x−1−lnx,it leads to(51)D0=wAA∗−wBB∗−wB∗IBI∗−wAI∗A∗I,D1=wAA∗−wBB∗−wB∗ABA∗,D2=wB∗SBS∗−wB∗AA∗B−wSA∗S∗A.
To verify the above expression, the definition ofwx and properties of logarithms were used to simplify the expression on the right hand of Di,i=0,1,2. Then applying the proposition A1 of [28] gives(52)D0≤−wBB∗−wB∗IBI∗−wAI∗A∗I≤0,D1≤−wBB∗−wB∗ABA∗≤0,D2≤−wB∗AA∗B−wSA∗S∗A≤0.
It is evidential thatdL1S,V,A,I,B/dt=0 only holds when S=S∗,V=V∗,A=A∗,I=I∗,B=B∗ and that E∗ is the only equilibrium state of model (1) on its line. Therefore, by Lyapunov–LaSalle asymptotic stability theorem [27], the positive equilibrium state E∗ is globally asymptotically stable in the positive region ℝ+5 whenever ℛC>1 and ρ=0.


## 2.4.1. Global Stability of Endemic Equilibrium

The nature of the endemic equilibriumE∗ indicates that for ℛC>1, then there exists a unique endemic equilibrium resulting in an unstable disease-free equilibrium. To establish the global stability of E∗, here we explore a special case (simpler case) of the model (1), where ρ=0. Let Ω0=S,V,A,I∈Ω:A=I=0 be the stable manifold of E0. Therefore, we claim the following result.Theorem 2.
The unique endemic equilibrium of the model (1) is globally asymptotically stable in Ω/Ω0 whenever ℛC>1 and ρ=0.Proof.
At the endemic equilibrium state of model (1) with ρ=0, the following equalities(29)Λ−βB∗S∗−μ+νS∗+θV∗=0,(30)νS∗−θ+μV∗=0,(31)βB∗S∗−μ+αA∗=0,(32)αA∗−c+μI∗=0,(33)ηA∗+I∗−ε+δB∗=0,hold, which further leads to(34)μ+ν=Λ−βB∗S∗+θV∗S∗,θ+μ=νS∗V∗,μ+α=βB∗S∗A∗,c+μ=αA∗I∗,ε+δ=ηA∗+I∗B∗.
Substituting equation (34) into model (1) result in model (35), in which its stability is the same as model (1):(35)dSdt=SΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1,dVdt=νS∗VV∗SV∗VS∗−1,dAdt=βB∗S∗AA∗BA∗SB∗S∗A−1,dIdt=αA∗II∗AI∗IA∗−1,dBdt=ηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Now, we define a Lyapunov function as(36)L1=S−S∗+S∗lnSS∗+V−V∗+V∗lnVV∗+A−A∗+A∗lnAA∗+k1I−I∗+I∗lnII∗+k2B−B∗+B∗lnBB∗,where k1 and k2 are positive constants to be determined, and thus k1>0,k2>0. The time derivative of the Lyapunov function along with the solution of model (35) yields the following equation:(37)dL1dt=1−S∗SSΛS∗S∗S−1+θV∗S∗VS∗V∗S−1−βB∗BB∗−1+1−V∗VνS∗VV∗SV∗VS∗−1+1−A∗AβB∗S∗AA∗BA∗SB∗S∗A−1+k11−I∗IαA∗II∗AI∗IA∗−1+k21−B∗BηBB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
Simplifying the above expression and rearranging lead to(38)dL1dt=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
We now let(39)L1=ΛS−S∗1S−1S∗+θS−S∗VS−V∗S∗+νV−V∗SV−S∗V∗,(40)L2=−S−S∗βB∗BB∗−1+A−A∗βB∗S∗A∗BA∗SB∗S∗A−1+k1I−I∗αA∗I∗AI∗IA∗−1+k2B−B∗ηB∗A∗B∗AA∗B−1+I∗IB∗I∗B−1.
MakingΛ and ν the subject in equations (29) and (30), respectively, and substituting into (39) with series of simplification give(41)L1=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗V.
Further elaboration onL2 results in(42)L2=k2ηA∗1+AA∗−BB∗−AB∗BA∗+I∗1+II∗−BB∗−B∗IBI∗+αk1A∗AA∗−II∗−AI∗IA∗+1+βB∗S∗BB∗−AA∗+SS∗−BSA∗B∗S∗A.
To obtaink1,k2, we make the following assertion by equating the constant coefficients of the terms A/A∗,I/I∗ and B/B∗ in L2 to zero, thus(43)k2ηA∗+αk1A∗−βB∗S∗=0,(44)k2ηI∗−αk1A∗=0,(45)βB∗S∗−k2ηA∗+I∗=0.
Solving (44) and (45), we obtain k1=βB∗S∗I∗/αA∗A∗+I∗ and k2=βB∗S∗/ηA∗+I∗. Replacing k1,k2 into (42) and simplifying yield(46)L2=βB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Hence, adding equations (41) and (46), we have(47)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗2+AA∗−BB∗−B∗IBI∗−AI∗A∗I+βB∗S∗A∗A∗+I∗1+AA∗−BB∗−B∗ABA∗+βBS∗1+B∗SBS∗−B∗AA∗B−SA∗S∗A,(48)dL1dt=μS∗2−SS∗−S∗S+μV∗3−S∗S−SV∗VS∗−VV∗+θV∗2−VS∗SV∗−V∗SS∗VβB∗S∗I∗A∗+I∗D0+βB∗S∗A∗A∗+I∗D1+βBS∗D2,where(49)D0=2+AA∗−BB∗−B∗IBI∗−AI∗A∗I,D1=1+AA∗−BB∗−B∗ABA∗,D2=1+B∗SBS∗−B∗AA∗B−SA∗S∗A.
Thus, ifD0≤0,D1≤0,D2≤0, and from the fact that the arithmetic mean is greater than or equal to the geometric mean, then it follows that dL1/dt≤0, provided the coefficients μS∗,μV∗,θV∗,βB∗S∗I∗/A∗+I∗,βB∗S∗A∗/A∗+I∗,βBS∗ are nonnegative. With the claim that(50)wx=x−1−lnx,it leads to(51)D0=wAA∗−wBB∗−wB∗IBI∗−wAI∗A∗I,D1=wAA∗−wBB∗−wB∗ABA∗,D2=wB∗SBS∗−wB∗AA∗B−wSA∗S∗A.
To verify the above expression, the definition ofwx and properties of logarithms were used to simplify the expression on the right hand of Di,i=0,1,2. Then applying the proposition A1 of [28] gives(52)D0≤−wBB∗−wB∗IBI∗−wAI∗A∗I≤0,D1≤−wBB∗−wB∗ABA∗≤0,D2≤−wB∗AA∗B−wSA∗S∗A≤0.
It is evidential thatdL1S,V,A,I,B/dt=0 only holds when S=S∗,V=V∗,A=A∗,I=I∗,B=B∗ and that E∗ is the only equilibrium state of model (1) on its line. Therefore, by Lyapunov–LaSalle asymptotic stability theorem [27], the positive equilibrium state E∗ is globally asymptotically stable in the positive region ℝ+5 whenever ℛC>1 and ρ=0.


## 3. Sensitivity Analysis

Here, we apply the concept of sensitivity analysis to ascertain the relative importance of each generic parameter in the basic reproduction number and the disease-induced state variables,E∗. Secondly, with parameter values from published papers and some realistic assumptions, we obtain both analytic and numerical values of the various parameters in ℛC. The analytic expressions obtained can be used to shed some light on how to control the onset of the disease in variant localities, if and only if the dynamics follow model (1). Parameter values were obtained from [19, 20], which were then converted into days; this is to enable us to see the daily trajectories of the model. The estimated parameters were calculated based on the life expectancy of dairy cattle and the natural decay rate of the bacteria from dust/soil (which ranges from 120 to 150 days [9]), and finally, we used the gestation period of cattle to denote the rate of new inflow of susceptible cattle/recruited cattle. All parameters used are given in Tables 1 and 2, respectively.Table 1
Model parameter descriptions and values.
Parameter
Description
Range
Baseline value
Reference
β
Environmental transmission rate (indirect)
0.0030–0.0943
0.0943 day−1/cattle
[19]
Λ
Birth/recruitment rate
0.0034–0.0036
0.0036 day−1
Estimated
η
Proportion of bacteria shed through mucus/faeces filling the environment
0.0021–0.0400
0.0400 day−1/cattle
[19, 20]
ρ
Transition rate fromA to S
0.0774–0.1208
0.1000 day−1
[19, 20]
α
Transition rateA to I
0.0001–0.0117
0.0029 day−1
[19, 20]
θ
Loss of vaccine immunity
0.0027–2.0000
0.0027 day−1
[20]
μ
Natural death rate of cattle
0.0001–0.0002
0.0002 day−1
Estimated
ν
Vaccination rate
0.1000–0.9000 unitless
0.1000
Assumed
ε
Environmental hygiene
0.1000–0.9000 unitless
0.1000
Assumed
Δ
Natural decay of bacteria from the dust/soil
0.0067–0.0083
0.0083 day−1
Estimated
c
Culling/mortality rate
0.0008–0.0286
0.0286 day−1
[20]Table 2
The elasticity indices ofℛC=1.5193 and the relative change of A∗,I∗, and B∗ at the endemic equilibrium points, to the parameters in the model.
Parameter
Baseline value
Elasticity index,ℛC
A
∗
I
∗
B
∗
β
0.0943
+1.0000
+1.9255
+1.9255
+1.9255
Λ
0.0036
+1.0000
+2.9255
+2.9255
+2.9255
η
20.0400
+1.0000
+1.9255
+1.9255
+2.9255
ρ
20.0097
−0.7578
−1.4592
−1.4592
−1.4592
α
10.0029
−0.1351
−1.1956
−0.1956
−1.1041
θ
10.0027
+0.9048
+1.7422
+1.7422
+1.7422
μ
10.0002
−0.9492
−1.8923
−1.8992
−1.8929
ν
10.1000
−0.9718
−1.8713
−1.8713
−1.8713
ε
10.1000
−0.9234
−1.7779
−1.7779
−2.7013
δ
20.0083
−0.0766
−0.1476
−0.1476
−0.2242
c
0.0286
−0.0908
−0.1749
−1.1680
−0.2658


### 3.1. The Impact of Transition Rateρ on the Reproduction Number

Taking the partial derivative of equation (17) with respect to ρ gives the relative importance of the discontinuation rate ρ (the feedback term) and thus the rate at which asymptomatic cattle build antibiotics and reenter the susceptible class:(53)∂ℛ0∂ρ=−Λα+μβηα+μ+ρ2δμ2<0.Hence, the basic reproduction numberℛ0 without culling, vaccination, and environmental hygiene is a decreasing factor of the parameter ρ as shown in Subsection 2.4.1. This implies that the basic reproduction number decreases with increase in the value of ρ. Consequently, taking the limit of ℛ0 as ρ becomes large gives(54)limρ⟶∞ℛ0=0.This means that any attempt to increaseρ will result in reducing the number of secondary infections, thereby eradicating the disease in the absence of other control measures, but a low rate of ρ will produce a high rate of asymptomatic cattle which may require effective control measures as shown in Section 4. From this, we infer that the development of any daily antibody boosters for dairy cattle (livestock) will help curtail the reproduction losses caused by C. burnetii and other bacterial diseases.


### 3.2. Calculation of Sensitivity/Elasticity Indices

Fundamentally, the onset of disease transmission is directly associated with the basic/control reproduction numberℛC, and the prevalence of the disease is associated with the endemic equilibrium point E∗ [29]. The sensitivity indices for the state variables are calculated as follows:(55)ϒpiXi=∂Xi∂pi×piXi,where Xi represent the equilibrium points at the endemic state and pi are the corresponding parameters in Xi [30]. For simplicity, the numerical values of the endemic state variables are shown in Table 2 without analytic expression due to the complex nature of the analytic result. Similarly, the sensitivity/elasticity indices of ℛC to the parameters in the model are also defined as follows:(56)ϒpiℛC=∂ℛC∂pi×piℛC,where the meaning of ℛC and pi remains unchanged. Therefore(57)ϒβℛC=1,ϒΛℛC=1,ϒηℛC=1,ϒρℛC=−ρα+μ+ρ,ϒαℛC=−αc−ρα+c+μα+μ+ρ,ϒνℛC=−νμ+ν+θ,ϒθℛC=νθμ+θμ+ν+θ,ϒεℛC=−εδ+ε,ϒδℛC=−δδ+ε,ϒcℛC=−αcc+μα+c+μ,ϒμℛC=−μ1α+μ+ρ+1μ+ν+θ+1c+μ−1μ+θ+1μ−1α+c+μ.Remark 1.
ϒ
α
ℛ
C
<
0 whenever c>ρ, this implies that an increase in culling/isolation of symptomatic cattle from farms during an outbreak will reduce the impact of α, thereby reducing the number of secondary infections and bacterial load in the environment as shown in Table 2. This further suggests that in choosing between culling and antibiotic booster during an outbreak, it is prudent to focus on culling/isolation.
In Table2, we obtain numerical values which indicate the relative importance of parameters in ℛC and the state variables at the endemic stage when other parameters are kept constant (depicting the local sensitivity of parameters in model (1)). The result reveals that β,λ,η and θ influence the stability of the disease and are also closely related. Thus, a 1% increase in β,λ,η, and θ will have a corresponding percentage increase in the basic reproduction number. Among the control measures implemented, the elasticity index in Table 2 indicates that the most efficient method in the infected classes will be vaccination ν and the most efficient method of reducing the bacteria from the environment may be a good environmental hygiene ε. The total efficacy of these control measures will be shown in Section 4. From Table 2, it also shows that a 0.1 effort in keeping the environment clean/burying of placenta discharge after birth or abortion may result in a 2.70 reduction in bacterial load in the affected herds’ environment and a corresponding reduction in the number of asymptomatic and symptomatic by 1.78, respectively. Interestingly, we notice that a 1% change in natural death may reduce the reproduction number by 0.95, but natural death cannot be used as a control measure; it, therefore, presupposes that the disease could be established as long as cattle live longer in an environment that is induced by the bacteria when there are no control measures; hence, vaccinating cattle and good environmental hygiene in farms will help mitigate the intensity of the organism on farms. It is well noticing that the analysis done so far represents the local sensitivity analysis of parameters in the model, which shows a close relationship between some of the parameters; hence, to see the global effect of these parameters on the model, we adopt the Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) which is demonstrated in the next Subsection 3.2.1.


#### 3.2.1. Global Sensitivity Analysis

In carrying out the global sensitivity analysis on the control reproduction numberℛC, we use the analytic expression in equation (16) and the range of values in Table 1. The objective is to ascertain the most influential parameters in ℛC and also give some insightful ideas from the PRCC plot, by using a sample size of 1000.The scatter plot in Figure2(a) shows a positive relation between ℛC and the parameters β,Λ,η,θ,α and a negative relation between ℛC and the parameters ν and ρ. The positive relation indicates that a high rate of either of these parameters β,Λ,η,θ,α will produce a high transmission rate during an outbreak, while the negative relation for ν and ρ indicates that an increase in parameter value of ν and ρ will help reduce the severity of the disease-induced rate. Figure 2(b) gives sharply defined simulation results of the scatter plots in Figure 2(a) and that of the numerical signs in Table 2, but we notice a reverse sign for α, which is as a result of a decrease in the culling/isolation of symptomatic dairy cattle, as indicated in Remark 1; therefore, as shown in Figure 2(b) that, to have a positive impact of c in combination of other control measures so as to reduce α, c should be greater than 0.12, and thus it is represented by the red line as cvital. Figure 2(b) also indicates that the most influential parameters that promote the disease spread are β,η, and α (thus with low culling/isolation rate) and the most influential control parameter is ε (environmental hygiene). Therefore, any measures to reduce the relative relevance of β,η, and α (with an increase in culling/isolation rate) will reduce the spreading rate of the bacteria in dairy cattle herds. Figure 2(c) shows the saturation effect of vaccination and environmental hygiene, and thus having ℛC<1 and ν,ε>0.5 may lead to a Coxiella-free environment.


## 3.1. The Impact of Transition Rateρ on the Reproduction Number

Taking the partial derivative of equation (17) with respect to ρ gives the relative importance of the discontinuation rate ρ (the feedback term) and thus the rate at which asymptomatic cattle build antibiotics and reenter the susceptible class:(53)∂ℛ0∂ρ=−Λα+μβηα+μ+ρ2δμ2<0.Hence, the basic reproduction numberℛ0 without culling, vaccination, and environmental hygiene is a decreasing factor of the parameter ρ as shown in Subsection 2.4.1. This implies that the basic reproduction number decreases with increase in the value of ρ. Consequently, taking the limit of ℛ0 as ρ becomes large gives(54)limρ⟶∞ℛ0=0.This means that any attempt to increaseρ will result in reducing the number of secondary infections, thereby eradicating the disease in the absence of other control measures, but a low rate of ρ will produce a high rate of asymptomatic cattle which may require effective control measures as shown in Section 4. From this, we infer that the development of any daily antibody boosters for dairy cattle (livestock) will help curtail the reproduction losses caused by C. burnetii and other bacterial diseases.


## 3.2. Calculation of Sensitivity/Elasticity Indices

Fundamentally, the onset of disease transmission is directly associated with the basic/control reproduction numberℛC, and the prevalence of the disease is associated with the endemic equilibrium point E∗ [29]. The sensitivity indices for the state variables are calculated as follows:(55)ϒpiXi=∂Xi∂pi×piXi,where Xi represent the equilibrium points at the endemic state and pi are the corresponding parameters in Xi [30]. For simplicity, the numerical values of the endemic state variables are shown in Table 2 without analytic expression due to the complex nature of the analytic result. Similarly, the sensitivity/elasticity indices of ℛC to the parameters in the model are also defined as follows:(56)ϒpiℛC=∂ℛC∂pi×piℛC,where the meaning of ℛC and pi remains unchanged. Therefore(57)ϒβℛC=1,ϒΛℛC=1,ϒηℛC=1,ϒρℛC=−ρα+μ+ρ,ϒαℛC=−αc−ρα+c+μα+μ+ρ,ϒνℛC=−νμ+ν+θ,ϒθℛC=νθμ+θμ+ν+θ,ϒεℛC=−εδ+ε,ϒδℛC=−δδ+ε,ϒcℛC=−αcc+μα+c+μ,ϒμℛC=−μ1α+μ+ρ+1μ+ν+θ+1c+μ−1μ+θ+1μ−1α+c+μ.Remark 1.
ϒ
α
ℛ
C
<
0 whenever c>ρ, this implies that an increase in culling/isolation of symptomatic cattle from farms during an outbreak will reduce the impact of α, thereby reducing the number of secondary infections and bacterial load in the environment as shown in Table 2. This further suggests that in choosing between culling and antibiotic booster during an outbreak, it is prudent to focus on culling/isolation.
In Table2, we obtain numerical values which indicate the relative importance of parameters in ℛC and the state variables at the endemic stage when other parameters are kept constant (depicting the local sensitivity of parameters in model (1)). The result reveals that β,λ,η and θ influence the stability of the disease and are also closely related. Thus, a 1% increase in β,λ,η, and θ will have a corresponding percentage increase in the basic reproduction number. Among the control measures implemented, the elasticity index in Table 2 indicates that the most efficient method in the infected classes will be vaccination ν and the most efficient method of reducing the bacteria from the environment may be a good environmental hygiene ε. The total efficacy of these control measures will be shown in Section 4. From Table 2, it also shows that a 0.1 effort in keeping the environment clean/burying of placenta discharge after birth or abortion may result in a 2.70 reduction in bacterial load in the affected herds’ environment and a corresponding reduction in the number of asymptomatic and symptomatic by 1.78, respectively. Interestingly, we notice that a 1% change in natural death may reduce the reproduction number by 0.95, but natural death cannot be used as a control measure; it, therefore, presupposes that the disease could be established as long as cattle live longer in an environment that is induced by the bacteria when there are no control measures; hence, vaccinating cattle and good environmental hygiene in farms will help mitigate the intensity of the organism on farms. It is well noticing that the analysis done so far represents the local sensitivity analysis of parameters in the model, which shows a close relationship between some of the parameters; hence, to see the global effect of these parameters on the model, we adopt the Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) which is demonstrated in the next Subsection 3.2.1.


### 3.2.1. Global Sensitivity Analysis

In carrying out the global sensitivity analysis on the control reproduction numberℛC, we use the analytic expression in equation (16) and the range of values in Table 1. The objective is to ascertain the most influential parameters in ℛC and also give some insightful ideas from the PRCC plot, by using a sample size of 1000.The scatter plot in Figure2(a) shows a positive relation between ℛC and the parameters β,Λ,η,θ,α and a negative relation between ℛC and the parameters ν and ρ. The positive relation indicates that a high rate of either of these parameters β,Λ,η,θ,α will produce a high transmission rate during an outbreak, while the negative relation for ν and ρ indicates that an increase in parameter value of ν and ρ will help reduce the severity of the disease-induced rate. Figure 2(b) gives sharply defined simulation results of the scatter plots in Figure 2(a) and that of the numerical signs in Table 2, but we notice a reverse sign for α, which is as a result of a decrease in the culling/isolation of symptomatic dairy cattle, as indicated in Remark 1; therefore, as shown in Figure 2(b) that, to have a positive impact of c in combination of other control measures so as to reduce α, c should be greater than 0.12, and thus it is represented by the red line as cvital. Figure 2(b) also indicates that the most influential parameters that promote the disease spread are β,η, and α (thus with low culling/isolation rate) and the most influential control parameter is ε (environmental hygiene). Therefore, any measures to reduce the relative relevance of β,η, and α (with an increase in culling/isolation rate) will reduce the spreading rate of the bacteria in dairy cattle herds. Figure 2(c) shows the saturation effect of vaccination and environmental hygiene, and thus having ℛC<1 and ν,ε>0.5 may lead to a Coxiella-free environment.


## 3.2.1. Global Sensitivity Analysis

In carrying out the global sensitivity analysis on the control reproduction numberℛC, we use the analytic expression in equation (16) and the range of values in Table 1. The objective is to ascertain the most influential parameters in ℛC and also give some insightful ideas from the PRCC plot, by using a sample size of 1000.The scatter plot in Figure2(a) shows a positive relation between ℛC and the parameters β,Λ,η,θ,α and a negative relation between ℛC and the parameters ν and ρ. The positive relation indicates that a high rate of either of these parameters β,Λ,η,θ,α will produce a high transmission rate during an outbreak, while the negative relation for ν and ρ indicates that an increase in parameter value of ν and ρ will help reduce the severity of the disease-induced rate. Figure 2(b) gives sharply defined simulation results of the scatter plots in Figure 2(a) and that of the numerical signs in Table 2, but we notice a reverse sign for α, which is as a result of a decrease in the culling/isolation of symptomatic dairy cattle, as indicated in Remark 1; therefore, as shown in Figure 2(b) that, to have a positive impact of c in combination of other control measures so as to reduce α, c should be greater than 0.12, and thus it is represented by the red line as cvital. Figure 2(b) also indicates that the most influential parameters that promote the disease spread are β,η, and α (thus with low culling/isolation rate) and the most influential control parameter is ε (environmental hygiene). Therefore, any measures to reduce the relative relevance of β,η, and α (with an increase in culling/isolation rate) will reduce the spreading rate of the bacteria in dairy cattle herds. Figure 2(c) shows the saturation effect of vaccination and environmental hygiene, and thus having ℛC<1 and ν,ε>0.5 may lead to a Coxiella-free environment.


## 4. Optimal Control Problem

Whenever the number of unknown parameters in a model is larger than the number of state equations, then the model in question usually allows for more than one solution. The optimization problem for such situations in continuous models is known as an optimal control problem. Here, we investigate the impact of time-variant controls: culling, vaccination, and environmental hygiene (cleaning/burying of placenta discharge after birth or abortion). This is to obtain the best control measure in preventing an outbreak of the disease in a situation where the cause of new infections is solely environmental. Thus, we seek to find optimal trajectories that show the relative importance of applying vaccination, culling, and environmental hygiene in controlling Q fever in an endemic setting. To see the impact of these controls, we modify model (1), thus setting ρ=0 (no natural resistance). The objective function is to reduce the numbers of asymptomatic cattle, symptomatic cattle, and bacterial load through vaccination ν, culling c, and environmental hygiene ε, respectively, and the cost associated with the implementation of the controls. Thus u1t∈Uad is the percentage of susceptible and asymptomatic cattle who are vaccinated per unit time, u2t∈Uad is the percentage of symptomatic cattle who are culled per unit time, and u3∈Uad is admissible measure for keeping the environment clean. Here Uad=uiuitare all  Lebesgue measurable and0≤uit≤umax,t∈0,tf, where i=1,2,3 and tf is the final time for running the control measures. Now the objective function that minimizes our desirable control problem is given as(58)Juit=minui,i=1,2,3∫0tfW1At+W2It+W3Bt+C12u12t+C22u22t+C32u32tdt,subject to(59)dStdt=Λ−1−u3tβBtSt−μ+u1tSt+θVt,dVtdt=u1tSt−θ+μVt,dAtdt=1−u3tβBtSt−μ+αAt,dItdt=αAt−u2t+μIt,dBtdt=ηAt+It−u3t+δBt,with initial conditions(60)S0=S0>0,V0=V0>0,A0=A0≥0,I0=I0≥0,B0=B0≥0,where 1−u3tβ is the reduction in transmission due to environmental hygiene, Wi,i=1,2,3 are small positive constant to keep balance in the size of St, At, It, and Bt, respectively [31]. The square term of the controls reflects the nonlinearity in the cost of controls, and the half-term minimizes the effect of applying the controls [32]. The positive constants C1,C2,C3 which are associated with uit are weight values such that 0<C1,C2<N, and 0<C3<B. The relative weights given to the positive constants associated with the control terms indicate greater or lower importance placed on minimizing the cost of a control measure [33]. Now to obtain the optimal solutions, we defined a Lagrangian L for the control problem, given as follows:(61)LS,A,I,B,u1,u2,u3=W1At+W2It+W3Bt+C12u12t+C22u22t+C32u32t.We now seek the minimal value of the Lagrangian function. Hence, we obtain a HamiltonianH for the control problem using Pontryagin’s maximum principle:(62)HS,A,I,B,ui,λj,t=L+λ1tdStdt+λ2tdVtdt+λ3tdAtdt+λ4tdItdt+λ5tdBtdt,where λj,j=1,2,3,4,5 are the adjoint variables associated with the states S,V,A,I, and B. Finally, taking the partial derivatives of H with respect to the state variables and the controls, we obtain the following relation for the adjoint and the optimal controls:(63)dλ1dt=−∂H∂S,dλ2dt=−∂H∂V,dλ3dt=−∂H∂A,dλ4dt=−∂H∂I,dλ5dt=−∂H∂B,u1∗=∂H∂u1=0,u2∗=∂H∂u2=0,u3∗=∂H∂u3=0,with transversality conditions (or boundary conditions) [31], λjtf=0,j=1,2,3,4,5. Therefore, the adjoint system is given as(64)dλ1dt=λ1μ+u1−Bβu3−1−λ2u1+Bβλ3u3−1,dλ2dt=λ2μ+θ−λ1θ,dλ3dt=λ3α+μ−αλ4−ηλ5−W1,dλ4dt=λ4μ+u2−W2−ηλ5,dλ5dt=λ5δ+u3−W3−Sβλ1u3−1+Sβλ3u3−1.By the optimality conditions, we have(65)∂H∂u1≔C1u1−Sλ1+Sλ2=0,∂H∂u2≔C2u2−Iλ4=0,∂H∂u3≔C3u3−Bλ5+BSβλ1−BSβλ3=0,⟹u1∗=Sλ1−Sλ2C1,u2∗=Iλ4C2,u3∗=Bλ5−BSβλ1+BSβλ3C3.Using the optimal solution in the control space [31], we have(66)u1∗=0,ifS∗λ1−S∗λ2C1≤0,S∗λ1−S∗λ2C1,if0<S∗λ1−S∗λ2C1<0.9,0.9,ifS∗λ1−S∗λ2C1≥0.9,u2∗=0,ifλ4I∗C2≤0,λ4I∗C2,if0<λ4I∗C2<1,1,ifλ4I∗C2≥1,u3∗=0,ifB∗λ5−B∗S∗βλ1+B∗S∗βλ3C3≤0,B∗λ5−B∗S∗βλ1+B∗S∗βλ3C3,if0<B∗λ5−B∗S∗βλ1+B∗S∗βλ3C3<0.8,0.8,ifB∗λ5−B∗S∗βλ1+B∗S∗βλ3C3≥0.8.Therefore, the compact notation (characterization) of the admissible control is then written as(67)u1∗=minmax0,λ1S∗−B∗S∗β−S∗λ2+B∗S∗βλ3C1,0.9,u2∗=minmax0,λ4I∗C2,1,u3∗=minmax0,λ5B∗−λ1βB∗S∗+λ3βB∗S∗C3,0.8.The control terms in the bound of 0.9 and 0.8 satisfy the inability to vaccinate all susceptible cattle at a time and a 0.2 failure in eliminating all bacteria in a particular farm at a given time, note that the control bounds are assumed for a realistic illustrative purposes and may not be practically practised, but the variation of these bounds will not narrow the essence of the studies. Substituting (67) into (59), (62), and (64), respectively, we obtain an optimality system with the corresponding Hamiltonian system H∗ which is then solved numerically to obtain optimal trajectories for our dynamical system.


### 4.1. Optimal Control Trajectories

In order to obtain the optimal trajectories for model (59), we employ the forward-backward sweep method, as described in [34]. We begin with an initial guess for the controls. The state equation (59) is solved forward in time while the adjoint function (64) is solved backward in time. The controls are then updated for a period of time until convergence is achieved. For better numerical output, we firstly made use of the following initial state values: S0=100,V0=2,A0=3,I0=1,B0=0.0025 and the stated parameter values in Table 1. The initial bacterial load is assumed to be small, the reason being that per the numerical simulations, the bacterial population grows more and more rapidly as depicted in Figures 3(f) and 4(f), respectively, by the red line, hence, the need to start with small initial values for the bacterial load so as to get a suitable numerical output for that of the initial conditions chosen for the animal population. The corresponding weight in the objective function (58) is taken to be W1=2,W2=4,W3=2,C1=1,C2=1,C3=0.2. The figures below depict the relative importance of finding optimal strategy/strategies in controlling the disease spread in cattle (animal population). Table 3 shows the total cost associated with the objective function (58) after a simulation time of 100 days. Secondly, we look at the control trajectories and the accompanying cost by changing the associated weights of control, cost of control implantation, and initial state values and thus W1=1,W2=1,W3=2,C1=1,C2=1,C3=0.5 and S0=1000,V0=10,A0=30,I0=50, B0=0.25 as depicted in Figure 4 and Table 3, item 3.Figure 3
Optimal trajectory results with and without controls. All initial state values are assumed for illustrative purposes. (a) Control profile. (b) Susceptible cattle. (c) Vaccinated cattle. (d) Asymptomatic cattle. (e) Symptomatic cattle. (f) Bacterial load.
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
Optimal trajectory results with and without controls, with a change in initial state value data. All initial state values are assumed for illustrative purposes. (a) Control profile. (b) Susceptible cattle. (c) Vaccinated cattle. (d) Asymptomatic cattle. (e) Symptomatic cattle. (f) Bacterial load.
(a)
(b)
(c)
(d)
(e)
(f)Table 3
The total costJ in regard to different control strategies associated with model (1) with ρ=0.
Item
Strategies
Associated costJ
1
Without optimal control,u1=0,u2=0,u3=0
28,977.63
2
With optimal control,u1≠0,u2≠0,u3≠0
623.04
3
With optimal control,u1≠0,u2≠0,u3≠0 with S0=1000
6,537.69
4
Culling only,u1=0,u2≠0,u3=0
23,057.75
5
Without environmental hygiene,u1≠0,u2≠0,u3=0
1,616.53
6
Without culling,u1≠0,u2=0,u3≠0
1,191.37
7
Without vaccination,u1=0,u2≠0,u3≠0
10,722.43Figure3(a) shows the control profile for the various bounds employed in model (59), which indicates that a bound of 0.8 for environmental hygiene is maintained for the entire simulation time, but a drop in the control bound for both vaccination and culling. Thus, the bound for culling and vaccination switches to 0.32 after 10 days. This indicates during the initial onset of an outbreak, effort should be made to cull all symptomatic cattle or isolate them from other susceptible (healthy) cattle for the first 10 days, after which less culling approach can be adopted, and also for first 10 days, at least 90% of cattle should receive vaccination during the initial stage of the disease outbreak. Figures 3(b) and 3(c) show that applying vaccination, culling and environmental hygiene simultaneously increases the rate of healthy cattle, thereby reducing the risk of infection in the susceptible cattle population and also reducing the number of asymptomatic cattle, symptomatic cattle, and bacterial load in the environment as shown in Figures 3(d)–3(f), respectively. We also notice that the gradual use of these controls in the long run is more cost-effective as shown in Table 3, item 2.In Figure4(a), we changed the initial population size, which indicates that an increase in population size requires that the initial control bound for vaccination be kept for 80 days, and thus an increase in population size of dairy cattle herds requires high control effort and high resources as depicted in item 3 of Table 3. The control trajectory for the vaccination bound and Figure 4(c) further indicate that for optimal eradication or reduction of the disease in a population of 1000 cattle (livestock), at least 80% of the animals should be vaccinated within the first 20–80 days. We also noticed that the use of vaccination, culling, and environmental hygiene has similar control efficacy in the disease compartments irrespective of the population size as shown in Figures 4(b)–4(f), respectively.In Figure5, we considered variations that may arise in the control measures, thus sticking to one or more controls. In Figure 5(a), we assumed that u2=0, that is a situation where culling is not applied, and it reveals that symptomatic cattle will always be in the population but stays below the trajectory when all controls are zero (see Figure 5(b)). Therefore, we notice that in a situation where a farm owner chooses to ignore culling and decide to stick to only vaccination and environmental hygiene during an outbreak of the disease may not lead to the eradication of the disease; hence, the combination of vaccination and environmental hygiene is not a good control measure in reducing the number of symptomatic cattle (see Figure 6(d)) for the implementation of culling only as a control measure on the symptomatic compartment. Figure 5(c) shows that in the absence of vaccination, the control bound for culling should remain at 1, which means that during the cause of culling, the aim should be focused on destroying all infected farms or animals which is sometimes practically hard to achieve and also comes with high control cost as shown in Table 3, item 7. Therefore, it clearly shows in Figures 5(d)–5(f) that effective vaccination plays a major role in reducing the disease-induced rate in dairy cattle herds depending only on culling. It can also be seen that without vaccination, the trajectories of Figures 5(d)–5(f) of the optimal control stay very close to the trajectories without optimal control.Figure 5
Optimal trajectory results marked by diversity in control procedures. All initial state values are assumed for illustrative purposes. (a) Without culling control profile. (b) Without culling in the objective function. (c) Without vaccinated control profile. (d) Without vaccination in the objective function. (e) Asymptomatic cattle (without vaccination in the objective function). (f) Bacterial load (without vaccination in the objective function).
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Optimal trajectory results marked by diversity in control procedures. All initial state values are assumed for illustrative purposes. (a) Control values (no environmental hygiene). (b) Bacterial load (no environmental hygiene). (c) Asymptomatic cattle (no environmental hygiene). (d) Symptomatic cattle (culling only).
(a)
(b)
(c)
(d)Figure6(a) shows the control profile in the absence of environmental hygiene; the figure depicts that with only vaccination and culling as control measures, the bound of vaccination is kept at 0.9 for 78 days before switching to 0.8, and the bound of culling switches to 0.45 as compared to Figure 3(a). Hence, Figure 6(b) depicts the noneffectiveness of vaccination and culling on the bacterial load in the environment until day 47, but due to vaccination, the number of asymptomatic cattle reduces as shown in Figure 6(c). Figure 6(d) shows the effect of culling on the number of symptomatic cattle, which indicates that early detection of symptomatic cattle on the farm will help in mitigating the disease within the first 5 days.


## 4.1. Optimal Control Trajectories

In order to obtain the optimal trajectories for model (59), we employ the forward-backward sweep method, as described in [34]. We begin with an initial guess for the controls. The state equation (59) is solved forward in time while the adjoint function (64) is solved backward in time. The controls are then updated for a period of time until convergence is achieved. For better numerical output, we firstly made use of the following initial state values: S0=100,V0=2,A0=3,I0=1,B0=0.0025 and the stated parameter values in Table 1. The initial bacterial load is assumed to be small, the reason being that per the numerical simulations, the bacterial population grows more and more rapidly as depicted in Figures 3(f) and 4(f), respectively, by the red line, hence, the need to start with small initial values for the bacterial load so as to get a suitable numerical output for that of the initial conditions chosen for the animal population. The corresponding weight in the objective function (58) is taken to be W1=2,W2=4,W3=2,C1=1,C2=1,C3=0.2. The figures below depict the relative importance of finding optimal strategy/strategies in controlling the disease spread in cattle (animal population). Table 3 shows the total cost associated with the objective function (58) after a simulation time of 100 days. Secondly, we look at the control trajectories and the accompanying cost by changing the associated weights of control, cost of control implantation, and initial state values and thus W1=1,W2=1,W3=2,C1=1,C2=1,C3=0.5 and S0=1000,V0=10,A0=30,I0=50, B0=0.25 as depicted in Figure 4 and Table 3, item 3.Figure 3
Optimal trajectory results with and without controls. All initial state values are assumed for illustrative purposes. (a) Control profile. (b) Susceptible cattle. (c) Vaccinated cattle. (d) Asymptomatic cattle. (e) Symptomatic cattle. (f) Bacterial load.
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
Optimal trajectory results with and without controls, with a change in initial state value data. All initial state values are assumed for illustrative purposes. (a) Control profile. (b) Susceptible cattle. (c) Vaccinated cattle. (d) Asymptomatic cattle. (e) Symptomatic cattle. (f) Bacterial load.
(a)
(b)
(c)
(d)
(e)
(f)Table 3
The total costJ in regard to different control strategies associated with model (1) with ρ=0.
Item
Strategies
Associated costJ
1
Without optimal control,u1=0,u2=0,u3=0
28,977.63
2
With optimal control,u1≠0,u2≠0,u3≠0
623.04
3
With optimal control,u1≠0,u2≠0,u3≠0 with S0=1000
6,537.69
4
Culling only,u1=0,u2≠0,u3=0
23,057.75
5
Without environmental hygiene,u1≠0,u2≠0,u3=0
1,616.53
6
Without culling,u1≠0,u2=0,u3≠0
1,191.37
7
Without vaccination,u1=0,u2≠0,u3≠0
10,722.43Figure3(a) shows the control profile for the various bounds employed in model (59), which indicates that a bound of 0.8 for environmental hygiene is maintained for the entire simulation time, but a drop in the control bound for both vaccination and culling. Thus, the bound for culling and vaccination switches to 0.32 after 10 days. This indicates during the initial onset of an outbreak, effort should be made to cull all symptomatic cattle or isolate them from other susceptible (healthy) cattle for the first 10 days, after which less culling approach can be adopted, and also for first 10 days, at least 90% of cattle should receive vaccination during the initial stage of the disease outbreak. Figures 3(b) and 3(c) show that applying vaccination, culling and environmental hygiene simultaneously increases the rate of healthy cattle, thereby reducing the risk of infection in the susceptible cattle population and also reducing the number of asymptomatic cattle, symptomatic cattle, and bacterial load in the environment as shown in Figures 3(d)–3(f), respectively. We also notice that the gradual use of these controls in the long run is more cost-effective as shown in Table 3, item 2.In Figure4(a), we changed the initial population size, which indicates that an increase in population size requires that the initial control bound for vaccination be kept for 80 days, and thus an increase in population size of dairy cattle herds requires high control effort and high resources as depicted in item 3 of Table 3. The control trajectory for the vaccination bound and Figure 4(c) further indicate that for optimal eradication or reduction of the disease in a population of 1000 cattle (livestock), at least 80% of the animals should be vaccinated within the first 20–80 days. We also noticed that the use of vaccination, culling, and environmental hygiene has similar control efficacy in the disease compartments irrespective of the population size as shown in Figures 4(b)–4(f), respectively.In Figure5, we considered variations that may arise in the control measures, thus sticking to one or more controls. In Figure 5(a), we assumed that u2=0, that is a situation where culling is not applied, and it reveals that symptomatic cattle will always be in the population but stays below the trajectory when all controls are zero (see Figure 5(b)). Therefore, we notice that in a situation where a farm owner chooses to ignore culling and decide to stick to only vaccination and environmental hygiene during an outbreak of the disease may not lead to the eradication of the disease; hence, the combination of vaccination and environmental hygiene is not a good control measure in reducing the number of symptomatic cattle (see Figure 6(d)) for the implementation of culling only as a control measure on the symptomatic compartment. Figure 5(c) shows that in the absence of vaccination, the control bound for culling should remain at 1, which means that during the cause of culling, the aim should be focused on destroying all infected farms or animals which is sometimes practically hard to achieve and also comes with high control cost as shown in Table 3, item 7. Therefore, it clearly shows in Figures 5(d)–5(f) that effective vaccination plays a major role in reducing the disease-induced rate in dairy cattle herds depending only on culling. It can also be seen that without vaccination, the trajectories of Figures 5(d)–5(f) of the optimal control stay very close to the trajectories without optimal control.Figure 5
Optimal trajectory results marked by diversity in control procedures. All initial state values are assumed for illustrative purposes. (a) Without culling control profile. (b) Without culling in the objective function. (c) Without vaccinated control profile. (d) Without vaccination in the objective function. (e) Asymptomatic cattle (without vaccination in the objective function). (f) Bacterial load (without vaccination in the objective function).
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Optimal trajectory results marked by diversity in control procedures. All initial state values are assumed for illustrative purposes. (a) Control values (no environmental hygiene). (b) Bacterial load (no environmental hygiene). (c) Asymptomatic cattle (no environmental hygiene). (d) Symptomatic cattle (culling only).
(a)
(b)
(c)
(d)Figure6(a) shows the control profile in the absence of environmental hygiene; the figure depicts that with only vaccination and culling as control measures, the bound of vaccination is kept at 0.9 for 78 days before switching to 0.8, and the bound of culling switches to 0.45 as compared to Figure 3(a). Hence, Figure 6(b) depicts the noneffectiveness of vaccination and culling on the bacterial load in the environment until day 47, but due to vaccination, the number of asymptomatic cattle reduces as shown in Figure 6(c). Figure 6(d) shows the effect of culling on the number of symptomatic cattle, which indicates that early detection of symptomatic cattle on the farm will help in mitigating the disease within the first 5 days.


## 5. Concluding Remarks

We presented a theoretical analysis of Q fever within dairy cattle herds, which may lead to new research on the disease dynamics, such as the effect of direct transmission, seasonal transmission, the dynamical role of ticks in maintaining the circle of Q fever in livestock, the role of media in averting the disease in humans, and the dynamical spread of the disease between livestock and humans. Here we analytically obtained the control states and numerically simulated it for different outcomes. We notice that culling only as a control measure produces a higher control cost, which also does not completely eradicate the number of symptomatic cattle due to the asymptomatic nature of the disease. We obtained a cost profile for the various control measures as shown in Table3 which suggests that the most effective method in getting a lower control cost, in the long run, is the continuous application of vaccination and good environmental hygiene together with culling or isolation of symptomatic cattle. The numerical simulation also reveals that irrespective of the control employed, the disease will always be present in the animal population, but vaccinating and cleaning/burying of placenta discharge after birth or abortion will help to mitigate the disease in the animal population. In totality, our work has revealed some global dynamics of Q fever in dairy cattle herds and the impact of sensitivity analysis and optimal control strategies in managing Q fever disease in dairy cattle herds (livestock). We also analytically obtained the local and global stability analysis of the disease by using matrix-theoretic method and a Lyapunov functional. The model studied here did not assume a complex model, but the results can provide some ideas for further studies and the prevention of the disease.



---

*Source: 6820608-2020-01-31.xml*