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

: 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 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 ﬁnd 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.


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].
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 R 0 . In Section 4, the numerical simulations of the mathematical model are presented. Section 5 contains the conclusions.

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(t) is divided into three subpopulations: susceptible cats (S(t)), infected cats (I(t)), and vaccinated/recovered cats (V R (t)). Thus, the model is a first order nonlinear system of ordinary differential equations given byṠ 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(t) = S(t) + I(t) + V R (t) = 1. The mathematical diagram is depicted in Figure 1. We study the dynamics of the model (1) in the restricted region: Now, we introduce a time delay τ > 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 becomeṡ with S(t) + I(t) + V R (t) = 1. Since the populations are normalized, we can study the following reduced model where V R (t) = 1 − S(t) − I(t). In addition, the system (3) satisfies the initial conditions given by with ξ(s) continuous function defined from the interval [−τ, 0] to R + , and with the norm ξ = sup −τ≤s≤0 |ξ(s)|. 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Ẇ where with the initial conditions given in (4). We use the following theorem: Suppose Ω is an open set in R × C [−τ, 0], R n + , f : Ω → R n is continuous, and f (t, ϕ) is Lipschitizian in ϕ in each compact set in Ω. If (σ, ϕ) ∈ Ω then there is a unique solution of model (5) through (σ, ϕ).
Let S = C [−τ, 0], R 3 + be the Banach space of continuous functions mapping the interval [−τ, 0] into R 3 + which is equipped with the sup-norm. The function f (t, W (t)) given by in (5) is defined in f : [0, ∞) × C [−τ, 0], R 3 + → R 3 is continuous and satisfies a local Lipschitz condition. Indeed, for ϕ 1 , Then, by Theorem 1 for any ϕ ∈ S there exists a unique solution W (t, ϕ) of the system (3) such that W 0 = ϕ. 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).

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 ≥ 0 and that they are bounded. These features are important from a realistic viewpoint. Proof. Using the first Equation of (3), we can see thaṫ Therefore, it follows that for all t ≥ 0. On the other hand, suppose that there exists a t 1 > 0 such that I(t 1 ) = 0,İ(t 1 ) ≤ 0, and I(t) > 0 for all t ∈ (0, t 1 ). Thus, from the second equation of model (3), one obtains that However, O(t) > 0 for all t ∈ (0, t 1 ). If this is not the case, there exists a t 2 > 0 such that t 2 < t 1 ,Ȯ(t 2 ) ≤ 0, O(t 2 ) = 0 and O(t) > 0 for t ∈ (0, t 2 ). From the third equation of system (3), it follows that 0 ≥ kI(t 2 − τ). Since τ > 0 and using the continuity of solution, then one gets that 0 ≥ lim τ→0 I(t 2 − τ) = I(t 2 ). This is a contradiction because I(t) > 0 for all t ∈ (0, t 1 ). Thus, using the continuity of O(t), it follows that O(t) > 0, with t ∈ (0, t 1 ]. Therefore, (6) is false and thus I(t) ≥ 0 for t ≥ 0. From the third equation, it is concluded that O(t) ≥ 0 for t ≥ 0. Now, the subpopulation S(t) is bounded by µ µ+γ , because by the standard comparison theorem [55], we can show that .
. Thus, the set ). Hence, the region O attracts all solutions 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 O.

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 The local stability of F * 0 is determined by the eigenvalues of the Jacobian of the system (1) at F * 0 : Clearly, the eigenvalues are given by λ = −µ − γ and the roots of the polynomial We can write it as 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 Another option to study the local stability is to use the next generation matrix method [56]. The system (3) with τ = 0 can be written aṡ where j = 1, 2, 3 and We can see that When

5.
If kβµ µ + γ < α µ 0 , then the eigenvalues of Jacobian matrix of model (3) evaluated at F * 0 with τ = 0 have negative real parts. Accordingly, the next generation matrix is given by The characteristic polynomial of the next generation matrix FV −1 is given by Thus, the basic reproductive number is the dominant eigenvalue or spectral radius of the next generation matrix FV −1 Thus, one obtains the following theorem, Proof. By Theorem 2 given in [57], since conditions A1-A5 are satisfied (see conditions 1-5 above).

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 Thus, one gets the transcendental characteristic equation for the disease-free equilibrium Proof. We define from Equation (10) Obviously, p is a continuous function. Moreover, Then, if R 0 > 1, there is a positive real root, and then the disease free equilibrium is locally unstable. Let us consider the case when R 0 < 1. Notice that when τ = 0, one obtains Hence, all the roots of g(λ) have negative real parts by the Routh-Hurwitz criterion. We can see that p(λ) does not have nonnegative real solutions because in p(λ), it is increasing when λ ≥ 0. Thus, if p(λ) has roots with nonnegative real parts, they must be complex and should have been obtained from a pair of complex conjugate roots. Therefore, p(λ) must have a pair of purely imaginary solutions.
Separating the real and imaginary part, one obtains Adding up the squares of the previous equations, one obtains that Thus, all the coefficients of Equation (11) are positive since the arithmetic mean is smaller than the geometric mean, and c + b = 1 − R 0 > 0. Consequently, Equation (11)) does not have positive real roots. Then, there is no ω such that iω 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 τ > 0. This implies that disease-free equilibrium F * 0 is locally asymptotically stable in O if 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 R 0 < 1, the disease-free equilibrium F * 0 is globally asymptotically stable (GAS). We have the following theorem.
Therefore, L ≤ 0 if R 0 ≤ 1. By LaSalle's invariance principle, this implies that the disease-free equilibrium F * 0 is globally asymptotically stable in O. This proves the theorem.

