Modeling the Differential Effect of Prescribed Fire on Multi-Vector Tick-Borne Tularemia Disease

1. Introduction

Tick-borne diseases pose a growing public health concern, with their incidence on the rise in various regions around the world; generally tick-borne diseases, doubled from 2006 to 2018 [18,41]. For instance in the United States over 48,000 cases of tick-borne illnesses were reported in 2016 with close to 51,000 cases reported in 2019 [19]. These figures include reported cases from tick-borne diseases such as Lyme disease, anaplasmosis, Rocky Mountain spotted fever (RMSF), Babesiosis, ehrlichia chaffeensis ehrlichiosis, and Tularemia among others. The number of human Tularemia reported cases has been increasing in recent times in several countries from countries in North America (Canada, Mexico, United States), to several Nordic countries (Denmark, Norway, Sweden, Finland, and Iceland), and some Asian countries [14,87]. For instance, in Europe close to 1,500 cases of Tularemia disease were reported in 2019 with over 1,000 cases reported in 2016 [1]. Francisella tularensis is mostly distributed in the Northern hemisphere and is not normally found in the southern hemisphere or the tropics [1].

Tularemia has recently become a significant re-emerging disease in the world. Its low infectious dose, easy, high dissemination with aerosols, and the ability to induce fatal disease make Francisella tularensis a potential agent of biological warfare. This is why the Centers for Disease Control and Prevention (CDC) classified it as a category A biological weapon (Dennis et al. 2001; Maurin 2015) [87].

Tularemia is a zoonotic disease caused by Francisella tularensis bacteria [1,44,50,87], a gram-negative coccobacillus-shaped bacterium; it is also known as deer fly fever, Ohara’s fever, rabbit fever and Pahvant Valley plague [42,87]. There are four subspecies of Francisella tularensis namely tularensis (type A), holarctica (type B), novicida and mediasiatica [50,87]. The most virulent subspecies are tularensis and holarctica [1,32,50,53,62]. The subspecies tularensis, the most virulent of the two virulent subspecies, only occur in North America, while the subspecies holarctica is the most widespread of the subspecies, mediasiatica is present in central Asia, and novicida is the least virulent [1]. The subspecies tularensis and holarctica are the major etiological agents of tularemia in humans [87]. While the subspecies novicida is rarely associated with human infections [74,87], the subspecies mediasiatica have never been documented in published literature to cause infection in humans [87].

Transmission of francisella tularensis to humans occur via multiple routes, such as consumption of contaminated food or water, handling of infected animals or bites from haematophagous arthropod vectors (such as ticks, deer flies, or mosquitoes) [1,24,45]; other routes include contact with aquatic environment, and inhalation of aerosols [8,48,68,87]. Once an individual is infected with francisella tularensis, the incubation period ranges between 1–21 days [17,22,26,75,76,87]. The symptoms of tularemia include fever, fatigue, chills, and headache. The clinical manifestation of tularemia in humans depends on the site and route of acquisition of the infection[27,87]. There are six clinical forms of infection in humans, these include ulcero-glandular, glandular, oculo-glandular, oro-pharyngeal, pneumonic, and typhoidal forms [23,27,87].

Ulcero-glandular tularemia develops as a result of direct contact with an infected animal or a vector bite such as tick or deer fly; it leads to having skin lesions and lymphadenopathy. Although glandular and ulceroglandular tularemia are similar; however, glandular differs from ulceroglandular with the presence of regional lymphadenopathy with no detectable skin lesion. The oculo-glandular tularemia often develops via direct contact with contaminated water or splash of infected animal’s body fluids into the conjunctiva. Oro-pharyngeal result from ingestion of the bacterium and lead to symptoms such as pharyngitis, fever, and cervical lymphadenitis appear [70]. Pneumonic tularemia develops from the inhalation of infectious aerosols, while typhoidal tularemia develops from ingesting contaminated food and water. Moreover, pneumonic and typhoidal tularemia are considered systemic forms as they develop by the spread of bacteria through blood circulation as a systemic disease [87]. All six clinical forms of tularemia can cause secondary pleuropneumonia, meningitis, and sepsis, shock and death in infected individuals [4]. Tularemia is often diagnosed through the use of culture, serology, or molecular methods. The disease is frequently treated using quinolones, tetracyclines, or aminoglycosides. There are no licensed vaccine available for the prevention of tularemia [79,87].

The bacterium has a wide host range including different vertebrate groups as well as invertebrates [50,71]. However, rodents, hares, and rabbits are the principal vertebrate hosts [1], while ticks are the principal enzootic vector and reservoir and major sources of infection in human [1,61,87]. According to the Centers for Disease Control and Prevention (CDC), in the United States, tularemia disease can be transmitted to humans by three tick species namely Amblyomma americanum (lone Star tick), Dermacentor variabilis (american Dog tick) and Dermacentor andersoni (wood tick) [24]. In Europe, tularemia is transmitted by Dermacentor reticulatus [1]. Ticks goes through four developmental stages from eggs to larva to nymph, and then adult [6]. During each life stage, the ticks needs a blood meal to develop into the next stage of their life cycle [12].

Fire is considered a basic ecological process that keeps a variety of vegetative communities intact [47]. Prescribed fire or controlled burn is a common and essential land management technique in ecosystems that depend on fire, such as grasslands, open pine forests, and wetlands [31,83]. Prescribed fire are fire set by a group of professionals under specific weather conditions to restore health to ecosystems that depend on fire [13,81]. Prescribed fires are commonly used to control tick population in different environments by killing the ticks directly and destroying their leaf litter habitats [39]. Given the increase in recent decades of tick-borne diseases and the discovery of several new pathogens [59,73], it is necessity to identify cost-effective and practical methods for reducing the risk of tick-borne diseases. Gleim et al. [39], found that extended prescribed fires not only significantly decreased tick abundance but also altered the composition of tick species.

Our goal in this is to investigate the differential effect of fire on the tick -borne disease via a mathematical model. However, fewer modeling work have been done to consider the effect of prescribed fire on tick-borne diseases [35,41]. Guo and Agusto [41], investigated the effect of high and low intensity of prescribed fire on a single-vector tick-borne disease. The study presented a compartmental model for ticks carrying Lyme disease and incorporated the effects of prescribed fire using an impulsive system. They recommended that high-intensity prescribed burns done annually resulted in the most significant reduction in infectious nymphs, the primary carriers of Lyme disease [41]. Fulk et al. [35] also explored the effect of prescribed burns on the tick populations and the spread of Lyme disease by developing a spatial stage-structured tick-host model to simulate the impact of controlled burning on tick populations. The numerical simulations explained the effects of different numbers of burns and the time between burns on tick populations. However, it was again found that consistent prescribed burning at high intensity was the most effective control method for tick populations [35]. Furthermore, Fulk and Agusto [34] showed significant increase in the number of Ehrlichiosis infected Amblyomma americanum ticks at temperature increases in the absence of prescribed burning. However with prescribed burning, they observed a reduction in the prevalence of infectious ticks even as temperature increases to level such as $2^0$C and $4^0$C.

In this research work, we develop a compartmental model for a two-ticks two-hosts system with both tick species carrying Francisella tularensis and we investigate the differential effect of fire on the ticks and the prevalence of tularemia disease on the hosts. We are not aware of any studies considering the transmission of a single pathogen by multiple vectors. Although several studies have been done on a single-vector single-pathogen system [2,36,37,56,64,78,86] and co-infection of multiple pathogens in a single tick vector [54,77,80]. Also, some studies have considered single vector-multiple host system. For instance, Occhibove et al. [66] considered a single-vector, multi-host model for two tick species (Ixodes ricinus and Ixodes trianguliceps) individually infected with Borrelia burgdorferi and Babesia microti; the model further incorporated the ecological relationships with non-host species. The study aimed to understand the dynamics of tick-borne pathogens in pathosystems that differ in vector generalism. The study found that the degree of vector generalism affected pathogen transmission with different dilution outcomes. Norman et al. [65] also considered a two-hosts one-tick system with one host viraemic and the other not. They found that the non-viraemic hosts could either amplify the tick population leading to the persistence of the virus, or they may dilute the infection and cause it to die out.

The remainder of the work in this paper is organized as follows. In Section 2, we formulate our baseline tick/tularemia disease model, compute the model basic reproduction number, and carry out a global stability analysis including sensitivity analysis to determine the parameter with the most impact on the response function (basic reproduction number). In Section 3, we describe the tick model with the effect of prescribed fire using an impulsive system of ordinary differential equations and discuss prescribed fire related parameters estimated from literature. In Section 4, we carry out a global sensitivity analysis of the developed models. In Section 5, we present some simulation results, and in Section 6 we discuss our findings and close with conclusions.

