Mathematical Modeling of Canine and Human Rabies

1. Introduction

Rabies is a fatal viral zoonosis primarily affecting mammals, including humans. It is responsible for nearly 59,000 human deaths per year worldwide, the majority occurring in Africa and Asia, with dogs being responsible for over 99% of human transmissions [1,2]. Prevention mainly relies on vaccination of reservoir animals, particularly dogs, and post-exposure prophylaxis in humans. Numerous studies have confirmed that mass dog vaccination is the most effective tool for interrupting transmission, often more cost-effective than solely strengthening post-exposure prophylaxis [3,4].

However, cross-border mobility of dogs complicates control strategies, especially in border regions where inter-state coordination is often weak. For example, studies conducted at the border between Chad and Cameroon have shown the importance of considering canine migratory flows to understand the spatial dynamics of rabies [5,6]. In the Serengeti region, analyses have highlighted that regular reintroductions of the disease from neighboring areas can compromise local efforts, making large-scale vaccination necessary [7].

In several regions of the world, notably in Africa and Asia, vaccination efforts are fragmented. Models applied to African urban contexts [8] or to China [9,10] have shown that vaccination coverage must exceed a critical threshold (often ≥ 70 %) to guarantee elimination, and that high turnover in the dog population leads to a rapid decline in coverage if campaigns are not repeated annually [11].

Mathematical models have improved the understanding of rabies transmission, particularly compartmental models of the SIR [12] or SEIR type adapted to animal and human contexts [13]. These works also highlight the importance of integrating the dynamics of canine and human populations within a realistic spatial framework. However, few studies have explored the impact of cross-border dog exchanges in a multi-country setting, even though this approach is crucial for achieving the goals of the WHO’s "Zero human deaths from dog-mediated rabies by 2030" strategy [1].

The aim of this study is to propose a multi-country mathematical model for canine rabies, accounting for dog movements between two neighboring countries to better understand the impact of cross-border mobility on disease dynamics and to evaluate the effectiveness of different control strategies.

Our specific objectives include, among others:

Developing a compartmental model structured by country, integrating bidirectional dog movements between the two countries;

Calculating and analyzing the basic reproduction number ${\mathcal{R}}_0$ for this multi-country system;

Studying the stability of the rabies-free equilibrium and the conditions for endemic persistence in each country;

Performing a sensitivity analysis to identify the most influential parameters on rabies spread;

Simulating the effect of different targeted vaccination and movement control strategies.

2. Model Formulation

To describe the dynamics of rabies in two interconnected countries, we construct a compartmental model structured into two subpopulations (human and canine) for each country denoted $i=1,2$ (home country) and $j=1,2$ (neighboring country) [14]. The model is inspired by previous work on canine rabies dynamics and cross-species transmission between dogs and humans [15,16,17].

For each country i,

• The canine population $N_c^i$ is divided into compartments: susceptible ($S_c^i$), exposed ($E_c^i$), infectious ($I_c^i$), vaccinated ($V_c^i$).
• The human population $N_h^i$ is divided into susceptible ($S_h^i$), exposed ($E_h^i$), and infected ($I_h^i$).
• Vaccination is applied only to dogs.
• The model is formulated for two interconnected countries $i=1,2$.
• Dog movements between countries are negligible.

2.1. Flow Diagram

Figure 1. illustrates the flow diagram of the model, showing transitions between compartments for canine and human populations in both countries. Arrows indicate transmission pathways, infection progression, vaccination, as well as demographic movements.

Figure 1. illustrates the flow diagram of the model, showing transitions between compartments for canine and human populations in both countries. Arrows indicate transmission pathways, infection progression, vaccination, as well as demographic movements.

Table 1. Model variables for each country $i=1,2$.

Variable	Description
$S_c^i$	Number of susceptible dogs in country i
$E_c^i$	Number of exposed dogs in country i
$I_{d_i}$	Number of infected dogs in country i
$V_c^i$	Number of vaccinated dogs in country i
$S_h^i$	Number of susceptible humans in country i
$E_h^i$	Number of exposed humans in country i
$T_h^i$	Number of treated humans in country i
$I_h^i$	Number of infected humans in country i
$R_h^i$	Number of recovered humans in country i

Table 1. Model variables for each country $i=1,2$.

Variable	Description
$S_c^i$	Number of susceptible dogs in country i
$E_c^i$	Number of exposed dogs in country i
$I_{d_i}$	Number of infected dogs in country i
$V_c^i$	Number of vaccinated dogs in country i
$S_h^i$	Number of susceptible humans in country i
$E_h^i$	Number of exposed humans in country i
$T_h^i$	Number of treated humans in country i
$I_h^i$	Number of infected humans in country i
$R_h^i$	Number of recovered humans in country i

Table 2. Model parameters.

Parameter	Description
${\Lambda }_h^i$	Human recruitment rate in country i
${\Lambda }_c^i$	Dog recruitment rate in country i
${\beta }_{ch}^i$	Dog-to-human transmission rate
${\beta }_{cc}^i$	Dog-to-dog transmission rate
${\mu }_c^i$	Natural mortality rate of dogs
${\alpha }^i$	Vaccine efficacy
${\nu }_c^i$	Canine vaccination rate
${\delta }_c^i$	Disease-induced mortality rate in dogs
${\theta }^i$	Probability of receiving post-exposure treatment
${\mu }_h^i$	Natural mortality rate of humans
${\gamma }^i$	Progression rate to infection in humans
${\delta }_h^i$	Disease-induced mortality rate in humans
$m_h^X,m_c^X$	Migration rates (humans, dogs) between countries

