Mathematical Modeling of Toxoplasmosis Considering a Time Delay in the Infectivity of Oocysts

Abstract

In this paper, we study the effect of the introduction of a time delay on the dynamics of toxoplasmosis. This time delay is the elapsed time from when oocysts become present in the environment and when they become infectious. We construct a mathematical model that includes cats and oocysts in the environment. We include the effect of oocysts, since they are crucial for the dynamics of toxoplasmosis. The likelihood of the acquisition of Toxoplasma gondii infection depends on the environmental load of the parasite. Furthermore, the model considers the possibility of vaccination of the feline host. In the mathematical model, we consider directly the infection of cats through the oocysts shed by other cats. We prove that the basic reproduction number ${\mathcal{R}}_0$ is a secondary parameter that determines the global dynamics of toxoplasmosis in cat populations. We study the effect of the time delay on the stability of the steady states. We find that the time delay cannot change the stability of the endemic state, which is an important result from the biological point of view. Numerical simulations are performed to support the theoretical results and obtain further insight into understanding toxoplasmosis dynamics in cat populations.

1. Introduction

T. gondii is an apicomplexan parasite with a worldwide distribution. Cats are the final host, and humans (and other warm-blooded vertebrates) are intermediate hosts [1,2]. Following a period of asexual reproduction by tachyzoite forms, the parasite enters a latent phase in the bradyzoite stage that persists for the host’s lifetime in pseudocysts, macrophages and neurons of various tissues, notably in the brain [1,2]. T. gondii encysts in the brain, where it can cause inflammation and inhibit apoptosis [2]. T. gondii is present in a variety of animals worldwide and mostly in cats [3,4,5]. The T. gondii parasite in cats goes through all stages of its life cycle and does not affect the cat’s life. Cats are immune to toxoplasma, and it is estimated that approximately 20 million oocysts, between 4 and 13 days after the infection, can be cast [6]. T. gondii can be vertically transmitted [7]. Cats are one of the Western world’s most popular pets [8].

T. gondii infection is highly prevalent around the world. The seroprevalence in cats and humans over different regions is often around 30 to 40 percent [9]. About one-third of the world human population has antibodies to T. gondii [2]. In 2013, the burden of congenital toxoplasmosis was estimated to be 1.2 million [10]. Eye issues from congenital infection are difficult to be identified at birth but are present in 20%–80% of congenitally infected persons by adulthood [11]. In addition, toxoplasmosis is considered to be a leading cause of death attributed to food-borne illness in the United States [11]. The testing of all pregnant women for T. gondii infection is routine in some European countries, but the cost–benefit of such mass screening has been debated [12].

Cats are the main reason for the maintenance of T. gondii in the environment [9]. Domestic cats play a crucial part in the spread of T. gondii by excreting environmentally resistant oocysts that may infect humans and other warm-blooded animals [13]. The oocysts released in the feces of infected cats can contaminate the environment, including water, vegetables and other kinds of food [14,15]. Thus, the T. gondii can be ingested orally. Gamma irradiation can help to decontaminate food contaminated with oocysts [15]. The contamination of drinking water with T. gondii oocysts has been documented [14]. Previous analysis shown that the feces of cats shedding T. gondii may contain $2.5\times {10}^6$ oocysts/gr and a cat may shed as many as 20 million oocysts per day [16]. Oocysts of T. gondii can be easily dispersed in the environment. The challenge of controlling the population of cats has been mentioned in several studies [17]. Studies have found population densities of approximately 500–1300 cats/square mile in urban areas [17,18]. Thus, controlling toxoplamosis is related to cat populations, and its study is justified.

Much scientific literature has documented the public health importance of parasites since they can impact population dynamics in different ways [19]. Mathematical models are simplifications of real-world processes that allow us to obtain better insight into the dynamics of these processes. Mathematical models have been used extensively to investigate the dynamics of diseases in populations [20,21,22]. Previous mathematical models have been used to study the dynamics of toxoplasmosis in human populations [13,23,24,25,26]. In [13], the authors studied the risks of exposure to T. gondii oocysts for humans and livestock living on farms. In addition, the authors discussed the role of young cats in the maintenance of environmental contamination by T. gondii oocysts on farms. The study included probabilistic factors in order to deal with some uncertainties related to the T. gondii infection dynamic in cat populations. An excellent systematic review of existing mathematical models for the transmission of T. gondii is presented in [10].

Mathematical models that consider time delays have been proposed and studied previously [27,28,29,30]. These models are useful to gain insight on epidemics and a variety of viruses [31,32,33,34,35,36]. In this paper, we extend the mathematical model presented in [37]. We include in the model a time delay since when oocysts are shed by cats, they are not immediately infective. In fact, the Toxoplasma parasite does not become infectious until one to five days after it is shed in a cat’s feces [11]. This is due to the fact that oocysts released from cat feces are not infective until sporulation, a process that requires external environmental exposure of 24–48 h [38]. Moreover, sporulation is affected by the temperature and, therefore, introduces uncertainty in the time delay [39]. Mathematical models that consider differential equations with delays are introduced in order to account for incubation periods, maturation times, age structure, etc. [40]. However, the resulting models present difficulties, which often require rather complicated tools for their analysis [40].

The mathematical model also considers a vaccination program. The model divides the cat population into three subpopulations: susceptibles $S\left(t\right)$, infected $I\left(t\right)$ and vaccinated/recovered $V_R\left(t\right)$. The model also considers the subpopulation of T. gondii oocysts since this form can transmit toxoplasmosis. The model is based on a system of differential equations, which includes different parameters that can impact dynamics of the transmission of toxoplasmosis in cat populations. The inclusion of a vaccination program allows the assessment of control strategies of health institutions. Vaccinated cats have lifelong immunity [41,42]. It is important to remark that in some regions where there are no cats, the prevalence of toxoplasmosis is null [4]. Previous articles studied and incorporated the vaccination of cats [37,43,44,45,46,47]. The inclusion of vaccines in the model allows the study of the impact of control strategies. In [46], the authors proposed a model to study when a vaccination program is cost effective.

We study the stability of the equilibrium points of the system and compute the basic reproductive number ${\mathcal{R}}_0$, which is a crucial threshold parameter for the dynamics of many diseases at the population level. We investigate the impact that the time delay has on the dynamics of toxoplasmosis and how modifying it can affect toxoplasmosis prevalence. Numerical simulations of the mathematical model are included to support the theoretical results and obtain insight into the spread of the disease in cat populations. The simulations are useful to assess prevention and control strategies. For instance, in [44] a computer simulation model that includes the vaccination of cats is presented.

The organization of this paper is the following: In Section 2, we present the mathematical model. In Section 3, we study the stability of the equilibrium points and calculate the basic reproduction number ${\mathcal{R}}_0$. In Section 4, the numerical simulations of the mathematical model are presented. Section 5 contains the conclusions.

2. Mathematical Model

