Stability analysis and numerical simulation of rabies spread model with delay effects

1.   Introduction

Rabies is a viral disease that attacks the nervous system and can be fatal once symptoms appear ^[1]. Rabies is usually spread through bites or scratches from wild animals. It can also spread through direct contact if wounds or mucosa (e.g., mouth and eyes) get contaminated with the saliva of the infected animals ^[1,2]. Rabies causes around 59,000 deaths worldwide each year, where cases mostly occur in Asia and Africa ^[1,3].

Dogs are the main carrier of rabies to humans. It can also be found in other animals, like raccoons, foxes, and cats ^[1,4]. Rabies usually needs laboratory testing to be diagnosed even though rabid animals may act strangely, such as being aggressive and tending to bite people or other animals, along with excessive drooling ^[5]. After biting, rabies needs an incubation period that can take up to 2–3 months and may vary from weeks to years, depending on certain factors ^[1]. There are no approved methods of treatment for rabies once symptoms appear. If someone has been exposed to rabies, either by bites or scratches, they should contact a healthcare professional as soon as possible ^[4]. In addition, vaccination of dogs can be chosen as one way to prevent rabies at the source. Vaccinations can also be done to people after or before rabies exposure ^[1]. General awareness and education on preventing spread of rabies, especially from domestic dogs, are also important to reduce incidence and deaths caused by rabies ^[1,6].

Mathematical modeling for epidemiology has given insights on the dynamic of rabies spread and measures on how to possibly reduce or eradicate the disease. Some researchers have studied mathematically the spread of rabies, especially within the Asian region. For example, Zhang et al. ^[7] has developed a model to analyze the rabies spread in China. The model observed the spread of rabies among dog and human populations, and the research also compared the effects of culling and vaccination on dogs to the spread of rabies. They found that the most effective methods to control the spread of rabies was by controlling birth of dogs and increasing vaccination coverage for dogs. Chen et al. ^[3], based on the constructed rabies spread model among human and dog populations, found that the immigration of dogs may create a disease endemic even if the disease dies out. They stated that the migration of dogs should be better monitored and always under surveillance. Tohma et al. ^[8] observed inter-island transmission of rabies in the Philippines. They stated that to control the spread of rabies, it was important to acknowledge that the inter-island transmission can occur because rabies can become endemic when the virus is introduced to an island that was previously rabies-free. Continuous rabies control programs in the Philippines, such as controlling transportation of dogs, should be implemented to prevent rabies spread. Huang and Li ^[9], also studying the case of rabies spread in China, developed the model of rabies spread among humans and domesticated and wild dog populations. The model was analyzed by fitting data obtained from literatures and officials in China, as well as studying the effectivity of different suppression methods. It was observed that relative suppression measures were including controlling birth of domesticated and wild dogs along with increasing immunity among domesticated dogs. Pantha et al. ^[10], who observed rabies among the populations of human, dog, and jackal, showed that jackals and dogs both played important roles in the spread of rabies in Nepal, even though dogs played a greater role in the spread. They also observed that interspecies transmission might occur between the jackal and dog populations. Some researchers also highlighted importance on implementation of spread control on dog populations. For instance, Asamoah et al. ^[11] obtained that, using optimal control, the deaths could be eradicated through mass vaccination of susceptible dogs and continuous use of pre and postexposure prophylaxis on humans. Carroll et al. ^[12] discovered that mass vaccination and population control for dogs are some of the effective measures to control the spread. Bornaa et al. ^[13] found that the efforts to control the spread of rabies should be focused more on the dogs than humans, which they stated that the disease can be controlled through reducing contact with infectious dogs, increasing vaccination, screening of recruited dogs, and culling of infectious dogs. Similarly, Renald et al. ^[14] found that the efforts of controlling the spread of rabies should be focused more on stray dogs.

Another way rabies transmission can be represented mathematically is with a system of delay differential equations. Time delay is a natural property of a system where the response of an action has a delayed effect ^[15]. A delay differential equation presents the population growth at time t that depends on the population at the past time ^[16]. There have been researchers that analyzed the model of delay differential equations, such as Song and Xu ^[17], who observed the existence of Hopf bifurcation in a model of a two neural-network system with multiple delays, and Kasbawati et al. ^[18] observed the pathogen-immune system interaction with delay effects. An example of research of rabies spread modeling by incorporating delay is by Abdulmajid and Hassan ^[2], where they formulated a SIV (susceptible-infected-vaccinated) epidemic model with delay to assess the effects of controls and time delay of incubation period on the transmission of rabies within human and dog populations. It was found that an increase in dog vaccination and decrease in puppies born were effective measures in eradicating rabies. In this research, we also formulated and analyzed a rabies spread model among human and dog populations as a system of delay differential equations with delay effects. In addition to the incubation period, vaccine activation time can be considered as a time delay that affects the rabies spread model. Our major research objectives are to analyze the stability of the delay system and analyze the effects of two delays on the stability of the system, as well as finding the best possible methods to reduce or eradicate the spread of rabies. The proposed model is believed to complement previous research on the complex nonlinear dynamics of rabies transmission.

2.   Methodology

The constructed model consists of two populations, human and dog populations, with both living within the same environment. The human population consists of three subpopulations: A susceptible human population (${S}_{H}$), infected human population (${I}_{H}$), and vaccinated human population (${V}_{H}$). The dog population also consists of three subpopulations: A susceptible dog population (${S}_{D}$), infected dog population (${I}_{D}$), and vaccinated dog population (${V}_{D}$). The model is built based on the following major assumptions:

1) No migrations happen to and from the system. The population size depends only on birth and death.

2) The transmission of rabies is only from rabid dogs. No transmission of rabies between susceptible and infected humans, nor is there any transmission between susceptible dogs and infected humans.

3) Vaccination is assumed to give a long-time immunity. Therefore, there is no immunity loss for vaccinated populations.

4) Exposure to the infected dogs does not instantly breed new infections. Instead, there is an incubation period.

5) Vaccination does not instantly give immunity. There is a period required for the vaccine to be able to give immunity.

6) Within the incubation and vaccination period, there is a chance that natural death may happen. Here, natural death is assumed to be an event with random arrival.

