Mathematical Modeling of Toxoplasmosis in Cats with Two Time Delays under Environmental Effects

Abstract

In this paper, we construct a more realistic mathematical model to study toxoplasmosis dynamics. The model considers two discrete time delays. The first delay is related to the latent phase, which is the time lag between when a susceptible cat has effective contact with an oocyst and when it begins to produce oocysts. The second discrete time delay is the time that elapses from when the oocysts become present in the environment to when they are able to infect. The main aim in this paper is to find the conditions under which the toxoplasmosis can disappear from the cat population and to study whether the time delays can affect the qualitative properties of the model. Thus, we investigate the impact of the combination of two discrete time delays on the toxoplasmosis dynamics. Using dynamical systems theory, we are able to find the basic reproduction number ${\mathcal{R}}_0^d$ that determines the global long-term dynamics of the toxoplasmosis. We prove that, if ${\mathcal{R}}_0^d<1$, the toxoplasmosis will be eradicated and that the toxoplasmosis-free equilibrium is globally stable. We design a Lyapunov function in order to prove the global stability of the toxoplasmosis-free equilibrium. We also prove that, if the threshold parameter ${\mathcal{R}}_0^d$ is greater than one, then there is only one toxoplasmosis-endemic equilibrium point, but the stability of this point is not theoretically proven. However, we obtained partial theoretical results and performed numerical simulations that suggest that, if ${\mathcal{R}}_0^d>1$, then the toxoplasmosis-endemic equilibrium point is globally stable. In addition, other numerical simulations were performed in order to help to support the theoretical stability results.

1. Introduction

T. gondii is a protozoan parasite that attacks and parasitizes animals including humans. Its only definitive hosts are members of the family Felidae. However, other species are intermediate hosts [1,2,3]. In humans, immunocompromised and congenitally infected persons can suffer from severe and lethal effects [4,5]. In individuals who are immunocompetent, toxoplasmosis infection usually only causes flu symptoms and chorioretinitis [4,6,7]. However, severe acute toxoplasmosis has been reported in immunocompetent persons [8]. In addition, it has been suggested that T. gondii infection might cause autoimmune diseases and neurological disorders [9,10,11,12,13]. It has been found that the genotype of the strain of T. gondii affects the severity of the infection [14,15,16].

T. gondii is one of the most significant foodborne pathogens [4,17,18,19]. In the USA alone, there are more than 40 million people with T. gondii [4,17]. It has been estimated that more than 60% of some communities around the world are affected. Prevalence is higher in regions that are hot, humid, and at lower altitudes [4,20]. For instance, in tropical South America, T. gondii has higher genetic diversity [16]. It has been estimated that 32% of the total consequences of toxoplasmosis are due to congenital toxoplasma infections [21]. Toxoplasmosis exists in all the countries around the world and it has been estimated that there are 190,000 congenital toxoplasmosis cases annually (1.5 per 1000 live births) [22,23]. However, the estimation related to acquired toxoplasmosis is uncertain. In the United States, it has been estimated that there are 58 congenital toxoplasmosis cases per 100,000 individuals [22,24]. It has also been estimated that, in the United States, the annual cost associated with T. gondii is 10,964 quality-adjusted life-years (ranking it third among food-related pathogens) and the cost of illness is USD 2.9 billion (ranked second) [22,25]. High associated costs are present in South America and in some low-income countries [22,23].

There are different routes of transmission of T. gondii; ingestion of sporulated oocysts, ingestion of viable tissue cysts, and congenital infection. For instance, the number of congenital toxoplasmosis infections has been estimated to be 1.2 million [26]. Oocyst contamination in fruits and greens contributes significantly more to total estimated T. gondii infections than bradyzoite-contaminated foods [19]. The infected cats are able to shed oocysts in the environment and humans can be infected after consumption of undercooked and raw meat [19,20,27,28,29,30]. It has been estimated that the feces of cats shedding T. gondii have $2.5\times {10}^6$ oocysts/gr and a cat sheds around 20 million oocysts per day [31]. The control of toxoplasmosis and the definitive host is a challenging problem [30,32,33,34,35]. Current health interventions are weak since they only focus on acute/reactivated forms of the disease [5,30]. Thus, studies related to toxoplasmosis interventions are crucial and worthwhile.

Animals usually do not suffer from symptoms related to toxoplasmosis infection [36,37,38,39]. However, they suffer from congenital transmission and its consequences, such as offspring with abnormalities [30]. Documentation of the consequences of T. gondii infection in animals has been presented [4,36,37,38,40,41]. Toxoplasmosis infections in many different animals have important consequences from different viewpoints [36,41,42,43,44]. Also, economic and public health impacts related to toxoplasmosis have been mentioned in the literature [30,45,46].

There have been several previous studies devoted to investigating toxoplasmosis dynamics [26,47]. Some of these works have used mathematical models with a variety of approaches [26]. Few mathematical models have focused on investigating cat–human transmission under different assumptions related to demographic factors [48,49]. However, there are other, more complex models that have been employed to study the dynamics of toxoplasmosis using other variables and hosts [47,50,51,52,53]. In addition, some researchers have investigated and explored the effect of cat and swine vaccinations [32,33,47,50,51,54]. Also, within-host mathematical models have been proposed to study the T. gondii dynamics within the host [55]. In [26], an excellent overview of mathematical models for toxoplasmosis dynamics is provided.

In this paper, we study the effects of the combination of two time delays on the dynamics of toxoplasmosis, taking into account the cat population and the number of oocysts in the environment. The introduction of two time delays allow us to study a more realistic scenario despite the typical limitations of mathematical models for real-world processes. One time delay is related to the cats’ exposed stage and the second delay is related to the time it takes an oocyst to become infectious after being shed by the cat. In particular, the first discrete time delay occurs due to the fact that, when a susceptible cat has effective contact with an oocyst, that cat starts an incubation period before being able to start shedding oocysts. Thus, the time delay is the incubation period of the cat. In order to study the aforementioned effects, we construct a mathematical model that uses a nonlinear system of delay differential equations. One key aspect of the model is the fact that, once the cats have effective contact with the oocysts, they are not able to produce oocysts instantaneously [2,6,20]. The inclusion of oocysts is crucial since they are highly responsible for the persistence of the T. gondii [2,6,20,27,29,56,57]. The constructed model encompasses the fact that cats can get the T. gondii parasite through effective contact with an oocyst. The model is more realistic than the mathematical model presented in [47]. This previous model does not include any time delay related to more realistic scenarios. There is a delay in the shedding of oocysts that lasts between 3 and 30 days after ingestion of tissue cysts [6,58,59]. The constructed model also includes the cats due to their importance in the reproduction of T. gondii and therefore for the toxoplasmosis disease [60]. Public health policies should take into account many factors to control the spread of toxoplasmosis; for instance, domestic cats and other intermediate hosts [43]. There are three major genotypes (type one, type two, and type three) of T gondii. In Europe and the USA, type two is the most prevalent for congenital toxoplasmosis. Their pathogenicity differs, but if we assume an average for all the types, then the model proposed in this study will be applicable. However, if we differentiate the genotypes and their features regarding pathogenicity then a more complex model will be required [61,62,63].

