Gut mucosal microbiota profiles linked to development of positional-specific human colorectal cancer

1.   Introduction

Ehrlichiosis, formally known as human monocytic ehrlichiosis, was recognized as a human disease in the late 1980s; however, official reporting by the CDC on the disease began in 2000 ^[1,2]. Canine ehrlichiosis, a disease caused by the Ehrlichia bacteria, has received particular attention due to its influence on the health of domesticated dogs. These illnesses, which are transmitted to people by tick bites, pose a substantial danger to world health due to their potential for broad transmission and complex ecological dynamics. Ehrlichia ewingii, Ehrlichia muris eauclairensis, and Ehrlichia chaffeensis are bacterial species belonging to the class Ehrlichiosis. Among these, Ehrlichia chaffeensis is the most prevalent in the United States ^[2]. The bacterium responsible for ehrlichiosis is transmitted to humans primarily by Amblyomma americanum, which is of significant concern in the United States not only due to its extremely aggressive and somewhat indiscriminate biting behavior throughout all three life stages, but also because of its potential to transmit several bacteria known to cause diseases in humans ^[3].

Amblyomma americanum, commonly known as the lone star tick, is widely distributed across the eastern and central United States, with a range that has expanded due to climate change and habitat alterations ^[4,5]. It thrives in wooded areas, grasslands, and dense brush, where it can easily find hosts. The tick's life cycle consists of four stages, egg, larva, nymph, and adult, each requiring a blood meal to progress to the next stage ^[6]. A. americanum feeds on a variety of animals, including mammals, birds, and reptiles. Larvae and nymphs feed on small mammals such as rodents and rabbits ^[7], while adult ticks feed on large mammals such as foxes and coyotes but they primarily rely on white-tailed deer as their main hosts ^[8]. Coyotes, foxes, and owl are predatory animals that feed on small animals such as rabbits and rodents ^[9].

Extensive research has been conducted on the role of host abundance and tick ecology in disease transmission (see ^{[10,11,12,13]} and references therein); however, the impact of predation — both lethal (direct consumption of hosts) and non-lethal (behavioral and physiological changes in host populations due to predation risk) — remains less explored in mathematical models. Fear of predation can lead to a range of behavioral, physiological, and ecological responses in prey populations. The effect of these responses on the prey population and disease dynamics can be as significant as direct predation itself ^[14]. For instance, the fear of predators can alter host behavior, leading to reduced foraging and movement, which may decrease tick-host contact rates ^[15,16]. Chronic stress from predation risk can suppress host immune responses, potentially increasing their susceptibility to infection ^[17].

Changes in habitat use due to predator presence may modify tick-host encounter rates, either reducing or amplifying disease transmission ^[16]. Prey-predator interactions are fundamental drivers of ecosystem dynamics, with cascading effects that can ripple through entire food webs ^[18,19]. Wang et al. ^[20] explore the fear effect in prey-predator interactions, revealing how the perception of predation risk can influence prey behavior and, by extension, disease transmission dynamics. This insight broadens our understanding of the non-lethal effects of predation.

Some recent studies have shed light on how the non-lethal effects of predation affect disease dynamics. In this context, Hossain et al. ^[21] observed that fear can eliminate chaotic oscillations of their proposed system by either making the endemic equilibrium regular (stable equilibium point or stable limit cycle) or moving it towards the disease-free state. Recently, Zhang ^[22] investigated the effect of the fear factor on an eco-epidemiological system with a standard incidence rate and found that an increase in the level of fear changed a disease-free equilibrium from unstable to stable. Moreover, the study by Pribadi and Chalasani ^[23] delves into fear conditioning in invertebrates, opening up new avenues for understanding the psychological and behavioral aspects of prey species when confronted with predation threats. However, the intricate relationship between prey-predator interactions, fear, and tick-borne diseases has not been comprehensively explored. This gap in our knowledge highlights the critical need for further investigation on how ecological factors like prey-predator interactions and fear responses influence disease dynamics. To effectively mitigate tick-borne diseases, we must understand how prey-predator interactions, mediated by fear, impact disease prevalence in ticks.

The goal of this study is to understand the complex interplay between lethal and non-lethal predation on the dynamics of tick-borne diseases with the aim of addressing the following specific research questions:

1. What are the cascading effects of fear within ecosystems, and how do they influence tick abundance, prey behavior, and ultimately, disease transmission dynamics?

2. Can fear-induced changes in host behavior and predator dynamics serve as effective tools in managing tick-borne disease risk within a given ecosystem?

The rest of the paper is organized as follows: in Section 2, we begin by formulating the ticks-prey-predator dynamics with disease.

In order to capture the intricacies of Amblyomma Americanum ticks' life cycle, we classify the ticks into their different life stages and explore stability analysis and sensitivity analyzes of the resulting model. Moving into Section 3, we use a reduced version of the disease model following some stated assumption and subsequently carry out the stability. In Section 5, we carry out extensive simulations using results from the sensitivity analysis. These simulations offer a glimpse into the practical implications of our findings, bridging the gap between theoretical insights and practical applications. We discuss and summarize our findings in Section 6.

2.   Model Formulation

To formulate the ticks-prey-predator dynamics with disease, we adapt some of the approaches used in earlier modeling of ticks dynamics in ^[24,25,26]. We assume that transovarial transmission is very unlikely for ehrlichiosis, so all eggs laid are without infection. The eggs hatch into susceptible larvae ($L_{S}$) which progresses into infectious larvae ($L_{I}$) after acquiring infection from feeding on an infected prey ($M_I$) at the rate $\lambda_{LM}$ expressed as

where $\beta_{LM}$ is the effective transmission probability and $M$ is the total prey population (susceptible and infected). The susceptible and infected larvae mature at the rate $\sigma_L$ or die in a density dependent manner $(\mu_L (L_S, L_I, M)$) express as

Thus, we have the following equations for susceptible and infected larvae:

