# Dynamics of Rabies Epidemics in Vampire Bats

Authors: Liang Tian; Juping Zhang
Journal: Complexity
 (2020)

---

## Abstract

In order to study the transmission of rabies epidemics in vampire bats, we propose a mathematical model for vampire bat rabies virus. A thresholdR0 is identified which determines the outcome of the disease. If R0<1, the disease-free equilibrium is globally asymptotically stable, and if R0>1, the endemic equilibrium is globally asymptotically stable with certain conditions. Through the numerical simulation, the correctness of the theoretical results is verified. We carry out the sensitivity analysis of the parameters which provide a theoretical basis for preventing and controlling the transmission of bat rabies.


---

## Body

## 1. Introduction

Rabies is a zoonotic disease caused by acute infectious viruses. Most rabies models have focused on disease dynamics and control within terrestrial mammals (e.g., dogs and foxes). Hou et al. [1] studied the prevalence and control of canine rabies in Guangdong Province, China. Asamoah et al. [2] studied an optimal control model of rabies transmission dynamics in dogs and the best way of reducing the death rate of rabies in humans. In fact, vampire bats are a host of zoonotic and potential zoonotic pathogens [3] and are also effective vectors for rabies. Most vampire bats transmit the virus through bites via infected saliva, but when the population density is high, they also spread the virus through aerosols. The habits of vampire bats are particularly exceptive. In the case of Mexican vampire bats, they mainly feed on the blood of animals. The unique tooth structure provides conditions for them to open the skin, absorb blood, and diffuse the rabies virus. The latent period of bat rabies is about 21–209 days. The longer the latent period, the more likely infected bats will survive long enough to enter hibernation and provide infectious contacts in the transmission season [4]. The target of bat rabies transmission is mainly large-scale warm-blooded animals such as cattle and sheep [5].VBRV has been spreading for more than 100 years [5]. Estimating from Brazil alone suggested that across Latin America at least thousands of cattle were infected with rabies each year, which became a major burden for the livestock industry [6]. Furthermore, human deaths related to vampire bats had been on the rise, with the earliest description dating back to the Spanish conquest of the Americas in the 16th century [7]. Chile reported the first case of an insect-eating bat carrying rabies virus infecting a human [8]. Martinez-Burnes et al. [9] experimentally confirmed the presence of rabies virus in vampire bats. Therefore, VBRV attracted close attention and in-depth study by researchers. They had different concerns about vampire bats and conducted different research. Obregón-morales et al. [10] demonstrated that rousette bats exposed to rabies had an innate low susceptibility to bat rabies. Based on five-year observations of large brown bats, George et al. [4] established a bat rabies seasonal transmission model and analyzed key parameters that determine the seasonality and maintenance of bat rabies. Blackwood et al. [11] used longitudinal study data from vampire bats to proposed four bat rabies models and developed a parameterized and evaluated stochastic performance model to assess the determinants of virus persistence in vampire bats. Aguilar-setien et al. [12] studied the saliva excretion of rabies virus by vampire bats and indicated that bats that were infected with rabies virus and survived will expel the virus through saliva in the early stage of infection. Kobayashi et al. [13] confirmed that there is a considerable diversity in bat rabies variants through sequence analysis of the Brazilian bat rabies virus, which provided conditions for the epidemic of bat rabies. Bat population was monitored for rabies by fluorescent antibody testing and simplified fluorescent inhibition microtesting by Almeida et al. [14], who indicated an increase in rabies virus circulation among bats and risk of a rabies outbreak. Streicker et al. [15] highlighted the importance of international viral dispersal for shaping the burden of rabies in Costa Rica and suggested a Central American corridor of rabies virus invasions between continents. Although researchers have carried out a lot of studies on bat rabies from different aspects, there are relatively few studies on relevant mathematical analysis from the perspective of the dynamic model. The transmission mechanism of bat rabies cannot be intuitively and objectively reflected by the mathematical model, and the key factors affecting the transmission cannot be analyzed by sensitivity analysis. In this paper, we will use mathematical derivation to study the bat rabies model and provide a theoretical basis.The remaining part of this paper is organized as follows. In Section2, we first propose the model of rabies epidemics in vampire bats. The next five sections are devoted to studying the stability of the established model. In Section 3, we derive the threshold R0. It is shown that if R0<1 then the disease-free equilibrium is globally asymptotically stable. In Sections 4–7, the model has a unique endemic equilibrium, which asymptotically attracts all nontrivial solutions under certain conditions. These theoretical results are supported with numerical simulations. Sensitivity analysis of R0 on various parameters is carried out at the end of Section 8. The paper concludes with a brief conclusion.