Mathematical models based on systems of delay differential equations have been used to study different kind of infectious diseases, including between hosts and within hosts [64,65,66,67,68,69,70,71,72,73]. The reason for including time delays is the fact that many processes do not have instantaneous effects in the real world. However, there are some drawbacks of employing mathematical models with delay differential equations. For instance, their analysis is usually more challenging and cumbersome. Furthermore, solving delay differential equations is more difficult than ordinary differential equations [74,75,76,77]. The inclusion of delays in the differential equations can have a great impact on the stability of the steady states and oscillatory behavior can appear with some particular values of the time delays [78,79,80,81,82,83,84,85]. Moreover, the inclusion of multiple time delays makes the mathematical analysis more complex and oftentimes intractable [86,87,88,89,90,91]. In particular, there are several works that have studied mathematical models that have two time delays [92,93,94]. For instance in [95], the authors used two time delays in a model to study the control of disease spreading in networks.

In summary, in our study, one main objective is to investigate the effect of the combination of two time delays on the toxoplasmosis infection dynamics. The importance of the combination of two time delays resides in the fact that, if a stability switch occurs in the model with only one time delay, then it might be possible that the result does not hold for the model with two time delays (see [81,96]). Moreover, the model with two time delays is closer to reality and also provides helpful insight into the mechanism of more complex dynamics [92]. Therefore, we investigate if the time delays are able to produce different behaviors and complexity in the toxoplasmosis dynamics. We explore the stability of the steady states of the model. We prove that, if ${\mathcal{R}}_0^d$ is greater than one, then there is only one toxoplasmosis-endemic equilibrium point, but due to the presence of two discrete time delays, we are not able to prove its stability. However, we obtained partial theoretical results and performed numerical simulations that suggest that, if ${\mathcal{R}}_0^d>1$, then the point is globally stable. We also performed in silico simulations to support the theoretical results. In addition, we simulated some scenarios to gain insight into the effect of the two time delays on the toxoplasmosis dynamics.

2. Mathematical Model

We designed a new mathematical model to explore and investigate the dynamics of toxoplasmosis in a population of cats. The proposed model includes two discrete time delays, which makes the model closer to reality. The first delay is related to the time that it takes for a cat to become infected after infectious contact with the oocysts. In particular, the first discrete time delay is the incubation period before the cats start shedding oocysts. Thus, the term $\beta S\left(t\right)O\left(t\right)$ represents that the cats have effective contact with the oocysts at time t and the term $-\beta S(t-{\tau }_1)O(t-{\tau }_1)e^{-\mu {\tau }_1}$ represents the cats leaving the incubation phase (after a time ${\tau }_1$) and becoming able to shed oocysts [69,97]. The second time delay is related to the time that it takes the oocysts to sporulate and become infectious (somewhere between 1 and 5 days) after they are shed in the environment by an infected cat [6,20,56,58]. The mathematical model includes vaccinated cats and infected cats [32,47,98]. The cats get infected by contact with oocysts. These infections depend on the amount of T. gondii in the environment [54]. The model uses a system of delay differential equations. Further details regarding the model can be seen in [47].

The model has five subpopulations: susceptible $S\left(t\right)$, latent $E\left(t\right)$, infected $I\left(t\right)$, vaccinated–recovered $V_R\left(t\right)$, and oocysts $O\left(t\right)$. The vaccine provides life immunity [32,33,99,100]. A susceptible cat goes to the latent stage after contact with an oocyst. The infected cats stay in this stage for a mean time of $\frac{1}{\alpha }$. The change in the number of oocysts in the environment is proportional to the infectious cats $I\left(t\right)$. The oocysts in the environment suffer from natural degradation.

The mathematical model is given in terms of the following nonlinear delay differential equation system

$$\begin{array}{ccc}\dot{S}\left(t\right)& \hfill =& \mu (1-I(t\left)\right)-\beta S\left(t\right)O\left(t\right)-(\mu +\gamma )S\left(t\right),\hfill \\ \dot{E}\left(t\right)& \hfill =& \beta S\left(t\right)O\left(t\right)-\beta S(t-{\tau }_1)O(t-{\tau }_1)e^{-\mu {\tau }_1}-\mu E\left(t\right),\hfill \\ \dot{I}\left(t\right)& \hfill =& \beta S(t-{\tau }_1)O(t-{\tau }_1)e^{-\mu {\tau }_1}-\alpha I\left(t\right),\hfill \\ \dot{V_R}\left(t\right)& \hfill =& \alpha I\left(t\right)+\gamma S\left(t\right)-\mu V_R\left(t\right),\hfill \\ \dot{O}\left(t\right)& \hfill =& kI(t-{\tau }_2)-{\mu }_{0}O\left(t\right).\hfill \end{array}$$

(1)

We assume that the total population of cats is given by $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+V_R\left(t\right)=1$ (without loss of generality). Therefore, since the state variable $V_R$ is not included explicitly in the first three parts of model (1), we can analyze the model without this variable. The graphical representation of mathematical model (1) can be observed in Figure 1. Note that vertical transmission is included in the cat population [28,40,101,102,103,104].

For the sake of simplicity, model (1) can be studied using the following simplified system

$$\begin{array}{ccc}\dot{S}\left(t\right)& \hfill =& \mu (1-I(t\left)\right)-\beta S\left(t\right)O\left(t\right)-(\mu +\gamma )S\left(t\right),\hfill \\ \dot{E}\left(t\right)& \hfill =& \beta S\left(t\right)O\left(t\right)-\beta S(t-{\tau }_1)O(t-{\tau }_1)e^{-\mu {\tau }_1}-\mu E\left(t\right),\hfill \\ \dot{I}\left(t\right)& \hfill =& \beta S(t-{\tau }_1)O(t-{\tau }_1)e^{-\mu {\tau }_1}-\alpha I\left(t\right),\hfill \\ \dot{O}\left(t\right)& \hfill =& kI(t-{\tau }_2)-{\mu }_{0}O\left(t\right),\hfill \end{array}$$

(2)

where the parameters $\mu ,\gamma ,{\tau }_1,{\tau }_2,\beta ,\alpha ,k$ are all positive. Consider

$${\mathbb{R}}_{+}^4=\left\{(x_1,x_2,x_3,x_4):x_i\ge 0,i=1,\cdots ,4\right\}.$$

In order to investigate the existence and uniqueness of model (2), we denote by $C\left[[-\tau ,0],{\mathbb{R}}_{+}^4\right]$ (with $\tau =max\left\{{\tau }_1,{\tau }_2\right\}$) the Banach space of continuous functions $\zeta$ defined in $\zeta :[-\tau ,0]\to {\mathbb{R}}_{+}^4$ such that $\zeta \left(t\right)=\left({\zeta }_1\left(t\right),{\zeta }_2\left(t\right),{\zeta }_3\left(t\right),{\zeta }_4\left(t\right)\right)$ and a norm defined by (see [105])

$$\begin{array}{c}\hfill ∥\zeta ∥=\underset{-\tau \le s\le 0}{sup}\left\{\left|{\zeta }_1\left(s\right)\right|,\left|{\zeta }_2\left(s\right)\right|,\left|{\zeta }_3\left(s\right)\right|\left|{\zeta }_4\left(s\right)\right|\right\}.\end{array}$$

Thus, for biological significance, we further assume that the initial conditions of model (2) have the following structure:

$$\begin{array}{l}S\left(s\right)={\zeta }_1\left(s\right)>0,E\left(s\right)={\zeta }_2\left(s\right)\ge 0,I\left(s\right)={\zeta }_3\left(s\right)\ge 0,O\left(s\right)={\zeta }_4\left(s\right)\ge 0,\forall s\in [-\tau ,0],\\ S\left(0\right)=S_0={\zeta }_1\left(0\right),E\left(0\right)=E_0={\zeta }_2\left(0\right),I\left(0\right)=I_0={\zeta }_3\left(0\right),O\left(0\right)=O_0={\zeta }_4\left(0\right).\end{array}$$