Table 2. Model parameters.

Parameter	Description
${\Lambda }_h^i$	Human recruitment rate in country i
${\Lambda }_c^i$	Dog recruitment rate in country i
${\beta }_{ch}^i$	Dog-to-human transmission rate
${\beta }_{cc}^i$	Dog-to-dog transmission rate
${\mu }_c^i$	Natural mortality rate of dogs
${\alpha }^i$	Vaccine efficacy
${\nu }_c^i$	Canine vaccination rate
${\delta }_c^i$	Disease-induced mortality rate in dogs
${\theta }^i$	Probability of receiving post-exposure treatment
${\mu }_h^i$	Natural mortality rate of humans
${\gamma }^i$	Progression rate to infection in humans
${\delta }_h^i$	Disease-induced mortality rate in humans
$m_h^X,m_c^X$	Migration rates (humans, dogs) between countries

The systems of equations are given by:

Dogs

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{S}_c^i}{dt}}& ={\Lambda }_c^i-{\beta }_{cc}{S}_c^i{\displaystyle \frac{I_c^i}{N_c^i}}-{\nu }_c^i{S}_c^i-{\mu }_c^i{S}_c^i+m_c^S{S}_c^j-m_c^S{S}_c^i,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{V}_c^i}{dt}}& ={\nu }_c^i{S}_c^i-{\mu }_c^i{V}_c^i+m_c^V{V}_{c_j}-m_c^V{V}_c^i\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{I}_c^i}{dt}}& ={\beta }_{cc}{S}_c^i{\displaystyle \frac{I_c^i}{N_c^i}}-({\mu }_c^i+{\delta }_c^i)I_c^i+m_c^I{I}_c^j-m_c^I{I}_c^i,\hfill \end{array}$$

(1c)

where

$$N_c^i=S_c^i+V_c^i+I_c^i.$$

(1d)

Humans

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{S}_h^i}{dt}}& ={\Lambda }_h^i-{\beta }_{ch}^i{S}_h^i{\displaystyle \frac{I_c^i}{N_c^i}}-{\mu }_h^i{S}_h^i+m_h^S{S}_h^j-m_h^S{S}_h^i,\hfill \end{array}$$

(1e)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{E}_h^i}{dt}}& ={\beta }_{ch}^i{S}_h^i{\displaystyle \frac{I_c^i}{N_c^i}}-(1-{\theta }^i){\gamma }^i{E}_h^i-{\mu }_h^i{E}_h^i+m_h^E{E}_h^j-m_h^E{E}_h^i,\hfill \end{array}$$

(1f)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{T}_h^i}{dt}}& ={\theta }^i{E}_h^i-{\alpha }^i{\theta }^i{T}_h^i-{\mu }_h^i{T}_h^i+m_h^T{T}_h^j-m_h^T{T}_h^i,\hfill \end{array}$$

(1g)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{I}_h^i}{dt}}& =(1-{\theta }^i){\gamma }^i{E}_h^i-({\delta }_h^i+{\mu }_h^i)I_h^i+m_h^I{I}_h^j-m_h^I{I}_h^i,\hfill \end{array}$$

(1h)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{R}_h^i}{dt}}& ={\alpha }^i{\theta }^i{T}_h^i-{\mu }_h^i{R}_h^i+m_h^R{R}_h^j-m_h^R{R}_h^i.\hfill \end{array}$$

(1i)

with initial conditions

$$S_h^i\left(0\right),S_c^i\left(0\right)>0,\mathrm{and}{E}_h^i\left(0\right),I_h^i\left(0\right),I_c^i\left(0\right),V_c^i\left(0\right),T_h^{i}0\left(t\right),R_h^i\left(0\right)\ge 0.$$

(1j)

3. Model Analysis

3.1. Vector Form

Before proceeding further, let us specify some vector and matrix notations that will be used subsequently.

Vectors are assumed to be column vectors but are written without specific orientation unless otherwise indicated [14]. If $\mathbf{X}=(X_1,\dots ,X_n)\in {\mathbb{R}}^n$, then $\mathbf{X}\ge \mathbf{0}$ means that $X_i\ge 0$, $\mathbf{X}>\mathbf{0}$ means that $\mathbf{X}\ge \mathbf{0}$ and there exists at least one i such that $x_i>0$; finally, $\mathbf{X}\gg \mathbf{0}$ means that $X_i>0$ for all $i=1,\dots ,n$. For $\mathbf{X},\mathbf{Y}\in {\mathbb{R}}^n$, we have $\mathbf{X}\ge \mathbf{Y}$, $\mathbf{X}>\mathbf{Y}$ and $\mathbf{X}\gg \mathbf{Y}$ if, respectively, $\mathbf{X}-\mathbf{Y}\ge \mathbf{0}$, $\mathbf{X}-\mathbf{Y}>\mathbf{0}$ and $\mathbf{X}-\mathbf{Y}\gg \mathbf{0}$. This same notation applies to matrices.

Writing system (Section 2) in vector form is particularly useful for the remainder of the work. We therefore introduce some notations here.

For any variable $X\in \{S_c,V_c,I_c,S_h{E}_h,T_h,I_h,R_h\}$, we denote $\mathbf{X}=(X_1,\dots ,X_n)$, and we define:

$$\mathbf{Z}=({\mathbf{S}}_c,{\mathbf{V}}_c,{\mathbf{I}}_c,{\mathbf{S}}_h,{\mathbf{E}}_h,{\mathbf{T}}_h,{\mathbf{I}}_h,{\mathbf{R}}_h)$$