## 2. Model Formulation

According to the transmission mechanism of rabies virus in bat population and the immunity of bats to rabies virus, the transfer diagram of the model is given by Figure1.Figure 1
The transfer diagram of the SEIVS model.Therefore, we formulate the dynamics model of rabies transmission in vampire bats:(1)dSdt=bN+εV−βSIN−μS,dEdt=βSIN−μE−τ1E,dIdt=ατ1E−τ2I−μI−δI,dVdt=1−ατ1E+τ2I−μV−εV.A population sizeNt is partitioned into subclasses of individuals who are susceptible, exposed (infected but not yet infectious), infectious, and immune, with sizes denoted by St, Et, It, and Vt, where Nt=St+Et+It+Vt. The meanings of related parameters in (1) and their values are shown in the Table 1 (time scale is week).Adding the four equations of system (1),(2)dNdt=bN−μN−δI.


## 3. Threshold and Stability of Disease-Free Equilibrium

Lets=S/N,e=E/N,i=I/N, and v=V/N denote the fractions of the classes S,E,I, and V in the population, respectively. It is easy to verify that s,e,i, and v satisfy the system of differential equations:(3)dsdt=b+εv−βsi−bs+δsi,dedt=βsi−τ1e−be+δei,didt=ατ1e−τ2i−δi−bi+δi2,dvdt=1−ατ1e+τ2i−εv−bv+δvi,and obviously s+e+i+v=1. Note that the total population size Nt does not appear in (3). Also, observe that the variable v does not appear in the first three equations of (3). This allows us to attack (3) by studying the system:(4)dsdt=b+ε1−s−e−i−βsi−bs+δsi,dedt=βsi−τ1e−be+δei,didt=ατ1e−τ2i−δi−bi+δi2.From biological considerations, we study (4) in the closed set:(5)Γ=s,e,i∈ℝ+30≤s,e,i≤1,0≤s+e+i≤1,where ℝ+3 denotes the nonnegative cone of ℝ3. It can be verified that Γ is positively invariant set with respect to (4).According to system (4), it is easy to get the disease-free equilibrium of E01,0,0. At the disease-free equilibrium of E0, the next-generation matrix theory [20] is used to derive(6)F=0β00,V=τ1+b0−ατ1τ2+δ+b.Let(7)R0=ατ1βb+τ1b+τ2+δ.Following [21], R0 is called the threshold. Therefore, when R0<1, the disease-free equilibrium E0 is locally asymptotically stable in Γ. In the following result, the global asymptotic stability of E0 should be understood in the sense of Lyapunov.Lemma 1 (see [22]).
LetΛ=e,i∈ℝ+20≤e+i≤1 and he,i=a1−b1e+c1−b1i+b1. Therefore, for any positive constants a1,b1, and c1, maxe,i∈Λhe,i=maxa1,b1,c1.Theorem 1.
The disease-free equilibriumE0=1,0,0 of (4) is globally asymptotically stable in Γ if R0<1.Proof.
LetL=ατ1e+b+τ1i. Then,(8)L′=ατ1βsi+ατ1δei−b+τ1b+τ2+δi+b+τ1δi2≤ατ1β1−e−ii+ατ1δei−b+τ1b+τ2+δi+b+τ1δi2=he,i−b+τ1b+τ2+δ,where(9)he,i=ατ1β1−e−i+ατ1δe+b+τ1=ατ1δ−ατ1βe+b+τ1δ−ατ1βi+ατ1β.
Applying Lemma1 to he,i leads to(10)L′=≤maxατ1δ,ατ1β,b+τ1δ−b+τ1b+τ2+δ≤0,ifR0<1.
IfL′=0 and i≠0, then(11)maxατ1δ,ατ1β,b+τ1δ=b+τ1b+τ2+δ.
Thus,L′=0 only if either i=0 or R0=1, that is, s+e+i=1. The maximum invariant set in s,e,i∈ΓL′=0 is the singleton E0. The global stability of E0 when R0<1 follows from LaSalle’s Invariance Principle [23].


