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

: In this article, a delay differential equations model is constructed to observe the spread of rabies among human and dog populations by considering two delay effects on incubation period and vaccine efficacy. Other parameters that affect the spread of rabies are also analyzed. Using the basic reproduction number, it is shown that dog populations and the two delays gives a significant effect on the spread of rabies among human and dog populations. The existence of two delays causes the system to experience Transcritical bifurcation instead of Hopf bifurcation. The numerical simulation shows that depending only on one control method is not enough to reduce or eradicate rabies within the dog populations; instead, it requires several combined strategies, such as increasing dog vaccinations, reducing contact with infected dogs, and controlling puppies’ birth. The spread within the human population will be reduced if the spread within the dog population is reduced.


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.

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 (  ), infected human population (  ), and vaccinated human population (  ).The dog population also consists of three subpopulations: A susceptible dog population (  ), infected dog population (  ), and vaccinated dog population (  ).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 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   (),   (), and   () are given as   (0) ≥ 0,   (0) ≥ 0, and   (0) ≥ 0, respectively, and historical functions of   (),   (), and   () are given as  1 () ≥ 0,  2 () ≥ 0, and  3 () ≥ 0, respectively.The functions of  1 () and  2 () are continuous functions within the interval [− 2 , 0], and  3 () is a continuous function within the interval [− 1 , 0].

Boundedness of the solution
Now, suppose that   =   +   +   and   =   +   +   , where   is the total human population and   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   tends to as  → ∞.Similarly, from (10), it can be seen that as  → ∞ .Therefore, the feasible region of the solutions of human and dog populations are given as } .Therefore, the model ( 1)-( 6) is mathematically well-posed and epidemiologically meaningful.

Disease-free equilibrium
From (15), it can be seen that   * = 0 or   * ≠ 0, which yields two different equilibria.For the case of   * = 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.

Basic reproduction number
Basic reproduction number ( 0 ) 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  0 here is by using next generation method, as given by [20].First, taking the infectious compartments from ( 2) and ( 5) by letting  1 =  ̇ and  2 =  ̇ gives The next generation matrix is given by   = − −1 , where the matrix  describes the events of producing new infections and the matrix  describes the events of infectious population transitions.Matrices  and  are constructed such that (17) can be written as  ̇= ( + ), where  = (  ,   ) .By taking the derivatives of  1 and  2 with respect to   and   and evaluating the derivatives around  0 , the matrices  and  are given as From ( 18) and ( 19), the next generation matrix is given as Therefore, by taking the spectral radius of   , the basic reproduction number is given as From ( 21), it can be seen the value of  0 depends only on dog parameters, as well as both delays  1 and  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.