as the vector of state variables.

We also introduce the following notations:

$$\Lambda =({\Lambda }_1,\dots ,{\Lambda }_n),\mathbf{\nu }=\mathsf{diag}({\nu }_1,\dots ,{\nu }_n),\mathbf{\delta }=\mathsf{diag}({\delta }_1,\dots ,{\delta }_n),\mathbf{\gamma }=\mathsf{diag}({\gamma }_1,\dots ,{\gamma }_n),$$

$$\mathbf{\theta }=\mathsf{diag}({\theta }_1,\dots ,{\theta }_n),\mathbf{\epsilon }=\mathsf{diag}({\epsilon }_1,\dots ,{\epsilon }_n),\mathbf{\alpha }=\mathsf{diag}({\alpha }_1,\dots ,{\alpha }_n),\mathbf{\mu }=\mathsf{diag}({\mu }_1,\dots ,{\mu }_n),$$

$${\mathbf{\beta }}_{cc}=\mathsf{diag}({\beta }_{cc}^1,\dots ,{\beta }_{cc}^n),{\mathbf{\beta }}_{ch}=\mathsf{diag}({\beta }_{ch}^1,\dots ,{\beta }_{ch}^n).$$

Finally, for each compartment $X\in \{S_c,V_c,I_c,S_h,E_h,I_h,R_h\}$, the associated movement matrix is given by:

$${\mathcal{M}}^X=\left(\begin{array}{cccc}-\sum _{k=1}^n{m}_{k1}^X& m_{12}^X& \dots & m_{1n}^X\\ m_{21}^X& -\sum _{k=1}^n{m}_{k2}^X& \dots & m_{2n}^X\\ ⋮& & \ddots & ⋮\\ m_{n1}^X& m_{n2}^X& \dots & -\sum _{k=1}^n{m}_{kn}^X\end{array}\right).$$

(2)

Movement matrices of the form (2) possess many useful properties for the analysis of metapopulation systems [14,18], which we will use later in the analysis.

We denote ∘ as the Hadamard product (element-wise product).

The vector form of system (Section 2) is then written as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{S}}_c& ={\Lambda }_c-{\mathbf{\beta }}_{cc}{\mathbf{S}}_c{\displaystyle \frac{{\mathbf{I}}_c}{{\mathbf{N}}_c}}-{\mathbf{\nu }}_c{\mathbf{S}}_c-{\mathbf{\mu }}_c{\mathbf{S}}_c+{\mathcal{M}}^{S_c}{\mathbf{S}}_c,\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{V}}_c& ={\mathbf{\nu }}_c{\mathbf{S}}_c-{\mathbf{\mu }}_c{\mathbf{V}}_c+{\mathcal{M}}^{V_c}{\mathbf{V}}_c,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{I}}_c& ={\mathbf{\beta }}_{cc}{\mathbf{S}}_c{\displaystyle \frac{{\mathbf{I}}_c}{{\mathbf{N}}_c}}-({\mathbf{\mu }}_c+{\mathbf{\delta }}_c){\mathbf{I}}_c+{\mathcal{M}}^{I_c}{\mathbf{I}}_c,\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{S}}_h& ={\Lambda }_h-{\mathbf{\beta }}_{ch}{\displaystyle \frac{{\mathbf{I}}_c}{{\mathbf{N}}_c}}\circ {\mathbf{S}}_h-\mathbf{\mu }{\mathbf{S}}_h+{\mathcal{M}}^S{\mathbf{S}}_h,\hfill \end{array}$$

(3d)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{E}}_h& ={\mathbf{\beta }}_{ch}{\displaystyle \frac{{\mathbf{I}}_c}{{\mathbf{N}}_c}}\circ {\mathbf{S}}_h-(1-\mathbf{\theta })\mathbf{\gamma }{\mathbf{E}}_h-\mathbf{\mu }{\mathbf{E}}_h+{\mathcal{M}}^E{\mathbf{E}}_h,\hfill \end{array}$$

(3e)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{T}}_h& =\mathbf{\theta }{\mathbf{E}}_h\mathbf{\alpha }\mathbf{\theta }{\mathbf{T}}_h-\mathbf{\mu }{\mathbf{T}}_h+{\mathcal{M}}^T{\mathbf{T}}_h,\hfill \end{array}$$

(3f)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{I}}_h& =(1-\mathbf{\theta })\mathbf{\gamma }{\mathbf{E}}_h-({\mathbf{\delta }}_h+\mathbf{\mu }){\mathbf{I}}_h+{\mathcal{M}}^I{\mathbf{I}}_h,\hfill \end{array}$$

(3g)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{R}}_h& =\mathbf{\alpha }\mathbf{\theta }{\mathbf{T}}_h-\mathbf{\mu }{\mathbf{R}}_h+{\mathcal{M}}^T{\mathbf{R}}_h.\hfill \end{array}$$

(3h)

3.2. Positivity

Theorem 3.1. 

The components $S_h^i\left(t\right),S_c^i\left(t\right),E_h^i\left(t\right),I_h^i\left(t\right),I_c^i\left(t\right),V_c^i\left(t\right),T_h^i\left(t\right),R_h^i\left(t\right)$ are positive for all time t.

Proof. 