## 4. Existence and Uniqueness of Endemic EquilibriumE∗

From Theorem1, system (4) admits a disease-free equilibrium E01,0,0 if R0<1 and E0 is unstable when R0>1. Furthermore, we remark that R0>1 implies β>δ. If R0>1, system (4) has a unique endemic equilibrium E∗s∗,e∗,i∗ and satisfies(12)b+ε1−s∗−e∗−i∗−βs∗i∗−bs∗+δs∗i∗=0,βs∗i∗−b+τ1e∗+δe∗i∗=0,ατ1e∗−b+τ2+δi∗+δi∗2=0,and thus s∗>0,e∗>0, and i∗>0. Adding the abovementioned equations leads to(13)b+ε−δi∗1−s∗−e∗−i∗=τ2i∗+1−ατ1e∗,which gives the following range of i∗:(14)0<i∗<min1,b+εδ.According to the first and third equations of system (12), we have(15)e∗=b+τ2+δ−δi∗i∗ατ1,s∗=b+ε1−i∗−εe∗b+ε+β−δi∗.Using (15), we eliminate s∗ and e∗ in the second equation of (12) and derive that i∗ satisfies(16)hi∗=a0i∗3+a1i∗2+a2i∗+a3=0,where(17)a0=δ−βδ2<0,a1=δεβ−b+εδ+β−δτ1+τ2+δ+2b,a2=ατ1βδ−βR0−ε+b+ετ1+τ2+δ+2bδ−b+τ2+δβε,a3=ατ1βb+ε1−1R0>0.Therefore,h0=a3>0 and h+∞=−∞<0. There is at least one positive equilibrium in the interval 0,+∞ for Intermediate Value Theorem.We see thati∗ satisfies(18)fi∗=gi∗,where(19)fi∗=1−δτ1i∗1−δb+τ2+δi∗1+β−δb+εi∗,gi∗=R0+εβδb+εb+τ1b+τ2+δi∗2−εβατ1+τ2+δ+bb+εb+τ1b+τ2+δi∗.WhenR0>1, assume that b>δ, discussing the number of intersections of fi∗ and gi∗ at [0, 1]. Obviously, g0=R0>1=f0, and the symmetry axis of gi∗ is greater than 1. Suppose that f1<g1, that is,(20)b+τ2βε+δ2+τ12+τ1bατ1βb+b+τ2βδ>1.Then,fi∗ and gi∗ have only one intersection in [0, 1]. Therefore, we see that (4) has a unique endemic equilibrium of E∗ at [0, 1] (see Figure 2).Figure 2
The existence and uniqueness ofi∗ in the interval [0, 1].


## 5. Local Asymptotic Stability of the Endemic EquilibriumE∗

To show the asymptotic stability of the endemic equilibriumE∗, we use Routh–Hurwitz conditions to proof it.Theorem 2.
IfR0>1, b>δ, and τ2+bβε+δ2+τ12+τ1b/ατ1βb+τ2+bβδ>1, then (4) has a unique endemic equilibrium E∗ in Γ and E∗ is locally asymptotically stable.Proof.
The Jacobian matrix of (4) at the endemic equilibrium E∗=s∗,e∗,i∗∈Γ is(21)JE∗=−b−ε−βi∗+δi∗−ε−ε−βs∗+δs∗βi∗−τ1−b+δi∗βs∗+δe∗0ατ1−τ2−δ−b+2δi∗.
Therefore, the characteristic equation corresponding to the Jacobian matrix is(22)λ3+b1λ2+b2λ+b3=0,where(23)b1=3b+τ1+τ2+ε+δ+βi∗−4δi∗,b2=τ2+δ+b−2δi∗2b+ε+τ1+βi∗−2δi∗+b+ε+βi∗−δi∗τ1+b−δi∗,b3=b+ε+βi∗−δi∗τ1+b−δi∗τ2+δ+b−2δi∗.
Using Routh–Hurwitz conditions, we have the following results ofb1>0 and b1b2−b3>0. This verifies local asymptotic stability of the endemic equilibrium E∗.