(3)

Moreover, for the continuity of the initial conditions, it is required that

$$\begin{array}{c}\hfill E\left(0\right)={\int }_{-\tau }^0{\zeta }_1\left(w\right){\zeta }_4\left(w\right)e^{\mu w}dw.\end{array}$$

The initial value of $E\left(0\right)$ must be related to the initial values of the variables S and O such that they satisfy the conditions in Equation (3), which is an integral constraint of the function ${\zeta }_2:[-\tau ,0]\to [0,\infty ],$ and this is important so that the positivity of the solutions is preserved [106]. From the theory of functional differential equations [107], we can apply the theorem that comes next.

Theorem 1.

Suppose $\mathcal{O}$ is an open set in $\mathcal{R}\times C\left([-\tau ,0],{\mathbb{R}}_{+}^4\right),$ $f:\mathsf{\Omega }\to {\mathcal{R}}^4$ is continuous, and $f(t,\zeta )$ is Lipschitzian in ζ in each compact set in $\mathsf{\Omega }.$ If $(\sigma ,\phi )\in \mathsf{\Omega }$, then there is a unique solution to model (2) through $(\sigma ,\phi ).$

Proof. 

See [107] (p. 44). □

Therefore, it follows that model (2) admits a unique solution provided that the conditions shown in Equation (3) hold.

Next, the positivity and boundedness for model (2) are established with the subsequent theorem.

Theorem 2.

The solution $X\left(t\right)=\left(S\right(t),E(t),I(t),O(t\left)\right)$ to Equation (2) is positive and uniformly bounded on $[0,+\infty )$ if Equation (3) holds.

Proof. 

From Equation (2), one can find that

$$\begin{array}{c}\hfill \dot{S}\left(t\right)=\mu (1-I\left(t\right))-\beta S\left(t\right)O\left(t\right)-\gamma S\left(t\right)\ge -\beta S\left(t\right)O\left(t\right)-\gamma S\left(t\right).\end{array}$$

Thus,

$$\begin{array}{c}\hfill S\left(t\right)\ge S\left(0\right)exp\left(-{\int }_0^t(\beta O\left(s\right)+\gamma )ds\right)>0,\forall t\ge 0.\end{array}$$

To show that the other variables preserve positivity over time, we proceed using the method of steps [106]. With this order of ideas, let $t\in [0,\tau ].$ With the initial conditions shown in Equation (3), it follows that $S(t-\tau )O(t-\tau )\ge 0$ in this interval. Then, from the third equation of system (2), and using the comparison criterion for differential equations, it follows that $I\left(t\right)\ge 0$ for all $t\in (0,\tau ].$ Analogously, it is verified that $O\left(t\right)>0$ with $t\in (0,\tau ],$ and, therefore, the non-negativity for $E\left(t\right)$ is a consequence of the simple fact that the solution of the differential equation related to $E\left(t\right)$ can be written as

$$\begin{array}{c}\hfill E\left(t\right)=e^{-\mu t}\left(E\left(0\right)+{\int }_{t-\tau }^{t}S\left(w\right)O\left(w\right)e^{\mu w}dw\right),\forall t\in (0,\tau ].\end{array}$$

The above reasoning can be applied for variables $E,I,O$ in the interval $(\tau ,2\tau ]$ and successively in intervals of the form $(n\tau ,(n+1\left)\tau \right],$ for $n\in {\mathbb{N}}_0$. Next, from system (2), one gets

$$\begin{array}{c}\hfill \dot{S}\left(t\right)+\dot{E}\left(t\right)+\dot{I}\left(t\right)+\dot{O}\left(t\right)\le \mu +k-M(S\left(t\right)+E\left(t\right)+I\left(t\right)+O\left(t\right)),\end{array}$$

where $M<min\left\{\mu +\gamma ,\mu ,\alpha +\mu ,{\mu }_0\right\}.$ As a result, one obtains

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; sup}(S\left(t\right)+E\left(t\right)+I\left(t\right)+O\left(t\right))\le {\displaystyle \frac{\mu +k}{M}}.\end{array}$$

Hence, we study the qualitative behavior of model (2) in the following region:

$$\mathcal{R}=\left\{(S,E,I,O)\in {\mathbb{R}}_{+}^4:0<S\left(t\right)+E\left(t\right)+I\left(t\right)+O\left(t\right)\le {\displaystyle \frac{\mu +k}{M}}\right\},$$

where ${\mathbb{R}}_{+}^4$ denotes the non-negative cone of ${\mathbb{R}}^4$. Thus, the set $\mathcal{R}$ is positively invariant. □

Now, we can focus on the stability analysis of model (2) in the restricted region $\mathcal{R}$.

3. Stability Analysis of the Model

The qualitative analysis of the mathematical model from Equation (2) can be undertaken by using linearizations and finding suitable Lyapunov functions [105,108]. Before we begin the stability analysis, we would like to prove the positivity and boundedness of the solutions to model (2).

3.1. Disease-Free Steady State for the Model without Delay

The steady states of model (2) without delay are the same as for the model with the two time delays ${\tau }_1$ and ${\tau }_2$. The steady states are crucial for the study of dynamical systems [108,109,110].

We can start analyzing the stability of the steady states of the model without considering the delays. This analysis has been performed previously and the details can be seen in [47,111,112]. The model from Equation (2) without delay has a toxoplasmosis-free equilibrium point $F^{*}$ and a toxoplasmosis-endemic equilibrium point. The toxoplasmosis-free equilibrium point is given by

$$F^{*}=\left(\frac{\mu }{\mu +\gamma },0,0,0\right).$$

(4)

It has also been found that the basic reproduction number for model (2) without any time delays is given by [47,52,111]:

$${\mathcal{R}}_0=\sqrt{\frac{k\beta \mu }{\alpha {\mu }_0\left(\gamma +\mu \right)}}.$$

(5)

If ${\mathcal{R}}_0<1$, the toxoplasmosis-free equilibrium point is globally asymptotically stable. On the other hand, the toxoplasmosis-free equilibrium point $F^{*}$ is unstable when ${\mathcal{R}}_0>1$ [47,52,53].

3.2. Endemic Steady State

Model (2) also has a unique toxoplasmosis-endemic equilibrium point. It is important to determine this endemic point for the stability analysis of model (2). Indeed, the toxoplasmosis-endemic equilibrium point is given by $E{P}^{*}=\left(S^{*},E^{*},I^{*},O^{*}\right)$, where

$$\begin{array}{l}{S}^{*}=\frac{\alpha {\mu }_0}{\beta k{\mathrm{e}}^{-\mu {\tau }_1}},\\ E^{*}=\frac{\alpha \left(-\mu k\beta {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \gamma {\mu }_0+\alpha {\mu }_0\mu \right)\left({\mathrm{e}}^{-\mu {\tau }_1}-1\right)}{\beta k{\mathrm{e}}^{-\mu {\tau }_1}\left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)\mu }=\frac{{\alpha }^2{\mu }_0(\gamma +\mu )\left(1-{\mathrm{e}}^{-\mu {\tau }_1}\right)\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)}{\beta k{\mathrm{e}}^{-\mu {\tau }_1}\left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)\mu },\\ I^{*}=-\frac{-\mu k\beta {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \gamma {\mu }_0+\alpha {\mu }_0\mu }{\beta \left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)k}=\frac{\alpha {\mu }_0(\gamma +\mu )\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)}{\beta \left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)k},\\ O^{*}=-\frac{-\mu k\beta {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \gamma {\mu }_0+\alpha {\mu }_0\mu }{\beta {\mu }_0\left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)}=\frac{\alpha (\gamma +\mu )\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)}{\beta \left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)},\end{array}$$