Equation (1e) implies ${\displaystyle \frac{d{S}_h^i}{dt}}+({\beta }_{ch}^i{\displaystyle \frac{I_c^i}{N_c^i}}+{\mu }_h^i)S_h^i={\Lambda }_h^i+m_h^S{S}_h^j-m_h^S{S}_h^i$. Its vector form gives ${\displaystyle \frac{d{\mathbf{S}}_h}{dt}}+({\mathbf{\beta }}_{ch}{\displaystyle \frac{{\mathbf{I}}_c}{{\mathbf{N}}_c}}+{\mathbf{\mu }}_h-{\mathcal{M}}_c^S){\mathbf{S}}_h={\Lambda }_h$, where ${\mathcal{M}}_c^S$ denotes the movement matrix of susceptible canines. The solution to this differential equation is

$${\mathbf{S}}_h\left(t\right)=e^{-{\int }_0^{t}f\left(\tau \right)d\tau }\left[{\mathbf{S}}_h\left(0\right)+\Lambda {\int }_0^t{e}^{-{\int }_0^{\tau }f\left(s\right)ds}\ge 0,\right]$$

with $g\left(t\right)={\mathbf{\beta }}_{ch}{\displaystyle \frac{{\mathbf{I}}_c\left(t\right)}{{\mathbf{N}}_c\left(t\right)}}+{\mathbf{\mu }}_h-{\mathcal{M}}_c^S$.

Similarly, the solution to Equation (1a) is given by

$${\mathbf{S}}_c\left(t\right)=e^{-{\int }_0^{t}g\left(\tau \right)d\tau }\left[{\mathbf{S}}_c\left(0\right)+\Lambda {\int }_0^t{e}^{-{\int }_0^{\tau }g\left(s\right)ds}\ge 0\right],$$

with $g\left(t\right)={\mathbf{\beta }}_{cc}{\displaystyle \frac{{\mathbf{I}}_c\left(t\right)}{{\mathbf{N}}_c\left(t\right)}}+{\mathbf{\nu }}_c+{\mathbf{\mu }}_c-{\mathcal{M}}_c^S$. □

Theorem 3.2. 

The total human $N_h^i$ and canine $N_c^i$ populations are bounded for all positive time t.

Proof. 

Define the total human population of country i: $N_h^i=S_h^i+E_h^i+T_h^i+I_h^i+R_h^i$.

Summing equations (1e)-(1i), we obtain:

$$\begin{array}{cc}\hfill N_h^{i^{\prime }}=& {\Lambda }_h^i-{\mu }_h^i{N}_h^i-{\delta }_h^i{I}_h^i\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \le & {\Lambda }_h^i-{\mu }_h^i{N}_h^i,\sin \mathrm{ce}{\delta }_h^i{I}_h^i\ge 0.\hfill \end{array}$$

(5)

This differential inequality is linear, and its solution is

$$\begin{array}{cc}\hfill N_h^i\left(t\right)\le & e^{-{\mu }_h^{i}t}\left[N_h\left(0\right)+{\Lambda }^i{\int }_0^t{e}^{-{\mu }_h^i\tau }d\tau \right]\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill =& e^{-{\mu }_h^{i}t}\left[N_h\left(0\right)-{\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}\right]+{\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}.\hfill \end{array}$$

(7)

Thus:

− If $N_h^i\left(t\right)\le {\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}$, then $0<N_h^i\left(t\right)\le {\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}$.

− If $N_h^i\left(t\right)>{\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}$, then $N_h^i\left(t\right)\le N_h^i\left(0\right),\forall t\ge 0$.

Therefore:

$$0\le N_h^i\left(t\right)\le max\left(N_h^i\left(0\right),{\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}}\right).$$

(8)

The same reasoning applies to the total canine population:

$$0\le N_c^i\left(t\right)\le max\left(N_c^i\left(0\right),{\displaystyle \frac{{\Lambda }_c^i}{{\mu }_c^i}}\right).$$

(9)

Thus, all human and canine populations are not only positive but also uniformly bounded in time. □

3.3. Isolated Case

This section considers the initial model without movement.

3.3.1. Disease-Free Equilibrium

The analysis of equilibrium points of a dynamical system allows for the identification of stationary states in which populations no longer vary over time.

Then, by setting the differential equations to zero and assuming the absence of infected or exposed individuals $(E_h^i=I_h^i=I_c^i=T_h^i=R_h^i=0)$, we obtain:

$$S_c^{i*}={\displaystyle \frac{{\Lambda }_c^i}{{\mu }_c^i+{\nu }_c^i}},V_c^{i*}={\displaystyle \frac{{\nu }_c^i{\Lambda }_c^i}{{\mu }_c^i({\mu }_c^i+{\nu }_c^i)}},S_h^{i*}={\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}},E_h^{i*}=I_h^{i*}=T_h^{i*}=R_h^{i*}=0$$

Thus, we have two disease-free equilibrium points:

$$E_c^{*}=\left({\displaystyle \frac{{\Lambda }_c^i}{{\mu }_c^i+{\nu }_c^i}},{\displaystyle \frac{{\nu }_c^i{\Lambda }_c^i}{{\mu }_c^i({\mu }_c^i+{\nu }_c^i)}},0\right)\mathrm{and}{E}_h^{*}=\left({\displaystyle \frac{{\Lambda }_h^i}{{\mu }_h^i}},0,0,0,0\right).$$

3.3.2. Basic Reproduction Number in an Isolated Location

To calculate the basic reproduction number,

We use the Next Generation Matrix method by van den Driessche and Watmough [19].

We consider the following infectious variables for country i, $(I_c^i,E_h^i)$.

The dynamics of these variables are decomposed into two functions:

• ${\mathcal{F}}_c$: rate of appearance of new infections in each compartment and
• ${\mathcal{V}}_c$: rate of transition between compartments (exits and entries due to progression or recovery),

