Dynamics and Stability Analysis of a Brucellosis Model with Two Discrete Delays

We present a mathematical model for brucellosis transmission that incorporates two discrete delays and culling of infected animals displaying signs of brucellosis infection. The first delay represents the incubation period while the second account for the time needed to detect and cull infectious animals. Feasibility and stability of the model steady states have been determined analytically and numerically. Further, the occurrence of Hopf bifurcation has been established. Overall the findings from the study, both analytical and numerical, suggest that the two delays can destabilize the system and periodic solutions can arise through Hopf bifurcation.


Introduction
Brucellosis is one of the neglected zoonotic diseases that remains a major public health problem world over, especially in Middle Eastern countries, southern Europe and North Africa, countries in South and Central Asia, sub-Saharan Africa, Mexico, the Caribbean, and countries in South and Central America [1], with an annual occurrence of more than 500 000 cases [2].
In animals, brucellosis is usually transmitted through direct contact between a susceptible and an infectious animal or indirectly, i.e., when a susceptible animal ingest contaminated forages or the excrement containing large quantities of bacteria, generally discharged by infected animals [3].In humans, however, the majority of the infections result from direct or indirect exposure to infected animals or consumption of raw animal products such as unpasteurized milk or cheese [4].Since human-to-human transmission of the disease is extremely rare [5], the ultimate management of human brucellosis can be achieved through effective control of brucellosis in livestock.Some researchers postulate that eradication of brucellosis in animals can be attained by combining vaccination with test-and-slaughter programs [1].
Mathematical models have proved to be essential guiding tools for epidemiologists, biologists, and policy makers.Models can provide solutions to phenomena which are difficult to measure practically.Recently, a number of mathematical models have been proposed to study the spread and control of brucellosis (see, for example, [3,[6][7][8][9][10][11][12][13][14][15] and references therein).A limitation of these previous studies, however, is the noninclusion of the time taken before an infectious animal is detected and culled, despite the fact that in many countries where the disease is endemic lack of financial and human resources often results in delay of detection and culling of infectious animals.The size of this delay may play an important role in minimizing the spread of the disease in the community.
It is therefore essential to gain a better and more comprehensive understanding of the effects of time delay on brucellosis transmission and control.Prior studies have shown that epidemic models with time delay often exhibit periodic solutions and as a consequence understanding the nature of these periodic outbreaks plays a crucial role in designing policies that can successfully control the disease.In fact, a recent analysis of brucellosis dataset in countries where the disease is endemic has shown that the disease incidences exhibit a strong periodic behavior with mortality and morbidity of the disease concentrated in a few months each year [16,17].Understanding the impact of such seasonal variations is crucial in managing the spread of the disease in the community.
Our main goal in this study is to explore the dynamics and stability analysis of a brucellosis model with two discrete delays.Hence we formulated a mathematical model that incorporates two discrete delays.The first delay represents the incubation period while the second delay accounts for the time taken to detect and cull infectious animals.In addition, we subdivide the total animal population into classes of susceptible, vaccinated, infectious undetected, and infectious detected animals.In certain situations immediate slaughter of detected animals may not be feasible and more often these animals are isolated from the rest.However, due to lack of financial and human resources, in addition to lack of knowledge and attitude of farmers, isolation of detected animals has not been a successful practice in most developing nations where animal infections are rampant.Thus in our modelling process we assume that a proportion of detected animals that are not immediately culled are also responsible for disease transmission.Utilizing both analytical and numerical results we have demonstrated that the two delays can destabilize the system and lead to Hopf bifurcation.
The article is organised as follows.The model description is given in Section 2. Analytical and numerical results are given in Sections 3 and 4, respectively.We end with Section 5 of conclusions.

Mathematical Model
We subdivide the total animal population () into compartments of susceptible (), vaccinated (), undetected infectious animals  1 (), and infectious detected and unculled  2 ().Although brucellosis can be transmitted indirectly (environmental transmission), prior studies [6,8] suggest that indirect transmission plays a relatively small role on the spread of brucellosis, and as such we have ignored this aspect in our study.Brucellosis dynamics in this study are governed by the following autonomous system: where  is the recruitment rate through birth,  is the natural death rate,  is the disease direct transmission rate,  is fraction of detected animals that have been culled,  is the vaccination rate,  is the vaccination waning rate,  1 represents the latency delay,  2 is the delay in detection,  is the modification factor,  is the rate at which animals are detected and quarantined,  is the culling rate of detected animals, and  is disease-induced death.

Stability of the Disease-Free Equilibrium.
In this section, we study the local and global stability of the disease-free equilibrium.
Theorem 1.The disease-free equilibrium E 0 of model ( 1) is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.
Proof.To study the local stability of the disease-free equilibrium E 0 , we linearized system (1) about this point and obtained the characteristic equation, given by the following determinant: where  is the eigenvalue.From ( 5) the characteristic equation is Clearly, − and −( +  + ) are eigenvalues and the other two can be determined from the following equation: Through direct calculation one can easily verify that (,  1 ,  2 ) is an analytic function and it follows that Now we proceed to investigate the distribution of the solutions of (7) in the following cases: (a) If R 0 < 1, then (0,  1 ,  2 ) > 0. Since the derivative    (,  1 ,  2 ) > 0 for  ≥ 0,  1 > 0 and  2 > 0, (7) has no zero root and positive real roots for all positive  1 and  2 .Now we assume that the solution of (7) does not have any purely imaginary roots  = , ( > 0) for some  1 > 0,  2 = 0. Then by computation,  must be positive real root of If R 0 < 1, (10) has no positive roots.Hence (7) does not have any purely imaginary roots.We can easily see that the roots of (, 0, 0) = 0 all have negative real parts when R 0 < 1.By the implicit function theorem and the continuity of (,  1 ,  2 ), we know that all roots of (7) have negative real parts for positive  1 and  2 = 0 which implies that E 0 is stable.
Since the derivative   (,  1 ,  2 ) > 0 for  ≥ 0,  1 > 0 and  2 > 0, (7) has a simple zero and no positive root for all positive  1 and  2 .By the same argument in case (), we can obtain that all roots of (7) have negative real parts for positive  1 and  2 = 0 except a zero root.Thus E 0 is a degenerate equilibrium of codimension and is stable except in one direction.
Theorem 2. The disease-free equilibrium of model ( 1) is globally asymptotically stable when R 0 ≤ 1 and unstable when R 0 > 1.

Disease Persistence.
System (1) is said to be uniformly persistent if there exists a constant  0 > 0 such that any solution ((), (),  1 (),  2 ()) of (1) satisfies Now we give a result on the uniform persistence of system (1).To proceed we introduce the following notation and terminology.Denote by (),  ≥ 0 the family of solution operators corresponding to (1).-limit set () of  consists of  ∈  such that there exists a sequence   → ∞ as  → ∞ with (  ) →  as  → ∞.
This completes the proof of Theorem 3. 1) admits a unique endemic equilibrium.