2. Model Formulation

According to the Centers for Disease Control and Prevention, tularemia disease humans become infected when they are bitten by vector bite such as ticks or deer flies or by handling infected or deceased animals, or by the consumption of contaminated food or water. In this section, we formulate a model using non-linear ordinary differential equation that accounts for the interplay between humans (a dead end host, there is no evidence of human-to-human transmission), the American dog tick (referred to as tick 1), the lone star tick (referred to as tick 2) and rodents (pathogen reservoir). This model incorporates various compartments that depict the epidemiological state of each species in the system.

In the model, the human population is divided into susceptible $S_H\left(t\right),$ exposed $E_H\left(t\right),$ asymptomatic $A_H\left(t\right),$ infected $I_H\left(t\right)$ and recovered $R_H\left(t\right).$ The total human population is denoted by $N_H\left(t\right)$ and is defined as

$$N_H=S_H+E_H+A_H+I_H+R_H.$$

The total rodent population is denoted by $N_M\left(t\right)$ and is defined as

$$N_M=S_M+I_M.$$

The tick population is divided by life stages, which include eggs, larvae, nymphs, and adults. Additionally, it is separated into susceptible and infected groups for larvae ($S_{L1}\left(t\right)$ and $I_{L1}\left(t\right)$) and ($S_{L2}\left(t\right)$ and $I_{L2}\left(t\right)$), nymphs ($S_{N1}\left(t\right)$ and $I_{N1}\left(t\right)$) and ($S_{N2}\left(t\right)$ and $I_{N2}\left(t\right)$), and adults ($S_{A1}\left(t\right)$ and $I_{A1}\left(t\right)$) and $S_{A2}\left(t\right)$ and $I_{A2}\left(t\right)$ for tick 1 and tick 2, respectively. Since ticks must feed on blood to get infected, and there is no vertical transmission of the disease from parent ticks to their offspring, all tick eggs remain infection free ($S_{E1}\left(t\right)$ and $S_{E2}\left(t\right)$).

For simplicity, we assume that human, tick and rodent populations are mixing homogeneously. The susceptible human sub-population are recruited at a rate ${\pi }_H.$ The rate at which a susceptible human progress to the exposed class is denoted by ${\lambda }_H$ and is defined as

$${\lambda }_H={\displaystyle \frac{{\beta }_{TH1}\left(I_{L1}+I_{N1}+I_{A1}\right)+{\beta }_{TH2}\left(I_{L2}+I_{N2}+I_{A2}\right)}{N_H}}.$$

This quantity ${\lambda }_H$ is the force of infection for the humans, where ${\beta }_{TH1}$ and ${\beta }_{TH2}$ are the probabilities of transmission from tick 1 to human and tick 2 to human, respectively.

The incubation period for tularemia in humans ranges 1–21 days [17,22,26,75,76,87]. Furthermore, 4–19% of the infected human cases are asymptomatic [9,42]. Thus, the exposed human subsequently becomes either asymptomatic at a rate $p{\sigma }_H$ or infected at a rate $(1-p){\sigma }_H$ with $p<1$. Tularemia in humans is usually treated with antibiotics such as streptomycin, gentamicin, doxycycline, and ciprofloxacin [21]. Depending on the stage of illness and the type of medication used, treatment can last as long as 10–21 days with many of the patients completely recovering from the illness [21,49]. Hence, we assume that individuals who are infected as well as those who are asymptomatic both recover at the rate ${\gamma }_H$ according to the reciprocal of the length of the treatment days. The human population is decreased by natural death at the rate denoted by ${\mu }_H$ and by disease induced mortality denoted by ${\delta }_H$. The disease induced mortality rate for F. tularensis tularensis infections is about 30–35% and about 5–15% for F. tularensis holarctica infections when left untreated [49,60,76]. When treated with the appropriate antibiotics, these figure reduces to 1–3% [49,60]. The mortality rate in patients with typhoidal tularemia is higher than those of other forms of tularemia [60]. The case fatality rate of typhoidal tularemia if untreated is approximately 35%; it is the most dangerous form of tularemia infections [60]. Untreated ulceroglandular tularemia infection has a case fatality rate of 5% [60]. Permanent immunity usually develops after a single episode of tularemia [60].

For the rodent population, the recruited at a rate into the susceptible sub-population is denoted by ${\pi }_M$. The rodents become infected following a bite from infected ticks of type 1 and 2 at the rate

$${\lambda }_M={\displaystyle \frac{{\beta }_{TM1}\left(I_{L1}+I_{N1}+I_{A1}\right)+{\beta }_{TM2}\left(I_{L2}+I_{N2}+I_{A2}\right)}{N_H}}.$$

where ${\beta }_{TM1}$ and ${\beta }_{TM2}$ are the probabilities of transmission from tick 1 to rodents and tick 2 to rodents, respectively. Rabbits, hares, and rodents often die in large numbers during outbreaks due to their high susceptibility to the bacteria [20]. The rodent population therefore is decreased by disease induced mortality denoted by ${\delta }_M$ and by natural death at the rate by ${\mu }_M$.

For ticks, we assume that every adult tick can reproduce, regardless of whether they are susceptible or already infected, and that there is no vertical transmission of the bacteria to the eggs from infected females. The eggs develop into the larvae at rates ${\tau }_{E1}$ and ${\tau }_{E2}$, with a certain proportion dying naturally at rates ${\mu }_{E1}$ and ${\mu }_{E2}$ for tick 1 and tick 2, respectively. The susceptible larvae of both tick 1 and tick 2 become infected larvae at a rate of ${\lambda }_{T1}$ and ${\lambda }_{T2},$ respectively or remain in the susceptible larval compartment. Thus, the force of infection in tick 1 and tick 2 (or the rate at which the ticks become infected after feeding on infected rodents) is given as

$${\lambda }_{T1}={\displaystyle \frac{{\beta }_{M1}{I}_H}{N_M}}\mathrm{and}{\lambda }_{T2}={\displaystyle \frac{{\beta }_{M2}{I}_M}{N_M}},$$

where ${\beta }_{M1}$ and ${\beta }_{M2}$ represents the probabilities that infections will occur when tick 1 or tick 2 bites an infected rodents. The populations of both susceptible and infected larvae are influenced by the natural mortality rate of ${\mu }_{L1}$ for tick 1 and ${\mu }_{L2}$ for tick 2. Regardless of their infection status, both susceptible and infected larvae from both tick species mature into their respective nymph stages at the rates ${\tau }_{L1}$ and ${\tau }_{L2}.$ Following a subsequent infectious blood meal, susceptible nymphs in tick 1 and tick 2 become infected nymph at a rate of ${\lambda }_{T1}$ and ${\lambda }_{T2},$ respectively. The populations of both susceptible and infected nymphs from both tick species are reduced due to natural death rate of ${\mu }_{N1}$ and ${\mu }_{N2},$ and they mature into adults at the rate ${\tau }_{N1}$ and ${\tau }_{N2},$ respectively.

Given the assumptions above, the following nonlinear equations are given for the transmission dynamics of tularemia in two tick species:

$$\left.\begin{array}{ccc}{\displaystyle \frac{d{S}_H}{dt}}\hfill & =& {\pi }_H-{\lambda }_H{S}_H-{\mu }_H{S}_H\hfill \\ {\displaystyle \frac{d{E}_H}{dt}}\hfill & =& {\lambda }_H{S}_H-({\sigma }_H+{\mu }_H)E_H\hfill \\ {\displaystyle \frac{d{A}_H}{dt}}\hfill & =& p{\sigma }_H{E}_H-({\gamma }_H+{\mu }_H)A_H\left[5pt\right]\hfill \\ {\displaystyle \frac{d{I}_H}{dt}}\hfill & =& (1-p){\sigma }_H{E}_H-({\gamma }_H+{\mu }_H+{\delta }_H)I_H\hfill \\ {\displaystyle \frac{d{R}_H}{dt}}\hfill & =& {\gamma }_H{A}_H+{\gamma }_H{I}_H-{\mu }_H{R}_H\hfill \\ {\displaystyle \frac{d{S}_M}{dt}}\hfill & =& {\pi }_m-{\lambda }_M{S}_M-{\mu }_M{S}_M\hfill \\ {\displaystyle \frac{d{I}_M}{dt}}\hfill & =& {\lambda }_M{S}_M-({\mu }_M+{\sigma }_H)I_M\hfill \\ {\displaystyle \frac{d{S}_{E_1}}{dt}}\hfill & =& {\pi }_{T1}\left(1-{\displaystyle \frac{S_{E_1}}{K_1}}\right)(S_{A_1}+I_{A_1})-({\tau }_{E1}+{\mu }_{E1})S_{E1}\hfill \\ {\displaystyle \frac{d{S}_{L1}}{dt}}\hfill & =& {\tau }_{E_1}{S}_{E1}-{\lambda }_{T1}{S}_{L1}-({\tau }_{L1}+{\mu }_{L1})S_{L1}\hfill \\ {\displaystyle \frac{d{I}_{L_1}}{dt}}\hfill & =& {\lambda }_{T1}{S}_{L1}-({\tau }_{L1}+{\mu }_{L1})I_{L1}\hfill \\ {\displaystyle \frac{d{S}_{N1}}{dt}}\hfill & =& {\tau }_{L1}{S}_{L1}-{\lambda }_{T1}{S}_{N1}-({\tau }_{N1}+{\mu }_{N1})S_{N1}\hfill \\ {\displaystyle \frac{d{I}_{N1}}{dt}}\hfill & =& {\tau }_{L1}{I}_{L1}+{\lambda }_{T1}{S}_{N1}-({\tau }_{N1}+{\mu }_{N1})I_{N1}\hfill \\ {\displaystyle \frac{d{S}_{A_1}}{dt}}\hfill & =& {\tau }_{N1}{S}_{N1}-{\lambda }_{T1}{S}_{A1}-{\mu }_{A1}{S}_{A1}\hfill \\ {\displaystyle \frac{d{I}_{A1}}{dt}}\hfill & =& {\tau }_{N1}{I}_{N_1}+{\lambda }_{TI}{S}_{A1}-{\mu }_{A1}{I}_{A1}\hfill \\ {\displaystyle \frac{d{S}_{E_2}}{dt}}\hfill & =& {\pi }_{T2}\left(1-{\displaystyle \frac{S_{E_2}}{K_2}}\right)(S_{A_2}+I_{A_2})-\left({\tau }_{E2}+{\mu }_{E2}\right)S_{E2}\hfill \\ {\displaystyle \frac{d{S}_{L2}}{dt}}\hfill & =& {\tau }_{E2}{S}_{E2}-{\lambda }_{T2}{S}_{L2}-({\tau }_{L2}+{\mu }_{L2})S_{L2}\hfill \\ {\displaystyle \frac{d{I}_{L2}}{dt}}\hfill & =& {\lambda }_{T2}{S}_{L2}-({\tau }_{L2}+{\mu }_{L2})I_{L2}\hfill \\ {\displaystyle \frac{d{S}_{N2}}{dt}}\hfill & =& {\tau }_{L2}{S}_{L2}-{\lambda }_{T2}{S}_{N2}-({\tau }_{N2}+{\mu }_{N2})S_{N2}\hfill \\ {\displaystyle \frac{d{I}_{N2}}{dt}}\hfill & =& {\tau }_{L2}{I}_{L2}+{\lambda }_{T2}{S}_{N2}-({\tau }_{N2}+{\mu }_{N2})I_{N2}\hfill \\ {\displaystyle \frac{d{S}_{A2}}{dt}}\hfill & =& {\tau }_{N2}{S}_{N2}-{\lambda }_{T2}{S}_{A2}-{\mu }_{A2}{S}_{A2}\hfill \\ {\displaystyle \frac{d{I}_{A2}}{dt}}\hfill & =& {\tau }_{N2}{I}_{N2}+{\lambda }_{T2}{S}_{A2}-{\mu }_{A2}{I}_{A2}\hfill \end{array}\right\}$$

(1)

The flow diagram for tularemia transmission is shown in Figure 1. The corresponding parameters and variables are described in Table 1 and Table 2.

2.1. Analysis of the Model

2.1.1. Basic Qualitative Properties

Positivity and Boundedness of Solutions

For the tularemia model (1) to be epidemiologically meaningful, it is important to prove that all its state variables are non-negative for all time. In other words, solutions of the model system (1) with non-negative initial data will remain non-negative for all time $t>0$.

Lemma 1. 

Let the initial data $F\left(0\right)\ge 0$ , where $F\left(t\right)=(S_H\left(t\right),E_H\left(t\right),A_H\left(t\right),I_H\left(t\right),R_H\left(t\right),$$S_M\left(t\right),I_M\left(t\right),S_{Ei}\left(t\right),S_{Li}\left(t\right),I_{Li}\left(t\right),S_{Ni}\left(t\right),I_{Ni}\left(t\right),S_{Ai}\left(t\right),I_{Ai}\left(t\right))$, where $i=1,2$. Then the solutions $F\left(t\right)$ of the tularemia model (1) are non-negative for all $t>0$. Furthermore

$$\underset{t\to \infty }{lim\; sup}{N}_H\left(t\right)\le \frac{{\pi }_H}{{\mu }_H},\underset{t\to \infty }{lim\; sup}{N}_M\left(t\right)\le \frac{{\pi }_M}{{\mu }_M},and\underset{t\to \infty }{lim\; sup}{N}_{Ti}\left(t\right)\le \frac{{\pi }_{Ti}}{{\mu }_{Ti}}.$$

where

$$N_H\left(t\right)=S_H\left(t\right),E_H\left(t\right),A_H\left(t\right),I_H\left(t\right),R_H\left(t\right),N_M\left(t\right)=S_M\left(t\right),I_M\left(t\right),$$

and

$$N_{Ti}\left(t\right)=S_{Ei}\left(t\right),S_{Li}\left(t\right),I_{Li}\left(t\right),S_{Ni}\left(t\right),I_{Ni}\left(t\right),S_{Ai}\left(t\right),I_{Ai}\left(t\right).$$

The proof of Lemma 1 is given in Appendix A.

Invariant Regions

The tularemia model (1) will be analyzed in a biologically-feasible region as follows. Consider the feasible region

$$\Omega ={\Omega }_H\cup {\Omega }_M\cup {\Omega }_{T1}\cup {\Omega }_{T2}\subset {\mathbb{R}}_{+}^5\times {\mathbb{R}}_{+}^2\times {\mathbb{R}}_{+}^7\times {\mathbb{R}}_{+}^7,$$

with,

$${\Omega }_H=\left\{(S_H\left(t\right),E_H\left(t\right),A_H\left(t\right),I_H\left(t\right),R_H\left(t\right))\in {\mathbb{R}}_{+}^5:N_H\left(t\right)\le \frac{{\pi }_H}{{\mu }_H}\right\},$$

$${\Omega }_M=\left\{(S_M\left(t\right),I_M\left(t\right))\in {\mathbb{R}}_{+}^2:N_M\left(t\right)\le \frac{{\pi }_M}{{\mu }_M}\right\},$$

and

$${\displaystyle {\Omega }_{Ti}=\left\{(S_{Ei}\left(t\right),S_{Li}\left(t\right),I_{Li}\left(t\right),S_{Ni}\left(t\right),I_{Ni}\left(t\right),S_{Ai}\left(t\right),I_{Ai}\left(t\right))\in {\mathbb{R}}_{+}^7:N_{Ti}\left(t\right)\le \frac{{\pi }_{Ti}}{{\mu }_{Ti}}\right\}.}$$

Lemma 2. 

The region $\Omega ={\Omega }_H\cup {\Omega }_M\cup {\Omega }_{T1}\cup {\Omega }_{T2}\subset {\mathbb{R}}_{+}^5\times {\mathbb{R}}_{+}^2\times {\mathbb{R}}_{+}^7\times {\mathbb{R}}_{+}^7$ is positively-invariant for the model (1) with non-negative initial conditions in ${\mathbb{R}}_{+}^{21}$.

The prove of Lemma 2 is given in Appendix B. In the next section, the conditions for the existence and stability of the equilibria of the model are stated.

2.1.2. Stability of Disease-Free Equilibrium (DFE)

The tularemia model has a disease-free equilibrium (DFE). The DFE is obtained by setting the right-hand sides and infected variables of the equation (2) to zero. The system has four equilibria ${\mathcal{E}}_0=({\mathcal{E}}_0^1,{\mathcal{E}}_0^2,{\mathcal{E}}_0^3,{\mathcal{E}}_0^4)$. The equilibrium ${\mathcal{E}}_0^1$, involves humans and rodents only; while ${\mathcal{E}}_0^2$ involves, humans, rodents and Demacenta variablis only; the equilirium ${\mathcal{E}}_0^3$ involves humans, rodents and Amblyomma americanum only. Lastly, the ${\mathcal{E}}_0^4$ involves humans, rodents and the two tick species. The expression for these equiliria are given as