which are given by:

$${\mathcal{F}}_c={\beta }_{cc}^i{S}_c^i{\displaystyle \frac{I_c^i}{N_c^i}}\mathrm{and}{\mathcal{V}}_c=({\mu }_c^i+{\delta }_c^i)I_c^i.$$

We compute the Jacobian of F and V at the disease-free equilibrium (DFE), i.e., when $I_c^i=E_h^i=I_h^i=0$.

The Jacobian of $F_c$ and $V_c$ at the DFE are:

$$F_c={\beta }_{cc}^i{\displaystyle \frac{S_c^{i*}}{N_c^{i*}}},\mathrm{and}{V}_c={\mu }_c^i+{\delta }_c^i\Rightarrow V_c^{-1}={\displaystyle \frac{1}{{\mu }_c^i+{\delta }_c^i}}$$

We define the next generation matrix $F_c{V}_c^{-1}$:

$$F_c{V}_c^{-1}={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)}}$$

The basic reproduction number ${\mathcal{R}}_0$ is the spectral radius (dominant eigenvalue) of $F_c{V}_c^{-1}$. Here, we have:

$${\mathcal{R}}_0^i={\mathcal{R}}_{0c}^i={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c}{({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)}}$$

(10)

Remark 3.1. 

This model assumes that humans do not retransmit rabies, which is biologically correct (rabies is not transmissible between humans). This is why $I_h^i$ is not a source of new infections, and humans do not fuel the epidemic reproduction dynamics.

3.3.3. Local Stability Analysis

Theorem 3.3. 

The disease-free equilibrium point is locally asymptotically stable when ${\mathcal{R}}_0^i<1$ and unstable otherwise.

Proof. 

The Jacobian matrix evaluated around the DFE is given by

$$J_{E_0}=\left[\begin{array}{ccc}-{\nu }_c^i-{\mu }_c^i& 0& -{\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{{\mu }_c^i+{\nu }_c^i}}\\ {\nu }_c^i& -{\mu }_c^i& 0\\ 0& 0& {\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{{\mu }_c^i+{\nu }_c^i}}-({\mu }_c^i+{\delta }_c^i)\end{array}\right]$$

The stability of the DFE depends on the sign of the eigenvalues of this matrix. We can see:

The 2×2 submatrix in the top left is triangular, so its eigenvalues are:

$${\lambda }_1=-{\nu }_c^i-{\mu }_c^i,{\lambda }_2=-{\nu }_c^i-{\mu }_c^i<0$$

The last eigenvalue (related to $I_c^i$ ) is:

$${\lambda }_3={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{{\mu }_c^i+{\nu }_c^i}}-({\mu }_c^i+{\delta }_c^i)=({\mu }_c^i+{\delta }_c^i)\left[{\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{({\mu }_c^i+{\delta }_c^i)({\mu }_c^i+{\delta }_c^i)}}-1\right]=({\mu }_c^i+{\delta }_c^i)\left[{\mathcal{R}}_0^i-1\right].$$

Then, if ${\mathcal{R}}_0^i<1$, ${\lambda }_3<0$. So all eigenvalues have negative real parts. Therefore, the disease-free equilibrium $E_0$ is locally asymptotically stable.

But, if ${\mathcal{R}}_0^i>1$, ${\lambda }_3>0$. And since one of the eigenvalues is positive, then the disease-free equilibrium $E_0$ is unstable. □

3.3.4. Global Stability of the Disease-Free Equilibrium

Theorem 3.4. 

The disease-free equilibrium $E_0$ is globally asymptotically stable if ${\mathcal{R}}_0^i<1$.

Proof. 

Let us verify the conditions of the Castillo-Chavez theorem.

• Let $z_1=(S_c^i,V_c^i)$ be the class of uninfected individuals, and $z_2=\left(I_c^i\right)$ the class of infected individuals.
  For $z_2=0$, the system becomes:

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d{S}_c^i}{dt}}& ={\Lambda }_c^i-{\beta }_{cc}{S}_c^i{\displaystyle \frac{I_c^i}{N_c^i}}-{\nu }_c^i{S}_c^i-{\mu }_c^i{S}_c^i,\hfill \end{array}$$

  (11)

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d{V}_c^i}{dt}}& ={\nu }_c^i{S}_c^i-{\mu }_c^i{V}_c^i.\hfill \end{array}$$

  (12)
  Solving these equations:

  $$S_c\left(t\right)=e^{-({\nu }_c^i+{\mu }_c^i)t}\left(S_H\left(0\right)-{\displaystyle \frac{{\Lambda }_H}{d}}+{\Lambda }_c^i{\int }_0^t{e}^{({\nu }_c^i+{\mu }_c^i)\tau }d\tau \right)$$

  i.e.,

  $$S_c\left(t\right)=e^{-({\nu }_c^i+{\mu }_c^i)t}\left(S_c^i\left(0\right)-{\displaystyle \frac{{\Lambda }_c^i}{{\nu }_c^i+{\mu }_c^i}}\right)+{\displaystyle \frac{{\Lambda }_c^i}{{\nu }_c^i+{\mu }_c^i}}$$

  and

  $$V_c\left(t\right)=e^{-{\mu }_c^{i}t}\left(V_c^i\left(0\right)-{\displaystyle \frac{{\nu }_c^i{\Lambda }_c^i}{{\mu }_c^i({\nu }_c^i+{\mu }_c^i)}}\right)+{\displaystyle \frac{{\nu }_c^i{\Lambda }_c^i}{{\mu }_c^i({\nu }_c^i+{\mu }_c^i)}}.$$
  Taking the limit as $t\to \infty$, we obtain:

  $$\underset{t\to \infty }{lim}{z}_1=\left({\displaystyle \frac{{\Lambda }_c^i}{{\nu }_c^i+{\mu }_c^i}},{\displaystyle \frac{{\nu }_c^i{\Lambda }_c^i}{{\mu }_c^i({\nu }_c^i+{\mu }_c^i)}}\right)=z_1^{★}.$$

  (13)
  So $z_1^{★}$ is globally stable when $z_2=0$.
