A Mathematical Model for Within-host Toxoplasma Gondii Invasion Dynamics

Toxoplasma gondii (T. gondii) is a protozoan parasite that infects a wide range of intermediate hosts, including all mammals and birds. Up to 20% of the human population in the US and 30% in the world are chronically infected. This paper presents a mathematical model to describe intra-host dynamics of T. gondii infection. The model considers the invasion process, egress kinetics, interconversion between fast-replicating tachyzoite stage and slowly replicating bradyzoite stage, as well as the host's immune response. Analytical and numerical studies of the model can help to understand the influences of various parameters to the transient and steady-state dynamics of the disease infection. A compartmental model representing the dynamics of T. gondii.


Introduction
Toxoplasma gondii, often referred to as T. gondii, is a parasite that is able to infect a wide range of hosts, including all mammals and birds [8].Up to one third of the world's human population and about 20% of the population in the US are estimated to carry a Toxoplasma infection [1].People can be infected by eating infected meat, by drinking water contaminated with the parasite, or by transmission from mother to fetus.During acute Toxoplasma infection, the patient typically exhibits a mild flu-like illness (swollen lymph nodes, or muscle aches and pains that last for a month or more) or no illness at all.However, people with a weakened immune system, such as AIDS patients or recipients of chemotherapy and organ transplant, may develop serious inflammation in the brain or the eyes.In most immunocompetent patients, the infection enters a latent phase, during which tissue cysts may form in the brain and muscle.Recent studies show that latent Toxoplasmosis may have significant effects on human behavior and may lead to neuropsychiatric disorders, e.g.schizophrenia.In addition, infection acquired during pregnancy may spread and cause severe damage to the fetus [1].
T. gondii has a complex life cycle, as seen in Figure 1.The parasite uses the feline to reproduce sexually.When the cat becomes infected, it sheds oocysts, which infect the environment.These oocysts can be ingested by mammals and birds which then become infected with the parasite [8].Eating another organism that is infected can also infect the secondary hosts.Within a host, T. gondii exists in two interconvertable stages: bradyzoites and tachyzoites.Bradyzoites have the slow-growing and encysted form whereas tachyzoites are the fast-replicating parasites.Tachyzoites disseminate within the host and lead to the acute phase of infection.After the bradyzoite-containing cysts are ingested by the host, the walls of these cysts are digested inside the host's stomach.Bradyzoites, which are resistant to gastric conditions in the stomach, will subsequently invade the host's epithelial cells of the small intestine and convert into tachyzoites there.While most of tachyzoites in immunocompetent hosts are eliminated by the innate and adaptive immune responses, some tachyzoites differentiate into the dormant bradyzoite stage inside host cells [8].The differentiation of tachyzoites into the bradyzoite stage plays an essential role in the development of tissue cysts, which allows life-long persistence of the parasites in the host.Reactivation of bradyzoites back to tachyzoites can lead to life threatening infection.The interconversion between tachyzoites and bradyzoites can be influenced by many in vivo and in vitro factors.
In this work, we aim to develop a mathematical model to understand the nonlinear, complex interactions between T. gondii invasion dynamics and host immune response.We develop a compartmental model to describe the invasion dynamics of T. gondii ; see Figure 2. The model has 7 state variables: the population size of uninfected cells, X; the population of cells infected with tachyzoites, Y T ; the population of cells containing early-stage bradyzoites, Y B ; the population of cells containing encysted bradyzoites, Y C ; the population of free tachyzoites, P T ; the population of free bradyzoites, P B ; and the effector cells of the host's immune response, Z.