$$\begin{array}{ccc}\hfill {\mathcal{E}}_0^1& =& (S_H^{*},S_M^{*},0,0,0,0,0,0,0,0\hfill \\ & =& \left({\displaystyle \frac{{\pi }_H}{{\mu }_H}},{\displaystyle \frac{{\pi }_M}{{\mu }_M}},0,0,0,0,0,0,0,0\right),\hfill \\ \hfill {\mathcal{E}}_0^2& =& \left(S_H^{*},S_M^{*},S_{E1}^{*},S_{L1}^{*},S_{N1}^{*},S_{A1}^{*},0,0,0,0\right)\hfill \\ \\ \hfill S_H^{*}& =& {\displaystyle \frac{{\pi }_H}{{\mu }_H}},S_M^{*}={\displaystyle \frac{{\pi }_M}{{\mu }_M}},\hfill \\ \hfill S_{E1}^{*}& =& {\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{{\tau }_{N1}{\tau }_{L1}{\tau }_{E1}{\pi }_{T1}}},S_{L1}^{*}={\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{\tau }_{N1}{\tau }_{L1}{\pi }_{T1}}},\hfill \\ \hfill S_{N1}^{*}& =& {\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{n}_3{\pi }_{T1}{\tau }_{N1}}},S_{A1}^{*}={\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{n}_3{\pi }_{T1}{\mu }_{A1}}},\hfill \\ \\ \hfill {\mathcal{E}}_0^3& =& (S_H^{*},S_M^{*},0,0,0,0,S_{E2}^{*},S_{L2}^{*},S_{N2}^{*},S_{A2}^{*})\hfill \\ \hfill S_H^{*}& =& {\displaystyle \frac{{\pi }_H}{{\mu }_H}},S_M^{*}={\displaystyle \frac{{\pi }_M}{{\mu }_M}},\hfill \\ \hfill S_{E2}^{*}& =& {\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{{\tau }_{N2}{\tau }_{L2}{\tau }_{E2}{\pi }_{T2}}},S_{L2}^{*}={\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{\tau }_{N2}{\tau }_{L2}{\pi }_{T2}}},\hfill \\ \hfill S_{N2}^{*}& =& {\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{n}_6{\pi }_{T2}{\tau }_{N2}}},S_{A2}^{*}={\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{n}_6{\pi }_{T2}{\mu }_{A2}}},\hfill \\ \hfill {\mathcal{E}}_0^4& =& \left(S_H^{*},S_M^{*},S_{E1}^{*},S_{L1}^{*},S_{N1}^{*},S_{A1}^{*},S_{E2}^{*},S_{L2}^{*},S_{N2}^{*},S_{A2}^{*}\right)\hfill \\ \\ \hfill S_H^{*}& =& {\displaystyle \frac{{\pi }_H}{{\mu }_H}},S_M^{*}={\displaystyle \frac{{\pi }_M}{{\mu }_M}},\hfill \\ \hfill S_{E1}^{*}& =& {\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{{\tau }_{N1}{\tau }_{L1}{\tau }_{E1}{\pi }_{T1}}},S_{L1}^{*}={\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{\tau }_{N1}{\tau }_{L1}{\pi }_{T1}}}\hfill \\ \hfill S_{N1}^{*}& =& {\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{n}_3{\pi }_{T1}{\tau }_{N1}}},S_{A1}^{*}={\displaystyle \frac{K_1({\pi }_{T1}{\tau }_{E1}{\tau }_{L1}{\tau }_{N1}-n_1{n}_2{n}_3{\mu }_{A1})}{n_2{n}_3{\pi }_{T1}{\mu }_{A1}}}\hfill \\ \hfill S_{E2}^{*}& =& {\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{{\tau }_{N2}{\tau }_{L2}{\tau }_{E2}{\pi }_{T2}}},S_{L2}^{*}={\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{\tau }_{N2}{\tau }_{L2}{\pi }_{T}2}}\hfill \\ \hfill S_{N2}^{*}& =& {\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{n}_6{\pi }_{T2}{\tau }_{N2}}},S_{A2}^{*}={\displaystyle \frac{K_2({\pi }_{T2}{\tau }_{E2}{\tau }_{L2}{\tau }_{N2}-n_4{n}_5{n}_6{\mu }_{A2})}{n_5{n}_6{\pi }_{T2}{\mu }_{A2}}},\hfill \end{array}$$

where $n_1={\tau }_{E1}+{\mu }_{E1}),n_2={\tau }_{L1}+{\mu }_{L1}),n_3={\tau }_{N1}+{\mu }_{N1},n_4={\tau }_{E2}+{\mu }_{E2},n_5={\tau }_{L2}+{\mu }_{L2},n_6={\tau }_{N2}+{\mu }_{N2}.$ We will focus on the the stability of ${\mathcal{E}}_0^4$, the stability of the other equiliria ${\mathcal{E}}_0=({\mathcal{E}}_0^1,{\mathcal{E}}_0^2,{\mathcal{E}}_0^3)$ can be determined using the same approach. Thus, the stability of ${\mathcal{E}}_0^4$ can be established using the next generation operator method on system (1). Taking $E_H,A_H,I_H,R_H,I_M,I_{L1},$$I_{N1},I_{A1},I_{L2},I_{N2}$, and $I_{A2}$ as the infected compartments and then using the aforementioned notation, the Jacobian F and V matrices for new infectious terms and the remaining transfer terms, respectively, are defined as:

$$F=\left(\begin{array}{cccccccccc}0& 0& 0& 0& {\beta }_{TH1}& {\beta }_{TH1}& {\beta }_{TH1}& {\beta }_{TH2}& {\beta }_{TH2}& {\beta }_{TH2}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\beta }_{TM1}& {\beta }_{TM1}& {\beta }_{TM1}& {\beta }_{TM2}& {\beta }_{TM2}& {\beta }_{TM2}\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M1}{S}_{L1}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M1}{S}_{N1}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M1}{S}_{A1}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M2}{S}_{L2}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M2}{S}_{N2}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{{\beta }_{M2}{S}_{A2}^{*}}{S_M^{*}}}& 0& 0& 0& 0& 0& 0\end{array}\right)$$

$$V=\left(\begin{array}{cccccccccc}{k}_1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -(1-p){\sigma }_H& k2& 0& 0& 0& 0& 0& 0& 0& 0\\ -p{\sigma }_H& 0& k_3& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& k_4& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& k_5& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& -{\tau }_{L1}& k_6& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -{\tau }_{N1}& {\mu }_{A1}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& k_7& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -{\tau }_{L2}& k_8& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -{\tau }_{N2}& {\mu }_{A2}\end{array}\right)$$

where, $k_1={\sigma }_H+{\mu }_H,k_2={\gamma }_H+{\mu }_H,k_3={\gamma }_H+{\mu }_H+{\delta }_H,k_4={\mu }_M]+{\delta }_M],k_5={\tau }_{L1}+{\mu }_{L1},k_6={\tau }_{N1}+{\mu }_{N1},k_7={\tau }_{L2}+{\mu }_{L2},k_8={\tau }_{N2}+{\mu }_{N2}$. Therefore, using the definition of the reproduction number ${\mathcal{R}}_0$ from [82] we have

$$\begin{array}{ccc}\hfill {\mathcal{R}}_0& =& \rho \left(F{V}^{-1}\right)=\sqrt{{\mathcal{R}}_{M1}{\mathcal{R}}_{T1}+{\mathcal{R}}_{M2}{\mathcal{R}}_{T2}}\hfill \end{array}$$

where $\rho$ is the spectral radius and

$$\begin{array}{ccc}\hfill {\mathcal{R}}_{M1}& =& {\displaystyle \frac{{\beta }_{TM1}}{k_4}},{\mathcal{R}}_{M2}={\displaystyle \frac{{\beta }_{TM2}}{k_4}},\hfill \\ \hfill {\mathcal{R}}_{T1}& =& {\displaystyle \frac{{\beta }_{M1}[S_{L1}^{*}(k_6{\mu }_{A1}+{\mu }_{A1}{\tau }_{L1}+{\tau }_{L1}{\tau }_{N1})+S_{N1}^{*}{k}_5({\mu }_{A1}+{\tau }_{N1})+S_{A1}^{*}{k}_5{k}_6]}{S_M^{*}{k}_5{k}_6{\mu }_{A1}}}\hfill \\ \hfill {\mathcal{R}}_{T2}& =& {\displaystyle \frac{{\beta }_{M2}[S_{L2}^{*}(k_8{\mu }_{A2}+{\mu }_{A2}{\tau }_{L2}+{\tau }_{L2}{\tau }_{N2})+S_{N2}^{*}{k}_7({\mu }_{A2}+{\tau }_{N2})+S_{A2}^{*}{k}_7{k}_8]}{S_M^{*}{k}_7{k}_8{\mu }_{A2}}}\hfill \end{array}$$