In this section, the mathematical model for the transmission of toxoplasmosis in cat populations is introduced. The model includes a time delay that represents the elapsed time from when the oocysts become present in the environment and when they are able to infect. The mathematical model also includes the possibility of vaccination of the cats. The constructed mathematical model considers the infection of cats through the oocysts shed by cats. The cats are the main host of the T. gondii parasite, and they shed them in the environment [9]. The model constructed includes direct contact between cats and oocysts. Contamination of the environment by oocysts is common [48]. As expected, the likelihood of infection depends on the environmental load of T. gondii [44]. Thus, the rate of infection is modeled using the amount of oocysts in the environment, which depends on the population of infected cats [12].

The mathematical model is based on a system of nonlinear ordinary differential equations. The model assumes that the oocysts suffer from natural decay when there are no infected cats. This is a common assumption in modeling biological processes. Other options, where there is a carrying capacity of the environment or a fixed lifetime, can be considered. Sporulated oocysts can survive for a long time under a variety of environmental conditions. In fact, oocysts can survive in moist soils for extended periods of time [12].

The model is based on the following assumptions:

•

The population of cats $N\left(t\right)$ is divided into three subpopulations: susceptible cats ($S\left(t\right)$), infected cats ($I\left(t\right)$), and vaccinated/recovered cats ($V_R\left(t\right)$). The cats in this last subpopulation have lifelong immunity.

•

The variable $O\left(t\right)$ denotes the number of oocysts in the environment.

•

The population of cats is constant, since the birth and death rates are equal to $\mu$.

•

Susceptible cats move to the infected subpopulation $I\left(t\right)$ if they have effective contact with a oocyst (at rate $\beta$).

•

The period from when the oocysts are shed by the cats until they are infective is a time delay of $\tau$.

•

Susceptible cats are vaccinated at a rate $\gamma$. An infected cat transits to the vaccinated/recovered subpopulation $V_R\left(t\right)$ at a rate of $\alpha$.

•

The increase in oocysts is proportional to the number of infectious cats.

•

${\mu }_0$ is the removal rate of oocysts from the environment.

•

Vertical transmission occurs in the cats. [12,49,50,51,52,53].

•

Vaccinated cats do not shed oocysts and have lifetime immunity.

•

Homogeneous mixing.

Thus, the model is a first order nonlinear system of ordinary differential equations given by

$$\begin{array}{ccc}\dot{S}\left(t\right)& \hfill =& \mu V_R\left(t\right)-\beta S\left(t\right)O\left(t\right)-\gamma S\left(t\right),\hfill \\ \dot{I}\left(t\right)& \hfill =& \beta S\left(t\right)O\left(t\right)-\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\left(t\right)-{\mu }_{0}O\left(t\right),\hfill \end{array}$$

(1)

where $k>0$ is the amount of oocysts shed per infected cat.

Without loss of generality, the population of cats is assumed to be $N\left(t\right)=S\left(t\right)+I\left(t\right)+V_R\left(t\right)=1$. The mathematical diagram is depicted in Figure 1. We study the dynamics of the model (1) in the restricted region:

$$\Omega =\left\{(S,I,V_R,O)\in {\mathbb{R}}_{+}^4/0<S\le \frac{\mu }{\mu +\gamma },0\le I\le {\displaystyle \frac{\beta k\mu }{{\mu }_0\alpha (\mu +\gamma )}},S+I+V_R=1,0\le O\le \frac{k}{{\mu }_0}\right\},$$

where ${\mathbb{R}}_{+}^4$ denotes the nonnegative cone of ${\mathbb{R}}^4$.

Now, we introduce a time delay $\tau >0$ in the mathematical model (1) in order to take into account that the Toxoplasma parasite does not become infectious until one to five days after it is shed in a cat’s feces [11]. Thus the mathematical model with delay becomes

$$\begin{array}{ccc}\dot{S}\left(t\right)& \hfill =& \mu V_R\left(t\right)-\beta S\left(t\right)O\left(t\right)-\gamma S\left(t\right),\hfill \\ \dot{I}\left(t\right)& \hfill =& \beta S\left(t\right)O\left(t\right)-\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 )-{\mu }_{0}O\left(t\right),\hfill \end{array}$$

(2)

with $S\left(t\right)+I\left(t\right)+V_R\left(t\right)=1.$ Since the populations are normalized, we can study the following reduced model

$$\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{I}\left(t\right)& \hfill =& \beta S\left(t\right)O\left(t\right)-\alpha I\left(t\right),\hfill \\ \dot{O}\left(t\right)& \hfill =& kI(t-\tau )-{\mu }_{0}O\left(t\right),\hfill \end{array}$$

(3)

where $V_R\left(t\right)=1-S\left(t\right)-I\left(t\right)$. In addition, the system (3) satisfies the initial conditions given by

$$\begin{array}{c}\hfill S\left(0\right)>0,O\left(0\right)\ge 0,I\left(s\right)=\xi \left(s\right)\ge 0,s\in [-\tau ,0],\end{array}$$

(4)

with $\xi \left(s\right)$ continuous function defined from the interval $[-\tau ,0]$ to ${\mathbb{R}}_{+}$, and with the norm $\parallel \xi \parallel =\underset{-\tau \le s\le 0}{sup}\left|\xi \left(s\right)\right|.$ These initial conditions mean that, initially, there are susceptible cats and oocysts in the environment (or none). In addition, before $t=0$, the infected population of cats is nonnegative. We will see later that either the populations of oocysts or infected cats should be positive in order to have the possibility of an epidemic. Now, the model (3) can be written as

$$\begin{array}{c}\hfill \dot{\mathcal{W}}\left(t\right)=f(t,\mathcal{W}\left(t\right)),\end{array}$$

(5)

where

$$\begin{array}{c}\hfill \mathcal{W}\left(t\right)=\left(\begin{array}{c}S\left(t\right)\\ I\left(t\right)\\ O\left(t\right)\end{array}\right),f(t,\mathcal{W})=\left(\begin{array}{c}\mu (1-I(t\left)\right)-\beta S\left(t\right)O\left(t\right)-(\mu +\gamma )S\left(t\right),\beta S\left(t\right)O\left(t\right)-\alpha I\left(t\right)\\ \beta S\left(t\right)O\left(t\right)-\alpha I\left(t\right)\\ I(t-\tau )-{\mu }_{0}O\left(t\right),\end{array}\right)\end{array}$$

with the initial conditions given in (4). We use the following theorem:

Theorem 1.

Suppose Ω is an open set in $\mathcal{R}\times C\left([-\tau ,0],{\mathbb{R}}_{+}^n\right),$ $f:\Omega \to {\mathcal{R}}^n$ is continuous, and $f(t,\phi )$ is Lipschitizian in φ in each compact set in $\Omega .$ If $(\sigma ,\phi )\in \Omega$ then there is a unique solution of model (5) through $(\sigma ,\phi ).$

Proof. 

See [54] (p. 44). □