Model Development
We assume uninfected cells are generated at a constant rate of λ and assume the average life time of an uninfected cell is 1/d.In the absence of an infection, the population dynamics of host cells is given by Ẋ = λ − dX.Under this simple population dynamics model, the number of uninfected cells converges to the equilibrium X 0 = λ/d.Free parasites infect uninfected cells at a rate proportional to the product of their abundance: β P T XP T for tachyzoites and β P B XP B for bradyzoites.The rate constants, β P T and β P B , are lumped parameters, describing the efficacy of the invasion process and depending on the rate at which the parasites find uninfected cells, the rate of parasite entry, and the probability of successful infection.
Infected cells produce free parasites at a rate proportional to their abundance: k T y T for tachyzoites and k C y C for encysted bradyzoites.Average life time of a cell infected with tachyzoites is 1/a T and that of a cell infected with encysted bradyzoites is 1/a C .Since tachyzoites replicate much faster than bradyzoites, one expect a C be much smaller than a T .Free parasites are removed from the system at a rate u T P T for tachyzoites and at rate u B P B for bradyzoites.Tachyzoites in a host cell can spontaneously convert to bradyzoites at a rate r T y T .To account for the reactivation process, we assume early-stage bradyzoites may convert to tachyzoites at a rate r B y B .We also consider a simple model of immune response.Much work has been done on immune response from Toxoplasma infection and many important mechanisms have been identified [11].Here, we introduce a variable Z to represent the overall effector cells without consideration of specific immune mechanisms.We assume the effector cells act on host cells infected with tachyzoites in a predator-prey manner.
Combining the above processes leads to the following system of equations: where ρ is the production rate of the effector cells, δ is the removal rate of the immune system response.The first term in the immune equation represents the activation process in response to the detection of infected cells whereas the second term in the immune equation represents natural decay of the immune effector.The activation process takes the form of Holling's Type II predator-prey relation [15,18,3].When the number of infected cells is small, the level of immune response is low.Then, the immune response increases at a great rate and saturates when the number of parasites is sufficiently large.Kafsack et al. [14] developed a mathematical model to interpret kinetics data collected for T. gondii invasion.Their results show that T. gondii invasion dynamics including contact, attaching, penetrating, and invasion occur within a few minutes.On the other hand, previous experiments show that replication and stage conversion dynamics take place in hours [13,8,20].We assume that the kinetics of free parasites are significantly faster than kinetics of other processes.Thus, we can adopt the common quasistatic approximation and assume the free parasites are in equilibrium.It follows from Equations ( 5)-( 6) that Substituting the above equations into Equations ( 1)-( 4), (7) leads to the following reduced model: where 3 Results

Disease-Free Equilibrium (DFE)
The structure of system ( 10)-( 14) implies there exists a unique non-negative DFE solution.
Denote this equilibrium solution by The stability of E 0 can be established using the next generation operator method on the system (10)- (14).We take, Y T , Y C , Y B , as our infected compartments, then using the notation in [19], the Jacobian matrices F and V for the new infection terms and the remaining transfer terms are respectively given by, It follows that the basic reproduction number of the system ( 10)-( 14), denoted by R 0 , is given by where ρ is the spectral radius.Further, using Theorem 2 in [19], the following result is established.
The basic reproduction number (R 0 ) measures the average number of new infections generated by a single infected individual in a completely susceptible population [2,6,12,19].Thus, Lemma 1 implies that T. gondii can be eliminated from within the host (when R 0 < 1) if the initial sizes of the sub-populations are in the basin of attraction of the DFE, E 0 .
Proof.The proof is based on using a comparison theorem.The equations for the infected components in ( 10)-( 14) can be written in terms of where, M = X * − X (t), the matrices F and V are as given in Section 3.1, and Q is the non-negative matrix given by Using the fact that the eigenvalues of the matrix F − V all have negative real parts (see the local stability result given in Lemma 1, where ρ(F V −1 ) < 1 if R 0 < 1 which is equivalent to F − V having eigenvalues with negative real parts when R 0 < 1 [19]), it follows that the linearized differential inequality system ( 19) is stable whenever Thus, by comparison theorem [16,9], The above result shows that T. gondii will be eliminated from within the host if the threshold quantity R 0 can be brought to a value less than unity.

Endemic Equilibrium (EE)
Let ) be any arbitrary equilibrium of the model ( 10)-( 14).Conditions for the existence of equilibria for which T. gondii is endemic within the host (where at least one of the infected variables is non-zero) can be obtained as follows.Let, be the associated force of infection.To determine the existence of the endemic equilibrium, we consider first the case where immunity is not present (i.e Z = 0).Setting the right-hand sides of the model to zero gives the following expressions (in terms of λ * 1 and λ * 2 at steady state): Substituting the expression in (21) into the expression in (20) gives It follows that x > 0 if and only if R 0 > 1.This result is summarized below: Theorem 3 The model ( 10)-( 14) with Z = 0 has a unique endemic equilibrium whenever R 0 > 1.
Next we consider the case where the immune system responds throughout infection period (i.e Z = 0).To determine the existence of the endemic equilibrium, setting the right-hand sides of the model to zero gives the following expressions: Substituting the expression in ( 22) into the expression in (20) gives It follows that x > 0 if and only if R 0 > 1 and ρ > δ.This result is summarized below: Theorem 4 The model ( 10)-( 14) with immune response (Z(0) > 0) has a unique endemic equilibrium whenever R 0 > 1 and ρ > δ.
Thus, to obtain a unique endemic equilibrium, Theorem 4 implies that in the presence of immune response, ρ, the maximum attack rate, must be greater than δ, the removal rate of the immune system response for this case to happen.