(6)

and where we define

$${\mathcal{R}}_0^d=\sqrt{\frac{\mu k\beta {\mathrm{e}}^{-\mu {\tau }_1}}{{\mu }_0\left(\gamma +\mu \right)\alpha }}.$$

(7)

Thus, this toxoplasmosis-endemic equilibrium makes biological sense if ${\mathcal{R}}_0^d>1$. Moreover, note that $E{P}^{*}$ is biologically feasible iff ${\mathcal{R}}_0^d\ge 1$. In particular, if ${\mathcal{R}}_0^d=1$, the point $E{P}^{*}$ collides with $F^{*}$. Now, we proceed to study model (2), which contains two discrete time delays ${\tau }_1$ and ${\tau }_2$.

3.3. Toxoplasmosis-Free Steady State Analysis

The steady states of model (2) can be found by setting the right-hand side as equal to zero and setting any state variable with delay as a state variable without delay. To study the stability of the toxoplasmosis-free steady state, we linearize system (2). The following is the characteristic equation:

$$|J_0+J_{{\tau }_1}{e}^{-\lambda {\tau }_1}+J_{{\tau }_2}{e}^{-\lambda {\tau }_2}-\lambda I|=0,$$

where

$$J_0=\left[\begin{array}{cccc}-\gamma -\mu & 0& -\mu & -\frac{\beta \mu }{\gamma +\mu }\\ 0& -\mu & 0& \frac{\beta \mu }{\gamma +\mu }\\ 0& 0& -\alpha & 0\\ 0& 0& 0& -{\mu }_0\end{array}\right],J_{{\tau }_1}=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& -\frac{\beta \mu {\mathrm{e}}^{-\mu {\tau }_1}}{\gamma +\mu }\\ 0& 0& 0& \frac{\beta \mu {\mathrm{e}}^{-\mu {\tau }_1}}{\gamma +\mu }\\ 0& 0& 0& 0\end{array}\right],J_{{\tau }_2}=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& k& 0\end{array}\right]$$

As a result, the transcendental characteristic equation is given by

$$|J_0+J_{{\tau }_1}{e}^{-\lambda {\tau }_1}+J_{{\tau }_2}{e}^{-\lambda {\tau }_2}-\lambda I|=\left(\lambda +\mu \right)\left(\lambda +\gamma +\mu \right)T\left(\lambda \right)=0,$$

(8)

where,

$$\begin{array}{ll}T\left(\lambda \right)& =\left(\gamma +\mu \right){\lambda }^2+\left(\gamma \alpha +\alpha \mu +\gamma {\mu }_0+\mu {\mu }_0\right)\lambda -{\mathrm{e}}^{-\lambda {\tau }_1}\beta \mu {\mathrm{e}}^{-\mu {\tau }_1}{\mathrm{e}}^{-\lambda {\tau }_2}k+\alpha \gamma {\mu }_0+\alpha {\mu }_0\mu \\ & =\left(\gamma +\mu \right){\lambda }^2+\left(\gamma +\mu \right)\left({\mu }_0+\alpha \right)\lambda -{\mathrm{e}}^{-\lambda {\tau }_1}\beta \mu {\mathrm{e}}^{-\mu {\tau }_1}{\mathrm{e}}^{-\lambda {\tau }_2}k+\alpha {\mu }_0\left(\gamma +\mu \right)\\ & =\left(\gamma +\mu \right)\left({\lambda }^2+\left({\mu }_0+\alpha \right)\lambda -\alpha {\mu }_0{\left({\mathcal{R}}_0^d\right)}^2{e}^{-\lambda ({\tau }_1+{\tau }_2)}+\alpha {\mu }_0\right).\end{array}$$

It is clear from Equation (8) that two eigenvalues are ${\lambda }_1=-\mu$ and ${\lambda }_2=-\mu -\gamma$. The other eigenvalues are given by the roots of $T\left(\lambda \right)$. We propose the following theorem.

Theorem 3.

The toxoplasmosis-free equilibrium $F^{*}$ of model (2) is locally asymptotically stable if ${\mathcal{R}}_0^d<1$ but unstable if ${\mathcal{R}}_0^d>1$.

Proof. 

Suppose that one eigenvalue $\lambda$ has a non-negative real part. Then, $T\left(\lambda \right)=0$ can be written as

$${\lambda }^2+\lambda (\alpha +{\mu }_0)+\alpha {\mu }_0=\alpha {\mu }_0{\left({\mathcal{R}}_0^d\right)}^2{\mathrm{e}}^{-\lambda \left({\tau }_1+{\tau }_2\right)}.$$

(9)

The modulus of the left-hand side of Equation (9) satisfies

$$|{\lambda }^2+\lambda (\alpha +{\mu }_0)+\alpha {\mu }_0| \ge \alpha {\mu }_0.$$

(10)

On the other hand, the modulus of the right-hand side of Equation (9) satisfies

$$|\alpha {\mu }_0{\left({\mathcal{R}}_0^d\right)}^2{\mathrm{e}}^{-\lambda \left({\tau }_1+{\tau }_2\right)}| <\alpha {\mu }_0{\left({\mathcal{R}}_0^d\right)}^2<\alpha {\mu }_0,$$

(11)

if ${\mathcal{R}}_0^d<1$. Then, Equation (9) with inequalities (10) and (11) lead to a contradiction. It follows that all the roots of the characteristic Equation (8) have a negative real part. Thus, the toxoplasmosis-free equilibrium point is locally asymptotically stable whenever ${\mathcal{R}}_0^d<1$. □

Now that we have obtained the local stability of the equilibrium point $F^{*}$, we can proceed to investigate its global stability.

3.4. Global Stability Analysis of the Toxoplasmosis-Free Equilibrium $F^{*}$

We can try to investigate the global stability of $F^{*}$ by using the Lyapunov theorem or the LaSalle variance principle [113]. Let us propose the following theorem.

Theorem 4.

The toxoplasmosis-free equilibrium point $F^{*}$ of system (2) is globally asymptotically stable in $\mathcal{R}$ if ${\mathcal{R}}_0^d<1.$

Proof. 

Denoting $y_t$ as the translation of the solution to system (2), we have that $y_t=(S(t+\xi ),I(t+\xi ),V(t+\xi ))$ where $\xi \in [-\tau ,0]$. Let us propose the following Lyapunov function:

$$\begin{array}{c}\hfill \mathcal{L}(I\left(t\right),S\left(t\right),O\left(t\right))=\beta {\int }_{t-{\tau }_1}^{t}S\left(\xi \right)O\left(\xi \right)e^{-\mu {\tau }_1}d\xi +\alpha {\int }_{t-{\tau }_2}^{t}I\left(\xi \right)d\xi +I\left(t\right)+\frac{\alpha }{k}O\left(t\right).\end{array}$$

(12)

Note that, at $F^{*}$, one has that $S^{*}=\frac{\mu }{\mu +\gamma }, I^{*}=0, O^{*}=0$ and therefore $\mathcal{L}\left(F^{*}\right)=0$. The Lyapunov function $\mathcal{L}$ satisfies

$$\begin{array}{c}\hfill \mathcal{L}\left(I\right(t),S(t),O(t\left)\right)>0,\end{array}$$

(13)

at any point different from $F^{*}$. The time derivative of $\mathcal{L}\left(I\right(t),S(t),O(t\left)\right)$ through the trajectories of model (2) is given by