Let $\mathcal{S}=C\left([-\tau ,0],{\mathbb{R}}_{+}^3\right)$ be the Banach space of continuous functions mapping the interval $[-\tau ,0]$ into ${\mathbb{R}}_{+}^3$ which is equipped with the sup-norm. The function $f(t,\mathcal{W}(t\left)\right)$ given by in (5) is defined in $f:[0,\infty )\times C\left([-\tau ,0],{\mathbb{R}}_{+}^3\right)\to {\mathbb{R}}^3$ is continuous and satisfies a local Lipschitz condition. Indeed, for ${\phi }_1,{\phi }_2\in C_0\subset \mathcal{S},$ $C_0$ compact set, ${\phi }_1={\left(S_1\left(t\right),I_1\left(t\right),O_1\left(t\right)\right)}^t,$ ${\phi }_2={\left(S_2\left(t\right),I_2\left(t\right),O_2\left(t\right)\right)}^t$ and $0<a<b$ such that $t\in [a,b]$, we obtain that

$$\begin{array}{c}\hfill ∥f(t,{\phi }_1)-f(t,{\phi }_2)∥\le K∥{\phi }_1-{\phi }_2∥,\end{array}$$

where $K=\underset{t\in [a,b]}{max}\left\{{\mu }_0+2\beta \left|S_1\left(t\right)\right|,\mu +\gamma 2\beta \left|O_2\left(t\right)\right|,\mu +\alpha +k\right\}.$ Then, by Theorem 1 for any $\phi \in \mathcal{S}$ there exists a unique solution $\mathcal{W}(t,\phi )$ of the system (3) such that ${\mathcal{W}}_0=\phi .$ Next, from the biological point of view, one characterization of epidemiological models represented by differential equations must be that their solutions are positive for all time. The next section is devoted to the stability analysis of the mathematical model (3).

3. Stability Analysis of the Model

In this section, we perform the stability analysis of the mathematical model (3). First, we prove that the solutions are positive for all $t\ge 0$ and that they are bounded. These features are important from a realistic viewpoint.

Theorem 2.

If the parameters of model (3) are all positive and the initial conditions given by (4) are satisfied, then the solutions of the model (3) given by $\left(S\left(t\right),I\left(t\right),O\left(t\right)\right)$ remain positive and uniformly bounded in $[0,+\infty )$.

Proof. 

Using the first Equation of (3), we can see that

$$\begin{array}{c}\hfill \dot{S}\left(t\right)=\mu (1-S\left(t\right)-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}$$

Therefore, it follows that

$$\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,\end{array}$$

for all $t\ge 0.$ On the other hand, suppose that there exists a $t_1>0$ such that $I\left(t_1\right)=0,$ $\dot{I}\left(t_1\right)\le 0,$ and $I\left(t\right)>0$ for all $t\in (0,t_1).$ Thus, from the second equation of model (3), one obtains that

$$\begin{array}{c}\hfill 0\ge \dot{I}\left(t_1\right)=\beta S\left(t_1\right)O\left(t_1\right).\end{array}$$

(6)

However, $O\left(t\right)>0$ for all $t\in (0,t_1).$ If this is not the case, there exists a $t_2>0$ such that $t_2<t_1$, $\dot{O}\left(t_2\right)\le 0,$ $O\left(t_2\right)=0$ and $O\left(t\right)>0$ for $t\in (0,t_2).$ From the third equation of system (3), it follows that $0\ge kI(t_2-\tau ).$ Since $\tau >0$ and using the continuity of solution, then one gets that $0\ge \underset{\tau \to 0}{lim}I(t_2-\tau )=I\left(t_2\right).$ This is a contradiction because $I\left(t\right)>0$ for all $t\in (0,t_1).$ Thus, using the continuity of $O\left(t\right),$ it follows that $O\left(t\right)>0,$ with $t\in (0,t_1].$ Therefore, (6) is false and thus $I\left(t\right)\ge 0$ for $t\ge 0.$ From the third equation, it is concluded that $O\left(t\right)\ge 0$ for $t\ge 0.$

Now, the subpopulation $S\left(t\right)$ is bounded by $\frac{\mu }{\mu +\gamma }$, because by the standard comparison theorem [55], we can show that

$$S\left(t\right)\le S\left(0\right)exp(-(\mu +\gamma )t)+\frac{\mu }{\mu +\gamma }(1-exp(-(\mu +\gamma )t)).$$

In particular, if $S\left(0\right)\le \frac{\mu }{\mu +\gamma }$, then $S\left(t\right)\le \frac{\mu }{\mu +\gamma }$. Analogously, one gets that $O\left(t\right)\le \frac{k}{{\mu }_0}$ if $O\left(0\right)\le \frac{k}{{\mu }_0},$ and $I\left(t\right)\le {\displaystyle \frac{\beta \mu k}{\alpha {\mu }_0(\mu +\gamma )}}$ if $I\left(0\right)\le {\displaystyle \frac{\beta \mu k}{\alpha {\mu }_0(\mu +\gamma )}}.$ Thus, the set

$$\begin{array}{c}\hfill \mathcal{O}=\left\{(S,I,O)\in {\mathbb{R}}_{+}^3/0<S\le \frac{\mu }{\mu +\gamma },0\le I\le {\displaystyle \frac{\beta k\mu }{{\mu }_0\alpha (\mu +\gamma )}},0\le O\le \frac{k}{{\mu }_0}\right\},\end{array}$$

is positively invariant. However, if $S\left(0\right)>\frac{\mu }{\mu +\gamma }$, then either the solutions enters $\mathcal{O}$ in finite time or $S\left(t\right)$ approaches $\frac{\mu }{\mu +\gamma }$ asymptotically (similarly for $O\left(t\right)$ if $O\left(0\right)>\frac{k}{{\mu }_0},$ and an infectious mouse if $I\left(t\right)>{\displaystyle \frac{\beta \mu k}{\alpha {\mu }_0(\mu +\gamma )}}$). Hence, the region $\mathcal{O}$ attracts all solutions in ${\mathbb{R}}_{+}^3$. □

Consequently, it is sufficient to study the dynamics of the solutions of model (3) in this set, i.e., system (3) is mathematically well posed in $\mathcal{O}.$

3.1. Disease-Free Equilibrium Point

Steady states of epidemiological models are important since they provide insightful information related to the dynamics of the disease. We will investigate the impact of a time delay on the stability of the steady states of the system (2). A first step is to analyze the stability of the steady states of the model without a time delay. Later, we study the stability when the time delay is introduced. We also will determine the conditions such that the time delay changes the stability of the steady states.

The reduced model (3) without delay has a disease-free equilibrium point ($F_0^{*}$) and a endemic point ($EE$). The disease-free point

$$F_0^{*}=\left(\frac{\mu }{\mu +\gamma },0,0\right)\in \mathcal{O}.$$

(7)

The local stability of $F_0^{*}$ is determined by the eigenvalues of the Jacobian of the system (1) at $F_0^{*}$:

$$J\left(F_0^{*}\right)=\left(\begin{array}{ccc}-\mu -\gamma & -\mu & -\frac{\beta \mu }{\mu +\gamma }\\ 0& -\alpha & \frac{\beta \mu }{\mu +\gamma }\\ 0& k& -{\mu }_0\end{array}\right).$$