The rabies spread dynamic between human and dog populations is fully illustrated in Figure 1. As shown in Figure 1, the susceptible human population is increased by average human birth per year (${A}_{H}$). It is decreased by natural death at rate ${m}_{H}$, infection by contact with infected dogs at rate ${\beta }_{HD}{S}_{H}{{I}_{D}}_{{\tau }_{1}}{e}^{-{m}_{D}{\tau }_{1}}$, and vaccination at rate ${k}_{H}{{S}_{H}}_{{\tau }_{2}}{e}^{-{m}_{H}{\tau }_{2}}$. Here, ${\tau }_{1}$ is the incubation delay and ${\tau }_{2}$ is the vaccination delay. Moreover, ${e}^{-{m}_{D}{\tau }_{1}}$ is the probability of rabid dogs surviving natural death within the time period $[0, {\tau }_{1}]$, ${e}^{-{m}_{H}{\tau }_{2}}$ is the probability of susceptible humans surviving natural death within the time period $[0, {\tau }_{2}]$, ${\beta }_{HD}$ is the transmission rate between susceptible humans and infected dogs, and ${k}_{H}$ is the vaccination rate of humans. The infected human population is increased through transmission that occurs due to contact between susceptible humans and infected dogs at rate ${\beta }_{HD}{S}_{H}{{I}_{D}}_{{\tau }_{1}}{e}^{-{m}_{D}{\tau }_{1}}$, and decreased by natural and rabies-related deaths at rate ${m}_{H}$ and ${\mu }_{H}$, respectively. The vaccinated human population is increased by vaccination of susceptible humans at rate ${k}_{H}{{S}_{H}}_{{\tau }_{2}}{e}^{-{m}_{H}{\tau }_{2}}$ and decreased by natural death at rate ${m}_{H}$.

The susceptible dog population is increased by average dog birth per year (${A}_{D}$). It is decreased by natural death at rate ${m}_{D}$, infection by contact with infected dogs at rate ${\beta }_{DD}{S}_{D}{{I}_{D}}_{{\tau }_{1}}{e}^{-{m}_{D}{\tau }_{1}}$, and vaccination at rate ${k}_{D}{{S}_{D}}_{{\tau }_{2}}{e}^{-{m}_{D}{\tau }_{2}}$. Here, the ${e}^{-{m}_{D}{\tau }_{2}}$ is the probability of susceptible dogs surviving natural death within the time period $[0, {\tau }_{2}]$, ${\beta }_{DD}$ is the transmission rate between susceptible dogs and infected dogs, and ${k}_{D}$ is the vaccination rate of dogs. The infected dog population is increased through infection by contact between susceptible and infected dogs at rate ${\beta }_{DD}{S}_{D}{{I}_{D}}_{{\tau }_{1}}{e}^{-{m}_{D}{\tau }_{1}}$, and decreased by natural and rabies-related deaths at rate ${m}_{D}$ and ${\mu }_{D}$, respectively. The vaccinated dog population is increased by vaccination of susceptible dogs at rate ${k}_{D}{{S}_{D}}_{{\tau }_{2}}{e}^{-{m}_{D}{\tau }_{2}}$ and decreased by natural death at rate ${m}_{D}$. Table 1 summarizes the description of all variables and model parameters

Based on the transmission diagram in Figure 1 and the assumptions described earlier, the rate of changes of all subpopulations can be written as a system of delay differential equations as follows:

with initial values of ${I}_{H}\left(t\right)$, ${V}_{H}\left(t\right)$, and ${V}_{D}\left(t\right)$ are given as ${I}_{H}\left(0\right)\ge 0$, ${V}_{H}\left(0\right)\ge 0$, and ${V}_{D}\left(0\right)\ge 0$, respectively, and historical functions of ${S}_{H}\left(t\right)$, ${S}_{D}\left(t\right)$, and ${I}_{D}\left(t\right)$ are given as ${\psi }_{1}\left(t\right)\ge 0$, ${\psi }_{2}\left(t\right)\ge 0$, and ${\psi }_{3}\left(t\right)\ge 0$, respectively. The functions of ${\psi }_{1}\left(t\right)$ and ${\psi }_{2}\left(t\right)$ are continuous functions within the interval $\left[-{\tau }_{2}, 0\right]$, and ${\psi }_{3}\left(t\right)$ is a continuous function within the interval $\left[-{\tau }_{1}, 0\right]$.

3.   Model analysis

3.1.   Basic properties of the model

Since the model (1)–(6) is an epidemiology model, it is important that the solution exists and bounded within a feasible region. This section will be exploring on the positivity and the boundedness of the solutions of (1)–(6). Suppose the model (1)–(6) is divided into two regions; ${\mathrm{\Omega }}_{H} = \left\{{S}_{H}, {I}_{H}, {V}_{H}\right\}\in {\mathbb{R}}_{+}^{3}$ and ${\mathrm{\Omega }}_{D} = \left\{{S}_{D}, {I}_{D}, {V}_{D}\right\}\in {\mathbb{R}}_{+}^{3}$. Thus, $\mathrm{\Omega } = {\mathrm{\Omega }}_{H}\times {\mathrm{\Omega }}_{D}$.

3.1.1.   Positivity of the solution

Suppose that ${S}_{H}$ is a decreasing function. Since ${A}_{H}\ge 0$, then (1) yields

Similarly, suppose that ${S}_{D}$ is a decreasing function. Since ${A}_{D}\ge 0$, then (4) yields

Therefore, ${S}_{H}$ and ${S}_{D}$ are positive for all $t > 0$. As ${S}_{H} > 0$ and ${S}_{D} > 0$, then from (3) and (6), it is obtained that ${\dot{V}}_{H}\ge -{m}_{H}{V}_{H}$ and ${\dot{V}}_{D}\ge -{m}_{D}{V}_{D}$, which yield ${V}_{H}\ge {V}_{H}\left(0\right){e}^{{m}_{H}t} > 0$ and ${V}_{D}\ge {V}_{D}\left(0\right){e}^{{m}_{D}t} > 0$, respectively. Therefore, ${V}_{H}$ and ${V}_{D}$ are also positive for all $t > 0$.