• Now consider the system for infected individuals:

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d{I}_c^i}{dt}}& ={\beta }_{cc}{S}_c^i{\displaystyle \frac{I_c^i}{N_c^i}}-({\mu }_c^i+{\delta }_c^i)I_c^i=g_i.\hfill \end{array}$$

  (14)
  Let $A={\displaystyle \frac{\partial \left(g_i\right)}{\partial z_2}}\left(E_0\right)$, we have:

  $$A=\left(\begin{array}{c}{\beta }_{cc}^i{\displaystyle \frac{{\mu }_c^i}{{\mu }^i+{\nu }^i}}-{\delta }_c^i-{\mu }_c^i\end{array}\right).$$
  Then

  $$\begin{array}{cc}\hfill \widehat{G}& =A{z}_2-g_i\hfill \end{array}$$

  (15)

  $$\begin{array}{cc}\hfill & ={\beta }_{cc}^i{\displaystyle \frac{{\mu }_c^i}{{\mu }^i+{\nu }^i}}{I}_c^i-({\delta }_c^i+{\mu }_c^i)I_c^i-{\beta }_{cc}^i{\displaystyle \frac{S_c^i}{N_c^i}}+({\delta }_c^i+{\mu }_c^i)I_c^i\hfill \end{array}$$

  (16)

  $$\begin{array}{cc}\hfill & ={\beta }_{cc}^i{\displaystyle \frac{{\mu }_c^i}{{\mu }^i+{\nu }^i}}{I}_c^i-{\beta }_{cc}^i{\displaystyle \frac{S_c^i}{N_c^i}}{I}_c^i\hfill \end{array}$$

  (17)

  $$\begin{array}{cc}\hfill & ={\beta }_{cc}^i{I}_c^i\left({\displaystyle \frac{{\mu }_c^i}{{\mu }^i+{\nu }^i}}-{\displaystyle \frac{S_c^i}{N_c^i}}\right)\ge 0\hfill \end{array}$$

  (18)

  because from (13) and (9)

  $${\displaystyle \frac{S_c^i}{N_c^i}}\le {\displaystyle \frac{{\mu }_c^i}{{\mu }^i+{\nu }^i}}.$$

Thus, $E_0$ is globally asymptotically stable when ${\mathcal{R}}_0<1$. □

3.4. Analytical Calculation of the Critical Vaccination Threshold

If

$${\mathcal{R}}_0={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)}},$$

the condition ${\mathcal{R}}_0<1$ implies:

$$({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)>{\beta }_{cc}^i{\mu }_c^i\Rightarrow {\nu }_c^i>{\mu }_c^i\left({\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}-1\right).$$

We define:

$${\nu }_c^{\mathrm{crit}}=max\left(0,{\mu }_c^i\left({\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}-1\right)\right).$$

• If ${\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}\le 1$, then ${\nu }_c^{\mathrm{crit}}\le 0$: no vaccination effort required to ensure ${\mathcal{R}}_0<1$ (theoretical viewpoint).
• If ${\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}>1$, then ${\nu }_c^{\mathrm{crit}}$ gives the minimal vaccination rate to achieve in steady state to obtain $R_0<1$.

Practical remark: this threshold is a theoretical guide; in reality, one must account for actual coverage, vaccine efficacies, parameter uncertainties, and logistics.

Corollary 3.1.(Critical canine vaccination threshold) From the expression

$${\mathcal{R}}_0={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)}},$$

the condition ${\mathcal{R}}_0<1$ provides a minimal threshold for the vaccination rate ${\nu }_c^i$, denoted ${\nu }_c^{\mathrm{crit}}$, such that if ${\nu }_c^i>{\nu }_c^{\mathrm{crit}}$ then ${\mathcal{R}}_0<1$.

Explicit formula:

$${\nu }_c^{\mathrm{crit}}=max\left\{0,{\mu }_c^i\left({\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}-1\right)\right\}.$$

Proof. 

Condition on ${\nu }_c^i$ for ${\mathcal{R}}_0<1$. We start from

$$R_0={\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)}}<1.$$

-