Clearly, the eigenvalues are given by $\lambda =-\mu -\gamma$ and the roots of the polynomial $p\left(\lambda \right)=(\lambda +\alpha )(\lambda +{\mu }_0)-\frac{k\beta \mu }{\mu +\gamma }$. We can write it as

$${\lambda }^2+a_1\lambda +b_1=0,$$

where $a_1=\alpha +{\mu }_0$ and $b_1=\alpha {\mu }_0-k\beta {\displaystyle \frac{\mu }{\mu +\gamma }}$. Clearly $a_1$ is positive and $b_1$ is positive if $k\beta {\displaystyle \frac{\mu }{\mu +\gamma }}<\alpha {\mu }_0$. Using the Routh–Hurwitz criterion, we have that the real part of the roots are negative under those conditions. Thus, we can define the basic reproduction number as

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

(8)

Another option to study the local stability is to use the next generation matrix method [56]. The system (3) with $\tau =0$ can be written as

$$\begin{array}{c}\hfill {\dot{\mathcal{W}}}_j\left(t\right)={\mathbb{F}}_j\left(\mathcal{W}\right)-{\mathbb{V}}_j^{-}\left(\mathcal{W}\right)+{\mathbb{V}}_j^{+}\left(\mathcal{W}\right),\end{array}$$

where $j=1,2,3$ and

$$\begin{array}{cc}\hfill & \mathbb{F}\left(\mathcal{W}\right)=\left(\begin{array}{c}{\mathbb{F}}_1\left(\mathcal{W}\right)\\ {\mathbb{F}}_2\left(\mathcal{W}\right)\\ {\mathbb{F}}_3\left(\mathcal{W}\right)\end{array}\right)=\left(\begin{array}{c}\beta S\left(t\right)O\left(t\right)\\ 0\\ 0\end{array}\right),{\mathbb{V}}^{+}\left(\mathcal{W}\right)=\left(\begin{array}{c}{\mathbb{V}}_1^{+}\left(\mathcal{W}\right)\\ {\mathbb{V}}_2^{+}\left(\mathcal{W}\right)\\ {\mathbb{V}}_3^{+}\left(\mathcal{W}\right)\end{array}\right)=\left(\begin{array}{c}0\\ kI\left(t\right)\\ \mu \end{array}\right),\hfill \\ & {\mathbb{V}}^{-}\left(\mathcal{W}\right)=\left(\begin{array}{c}{\mathbb{V}}_1^{-}\left(\mathcal{W}\right)\\ {\mathbb{V}}_2^{-}\left(\mathcal{W}\right)\\ {\mathbb{V}}_3^{-}\left(\mathcal{W}\right)\end{array}\right)=\left(\begin{array}{c}\alpha I\left(t\right)\\ {\mu }_{0}O\left(t\right)\\ \mu I\left(t\right)+\beta S\left(t\right)O\left(t\right)+(\mu +\gamma )S\left(t\right)\end{array}\right).\hfill \end{array}$$

We can see that

1.

For $\mathcal{W}\left(t\right)\ge 0,$ it follows that ${\mathbb{F}}_j\left(\mathcal{W}\right),{\mathbb{V}}_j^{-}\left(\mathcal{W}\right),{\mathbb{V}}_j^{+}\left(\mathcal{W}\right)\ge 0,$ for $j=1,2,3.$

2.

When $I\left(t\right)=0$ then ${\mathbb{V}}_1^{-}\left(\mathcal{W}\right)=0.$ Notice that when $\mathcal{W}=F_0^{*}={\left(\frac{\mu }{\mu +\gamma },0,0\right)}^t$ then ${\mathbb{V}}_1^{-}\left(F_0^{*}\right)=0.$

3.

${\mathbb{F}}_j\left(\mathcal{W}\right)=0,$ for $j>1.$

4.

If $\mathcal{W}=F_0^{*}={\left(\frac{\mu }{\mu +\gamma },0,0\right)}^t,$ then ${\mathbb{F}}_1\left(F_0^{*}\right)=0,$ and ${\mathbb{V}}_1^{+}\left(F_0^{*}\right)=0.$

5.

If ${\displaystyle \frac{k\beta \mu }{\mu +\gamma }}<\alpha {\mu }_0$, then the eigenvalues of Jacobian matrix of model (3) evaluated at $F_0^{*}$ with $\tau =0$ have negative real parts.

Accordingly, the next generation matrix is given by

$$\begin{array}{c}\hfill F{V}^{-1}=\left(\begin{array}{cc}{\displaystyle \frac{\beta \mu k}{(\mu +\gamma )\alpha {\mu }_0}}& {\displaystyle \frac{\beta \mu }{(\mu +\gamma ){\mu }_0}}\\ \\ 0& 0\end{array}\right).\end{array}$$

The characteristic polynomial of the next generation matrix $F{V}^{-1}$ is given by

$$\begin{array}{c}\hfill {\lambda }^2-{\displaystyle \frac{\beta \mu k\lambda }{(\mu +\gamma )\alpha {\mu }_0}}=0.\end{array}$$

Thus, the basic reproductive number is the dominant eigenvalue or spectral radius of the next generation matrix $F{V}^{-1}$

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

(9)

Thus, one obtains the following theorem,

Theorem 3.

The disease-free equilibrium $F_0^{*}$ of model (3) with $\tau =0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$, but unstable if ${\mathcal{R}}_0>1$.

Proof. 

By Theorem 2 given in [57], since conditions A1–A5 are satisfied (see conditions 1–5 above). □

3.2. Disease-Free Equilibrium Point Analysis for the Delay Model

The steady states of the model with time delay are the same as the model without delay. To study the stability of the steady states, we linearize the system (3). The characteristic equation is given by

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

where

$$J_0=\left(\begin{array}{ccc}-\mu -\gamma & -\mu & -\frac{\beta \mu }{\mu +\gamma }\\ 0& -\alpha & \frac{\beta \mu }{\mu +\gamma }\\ 0& 0& -{\mu }_0\end{array}\right),J_{\tau }=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& k& 0\end{array}\right).$$

Thus, one gets the transcendental characteristic equation for the disease-free equilibrium

$$\left|\lambda I-J_0-J_{\tau }{e}^{-\lambda \tau }\right|=(\lambda +\mu +\gamma )({\lambda }^2+a\lambda +b{e}^{-\lambda \tau }+c)=0,$$

(10)

where $a=\alpha +{\mu }_0,b=-{\displaystyle \frac{k\beta \mu }{\mu +\gamma }},$ and $c=\alpha {\mu }_0$.

Theorem 4.

The disease-free equilibrium $F_0^{*}$ of model (3) is locally asymptotically stable if ${\mathcal{R}}_0<1$, but unstable if ${\mathcal{R}}_0>1$.

Proof. 

We define from Equation(10) the following function $p\left(\lambda \right)={\lambda }^2+a\lambda +b{e}^{-\lambda \tau }+c$. Obviously, p is a continuous function. Moreover,

$$p\left(0\right)=b+c=-{\displaystyle \frac{k\beta \mu }{\mu +\gamma }}+\alpha {\mu }_0=1-{\mathcal{R}}_0,$$

and