Furthermore, for (2) and (6), it is analyzed by using the method in ^[2]. First, it should be noted that ${I}_{H}\left(\theta \right)$ and ${I}_{D}\left(\theta \right)$ are positive for $\theta \in \left[-{\tau }_{1}, 0\right]$. Suppose that, on the contrary, there exists ${t}^{*} > 0$ such that ${I}_{H}\left({t}^{*}\right) = 0$ and ${I}_{D}\left({t}^{*}\right) = 0$, $\dot{{I}_{H}}\left({t}^{*}\right)\le 0$ and $\dot{{I}_{D}}\left({t}^{*}\right)\le 0$, and ${S}_{H} > 0$, ${I}_{H} > 0$, ${V}_{H} > 0$, ${S}_{D} > 0$, ${I}_{D} > 0$, and ${V}_{D} > 0$, for $t\in [0, {t}^{*})$. From (2) and (6), it yields

which contradicts to the initial assumption. Therefore, ${I}_{H}$ and ${I}_{D}$ are also positive for all $t > 0$.

3.1.2.   Boundedness of the solution

Now, suppose that ${N}_{H} = {S}_{H}+{I}_{H}+{V}_{H}$ and ${N}_{D} = {S}_{D}+{I}_{D}+{V}_{D}$, where ${N}_{H}$ is the total human population and ${N}_{D}$ is the total dog population. Adding (1)–(3) and (4)–(6) results in

As performed by ^[2,11,13], the solution of (7) and (8) can be written as

It can be seen from (9) that ${N}_{H}$ tends to $\frac{{A}_{H}}{{m}_{H}}$ as $t\to \infty$. Similarly, from (10), it can be seen that ${N}_{D}$ tends to $\frac{{A}_{D}}{{m}_{D}}$ as $t\to \infty$. Therefore, the feasible region of the solutions of human and dog populations are given as ${\mathrm{\Omega }}_{H} = \left\{\left({S}_{H}, {I}_{H}, {V}_{H}\right)\in {\mathbb{R}}_{+}^{3}, {N}_{H}\le \frac{{A}_{H}}{{m}_{H}}\right\}$ and ${\mathrm{\Omega }}_{D} = \left\{\left({S}_{D}, {I}_{D}, {V}_{D}\right)\in {\mathbb{R}}_{+}^{3}, {N}_{D}\le \frac{{A}_{D}}{{m}_{D}}\right\}$. This gives the feasible region of the solutions of the system (1)–(6) as $\mathrm{\Omega } = \left\{\left({S}_{H}, {I}_{H}, {V}_{H}, {S}_{D}, {I}_{D}, {V}_{D}\right)\in {\mathbb{R}}_{+}^{6}, {N}_{H}\le \frac{{A}_{H}}{{m}_{H}}, {N}_{D}\le \frac{{A}_{D}}{{m}_{D}}\right\}$. Therefore, the model (1)–(6) is mathematically well-posed and epidemiologically meaningful.

3.2.   Existence of equilibrium solutions

Suppose ${\boldsymbol{X}}^{*} = \left({S}_{H}^{*}, {I}_{H}^{*}, {V}_{H}^{*}, {S}_{D}^{*}, {I}_{D}^{*}, {V}_{D}^{*}\right)\in {\mathbb{R}}_{6}^{+}$ is an arbitrary equilibrium of the system on (1)–(6). At the equilibrium, it applies that ${S}_{H} = {{S}_{H}}_{{\tau }_{2}} = {S}_{H}^{*},$ ${I}_{H} = {I}_{H}^{*},$ ${V}_{H} = {V}_{H}^{*},$ ${S}_{D} = {{S}_{D}}_{{\tau }_{2}} = {S}_{D}^{*},$ ${I}_{D} = {{I}_{D}}_{{\tau }_{1}} = {I}_{D}^{*},$ ${V}_{D} = {V}_{D}^{*},$ and ${\dot{S}}_{H} = {\dot{I}}_{H} = {\dot{V}}_{H} = {\dot{S}}_{D} = {\dot{I}}_{D} = {\dot{V}}_{D} = 0.$ Applying these conditions to (1)–(6) results in

3.2.1.   Disease-free equilibrium

From (15), it can be seen that ${I}_{D}^{*} = 0$ or ${I}_{D}^{*}\ne 0$, which yields two different equilibria. For the case of ${I}_{D}^{*} = 0$, after substituting and solving (11)–(16), it gives the disease-free equilibrium given as

As can be seen, the non-existence of the infected dog population within the equilibrium also causes the non-existence of infected human population within the equilibrium. This suggests that control of spread of rabies in the dog population is sufficient to suppress the spread of rabies for both human and dog populations.

3.2.2.   Basic reproduction number

Basic reproduction number $\left({R}_{0}\right)$ is the number of secondary infections caused by a primary infection when it is introduced to a susceptible population within the infectious period of the primary infection ^[19]. The method of finding ${R}_{0}$ here is by using next generation method, as given by ^[20]. First, taking the infectious compartments from (2) and (5) by letting ${f}_{1} = {\dot{I}}_{H}$ and ${f}_{2} = {\dot{I}}_{D}$ gives

The next generation matrix is given by ${K}_{L} = -T{\Sigma }^{-1}$, where the matrix $T$ describes the events of producing new infections and the matrix $\Sigma$ describes the events of infectious population transitions. Matrices $T$ and $\Sigma$ are constructed such that (17) can be written as $\dot{\boldsymbol{F}} = \left(T+\Sigma \right)\boldsymbol{F},$ where $\boldsymbol{F} = \left({I}_{H}, {I}_{D}\right)$. By taking the derivatives of ${f}_{1}$ and ${f}_{2}$ with respect to ${I}_{H}$ and ${I}_{D}$ and evaluating the derivatives around ${E}_{0}$, the matrices $T$ and $\Sigma$ are given as

From (18) and (19), the next generation matrix is given as

Therefore, by taking the spectral radius of ${K}_{L}$, the basic reproduction number is given as

From (21), it can be seen the value of ${R}_{0}$ depends only on dog parameters, as well as both delays ${\tau }_{1}$ and ${\tau }_{2}$. Therefore, the efforts to control or eradicate the spread of rabies must be focused on the dog population. This supports the statement within Section 3.2.1 that the control of rabies spread within dog population is sufficient to control the spread within human population as well.

3.2.3.   Endemic equilibrium

For the case of ${I}_{D}^{*}\ne 0$, after substituting and solving (11)–(16), it gives the endemic equilibrium given as ${E}_{\varepsilon } = \left({S}_{H}^{\varepsilon }, {I}_{H}^{\varepsilon }, {V}_{H}^{\varepsilon }, {S}_{D}^{\varepsilon }, {I}_{D}^{\varepsilon }, {V}_{D}^{\varepsilon }\right)$, with