$$\begin{array}{cc}\hfill {\mathcal{L}}^{\prime }=& \beta e^{-\mu {\tau }_1}\left[S\left(t\right)O\left(t\right)-S(t-\tau )O(t-\tau )\right]+\alpha (I\left(t\right)-I(t-{\tau }_2))+I^{\prime }\left(t\right)+\frac{\alpha }{k}{O}^{\prime }\left(t\right)\hfill \\ \hfill =& \beta e^{-\mu {\tau }_1}\left[S\left(t\right)O\left(t\right)-S(t-{\tau }_1)O(t-{\tau }_1)\right]+\alpha (I\left(t\right)-I(t-{\tau }_2))+\alpha I(t-{\tau }_2)-{\displaystyle \frac{\alpha {\mu }_0}{k}}O\left(t\right)\hfill \\ \hfill =& \beta e^{-\mu {\tau }_1}S\left(t\right)O\left(t\right)-{\displaystyle \frac{\alpha {\mu }_0}{k}}O\left(t\right)\hfill \\ \hfill \le & \beta {\displaystyle \frac{\mu }{\mu +\gamma }}O\left(t\right)e^{-\mu {\tau }_1}-{\displaystyle \frac{\alpha {\mu }_0}{k}}O\left(t\right)\hfill \\ \hfill =& {\displaystyle \frac{\alpha {\mu }_0}{k}}O\left(t\right)\left({\displaystyle \frac{\mu \beta k{e}^{-\mu {\tau }_1}}{\alpha {\mu }_0(\mu +\gamma )}}-1\right)\hfill \\ \hfill =& {\displaystyle \frac{\alpha {\mu }_0}{k}}O\left(t\right)({\left({\mathcal{R}}_0^d\right)}^2-1),\hfill \end{array}$$

since, from first equation of system (2), one can deduce that $S\left(t\right)\le \frac{\mu }{\mu +\gamma }$. Thus, one gets

$${\mathcal{L}}^{\prime }\le \frac{\alpha {\mu }_0}{k}O\left(t\right)({\left({\mathcal{R}}_0^d\right)}^2-1).$$

Therefore, ${\mathcal{L}}^{\prime }\le 0$ if ${\mathcal{R}}_0^d\le 1$. Using LaSalle’s invariance principle, one finds that $F^{*}$ is globally asymptotically stable (GAS) in $\mathcal{R}$. □

3.5. Stability Analysis of the Toxoplasmosis-Endemic Steady State

The characteristic equation for the toxoplasmosis-endemic steady state is given by

$$|J_E+J_{{\tau }_1}{e}^{-\lambda {\tau }_1}+J_{{\tau }_2}{e}^{-\lambda {\tau }_2}-\lambda I|=0,$$

(14)

where

$$J_E=\left[\begin{array}{cccc}-\beta O^{*}-\gamma -\mu & 0& -\mu & -\beta S^{*}\\ \beta O^{*}& -\mu & 0& \beta S^{*}\\ 0& 0& -\alpha & 0\\ 0& 0& 0& -{\mu }_0\end{array}\right],J_{{\tau }_1}=\left[\begin{array}{cccc}0& 0& 0& 0\\ -\beta O^{*}{\mathrm{e}}^{-\mu {\tau }_1}& 0& 0& -\beta S^{*}{\mathrm{e}}^{-\mu {\tau }_1}\\ \beta O^{*}{\mathrm{e}}^{-\mu {\tau }_1}& 0& 0& \beta S^{*}{\mathrm{e}}^{-\mu {\tau }_1}\\ 0& 0& 0& 0\end{array}\right],\mathrm{and}{J}_{{\tau }_2}=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& k& 0\end{array}\right].$$

The explicit form of Equation (14) is

$$\begin{array}{cc}\hfill & -\left(\mu +\lambda \right)(-S^{*}{\mathrm{e}}^{-\lambda {\tau }_1-\lambda {\tau }_2-\mu {\tau }_1}\left(\gamma \beta k+\beta k\lambda +\beta k\mu \right)+O^{*}{\mathrm{e}}^{-{\tau }_1\left(\mu +\lambda \right)}\beta \lambda \mu \hfill \\ & +O^{*}{\mathrm{e}}^{-{\tau }_1\left(\mu +\lambda \right)}\beta \mu {\mu }_0+O^{*}\alpha \beta \lambda +O^{*}\alpha \beta {\mu }_0+O^{*}\beta {\lambda }^2+O^{*}\beta \lambda {\mu }_0+\gamma \alpha \lambda +\alpha \gamma {\mu }_0\hfill \\ & +{\lambda }^2\alpha +\alpha \lambda \mu +\alpha \lambda {\mu }_0+\alpha {\mu }_0\mu +\gamma {\lambda }^2+\gamma \lambda {\mu }_0+{\lambda }^3+{\lambda }^2\mu +{\lambda }^2{\mu }_0+\lambda \mu {\mu }_0)=0.\hfill \end{array}$$

Clearly, one eigenvalue is $\lambda =-\mu$. The other eigenvalues are the roots of the following equation

$$k{S}^{*}\beta \left(\mu +\gamma +\lambda \right){\mathrm{e}}^{\left(-{\tau }_1-{\tau }_2\right)\lambda -\mu {\tau }_1}-\left({\mu }_0+\lambda \right)\left({\mathrm{e}}^{-{\tau }_1\left(\mu +\lambda \right)}\beta \mu O^{*}+\left(\alpha +\lambda \right)\left(\beta O^{*}+\gamma +\lambda +\mu \right)\right)=0.$$

(15)

Let us assume that a root $\lambda$ of this equation has a non-negative real part. Thus, Equation (15) can be rewritten as

$$-\left({\mu }_0+\lambda \right){\mathrm{e}}^{-{\tau }_1\left(\mu +\lambda \right)}\beta \mu O^{*}=-k{S}^{*}\beta \left(\mu +\gamma +\lambda \right){\mathrm{e}}^{\left(-{\tau }_1-{\tau }_2\right)\lambda -\mu {\tau }_1}+\left({\mu }_0+\lambda \right)\left(\alpha +\lambda \right)\left(\beta O^{*}+\gamma +\lambda +\mu \right).$$

The modulus of the left-hand side of this last equation satisfies

$$\begin{array}{cc}\hfill ∥-\left({\mu }_0+\lambda \right){\mathrm{e}}^{-{\tau }_1\left(\mu +\lambda \right)}\beta \mu O^{*}∥& =∥{\mu }_0+\lambda ∥e^{-\mu {\tau }_1}\beta \mu O^{*}∥e^{-{\tau }_1\lambda }∥\hfill \\ \hfill & =∥{\mu }_0+\lambda ∥e^{-\mu {\tau }_1}\beta \mu ∥e^{-{\tau }_1\lambda }∥\frac{\alpha (\gamma +\mu )\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)}{\beta \left(\mu {\mathrm{e}}^{-\mu {\tau }_1}+\alpha \right)}\hfill \\ \hfill & <∥{\mu }_0+\lambda ∥\alpha (\gamma +\mu )\left({\left({\mathcal{R}}_0^d\right)}^2-1\right).\hfill \end{array}$$

Now, the modulus of the right-hand side of Equation (15) satisfies