Regardless of whether or not larvae are susceptible or infectious ($L_{S}$ and $L_{I}$, respectively), they molt into nymphs at the rate $\sigma_L$ once fully engorged from a blood meal, or die at the density dependent rate, $\mu_N (N_S, N_I, M, P)$ given as

The nymph life-stage is the least discriminant when it comes to host preference and as such, we incorporate the potential for nymphs to feed on both prey and predator which are small and large hosts respectively. At this stage the susceptible nymphs may become infected if feeding on an infectious host. The force of infection, $\lambda_{NM}$, is defined as:

where $\beta_{NM}$ is the probability that a susceptible nymph becomes infected after feeding on an infectious host. $M$, and $P$ are the number of susceptible and infectious prey and predator hosts. The susceptible and infectious nymphs mature to become adult ticks at the rates $\sigma_N$. The equations describing the rates of change for susceptible and infectious nymph populations ($N_{S}$ and $N_{I}$, respectively) are:

The susceptible and infectious adult ticks quest for hosts, especially larger hosts like the predator population. At this point, there is the possibility of a susceptible tick becoming infected if feeding on an infectious host. The force of infection for feeding adult ticks is:

where $\beta_{AP}$ is the probability that a susceptible adult becomes infected after feeding on an infectious host. The susceptible and infectious populations are reduced by natural death at the rates ($\mu_A (A_S, A_I, P) A_S$), expressed as

With these, we have the following equations for susceptible and infectious feeding adults, given as $A_{S}$ and $A_{I}$:

Note that once a larvae (or nymph or adult) becomes infectious, it remains so for the rest of its life.

Like the expression for the susceptible and infected tick populations we can express the susceptible and infected in the prey and predator systems as

where the forces of infection for the prey and predator $\lambda_{MT}$ and $\lambda_{PA}$ are express as:

where $\beta_{MT}$ and $\beta_{PA}$ are the probability that a susceptible prey and predator becomes infected after been fed on by an infectious tick.

Furthermore, we assume that the host do not recover from the infection, see ^[27,28]. With these assumptions we have the following system of equations for the eggs and the larvae, nymphs, and adults:

The model flow diagram is given in Figure 1.

2.1.   Analysis of the ticks-prey-predator model (2.1)

2.1.1.   Positivity and boundedness of solutions

Following the results in subsection 3.1 and the approach in Guo and Agusto ^[29] we can prove the following lemmas for model positivity and boundedness of solutions.

Lemma 1. Let the initial data $F(0) \geq 0$, where $F(t) = (E(t), L_S(t), L_I(t), N_S(t), N_I(t),$ $A_S(t), A_I(t), M_S(t), M_I(t), P_S(t), P_I(t)$). Then the solutions $F(t)$ of the system are non-negative for all $t > 0$.

Lemma 2. The region $\Omega = \Omega_T \cup \Omega_{MP} \subset \mathbb{R}^{7}_+\times \mathbb{R}^{4}_+$ is positively-invariant for the system with non-negative initial conditions in $\mathbb{R}^{11}_+$.

where,

$N_T = L_S(t)+L_I(t)+N_S(t)+N_I(t)+A_S(t)+ A_I$, and $E(t) < K_E$.

where $r = \dfrac{K_M(\pi_M-\mu_M+\eta)^2}{4 \pi_M} > 0$.

Proof. The first seven equations of model (2.1) give the following after summing the equations representing the larvae, nymphs, and adult stages

where $\mu_T = \min \{\mu_L(L_{S}, L_{I}, M), \mu_N(N_{S}, N_{I}, M, P), \mu_A(A_{S}, A_{I}, P)\},$ Since $K_E$ is the carrying capacity, it follows that $E\leq K_E$. Hence, equation (2.3) becomes

Thus,

If $N_T(0) = \dfrac{\sigma_EK_E}{\mu_T}$, then $N_T(t) = \dfrac{\sigma_E K_E}{\mu_T}$. Hence, equation (2.4) implies that $N_T(t)$ is bounded. Furthermore, in Section 3.1 below, we show that the domain $\Omega_{MP}$ is bounded. Thus, all solutions starting in the region $\Omega$ remain in $\Omega$. Hence, the region is positively-invariant and hence, the region $\Omega$ attracts all solutions in $\mathbb{R}^{11}_+$. □

2.1.2.   Stability of disease-free equilibrium (DFE)

The tick-prey-predator model has a disease-free equilibrium (DFE) denoted by $\mathcal{E}_0$. The DFE is obtained by setting the right-hand sides of the equations in the model (2.1) to zero, which is given by

The analytical expression of the DFE is difficult to obtain due to the nonlinearity of the density-dependent fecundity and mortality functions in the tick population, see subseection 3.1.2 above. To establish the stability of disease-free equilibrium ($\mathcal{E}_0$) we use the reproduction number $\mathcal{R}_0$ which is obtained using the next generation operator method on system (2.1). Taking $L_I, N_I, A_I, M_I$, and $P_I$ as the infected compartments and then using the notation in ^[30], the Jacobian $F$ and $V$ matrices for new infectious terms and the remaining transfer terms are defined as:

where $k_1 = \dfrac{\alpha_MP^*}{a+M^*} + \mu_M, \; k_2 = \sigma_L + 0.25 + 0.015\; {\rm ln}\bigg(\dfrac{1.01}{M^*_S}\bigg), \; k_3 = \sigma_N + 0.45 + 0.02 \; {\rm ln}\bigg(\dfrac{1.01}{M^*_S + P^*_S}\bigg),$ and $\mu_{A_I}(P^*) = 0.35 + 0.03\; {\rm ln}\bigg(\dfrac{1.01}{P^*_S}\bigg)$.

Therefore, the reproduction is given as ${\mathcal{R}}_0 = \rho({F}{V}^{-1})$, where $\rho$ is the spectra radius is given as

It is worth mentioning that $\mathcal{R}_0 > 0$ provided

Thus, the following result is established using Theorem 2 in ^[30].