where the endemic equilibrium exists for ${R}_{0} > 1$.

3.3.   Local stability of equilibria

Analyzing local stability can be done by using linearization method. Suppose $\boldsymbol{X} = \left({S}_{H}\left(t\right), {I}_{H}\left(t\right), {V}_{H}\left(t\right), {S}_{D}\left(t\right), {I}_{D}\left(t\right), {V}_{D}\left(t\right)\right)$ and ${\boldsymbol{X}}_{{\tau }_{i}} = \left({{S}_{H}}_{{\tau }_{i}}, {{I}_{H}}_{{\tau }_{i}}, {{V}_{H}}_{{\tau }_{i}}, {{S}_{D}}_{{\tau }_{i}}, {{I}_{D}}_{{\tau }_{i}}, {{V}_{D}}_{{\tau }_{i}}\right)$ for $i = 1, 2$. Then, by using Taylor series expansion around ${\boldsymbol{X}}^{*}$ on (1) – (6), it results in

where ${J}_{0}$ is the Jacobian matrix with respect to $\boldsymbol{X}$, ${J}_{1}$ is the Jacobian matrix with respect to ${\boldsymbol{X}}_{{\tau }_{1}}$, and ${J}_{2}$ is the Jacobian matrix with respect to ${\boldsymbol{X}}_{{\tau }_{2}}$. The matrices ${J}_{0}, {J}_{1}, {J}_{2}$ are described as follows

with

Since the system is now in a linearized form, therefore the solution must in the form $\boldsymbol{X}\left(t\right) = \boldsymbol{C}{e}^{\lambda t}$, where $\boldsymbol{C}$ is a non-zero vector. From (22), the characteristic equation is given by

Applying ${J}_{0}, {J}_{1},$ and ${J}_{2}$ to (23) yields

with

Solving (24) results in

3.3.1.   Local stability of disease-free equilibrium

By substituting ${E}_{0}$ to (25), the characteristic equation becomes

From (26), some of the eigenvalues obtained are ${\lambda }_{1} = -\left({m}_{H}+{\mu }_{H}\right)$, ${\lambda }_{2} = -{m}_{H}$, and ${\lambda }_{3} = -{m}_{D}$, which are all real negatives. The other eigenvalues are given by the following equations:

which will be analyzed by considering several cases of ${\tau }_{1}$ and ${\tau }_{2}$.

Suppose ${\tau }_{1} = 0$ and ${\tau }_{2} = 0$, then the other eigenvalues are given as ${\lambda }_{4} = \left({m}_{D}+{\mu }_{D}\right)\left({R}_{0}-1\right)$, ${\lambda }_{5} = -\left({m}_{H}+{k}_{H}\right)$, and ${\lambda }_{6} = -\left({m}_{D}+{k}_{D}\right)$, which are all real negative eigenvalues if ${R}_{0} < 1$. Therefore, for ${\tau }_{1} = {\tau }_{2} = 0$, the disease-free equilibrium is locally asymptotically stable if ${R}_{0} < 1$.

Suppose ${\tau }_{1} > 0$ and ${\tau }_{2} = 0$. From (28) and (29), it yields ${\lambda }_{4} = -\left({m}_{H}+{k}_{H}\right)$ and ${\lambda }_{5} = -\left({m}_{D}+{k}_{D}\right)$, which are all real negative eigenvalues. For (27), it is assumed that the system experiences stability switch by the existence of a pair of imaginary eigenvalues $\lambda = \pm i\omega$, with $\omega \in {\mathbb{R}}^{+}$. Without losing generality, substituting $\lambda = i\omega$ to (27), and separating real and imaginary parts results in

Solving (30) and (31) by squaring and adding gives

or

In terms of ${R}_{0}$, we have $\omega = \left({m}_{D}+{\mu }_{D}\right)\sqrt{\left({R}_{0}-1\right)\left({R}_{0}+1\right)}$, where there are no $\omega \in {\mathbb{R}}^{+}$ if ${R}_{0} < 1$. Therefore, the roots of (27) must be real. It is left to verify that the roots of (27) are real and negative. Suppose that the left-hand side of (27) can be written as $g\left(\lambda \right) = \lambda +\left({m}_{D}+{\mu }_{D}\right)\left(1-{R}_{0}{e}^{-\lambda {\tau }_{1}}\right)$. It can be seen that $g\left(0\right) = \left({m}_{D}+{\mu }_{D}\right)\left(1-{R}_{0}\right) > 0$, if ${R}_{0} < 1$, and $\underset{\lambda \to -\infty }{\mathrm{lim}}g\left(\lambda \right) = -\infty .$ Since $g\left(\lambda \right)$ is an increasing function, given that $g\text{'}\left(\lambda \right) > 0$, then the roots of (27) must be real and negative when ${R}_{0} < 1$. Therefore, for ${\tau }_{1} > 0$ and ${\tau }_{2} = 0$, the disease-free equilibrium is locally asymptotically stable if ${R}_{0} < 1$.

Suppose ${\tau }_{1} = 0$ and ${\tau }_{2} > 0$. From (27), it yields ${\lambda }_{4} = \left({m}_{D}+{\mu }_{D}\right)\left({R}_{0}-1\right)$, which is a real negative eigenvalue if ${R}_{0} < 1$. For (28) and (29), the equations can be analyzed by using Lambert $W$ function, as illustrated by ^[21,22]. Solving (28) yields

Since $-{k}_{H}{\tau }_{2} < 0$, then the type of roots of Eq (28), based on Eq (32), are classified as follows:

1) two negative real roots; ${\lambda }_{\mathrm{5, 1}} < {\lambda }_{\mathrm{5, 2}} < -{m}_{H}$, if $0 < {\tau }_{2} < \frac{1}{{k}_{H}e}$,

2) one negative real root; ${\lambda }_{5} = -\frac{1}{{\tau }_{2}}-{m}_{H}$, if ${\tau }_{2} = \frac{1}{{k}_{H}e}$, and

3) no real roots, if ${\tau }_{2} > \frac{1}{{k}_{H}e}$.