$$\underset{\lambda \to \infty }{lim}p\left(\lambda \right)=+\infty .$$

Then, if ${\mathcal{R}}_0>1$, there is a positive real root, and then the disease free equilibrium is locally unstable. Let us consider the case when ${\mathcal{R}}_0<1$. Notice that when $\tau =0$, one obtains

$$g\left(\lambda \right)={\lambda }^2+a\lambda +1-{\mathcal{R}}_0.$$

Hence, all the roots of $g\left(\lambda \right)$ have negative real parts by the Routh–Hurwitz criterion.

We can see that $p\left(\lambda \right)$ does not have nonnegative real solutions because in $p\left(\lambda \right)$, it is increasing when $\lambda \ge 0$. Thus, if $p\left(\lambda \right)$ has roots with nonnegative real parts, they must be complex and should have been obtained from a pair of complex conjugate roots. Therefore, $p\left(\lambda \right)$ must have a pair of purely imaginary solutions.

Suppose that $\lambda =i\omega (\omega >0)$ is a root of $p\left(\lambda \right)$. Then, it must satisfy

$$-{\omega }^2+ai\omega +b{e}^{-i\omega \tau }+c=0.$$

Separating the real and imaginary part, one obtains

$$-{\omega }^2+bcos\left(\omega \tau \right)+c=0\Rightarrow -{\omega }^2+c=-bcos\left(\omega \tau \right),$$

and

$$a\omega -bsin\left(\omega \tau \right)=0\Rightarrow a\omega =bsin\left(\omega \tau \right).$$

Adding up the squares of the previous equations, one obtains that

$${\omega }^4+(a^2-2c){\omega }^2+(c^2-b^2)=0.$$

(11)

Thus, all the coefficients of Equation (11) are positive since the arithmetic mean is smaller than the geometric mean, and $c+b=1-{\mathcal{R}}_0>0.$ Consequently, Equation (11)) does not have positive real roots. Then, there is no $\omega$ such that $i\omega$ is a solution of Equation (10). Therefore, it follows from Rouché’s theorem [58], that the real parts of all the eigenvalues of the characteristic Equation (10) are negative for all $\tau >0$. This implies that disease-free equilibrium $F_0^{*}$ is locally asymptotically stable in $\mathcal{O}$ if ${\mathcal{R}}_0<1$ for the model (3). □

3.3. Global Stability of the Disease-Free Equilibrium $F_0^{*}$ for the Delay Model

Now, we proceed to study the global stability of the disease-free equilibrium $F_0^{*}$. This means that the extinction of toxoplasmosis is independent of the initial conditions of the subpopulations. We want to prove that if ${\mathcal{R}}_0<1,$ the disease-free equilibrium $F_0^{*}$ is globally asymptotically stable (GAS). We have the following theorem.

Theorem 5.

The disease-free equilibrium point $F_0^{*}$ of system (3) is globally asymptotically stable in $\mathcal{O}$ if ${\mathcal{R}}_0\le 1.$

Proof. 

To proceed with the global stability of the disease-free equilibrium $F_0^{*}$, we denote by $y_t$ the translation of the solution of the system (3), by $y_t=(S(t+\xi ),I(t+\xi ),V(t+\xi ))$ where $\xi \in [-\tau ,0]$. Consider the following Lyapunov functional:

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

(12)

At the disease-free equilibrium point $F_0^{*}$, we have that $S^{*}=\frac{\mu }{\mu +\gamma },I^{*}=0,O^{*}=0$, and $\mathcal{L}\left(F_0^{*}\right)=0$. The functional $\mathcal{L}$ also 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 than the disease-free equilibrium point $F_0^{*}$. Now, computing the time derivative of $\mathcal{L}\left(I\right(t),S(t),O(t\left)\right)$ along the trajectories of model (3), one obtains that

$$\begin{array}{cc}\hfill {\mathcal{L}}^{\prime }=& k\left[I\left(t\right)-I(t-\tau )\right]+O^{\prime }\left(t\right)=kI\left(t\right)-kI(t-\tau )+kI(t-\tau )-{\mu }_{0}O\left(t\right)\hfill \\ \hfill =& kI\left(t\right)-{\mu }_{0}O\left(t\right).\hfill \end{array}$$

However, we know from the restricted region $\mathcal{O}$ that $O\left(t\right)\le {\displaystyle \frac{k}{{\mu }_0}}$ and $0\le I\left(t\right)\le {\displaystyle \frac{\beta k\mu }{{\mu }_0\alpha (\mu +\gamma )}}={\mathcal{R}}_0$. Hence, $O\left(t\right){\mu }_o\le k$ and $kI\left(t\right)\le k{\mathcal{R}}_0$. Then,

$$kI\left(t\right)-O\left(t\right){\mu }_0\le k{\mathcal{R}}_0-k=k({\mathcal{R}}_0-1).$$

Thus, one obtains

$${\mathcal{L}}^{\prime }\le k({\mathcal{R}}_0-1).$$

Therefore, ${\mathcal{L}}^{\prime }\le 0$ if ${\mathcal{R}}_0\le 1$. By LaSalle’s invariance principle, this implies that the disease-free equilibrium $F_0^{*}$ is globally asymptotically stable in $\mathcal{O}$. This proves the theorem. □

3.4. Local Stability Analysis of the Endemic Point Equilibrium

The endemic equilibrium for the model (3) is

$$\begin{array}{c}\hfill E^{*}=\left(S^{*},I^{*},O^{*}\right),\end{array}$$

(14)

where

$$E^{*}=\left(S^{*}={\displaystyle \frac{\alpha {\mu }_0}{\beta k}},I^{*}={\displaystyle \frac{\mu \beta k-(\gamma +\mu )\alpha {\mu }_0}{\beta k(\mu +\alpha )}},O^{*}={\displaystyle \frac{\mu \beta k-(\gamma +\mu )\alpha {\mu }_0}{{\mu }_0\beta (\mu +\alpha )}}\right)\in \mathcal{O}.$$

(15)

This endemic equilibrium exists (positive components) only if ${\mathcal{R}}_0-1>0$. Now, we proceed to study the stability of the endemic point $E^{*}$ for the delayed model (3). Linearizing this model, one obtains the following transcendental characteristic equation

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

(16)

where

$$J_E=\left(\begin{array}{ccc}-\beta O^{*}-\mu -\gamma -\lambda & -\mu & -\beta S^{*}\\ \beta O^{*}& -\alpha -\lambda & \beta S^{*}\\ 0& 0& -{\mu }_0-\lambda \end{array}\right),\mathrm{and}{J}_{\tau }=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& k& 0\end{array}\right).$$

Thus, one obtains the characteristic equation

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

(17)

We can rewrite this characteristic equation as

$$\begin{array}{c}\hfill {\lambda }^3+A{\lambda }^2+B\lambda +C=e^{-\lambda \tau }(T_1\lambda +T_2),\end{array}$$

(18)

where