Lemma 3. The DFE of ticks-prey-predator model (2.1) with disease, given by ${{\cal E}}_{0}$, is locally asymptotically stable (LAS) if ${{\cal R}}_0 < 1$, and unstable if ${{\cal R}}_0 > 1$.

The basic reproduction number (${{\cal R}}_0$) measures the average number of new infections generated by a single infected individual (tick, prey or predator) in a completely susceptible population ^{[30,31,32,33]}. Thus, Lemma 3 implies that malaria can be eliminated from the human population (when ${{\cal R}}_0 < 1)$ if the initial sizes of the sub-populations are in the basin of attraction of the DFE, ${{\cal E}}_{0}$.

3.   Ticks-prey-predator model without disease

To formulate the ticks-prey-predator dynamics, we classified the Amblyomma Americanum ticks into the following life stages: egg $E$, larvae $L$, nymph $N$, and adult $A$ while $P$ and $M$ are predator and prey populations respectively.

The number of eggs are determined, in part, by the number of egg-laying adults. The fecundity of egg-laying adults is reduced in a density-dependent manner by the factor $\pi_E (A, P)$, with a female having the capacity to lay a maximum of $f$ = 6000 eggs if $\pi_E (A, P) = 1$. The function is given as:

where $P$ is the population of the predator. The egg hatch into larvae at the rate $\sigma_E$ or the decay at a fixed daily, per-capita rate ($\mu_E$). The equation for the eggs ($E$) is given by:

The larvae population ($L$) increases depending on the egg developmental rate or hatching rate ($\sigma_E$). The population decreases via maturation by a fixed rate $\sigma_L$ as well as via death in a density dependent manner described as:

Thus, the equation for the larvae ($L$) population is given by:

The population of the nymphs increase following the maturation of the larvae. Nymphs that survive feeding on a host then become engorged, once engorged, they either molt into adult ticks at the rate $\sigma_N$, or die in a density-dependent manner described as:

where $M$ is population of the prey. Hence, the equations describing the rates of change for the nymph population ($N$) is given as:

Once attached to a host, the adult ticks can take up to several days to finish feeding. They may also die in a density dependent manner:

All of these components lead to the following equations for adult ticks, given as:

We assume that the prey population grew logistically; however, the population is also limited due to the fear of the predator at the rate $\bigg(\dfrac{\pi_M}{1+\tau P}\bigg)$, where $\pi_M$ is the preys birth rate and $\tau$ regulates the strength of fear incited by the predator. This decreasing fear function $\left(\frac{1}{1+\tau P}\right)$ was first introduced by Wang et al. ^[20]. The function has gained widespread acceptance within the research community, as evidenced by its adoption and application in various studies. For more information, see ^[21,34,35] and references therein. Biologically meaningful properties are embedded within the fear function, which dictate its behavior within the model framework.

$({\rm{i}})$ In the absence of fear ($\tau = 0$), the fear function yields:

$({\rm{ii}})$ When there are no predators ($P = 0$), the fear function becomes:

$({\rm{iii}})$ As the fear level increases ($\tau\rightarrow \infty$), the fear function approaches zero:

$({\rm{iv}})$ When the predator population is abundant ($P\rightarrow \infty$), the fear function tends to zero, causing the birthrate or disease transmission rate to decline:

The properties of the fear function represent biological and epidemiological factors such as fear intensity, predator presence, and population dynamics, all of which impact disease transmission.

The prey population is reduced by predation by a Holling type 2 functional response at the rate $\bigg(\dfrac{\alpha_M M}{a+M}\bigg)$, where $\alpha_M$ is the attack rate of the predator. We note that $\mu_M$ represents the death rate of the prey that is not influenced by fear of predators and thus is isolated from the logistic equation for the prey population. This novel phenomenon of splitting the natural death rate into two in prey dynamics involving fear of predation was recently introduced by Antwi-Fordjour et al. ^[36]. Thus, the rate of change of the prey population is expressed as:

The predator population increases after predation of the prey population at the rate $\bigg(\dfrac{\pi_P M}{a+M}\bigg)$, where $\pi_P$ is the energy conversion rate due to predation of the prey. The energy conversion rate from prey to predator measures how efficiently a predator converts the energy from its prey into its own growth and reproduction. The predator population is reduced due to mortality at the rate $\mu_P$. Thus the equation of the predator population is given as:

Combining the above description and assumptions we have the following system of equations for the ticks-prey-predator system:

The model flow diagram is given in Figure 2 and the description of the parameters and the values used in simulations are given in Table 1.

3.1.   Quantitative analysis of ticks-prey-predator model (3.1)

3.1.1.   Positivity and boundedness of solutions

We shall investigate the positivity and boundedness of the prey and predator populations of system (3.1). From an ecological standpoint, the biomass of a population cannot exhibit negativity at any point in time, thereby necessitating the persistence of non-negative solutions. Boundedness guarantees that these populations will not undergo exponential growth over an extended duration.

Theorem 1. For all $t\geq 0$, the solutions $(M(t), P(t))$ of the system (3.1) with initial conditions $M(0)\geq 0$ and $P(0)\geq 0$ remain non-negative.

Proof. Since the functions on the right-hand side of system (3.1) are continuous and locally Lipschitz in the positive octant $\mathbb{R}^2_+$, the system has a unique solution that exists for all $t\in [0, \infty)$.

The predator ($P$) and prey ($M$) equations of system (3.1) can be rewritten as

where

We apply the integrating factor technique to system (3.2), a fundamental method for solving differential equations, to find the solutions of the system (3.1);

From the expressions above, it is clear that $M(t)$ and $P(t)$ remain non-negative for all future time if they initiate from an interior initial point of

□

Theorem 2. The solutions $(M(t), P(t))$ of system (3.1) with initial conditions $M(0) > 0, P(0) > 0$ are ultimately bounded.

Proof. Consider the function $W(M(t), P(t)) = M(t)+P(t)$. Thus