## 6. Uniform Persistence

To prove the global stability ofE∗, in this section, we investigate the uniform persistence of (4) with the uniform persistence theory of Hale and Waltman [24].Definition 1.
LetX be a complete metric space. Suppose that X0 is an open set in X, ∂X0⊂X, X0∩∂X0=∅. Assume that Tt is a C0-semigroup of X satisfying(24)Tt:X0⟶X0,Tt:∂X0⟶∂X0.
LetT∂t=Tt∂X0, and let A∂ be the global attractor for T∂t.Definition 2.
The semigroupTt is said to be point dissipative in X if there is a bound nonempty set B in X such that, for any x∈X, there is t0=t0x,B such that Ttx∈B for t≥t0.Definition 3.
A nonempty invariant subsetM of X is called an isolated invariant set if it is the maximal invariant set of a neighborhood of itself.
The orbitγ+x through x is defined as γ+=∪t≥0Ttx. The ω-limit set is defined as(25)ωx=∩τ≥0C1∪t≥τTtx,where C1 is a closure.
Thestable set of a compact invariant set A is denoted by Ws and is defined as(26)WsA=xx∈X,ωx≠∅,ωx⊂A.
The particular invariant set of interest is(27)A∂˜=∪x∈A∂ωx.Definition 4.
A
∂
˜ is isolated if there exists a covering M=∪i−1kMi of A∂˜ by pairwise disjoint, compact, and isolated invariant sets M1,M1,…,Mk for T∂ such that each Mi is also an isolated invariant set for T. M is called an isolated covering.Definition 5.
A
∂
˜ will be called acyclic if there exists some isolated covering M=∪i−1kMi of A∂˜ such that no subset of Mi’s forms a cycle.Lemma 2.
Suppose thatTt satisfies (4) and that the following conditions are valid:(i)
There ist0>0 such that Tt is compact for t>t0(ii)
Tt is point dissipative in X(iii)
A∂˜=∪x∈A∂ωx is isolated and has an acyclic covering M, where M=∪i−1kMi(iv)
WsMi∩X0=∅, i=1,2,…,n
Then,Tt is uniformly persistent.Theorem 3.
IfR0>1, system (4) is uniformly persistent.Proof.
LetC−τ,0,ℝ+3 be a Banach space which consists of −τ,0 to ℝ+3 continuous functions and let C˜=φ1,φ2,φ3∈−τ,0,ℝ+3,φ1θ≠0,φ2θ=φ3θ≡0,θ∈−τ,0. We assume X0=C−τ,0,ℝ+3, ∂X0=C˜. Then, (4) has a constant solution, such that E˜⊂C˜, where E˜=φ1,φ2,φ3∈−τ,0,ℝ+3,φ1θ=b+σ1−s,φ2θ=φ3θ≡0,θ∈−τ,0.
Defineb+σ1−s=s0. Using Definition 1, (4) satisfies conditions (i) and (ii) of Lemma 2. Furthermore, X0 and ∂X0 are invariant sets; then, (24) is also established. e′tφ1,φ2,φ3∈C˜≡0,i′tφ1,φ2,φ3∈C˜≡0 deduce that etφ1,φ2,φ3∈C˜≡0 and itφ1,φ2,φ3∈C˜≡0 are satisfied for any t≥0. Therefore, s′tφ1,φ2,φ3∈C˜=b+σ1−s, that is, stφ1,φ2,φ3∈C˜=s0. As a result, all the points in C˜ tend to E˜, that is, C˜=WsE˜. Obviously, A∂˜=E˜ is isolated. We know that the flow of system (4) is an acyclic of A∂˜; thus, condition (iii) of Lemma 2 is established.
Next we proveWsE˜∩X0=∅. Using the counterevidence, we assume WsE˜∩X0≠∅; there is a positive solution for system (4) st,et,it satisfying limt⟶+∞st,et,it=s0,0,0. Let σ be any sufficiently small positive constant; there is η=ησ, such that(28)st>s0−σ>0,et<σ,it<σ,t≥η.
From the second equation of system (4), we know that(29)e′t>βs0−σit+δetit−τ1+bet.
Thus, we can construct the following auxiliary system:(30)m1′t=βs0−σm2t+δm1tm2t−τ1+bm1t,m2′t=ατ1m1t+δm22t−τ2+δ+bm2t,m1t≡et,m2t≡it,t∈η,η+σ,which obtains et≥m1t and it≥m2t, using comparison principle for any t≥η. If R0>1, when b>δ, the unique endemic equilibrium m1∗,m2∗ of the auxiliary system has(31)limt⟶+∞m1t=Q1δ+Q2s0−σβ2b+2τ1−Q1δ−Q2δ≔m1∗>σ,limt⟶+∞m2t=Q1δ+Q2≔m2∗>σ,where Q1=τ1+τ2+δ+2b and Q2=Q12δ2+4ατ1βδ2s0−σ−1/R01/2. Thus, liminft⟶+∞et≥m1∗>σ and liminft⟶+∞it≥m2∗>σ for t>η. This result contradicts the boundedness of the solution. Then, all of the conditions of Lemma 2 are valid. This establishes the theorem.