$$\begin{array}{cc}\hfill A& =\beta O^{*}+\alpha +\gamma +\mu +{\mu }_0,\hfill \\ \hfill B& =O^{*}\alpha \beta +\beta O^{*}\mu +\beta O^{*}{\mu }_0+\gamma \alpha +\gamma {\mu }_0+\alpha \mu +\alpha {\mu }_0+\mu {\mu }_0\hfill \\ \hfill & =O^{*}\beta (\alpha +\mu )+\beta O^{*}{\mu }_0+\alpha (\mu +\gamma )+\gamma {\mu }_0+\alpha {\mu }_0+\mu {\mu }_0\hfill \\ \hfill & =(\gamma +\mu )\alpha ({\mathcal{R}}_0-1)+\alpha (\mu +\gamma )+\beta O^{*}{\mu }_0+\gamma {\mu }_0+\alpha {\mu }_0+\mu {\mu }_0\hfill \\ \hfill & =\alpha (\mu +\gamma ){\mathcal{R}}_0+\beta O^{*}{\mu }_0+\gamma {\mu }_0+\alpha {\mu }_0+\mu {\mu }_0,\hfill \\ \hfill C& =O^{*}\alpha \beta {\mu }_0+O^{*}\beta \mu {\mu }_0+\alpha \gamma {\mu }_0+\alpha {\mu }_0\mu =O^{*}\beta {\mu }_0(\alpha +\mu )+\alpha {\mu }_0(\gamma +\mu )\hfill \\ \hfill & =(\gamma +\mu )\alpha ({\mathcal{R}}_0-1){\mu }_0+\alpha {\mu }_0(\gamma +\mu )=(\gamma +\mu )\alpha {\mathcal{R}}_0{\mu }_0,\hfill \\ \hfill T_1& =S^{*}\beta k=\alpha {\mu }_0,\hfill \\ \hfill T_2& =S^{*}\beta k(\gamma +\mu )=\alpha {\mu }_0(\gamma +\mu ),\hfill \end{array}$$

where $A,B,C,T_1,T_2\ge 0$. Rewriting Equation (18), one obtains

$$\begin{array}{cc}\hfill {\lambda }^3+A{\lambda }^2+\left(\alpha (\mu +\gamma ){\mathcal{R}}_0+\beta O^{*}{\mu }_0+\gamma {\mu }_0+\alpha {\mu }_0+\mu {\mu }_0\right)\lambda & \hfill \\ \hfill =e^{-\lambda \tau }(T_1\lambda +T_2)-C& \\ \hfill =e^{-\lambda \tau }\alpha {\mu }_0\lambda +e^{-\lambda \tau }\alpha {\mu }_0(\gamma +\mu )-(\gamma +\mu )\alpha {\mathcal{R}}_0{\mu }_0& \\ \hfill =e^{-\lambda \tau }\alpha {\mu }_0\lambda +\alpha {\mu }_0(\gamma +\mu )[e^{-\lambda \tau }-{\mathcal{R}}_0].\end{array}$$

Now we can rewrite it as

$$\begin{array}{cc}\hfill {\lambda }^3+A{\lambda }^2+& \left(\alpha (\mu +\gamma ){\mathcal{R}}_0+\beta O^{*}{\mu }_0+\gamma {\mu }_0+\mu {\mu }_0\right)\lambda \hfill \\ \hfill =& e^{-\lambda \tau }\alpha {\mu }_0\lambda -\alpha {\mu }_0\lambda +\alpha {\mu }_0(\gamma +\mu )[e^{-\lambda \tau }-{\mathcal{R}}_0]\hfill \\ \hfill =& \alpha {\mu }_0\lambda [e^{-\lambda \tau }-1]+\alpha {\mu }_0(\gamma +\mu )[e^{-\lambda \tau }-{\mathcal{R}}_0].\hfill \end{array}$$

Thus, one obtains

$$\begin{array}{c}\hfill {\lambda }^3+\tilde{A}{\lambda }^2+\tilde{B}\lambda =\widetilde{T_1}\lambda +\widetilde{T_2},\end{array}$$

(19)

where

$$\begin{array}{cc}\hfill \tilde{A}& =A\ge 0,\hfill \\ \hfill \tilde{B}& =\alpha (\mu +\gamma ){\mathcal{R}}_0+\beta O^{*}{\mu }_0+\gamma {\mu }_0+\mu {\mu }_0\ge 0,\hfill \\ \hfill \widetilde{T_1}& =\alpha {\mu }_0[e^{-\lambda \tau }-1]\le 0,\hfill \\ \hfill \widetilde{T_2}& =\alpha {\mu }_0(\gamma +\mu )[e^{-\lambda \tau }-{\mathcal{R}}_0].\hfill \end{array}$$

Therefore, if ${\mathcal{R}}_0>1$, then $\widetilde{T_2}\le 0$. Thus, the left-hand side of Equation (19) is positive and the right-hand side is negative for all $\lambda \ge 0$. Therefore, for any $\lambda \ge 0$, Equation (18) cannot have real non-negative solutions.

The next step is to find the distribution of the roots of Equation (18) when $\tau >0$. To discard complex conjugate solutions with nonnegative real parts, let us assume that $\lambda =i\omega (\omega >0)$ is a root of Equation (18). This occurs if and only if $\omega$ satisfies the next equation:

$$\begin{array}{cc}\hfill -i{\omega }^3-A{\omega }^2+Bi\omega +C& =e^{-i\omega \tau }[i\omega T_1+T_2]=(cos\left(\omega \tau \right)-isin\left(\omega \tau \right))[i\omega T_1+T_2].\hfill \end{array}$$

Separating the imaginary and real part, one obtains

$$-{\omega }^3+B\omega =T_1\omega cos\left(\omega \tau \right)-T_{2}sin\left(\omega \tau \right),$$

(20)

and

$$-A{\omega }^2+C=T_{2}cos\left(\omega \tau \right)+T_1\omega sin\left(\omega \tau \right).$$

(21)

Squaring the previous equations, one obtains

$${\omega }^6-2B{\omega }^4+B^2{\omega }^2=T_1^2{\omega }^2{cos}^2\left(\omega \tau \right)-2T_1{T}_2\omega cos\left(\omega \tau \right)sin\left(\omega \tau \right)+T_2^2{sin}^2\left(\omega \tau \right),$$

and

$$A^2{\omega }^4-2AC{\omega }^2+C^2=T_2^2{cos}^2\left(\omega \tau \right)+2T_1{T}_2\omega cos\left(\omega \tau \right)sin\left(\omega \tau \right)+T_1^2{\omega }^2{sin}^2\left(\omega \tau \right).$$

Adding both equations, we obtain that

$${\omega }^6+(A^2-2B){\omega }^4+(B^2-2AC){\omega }^2+C^2=T_1^2{\omega }^2+T_2^2.$$

Rearranging,

$${\omega }^6+(A^2-2B){\omega }^4+(B^2-2AC-T_1^2){\omega }^2+C^2-T_2^2=0.$$

(22)

Now let $z={\omega }^2$, then

$$z^3+(A^2-2B)z^2+(B^2-2AC-T_1^2)z+C^2-T_2^2=0.$$

(23)

Thus