Assume $\pi_P < \alpha_M$, and $\eta\in \mathbb{R}_+$ such that $\eta \leq \mu_P$ and $\pi_M+\eta > \mu_M$. Due to the positivity of the solutions, it is obvious that $M(t)\leq K_M$ by comparison to the logistic equation.

Then,

where $r = \dfrac{K_M(\pi_M-\mu_M+\eta)^2}{4 \pi_M} > 0$. By solving the differential inequality above, we obtain the solution

Thus,

Therefore, the solutions $(M(t), P(t))$ to the system (3.1) with positive initial conditions are ultimately bounded. □

3.1.2.   Stability analysis of the equilibrium points

To find equilibrium points, we need to set the derivatives of all variables to zero and solve the resulting equations simultaneously.

Denote the equilibrium values as follows:

Since the equations are highly nonlinear in nature due to the density-dependent fecundity and mortality in the tick population, making it difficult to find analytical solutions. As a result, it may not be feasible to solve them explicitly. Instead, it may be necessary to employ numerical methods or utilize software tools to obtain approximate values for the equilibrium points.

However, it is biologically feasible to assume that the number of larvae, nymphs, and adult ticks are smaller than the number of eggs (that is, $A\! < N\! < L < \!K_E$), since most ticks die if they are unable to find a host to fed on ^[40]; from subsection 3.1 above we have that the prey and predator population are bounded, that is $|M|\leq C_1$ and $|P|\leq C_2$, where $C_1$ and $C_2$ are non-negative constants; then we can assume that the tick fecundity and mortality are constant rates (that is $\pi_E, \; \mu_E, \; \mu_L, \; \mu_N$, and $\mu_A$). Thus model (3.1) reduces to the following simple system

where $\pi_T = f\pi_E$, and $T$ denotes the total tick population from eggs to adult. The reduced system (3.3) has three non-negative steady states namely:

$({\rm{i}})$ Tick-only equilibrium $(T, M, P) = \left(\dfrac{\pi_T}{\mu_T}, 0, 0\right)$.

$({\rm{ii}})$ Predator-free equilibrium $(T, M, P) = \left(\dfrac{\pi_T}{\mu_T}, K_M \left(1-\frac{\mu _M}{\pi _M}\right), 0 \right)$ for $\pi_M > \mu_M$.

$({\rm{iii}})$ The coexistence equilibrium $(T^*, M^*, P^*) = \left(\dfrac{\pi_T}{\mu_T}, \dfrac{a \mu _P}{\pi _P-\mu _P}, \dfrac{-w_1\mp\sqrt{w_2}}{2 \alpha _M \sqrt{K_M} \tau \left(\pi _P-\mu _P\right)}\right)$, where $\pi_P > \mu_P, \; w_2 > w_1^2$, and

Thus, the following Theorem gives the conditions for the stability of equilibrium points obtained from system (3.3).

Theorem 3. Consider the system given by system (3.3). Then the following statements hold:

(i) The tick-only equilibrium is locally asymptotically stable if $\mu_M > \pi_M$,

(ii) The predator-free equilibrium point is locally asymptotically stable if $\mu_P > \frac{K_M \pi _P \left(\mu _M-\pi _M\right)}{K_M \mu _M-\pi _M \left(a+K_M\right)}$,

(iii) The coexistence equilibrium is locally asymptotically stable provided the signs of $w_3$ and $w_4$ given in (3.6) are positive.

Proof. The Jacobain matrix for the system (3.3) at any equilibrium point is given as

$({\rm{i}})$ To analyze the stability of the equilibrium where only ticks are present, we begin by evaluating the Jacobian matrix at the tick-only equilibrium point.

The eigenvalues associated with the Jacobian matrix $J^*$ are $\lambda_1 = -\mu_T < 0, \; \lambda_2 = \pi_M-\mu_M$ and $\lambda_3 = -\mu_P < 0$. Thus the tick-only equilibrium is locally asymptotically stable if $\mu_M > \pi_M$.

$({\rm{ii}})$ To investigate the stability of predator-free equilibrium, we evaluate the Jacobian matrix at the predator-free equilibrium point,

The eigenvalues associated with the Jacobian matrix $J^{**}$ are $\lambda_1 = -\mu_T < 0, \; \lambda_2 = -\pi_M+\mu_M < 0$ and $\lambda_3 = \frac{K_M \pi _P \left(\mu _M-\pi _M\right)}{K_M \mu _M-\pi _M \left(a+K_M\right)}-\mu _P$. Thus the predator-free equilibrium is locally asymptotically stable if $\mu _P > \frac{K_M \pi _P \left(\mu _M-\pi _M\right)}{K_M \mu _M-\pi _M \left(a+K_M\right)}$,

$({\rm{iii}})$ From the Jacobian matrix (3.4), we derive the following characteristic equation

where

note that $w_3 > 0$, if

and

Additionally, $w_4 > 0$ if

The eigenvalues associated with the equation in (3.6) are $\lambda_1 = -\mu_T < 0$ and other two roots (i.e., $\lambda_2$ and $\lambda_3$) of the quadratic equation in $\lambda^2 + w_3\lambda + w_4 = 0$. The other roots are given as $\lambda_2+\lambda_3 = w_3 > 0$ and $\lambda_2\lambda_3 = w_4 > 0$. Thus the stability of the coexistence equilibrium point depend on sign of $w_3$ and $w_4$.

□

4.   Sensitivity analysis for models (2.1) and (3.1)

Deterministic model outputs are governed by its input parameters, which exhibit some uncertainty in the process of their selection. Exploration of the ticks-prey-predator models (2.1) and (3.1) above requires determination of which parameters most strongly affect key model outputs (such as the number of infected individuals). A global sensitivity analysis was carried out 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 of the numerical simulations to variations in each parameter. LHS is a stratified sampling without replacement technique which allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter ^{[41,42,43,44]}. To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. PRCC measures the strength of the relationship between the model outcome and the parameters, stating the degree of the effect that each parameter has on the outcome ^{[41,42,43,44]}. It is a reliable sensitivity measure for nonlinear yet monotonic relationships between the model input and output, provided that the inputs have minimal or no correlation ^[42].