## 7. Global Asymptotic Stability of Endemic EquilibriumE∗

Theorem 4.
IfR0>1, suppose that τ2+bβε+δ2+τ12+τ1b/ατ1βb+τ2+bβδ>1 and b>δ; there is a unique endemic equilibrium E∗. When 2δ−b<ε<minδ,τ1/2, E∗ is globally asymptotically stable.Proof.
LetPs,e,i=diag1,e/i,e/i; then, P is C1 and nonsingular in Γ. Let f denote the vector of system (4); then, PfP−1=diag0,e′/e−i′/i,e′/e−i′/i. The second compound matrix [25] A2 associated with a general solution xt=st,et,it of system (4) takes the form(32)A2=−2b−τ1−ε−βi+2δiβs+δeε+βs−δsατ1−2b−ε−δ−τ2−βi+3δi−ε0βi−τ1−τ2−2b−δ+3δi.
Then,(33)PfP−1+PA2P−1=X+2δiβs+δeε+βs−δsieατ1eiY+3δi+e′e−i′i−ε0βiZ+3δi+e′e−i′i,where X=−2b−ε−τ1−βi,Y=−2b−ε−τ2−δ−βi,andZ=−τ1−τ2−δ−2b. Let(34)B=PfP−1+PA2P−1,that is, B=B11B12B21B22, where B11=X+2δi,B12=βs+δe,ε+βs−δsi/e, and(35)B21=ατ1ei0,B22=Y+3δi+e′e−i′i−εβiZ+3δi+e′e−i′i.
Letu,v,w be a vector in ℝ3; we can define a vector norm in ℝ3 as(36)u,v,w=maxu,v+w.
Letμ be the Lozinskiiˇ measure with respect to this norm; then, as described in [26], we choose(37)μB≤supg1,g2,where(38)g1=μ1B11+B12,g2=B21+μ1B22.B12 and B21 are matrix norms with respect to the L1 vector norm and μ1 denotes the Lozinskiiˇ measure with respect to the L1 norm. Then,(39)B21=ατ1ei,μ1B11=−2b−ε−τ1−βi+2δi,μ1B22=−2b−ε−τ2−δ+3δi+e′e−i′i,B12=maxβs+δeie,ε+βs−δsie.
Therefore,(40)g1=−2b−ε−τ1−βi+2δi+maxβs+δeie,ε+βs−δsie,g2=−2b−ε−τ2−δ+3δi+ατ1ei+e′e−i′i.
From system (4), it follows that e′/e=βsi/e−τ1−b+δi and ατ1e/i=i′/i+τ2+δ+b−δi. We have(41)g1=e′e−b−ε−βi+2δi+ε−δs+eie,ε≥δ,e′e−b−ε−βi+2δi,ε<δ.g2=e′e−b−ε+2δi,ε≥τ12,e′e−b−τ1+ε+2δi,ε<τ12.when 2δ−b<ε<minδ,τ1/2, according to the expression of g1 and g2, g1<g2 is verified; that is,(42)μB≤supg1,g2=e′e−b−ε+2δi.
Therefore, we have(43)lim supt⟶+∞supx0∈Γ1t∫0tμBds<1tlnetlne0+b+εi−1<0,with Mean Value Theorems for Integrals. This completes the proof.
Unfortunately, we only prove one case. It is difficult to compare the sizes ofg1 and g2. Hence, other cases are not discussed in this paper.