The expressions ${\mathcal{R}}_{M1}$ and ${\mathcal{R}}_{M2}$ are the average number of secondary infections in rodents due to infectious Demacenta variablis and Amblyomma americanum while ${\mathcal{R}}_{T1}$ and ${\mathcal{R}}_{T2}$ are the average number of infections in ticks 1 and 2 due to an infectious rodents. Furthermore, using Theorem 2 in [82], the following result is established.

Lemma 3. 

The disease-free equilibrium (DFE) of the tularemia model (1) is locally asymptotically stable (LAS) if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

The quantity ${\mathcal{R}}_0$ is basic reproduction number and it is defined as the average number of new infections that result from one infectious individual in a population that is fully susceptible. [5,30,46,82]. The epidemiological significance of Lemma 3 is that tularemia model (1) will be eliminated from within a herd if the reproduction number (${\mathcal{R}}_0$) can be brought to (and maintained at) a value less than unity.

3. Tick Model with Prescribed Fire

In this section, we consider the effect of fire on ticks and the rodent population, and the subsequent effect on the incidence and prevalence of the disease on humans. Note that fire is not explicitly incorporated into the tularemia model (1); rather, we consider the effect of fire on population size after the burns. To introduce the effect of prescribed fire into the tularemia disease model we have the following system of non-autonomous impulsive differential equations.

$$\left.\begin{array}{ccc}{\displaystyle \frac{d{S}_H}{dt}}\hfill & =& {\pi }_H-{\lambda }_H{S}_H-{\mu }_H{S}_H\hfill \\ {\displaystyle \frac{d{E}_H}{dt}}\hfill & =& {\lambda }_H{S}_H-({\sigma }_H+{\mu }_H)E_H\hfill \\ {\displaystyle \frac{d{A}_H}{dt}}\hfill & =& p{\sigma }_H{E}_H-({\gamma }_H+{\mu }_H)A_H\hfill \\ {\displaystyle \frac{d{I}_H}{dt}}\hfill & =& (1-p){\sigma }_H{E}_H-({\gamma }_H+{\mu }_H+{\delta }_H)I_H\hfill \\ {\displaystyle \frac{d{R}_H}{dt}}\hfill & =& {\gamma }_H{A}_H+{\gamma }_H{I}_H-{\mu }_H{R}_H\hfill \\ {\displaystyle \frac{d{S}_M}{dt}}\hfill & =& {\pi }_M-{\lambda }_M{S}_M-{\mu }_M{S}_M\hfill \\ {\displaystyle \frac{d{E}_M}{dt}}\hfill & =& {\lambda }_M{S}_M-({\mu }_M+{\sigma }_H)I_M\hfill \\ {\displaystyle \frac{d{S}_{E_1}}{dt}}\hfill & =& {\pi }_{T1}\left(1-{\displaystyle \frac{S_{E_1}}{K_1}}\right)(S_{A_1}+I_{A_1})-({\tau }_{E1}+{\mu }_{E1})S_{E1}\hfill \\ {\displaystyle \frac{d{S}_{L1}}{dt}}\hfill & =& {\tau }_{E_1}{S}_{E1}-{\lambda }_{T1}{S}_{L1}-({\tau }_{L1}+{\mu }_{L1})S_{L1}\hfill \\ {\displaystyle \frac{d{I}_{L_1}}{dt}}\hfill & =& {\lambda }_{T1}{S}_{L1}-({\tau }_{L1}+{\mu }_{L1})I_{L1}\hfill \\ {\displaystyle \frac{d{S}_{N1}}{dt}}\hfill & =& {\tau }_{L1}{S}_{L1}-{\lambda }_{T1}{S}_{N1}-({\tau }_{N1}+{\mu }_{N1})S_{N1}\hfill \\ {\displaystyle \frac{d{I}_{N1}}{dt}}\hfill & =& {\tau }_{L1}{I}_{L1}+{\lambda }_{T1}{S}_{N1}-({\tau }_{N1}+{\mu }_{N1})I_{N1}\hfill \\ {\displaystyle \frac{d{S}_{A_1}}{dt}}\hfill & =& {\tau }_{N1}{S}_{N1}-{\lambda }_{T1}{S}_{A1}-{\mu }_{A1}{S}_{A1}\hfill \\ {\displaystyle \frac{d{I}_{A1}}{dt}}\hfill & =& {\tau }_{N1}{I}_{N_1}+{\lambda }_{TI}{S}_{A1}-{\mu }_{A1}{I}_{A1}\hfill \\ {\displaystyle \frac{d{S}_{E_2}}{dt}}\hfill & =& {\pi }_{T2}\left(1-{\displaystyle \frac{S_{E_2}}{K_2}}\right)(S_{A_2}+I_{A_2})-\left({\tau }_{E2}+{\mu }_{E2}\right)S_{E2}\hfill \\ {\displaystyle \frac{d{S}_{L2}}{dt}}\hfill & =& {\tau }_{E2}{S}_{E2}-{\lambda }_{T2}{S}_{L2}-({\tau }_{L2}+{\mu }_{L2})S_{L2}\hfill \\ {\displaystyle \frac{d{I}_{L2}}{dt}}\hfill & =& {\lambda }_{T2}{S}_{L2}-({\tau }_{L2}+{\mu }_{L2})I_{L2}\hfill \\ {\displaystyle \frac{d{S}_{N2}}{dt}}\hfill & =& {\tau }_{L2}{S}_{L2}-{\lambda }_{T2}{S}_{N2}-({\tau }_{N2}+{\mu }_{N2})S_{N2}\hfill \\ {\displaystyle \frac{d{I}_{N2}}{dt}}\hfill & =& {\tau }_{L2}{I}_{L2}+{\lambda }_{T2}{S}_{N2}-({\tau }_{N2}+{\mu }_{N2})I_{N2}\hfill \\ {\displaystyle \frac{d{S}_{A2}}{dt}}\hfill & =& {\tau }_{N2}{S}_{N2}-{\lambda }_{T2}{S}_{A2}-{\mu }_{A2}{S}_{A2}\hfill \\ {\displaystyle \frac{d{I}_{A2}}{dt}}\hfill & =& {\tau }_{N2}{I}_{N2}+{\lambda }_{T2}{S}_{A2}-{\mu }_{A2}{I}_{A2}\hfill \end{array}\right\}{t}_n\ne nT,n\in \mathbb{N}$$

(2)

subject to the prescribed fire impulsive condition

$$\left.\begin{array}{c}{S}_H\left(n{T}^{+}\right)=S_H\left(n{T}^{-}\right),E_H\left(n{T}^{+}\right)=E_H\left(n{T}^{-}\right),A_H\left(n{T}^{+}\right)=A_H\left(n{T}^{-}\right),\hfill \\ I_H\left(n{T}^{+}\right)=I_H\left(n{T}^{-}\right),R_H\left(n{T}^{+}\right)=R_H\left(n{T}^{-}\right),\hfill \\ S_M\left(n{T}^{+}\right)=(1-{\nu }_M)S_M\left(n{T}^{-}\right),I_M\left(n{T}^{+}\right)=(1-{\nu }_M)l_M\left(n{T}^{-}\right),\hfill \\ S_{Ei}\left(n{T}^{+}\right)=(1-{\nu }_{Ei})S_{Ei}\left(n{T}^{-}\right)\hfill \\ S_{Li}\left(n{T}^{+}\right)=(1-{\nu }_{Li})S_{Li}\left(n{T}^{-}\right),I_{Li}\left(n{T}^{+}\right)=(1-{\nu }_{Li})S_{Li}\left(n{T}^{-}\right),\hfill \\ S_{Ni}\left(n{T}^{+}\right)=(1-{\nu }_{Ni})S_{Ni}\left(n{T}^{-}\right),I_{Ni}\left(n{T}^{+}\right)=(1-{\nu }_{Ni})I_{Ni}\left(n{T}^{-}\right),\hfill \\ S_{Ai}\left(n{T}^{+}\right)=(1-{\nu }_{Ai})S_{Ai}\left(n{T}^{-}\right),I_{Ai}\left(n{T}^{+}\right)=(1-{\nu }_{Ai})I_{Ai}\left(n{T}^{-}\right),\hfill \end{array}\right\}{t}_n=nT,n\in \mathbb{N},$$

(3)