We start by exploring first the sensitivity analysis of model(3.1) without disease and then model (2.1) with disease.

4.1.   Sensitivity analysis of model (3.1) without disease

For the ticks-prey-predator model (3.1), a total of 1,000 runs were then carried out for the model using parameter values sampled from the LHS matrix and the outcome were ranked using the PRCC method. The parameter values used are given in Tables 1 (with ranges varying from $\pm 20\%$ the stated baseline values). The response functions in this case are the sum of the ticks $L+N+A$, the preys ($M$), and the predators ($P$).

The outcomes of the global sensitivity analysis are shown in Figure 3 and given in Table 2.

The parameters with the most significant impact on the sum of ticks ($L+N+A$) are tick eggs maturation rate ($\sigma_E$), the egg carrying capacity ($K_E$), the predator birth and natural death rates ($\pi_P$, and $\mu_P$), the predator's attack rate ($\alpha_M$), and the half saturation constant ($a$). Details are shown in Figure 3(a).

The significant parameters for prey ($M$) response function are shown in Figure 3(b); these parameters are the carrying capacity ($K_M$) of the prey, the predator birth and natural death rates ($\pi_P$, and $\mu_P$), the predators attack rate ($\alpha_M$), and the half saturation constant ($a$).

For the predators ($P$) response function, the significant parameters are the predators birth and natural death rate ($\pi_P$ and $\mu_P$). Other parameters are the preys' birth and natural death rate ($\pi_M$ and $\mu_M$) and their carrying capacity ($K_M$). The predators attack rate ($\alpha_M$), the half saturation constant ($a$), and the fear of predators ($\tau$) are equally significant, see Figure 3(c).

Note that the fear factor $\tau$ have a variable sensitivity indices for the three response function. It has a negative impact on the sum of the infected ticks, and the predator, and a positive significant impact on the prey population; however, the magnitude of the index is small.

The sensitivity indices are either positive or negative indicating the direction of the impact on the response functions, and the values of the indices indicate the magnitude of the response functions, see Table 2.

4.2.   Sensitivity analysis of model (2.1) with disease

For ticks-prey-predator model (2.1) a total of 1,000 simulations (runs) of the model for the LHS matrix were also carried out, using the parameter values given in Tables 1 (with ranges also varying from $\pm 20\%$ the stated baseline values) and the same response functions namely the sum of infected ticks $L_I+N_I+A_I$, infected preys ($M_I$), and infected predators ($P_I$). The PRCC method was then used to rank the parameter after the simulation runs were completed. The outcomes of the global sensitivity analysis are shown in Figure 4 and given in Table 3.

Figure 4(a) depict the parameters with the most significant impact on the sum of infected ticks ($L_I+N_I+A_I$), these are tick eggs maturation rate ($\sigma_E$), the egg carrying capacity ($K_E$), the predators birth rate ($\pi_P$), the predators attack rate ($\alpha_M$), and the half saturation constant ($a$).

The impact of the model parameters on the infected preys ($M_I$) are shown in Figure 4(c). The impact of the parameters on the response function $M_I$ are small, however, some of the effect of these parameters significant on $M_I$. The significant parameters are the nymphs maturation rate ($\sigma_N$), the predator birth rates ($\pi_P$), and the half saturation constant ($a$).

Lastly for the infected predators ($P_I$) depicted in Figure 4(c), the significant parameters are the egg carrying capacity ($K_E$), the probability that a susceptible prey becomes infected from an infectious larvae tick ($\beta_{LM}$), the prey carrying capacity ($K_M$), preys birth and natural death rate ($\pi_M$ and $\mu_M$). Other parameters are the fear of the predator ($\tau$), predator biomass conversion efficiency ($\pi_P$), predators natural death rate ($\mu_P$), the predators attack rate ($\alpha_M$), and the half saturation constant ($a$).

The fear term $\tau$ as for model (3.1) have a variable sensitivity indices for the three response function. It has a negative impact on the sum of the infected ticks, and the predator, and a positive significant impact on the prey population; however, the magnitude of the index is small.

It is interesting to note that the infection transmission probabilities have mostly non-significant impact on all the response functions $L_I+N_I+A_I$, $M_I$, and $P_I$ except for the probability that a susceptible prey becomes infected from an infectious larvae tick ($\beta_{LM}$).

Parameters with the largest partial rank correlation coefficient (PRCC) values have the largest impact on each of the response functions. The parameters with positive PRCC values have a positive impact on the response function and lead to increase in their values. While the parameters with negative PRCC values will reduce the response function values.

5.   Numerical simulations

In this section, we simulate models (2.1) and (3.1) using initial conditions and parameter values that lead to unstable and stable dynamics, see Section 3.1. We start with the tick-prey-predator model (3.1) without disease and choose the unstable equilibrium due to the interesting dynamics exhibited by this equilibrium point. We also simulate the model (3.1) to show the impact of the attack rate and fear of the predator following the result of the sensitivity analysis in Section 4. Lastly, we simulate the ticks-prey-predator model (2.1) with disease to show the disease transmission dynamics given the interactions between the ticks, prey, and predator.

5.1.   Numerical simulations of model (3.1)

5.1.1.   Stable and unstable dynamics

We simulate the ticks-prey-predator model (3.1) using initial conditions and parameters values that lead to unstable and stable dynamics and focus on the prey-predator interaction. Figure 5(a) shows the dynamics when the equilibrium is unstable; while Figure 5(b) shows the dynamics when the equilibrium is stable. We see in Figure 5(a) an oscillatory dynamics as the prey and predator interact. This oscillatory behavior dissipates in a stable environment; see Figure 5(b). A stable attractor is observed in Figure 5(c). Here, the trajectory of the initial condition (IC) converges toward a stable limit cycle.