Therefore, ${\lambda }_{\mathrm{5, 1}}$ and ${\lambda }_{\mathrm{5, 2}}$ are real negative eigenvalues if ${\tau }_{2}\le \frac{1}{{k}_{H}e}$.

Now, suppose that ${\tau }_{2} > \frac{1}{{k}_{H}e}$. Then, the roots of (28) are complex. The pair of complex roots can be written as $\lambda = u\pm i\omega$, with $u, \omega \in \mathbb{ }\mathbb{R}$. Without losing generality, substituting $\lambda = u+i\omega$ to (28), and separating real and imaginary parts results in

Taking the ratio of (33) and (34), and substituting it back to (34) yields

which gives the relation of $\omega$ and ${\tau }_{2}$. Furthermore, squaring and adding (33) and (34), and substituting the result to (33) yields

with ${\omega }^{*} = \sqrt{{\left({k}_{H}{e}^{-{\tau }_{2}\left(u+{m}_{H}\right)}\right)}^{2}-{\left(u+{m}_{H}\right)}^{2}}$. As illustrated on Figure 2(a), for ${\tau }_{2} > \frac{1}{{k}_{H}e}$, the real part of the eigenvalues of Eq (28) always lies on the negative region.

Similarly, solving for (29) yields

Since $-{k}_{D}{\tau }_{2} < 0$, then the type of roots of (29), based on (37), are classified as follows:

1) two negative real roots; ${\lambda }_{\mathrm{6, 1}} < {\lambda }_{\mathrm{6, 2}} < -{m}_{D}$, if $0 < {\tau }_{2} < \frac{1}{{k}_{D}e}$,

2) one negative real root; ${\lambda }_{6} = -\frac{1}{{\tau }_{2}}-{m}_{D}$, if ${\tau }_{2} = \frac{1}{{k}_{D}e}$, and

3) no real roots, if ${\tau }_{2} > \frac{1}{{k}_{D}e}$.

Therefore, ${\lambda }_{\mathrm{6, 1}}$ and ${\lambda }_{\mathrm{6, 2}}$ are real negative eigenvalues if ${\tau }_{2}\le \frac{1}{{k}_{D}e}$.

Now, suppose that ${\tau }_{2} > \frac{1}{{k}_{D}e}$. Then, the roots of (29) are complex. The pair of complex roots can be written as $\lambda = u\pm i\omega$, with $u, \omega \in \mathbb{ }\mathbb{R}$. Without losing generality, substituting $\lambda = u+i\omega$ to (29), and separating real and imaginary parts results in

Taking the ratio of (38) and (39), and substituting it back to (39) yields

which gives the relation of $\omega$ and ${\tau }_{2}$. Moreover, squaring and adding (38) and (39), and substituting the result to (38) yields

with ${\omega }^{*} = \sqrt{{\left({k}_{D}{e}^{-{\tau }_{2}\left(u+{m}_{D}\right)}\right)}^{2}-{\left(u+{m}_{D}\right)}^{2}}$. As shown on Figure 3(a), for ${\tau }_{2} > \frac{1}{{k}_{D}e}$, the real part of the eigenvalues of (29) always lies on the negative region. Therefore, for ${\tau }_{1} = 0$ and ${\tau }_{2} > 0$, the disease-free equilibrium is locally asymptotically stable if ${R}_{0} < 1$.

Now, suppose ${\tau }_{1}, {\tau }_{2} > 0$. Since (27) only depends on ${\tau }_{1}$, then the state of ${\tau }_{2}$ can be rejected. Therefore, the analysis will give the same result as the second case where ${\tau }_{1} > 0$ and ${\tau }_{2} = 0$, where the roots of (27) are real negative eigenvalues with no stability switch.

Furthermore, since (28) and (29) only depend on ${\tau }_{2}$, then the state of ${\tau }_{1}$ can be rejected. Thus, the analysis will give the same result as the third case where ${\tau }_{1} = 0$ and ${\tau }_{2} > 0$, where the roots of (28) are real and negative for ${\tau }_{2}\le \frac{1}{{k}_{H}e}$ and complex with negative real part for ${\tau }_{2} > \frac{1}{{k}_{H}e}$, and the roots of (29) are real and negative for ${\tau }_{2}\le \frac{1}{{k}_{D}e}$ and complex with negative real part for ${\tau }_{2} > \frac{1}{{k}_{D}e}$. Therefore, for ${\tau }_{1}, {\tau }_{2} > 0$, the disease-free equilibrium is locally asymptotically stable if ${R}_{0} < 1$.

Therefore, based on the mathematical analysis for cases of ${\tau }_{1}$ and ${\tau }_{2}$, it can be concluded that for all ${\tau }_{1}\ge 0$ and ${\tau }_{2}\ge 0$, the disease-free equilibrium is locally asymptotically stable if ${R}_{0} < 1$, meaning that the rate of production of new infections is lower than the rate of loss of infected populations.

1.3.2.   Local stability of endemic equilibrium

By substituting ${E}_{\varepsilon }$ to (25), the characteristic equation becomes

From (42), some of the eigenvalues obtained are ${\lambda }_{1} = -\left({m}_{H}+{\mu }_{H}\right)$, ${\lambda }_{2} = -{m}_{H}$, and ${\lambda }_{3} = -{m}_{D}$, which are all real negatives. The other eigenvalues are given by the following equations:

which will be analyzed by considering several cases of ${\tau }_{1}$ and ${\tau }_{2}$.

Suppose ${\tau }_{1} = 0$ and ${\tau }_{2} = 0$. Then, the other eigenvalues are given as ${\lambda }_{4} = -\left({m}_{H}+\frac{{\beta }_{HD}}{{\beta }_{DD}}\left({m}_{D}+{k}_{D}\right)\left({R}_{0}-1\right)+{k}_{H}\right)$, ${\lambda }_{5} = \frac{-\left({m}_{D}+{k}_{D}\right){R}_{0}+\sqrt{\mathcal{D}}}{2}$, and ${\lambda }_{6} = \frac{-\left({m}_{D}+{k}_{D}\right){R}_{0}-\sqrt{\mathcal{D}}}{2}$, with $\mathcal{D} = {\left(\left({m}_{D}+{k}_{D}\right){R}_{0}\right)}^{2}-4\left({m}_{D}+{k}_{D}\right)\left({m}_{D}+{\mu }_{D}\right)\left({R}_{0}-1\right).$ Eigenvalue ${\lambda }_{4}$ will be a real negative eigenvalue if ${R}_{0} > 1$. Moreover, if $\mathcal{D} < 0$, then ${\lambda }_{5}$ and ${\lambda }_{6}$ will be a pair of complex eigenvalues with negative real part, and if $\mathcal{D}\ge 0$, then ${\lambda }_{5}$ and ${\lambda }_{6}$ will be real negative eigenvalues if ${R}_{0} > 1$. Therefore, for ${\tau }_{1} = {\tau }_{2} = 0$, the endemic equilibrium is locally asymptotically stable if ${R}_{0} > 1$.

