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
C and
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
exposed
asymptomatic
infected
and recovered
The total human population is denoted by
and is defined as
The total rodent population is denoted by
and is defined as
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 (
and
) and (
and
), nymphs (
and
) and (
and
), and adults (
and
) and
and
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 (
and
).
For simplicity, we assume that human, tick and rodent populations are mixing homogeneously. The susceptible human sub-population are recruited at a rate
The rate at which a susceptible human progress to the exposed class is denoted by
and is defined as
This quantity
is the force of infection for the humans, where
and
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
or infected at a rate
with
. 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
according to the reciprocal of the length of the treatment days. The human population is decreased by natural death at the rate denoted by
and by disease induced mortality denoted by
. 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
. The rodents become infected following a bite from infected ticks of type 1 and 2 at the rate
where
and
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
and by natural death at the rate by
.
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
and
, with a certain proportion dying naturally at rates
and
for tick 1 and tick 2, respectively. The susceptible larvae of both tick 1 and tick 2 become infected larvae at a rate of
and
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
where
and
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
for tick 1 and
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
and
Following a subsequent infectious blood meal, susceptible nymphs in tick 1 and tick 2 become infected nymph at a rate of
and
respectively. The populations of both susceptible and infected nymphs from both tick species are reduced due to natural death rate of
and
and they mature into adults at the rate
and
respectively.
Given the assumptions above, the following nonlinear equations are given for the transmission dynamics of tularemia in two tick species:
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
.
Lemma 1.
Let the initial data , where , where . Then the solutions of the tularemia model (1) are non-negative for all . Furthermore
Invariant Regions
The tularemia model (
1) will be analyzed in a biologically-feasible region as follows. Consider the feasible region
with,
and
Lemma 2. The region is positively-invariant for the model (1) with non-negative initial conditions in .
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
. The equilibrium
, involves humans and rodents only; while
involves, humans, rodents and
Demacenta variablis only; the equilirium
involves humans, rodents and
Amblyomma americanum only. Lastly, the
involves humans, rodents and the two tick species. The expression for these equiliria are given as
where
We will focus on the the stability of
, the stability of the other equiliria
can be determined using the same approach. Thus, the stability of
can be established using the next generation operator method on system (
1). Taking
, and
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:
where,
. Therefore, using the definition of the reproduction number
from [
82] we have
where
is the spectral radius and
The expressions
and
are the average number of secondary infections in rodents due to infectious
Demacenta variablis and
Amblyomma americanum while
and
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 and unstable if .
The quantity
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 (
) can be brought to (and maintained at) a value less than unity.
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 (
). 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
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
of the baseline values respectively and the reproduction number (
) 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 (
), the infected rodents (
), the sum of infected
Demacenta variablis ticks in all life stages (
), and the sum of infected
Amblyomma americanum ticks in all life stages (
).
4.1. Global Sensitivity Analysis for Tularaemia Model (1)
The outcome of the global sensitivity analysis for tularaemia model (
1) using the reproduction number (
) 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
. The parameters with the most impacts on (
), are rodent birth rate (
), rodent transmission probability (
) to tick 2 (
Amblyomma americanum), rodent death (
) and disease induced mortality (
) rate,
A. americanum carrying capacity (
), transmission probability (
) of tularemia to rodents from
A. americanum, and the maturation rate (
) of
A. americanum larvae and its adult death rate (
. Notice that the parameters with significant effects on
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 (
, disease progression rate (
), and the recovery rate (
). For the infected rodent response function, the significant parameters are rodent transmission probability (
) to
A. americanum and human disease induced death rate (
). For the response function related to the sum of infected
Demacenta variablis ticks in all life stages (
), the significant parameters are the mortality rate (
) of adult
D. variablis ticks, the maturation rate (
) of
D. variablis nymphs, the mortality (
) and maturation (
) rates of larvae, eggs in-viability (
) and maturation (
) rates for
D. variablis, and transmission probability (
) of
D. variablis to rodents infection; finally for this response function is the rodent disease-induced death rate (
). Lastly, for the sum of infected
A. americanum ticks (
) response function, the significant parameters are the carrying capacity (
) of
A. americanum, the transmission probability
) of
A. americanum to-rodent, the eggs maturation (
) rate, the larvae maturation rate (
), and the nymphs maturation (
) 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 (
) 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 (
) and, the sum of infected
Amblyomma americanum ticks in all life stages (
).
For the two response functions, the number of burns (pulses) and the time between the burns () were the most significant parameters with negative and positive PRCC values followed by and . The sign on the pulses (negative) implies a negative impact on two ticks species. The sign on 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.