where $t_n$ is the times that prescribed fire is implemented, which may be fixed or non-fixed; in this study we will consider the case with fixed times. The parameters ${\nu }_j,\mathrm{where}j=Ei,Li,Ni,Ai,M$ are the proportion of tick type $i=1,2$ and mice population that is reduced by the fire. In Section 3.1 below, we discuss how these parameters are estimated using data from literature.

3.1. Estimating the Reduction Proportion ${\nu }_j$ Due to Prescribed Fire

To estimate the proportion (${\nu }_{Ei}$, ${\nu }_{Li}$, ${\nu }_{Ni}$, ${\nu }_{Ai},i=1.2$, and ${\nu }_M$) reduced in the ticks and rodent populations due to fire, we used data from [38,67]. The study sites in [38] were located in Western Illinois University’s Alice L. Kibbe Field Station located in Warsaw, Hancock County, Illinois, USA. These sites consisted of several woodlands (oak-hickory, oak barrens), floodplain forests, restored tallgrass prairies and hill prairies. Low intensity burns were carried out at these sites; for low intensity burns, the heights of most flame were less than 1 meter and plant burnt were limited to the understory vegetation [85]. The last time the entire study site was burnt was in 2004 (B04) and two additional burns were later carried out in spring of 2014 (B14) and 2015 (B15).

Following the burns, ticks were collected by sweeping through the vegetation with a $1m^2$ flannel cloth attached to a bamboo stick; this is known as flagging method. This was repeated every two weeks when the vegetation was dry between 12:00 and 18:00 hours from 9 May 2015 to 30 October 2015 and 22 April 2016 to 4 November 2016. To ensure tick mortality, the flannel cloth was placed in a sealed bag and frozen for three days. The ticks were later removed and identified using taxonomic keys, and their DNA extracted for pathogen testing. A total of 23 Dermacentor variabilis, 54 Ixodes scapularis, and 2788 Amblyomma americanum ticks were collected during the two consecutive years 2015 and 2016. A total of 23 D. variabilis were collected, $n=17$ (74%) were adults, $n=2$ (9%) were nymphs, and $n=4$ (17%) were larvae. For I. scapularis, 4% ($n=2$) of the 54 collected, were adults, 22% ($n=12$) were nymphs, and 74% (n = 40) were larvae. In the case of Amblyoma americanum, the ticks collected in B04 made up 51% of the collection ($n=1433$), while those collected in B14 made up 37% ($n=1045$) of the collections, and those collected in B15 constituted 11% $(n=307$) of the collection. Of these ticks, 2% ($n=67$) were adults, 4% ($n=107$) were nymphs, and 93% ($n=2614$) were larvae. Although Ixodes scapularis were collected, we only estimate the parameters quantifying the reduction in the tick population due to fire using the data for Amblyoma americanum and Dermacentor variabilis since they are competent vector for the tuleramia.

Note that the study did not indicate if these were data for the pre-burn or post burn number of I. scapularis collected, nor did it provide data for the number of mice that were caught. So, we use the data for the rodent population provided in [67]. This study have data for six different rodent populations namely California kangaroo rat (Dipodomys californicus californicus), brush mouse (Peromyscus boylii), pin$\tilde{o}$n mouse (P. truei sequoiensis), deer mouse (P. maniculatus gambelii), dusky-footed woodrat (Neotoma fuscipes), and western harvest mouse (Reithrodontomys megalotis longicaudus). After the fire treatment, about half as many rodents were trapped at the treated sites compared with control sites. At the two study sites (MG and DBP) only six tick species were found namely, Ixodes pacificus, Ixodes jellisoni, Ixodes spinipalpis, Ixodes woodi, Dermacentor occidentalis and Dermacentor parumapertus. The data for the Dermacentor subspecies (occidentalis and parumapertus) were not use in estimating the proportions since we are not aware of the transmission efficiency of these vectors unlike Dermacentor andersoni [24] which we do not have data for.

Estimating parameters ${\nu }_{Li},{\nu }_{Ni},{\nu }_{Ai},{\nu }_M,i=1,2$: To estimate these parameters, we divide the numbers collected in each age group by the total number collected and subtract the proportion obtained from 1 to obtain the proportion reduced by fire. For instance, for Dermacentor variabilis, we have ${\nu }_{A1}=1-17/23=0.2609$ for the adult population. For nymphs it is ${\nu }_{N1}=1-2/23=0.9130$, and ${\nu }_{L1}=1-4/23=0.8261$ for larvae. The proportion for Amblyoma americanum are similarly estimated as ${\nu }_{A1}=1-67/2788$, ${\nu }_{N1}=1-107/2788$, and ${\nu }_{L1}=1-2614/2788$. Lastly, the proportion for the rodent population using the data in [67] is given as ${\nu }_M=1-17/23$.

4. Global Sensitivity Analysis

Sensitivity analysis procedure is often used to determine the impact and contribution of the model parameters to model outputs (such as the infected population)[7,25,43]. Results of the sensitivity analysis help to identify the best parameters to target for intervention or for future surveillance data gathering. We implement a global sensitivity analysis using Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC) to assess the impact of parameter uncertainty and the sensitivity of these key model outputs. The LHS method is a stratified sampling technique without replacement this allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter [11,57,58,72]. The PRCC on the other hand measures the strength of the relationship between the parameters and the model outcome, stating the degree of the impact each parameter has on the model outcome [11,57,58,72].

We start by generating the LHS matrices and assuming all the model parameters are uniformly distributed except for the parameters representing number of burns (pulses) and the time between burns ($\tau$). With these parameters we sampled their values from two pseudo-Poisson distributions that exclude zero as in [34,35]. For instance, Fulk et al. [34,35], sampled the values for pulses from a Poisson distribution centered on 10 (since they ran a scenario for 10 years) while $\tau$ was sampled from a Poisson distribution centered on 10. Once the initial distributions are created, “any zeros in either sample were changed to ones to avoid issues with implementing the burns" [34,35]. We then carry out a total of 1,000 simulations (runs) of the model for the LHS matrix, using the parameter values given in Table 2; the minimums and maximums for each parameter is set to $\pm 20\%$ of the baseline values respectively and the reproduction number (${\mathcal{R}}_0$) as the response function for tularemia model (1). We also considered other response functions at the end of the simulation for the tularemia models (1) and (2); these are the infected humans ($I_H$), the infected rodents ($I_M$), the sum of infected Demacenta variablis ticks in all life stages ($I_{L1}+I_{N1}+I_{A1}$), and the sum of infected Amblyomma americanum ticks in all life stages ($I_{L2}+I_{N2}+I_{A2}$).

4.1. Global Sensitivity Analysis for Tularaemia Model (1)

The outcome of the global sensitivity analysis for tularaemia model (1) using the reproduction number (${\mathcal{R}}_0$) is shown in Figure 2. The parameters with significant effect on the reproduction number are those parameters whose sensitivity index have significant p-values less than or equal to $0.05$. The parameters with the most impacts on (${\mathcal{R}}_0$), are rodent birth rate (${\pi }_M$), rodent transmission probability (${\beta }_{TM2}$) to tick 2 (Amblyomma americanum), rodent death (${\mu }_M$) and disease induced mortality (${\delta }_M$) rate, A. americanum carrying capacity ($K_2$), transmission probability (${\beta }_{M2}$) of tularemia to rodents from A. americanum, and the maturation rate (${\tau }_{E2}$) of A. americanum larvae and its adult death rate (${\mu }_{A2})$. Notice that the parameters with significant effects on ${\mathcal{R}}_0$ are related to the rodents and A. americanum.

For the other response functions (see Figure 3), the most significant parameters related to the infected humans, for instance are humans birth and death rates (${\pi }_H,{\mu }_H)$, disease progression rate (${\sigma }_H$), and the recovery rate (${\gamma }_H$). For the infected rodent response function, the significant parameters are rodent transmission probability (${\beta }_{TM2}$) to A. americanum and human disease induced death rate (${\delta }_H$). For the response function related to the sum of infected Demacenta variablis ticks in all life stages ($I_{L1}+I_{N1}+I_{A1}$), the significant parameters are the mortality rate (${\mu }_{A1}$) of adult D. variablis ticks, the maturation rate (${\tau }_{N1}$) of D. variablis nymphs, the mortality (${\mu }_{L1}$) and maturation (${\tau }_{L1}$) rates of larvae, eggs in-viability (${\mu }_{E1}$) and maturation (${\tau }_{E1}$) rates for D. variablis, and transmission probability (${\beta }_{M1}$) of D. variablis to rodents infection; finally for this response function is the rodent disease-induced death rate (${\delta }_M$). Lastly, for the sum of infected A. americanum ticks ($I_{L2}+I_{N2}+I_{A2}$) response function, the significant parameters are the carrying capacity ($K_2$) of A. americanum, the transmission probability $({\beta }_{M2}$) of A. americanum to-rodent, the eggs maturation (${\tau }_{E2}$) rate, the larvae maturation rate (${\tau }_{L2}$), and the nymphs maturation (${\tau }_{N2}$) rate.