Suppose ${\tau }_{1} > 0$ and ${\tau }_{2} = 0$. From (43), the eigenvalue is given by ${\lambda }_{4} = -\left({m}_{H}+\frac{{\beta }_{HD}}{{\beta }_{DD}}\left({m}_{D}+{k}_{D}\right)\left({R}_{0}-1\right)+{k}_{H}\right)$, which is a real negative eigenvalue if ${R}_{0} > 1$. Furthermore, the other eigenvalues are given by

For (45), assume that the system experiences stability switch by the existence of a pair of imaginary eigenvalues $\lambda = \pm i\omega$, with $\omega \in {\mathbb{R}}^{+}$. Without losing generality, substituting $\lambda = i\omega$ to (45), and separating real and imaginary parts results in

Squaring and adding (46) and (47) results in

Substituting $u = {\omega }^{2}$ to (48), gives the equation

which is a quadratic equation. Suppose that ${u}_{1}$ and ${u}_{2}$ are the roots of (49), then using Vieta's formula, it is obtained that ${u}_{1}+{u}_{2} = -{\left({m}_{D}+{k}_{D}\right)}^{2}{{R}_{0}}^{2}$ and ${u}_{1}\cdot {u}_{2} = {\left({m}_{D}+{k}_{D}\right)}^{2}{\left({m}_{D}+{\mu }_{D}\right)}^{2}\left({{R}_{0}}^{2}-1\right)$. It can be seen that ${u}_{1}+{u}_{2} < 0$ and ${u}_{1}\cdot {u}_{2} > 0$ if ${R}_{0} > 1$ implying that the roots of (49) are real negative roots. Therefore, there are no real roots of (48), which implies that the eigenvalues of (45) must be real negatives, with no stability switch. Therefore, for ${\tau }_{1} > 0$ and ${\tau }_{2} = 0$, the endemic equilibrium is locally asymptotically stable if ${R}_{0} > 1$.

Next, suppose ${\tau }_{1} = 0$ and ${\tau }_{2} > 0$. (43) becomes

(50) can also be analyzed by using Lambert $W$ function. Solving (50) yields

with ${B}_{1} = \frac{{\beta }_{HD}}{{\beta }_{DD}}\left({m}_{D}+{k}_{D}{e}^{-{m}_{D}{\tau }_{2}}\right)\left({R}_{0}-1\right)$ and ${B}_{2} = -{k}_{H}{\tau }_{2}{e}^{{B}_{1}{\tau }_{2}}$. Since ${B}_{2} < 0$, then the type of roots of (50), based on (51), are classified as follows:

1) two negative real roots; ${\lambda }_{\mathrm{4, 1}} < {\lambda }_{\mathrm{4, 2}} < -\left({m}_{H}+{B}_{1}\right),$ if $-{e}^{-1} < {B}_{2} < 0$,

2) one negative real root; ${\lambda }_{4} = -\frac{1}{{\tau }_{2}}-\left({m}_{H}+{B}_{1}\right)$, if ${B}_{2} = -{e}^{-1}$, and

3) no real roots, if ${B}_{2} < -{e}^{-1}$.

Therefore, ${\lambda }_{\mathrm{4, 1}}$ and ${\lambda }_{\mathrm{4, 2}}$ are real and negative if ${B}_{2}\ge -{e}^{-1}$.

Now, suppose that ${B}_{2} < -{e}^{-1}$, then the roots of (50) are complex. The pair of complex roots can be written as $\lambda = u\pm i\omega$, with $u, \omega \in \mathbb{ }\mathbb{R}$. Without losing generality, substituting $\lambda = u+i\omega$ to (50), and separating real and imaginary parts results in

Taking the ratio of (52) and (53), and substituting it back to (53) yields

which gives the relation of $\omega$ and ${\tau }_{2}$. Moreover, squaring and adding Eqs (52) and (53), and substituting the result to Eq (52) yields

with ${\omega }^{**} = \sqrt{{\left({k}_{H}{e}^{-{\tau }_{2}\left(u+{m}_{H}\right)}\right)}^{2}-{\left(u+{m}_{H}+{B}_{1}\right)}^{2}}$. Based on Figure 4(a), for ${B}_{2} < -{e}^{-1}$, or ${\tau }_{2} > {{\tau }_{2}}^{*}\approx 1.214$, the real part of the eigenvalues of (29) always lies on the negative region. Therefore, equation (44) becomes

For (56), assume that the system experiences stability switch by the existence of a pair of imaginary eigenvalues $\lambda = \pm i\omega$, with $\omega \in {\mathbb{R}}^{+}$. Without losing generality, substituting $\lambda = i\omega$ to (56), and separating real and imaginary parts results in

Squaring and adding (57) and (58) results in

with ${p}_{1} = {m}_{D}+{\mu }_{D}$, ${p}_{2} = {k}_{D}{e}^{-{m}_{D}{\tau }_{2}}$, and ${p}_{3} = {m}_{D}+{k}_{D}{e}^{-{m}_{D}{\tau }_{2}}$. Substituting $u = {\omega }^{2}$ to (59), gives the equation

which is a quadratic equation. Suppose that ${u}_{1}$ and ${u}_{2}$ are the roots of Eq (60), then using Vieta's formula, it is obtained that ${u}_{1}+{u}_{2} = -{p}_{3}\left({p}_{3}{{R}_{0}}^{2}-2{p}_{1}\left({R}_{0}-1\right)-2{p}_{2}{R}_{0}\right)$ and ${u}_{1}\cdot {u}_{2} = {{p}_{1}}^{2}{{p}_{3}}^{2}{\left({R}_{0}-1\right)}^{2}$. For positive roots to exist for Eq (60), it is required that, by using Vieta's formula, ${p}_{3}{{R}_{0}}^{2}-2{p}_{1}\left({R}_{0}-1\right)-2{p}_{2}{R}_{0} < 0$, which is implicit on ${\tau }_{2}$. Therefore, the critical value of ${\tau }_{2}$ is given by ${p}_{3}{{R}_{0}}^{2}-2{p}_{1}\left({R}_{0}-1\right)-2{p}_{2}{R}_{0} = 0$.