$$\begin{array}{cc}\hfill & ∥-k{S}^{*}\beta \left(\mu +\gamma +\lambda \right){\mathrm{e}}^{\left(-{\tau }_1-{\tau }_2\right)\lambda -\mu {\tau }_1}+\left({\mu }_0+\lambda \right)\left(\alpha +\lambda \right)\left(\beta O^{*}+\gamma +\lambda +\mu \right)∥\hfill \\ \hfill & >\left|∥\left({\mu }_0+\lambda \right)\left(\alpha +\lambda \right)\left(\beta O^{*}+\gamma +\lambda +\mu \right)∥-∥k{S}^{*}\beta \left(\mu +\gamma +\lambda \right){\mathrm{e}}^{\left(-{\tau }_1-{\tau }_2\right)\lambda -\mu {\tau }_1}∥\right|\hfill \\ \hfill & >∥\gamma +\lambda +\mu ∥\left|{\mu }_0\alpha -k{S}^{*}\beta e^{-\mu {\tau }_1}∥e^{\left(-{\tau }_1-{\tau }_2\right)\lambda }∥\right|\hfill \\ \hfill & >\left(\gamma +\mu \right)\left|{\mu }_0\alpha -k\frac{\alpha {\mu }_0}{\beta k{\mathrm{e}}^{-\mu {\tau }_1}}\beta e^{-\mu {\tau }_1}∥e^{\left(-{\tau }_1-{\tau }_2\right)\lambda }∥\right|=\left(\gamma +\mu \right){\mu }_0\alpha \left|1-∥e^{\left(-{\tau }_1-{\tau }_2\right)\lambda }∥\right|.\hfill \end{array}$$

Accordingly, one gets that

$$\begin{array}{c}\hfill ∥{\mu }_0+\lambda ∥\alpha (\gamma +\mu )\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)>\left(\gamma +\mu \right){\mu }_0\alpha \left|1-∥e^{\left(-{\tau }_1-{\tau }_2\right)\lambda }∥\right|\end{array}$$

i.e.,

$$\begin{array}{c}\hfill ∥{\mu }_0+\lambda ∥\left({\left({\mathcal{R}}_0^d\right)}^2-1\right)>{\mu }_0\left|1-∥e^{\left(-{\tau }_1-{\tau }_2\right)\lambda }∥\right|,\forall {\mathcal{R}}_0^d>1.\end{array}$$

Thus, as ${\mathcal{R}}_0^d\to 1^{+}$, we have a contradiction. Therefore, $\lambda$ cannot have a positive real part.

Now, Equation (15) can be rewritten as

$$\begin{array}{l}({\mu }_0+\lambda )(\alpha +\lambda )(\beta O^{*}+\gamma +\mu +\lambda )+\beta \mu O^{*}({\mu }_0+\lambda ){\mathrm{e}}^{-\mu {\tau }_1-\lambda {\tau }_1}\\ -\beta k{S}^{*}(\gamma +\mu +\lambda ){\mathrm{e}}^{-\mu {\tau }_1-\lambda {\tau }_1-\lambda {\tau }_2}=0.\end{array}$$

(16)

Equation (16) is a function of $\lambda ,{\tau }_1$, and ${\tau }_2$ with the following form

$$\begin{array}{c}\hfill P(\lambda ,{\tau }_1,{\tau }_2)=\left(P_0\left(\lambda \right)+P_1\left(\lambda \right)e^{-\lambda {\tau }_1}\right)+\left(P_2\left(\lambda \right)+P_3\left(\lambda \right)e^{-\lambda {\tau }_1}\right)e^{-\lambda {\tau }_2},\end{array}$$

(17)

where $P_0\left(\lambda \right)=({\mu }_0+\lambda )(\alpha +\lambda )(\beta O^{*}+\mu +\gamma +\lambda )$. It is clear that the degree of $P_0\left(\lambda \right)$ is three. The other polynomials $P_i\left(\lambda \right)$ are given by $P_1\left(\lambda \right)=\beta \mu ({\mu }_0+\lambda )O^{*}{\mathrm{e}}^{-\mu {\tau }_1}$, $P_2\left(\lambda \right)=0$, and $P_3\left(\lambda \right)=-\alpha {\mu }_0\left(\gamma +\lambda +\mu \right)$. Notice that

$$deg\left(P_0\left(\lambda \right)\right)\ge max\{deg\left(P_1\left(\lambda \right)\right),deg\left(P_2\left(\lambda \right)\right),deg\left(P_3\left(\lambda \right)\right)\},$$

and

$$\begin{array}{cc}\hfill \sum _{i=0}^3{P}_i\left(0\right)& ={\mu }_0\alpha (\beta O^{*}+\mu +\gamma )+\beta \mu {\mu }_0{O}^{*}{\mathrm{e}}^{-\mu {\tau }_1}-\alpha {\mu }_0\left(\gamma +\mu \right).\hfill \end{array}$$

The parameters of model (2) satisfy $\alpha >0,{\mu }_0>0,\mu >0,$ and $\beta >0$; then, the above expression will be zero if $O^{*}=0$. Note that if $O^{*}$ is a component of the endemic equilibrium point $E{P}^{*}$, then $O^{*}$ cannot be equal to zero. Hence, the above expression cannot be zero. Thus, one gets that

$$P_0\left(0\right)+P_1\left(0\right)+P_2\left(0\right)+P_3\left(0\right)\ne 0.$$

Therefore, $\lambda =0$ cannot be an eigenvalue or a root of Equation (14). Note that the polynomials $P_0,P_1,P_2,$ and $P_3$ do not have common zeros. On the other hand, the polynomials $P_0\left(\lambda \right),P_1\left(\lambda \right),P_2\left(\lambda \right),$ and $P_3\left(\lambda \right)$ hold with

$$\begin{array}{c}\hfill \underset{\lambda \to \infty }{lim}\left(\left|\frac{P_1\left(\lambda \right)}{P_0\left(\lambda \right)}\right|+\left|\frac{P_2\left(\lambda \right)}{P_0\left(\lambda \right)}\right|+\left|\frac{P_3\left(\lambda \right)}{P_0\left(\lambda \right)}\right|\right)<1.\end{array}$$

Therefore, Equation (17) is a characteristic equation (see [93]).

To study stability switching, we consider purely imaginary characteristic roots since $\lambda \ne 0$ and roots of a real function always occur as conjugate pairs. Let us assume that $\lambda =i\omega (\omega >0)$ is an eigenvalue. Substituting this into Equation (17), one obtains

$$\begin{array}{c}\hfill P(i\omega ,{\tau }_1,{\tau }_2)=\left(P_0\left(i\omega \right)+P_1\left(i\omega \right)e^{-i\omega {\tau }_1}\right)+\left(P_2\left(i\omega \right)+P_3\left(i\omega \right)e^{-i\omega {\tau }_1}\right)e^{-i\omega {\tau }_2}=0.\end{array}$$

(18)

Since $P_2\left(i\omega \right)=0$ and $∥{\mathrm{e}}^{-i\omega {\tau }_2}∥=1,$ it follows that

$$\left|\right|P_0\left(i\omega \right)+P_1\left(i\omega \right){\mathrm{e}}^{-i\omega {\tau }_1}\left|\right|=\left|\right|P_3\left(i\omega \right){\mathrm{e}}^{-i\omega {\tau }_1}\left|\right|=\left|\right|P_3\left(i\omega \right)\left|\right|.$$

(19)

Then, the above equation can be written as

$$(P_0\left(i\omega \right)+P_1\left(i\omega \right){\mathrm{e}}^{-i\omega {\tau }_1})(\overline{P_0}\left(i\omega \right)+\overline{P_1}\left(i\omega \right){\mathrm{e}}^{i\omega {\tau }_1})=\left(P_3\left(i\omega \right){\mathrm{e}}^{-i\omega {\tau }_1}\right)\left(\overline{P_3}\left(i\omega \right){\mathrm{e}}^{i\omega {\tau }_1}\right).$$