$$z^3+\alpha z^2+\beta z+\nu =0,$$

(24)

where $\alpha =A^2-2B$, $\beta =B^2-2AC-T_1^2$ and $\nu =C^2-T_2^2$.

Lemma 1.

If $\alpha \ge 0$, $\beta \ge 0$ and $\nu \ge 0$, then Equation (24) has no positive real roots.

Proof. 

Let $h\left(\xi \right)={\xi }^3+\alpha {\xi }^2+\beta \xi +\nu .$ Deriving with respect to $\xi$, one obtains that $h^{\prime }\left(\xi \right)=3{\xi }^2+2\alpha \xi +\beta$. Here $\xi ={\omega }^2>0$, $\alpha \ge 0$, $\beta \ge 0$ and $\nu \ge 0$. Thus, $h^{\prime }\left(\xi \right)>0$, so $h\left(\xi \right)$ is an increasing function of $\xi \ge 0$. Therefore, since $h\left(0\right)=\nu \ge 0$, then Equation (24) has no positive real roots. □

Therefore, there is no $\lambda =i\omega (\omega >0)$ such that it is a solution of Equation (18) if the previous conditions are satisfied. Accordingly, by Rouché’s theorem [58], the real parts of all the eigenvalues of Equation (18) are negative for all values of the time delay $\tau \ge 0$ under the previous conditions.

3.5. Bifurcation Analysis

In the previous subsection, we found conditions such that the endemic equilibrium $E^{*}$ is stable. If these conditions are not fulfilled, then the stability of the endemic equilibrium depends on the time delay $\tau$. Therefore, the endemic equilibrium point can lose stability, and periodic solutions can appear. For instance, if the parameter $\nu <0$, then the polynomial related to Equation (24) has at least one positive root $z_0$, and then the endemic equilibrium is no longer stable.

For the bifurcation analysis, we will use the time delay $\tau$ since we are interested in understanding the effect of a time delay in the infectivity of the oocysts. For instance, some chemicals can be developed such that the infectivity of the oocysts in the environment is delayed. It has been shown that Gamma irradiation of oocysts changes their infectivity [15].

Let us see how the roots of Equation (18) vary when the bifurcation parameter $\tau$ varies. Let $\lambda \left(\tau \right)=\eta \left(\tau \right)+i\omega \left(\tau \right)$ be one eigenvalue of Equation (18) such that for some particular value of the bifurcation parameter ${\tau }_0$, we have that the real part is zero, i.e., $\eta \left({\tau }_0\right)=0$, and the imaginary part $\omega \left({\tau }_0\right)={\omega }_0$. Since the roots come in pairs, we can assume that without loss of generality, ${\omega }_0>0$. We propose the following theorem.

Theorem 6.

Assume that ${\mathcal{R}}_0>1$, then the endemic equilibrium $E^{*}$ of the delay model (3) is asymptotically stable for all $\tau \ge 0$ and, therefore, Hopf bifurcation cannot occur.

Proof. 

We can proceed as before using Equation (18) to obtain Equation (23). Now computing $\alpha$, $\beta$ and $\nu$ of Equation (24), one obtains that

$$\begin{array}{cc}\hfill & \alpha ={\nu }^2+2(2\beta O^{*}+\mu )\nu +{\beta }^2{\left(O^{*}\right)}^2+{\mu }^2+{\alpha }^2+{\mu }_0^2>0,\hfill \\ \hfill & \beta ={\displaystyle \frac{\left(q+({\mathcal{R}}_0^2{\alpha }^3+2{\mathcal{R}}_0^2{\alpha }^2\nu +{\mu }_0^2({\mathcal{R}}_0^2-2{\mathcal{R}}_0+2)\alpha +2{\mathcal{R}}_0\nu {\mu }_0^2)\alpha \mu +{\mathcal{R}}_0^2{\alpha }^2\nu \left({\alpha }^2{\mu }_0^2\right)\right)(\nu +\mu )}{{(\alpha +\mu )}^2}}>0,\hfill \end{array}$$

with $q=({\mathcal{R}}_0^2{\alpha }^2+{\mu }_0^2)\left({\mu }^3+2(\alpha +\nu /2){\mu }_0\right),$ and

$$\nu ={\alpha }^2{(\nu +\mu )}^2{\mu }_0^2({\mathcal{R}}_0-1).$$

It follows that $\alpha =A^2-2B>0$, $\beta =B^2-2AC-T_1^2>0$ and $\nu ={\alpha }^2{(\nu +\mu )}^2{\mu }_0^2({\mathcal{R}}_0-1)>0$ if ${\mathcal{R}}_0>1$. Therefore, using Lemma 1, there is no $\lambda =i\omega (\omega >0)$ such that it is a solution of Equation (18). Accordingly, by Rouché’s theorem [58], the real parts of all the eigenvalues of Equation (18) are negative for all values of the time delay $\tau \ge 0$. Thus, Hopf bifurcation cannot occur since the next transversal condition is not satisfied for any $\tau$

$${\displaystyle \frac{d\mathfrak{Re}\left(\lambda \right(\tau \left)\right)}{d\tau }}{|}_{\tau ={\tau }_0}>0.$$

(25)

□

Thus, this Theorem indicates that if ${\mathcal{R}}_0>1$, there are no values of the parameter $\tau$ such that the system (3) exhibits a Hopf bifurcation. Summarizing the above results, we prove that the introduction of a time delay $\tau$ in the infectivity of the oocysts cannot cause the endemic state $E^{*}$ of the delayed model (3) to become unstable. Moreover, it cannot cause periodic solutions, due to Hopf bifurcation. In the next section, we perform numerical simulations that are in good agreement with this theoretical result.

4. Numerical Simulations

In this section, we present a variety of numerical simulations of the mathematical model (3) in order to explore the dynamics of toxoplasmosis and support the theoretical results obtained in the previous section. We explore a variety of scenarios with different time delays and combining different parameter values in order to obtain values of the basic reproduction number ${\mathcal{R}}_0$ greater and less than one since this is an important threshold for the outcome of the disease. These simulations provide insights related to the dynamics of toxoplasmosis at the population level. We vary the values of the infectivity of the oocysts and the time delay $\tau$ under the assumption that it can be modified by ecologically friendly products that affect the infectivity of the oocysts. However, modifying the time when oocysts become infectious might be difficult in the real world [59].

We use Matlab to compute the numerical simulations (on Intel(R) Core(TM) i7-8700 CPU 3.20 GHz and 16 GB RAM). In particular, we use the function (or solver) dde23, which numerically solves delay differential equations [60,61]. For all the numerical simulations, we use weeks as the time unit.

For most of the the numerical simulations, we use the base parameter values presented in

Table 1. Description of the parameters of the mathematical model (3) and values.

Parameter	Description	Value
$\mu$	Birth/Death rates (cats)	$1/260$ (1/weeks) [62]
$\alpha$	Shedding period	$1/2$ (1/weeks) [6]
${\mu }_0$	Clearance rate	$1/26$ (1/day) [6,44]
k	Oocysts per day (cat)	$20\times {10}^6$ (1/day) [16]
$\beta$	Transmission rate	varied
$\gamma$	Vaccination rate	varied