To verify the given critical value of ${\tau }_{2}$ causes Hopf bifurcation, we can check the transversality condition, as performed by ^[17]. The transversality condition is satisfied when $\mathrm{sign}\left[{\left.\frac{dRe\left(\lambda \right)}{d{\tau }_{2}}\right|}_{{\tau }_{2} = {{\tau }_{2}}^{*}}\right] = \mathrm{sign}\left[{\left.\frac{dF\left(u\right)}{du}\right|}_{u = {\omega }^{2}}\right]$, where $F\left(u\right)$ is described as the left-hand side of (60). Differentiating (56) with respect to ${\tau }_{2}$, and rearranging results in

Since for ${\tau }_{2} = {{\tau }_{2}}^{*}$, $\lambda = i\omega$, then by letting ${q}_{1} = {p}_{3}{R}_{0}-{p}_{2},$ ${q}_{2} = {p}_{1}{p}_{3}\left({R}_{0}-1\right)$, ${\upsilon }_{1} = {m}_{D}{p}_{2}-{m}_{D}{q}_{1}-{q}_{2},$ ${\upsilon }_{2} = {m}_{D}+{q}_{1},$ ${\upsilon }_{3} = {m}_{D}\left({p}_{1}{p}_{2}-{q}_{2}\right)$, then from (61), it yields

or

Taking the real part of (62) results in

which is known that $\mathrm{sign}\left\{Re\left[{\left(\frac{d\lambda }{d{\tau }_{2}}\right)}_{{\tau }_{2} = {{\tau }_{2}}^{*}}^{-1}\right]\right\} = \mathrm{sign}\left\{{\left.\frac{dRe\left(\lambda \right)}{d{\tau }_{2}}\right|}_{{\tau }_{2} = {{\tau }_{2}}^{*}}\right\}$. Furthermore, from Eq (60), it is obtained that

From (63) and (64), it can be concluded that the transversality condition is not satisfied. Therefore, the critical value of ${\tau }_{2}$, obtained from (60), is not the bifurcation point of the equilibrium. Therefore, the given critical value does not cause stability switch for the system. Moreover, for ${\tau }_{1} = 0$ and ${\tau }_{2} > 0$, the endemic equilibrium is locally asymptotically stable if ${R}_{0} > 1$.

Suppose ${\tau }_{1}, {\tau }_{2} > 0$. Eq (43) can be reduced to only depending on ${\tau }_{2}$, which results in Eq (50). Therefore, the analysis can be done as in previous case, which will give similar result where the eigenvalues will always be negative. Moreover, Eq (44) can be analyzed by assuming the existence of a pair of imaginary eigenvalues $\lambda = \pm i\omega$, with $\omega \in {\mathbb{R}}^{+}$.

By letting ${r}_{1} = \left({m}_{D}+{\mu }_{D}\right)\sin\omega {\tau }_{1}$, ${r}_{2} = {\omega }_{4}+\left({m}_{D}+{\mu }_{D}\right)\sin\omega {\tau }_{1}$, ${r}_{3} = \left({m}_{D}+{\mu }_{D}\right)\left(1-\cos\omega {\tau }_{1}\right)$, ${r}_{4} = \left({m}_{D}+{\mu }_{D}\right)\left({m}_{D}+{k}_{D}{e}^{-{m}_{D}{\tau }_{2}}\right)\left({R}_{0}-1\right)$, ${r}_{5} = {m}_{D}+{k}_{D}{e}^{-{m}_{D}{\tau }_{2}}$, ${s}_{1} = {\omega }^{2}+\omega {r}_{1}-{m}_{D}{r}_{3}-{r}_{4}$, and ${s}_{2} = \omega {r}_{5}{R}_{0}+{m}_{D}{r}_{1}+\omega {r}_{3}+\omega {k}_{D}{e}^{-{m}_{D}{\tau }_{2}}$, the critical point of ${\tau }_{2}$ is given as an implicit function based on ${\tau }_{1}$ and ${\tau }_{2}$, i.e.,

The only thing left is to verify the existence of Hopf bifurcation by the critical delay ${\tau }_{2}$. Based the transversality analysis done within the previous case, it can be concluded that since ${\tau }_{1}$ does not cause any Hopf bifurcation, then the given critical delay of ${\tau }_{2}$, that depends on ${\tau }_{1}$, does not cause any Hopf bifurcation. Therefore, for ${\tau }_{1}, {\tau }_{2} > 0$, the endemic equilibrium is locally asymptotically stable if ${R}_{0} > 1$.

Based on mathematical analysis for cases of ${\tau }_{1}$ and ${\tau }_{2}$, it can be concluded that for all ${\tau }_{1}\ge 0$ and ${\tau }_{2}\ge 0$, the endemic equilibrium is locally asymptotically stable if ${R}_{0} > 1$, meaning that the rate of production of new infections is higher than the rate of loss of infected populations. From the local stability analysis on both equilibria, the existence of two delays on the system does not cause the system to experience Hopf bifurcation as expected. Instead, the stability of both equilibria relies heavily on ${R}_{0}$, where the stable equilibrium shifted as ${R}_{0}$ passes 1. This is known as Transcritical bifurcation ^[23].

4.   Numerical simulation

The numerical simulation for this research is done using Python version 3.9.13. The numerical simulation is aimed to back up the stability analysis result that has been obtained previously for both equilibrium states, as well as to observe the effects on the variations of some parameters in reducing or eradicating the disease. Numerical simulation is obtained by first solving the system numerically using the Python software. The system is solved using fourth order Runge-Kutta that has been modified to solve the delay system on (1)–(6). The parameter values used for the simulation is given in Table 2.