After simplifying, one gets

$$\begin{array}{cc}\hfill {∥P_0\left(i\omega \right)∥}^2+{∥P_1\left(i\omega \right)∥}^2-{∥P_3\left(i\omega \right)∥}^2=& -2Re\left(P_0\left(i\omega \right)\overline{P_1}\left(i\omega \right)\right)cos\left(\omega {\tau }_1\right)\hfill \\ \hfill & +2Im\left(P_0\left(i\omega \right)\overline{P_1}\left(i\omega \right)\right)sin\left(\omega {\tau }_1\right).\hfill \end{array}$$

Hence,

$${∥P_0\left(i\omega \right)∥}^2+{∥P_1\left(i\omega \right)∥}^2-{∥P_3\left(i\omega \right)∥}^2=2A_1\left(i\omega \right)cos\left(\omega {\tau }_1\right)-2B_1\left(i\omega \right)sin\left(\omega {\tau }_1\right),$$

(20)

where

$$A_1\left(i\omega \right)=-Re\left(P_0\left(i\omega \right)\overline{P_1}\left(i\omega \right)\right)\mathrm{and}{B}_1\left(\omega \right)=-Im\left(P_0\left(i\omega \right)\overline{P_1}\left(i\omega \right)\right).$$

(21)

Suppose that there exists $\omega$ such that $A_1{\left(i\omega \right)}^2+B_1{\left(i\omega \right)}^2=0.$ Then

$$\begin{array}{l}{A}_1\left(i\omega \right)=B_1\left(i\omega \right)=0⟺P_0\left(i\omega \right)\overline{P_1}\left(i\omega \right)=0,\\ {∥P_0\left(i\omega \right)∥}^2+{∥P_1\left(i\omega \right)∥}^2={∥P_3\left(i\omega \right)∥}^2.\end{array}$$

(22)

On the other hand, we have that

$$\begin{array}{cc}\hfill P_0\left(i\omega \right)=& (-{\omega }^3+\alpha {\mu }_0\omega +(\alpha \omega +{\mu }_0\omega )(\beta O^{*}+\gamma +\mu ))i-(\alpha \omega +{\mu }_0\omega )\omega \hfill \\ \hfill & -({\omega }^2-{\mu }_0\alpha )(\beta O^{*}+\gamma +\mu ),\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill P_1\left(i\omega \right)& ={\mathrm{e}}^{-\mu {\tau }_1}\beta \mu O^{*}({\mu }_0+i\omega ),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \overline{P_1}\left(i\omega \right)& ={\mathrm{e}}^{-\mu {\tau }_1}\beta \mu O^{*}({\mu }_0-i\omega ).\hfill \end{array}$$

(24)

After simplification, we can write the following equation

$$\begin{array}{cc}\hfill P_0\overline{P_1}& =\beta \mu O^{*}{\mathrm{e}}^{-\mu {\tau }_1}({\mu }_0^2+{\omega }^2)(\beta O^{*}+\alpha +\gamma +\mu )\hfill \\ & +i\beta \mu \omega O^{*}{\mathrm{e}}^{-\mu {\tau }_1}({\mu }_0^2+{\omega }^2)(\alpha \beta O^{*}+\gamma \alpha +\alpha \mu -{\omega }^2)=0.\hfill \end{array}$$

Now, both the real part and the imaginary part of the above equation must be zero separately

$$\beta \mu O^{*}{\mathrm{e}}^{-\mu {\tau }_1}({\mu }_0^2+{\omega }^2)(\beta O^{*}+\alpha +\gamma +\mu )=0,$$

(25)

and

$$\beta \mu \omega O^{*}{\mathrm{e}}^{-\mu {\tau }_1}({\mu }_0^2+{\omega }^2)(\alpha \beta O^{*}+\gamma \alpha +\alpha \mu -{\omega }^2)=0.$$

(26)

Solving Equation (25) for $\omega$, one gets that $\omega =\pm {\mu }_{0}i$. Now, solving Equation (26), one obtains that $\omega =\pm {\mu }_{0}i$ or $\omega =\pm \sqrt{O^{*}\alpha \beta +\gamma \alpha +\alpha \mu }$. Therefore, the only option with which both Equations (25) and (26) are satisfied is to have $\omega =\pm {\mu }_{0}i$. Since $\omega$ is a positive real number, $P_0\overline{P_1}$ cannot be zero for any value of $\omega$ for all ${\tau }_1\in {\mathbb{R}}_{+}$.

Next, we should consider the case where $A_1{\left(\omega \right)}^2+B_1{\left(\omega \right)}^2>0$. In [93], it is shown that, if the inequality

$$\left||P_0{|}^2+|P_1{|}^2-|P_2{|}^2-{\left|P_3\right|}^2\right|\le 2\sqrt{A_1^2+B_1^2}$$

(27)

is not satisfied, then we can guarantee that there is no real positive $\omega$. Therefore, we would be able to conclude that the toxoplasmosis-endemic equilibrium point $E{P}^{*}$ cannot lose stability regardless of the values of the time delays ${\tau }_1$ and ${\tau }_2$. We were not able to prove that this inequality is not true for any positive parameter values and discrete time delays ${\tau }_1$ and ${\tau }_2$. However, we tested the inequality with a large number of values for the parameters and time delays and it was found that the inequality was not satisfied. Thus, these numerical results suggest that $E{P}^{*}$ is locally stable and no Hopf bifurcation occurs [105].

In the next section, we provide numerical simulations that give support to all the above mathematical theoretical results. Moreover, we describe numerical simulations that show the importance of the time discrete delay and the basic reproduction number ${\mathcal{R}}_0^d$ in the delayed model. In addition, we discuss simulations that suggest that the toxoplasmosis-endemic equilibrium point $E{P}^{*}$ cannot lose stability regardless of the values of the time delays ${\tau }_1$ and ${\tau }_2$.

4. Numerical Simulations

In order to investigate the dynamics of toxoplasmosis and corroborate the theoretical findings from the preceding sections, we undertook a range of in silico simulations of mathematical model (1). For the purpose of acquiring distinct values for ${\mathcal{R}}_0^d$ and ${\mathcal{R}}_0$, we investigated a range of scenarios with a variety of values for the time delays and the parameters. In particular, we were interested in the values ${\mathcal{R}}_0^d>1$ and ${\mathcal{R}}_0^d<1$ since these are the threshold values for the long-term behavior of the toxoplasmosis disease. The in silico simulations provide further insights, such as into the transient and long-term behavior of the toxoplasmosis disease in the cat population and the oocyst population. We varied the values of the transmission rate and infectivity of the oocysts. Thus, we could study different scenarios related to the dynamics.

We relied on the Matlab built-in function dde23 to compute the numerical solutions to the delay differential equations [77,114]. For the in silico simulations, we used the values presented in

Table 1. Explanation of the parameters of the delayed model (1) and their values.

Parameter	Description	Value
$\mu$	Birth/death rates (cats)	$1/260$ (1/weeks) [115]
$\alpha$	Shedding period	$1/2$ (1/weeks) [56]
${\mu }_0$	Clearance rate	$1/26$ (1/day) [54,56]
k	Oocysts per day (cat)	$20\times {10}^6$ (1/day) [116]
$\beta$	Transmission rate	Varied
$\gamma$	Vaccination rate	Varied