Next, we simulate the ticks-prey-predator model (3.1) for the unstable prey-predator equilibrium to see the effect of the prey-predator interaction on the tick dynamics. Figure 6(a) shows the prey-predator interaction as in Figure 5. In Figure 6(b) we plot the simulation result of larvae and prey together. The figure shows the larvae population following the trajectory to the prey population. Similarly, we plot the simulation results of the eggs, nymphs, adults, along with predator. We also observed in Figure 6(c) the nymphs and the adult population following, in this case, the trajectory of the predator population, but that is not the case with the eggs.

5.2.   Numerical simulations of model (2.1) with disease

We simulate the ticks-prey-predator model (2.1) with disease using initial conditions and parameters values that lead to unstable dynamics and show the disease transmission dynamics given the interactions between the ticks, prey, and predator. We observed oscillatory dynamics in susceptible and infected populations as prey and predator interact (see Figure 7). These oscillatory dynamics are similar to the dynamics observed in model (3.1) without disease. Figure 7(a) shows the dynamics among the susceptible ticks, and prey population while Figure 7(b) shows the dynamics of the infected ticks and the prey population.

We observed in Figure 7(a) the trajectory of the susceptible larvae population follow the trajectory of the susceptible prey population as the model without disease. While susceptible nymphs and adults follow the trajectories of susceptible predators, see Figure 7(c). For the infected larvae in Figure 7(b), we observed more infected larvae when the infected prey population are low and few infected larvae when there are more infected prey. We also observed in Figures 7(c) and 7(d) the trajectories of the nymphs and adult population either susceptible or infected follow the trajectory of the predator population (susceptible or infected).

5.2.1.   Simulation scenario based on sensitivity analysis

Next, we run some simulation scenario using the results from the sensitivity analysis varying the prey-predator parameters $\alpha_M$ (predator attack rate), $\tau$ (strength of fear of predators). The parameter $\alpha_M$ have significant negative impact on the response functions for the sum of infected ticks $L_I+N_I+A_I$, infected prey ($M_I$) and infected predator $P_I$; while the impact of $\tau$ is minimal and positive on the sum of infected ticks $L_I+N_I+A_I$, and infected prey ($M_I$). It is however, has negative impact on $P_I$. Although the impact of this parameter is low and insignificant, we include this in our simulation because of it's biological importance.

5.2.2.   Varying predator attack rate ($\alpha_M$)

We set the predator attack rate $\alpha_M$ to 0.5, 0.75, 0.95 and simulate model (3.1) to explore the effect of predator attack rate ($\alpha_M$) on the system. We observed from Figure 8(a) a significant decrease in the susceptible predator population ($P$) as the attack rate $\alpha_M$ increases with a slight increase in the prey population ($M$) during the times the prey population increases following a decrease in the predator population. Figure 8(b) show the infected populations. We observed a decrease in both prey and predator populations with increasing attack rate; however, there are more infected prey when the predator are low. Furthermore, we observed in Figures 8(c) and 8(d) the cascading effect of the predator attack on prey who are the hosts to the larvae and nymph tick populations; thus as the attack rate increases we observed a decrease in the tick populations namely the larvae, nymphs, and adults.

5.2.3.   Varying the fear of the predator ($\tau$)

Next, we vary the fear of the predator and set the fear parameter $\tau$ to 1.50, 2.50, 4.50, and simulate model (2.1) to explore the effect of the predator fear rate ($\tau$) on the system. We observed from Figure 9(a) a decrease in the susceptible predator population ($P_S$) as $\tau$ increases with a horizontal shift to the right in the prey population ($M_S$) when the prey population increases following a decrease in the predator population. Figure 9(b) show the infected populations ($P_I, M_I$). Both the prey and predator populations decrease with increasing attack rate; however, there are more infected prey when the predator are low. Figures 9(c) and 9(d) shows the cascading effect of the fear of predator on susceptible and infected tick populations; thus as the fear rate increases, we observed a decrease in larvae, nymphs, and adults tick populations; although this decrease is not as much as what we observed in Figure 8(b) above.

5.2.4.   Varying attack ($\alpha_M$) and fear ($\tau$) rates

Lastly, we vary together the predator attack rate ($\alpha_M$) and the fear of the predator ($\tau$) and simulate model (2.1) to explore the joint effect of the predator attack rate ($\alpha_M$) and prey fear rate ($\tau$) on the system. We observed a significant decrease in the susceptible predator population ($P_S$) in Figure 10(a) as $\tau$ increases with a horizontal shift to the right in the susceptible prey population ($M_S$) at the times the prey population increases following a decrease in the predator population. Figure 10(b) show the infected populations. We observed a decrease in both prey and predator populations with increasing attack rate; however, there are more infected prey when the predator when the prey are low. We also observed a strong cascading effect of the fear of predator on ticks in Figure 10(c) and 10(d); thus leading to a decrease in the larvae, nymphs, and adults tick populations as the attack and fear rates increases.

5.2.5.   Varying transmission probability ($\beta_{LM}$)

Next, we vary the probability that a susceptible prey becomes infected from an infectious larvae tick ($\beta_{LM}$) and set $\beta_{LM} = 0.018, \; 2\times0.018, \; 3\times0.0.18$ in Figure 11. We also vary the probability that a susceptible prey and predator becomes infected from an infectious nymphs and adult tick, and vice versa. We then double and triple $\beta_{LM}, \; \beta_{NM}, \; \beta_{NP}, \; \beta_{AP}, \; \beta_{PA}, \; \beta_{MT}$. Although these parameters have low sensitivities index from the sensitivity analysis results. We observed that as the transmission probabilities increases the susceptible populations (host and ticks) decreases and the infected populations increases.

5.2.6.   Summary of the simulation scenarios

