J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 Journal of the Nigerian Society of Physical Sciences Dynamics of Toxoplasmosis Disease in Cats population with vaccination Idris Babaji Muhammada, Salisu Usainib,∗ aDepartment of Mathematics ,Faculty of Science, Bauchi State University, Gadau P.M.B65,Gadau, Nigeria bDepartment of Mathematics , Kano University of Science and Technology, Wudil, P.M.B. 3244, Kano, Nigeria Abstract We extend the deterministic model for the dynamics of toxoplasmosis proposed by Arenas et al. in 2010, by separating vaccinated and recovered classes. The model exhibits two equilibrium points, the disease-free and endemic steady states. These points are both locally and globally stable asymptotically when the threshold parameter Rv is less than and greater than unity, respectively. The sensitivity analysis of the model parameters reveals that the vaccination parameter π is more sensitive to changes than any other parameter. Indeed, as expected the numerical simulations reveal that the higher the vaccination rate of susceptible cats the smaller the value of the threshold Rv (i.e., increase in π results in the decrease in Rv), leading to the eradication of toxoplasmosis in cats population. DOI:10.46481/jnsps.2021.141 Keywords: Toxoplasmosis, Vaccination, Cats, sensitivity index Article History : Received: 31 October 2020 Received in revised form: 24 November 2020 Accepted for publication: 05 January 2021 Published: 27 February 2021 c©2021 Journal of the Nigerian Society of Physical Sciences. All rights reserved. Communicated by: T. Latunde 1. Introduction The causative agent of toxoplasmosis is toxoplasma gondii (T. gondii), which is a protozoan of the order Eucoccidioroda and the family Sarcocystidae. T. gondii is characterized as an intracellular parasite that lives in the host cell by regulating vi- tal processes to acquire nutrients, guaranteeing its survival and thus evading the host immune system [12]. It is an infectious pathogen which enters its host through ingestion of either the oocysts, the trachyzoite, or tissue cysts (bradyzoites) from con- taminated water, soil, or infected meat. The three major types of disease reservoirs of T. godii are: Domestic and wild cats, soil and water contaminated with feces (Non-living reservoirs) and ∗Corresponding author tel. no: +2348064237334 Email address: salisu.usaini@kustwudil.edu.ng (Salisu Usaini ) bradyzoites in tissue cysts of intermediate hosts (animal reser- voirs) [25]. The intermediate hosts of T. gondii include sheep, goats, rodents, cattle, swine, chicken and birds. The modes of T. gondii transmissions include ingestion of oocysts that excrete in cats feces for which tissue cysts develop through exposure to cat litter or soil and water. This is the most common and well known modes of transmission to other animals. T. gondii can also be transmitted via consumption of T. gondii tissue cysts in raw or undercooked meats, unpasteurized milk and consump- tion of oocysts in foods infected by contaminated fomites. This is called the food born transmission. Another route of transmis- sion is congenital toxoplasmos (transplacental transmission) in which a mother infect her offspring in uterus through the blood during pregnancy. A pregnant women with acute infections can transmit to the fetus and cause severe illness such as men- tal retardation, blindness, and epilepsy [12, 25]. It was re- 17 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 18 ported in a recent study in dogs that sexual transmission of T. gondii in canine species is possible [25]. Treatment of tox- oplasmosis include sulfonamides against murine toxoplasmo- sis, combined therapy with sulfonamides and pyrimethamine as the standard therapy for toxoplasmosis in humans, spiramycin during pregnancy to reduce transmission of the parasite from mother to fetus, clindamycin especially for patients allergic to sulphonamides [12]. A number of mathematical epidemic models for the trans- mission dynamics of toxoplasmosis parasites with or without vaccination were proposed in the literature. In 2008, Aranda et al. [6], proposed a simple mathematical model to study the dy- namics of Toxoplasmosis disease in the population of Colom- bia. Numerical simulations of the model using available data reveals the effect of some (hygienic actions, educations pro- grams, more testing and treatments) strategies for the control of toxoplasmosis parasite. This model seems to be the first mathe- matical model for the transmission dynamics of Toxoplasmo- sis disease in the human population. In the following year, an epidemic model of toxoplasmosis disease in human and cat populations was presented in [9]. The analysis of the model indicated that toxoplasmosis disease persistence or extinction is determined by the threshold parameter, R0. Moreover, they showed by numerical simulations the importance of cats verti- cal transmission to the dynamics of the infection. In fact, the higher the vertical transmission the higher the proportions of infected cats and infected humans at the endemic state. Thus, the vertical transmission could be an important mechanism that favours the maintenance of the virus areas with low human den- sities. Another epidemiological model of toxoplasmosis in a cat population with a continuous vaccination schedule was pro- posed in [2]. Indirectly, the model considers the infection of prey through the oocyst shedding by cats. They proved that the global dynamics and disease outcome are solely determined by the basic reproduction number, R0. Furthermore, numerical simulations of the model suggest that the continuous vaccina- tion program is more effective than the removing of oocysts in the environment since the increase of the latter in the environ- ment to feasible values is not enough to reduce the threshold parameter, R0 to a value less than unity. The spatial spread of toxoplasmosis parasite was also con- sidered later on. In 2012, a spatial model for the spread of tox- oplasma gondii through a heterogeneous predatorprey system was proposed in [16]. They considered some relevant toy mod- els due to the complexity of the proposed model. Then they proved the existence and local stability of a persistent steady state for the underlying predatorprey model systems. A spa- tial mathematical model for the reproduction dynamics of the Toxoplasma gondii parasite in the definitive host Felis catus (cat) was presented in [10]. Both asexual and sexual repro- duction processes of the T. gondii parasite were incorporated in the model. Some numerical results showed that variations in the transition and loss rates do not produce significant changes in the reproduction, propagation and creation of new popula- tions. On the other hand, with either low or high consumption of oocysts from the environment by the cat, the involved pop- ulations are always reproduced. Then they spread by all over epithelial cells and subsequently are expelled to the environ- ment through the cat feces. Ferreira et al. [11], studied the dynamical behaviours of both deterministic and spatial models of toxoplasmosis disease in cat and human populations. They showed that the deterministic model exhibits a trans-critical bi- furcation and the system has no periodic orbits inside the pos- itive octant. The global asymptotic stability of the endemic equilibrium point of the model was proven in the first octant. Moreover, they carried out the local stability analysis for the spatially homogeneous equilibrium points of the reaction diffu- sion model. In the restricted region (the first octant), the global stability of both the disease-free and endemic states of such a model were established. In the last decade, a hybrid model for the spread of tox- oplasmosis between two cities of Colombia in cat and human populations was proposed by Peña et al. in [4]. The model was developed using System Dynamics (SD) and Geographic In- formation Systems (GIS). The analysis of the model indicates that the disease persistence continues in the first city, even after transported to another community. This demonstrated that these diseases should be treated through cooperative mechanisms be- tween communities. A deterministic mathematical model of interaction between Toxoplasma Gondii transmission dynamics and host immune responses was developed and presented in [1]. The local asymptotic stability analysis of the model equilibria and numerical simulations were carried out to understand the transmission dynamics and the impact of different functional responses. In fact, the Holling type II functional response en- hances the effect or cells of hosts immune response. In a recent development, an epidemic mathematical model of toxoplasmosis disease was presented in [15] to study the dy- namic behaviour of the parasite for controlling the epidemic. The existence and stability analysis of each equilibrium point of the model was carried out. They constructed a function of two important parameters of the model namely the controlled rate and the rate of infected births which influences the dis- ease spread. The bifurcation analysis of each region was car- ried out from each function of such parameters. Three regional conditions were obtained through this analysis, showing the dynamics of the toxoplasmosis epidemic of these two impor- tant parameters with each interpretation of the bifurcation re- gion. The widespread of co-infection by Plasmodium species and Toxoplasma gondii in humans mostly in the tropical re- gion motivated Ogunmiloro to study the co-infection dynam- ics of malaria-toxoplasmosis in the mainly human and feline susceptible host population in [18]. The analysis showed that malaria-toxoplasmosis infection free equilibrium is locally and globally stable if the basic reproduction number is less than unity. While the co-infection occurs when the basic reproduc- tion number is greater than unity. Furthermore, sensitivity anal- ysis of the model parameters indicated the need for proper op- timum medical strategies to reduce and contain the spread of malaria-toxoplasmosis co-infections. In 2010, a deterministic mathematical model for the trans- mission dynamics of toxoplasmosis disease in cats population under vaccination was proposed in [2]. Implicitly, their model considers the infection of prey through the oocyst shedding by 18 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 19 cats in the environment. The analysis of the model shows that the global dynamics and disease outcome are completely de- termined by the basic reproduction number, R0. Moreover, the numerical results reveal that the effectiveness of the continuous vaccination program is more than the removing of oocysts in the environment. This is for the fact that an increase of the lat- ter in the environment to feasible values is not enough to control disease persistence by reducing the threshold parameter, R0 to a value less than one. It is instructive to note that in their model the natural immunity and vaccine protection were assumed to be the same and provide permanent immunity. Therefore, they merged the vaccinated cats and non vaccinated recovered cats in the same class. We extend their model here by separating vaccinated cats from non vaccinated recovered cats to have two classes. This is done for the fact that there is increasing clin- ical evidence which differentiate between these two immuni- ties against toxoplasmosis in cats [5, 7, 8, 21]. As in [2], we show that the disease outcome, persistence or extinction is de- termined by the vaccinated reproduction number Rv. The sensi- tivity analysis of the model parameters indicates that the vacci- nation parameter π is more sensitive to changes than any other parameter of the model. Furthermore, numerical simulations show that the higher vaccination rate of susceptible cats reduces the value of the threshold parameter Rv the more. This is also similar to the result presented in [2] that vaccination program is more effective than the removal of oocysts in the environment. We organise the remaining part of the paper as follows: In Section 2, model description is provided, Section 3 presents the stability analysis of the model equilibria while sensitivity anal- ysis of the model parameters is presented in Section 4. Finally, we discuss the results and give a concluding remark in Section 5. 2. Model Formulation To extend the model proposed by Arenas et al. [2], we sepa- rate the compartment of vaccinated cats from that of recovered (vaccinated) cats. Thus we divide the total population of cats to four compartments of susceptible (S(t)), cats that may be- come infected after an effective contact with oocysts, infective (I(t) with T. gondii, vaccinated (V (t)), and non vaccinated re- covered R(t). Using the following notations and hypotheses as presented in [2], we have the required model (1). The rates of transfer between the four classes are shown in Figure 1. 2.1. Basic Assumptions i. O(t) repreaent the number of oocysts in the environment at time t, ii. The birth rate (δ) of cats is assumed to be equal to the cats natural death rate µ, iii. The non vaccinated susceptible cats transits to an infec- tive compartment after an effective contact with oocysts at the rate β. Thus, β depends on the probability of an effective contact with oocysts, τR μI μ0O μS �I αI μV μRγV πS δI βSO δ(S+V+R) S I V R O Figure 1: Schematic diagram of model (1). iv. Vaccinated susceptible cat S (t) transits to the vaccinated compartment (V (t)) at a rate π and infected cats to recov- ered compartment R(t) at a rate α, v. Number of oocysts, O(t) at time t determine the number of infected cats I(t), vii. µ0 is the death rate of the oocysts which can be modified by the removal of oocysts from the environment, vi. It is assumed that there is vertical transmission in the in- fected cats population (i.e., transmission from mother to fetus), viii. Vaccinated cats do not shed oocysts and acquire perma- nent immunity so that they move to recovered compart- ment R(t). The recovered cats can be re-infected with T.gondii but do not shed oocysts again and so, they move to susceptible class at the rate τ, ix. The susceptible cats S (t) have the same probability of be- coming infected (i.e., homogeneous mixing is assumed) x. κ > 0 is the rate of appearance of new oocysts from the infected cats into the environment, xi. The parameter γ is the progression rate from vaccinated class to the recovered compartment. dS dt = µV + µR −βS O −πS + τR, dI dt = βS O −αI, dV dt = πS −µV −γV, dR dt = αI −µR −τR + γV, dO dt = κI −µ0O, (1) with N(t) = S + I + V + R = 1. In the first equation of model (1), the rate of change of the susceptible cats population is increased by per capita birth rate 19 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 20 Table 1: parameter description Parameter Description Value Source β Effective contact rate 0.5254 Wu and Pei π Vaccination rate of susceptible 0.2 Arena et al. α HBV induced death rate 0.5 Arenas et al. µ0 Death rate of Oocysts 7 100000 Wu and Pei. µ Natural birth rate 0.652 Wu and Pei τ Immunity rate 0.452 assumed κ Rate of appearance of oocyst 120 Wu and Pei γ Progression rate of vaccinated cats 0.3 assumed δ of susceptible, vaccinated and recovered classes, the waning immunity after recovery at the rate τ. It decreases by vaccina- tion at the rate, π, natural death µ and effective contact at the rate β between the oocysts and the healthy cats. The rate of change of Infected cats is increased by effective contact of the susceptible cats and the oocysts at the rate β, with birth rate δ of the infected cats. While it decreases by contracting the disease at the rate α. The third equation of model 1 is increased by vac- cinating susceptible cats at the rate π and decreases by natural death µ and the rate at which they move to recovered class γ. Similarly, the rate of change of recovered cats is increased at the per capita rates of recovery of infected and vaccinated cats α and γ, respectively. It is reduced by natural death at the rate µ and wane of immunity from the recovered class at the rate τ. Finally, the rate of change of oocyst with time is increased by shedding the virus at the rate, κ and reduced by the natural death of oocysts at the rate µ0. Using R(t) = 1 − S − I − V , the system of equations (1) becomes. dS dt = (µ + τ)(1 − S − I) −τV −βS O −πS, dI dt = βS O −αI, dV dt = πS −µV −γV, dO dt = κI −µ0O (2) Theorem 1. Model (1) is positively invariant and attractive in the biologically feasible region Ω0 = { (S, I, V, O) ∈ R4+ : 0 < S ≤ (µ + τ)(µ + γ) (µ + γ)(µ + π + τ) + πτ , 0 6 O 6 κ µ 0 } . Proof. Using the first and fourth equations of system (2), we have dS dt ≤ (µ + τ)(µ + γ) − [(π + µ + τ)(µ + γ) + τπ]S and dO dt ≤ κ−µ0O. Then dS dt + [(π + µ + τ)(µ + γ) + τπ]S ≤ (µ + τ)(µ + γ), so that S (t) ≤ (µ + τ)(µ + γ) (π + µ + τ)(µ + γ) + τπ and O(t) ≤ O(0)e−µ0 t + κ µ 0 (1 − e−µ0 t), with O(t) ≤ κ µ 0 if O(0) < κ µ 0 as t approaches ∞. Thus Ω is positively invariant. Moreover, if S (0) > (µ + τ)(µ + γ) (π + µ + τ)(µ + γ) + τπ and O(0) > κ µ 0 then the solution either enters Ω in finite times or S (t), O(t) ap- proaches (µ+τ)(µ+γ)(π+µ+τ)(µ+γ)+τπ , κ µ 0 asymptotically [2]. Hence the re- gion attract all solutions in R4+. Thus it is sufficient to consider the dynamics of the model (2) in Ω since it is mathematically well-posed. 3. Model Analysis 3.1. Disease Free Equilibrium (DFE) In the absence of infection (i.e., when all the infected com- partments of the model (2) are empty), we have the disease free equilibrium point of the model denoted by P0. Then P0 = (S 0, I0, V0, O0) = ( (µ + τ)(µ + γ) (π + µ + τ)(µ + γ) + τπ , 0, πS 0 µ + γ , 0 ) . For us to establish the linear stability of equilibria, we need a threshold parameter usually called the basic reproduction num- ber. This threshold parameter is defined as the number of sec- ondary infections produced by one infective when introduce into the susceptible population [13]. The next generation op- erator method can be used to compute such a threshold param- eter as presented in [20]. Moreover, when there is a vaccination such a parameter is called a vaccinated reproduction number (see for instance, [22]). The method consists of two different matrices, the matrix of the new infection terms and that of the transition terms denoted by G and M, respectively. Thus 20 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 21 G = [ βS 0 0 0 ] , M = [ 0 α µ0 −κ ] . Then GM−1 = [ κβ(µ+τ)(µ+γ) αµ0 [(µ+π+τ)(µ+γ)+τπ] κβ(µ+τ)(µ+γ) (µ+π+τ)(µ+γ)+τπ 0 0 ] Hence, vaccinated reproductive number Rv is the largest eigen- value of the matrix GM−1. Therefore Rv = ρ(GM −1) = κβ(µ + τ)(µ + γ) αµ0[(µ + π + τ)(µ + γ) + τπ] . 3.1.1. Local Stability Analysis of DFE Theorem 2. The DFE, P0 of model (2) is locally asymptoti- cally stable on Ω if Rv < 1 and is unstable whenever Rv > 1. Proof. The Jacobian matrix of system (2) evaluated at the DF E is as follows J(P0) =  −(µ + π + τ) −(µ + τ) −τ −βS 0 0 −α 0 βS 0 π 0 −(µ + γ) 0 0 κ 0 −µ0  . Then, interchanging row-3 and row-4 and using row operations, the matrix J(P0) becomes J(P0) =  −(µ + π + τ) −(µ + τ) −τ −βS 0 0 −α 0 βS 0 0 κ 0 −µ0 π 0 −(µ + γ) 0  . Divide row 4 by π, then new row 4 is R4 = (π + τ + µ)R4 + R1 and so J(P0) =  −(µ + π + τ) −(µ + τ) −τ −βS 0 0 −α 0 βS 0 0 κ 0 −µ0 0 −(µ + τ) − [ (µ+γ)(µ+π+τ) π + τ ] −βS 0  . Thus the first two eigenvalues of J(P0) are λ1 = −(µ + π + τ) and λ2 = − [ (µ + γ)(µ + π + τ) π + τ ] . While the remaining two eigenvalues are for the following 2×2 sub-matrix of J(P0). J0(P0) = [ −α βS 0 κ −µ0 ] . Now, the trace and the determinant of this matrix are respec- tively given by tr(J(P0)) = −(α + µ0) < 0 and det(P0) = αµ0 − κβ(µ + τ)(µ + γ) (µ + γ)(µ + π + τ) + τπ = αµ0 ( 1 − κβ(µ + τ)(µ + γ) αµ0[(µ + γ)(µ + π + τ) + τπ] ) = αµ0(1 −Rv). Thus the determinant is greater than zero if Rv < 1. Hence, all the four eigenvalues of J(P0) have negative real parts and so, the DFE, P0 of the system (2) is locally asymptotically stable whenever Rv < 1. 3.1.2. Global Stability of the DFE Using Lyapunov Principle method in conjunction with LaSalle’s Invariant Principles [19], we establish the global asymptotic stability of the DFE as presented in the following theorem. Theorem 3. If Rv 6 1, then the DFE, P0 of model (2) is glob- ally asymptotically stable on Ω. Proof. We define a Lyapunov function by H(I, O) = I(t) + β(µ + τ)(µ + γ) µ0[(µ + π + τ)(µ + γ) + πτ] O(t). Then, the time derivative of V along the solutions of system (2) is Ḣ(I, O) = İ(t) + β(µ + τ)(µ + γ) µ0[(µ + π + τ)(µ + γ) + πτ] Ȯ(t) = βS O −αI + β(µ + τ)(µ + γ) µ0[(µ + π + τ)(µ + γ) + πτ] (κI −µ0O) = βO ( S − (µ + τ)(µ + γ) (µ + π + τ)(µ + γ) + πτ ) + α(Rv − 1)I = βO(S − S 0) + α(Rv − 1)I. Then, Ḣ 6 0 since S 6 S 0 and Rv ≤ 1, with V̇ = 0 if and only if S = S 0 and Rv = 1. Hence, V is a Lyapunov function on Ω and so, S → (µ+τ)(µ+γ)(µ+π+τ)(µ+γ)+πτ as t →∞. Therefore, the maximum invariant set in {(S, I, O, V ) ∈ Ω0 : Ḣ 6 0}, is a singleton set P0. It follows from LaSalle’s Invariant Principles [19] that P0 is globally asymptotically stabile when Rv 6 1. 3.2. Endemic Equilibrium (EE) In the presence of disease infection, the infected compart- ments (I, O) are non empty and so, the disease invades the pop- ulation. Let Pe = (S e, Ie, Oe, Ve) be an endemic equilibrium point of model (2), then setting the right hand side of model (2) to zero we obtain after algebraic manipulations. S e = αµ0 κβ , Ie = αµ0[(µ + π + τ) + τπ](Rv − 1) κβ(µ + γ)(α + µ + τ) , Ve = απµ0 κβ(µ + γ) , Oe = αµ0[(µ + π + τ) + τπ](Rv − 1) βµ0(µ + γ)(α + µ + τ) . Hence, Pe is the unique endemic equilibrium point of sys- tem (2) which exist if and only if Rv > 1. 21 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 22 3.2.1. Local stability analysis of the EEP Theorem 4. The EEP, Pe of model (2) is locally asymptotically stable in the interior of Ω if Rv > 1. Proof. Evaluating the Jacobian matrix of system (2) at Pe gives J(Pe) =  −a −(µ + τ) −αµ0 κ −τ b −α αµ0 κ 0 π 0 0 −µ 0 κ −µ0 0  , with a = (µ+π+τ+βOe) and b = βOe. Then we use elementary row operations to obtain J(Pe) =  −a −(µ + τ) −αµ0 κ −τ 0 −[ αab + (µ + τ)] αµ0 κ [ ab − 1] 0 0 −(µ + τ) −αµ0 κ −( aµ π + τ) 0 κ −µ0 0  . Thus the eigenvalues of J(Pe) are λ1 = −a λ2 = −( aµ π + τ), λ3 = −[ αa b + (µ + τ)], and λ4 = − µ0 κ (µ + τ + α). Therefore, all the four eigenvalues are real and negative and so, the EEP, Pe is locally asymptotically stable if Rv > 1. . 3.2.2. Global Stability of Endemic Equilibrium Theorem 5. If Rv > 1, then there exist an endemic equilibrium point Pe and it is globally asymptotically stable in the interior of Ω. Proof. : Given that Rv > 1, then the existence of EEP is guaran- teed (see Section 3.4). Consider the Lyapunov function below: L(S, I, O, V ) = S − S e − S e ln ( S S e ) + I − Ie − Ie ln ( I Ie ) + V − Ve − Ve ln ( V Ve ) + O − Oe − Oe ln ( O Oe ) , The time derivative of L along the solution curve of the system (2) yields L̇ = ( 1 − S e S ) Ṡ + ( (1 − Ve V ) V̇ + ( 1 − Ie I ) İ + ( 1 − Oe O ) Ȯ = ( 1 − S e S ) [(µ + τ)(1 − S − I) −τV −βS O −πS ] + ( 1 − Ie I ) (βS O −αI) + ( (1 − Ve V ) [πS − (µ + γ)V ] + ( 1 − Oe O ) (κI −µ0O) (3) But at endemic state, we have µ + τ = (µ + τ)(S e + Ie) + τVe + βS eOe + πS e αIe = βOeS e, πS e = (µ + γ)Ve, κIe = µ0Oe. (4) Then using (4) in (3), we obtain alter algebraic manipulations L̇ = (µ + τ)S e) ( 2 − S S e − S e S ) + πS e ( 2 − S S e − S e S ) + (µ + τ)Ie ( 1 − S e S ) − (µ + τ)I ( 1 − S e S ) + τVe ( 1 − S e S ) −τV ( 1 − S e S ) + βOeS e ( 1 − S e S ) + βOS e ( 1 − S S e ) + βOS ( 1 − Ie I ) + βOeS e ( 1 − I Ie ) + πS ( 1 − V Ve ) + πS e ( 1 − Ve V ) + κI ( 1 − Oe O ) + κIe ( 1 − O Oe ) . (5) Since S ≤ S e, I ≤ Ie, V ≤ Ve and O ≤ Oe, then (5) becomes L̇ ≤ (µ + τ)S e) ( 2 − S S e − S e S ) + πS e ( 2 − S S e − S e S ) + βOeS e ( 2 − S S e − S e S ) + βOeS e ( 2 − I Ie − Ie I ) + πS e ( 2 − V Ve − Ve V ) + κIe ( 2 − O Oe − Oe O ) = S e(µ + τ + π + βOe) ( 2 − S S e − S e S ) + βOeS e ( 2 − I Ie − Ie I ) + πS e ( 2 − V Ve − Ve V ) + κIe ( 2 − O Oe − Oe O ) . (6) It follows from arithmetic-geometric inequality that( 2 − S S e − S e S ) ≤ 0, ( 2 − I Ie − Ie I ) ≤ 0,( 2 − V Ve − Ve V ) ≤ 0, ( 2 − O Oe − Oe O ) ≤ 0, and so, L̇ ≤ 0. Then we conclude this theorem using similar argument as in the proof of Theorem 3. 4. Sensitivity analysis One can observe from Section 3, that the disease persis- tence or otherwise in a community is determined by the thresh- old parameter Rv. Moreover, the uniqueness properties of both disease-free (P0) and endemic (Pe) equilibria indicate that there is an exchange of stability between P0 and Pe when Rv = 1 (i.e., P0 undergoes a transcritical bifurcation at Rv = 1, see for instance [23, Theorem 4]). However, it is instructive to investi- gate the model parameter that has greatest effect on this thresh- old parameter, thereby having same effect on the toxoplasmosis disease persistence. To achieve that, we can apply sensitivity indices to measure the relative change of a state variable with a changing parameter [17, 23]. Then using the definition in [17], the normalised forward sensitivity index with respect to each of the parameters of Rv is 22 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 23 Figure 2: Sensitivity analysis of the reproductive number Rv, showing changes in the parameters (a) π, (b) β, (c) κ and (d) µ0. Parameter values used are as presented in Table 1. the following. Υ Rv κ = ∂Rv ∂κ × κ Rv = +1 Υ Rv β = ∂Rv ∂β × β Rv = +1 Υ Rv µ = ∂Rv ∂µ × µ Rv = πµ[αµ0(µ + γ) + τ][αµ0(µ + π + τ)(µ0 + γ) + πτ] (µ + τ)[αµ0(µ + γ)(µ + π + γ) + πτ]2 Υ Rv τ = ∂Rv ∂τ × τ Rv = πτ[αµ0(µ + γ) −µ][αµ0(µ + π + τ)(µ0 + γ) + πτ] (µ + τ)[αµ0(µ + γ)(µ + π + γ) + πτ]2 Υ Rv γ = ∂Rv ∂γ × γ Rv = γπτ[αµ0(µ + γ)(µ + π + τ) + πτ] (µ + γ)[αµ0(µ + γ)(µ + π + γ) + πτ]2 Υ Rv α = ∂Rv ∂α × α Rv = − κβµ0(µ + τ)(µ + γ)2(µ + π + τ) [αµ0(µ + γ)(µ + π + τ) + πτ]2 Υ Rv µ0 = ∂Rv ∂µ0 × µ0 Rv = − αµ0(µ + γ)(µ + π + τ) αµ0(µ + γ)(µ + π + τ) + πτ Υ Rv π = ∂Rv ∂π × κ Rv = − π[αµ0(µ + γ) + τ] αµ0(µ + γ)(µ + π + τ) + πτ (7) Figure 3: Sensitivity analysis of the reproductive number Rv, showing changes in the parameters (e) γ, (f) µ, (g) τ and (h) α. Parameter values used are as presented in Table 1. It can be observed from equation (7) that most of the expres- sions for the sensitivity indices are complex and so, we evaluate these indices at the parameter values presented in Table 1. The 23 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 24 Table 2: Sensitivity indices of Rv to parameters for the toxoplasmosis model, evaluated at the parameter values presented in Table 1. The parameters are ordered from most sensitive to least one. parameter sensitivity index 1. π −644.49655 2. β +1.00000 3. κ +1.00000 4. µ0 −1.00000 5. γ 0.94230 6. µ 0.58798 7. τ −0.58650 8. α −0.00003 resulting sensitivity indices of Rv to the eight model parameters from the most sensitive to the least one are shown in Table 2. We infer from Figures 2 and 3 that the threshold parameter Rv increases with increasing values of β,κ,µ and τ. While it decreases with increasing values of π,µ0 and α. As expected the vaccination parameter π is more sensitive to changes than any other parameter of the model. 5. Discussion and conclusion In this note, we extend the deterministic model for the trans- mission dynamics of toxoplasmosis in cats population with vac- cination by separating the vaccinated and the recovered com- partments. The stability analysis of the model equilibria are car- ried out and we show that the unique disease-free equilibrium as well as the unique endemic state of the model are asymptot- ically stable globally under certain condition. That is, the DFE, P0 is stable globally asymptotically when Rv < 1, while the EEP, Pe is globally asymptotically stable if Rv > 1. In order to complete our investigation, we use the sensitivity indices as pre- sented in [17] to determine the most sensitive model parameter. As expected, the vaccination parameter, π is the most sensitive to changes than any other parameter of the model. Moreover, considering the complex nature of the expressions for the sen- sitivity indices in equation (7), we evaluate these indices at the parameter values presented in Table 1. Then, we order them from the most sensitive to the least one as presented in Table 2. Furthermore, the numerical simulations in Figures 2 and 3 reveal how the threshold parameter Rv varies with varying val- ues of each of the model parameters. In fact, Rv increases with increasing values of the parameters β,κ,µ and τ and decreases with increasing values of π,µ0 and α. In conclusion, the sensitivity analysis indicates that the best way to contain the spread of toxoplasmosis is either to increase the vaccination rate of susceptible cats or to increase the death rate of the oocysts by environmental sanitation. It can be ob- served from Figure 2 that vaccination is the best control strat- egy than the removal of the oocysts. This is similar to the result presented in [2]. References [1] A. Gautam & A. Priyadarshi, “Mathematical Modelling of Toxoplasma Gondii and Host Immune Response”, AIP Conference Proceedings 1975 (2018) 030002, doi: 10.1063/1.5042172. [2] A. J. Arenas, G. González-Parra & R. J. “Villanueva Micó, Modeling toxoplasmosis spread in cat populations under vaccination”, Theor. Popul. Biol. 77 (2010) 227, https://doi.org/10.1016/j.tpb.2010.03.005. [3] C. Castillo-Chavez & B. Song, “Dynamical models of tuberculo- sis and their applications”, Math. Biosci. Engin. 1 (2004) 361, doi: 10.3934/mbe.2004.1.361. [4] C. Peña & K. G. Hermes Martinez, “Hybrid model of the spread of Toxo- plasmosis between two Town of Colombia”, Tecciencia 7 (2015) 1, DOI: http:/dx.doi.org/10.18180/tecciencia.2015.18.1. [5] C. Ramakrishnan, S. Maier, R. A. Walker, H. Rehrauer, et al., “An experimental genetically attenuated live vaccine to prevent transmis- sion of Toxoplasma gondii by cats”, Scientific Reports 9 (2019) 1479, https://doi.org/10.1038/s41598-018-37671-8. [6] D. F. Aranda, R. J. Villanueva, A. J. Arenas & G. C. Gonzalez- Parra, “Mathematical modeling of Toxoplasmosis disease in varying size populations”, Comp. Math. Appl. 56 (2008) 690, https://doi.org/10.1016/j.camwa.2008.01.008. [7] D. L. Zulpo, A. S. Sammi, J. R. dos Santos, et al., “Toxoplasma gondii: A study of oocyst re-shedding in domestic cats”, J. Vet. Parasitol. 249 (2018) 17. [8] E. A. Innes, P. M. Bartley, S. Maley, F. Katzer & D. Buxton, “Veterinary vaccines against Toxoplasma gondii, Mem Inst Oswaldo Cruz”, Rio de Janeiro 104 (2009) 246. [9] G. C. González-Parra, A. J. Arenas, D. F. Aranda, R. J. Villanueva & L. Jódar, “Dynamics of a model of Toxoplasmosis disease in human and cat populations”, Comp.Math. Appl. 57 (2009) 1692, https://doi.org/10.1016/j.camwa.2008.09.012. [10] J. A. León Marín & I. D. Gandica, “A mathematical model for the repro- duction dynamics of toxoplasma gondii”, J. Biol. Syst. 23 (2015) S91, doi: 10.1142/S0218339015400082. [11] J. D. Ferreira, L. M. Echeverry & C. A. Peña Rincon, “Stability and bifur- cation in epidemic models describing the transmission of toxoplasmosis in human and cat populations”, Math. Meth. Appl. Sci. 40 (2017) 5575, doi: 10.1002/mma.4410. [12] J. P. Dubey, “The history of Toxoplasma gondii: The first 100 years”, J. Eukaryot. Microbiol. 55 (2008) 467, doi: 10.1111/j.1550- 7408.2008.00345.x. [13] K. M. Addo, “SEIR Model for Dogs Rabies. A case study: Bango District, Ghana”, Master’s thesis, Kwame Nkurumah University of Science and Technology, Ghana (2012). [14] K. Sornsong, S. Naowarat & P. Thongjaem, “Mathematical model for the dynamics transmission of rabies with control measures”, Austr. J. Basics and Appl. Sci. 10 (2016) 169. [15] M. Hari & Y. Zulfahmi, “Bifurcation Analysis of Toxoplasmosis Epi- demic Control on Increased Controlled Rate of Suppressing the Rate of Infected Births”, Int. J. Comp. Sci. Appl. Math. 6 (2020) 1, doi: 10.12962/j24775401.v6i1.5978. [16] M. Langlais, M. Lélu, C. Avenet & E. Gilot-Fromont, “A simplified model system for Toxoplasma gondii spread within a heterogeneous environ- ment”, Nonlinear Dyn. 68 (2012) 381, doi:10.1007/s11071-011-0255-4. [17] N. Chitnis, J. M. Hyman & J. M. Cushing, “Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model”, Bull. Math. Biol. 70 (2008) 1272, doi: 10.1007/s11538-008-9299-0. [18] O. M. Ogunmiloro, “Mathematical Modeling of the Coinfection Dynam- ics of Malaria-Toxoplasmosis in the Tropics”, Biomet. Lett. 56 (2019) 139, doi: https://doi.org/10.2478/bile-2019-0013. [19] P. Glendinning, Stability, Instability and Choas: an introduction to the theory of nonlinear differential equations, Cambridge University Press, (1994). [20] P. Van den Driessche & J. Wanmough, “Reproduction numbers and sub- threshold endemic equilibria for compartimental models of disease trans- mition”, Mathematical Biosciences 180 (2000) 29. [21] R. Verma & P. Khanna, “Development of Toxoplasma gondii vaccine: A global challenge”, Human Vac. Imm. therap. 9 (2013) 291. [22] S.M. Garba, A.B. Gumel & M.R. Abu Bakar, “Backward bifurcations in dengue transmission dynamics”, Math. Biosci. 215 (2008) 11. 24 Muhammad & Usaini / J. Nig. Soc. Phys. Sci. 3 (2021) 17–25 25 [23] S. Usaini, U. T. Mustapha & S. M. Sabiu, “Modelling scholastic under- achievement as a contagious disease”, Math. Meth. Appl. Sci. Special Issue (2018) 1, https://doi.org/10.1002/mma.4924. [24] T. Loung & S. A. Davis, “Rabies vaccination targets for stray dogs popu- lation”, Frontiers in Vet. Sci. 4 (2017) 52. [25] T. Dubie, G. Terefe, M. Asaye & T. Sisay , “Toxoplasmosis: Epidemiol- ogy with the emphasis of its public health importance” 2 (2014) 097. 25