. The initial conditions were varied in order to obtain a more robust support for the theoretical results. In addition, the equilibrium points were computed for each scenario using the numerical simulations and we contrasted them with the theoretical ones.

4.1. Numerical Simulations for the Scenarios When ${\mathcal{R}}_0^d<1$

Here, we describe the in silico simulations where the basic reproduction number ${\mathcal{R}}_0^d<1$ and both time delays (${\tau }_1$ and ${\tau }_2$) were different from zero. We used the following initial condition: $S\left(0\right)=0.5,E\left(0\right)=0.3,I\left(0\right)=0.2,V_R\left(0\right)=0$, and $O\left(0\right)=1$. Figure 2 shows that the system reached the toxoplasmosis-free equilibrium $F^{*}$. In this case, we have that ${\mathcal{R}}_0>1$, but as can be seen, the crucial parameter is ${\mathcal{R}}_0^d$. It is important to mention that when there are no time delays one obtains that ${\mathcal{R}}_0^d={\mathcal{R}}_0$. Thus, it is possible that, for the toxoplasmosis model without time delays, the disease can be prevalent, but when time delays are taken into account, the disease is eradicated. This has a relevant impact from a public health point of view, since it shows the relevance of the time delays in the dynamics of toxoplasmosis. Figure 3 shows the same scenario but using much larger time delays (${\tau }_1=10$ and ${\tau }_2=10$). Again, the system reached the toxoplasmosis-free equilibrium $F^{*}$ since ${\mathcal{R}}_0^d<1$. Notice that the increase in the time delays did not affect the value of $F^{*}$. Finally, we increased the time delays to much larger values to show that long-term behavior would be maintained despite ${\mathcal{R}}_0^d$ being significantly lower. This can be observed in Figure 4.

4.2. Numerical Simulations for the Scenarios When ${\mathcal{R}}_0^d>1$

Here, we describe the simulations performed where the basic reproduction number ${\mathcal{R}}_0^d>1$ and both time delays (${\tau }_1$ and ${\tau }_2$) were different from zero. Figure 5 shows that the system reached the toxoplasmosis-endemic equilibrium $E{P}^{*}$. We increased the transmission rate $\beta$ in order to have ${\mathcal{R}}_0^d>1$. The system must approach the endemic steady state $E{P}^{*}$ in accordance with the theoretical results. This section’s second numerical simulation took a greater transmission rate $\beta$ into account. Again, as Figure 6 shows, the system reached the toxoplasmosis-endemic equilibrium $E{P}^{*}$ but with larger populations of infected cats and oocysts. This was expected since the toxoplasmosis-endemic equilibrium $E{P}^{*}$ depends on the time delay $\beta$ (see Equation (6)). Moreover, the theoretical toxoplasmosis-endemic equilibrium $E{P}^{*}$ agreed with the numerical endemic equilibrium. The next simulation considered very large time delays to study the possibility of the system losing stability. The results from the previous section suggest that this is not feasible from a mathematical theoretical viewpoint. Figure 7 shows that the system still reached the toxoplasmosis-endemic equilibrium $E{P}^{*}$ but with larger populations of infected cats and oocysts due to the large transmission rate. It is important to mention that, for many models based on delay differential equations, large time delays can give conditions such that periodic solutions arise. This is an important aspect to consider since, instead of having a steady state, oscillatory behavior could be obtained. It is important to understand this situation from an ecological viewpoint since incorrect public health interventions can be implemented. However, the numerical simulations from this section and the results from the previous section suggest that the mathematical model from Equation (1) cannot have periodic solutions, even for very large time delays. We tested further time delays, such as ${\tau }_1>1\times {10}^3$ and ${\tau }_2>1\times {10}^3$. Although those numerical results are not shown here, they were in good agreement with the previous results.

4.3. Numerical Tests for Hopf Bifurcation When ${\mathcal{R}}_0^d>1$

Here, we show that, for a set of parameter values that are similar to the previous numerical simulations, the inequality (27) is not satisfied. Therefore, a Hopf bifurcation cannot occur for the set of parameter values that we chose. Let us first rewrite Equation (27) as

$$\left||P_0{|}^2+|P_1{|}^2-|P_2{|}^2-{|P_3|}^2\right|-2\sqrt{A_1^2+B_1^2}\le 0.$$

(28)

If we vary $\omega$, we can compute the values of the left-hand side of the inequality (28). Figure 8 shows these values for a large range of values for $\omega$. We use the log scale for the values of the left-hand side of the inequality (28) and, as can be seen, the inequality (28) is not satisfied. It is important to point out that this numerical procedure is not a proof since there are infinitely many potential values for the parameters of model (1). However, these results and the previous numerical simulations suggest that a Hopf bifurcation cannot occur.

5. Conclusions

In order to explore the dynamics of toxoplasmosis, we designed a more realistic mathematical model in this work. The model considers two discrete time delays: the first time delay ${\tau }_1$ is related to the time it takes the susceptible cats to become infected after effective contact with an oocyst and the second time delay ${\tau }_2$ is the time that elapses between when the oocysts become present in the environment and when they are able to infect. The main achievement was finding the conditions under which the toxoplasmosis can disappear from the cat population. This can only occur if there are no oocysts in the environment. We also found conditions where the toxoplasmosis became endemic. In addition, we investigated the impact of the two time delays on the toxoplasmosis dynamics. Using dynamical systems theory, we were able to find the basic reproduction number ${\mathcal{R}}_0^d$ for the delay model, which differed from that for the model without delays. We demonstrated that it is a threshold that determines the global long-term dynamics of the toxoplasmosis. We proved that, if ${\mathcal{R}}_0^d<1$, the toxoplasmosis would be eradicated and the toxoplasmosis-free equilibrium would be globally stable. On the other hand, we proved that, if the threshold parameter ${\mathcal{R}}_0^d$ is greater than one, then there is only one toxoplasmosis-endemic equilibrium point. We obtained partial theoretical results and performed numerical simulations that suggested that, if ${\mathcal{R}}_0^d>1$, then the toxoplasmosis-endemic equilibrium point is globally stable. This is a crucial result since it means that, if the basic reproduction number ${\mathcal{R}}_0^d$ is greater than one, the system goes exactly to the toxoplasmosis-endemic steady state. Thus, a globally attractive toxoplasmosis-endemic equilibrium exists regardless of the time delay. Additional numerical simulations were performed and supported the theoretical stability results.

The results obtained in this paper allow us to discuss further the effect of vaccination programs on toxoplasmosis under a more realistic scenario. We provided mathematical analysis that is valid with any parameter values. It is important to mention that, in reality, the values of the parameters have uncertainty, as is common in biology, epidemiology, and public health. The proposed model can be calibrated to real data if we have enough data for some state variables, such as the number of infected cats and the density of oocysts in the environment. This calibration is possible since some values of the parameters are known approximately. However, the uniqueness of the calibration would require identifiability analysis [117,118]. Another aspect that we would like to point out is that the temperature impacts the time delays related to the infectiousness of the oocysts [57]. Therefore, due to climate change, the time delays may be modified in the future. In summary, this study provides further insight into toxoplasmosis dynamics. However, there are still a great number of open problems where even more realistic facts can be included in the mathematical models. For instance, including intermediate hosts, such as mice, would be an important aspect. This is one of the limitations of this work. As with any mathematical model related to biological processes, there are some limitations related to the complex real world. The developed model does not consider humans, who, in a very complex process, might affect the toxoplasmosis dynamics. The model also includes homogeneous mixing, which might not be suitable for some scenarios. Finally, future works that include additional intermediate hosts or new time delays, although challenging, could be useful.