Isolate ${\nu }_c^i$: Multiply by $({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i)>0$:

$${\beta }_{cc}^i{\mu }_c^i<({\mu }_c^i+{\nu }_c^i)({\mu }_c^i+{\delta }_c^i).$$

Divide by ${\mu }_c^i+{\delta }_c^i>0$:

$${\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{{\mu }_c^i+{\delta }_c^i}}<{\mu }_c^i+{\nu }_c^i.$$

Isolate ${\nu }_c^i$:

$${\nu }_c^i>{\displaystyle \frac{{\beta }_{cc}^i{\mu }_c^i}{{\mu }_c^i+{\delta }_c^i}}-{\mu }_c^i={\mu }_c^i\left({\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}-1\right).$$

-

Negative case: if ${\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}\le 1$, the right-hand side is $\le 0$, so a rate ${\nu }_c=0$ satisfies ${\mathcal{R}}_0<1$. Otherwise, a positive minimal vaccination effort is required.

-

Compact formulation:

$${\nu }_c^{\mathrm{crit}}=max\left\{0,{\mu }_c^i\left({\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}-1\right)\right\}.$$

□

3.4.1. Equivalent Expression for Vaccination Coverage p

We express p as the proportion of vaccinated dogs at equilibrium:

$$p={\displaystyle \frac{V_c^{i*}}{N_c^{i*}}}={\displaystyle \frac{{\nu }_c^i}{{\mu }_c^i+{\nu }_c^i}}.$$

Solving for ${\nu }_c^i$:

$${\nu }_c^i={\displaystyle \frac{p{\mu }_c^i}{1-p}}.$$

Substitute into ${\mathcal{R}}_0$:

$${\mathcal{R}}_0={\displaystyle \frac{{\beta }_{cc}^i}{{\mu }_c^i+{\delta }_c^i}}·{\displaystyle \frac{{\mu }_c^i}{1-p}}.$$

The condition ${\mathcal{R}}_0<1$ becomes:

$$1-p<{\displaystyle \frac{{\mu }_c^i+{\delta }_c^i}{{\beta }_{cc}^i}}\iff p>1-{\displaystyle \frac{{\mu }_c^i+{\delta }_c^i}{{\beta }_{cc}^i}}.$$

Critical coverage:

$$p^{\mathrm{crit}}=max\left\{0,1-{\displaystyle \frac{{\mu }_c^i+{\delta }_c^i}{{\beta }_{cc}^i}}\right\}.$$

We verify that this corresponds to ${\nu }_c^{\mathrm{crit}}$:

$${\nu }_c^{\mathrm{crit}}={\displaystyle \frac{p^{\mathrm{crit}}{\mu }_c^i}{1-p^{\mathrm{crit}}}}.$$

3.4.2. Generalization: Imperfect Vaccine

If the vaccine has an efficacy $0\le \epsilon \le 1$, the protected fraction is $\epsilon p$ and the remaining fraction is $1-\epsilon p$.

The condition becomes:

$$p^{\mathrm{crit}}=max\left\{0,1-{\displaystyle \frac{{\mu }_c^i+{\delta }_c^i}{{\epsilon }^i{\beta }_{cc}^i}}\right\}.$$

3.5. Disease-Free Equilibrium

Then, by setting the differential equations of system Section 3.1 to zero and assuming the absence of infected or exposed individuals $({\mathbf{E}}_h={\mathbf{I}}_h={\mathbf{I}}_c={\mathbf{T}}_h={\mathbf{R}}_h=\mathbf{0})$, and noting that the matrix $\mathbf{\mu }-\mathcal{M}$ is non-singular and thus invertible, we obtain:

$${\mathbf{S}}_c^{*}={\Lambda }_c{({\mathbf{\mu }}_c+{\mathbf{\nu }}_c-{\mathcal{M}}^{S_c})}^{-1},{\mathbf{V}}_c^{*}={\mathbf{\nu }}_c{\Lambda }_c{\left({\mathbf{\mu }}_c({\mathbf{\mu }}_c+{\mathbf{\nu }}_c)\right)}^{-1},{\mathbf{S}}_h^{*}={\Lambda }_h{\mathbf{\mu }}_h^{-1},{\mathbf{E}}_h^{*}={\mathbf{I}}_h^{*}={\mathbf{T}}_h^{*}={\mathbf{R}}_h^{*}=0$$

Thus, we have two disease-free equilibrium points:

$${\mathbf{E}}_c^{*}=\left({\Lambda }_c{({\mathbf{\mu }}_c+{\mathbf{\nu }}_c-{\mathcal{M}}^{S_c})}^{-1},{\mathbf{\nu }}_c{\Lambda }_c{\left({\mathbf{\mu }}_c({\mathbf{\mu }}_c+{\mathbf{\nu }}_c)\right)}^{-1},\mathbf{0}\right),{\mathbf{E}}_h^{*}=\left({\Lambda }_h{\mathbf{\mu }}_h^{-1},\mathbf{0},\mathbf{0},\mathbf{0},\mathbf{0}\right).$$

(19)

3.6. Basic Reproduction Number for the Global System

Using the same approach as in Section 3.3.2, we obtain:

$$F{V}^{-1}={\mathbf{\beta }}_{cc}{\mathbf{\mu }}_c{({\mathbf{\mu }}_c+{\mathbf{\nu }}_c-{\mathcal{M}}^{S_c})}^{-1}{({\mathbf{\mu }}_c+{\mathbf{\delta }}_c-{\mathcal{M}}^{I_c})}^{-1}.$$

Then, the basic reproduction number is:

$${\mathcal{R}}_0^{global}=\rho \left({\mathbf{\beta }}_{cc}{\mathbf{\mu }}_c{({\mathbf{\mu }}_c+{\mathbf{\nu }}_c-{\mathcal{M}}^{S_c})}^{-1}{({\mathbf{\mu }}_c+{\mathbf{\delta }}_c-{\mathcal{M}}^{I_c})}^{-1}\right).$$

(20)

3.7. Local Stability

First, from [19], we have the following result.

Theorem 3.5. 

$E_0^{★}$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Proof. 

We need to verify that assumptions A1-A5 of [19] are satisfied. Assumptions A1-A4 follow from the procedure used above to derive F and V in the calculation of ${\mathcal{R}}_0$. Therefore, all we need to verify is that the disease-free system at the DFE is locally asymptotically stable. In the absence of disease, (Section 3.1) is a linear system

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}{\mathbf{S}}_c& ={\Lambda }_c+({\mathcal{M}}^{S_c}-{\mathbf{\nu }}_c-{\mathbf{\mu }}_c){\mathbf{S}}_c+{\mathbf{\epsilon }}_c{\mathbf{R}}_c\hfill \\ \hfill {\displaystyle \frac{d}{dt}}{\mathbf{V}}_c& ={\mathbf{\nu }}_c{\mathbf{S}}_c+({\mathcal{M}}^{V_c}-{\mathbf{\mu }}_c){\mathbf{V}}_c.\hfill \end{array}$$