Discussion
Recall from (17) that the reproduction number R 0 is proportional to X 0 , the total number of host cells before invasion of the parasites.Assuming the invasion and reproduction kinetics of the parasites are the same in different organs, larger organs intend to have more cells and thus larger R 0 values, which will make them more suitable for the parasites to dwell.This is probably why T. gondii are most frequently found in brain, heart, and muscle in a host [8].To investigate the influences of various parameters to the diseased state, we will rewrite the endemic equilibria in (21) and ( 22) more explicitly without using the terms λ 1 and λ 2 .
The equilibrium without immune response in (21) can be rewritten as follows: It follows from (23) that, in endemic state, the number of healthy cells, X * , is less than the original number of cells, X 0 .The model predicts that the steady state of the disease reaches a dynamic balance between three different stages of the parasites: T is positively correlated to X 0 .Thus, steady state parasite loads are related to the size of the organ.Also, recall the reproduction rate of uninfected cells is λ = d X 0 .Assume the reproduction rate λ is a constant, it can be seen from ( 23) that, an organ with longer life expectance 1/d will have larger parasite loads.Since brain cells are permanent, this also explains why T. gondii is mostly like to be found in brain [8].
Similarly, the equilibrium with immune response in ( 22) can be rewritten as follows: It is interesting to note that the relationships between the three parasite loads, Y * T , Y * C , and Y * B , remain the same as those in the absence of immune response although the steady state value of Y * T is now determined by kinetics of immune response.We also note the parasite loads Y * T and Y * C do not depend on the original number of host cells X 0 whereas the load Y * B is proportional to X 0 .Moreover, when the reproduction rate of uninfected cells λ = d X 0 is kept at a constant, a larger life expectance 1/d will lead to a larger value of X * and thus a higher parasite load Y * B .Again, the analysis shows that T. gondii favors cells with long life expectance such as the brain.

Numerical Simulations
In order to investigate the effects of T. gondii infection, we first introduce model assumptions and estimate parameters of infection dynamics using experimental data available in the literature.Numerical simulations here use a mouse spleen as an example.We estimate a healthy spleen has X 0 = 10 8 cells.Assume that the life expectancy of spleen cells is 1 month, which leads to the death rate as d = 1.389×10 −3 h −1 .In the current model, we assume at most one parasitophorous vacuole (PV) can form within a host cell.We further assume parasites within the same PV are in the same stage and replicate simultaneously.Experimental data in Weiss and Kim [20] indicate that the doubling time of tachyzoites is 6 hours and that of bradyzoites is 24 hours.Tachyzoites in vivo often lyse the host cell after reproducing 2 or 3 times [4]; thus, we choose a T = 1/(18h) and k T = 8/(18h).Cysts of bradyzoite may contain more than 1000 parasites [7].We assume encysted brayzoites burst after reproducing 10 times and thus estimate a C = 1/(240h) and k C = 1024/(240h).After a parasite is released from a PV into the organ, we assume the parasite may interact with 10 host cells and the probability of invasion of individual host is 2%.It follows that β T = 8.889 × 10 −10 h −1 and β B = 8.533 × 10 −9 h −1 .Since tachyzoites can convert to bradyzoites after about 20 generations of reproduction [17], we estimate r T = 1/ (108h).Weiss and Kim [20] showed that 48 hours after bradyzoites invade a tissue, tachyzoites start to appear; thus, we estimate r B = 1/ (48h).The current model does not consider detailed immune mechanisms.Instead, we consider the effector cells of the immune system acting on tachyzoites in PV.Assume the interaction rate between effector cells and tachyzoite PVs as c T = 1.67 × 10 −8 h −1 .Assume the degradation rate of the immune effector cells to be δ = 1/ (48h).Let h = 10 5 to account for the memory effect of immune response and let ρ = 10/ (24h) be the response rate.
Substituting the above parameters into Equation ( 17) yields R 0 = 30.63,which indicates that an infected cell will, on average, infect 30.63 uninfected host cells per hour.Consider an immuno-incompetent host, in which Z = 0, it follows that the equilibrium diseased state is X * = 3.26 × 10 6 , Y * T = 2.07 × 10 6 , Y * B = 6.16 × 10 6 , and Y * C = 4.61 × 10 6 .The total number of cells is 1.61×10 7 , which is reduced to 16% of the original number, X 0 .In contrast, in an immunocompetent host, the long term solution is X * = 9.30 × 10 7 , Y * T = 5.26 × 10 3 , Y * B = 4.46 × 10 5 , Y * C = 1.17 × 10 4 and Z * = 1.07 × 10 8 .The total number of cells is 9.34 × 10 7 , only slightly less than the number for disease free state.Transient responses in the absence of immune response are shown in Figure 3 while transient responses in the presence of immune response are shown in Figure 4.
We further use the model to study reactivation.While the number of parasites in an immunocompetent host is greatly suppressed, the parasites can be reactivated when the host becomes immunoincompetent [5].Here, we consider reactivation due to temporary impairment in the immune system.We assume a host is first infected with T. gondii and, after 3000 hours, the host's immune system is temporarily impaired so that the production rate ρ is reduced to 10% of the nominal value.We assume the impairment lasts for 200 hours and then the host's immunity recovers to the original level.This scenario would be analogous to a host suffering from immunodepression for a few weeks before overcoming the secondary infection and allowing for a full immune response to the toxoplasma.The results for this particular situation are shown in figure 5.It is clear that this temporary decrease in  immune response leads to a significant increase of parasite load within the host.