. A variety of initial conditions are used in order to support the theoretical results regarding local and global stability. We compute the equilibrium points for each numerical simulation and compare them with the theoretical result.

4.1. Numerical Simulations without Time Delay

We assume a transmission rate $\beta$ such that the basic reproduction number ${\mathcal{R}}_0<1$. Figure 2 shows the dynamics of the subpopulations when there is no time delay ($\tau =0$). In this case, the system approaches the disease-free equilibrium $F^{*}$ as expected from the theoretical results. Figure 3 shows the dynamics of the subpopulations when ${\mathcal{R}}_0>1$ and there is no time delay ($\tau =0$). In this case, the system does not approach the disease-free equilibrium $F^{*}$ as the theoretical results indicate.

4.2. Numerical Simulations with Time Delay

First, we consider a transmission rate $\beta$ such that the basic reproduction number ${\mathcal{R}}_0<1$. Figure 4 shows that in this case, the system still approaches the disease-free equilibrium $F^{*}$ despite the introduction of a time delay on the infectivity of the oocysts. We add an additional numerical simulation with a large time delay ($\tau =50$) to observe that the solution still approaches the disease-free equilibrium $F^{*}$. Figure 5 shows the graph of the solution for this previous scenario. The theoretical results indicate that whenever ${\mathcal{R}}_0<1$ and for any initial condition, the disease will disappear due to the global stability. This is quite important since it guarantees that the disease can become extinct. Figure 6 shows the dynamics of the subpopulations when ${\mathcal{R}}_0>1$. Notice that the system does not approach the disease-free equilibrium $F^{*}$. This is due to an increase in the infectivity of the oocysts that also increases the transmission rate $\beta$.

4.3. Numerical Simulations without Hopf Bifurcation

Now, we consider a transmission rate $\beta$ such that the basic reproduction number ${\mathcal{R}}_0>1$, which allows the system to approach the endemic steadiness. Oftentimes, increasing the time delay can give conditions such that periodic solutions arise and stability is no longer present. This is crucial since we consider that the time delay $\tau$ can be modified by ecologically friendly products that affect the infectivity of the oocysts. However, by the theoretical results, we know that, even for large time delays, periodic solutions will not arise, and the stability of the endemic equilibrium will be maintained.

Figure 7 shows the dynamics when $\tau =300$. In this case, the system still approaches the endemic steady state $E^{*}$ despite the introduction of a very large time delay. The theoretical results indicate that whenever ${\mathcal{R}}_0>1$, the endemic steady state $E^{*}$ is stable. This is important since it means that despite the application of ecologically friendly chemicals to increase the delay of when the oocysts become infectious, the disease persists. Figure 8 shows the dynamics of the subpopulations when there is no time delay and ${\mathcal{R}}_0>1$. Finally, we take initial conditions that are very close to the disease-free equilibrium (far from the endemic point) and ${\mathcal{R}}_0>1$. In this case, Figure 9 shows that the system still approaches the endemic point despite a very long time delay of $\tau =300$. This result is in good agreement with the theoretical results.

4.4. Effect of the Time Delay

In this subsection, we perform additional numerical simulations to see the effect of the time delay. From the theoretical results, we obtain the conditions such that the solutions approach one of the steady states. However, the numerical simulations provide the transient behavior of the solutions. Figure 10 shows a comparison of the solutions of models (1) and (3), i.e., without delay and with delay. In this case, we consider that the basic reproduction number ${\mathcal{R}}_0>1$, so both solutions approach the endemic steady state $E^{*}$, and this is in good agreement with the theoretical results. The transient behavior of the solutions differs due to the time delay. It can be observed that the peaks of infection occur later on the delayed model and that both systems approach the endemic steady state $E^{*}$. Figure 11 shows the solutions when the time delay is increased to $\tau =30$. Again, the infection peak of the delayed model comes later than the model without delay. Moreover, the peak of infected cats of the delayed model is lower. Thus, the time delay acts as a damping factor on the epidemics when the basic reproduction number ${\mathcal{R}}_0>1$. One additional simulation is performed with a very large time delay $\tau =300$. Figure 12 shows that the delayed model has lower infection peaks but still approaches the correct steady state despite the long time delay, which is more challenging for the numerical scheme.

5. Conclusions

In this paper, we constructed a mathematical model to study toxoplasmosis dynamics. The model considers a time delay from when the oocysts become present in the environment to when they are able to infect. We investigated under what conditions the toxoplasmosis can be eradicated and the impact of the time delay on the dynamics. We calculated the basic reproduction number ${\mathcal{R}}_0$, which is an important parameter that determines the global dynamics of the toxoplasmosis. We found that if ${\mathcal{R}}_0<1$, the disease-free equilibrium is globally stable with or without a time delay, and toxoplasmosis can be eradicated. If ${\mathcal{R}}_0>1$, we found that there is only one endemic equilibrium and we proved that this is locally stable when there is no time delay. Moreover, we found conditions such that the endemic steady state cannot become unstable and Hopf bifurcation cannot arise. This means that under any conditions or any value of the time delay, the toxoplasmosis will persist over time if the basic reproduction number ${\mathcal{R}}_0>1$. Thus, we proved that the introduction of the time delay cannot destabilize the delayed model, and periodic solutions due to Hopf bifurcation cannot appear. Numerical simulations were performed, and they are in good agreement with the theoretical results.

The results presented in this study allow us to discuss the impact of vaccination programs against toxoplasmosis and the time delay related to the infectivity of the oocysts on the dynamics of toxoplasmosis in cat populations. The mathematical analysis provided in this article is applicable for any parameter values. The values of some parameters were taken from specialized works related to toxoplasmosis, despite there being some uncertainty with them as is common in many biological processes. The basic reproduction number found here determines the impact of the time delay on the dynamics of the toxoplasmosis on the population of cats. Based on the results, we observed that the time delay cannot affect the steady states of the biological system. Moreover, the time delay cannot change the stability of the endemic state, which is an important result from a biological point of view. It is important to remark that temperatures affect the time delay or when the oocysts become infective, and thus, with climate change, the time delay can be affected [39].

The constructed mathematical model also allows the evaluation of prevention and control strategies designed by health institutions against toxoplasmosis. Health authorities can implement a vaccination program for cats and/or use ecologically friendly products to reduce the amount of oocysts and/or the time when the oocysts become infective. Studies using optimal control are recommended to evaluate the feasibility of these strategies, taking into account economic and ecological factors [63,64].

Finally, we would like to mention some limitations of this work and possible future work. As with any mathematical model for biological processes, there are limitations due to the complexity of the real world. The proposed model does not consider other intermediate hosts, such as humans and mice. The classical well-known homogeneous mixing assumption used in most epidemiological mathematical models might not be appropriate for some regions, but further research is needed. Creating different models for rural and urban areas might be necessary. In this work, biological values for some parameters were used, but there are still unknowns about the crucial transmission rate. Further studies that consider additional hosts and the introduction of a latent class or additional time delays could be profitable, although challenging.