The PRCC index values of some of these parameters have positive signs while others have negative signs. The positive signs means that any increase in these parameters will lead to an increase in all the response functions. While the negative signs means that an increase in the parameters will lead to a decrease in the response functions. Hence, these parameters would be useful targets during mitigation efforts. Therefore, control strategies which targets these parameters with significant PRCC values will give the greatest impact on the model response functions.

4.2. Global Sensitivity Analysis for Tularaemia Model (2)

For the sensitivity analysis for tularaemia model (2), we included the number of burns (pulses) and the time between the burns ($\tau$) and drew from Poisson distributions for these parameters. We then used as response functions the sum of infected Demacenta variablis ticks in all life stages ($I_{L1}+I_{N1}+I_{A1}$) and, the sum of infected Amblyomma americanum ticks in all life stages ($I_{L2}+I_{N2}+I_{A2}$).

For the two response functions, the number of burns (pulses) and the time between the burns ($\tau$) were the most significant parameters with negative and positive PRCC values followed by ${\tau }_{L1}$ and ${\beta }_{TM2}$. The sign on the pulses (negative) implies a negative impact on two ticks species. The sign on $\tau$ on the other hand is positive pointing to a positive impact on the ticks. In the next section below we will explore the impact of these two parameters on our model.

Thus, as the frequency of the burn increases few ticks will remain. That is, as the time between the burn increases more ticks will be found.

5. Numerical Simulations: Effect of Fire

Our goal in this research work is to investigate the differential effect of fire on the two tick species and the prevalence of tularemia disease on the two host populations (humans and rodents). To address our research goal, we simulate the tularemia model (2) with prescribed and compare the outcome to the tularemia model (1) without prescribe fire. We focus on infected humans, rodents, and the nymph populations of the two tick species.

In Figure 5 we observed a substantial difference between the number of humans and rodents infected with tularemia when prescribed fire is implemented and when it is absent (see Figure 5a,b). Although we have the same number of infected humans and rodents in the first year before the application of fire but once prescribed fire in implemented we see a reduction in infected humans and rodent population. We observed similar dynamics in both infected Demacenta variablis and Amblyomma americanum ticks but the population of A. americanum rebound quicker than those of D. variablis. Thus showing the differential effect of prescribed fire on the two tick species.

5.1. Differential Effect of Fire on Infected New Cases

Next we consider the differential effect of prescribed fire on cumulative humans and rodents infected new cases. To account for the cumulative number of new cases, we consider for humans the transition to the exposed class ($E_H$) from the susceptible class ($S_H$) due to the individual tick species. Similarly for the rodents, we consider the transition into the infected class ($I_M$) due to the individual ticks (see the force of infection for humans in models (1) and (2)). We then assumed all the ticks parameter values are the same except for the proportion reduced by fire. Figure 6 reveals the outcome of the simulation clearly showing the differential effect of fire. Furthermore, we observed the difference in the cumulative number of new cases due to the individual ticks, with Amblyomma americanum producing more infections than D. variablis, since many of them are remaining after each burn. Observe that in the absence of fire, the number of new cases from the infectious ticks of the two species are the same.

Next we simulate the models to show the differential effect of prescribed fire on the incidence of the disease due to individual tick species. Here the parameter values in Table 2 are used for the simulations and we can further see the differential effect of prescribed fire on the incidence of the disease.

5.2. Frequency and Time between Burn

Next, we consider the time between each burn and explore the effect of the frequency of the burn over ten years. We consider different burning regimes, for instance, we burn annually for a period of ten years, once every other year, once every five years, and once every ten years. We found in Figure 8 more infected humans and rodents as duration between burns increases. Similar increase in the number of infectious nymph populations is observed as the duration between each burn increases. However, the frequently prescribed fire is implemented the fewer the infected populations (humans, rodents, and ticks).

5.3. Burn Environments

Next, we consider the effect of burning in different environments. Using burn data Gleim et al. [39], we consider two additional environments, namely unburned area surrounded with burned area (UBB) and burned area surrounded by unburned area (BUB). This study did not collect data for the rodent population, so we used our baseline data obtained from [41,67] for the rodent population. We then simulate the tularemia model (2) and compare the results with no burn, UBB, BUB, and our baseline burn using data from [41,67].

In Figure 9, we observed almost similar number of infected Demacenta variablis ticks in the UBB and BUB environments. The ticks obtained in these environment are fewer in number than the number obtained in the baseline and unburned environments. For Amblyomma americanum ticks on the other hand, slightly more ticks were present in the UBB environment compared to the BUB, the baseline, and the no burn environments. This dynamics was also observed with the infected humans unlike the rodents that remained the same except for the no burn environment.

6. Discussion and Conclusion

Discussion

In this paper we developed and analyzed a mathematical model for a tick-borne tularemia disease transmission dynamics in humans, rodents and two ticks species (Demacenta variablis and Amblyomma americanum) using the natural history of infection of the disease; the model incorporates prescribed fire. Tularemia is a zoonotic disease caused by Francisella tularensis bacteria [1,44,50,87], a gram-negative coccobacillus-shaped bacterium; There are four subspecies of Francisella tularensis namely tularensis (type A), holarctica (type B),novicida and mediasiatica [50,87]. The most virulent subspecies are tularensis and holarctica [1,32,50,53,62]. The subspecies tularensis and holarctica are the major etiological agents of tularemia in humans [87]. There are multiple transmission routes of the infection to humans such as consumption of contaminated food or water, handling of infected animals or bites from haematophagous arthropod vectors (such as ticks, deer flies, or mosquitoes) [1,24,45]; other routes include contact with aquatic environment, and inhalation of aerosols [8,48,68,87]. In this study we focus on transmission via the bites of ticks.

The goal of the study is to explore the differential effect of fire via prescribed burn on the two ticks species and how that affects the prevalence of the disease in humans and rodents. Prescribed fire or controlled burn is a common and essential land management technique used in ecosystems that depend on fire, such as grasslands, open pine forests, and wetlands [31,83]. Prescribed fire was used often used by the native Americans long before the era of fire suppression. These practises are been revived by several native tribes across the United States [3,29]. These cultural fires help to repair degraded soil and increase biodiversity while reducing the leave litters and overgrowth that contributes to wild fires [3,29,84]. The native plants, which are adapted to fire are not harm by these cultural burns, in fact many of these plants need fire for their seeds to germinate [3].

Quantitative analysis of the tularemia model (1) without fire indicate that the model has a disease-free equilibrium that is locally asymptotically stable when the reproduction number is less than one. The global sensitivity analysis of this model using LHS/PRCC method was carried out to determine the parameter with the most influence on several response functions like the reproduction number (${\mathcal{R}}_0$), the infected humans ($I_H$), the infected rodents ($I_M$), the sum of infected Demacenta variablis ticks in all life stages ($I_{L1}+I_{N1}+I_{A1}$), and the sum of infected Amblyomma americanum ticks in all life stages ($I_{L2}+I_{N2}+I_{A2}$). To implement the analysis, we used parameter values obtained from literature in most cases as well as assumed some where we could not find their values. Identification of these parameters are valuable for targeting interventions and perhaps for gathering data for other analyzes.

The parameters with the most significant impacts on the response functions are rodent birth rate (${\pi }_M$), rodent transmission probability (${\beta }_{TM2}$) to Amblyomma americanum, rodent death (${\mu }_M$) and disease induced mortality (${\delta }_M$) rates, A. americanum carrying capacity ($K_2$), transmission probability (${\beta }_{M2}$) of tularemia to rodents from A. americanum, and the maturation rate (${\tau }_{E2}$) of A. americanum larvae and its adult death rate (${\mu }_{A2})$ (see Figure 2 and Figure 3). Others include humans birth and death rates (${\pi }_H,{\mu }_H)$, disease progression rate (${\sigma }_H$), and the recovery rate (${\gamma }_H$), human disease induced death rate (${\delta }_H$). The significant parameters also include the mortality rate (${\mu }_{A1}$) of adult D. variablis ticks, the maturation rate (${\tau }_{N1}$) of D. variablis nymphs, the mortality (${\mu }_{L1}$) and maturation (${\tau }_{L1}$) rates of larvae, eggs in-viability (${\mu }_{E1}$) and maturation (${\tau }_{E1}$) rates for D. variablis, and transmission probability ( ${\beta }_{M1}$) of D. variablis to rodents infection, rodent disease-induced death rate (${\delta }_M$), the eggs maturation (${\tau }_{E2}$) rate, the larvae maturation rate (${\tau }_{L2}$), and the nymphs maturation (${\tau }_{N2}$) rate.