This is exactly (Section 3.5), so we know it has a unique equilibrium, the disease-free equilibrium (19). The Jacobian matrix of the system at any point takes the form

$$\left(\begin{array}{cc}{\mathcal{M}}^{S_c}-{\mathbf{\nu }}_c-\mathbf{\mu }& 0\\ {\mathbf{\nu }}_c& {\mathcal{M}}^{V_c}-{\mathbf{\mu }}_c\end{array}\right).$$

We saw in Section 3.1 that this matrix is invertible. According to [18], the spectral abscissa of this matrix is negative, so the disease-free equilibrium is always (locally) asymptotically stable. The result then follows. □

4. Sensitivity Analysis

We use the data from Table 3 for our numerical simulations.

Figure 2. PRCC of ${\mathcal{R}}_0^i$.

Figure 2. PRCC of ${\mathcal{R}}_0^i$.

Interpretation: This analysis reveals which parameters to target for effective interventions:

Parameters with high positive PRCC represent priority targets for reduction measures (e.g., decreasing transmission rates). Parameters with high negative PRCC represent intervention opportunities (e.g., increasing vaccination or treatment rates).

5. Numerical Simulations

To visualize the temporal dynamics of the model and understand the impact of the basic reproduction number (${\mathcal{R}}_0$), we performed numerical simulations for canine and human populations. Two contrasting scenarios were studied: one where ${\mathcal{R}}_0<1$, indicating the epidemic is theoretically controlled, and another where ${\mathcal{R}}_0>1$, a scenario in which the disease can spread within the population.

6. Discussion

Figure 3 shows the dynamics for ${\mathcal{R}}_0^c=1.827$. A rapid epidemic is observed: the infected compartment ($I_c$) experiences rapid exponential growth, peaking at a high level where a large portion of the population is simultaneously infected. Subsequently, the number of infected decreases as the pool of susceptibles ($S_c$) is depleted, either by infection or vaccination. This curve is characteristic of an uncontrolled epidemic.

Conversely, Figure 4 illustrates the scenario where ${\mathcal{R}}_0^c=0.818$. Here, the number of infected decreases monotonically from the start, without an epidemic peak. Each infected individual generates on average less than one new infection, leading to the natural and rapid extinction of the disease. The vaccinated population ($V_c$) remains significant, contributing to herd immunity. This result validates the theoretical threshold: if ${\mathcal{R}}_0<1$, the disease cannot be maintained in the population.

The impact of the two canine scenarios on public health is striking. Figure 5 shows that an uncontrolled canine epidemic (${\mathcal{R}}_0^c>1$) leads to a significant risk for humans (${\mathcal{R}}_0^h=1.089>1$), resulting in an increase in cases of exposure ($E_h$), infection ($I_h$), and requiring treatment ($T_h$).

In contrast, Figure 6 demonstrates the beneficial effect of controlling rabies at its animal source. When ${\mathcal{R}}_0^c<1$, transmission from dog to human is interrupted, as evidenced by the value ${\mathcal{R}}_0^h=0.096\ll 1$. The curves for the exposed, infected, and treated compartments remain consistently at negligible, or even zero, levels. The human population remains mostly, if not entirely, in the susceptible ($S_h$) or immune ($R_h$) compartment, confirming that the best strategy to protect humans is to massively vaccinate dogs.

These results highlight the critical relationship between disease control in the animal reservoir (the dog) and the risk to human health. They emphasize the importance of maintaining sufficient canine vaccination coverage to ensure that ${\mathcal{R}}_0^c$ remains sustainably below 1, in order to prevent both canine epidemics and human rabies cases.

7. Conclusion

This study developed and analyzed a metapopulation-type mathematical model for the dynamics of canine and human rabies in a cross-border context. The theoretical analysis allowed for the determination of the basic reproduction number, ${\mathcal{R}}_0$, whose threshold value of 1 governs the stability of the disease-free equilibrium. The results rigorously demonstrate that this equilibrium is locally and globally asymptotically stable when ${\mathcal{R}}_0<1$, guaranteeing epidemic extinction, and unstable otherwise, leading to endemic persistence. Sensitivity analysis identified the parameters most influential on $\mathcal{R}0$, such as the inter-canine transmission rate ($\beta cc$) and the vaccination rate (${\nu }_c$), thus providing priority targets for intervention strategies.

Numerical simulations concretely illustrated the crucial impact of canine vaccination. An insufficient vaccination rate, leading to ${\mathcal{R}}_0>1$, results in an explosive epidemic in the canine population, which immediately translates into an increased risk of transmission to humans. Conversely, vaccination coverage exceeding the critical threshold, ensuring ${\mathcal{R}}_0<1$, effectively extinguishes transmission in both dogs and humans. These results underscore the imperative necessity of cross-border coordination and massive, sustained dog vaccination campaigns as a central strategy to achieve the global goal of "Zero human deaths from dog-mediated rabies by 2030".