Endemic equilibrium
For the case of   * ≠ 0, after substituting and solving ( 11)-( 16), it gives the endemic equilibrium given as   = (   ,  where the endemic equilibrium exists for  0 > 1.
Then, by using Taylor series expansion around  * on (1) -( 6), it results in where  0 is the Jacobian matrix with respect to ,  1 is the Jacobian matrix with respect to   1 , and  2 is the Jacobian matrix with respect to   2 .The matrices  0 ,  1 ,  2 are described as follows Since the system is now in a linearized form, therefore the solution must in the form () =   , where  is a non-zero vector.From (22), the characteristic equation is given by Applying  0 ,  1 , and  2 to (23) yields with From ( 26), some of the eigenvalues obtained are  1 = −(  +   ) ,  2 = −  , and  3 = −  , which are all real negatives.The other eigenvalues are given by the following equations: which will be analyzed by considering several cases of  1 and  2 .
Furthermore, since (28) and (29) only depend on  2 , then the state of  1 can be rejected.Thus, the analysis will give the same result as the third case where  1 = 0 and  2 > 0, where the roots of (28) are real and negative for  Therefore, for  1 ,  2 > 0, the disease-free equilibrium is locally asymptotically stable if  0 < 1.
Therefore, based on the mathematical analysis for cases of  1 and  2 , it can be concluded that for all  1 ≥ 0 and  2 ≥ 0, the disease-free equilibrium is locally asymptotically stable if  0 < 1, meaning that the rate of production of new infections is lower than the rate of loss of infected populations.

Local stability of endemic equilibrium
By substituting   to (25) which will be analyzed by considering several cases of  1 and  2 .

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. Latency time for vaccination 1/10 year assumed Calculation of a basic reproduction number using values given by Table 2 results in  0 ≈ 0.335 < 1, which satisfies the stability condition of the disease-free equilibrium  0 .Moreover, for the endemic equilibrium   , the parameter values used for the numerical simulation are the same as those given in Table 2, except for   = 0.001,   = 0.01,   = 0.25, and   = 0.25.This is required to satisfy the existence and stability conditions of   , i.e.,  0 > 1 .Calculation on basic reproduction number using these modified parameter values results in  0 ≈ 6.055 > 1.
Figure 5 shows that all the solutions are asymptotically stable and tends to  0 as  → ∞ .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   as  → ∞, 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  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  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,   , and the average birth of puppies,   .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   = 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 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  1 and  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.Volume 9, Issue 2, 3399-3425.

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 diseasefree equilibrium, and the endemic equilibrium.The disease-free equilibrium is locally asymptotically stable for all positive delays if  0 < 1, and the endemic equilibrium is locally asymptotically stable for all positive delays if  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.

Figure 2 .
Figure 2. The numerical plots of (a) real part of the complex eigenvalues of (28), based on (36), and (b) imaginary part of the complex eigenvalues of (28), based on (35).

Figure 3 .
Figure 3.The numerical plots of (a) real part of the complex eigenvalues of (29), based on (41), and (b) imaginary part of the complex eigenvalues of (29), based on (40).
and complex with negative real part for  2 > 1    , and the roots of (29) are real and negative for  2 ≤ 1    and complex with negative real part for  2 > 1    .

Figure 1. Schematic diagram of the rabies spread model. The susceptible dog population is increased by average dog birth per year (𝐴 𝐷 ). It is decreased by natural death at rate 𝑚 𝐷 , infection by contact with infected dogs at rate 𝛽 𝐷𝐷 𝑆 𝐷 𝐼 𝐷 𝜏 1 𝑒 −𝑚 𝐷 𝜏 1 , and vaccination at rate 𝑘 𝐷 𝑆 𝐷 𝜏 2 𝑒 −𝑚 𝐷 𝜏 2 . Here, the 𝑒 −𝑚 𝐷 𝜏 2 is the probability of susceptible dogs surviving natural
, the susceptible human population is increased by average human birth per year (  ).It is decreased by natural death at rate   , infection by contact with infected dogs at rate        1  −   1 , and vaccination at rate      2  −   2 .Here,  1 is the incubation delay and  2 is the vaccination delay.Moreover,  −   1 is the probability of rabid dogs surviving natural death within the time period [0,  1 ] ,  −   2 is the probability of susceptible humans surviving natural death within the time period [0,  2 ],   is the transmission rate between susceptible humans and infected dogs, and   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        1  −   1 , and decreased by natural and rabies-related deaths at rate   and   , respectively.The vaccinated human population is increased by vaccination of susceptible humans at rate      2  −   2 and decreased by natural death at rate   .death within the time period [0,  2 ],   is the transmission rate between susceptible dogs and infected dogs, and   is the vaccination rate of dogs.The infected dog population is increased through infection by contact between susceptible and infected dogs at rate        1  −   1 , and decreased by natural and rabies-related deaths at rate   and   , respectively.The vaccinated dog population is increased by vaccination of susceptible dogs at rate      2  −   2 and decreased by natural death at rate   .Table 1 summarizes the description of all variables and model parameters

Table 2 .
Values of the model's parameters used for numerical simulation.