Existence of the Endemic Equilibrium
Proof.The endemic equilibrium E * = ( * ,  * ,  * 1 ,  * 2 ) of model ( 1) is determined by the following: From the last equation in (24) we have The first two equations in (24) give For  1 ̸ = 0, substituting (25) into the third equation in (24) gives Substituting ( 26) into (27) yields Direct calculations show where then function ( 1 ) is monotonic decreasing for  1 > 0, and then we can define the following function: Therefore, by monotonicity of a function ( 1 ) there exists a unique endemic equilibrium E * = ( * ,  * ,  * 1 ,  * 2 ) 3.6.Stability of the Endemic Equilibrium.In this section, we will investigate the local and global stability of the endemic equilibrium point.
Theorem 5.The endemic equilibrium E * of the system ( 1) is locally asymptotically stable if R 0 > 1 and conditions (36) are satisfied.
Proof.The characteristics equation of system (1) on the infected equilibrium E * is given by the following determinant: with After some algebraic manipulations one can show that the characteristic equation has the following form: with By the Routh-Hurwitz criterion, all roots of the characteristics equation (34) have negative real parts and the endemic equilibrium E * of system (1) is locally asymptotically stable if Now, we wish to explore if there is a possibility of having complex roots with positive real part for (a)  1 > 0,  2 = 0 and (b)  1 = 0,  2 > 0. We now proceed to explore the above cases as follows.
(a) If  1 > 0,  2 = 0, then the characteristics equation (34) becomes with Now we need to show that all roots of (37) have negative real parts for all  1 ∈ (0,  * ).To do so, we show that (37) does not have any purely imaginary roots for all  1 ∈ (0,  * ).We assume that  =  with  > 0 is a root of (37).Then  must satisfy the following system: Now, we square both sides of each equation above and add the resulting equations; thus  must be a positive root of where Let  =  2 ; then (40) becomes One can observe that if   ≥ 0, ( = 1, 2, 3, 4), then (42) has no positive roots.Therefore (37) does not have any purely imaginary roots for all  1 ∈ (0,  * ) so that all roots of the characteristic equation (37) have negative real parts and the endemic equilibrium E * of ( 1) is stable for all  1 ∈ (0,  * ).
(b) If  2 > 0, Using the same discussion as in the above case then (43) can be written as with It follows that all roots of (43) have negative real parts for  2 (0,  * 2 ) when   ≥ 0,  = 1, 2, 3, 4 and this implies that endemic equilibrium is locally asymptotically stable for  2 ∈ (0,  * 2 ).This completes the proof.We now explore the global stability of the endemic equilibrium.Theorem 6.If R 0 > 1, then E * is globally asymptotically stable.
Proof.Let us consider the Lyapunov function: Here, The derivatives of W 1 () are given by Substituting the appropriate differentials from (1) we have At endemic equilibrium, we have Using the above constants, we have ) . (52) The derivatives of W 2 are given by ) . ( The derivatives of W 3 () are given by ) . ( Combining the derivatives of W  (), for  = 1, 2, 3, we have for all () > 0 and () > 0, because the arithmetic mean is greater than or equal to the geometric mean.Further, note that () = 1 − () + ln () is always nonpositive for any function () > 0, and () = 0 if and only if () = 1.

Hopf Bifurcation Analysis.
In this section we determine criteria for Hopf bifurcation to occur using the time delay  1 and  2 as the bifurcation parameters to find the interval in which the infected equilibria is stable and unstable out of the same margins.Now to consider Hopf bifurcation we consider the cases (a)  1 =  10 > 0,  2 = 0 and (b)  2 =  20 > 0 and  1 = 0. Our analysis is as follows: (a) When  1 =  10 > 0 and  2 = 0 we need to show that Re( 10 )/ 1 > 0 differentiating (37) with respect to  1 we get This gives with In our case considering () =  4 +  1  3 +  2  2 +  3  +  4 = 0 defined in (42), and assuming  3 < 0 and  2 0 as the largest positive root, we have The above analysis can be summarized into the following theorem.
Theorem 8. Suppose that () R 0 > 1.If either ()  4 < 0 or ()  4 ≥ 0 and  3 < 0 is satisfied, and  0 is the largest positive simple root of (42) then the infected equilibrium E * of model ( 1) is locally asymptotically stable when  1 <  10 and unstable when  1 >  10 where when  1 =  10 , a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from E * as  1 passes through the critical value  10 .
(b) When  2 =  20 > 0 and  1 = 0 we also need to show that Re( 20 )/ 2 > 0 differentiating (43) with respect to  2 we get (4 From the analysis above, we can deduce that Hopf bifurcations may arise if conditions in Theorems 8 and 10 are satisfied.Thus, the introduction of time delay in system (1) can destabilize the system.

Numerical Results
In order to explore the behavior of system (1) and illustrate the stability of the equilibria solutions, we numerically solve system (1) using MATLAB and parameter values adopted from Table 1.
In Figure 1 we illustrate the effects of varying the delay ( 1 =  2 ) on the dynamics of system (1).Figures 1(a) and 1(b)  Figure 1: Stability of the infected and free-infected equilibrium of model system (1) showing plots of  1 () and  2 () with varying delay ( 1 =  2 ).The direction of the arrow depicts an increase in delay with a step size of 2.0 starting from 2.0 to 10.The blue patterns in both (a) and (b) highlight brucellosis dynamics when R 0 < 1 while the red pattern is for R 0 > 1.Initial population levels were assumed as follows: (0) = 1000 animals, (0) = 500 animals,  1 (0) = 500 animals, and  2 (0) = 0 animals.demonstrate that the system approaches the stable diseasefree or endemic equilibrium for R 0 < 1 and R 0 > 1, respectively.One should note that, according to Theorems 2 and 6, the stability of the model steady states does not depend on the value of the time delays, but rather on the basic reproduction number R 0 , only.In addition, we observe that the range of values for the two time delays does not lead to periodic solutions but an increase in both delays translate to an increase in the infectious population, both detected and undetected.Figure 2 depicts the numbers of infectious undetected and infectious detected animals over time with varying delays.The results clearly show that the incubation related delay ( 1 ) has more influence on shaping the dynamics of brucellosis compared to the culling related delay ( 2 ).More precisely, the incubation period delay significantly increases the infectious population (both detected and undetected) for 0 <  < 20 and thereafter its impact will be the same as that of detection ( 2 ).
In Figure 3 we illustrate the stability of the disease-free equilibrium E 0 with  1 = 30 and  2 = 5 (note that R 0 = 0.686281).As we can observe, for certain parameter values and initial population levels, system (1) exhibits some periodic oscillation.Precisely, we note that the infected population ( 1 () and  2 ()) oscillates with a reduced amplitude from the start till when  is approximately 400; thereafter the oscillations dies off the solutions converges to the diseasefree equilibrium.These simulation results demonstrate the occurrence of periodic solutions through Hopf bifurcation for delay values  1 = 30 and  2 = 5.In contrast, we can observe that there are no periodic oscillations for the uninfected populations () and ().
equilibrium for different pair of delay values ( 1 ,  2 ) and from the simulation results we can conclude that both delays do not have a huge influence on the stability of disease-free equilibrium.
In Figure 5 we observe that, for certain parameter values and initial population levels, system (1) may admit periodic oscillations when R 0 > 1.As we can observe, when R 0 = 3.77333 both the solutions of the infected and uninfected populations exhibit periodic oscillation for a certain period, before stability at endemic point is attained.
In Figure 6, we illustrate the dynamics for model system (1) with respect to the stability of endemic equilibrium for several pair of delay values ( 1 ,  2 ).The results confirm that the incubation related delay  1 has more influence on shaping the dynamic of brucellosis compared to the culling related delay  2 .
To explore influence of model parameters on the reproduction number R 0 , we perform a local sensitivity analysis of the basic reproduction number following the approach in [23].The local sensitivity analysis will be useful on identifying parameters with greatest influence to change R 0 .To this end, denoting by Φ the generic parameter of system (1), we evaluate the normalised sensitivity index: which indicates how sensitive R 0 is to a change of parameter Φ. Model parameters with positive index increase the value of R 0 whenever they are increased while those with a negative index decrease the value of R 0 whenever they are increased.We consider the parameter values in Table 1, and we set  = 0.015 in order to evaluate the normalized sensitivity index and the results are depicted in Figure 7. Here, we observe that parameters , , ,  have a positive correlation with R 0 , such that increasing these parameters will increase R 0 .However, it is the increase of  and  that has the greatest influence to change R 0 .Precisely, increasing either  or  by 50% will increase R 0 by 50%.We also note that increasing parameters, , , , , , and , will lower the reproduction number.Figure 5: Stability of E * equilibrium of model system (1) with R 0 = 3.77333.The time delay  1 was fixed to be 60 and  2 was fixed to be 1.We set the model parameters and variables as follows:  = 6.844 × 10 −6 animal −1 year −1 ,  = 0.2,  = 0.015 year −1 , (0) = 100 animals, (0) = 0 animals,  1 (0) = 10 animals, and  2 (0) = 0 animals while the other parameter values are as in Table 1.

Conclusion
Zoonotic brucellosis remains a major public health problem in many developing nations.This is mainly attributed to several challenges associated with effective disease control in these nations.The challenges for effective control of brucellosis in developing nations range from inadequate veterinary personnel and vaccines as well as the failure by farmers to adhere to some of the aforementioned brucellosis control and eradication program activities.Furthermore, these challenges often lead to delay in detection and culling of infectious animals.In this study, we developed and analysed a mathematical model for brucellosis infection that incorporates two discrete delays.The first delay accounts for the latent period and the second delay represents the time taken to detect infectious animals.We computed the basic reproduction number and demonstrated that it is an important threshold quantity for stability of equilibria.By constructing suitable Lyapunov functionals, it has been shown that the model has a globally asymptotically stable infection-free equilibrium whenever the reproduction is less than unity.Further, it has been demonstrated that whenever the model reproduction number is greater than unity then the model has a unique endemic equilibrium point which is globally asymptotically stable.Numerical simulations are carried out to illustrate the main results.
studies suggest that culling of both infected and susceptible animal may be more effective [24,25].The rationale is that by decreasing the host density, then the number of contacts per unit time between animals is low, thereby reducing disease transmission.In [25] it was demonstrated that culling of both susceptible and symptomatic animals only can be effective whenever the number of infected host is above a certain critical level [25].We expect to improve this study in our future work by developing brucellosis model(s) with time delay that will enable us to compare aforementioned aspects.

Figure 7 :
Figure 7: Sensitivity index for R 0 with respect to model parameters that define it.

Table 1 :
Model parameters and variables and their baseline values.