In this section we ran some simulation scenarios using some of the parameters with the most impact from the sensitivity analysis, specifically we varied the prey-predator parameters $\alpha_M$ (predator attack rate), $\tau$ (strength of fear of predators). We also varied the transmission probabilities$\beta_{LM}, \; \beta_{NM}, \; \beta_{NP}, \; \beta_{AP}, \; \beta_{PA}, \; \beta_{MT}$, although these parameters had minimal effect on the response functions used in the sensitivity analysis. We observed a cascading effect of the predator's attack and fear on the prey and the tick populations, since the prey are hosts of the larvae and nymphs. We summarized below the findings from the simulations

$({\rm{i}})$ Varying predator attack rate ($\alpha_M$): As the attack rate on the prey increased the tick populations namely the larvae, nymphs, and adults we decreased;

$({\rm{ii}})$ Varying the fear of the predator ($\tau$): As the fear rate increased, the larvae, nymphs, and adults tick populations decreased;

$({\rm{iii}})$ Varying transmission probability: As the transmission probabilities increased the infected populations increased leading to a decrease in the susceptible populations both host and ticks alike.

The implications of these results are discussed next in the section below.

6.   Discussion

In this paper we developed and analyzed two mathematical models involving ticks-prey-predator interactions model (2.1) with disease and model (3.1) without disease. We use Erhlichiosis as our model disease transmitted by ticks, see Erhlichiosis model in ^[24]. The quantitative analysis of model (2.1) with disease transmission, indicate that the disease-free equilibrium is locally asymptotically stable when the reproduction number is less than one and unstable otherwise.

To carry out the quantitative analysis of ticks-prey-predator model (3.1), we first simplify the model. The quantitative analysis of the simplified form using boundedness of the prey and predator indicates that the model has three equilibrium points: tick-only equilibrium, predator-free equilibrium, and coexistence ticks-prey-predator equilibrium. These three equilibrium points are locally asymptotically stable given the conditions in Theorem 3.

Next we carried out a global sensitivity analysis using LHS/PRCC method to determine the parameter with the most influence on the response functions obtained from the ticks-prey-predator model (3.1) without disease. The parameters with the most significant impact on the sum of ticks ($L+N+A$), the prey ($M$), and the predators ($P$) are the tick eggs maturation rate ($\sigma_E$), the egg carrying capacity ($K_E$). Other significant parameters are the preys' birth and natural death rate ($\pi_M$ and $\mu_M$) and their carrying capacity ($K_M$). The parameter denoting the fear of predators ($\tau$) is equally significant. The last sets of significant parameters are the predator birth and natural death rates ($\pi_P$, and $\mu_P$), the predator's attack rate ($\alpha_M$), and the half saturation constant ($a$).

We also carried out a global sensitivity analysis using LHS/PRCC method to determine the parameters of the ticks-prey-predator disease model (2.1) with the most influence on the response functions namely the sum of infected ticks ($L_I+N_I+A_I$), infected prey ($M_I$), and the infected predators ($P_I$). The most significant parameters for these functions shown in Figure 4 and given in Table 3. These parameters are tick eggs maturation rate ($\sigma_E$), the egg carrying capacity ($K_E$), the probability that a susceptible prey becomes infected from an infectious larvae tick ($\beta_{LM}$). The other significant parameters are the preys birth and natural death rate ($\pi_M$ and $\mu_M$), and their carrying capacity ($K_M$). The other equally significant parameters are the predator biomass conversion efficiency ($\pi_P$) or predators birth rate, the natural death rate ($\mu_P$), the attack rate ($\alpha_M$), and the half saturation constant ($a$). The other parameters are the fear of the predator ($\tau$).

The ticks-prey-predator interactions models (2.1) and (3.1) with and without disease have several significant parameters (like $\sigma_E, K_E, \pi_M, \mu_M, K_M, \pi_P, \mu_P, \alpha_M$, and $\tau$) that are common given the three response functions under consideration. The parameter $\beta_{LM}$ from model (2.1) is the only parameter not common to both models.

The fear term $\tau$ in models (2.1) and (3.1) have variable sensitivity indices for the three response functions. It has a negative impact on the sum of the infected ticks, and the predator, and a positive significant impact on the prey population, although their magnitude if small. This indicate that the effect of fear is negative on the ticks and predator and positive on the prey. This positive effect on the prey we expect is in the short term because studies have show that fear can affect the birth rate of prey population ^[45], limit foraging activities ^[46]. For instance, mule deer limit their forage due lion attacks ^[46]. Thus, in the long run that effect on the prey would be negative since the fear of the predator would prevent them from foraging for food and even mating leading to a decrease in their population over time ^[45,46,47].

In Figure 6, we focus on the simulation of the ticks-prey-predator model (3.1) for the unstable prey-predator equilibrium shown in Figure 5(a) to see the effect of the prey-predator interaction on the ticks dynamics. Figure 6(a) shows the prey-predator dynamics. In Figure 6(b), we see that the larvae population follows the trajectory of the prey population, this is due to a number of factors: for instance, the larvae feeds on the prey. Furthermore, the larvae mortality which is density dependent is also dependent on the prey population; hence, the larvae population decreases with the prey population and increases when the prey population increases. We also observed in Figure 6(c) that the nymphs and the adult follow the predator population trajectory since the egg laying rate, nymph, and adult mortality both depends on the predator population. Additionally, the nymph population also depends on the prey population. However, from our numerical simulations, we observed dominance by the predator population in terms of following the trajectory. This is quite interesting and worthy of further investigation to ascertain the true cause of dominance by the predator population. This observation substantiate the results obtained in ^[48]. Although this interactions is between coral reef fishes; however, the authors emphasized the importance of predation and competition within hierarchies between two coral reef fishes.