Local Stability Analysis of the Endemic Point Equilibrium
The endemic equilibrium for the model (3) is where This endemic equilibrium exists (positive components) only if 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 where Thus, one obtains the characteristic equation We can rewrite this characteristic equation as where where A, B, C, T 1 , T 2 ≥ 0. Rewriting Equation (18), one obtains

Now we can rewrite it as
Thus, one obtains Therefore, if R 0 > 1, thenT 2 ≤ 0. Thus, the left-hand side of Equation (19) is positive and the right-hand side is negative for all λ ≥ 0. Therefore, for any λ ≥ 0, Equation (18) cannot have real non-negative solutions.
The next step is to find the distribution of the roots of Equation (18) when τ > 0. To discard complex conjugate solutions with nonnegative real parts, let us assume that λ = i ω (ω > 0) is a root of Equation (18). This occurs if and only if ω satisfies the next equation: Separating the imaginary and real part, one obtains and Squaring the previous equations, one obtains and A 2 ω 4 − 2AC ω 2 + C 2 = T 2 2 cos 2 (ωτ) + 2T 1 T 2 ω cos(ωτ) sin(ωτ) + T 2 1 ω 2 sin 2 (ωτ).
Therefore, there is no λ = iω (ω > 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 τ ≥ 0 under the previous conditions.

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 τ. Therefore, the endemic equilibrium point can lose stability, and periodic solutions can appear. For instance, if the parameter ν < 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 τ 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 τ varies. Let λ(τ) = η(τ) + iω(τ) be one eigenvalue of Equation (18) such that for some particular value of the bifurcation parameter τ 0 , we have that the real part is zero, i.e., η(τ 0 ) = 0, and the imaginary part ω(τ 0 ) = ω 0 . Since the roots come in pairs, we can assume that without loss of generality, ω 0 > 0. We propose the following theorem. Theorem 6. Assume that R 0 > 1, then the endemic equilibrium E * of the delay model (3) is asymptotically stable for all τ ≥ 0 and, therefore, Hopf bifurcation cannot occur.
Proof. We can proceed as before using Equation (18) to obtain Equation (23). Now computing α, β and ν of Equation (24), one obtains that with q = (R 2 0 α 2 + µ 2 0 ) µ 3 + 2(α + ν/2)µ 0 , and Therefore, using Lemma 1, there is no λ = iω (ω > 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 τ ≥ 0. Thus, Hopf bifurcation cannot occur since the next transversal condition is not satisfied for any τ Thus, this Theorem indicates that if R 0 > 1, there are no values of the parameter τ such that the system (3) exhibits a Hopf bifurcation. Summarizing the above results, we prove that the introduction of a time delay τ 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.

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 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 τ 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. 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.

Numerical Simulations without Time Delay
We assume a transmission rate β such that the basic reproduction number R 0 < 1. Figure 2 shows the dynamics of the subpopulations when there is no time delay (τ = 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 R 0 > 1 and there is no time delay (τ = 0). In this case, the system does not approach the disease-free equilibrium F * as the theoretical results indicate.

Numerical Simulations with Time Delay
First, we consider a transmission rate β such that the basic reproduction number 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 (τ = 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 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 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 β.

Numerical Simulations without Hopf Bifurcation
Now, we consider a transmission rate β such that the basic reproduction number 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 τ 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 τ = 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 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 R 0 > 1. Finally, we take initial conditions that are very close to the disease-free equilibrium (far from the endemic point) and R 0 > 1. In this case, Figure 9 shows that the system still approaches the endemic point despite a very long time delay of τ = 300. This result is in good agreement with the theoretical results.

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 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 τ = 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 R 0 > 1. One additional simulation is performed with a very large time delay τ = 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.

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 R 0 , which is an important parameter that determines the global dynamics of the toxoplasmosis. We found that if R 0 < 1, the disease-free equilibrium is globally stable with or without a time delay, and toxoplasmosis can be eradicated. If 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 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.  Acknowledgments: The authors are grateful to the reviewers for their careful reading of this manuscript and their useful comments to improve the content of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.