## 8. Numerical Simulations

To verify the theoretical results, we perform numerical simulations of system (4). We assume β=0.2, δ=0.06, and τ2=0.1. The values of other parameters are shown in Table 1. Figure 3(a) shows the time-variation diagram of system (4) state variables. It is found that if R0<1, the curve tends to be stable with time, and the disease-free equilibrium E0 is globally asymptotically stable. Moreover, we take β=0.0106, δ=0.007, and τ2=0.015. As shown in Figure 3(b), it is found that if R0>1, the state variables of system (4) change with time, and the endemic equilibrium E∗ is globally asymptotically stable.Table 1
Model parameters.
Parameter
Description
Value
Source
b
Birth rate
0.5/52
[16]
1
/
ε
Mean duration of time inV class
18–20 weeks
[17, 18]
β
Coefficient of infectious transmissibility
Assumption
—
1
/
μ
Average lifespan
177 weeks
[16]
δ
Disease mortality rate
Assumption
—
1
/
τ
1
Mean duration of time inE class
3–30 weeks
[19]
1
/
τ
2
Mean duration of time inI class
Assumption
—
α
Probability of infection given exposure
0–0.89
[11]Figure 3
Time series diagram ofs, e, and i. (a) R0<1. (b) R0<1.
(a)
(b)Figure4 shows it will increase with the increase of the parameters α, β, τ1, and τ2 if R0<1. But the change of the parameter τ1 does not change significantly for the proportion of infected individuals. The larger the parameter τ2 is, the shorter the time for the proportion of infected individuals to reach stable state is. Figure 5 shows it will increase with the increase of the parameters α, β, τ1, and τ2 if R0>1. We can see that when the parameters values of α, β, τ1, and τ2 change slightly, this will cause a large change in it.Figure 4
The curve-trend diagram ofit with time, when the parameters α, β, τ1, and τ2 take different values.
(a)
(b)
(c)
(d)Figure 5
The curve-trend diagram ofit with time, when the parameters α, β, τ1, and τ2 take different values.
(a)
(b)
(c)
(d)We further simulate the effect of two different parameters on the thresholdR0 (Figure 6). When two parameters are higher simultaneously, the threshold R0 also increases. For two different parameters, the threshold changes differently.Figure 6
The effect of the parametersδ, β, τ1, and τ2 for R0.
(a)
(b)
(c)


## 9. Conclusion

This paper considers the propagation of bat rabies with latent period and establishes a dynamic model. The threshold of the model is calculated asR0, which proves that when R0<1, disease-free equilibrium E0 is globally asymptotically stable; when R0>1, under certain conditions, the system has a unique endemic equilibrium E∗. Under the restrictions, we use the geometric method [27] to prove the global stability of E∗. Finally, numerical simulations and parameter sensitivity analysis are carried out to verify the correctness of the theoretical results. This is consistent with [4]. That is, bats have an adequate latent period and a very low disease-related mortality rate of δ, which further leads to the maintenance of the virus and continuous transmission. Obviously, our results do not reflect changes of the seasonal transmission compared to [4]. In subsequent studies, we will consider the effects of bat population density constraints, seasonal birth impulses, and hibernation on virus transmission.



---

*Source: 7032451-2020-04-13.xml*