Considering the effect of varying the predators attack rate $(\alpha_M)$ is crucial in understanding disease outbreaks within the populations. In Figure 8, we investigate how changes in the predators attack rate impact disease transmission. By adjusting $\alpha_M$ from 0.5 to 0.95, we observed distinct outcomes. We observed in Figure 8(a) a decline in the populations of infected larvae, nymphs, and adult ticks as the predators attack rate on the prey increases. This decline indicates that the infected predators play a significant role in reducing the prevalence of tick-borne diseases by preying on the infected ticks. Consequently, the overall tick burden decreases, leading to a potential reduction in disease transmission to hosts. This observation from our numerical simulations substantiate experimental results in ecology where predators and competitors of vertebrates can in theory reduce the density of infected nymphs ^[18], and thus may have a cascading impacts on the tick-borne disease risk. However, in Figure 8(b), we observe a different trend. Here, despite the increase in predator activity, we notice a decrease in the predator population and a marginal rise in the prey population. This unexpected outcome suggests a possible trade-off between predator survival and prey control efficiency. The decline in the infected predator population might be attributed to the additional time and energy expended in intensifying attacks on infected prey individuals. As a result, while predator activity may initially reduce disease prevalence, sustained high attack rates could lead to a decline in predator fitness, potentially impacting long-term disease regulation dynamics.

The non-lethal effect of predation due to fear $(\tau)$ is further explored in Figure 9. We observed a decline and shift in the dynamics of the infected ticks populations, i.e., larvae, nymphs, and adult ticks as the strength of fear parameter increases, see Figure 9(a). This observation indicates that intensifying the strength of fear of infected predators among infected prey populations lead to a decrease in tick burden within the ecosystem, thus a cascading effects on the tick-borne disease risks. Also, in Figure 9(b), we observed a decrease and shift in the infected predator populations and a slight decrease in the population of the infected prey. These could be due to several factors. Intensifying the strength fear of predators may have resulted in reduced encounter rates between infected predators and infected prey, limiting opportunities for predation and thus impacting the infected predator populations. In a related study ^[49], researchers found similar results. They noticed that as the fear of predation increased, the predator population decreased, and in some cases, high levels of fear even caused the predators to become extinct. Additionally, an increase in the strength of fear $(\tau)$ could lead to changes in infected predator behavior, such as decreased hunting efficiency. Conversely, the slight decrease in infected prey populations may be as a result from a combination of reduced predation pressure and altered prey behavior in response to fear. The infected prey populations may exhibit avoidance behaviors or alter their habitat use patterns to reduce encounters with infected predators, leading to a decrease in infected prey densities over time ^[50].

In Figure 10, we carried out numerical simulations where we simultaneously increased both the infected predators' attack rate $(\alpha_M)$ and the strength of infected predators fear $(\tau)$ to observe their combined impact on the dynamics of the system. Our observations revealed a trend consistent with the findings from individual parameter variations. As we increased both $\alpha_M$ and $\tau$ simultaneously, we observed a decrease in the populations of infected larvae, nymphs, and adult ticks, see Figure 10(a). Similarly, we observed a decline and a shift in the predators population. Therefore, the combined effect of predators attack rate and the intensity of fear can have intricate consequences on populations of infected prey and predators, consequently influencing the risk of tick-borne diseases.

7.   Conclusion

To conclude, in this study we developed a deterministic model of ordinary differential equations to gain insight on the effect of prey-predator interaction and the cascading effect of fear of the predators on ticks population and disease prevalence. We found that tick host predator and fear of the predators reduces the tick populations and disease prevalence. We summarize the other results as follows:

$({\rm{i}})$ The results of the Sensitivity analyzes was implemented using three response or output functions namely the sum of infected ticks ($L_I+N_I+A_I$), infected predators ($M_I$), and the infected predators ($P_I$) to identify the parameters with the most impact on these functions in no particular order. The most significant common parameters among these functions are the eggs maturation rates ($\sigma_E$), and carrying capacity ($K_E$). Other significant parameters are the preys birth and natural death rates ($\pi_M, \mu_M)$. The last sets of significant common parameters are the predator birth and natural death rates ($\pi_P$, and $\mu_P$), the predator's attack rate ($\alpha_M$), and the half saturation constant ($a$). The probability that a susceptible prey becomes infected from an infectious larvae tick ($\beta_{LM}$) is also a significant parameter.

$({\rm{ii}})$ The trajectory of the susceptible larvae population follow the trajectory of the susceptible prey population. While susceptible nymphs and adults follow the trajectories of susceptible predator. The trajectories of the nymphs and adult population either susceptible or infected follow the trajectory of the predator population (susceptible or infected).

$({\rm{iii}})$ There are more infected larvae when infected prey population are low and few infected larvae when there are more infected prey. Similar dynamics was observed for the infected nymphs and adult ticks and infected predator population.

$({\rm{iv}})$ Prey predation and fear of predator have cascading effect on tick populations and disease prevalence. As prey attack rate increases, there is corresponding reduction in the prey and ticks population; furthermore as the fear of the predator increases we observed a reduction in the prey population which subsequently lead to a decrease in the ticks populations.

In this study, we have used mathematical models to investigate the cascading effect of prey-predator interaction and predator fear tick populations and disease prevalence. One of the limitations of this study is not including additional source of food for the predator. We know that often times, the predator have other sources of food aside from the prey populations, see ^[51] and references therein. Another limitation of the study is that we were not specific as to the kind of prey population we are interested in; however, tick larvae and nymphs are known to feed on small mammals while the adult ticks feed on much larger mammals including the predators which we capture in the two models in this study. Another major drawback was the availability of data. Although, we have off set this drawback by using parameter estimates from literature; we are confident in our conclusions since are results are backed up by the outcomes obtained from the sensitivity analysis we carried out.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This material is based upon work supported by the National Science Foundation (NSF) under Grant No. DMS-1928930, National Security Agency (NSA) under Grant No. H98230-23-1-0012 and the Sloan Foundation under Grant G-2020-14104 while the author participated in ADJOINT hosted by the Mathematical Sciences Research Institute in Berkeley, California, in [2023-2024]. FBA had additional support from the National Science Foundation under the EPSCOR Track 2 grant number 1920946.

Conflict of interest

The authors declare there is no conflict of interest.