The sensitivity analysis for tularemia model (2) with prescribed fire was also carried out using as response functions the sum of infected D. variablis ticks ($I_{L1}+I_{N1}+I_{A1}$) and the sum of infected A. americanum ($I_{L2}+I_{N2}+I_{A2}$) ticks at all life stages as the model without fire; the number of burns (pulses) and the time between the burns ($\tau$) were drawn from Poisson distributions. For the two response functions, the number of burns (pulses) and the time between the burns ($\tau$) were the most significant parameters with PRCC values having negative and positive signs meaning negative and positive impact on the two ticks species. Thus, as the frequency of prescribed fire increases fewer ticks will remain or will be found. Furthermore, as the time between the burn increases more ticks will be found since their population would have increased in the time before the implementation of the next burn which will lead to reduction of their population. This result aligns with the outcome obtained by Guo and Agusto [41], and Fulk et al. [34,35].

Arthropods have complex reactions to fire [10]; their responses to fire depend on the intensity and severity of the fire. Fire can lead to abundance or reduction in their population and composition, so much so that only certain arthropods will be able to live in habitats that experience frequent fires especially wildfires [10]. In Figure 6 and Figure 7 we examine the differential effect of fire on the incidence of tularemia due to the two species using data collected by Gleim et al. [39,40]. The data showed more A. americanum than D. variablis were collected after the different burn regimes. This abundance was evident in the cumulative number of new human and rodent infection due to the tick species. To clearly show the differential effect of fire, we assumed the same parameter values for the two species except for the proportion reduced due to fire. As expected A. americanum which was reduced the least produced more infection, see Figure 6. In Figure 7, we relaxed this assumption and simulated the model using the parameters in Table 2; this figure further emphasis the result observed in Figure 7 that A. americanum led to more infection than D. variablis due to its abundance after the burn.

The numerical simulations of the tularaemia model (2) depicted in Figure 5, offer valuable insights into the dynamic interplay between prescribed fire and the spread of tularemia. The observed decrease in the number of infected, particularly in tick species D. variablis and A. americanum, underscores the potential efficacy of prescribed fire in mitigating the prevalence of the disease. This result aligns with the results obtained by Guo and Agusto [41] which showed reduction in Lyme diseased with the use of prescribed fire. Guo and Agusto [41] further showed that high intensity burn reduces more ticks than low intensity burn. This current study further aligns with study by Fulk et al. [35] which showed reduction in the prevalence of Lyme disease in different spatial setting (grassland or wooded area) with the application of prescribed burn. Similar results was also obtained in the study by Fulk and Agusto [34] on the effects of rising temperature and prescribed fire on Amblyomma Americanum with ehrlichiosis. The study found that with significant increase in temperature as high a 2⁰C or 4⁰C, implementing prescribed burning can still lead to reduction in the prevalence of infectious questing nymphs and adults with Ehrlichiosis.

Using the results from the sensitivity analysis we explore in Figure 8 the frequency and time between burns. The parameters used for the proportion reduced as a result of fire was obtained from [41,67]. Considering the results from Figure 8, there is a considerable reduction in infected humans and rodents as well as the infectious nymphs generally as a result of the burn. However, the frequency of the burn plays a significant role in ensuring whether there is a significant reduction or not. As can be deduced from the figure, one burn every year ensures a significant reduction in the infected (humans, rodents, and ticks).The simulation results emphasize the importance of well-timed fire applications in controlling tularemia. This result confirms the results from Guo and Agusto [41] which explore the effect of fire frequency and the time between burns. As identified in the study, the duration between burns has a more significant effect on ticks with higher-intensity fires than with lower-intensity fires. However, fire intensity appears to have a larger influence on tick reduction than the duration of the burns, as burning fewer times at a higher intensity is more effective than burning more times at a lower intensity.

The investigation of the different burn environments in Figure 9 reveals intriguing infection patterns in different burn environments (UBB and BUB). The number of infected was almost similar especially for D. variabilis and A. Americanum. Gleim et al. [39,40] indicted that frequent and long-term prescribed fire significantly reduces tick abundance at sites with varying burn regimes such as UBB and BUB. The results of this study suggests that the spatial arrangement of burned and unburned areas does not matter on the overall impact of prescribed fire on ticks and tularemia dynamics. This result may be affected by the fact that the same proportion of rodent population was used in the simulation of both environments. Considering the study done by Gleim et al. [39], a slightly different number of nymphs and adults of both species were obtained in both environment. Gleim et al. [39] research showed the differential effect of fire on these tick species; A. americanum is more fire tolerant than D. variablis; more A. americanum larvae were found in the two environments over none found for D. variablis. Our simulation results also shows this differential effect of fire on the two tick species, fewer nymphs of D. variablis remained at the end of our simulation period compare to the A. americanum remained.

6.1. Conclusion

To conclude, in this study we developed a deterministic model of ordinary and impulsive differential equations to gain insight in the differential effect of prescribed fire on Demacenta variablis and Amblyomma americanum ticks and the prevalence of tularemina. We found that prescribed fire can differentially reduce the number of the two ticks with D. variablis been reduced the most compare to A. americanum. We summarize the other results as follows:

(i)

The results of the sensitivity analyzes using as response or output functions the reproduction number (${\mathcal{R}}_0$), the infected humans ($I_H$), the infected rodents ($I_M$), the sum of infected D. variablis ticks at all life stages ($I_{L1}+I_{N1}+I_{A1}$), and the sum of infected A. americanum ticks at all life stages ($I_{L2}+I_{N2}+I_{A2}$) to identify the parameters with the most impact on these functions in no particular order are: the most significant parameters related infected humans are the birth, death, and disease induced death rates (${\pi }_H,{\mu }_H$, ${\delta }_H$), human disease progression rate (${\sigma }_H$), and the recovery rate (${\gamma }_H$); the most significant parameters related to infected rodents are rodent birth, death, and disease induced death rates (${\pi }_M$, ${\mu }_M$, ${\delta }_M$), rodent transmission probability (${\beta }_{TM2}$) to A. americanum. The significant parameters related D. variablis are eggs maturation (${\tau }_{E1}$) and in-viability (${\mu }_{E1}$) rates, the larvae maturation (${\tau }_{L1}$) and mortality (${\mu }_{L1}$) rates, the nymphs maturation rate (${\tau }_{N1}$), and the adult mortality rate (${\mu }_{A1}$). The significant parameters related to A. americanum are its carrying capacity ($K_2$), the eggs maturation (${\tau }_{E2}$) rate, the larvae maturation rate (${\tau }_{L2}$), and the nymphs maturation (${\tau }_{N2}$) rate, its adult death rate (${\mu }_{A2})$, and the transmission probability (${\beta }_{M2}$) to rodents.

(ii)

Prescribed fire can reduce the number of ticks leading to a reduction in the number of tularemia infected humans, rodents and ticks.

(iii)

A. americanum produced more tulareamia infection in humans and rodents than D. variablis due to the differential effect of prescribed fire which leaves more A. americanum than D. variablis after a burn.

(iv)

As time between burn increases, more infected humans, rodents, and ticks increases. Frequent burning reduces the number of ticks and therefore infections.

(v)

The spatial arrangement of UBB and BUB areas may not matter as either arrangement led to fewer ticks and reduction in tularemia transmission with the implementation of prescribed fire.

In this study, we have used mathematical models to investigate the differential effect of prescribed burn on two tick species namely D. variablis and A. americanum and explored the effectiveness of prescribed fire in managing tick populations and lowering the prevalence of tularemia disease in the population (humans, rodents, and ticks). Our study is not without some limitations; the biggest drawbacks was the availability of data. For instance, there are no rodent data for the UBB and BUB environment. Despite this drawback, we have used the limited available data and we understand it might affect our results but we have confidence in our conclusions given the results from the sensitivity analysis we carried out.

Occhibove et al. [66] found that the degree of vector generalism affected pathogen transmission with different dilution outcomes. A. americanum is known as an aggressive and generalist hematophagous tick species [33,52,69]. In future studies, we will explore the differential effect of prescribed burn on the dilution or amplification tendencies of this tick species on the transmission dynamics of tick-borne diseases. Furthermore, we would also explore other ecological implications of our findings, by considering factors such as habitat fragmentation and species interactions in the presence of prescribed fire.