Conclusion
We have developed a mathematical framework to investigate intra-host dynamics of T. gondii.Assumptions and simplifications have been made about the biological processes including invasion, replication, and stage conversion.Parameters of the model have been estimated based on available experimental data.In the differential equation model we created, the effects of spreading parasites are examined.In the analysis, we first found the fixed points and analyzed the rate of infection R 0 .The first fixed point in the analytical model was the parasite-free equilibrium point.The second and third equilibrium points are identical in expressions for Y T , Y C , and Y B .It is interesting to note that the second and third equilibrium points vary in the values for the immune response fixed point.The immune Figure 5: Reactivation of parasites due to temporary impairment of the immune system: variation of healthy host cells (left) and variation of infected host cells (right).The infected host cells include those infected with tachyzoites (solid), cysted bradyzoites (dashed), and early-stage bradyzoites (dash dotted).regulated equilibrium only exists when ρ > δ.
It is also important to note the assumption made in analyzing this model.We have assumed the invasion dynamics of the free parasites are much faster than the replication and stage conversion, which leads to quasi-static simplification of the free parasites.This assumption is valid because the bursting of cells with tachyzoites release free tachyzoites, this release was modeled as a contribution of infection directly to uninfected hosts.The use of the Holling's Type II functional response models the way an immune system should respond: Very limited response to low numbers of invaders and a quickly growing response up to a threshold as the number of invaders increases.This type of response allows us to simplify the immune system into a mathematical model.
The critical value for R 0 found in this model indicates that our infection rate is dependent upon the initial number of uninfected hosts, the death rate of cells infected with bradyzoites and tachyzoites, the invasion rate of tachyzoites and bradyzoites contained within a host, and the conversion rate from tachyzoites to bradyzoites.Analyses of the reproduction number and the endemic solutions indicate that T. gondii favors large organs with long life expectance.This agrees with the experimental observation that T. gondii are commonly found in skeletal muscles, brain, and myocardium.
Numerical simulations show that the immune system plays a pivotal role in suppressing the growth of tachyzoites within host cells.This suppression, once the system reaches endemic behavior, allows for the body to exit the acute infection stage and begin the long-term, virtually symptom-free, state.Without immune response, the tachyzoites would be free to replicate and invade many different hosts.
Tachyzoites are rapidly dividing and responsible for the acute infection whereas the slowly replicating bradyzoites are located within tissue cysts, which protect the parasite from the host immune system and make it inaccessible to drugs [8].The differentiation of tachyzoites into bradyzoites is a response to the onset of protective immunity whereas the dormant bradyzoites are able to reconvert into tachyzoites to cause fatal infection in patients.
Therefore, stage conversion between tachyzoites and bradyzoites plays a pivotal role in the pathogenesis, transmission, persistence, and reactivation of the disease.
Future work of this system will include a more detailed description of immune response.While using the Holling's Type may accurately model the conceptual framework of an immune response, more evidence is needed to qualify this technique as accurate.A further expansion of this model might include a spatial array in which the disease can propagate.While our model tracks the disease throughout the spleen and assumes a homogenous distribution of cells throughout, the parasites are actually capable of starting in the stomach and invading the brain, muscles, and liver.Therefore, a spatial model may be able to describe the complicated dynamics that describe how the parasites can move throughout the body.

Figure 1 :
Figure 1: Diagram representing the life cycle of T. gondii ; Figure adopted from [10].

Figure 2 :
Figure 2: Compartmental Model description of the life-cycle of Toxoplasma gondii.

Figure 3 :
Figure 3: Response in the absence of immune response: variation of healthy host cells (left) and variation of infected host cells (right).The infected host cells include those infected with tachyzoites (solid), cysted bradyzoites (dashed), and early-stage bradyzoites (dash dotted).

Figure 4 :
Figure 4: Response in the presence of immune response: variation of healthy host cells (left) and variation of infected host cells (right).The infected host cells include those infected with tachyzoites (solid), cysted bradyzoites (dashed), and early-stage bradyzoites (dash dotted).