Calculation of a basic reproduction number using values given by Table 2 results in ${R}_{0}\approx 0.335 < 1$, which satisfies the stability condition of the disease-free equilibrium ${E}_{0}$. Moreover, for the endemic equilibrium ${E}_{\varepsilon }$, the parameter values used for the numerical simulation are the same as those given in Table 2, except for ${\beta }_{HD} = 0.001$, ${\beta }_{DD} = 0.01$, ${k}_{H} = 0.25$, and ${k}_{D} = 0.25$. This is required to satisfy the existence and stability conditions of ${E}_{\varepsilon }$, i.e., ${R}_{0} > 1$. Calculation on basic reproduction number using these modified parameter values results in ${R}_{0}\approx 6.055 > 1$.

Figure 5 shows that all the solutions are asymptotically stable and tends to ${E}_{0}$ as $t\to \infty$. From Figure 5(c), the infection on human populations gradually decreases, and after about 8 years, the disease dies out, without any spike on case numbers due to low contact among susceptible humans and rabid dogs. From Figure 5(a), the susceptible human population gradually decreases, indicating that a lot of susceptible subpopulation transition to vaccinated subpopulation due to high vaccination coverage. It can be seen in Figure 5(e) at the vaccinated human population keeps increasing, with slower growth each time, indicating that it will eventually reach an equilibrium state. From Figure 5(d), the infection on dog populations experiences an increase for the number of cases, albeit only less than 30 cases within the first years, and then the infection gradually decreases until about 10 years when the disease finally dies out. Furthermore, from Figure 5(f), it shows that the vaccinated dog population keeps increasing due to high vaccination coverage, corresponding to susceptible dog population that experiences decrement.

Figure 6 shows that all solutions are stable and tends to ${E}_{\varepsilon }$ as $t\to \infty$, albeit there is a bit of fluctuation within the first years of the epidemic. As seen in Figure 6(c), within the first years, the disease experiences a huge spike in the number of cases from 500 up to around 6,800 cases due to high contact rate among susceptible humans and rabid dogs. This correlates to the rapid decrement of susceptible human population, as can be seen on Figure 6(a). Moreover, in Figure 6(e), the vaccinated population grows, albeit not as rapid as on ${E}_{0}$, which only reached approximately 24,000 humans after 50 years. This is the impact of low vaccination coverage and high contact rate, which causes a high number of casualties due to the disease. From Figure 6(d), within the first years of the epidemic, there is a spike in the number of cases from 50 up to 1,300 cases within a year. This impacted the rapid decrement of the susceptible dog populations, as well as the rapid increase of the case numbers in human population. Albeit after this spike, the case number flattens, until it reaches an equilibrium state after about 5 years. Although the case flattens quicker than on ${E}_{0}$, the high number of cases causes high number of causalities because of the disease.

The effects of varying parameters and the impact to the number of infected humans and dogs are also observed, i.e., the effects of variation on dog vaccination rate, ${k}_{D}$, and the average birth of puppies, ${A}_{D}$. Unlike the methods presented in ^[24,25,26], where they use game theory strategies in vaccination, the vaccination factor here does not depend on the voluntary decision of each person, but is rather forced to achieve a certain percentage of vaccination, as the cost of vaccination is not considered here. From Figure 7, by increasing vaccination coverage, it can help reduce the number of cases during the spike of cases. However, it requires a high coverage of vaccination to completely eradicate the disease, up to ${k}_{D} = 2$, which is over 100% vaccination coverage. From Figure 8, even though reducing the average birth of puppies to 30 is sufficient to eliminate the disease, it does not help to reduce the number of cases spike during the first years of the epidemic.

Therefore, it may help to combine several control strategies to reduce the number of cases spiking, or to allow the disease to die out without having to experience a spike in the number of cases. For example, by combining the strategy to increase vaccination coverage: Limiting contact with rabid dogs, as well as reducing the number of puppies born, this can be done by raising awareness to pet owners to vaccinate their pets, especially dog owners. Educating the mass to detect signs of rabies in animals may help people to avoid interactions with possibly infected animals. Controlling wild dog populations is sufficient to help reduce the birth of wild puppies on the streets, which are more susceptible to catch the disease and spread them.

The success of vaccination in controlling the disease depends on the success on other control measures. This is different from using voluntary decision, as demonstrated in ^[24], where vaccination cost determines the decision of taking vaccination. This would require more pressure on other controls that do not require voluntary decision making.

As an addition, the effects of variation of ${\tau }_{1}$ and ${\tau }_{2}$ will also be observed here to see how both delay affects the solution of the system. From Figure 9, extending the delay of incubation allows the disease to experience a smaller spike of cases. However, extending the incubation period does not help with eradicating the disease. Therefore, if there exists a way to prolong incubation, this requires a combined control strategy to eliminate the disease. It can also be seen in Figure 10 that variation on vaccination latency also does not help with eradicating the disease, nor helps by reducing the number of cases spike. Therefore, any efforts of controlling vaccination latency should be accompanied by other control strategies to help eliminate the disease.

5.   Conclusions

In this article, a mathematical model of rabies disease spread among human and dog populations has been developed by incorporating two discrete delays on incubation and vaccination. Several key findings that were discussed here are about the local stability of the system, the strategies that can be implemented to reduce or eradicate diseases, and several potential developments that can be considered for future studies. The stability of the model is analyzed around two types of equilibrium: The disease-free equilibrium, and the endemic equilibrium. The disease-free equilibrium is locally asymptotically stable for all positive delays if ${R}_{0} < 1$, and the endemic equilibrium is locally asymptotically stable for all positive delays if ${R}_{0} > 1$. It is also discovered that both delays affect the spread of rabies based on the basic reproduction number. However, based on the local stability analysis, it was found that both delays do not affect the stability switch on the local stability of both equilibria. In order to reduce or eradicate the disease, applying one control strategy is not enough. Multiple control strategies are highly suggested in order to eradicate the disease. For example, by combining the strategy to increase dog vaccination coverage, limiting the contact with rabid dogs, and decreasing the number of births of puppies, they are sufficient to help in eradicating the disease. It is also found that there are factors that are overlooked within this research, which may become a potential development for future studies, such as vaccines losing their efficacy over time and incorporating diffusion factors that give a new insight of the model as a PDE system.

Use of AI tools declaration

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

Acknowledgments

The authors thank the anonymous reviewers for their valuable comments and suggestions that improved the quality of this paper.

Conflicts of interest

The authors declare no conflicts of interest regarding the publication of this article.