A fractional diffusion model of CD8+T cells response to parasitic infection in the brain

Toxoplasma gondii (T. gondii) is a parasitic pathogen that causes serious brain diseases in fetuses and patients with immunodeficiency, particularly AIDS patients. In the field of immunology, a large number of studies have shown that effector CD8 + T cells respond to T. gondii infection in the brain tissue through controlling the proliferation of intracellular parasites and killing infected brain cells. These protective mechanisms do not occur without T cell movement and searching for infected cells, as a fundamental  feature of the immune system. Following infection with a pathogen in a tissue, in their search for infected cells, CD8 + T cells can perform different stochastic searches, including Levy and Brownian random walks. Statistical analysis of CD8 + T cells in response to infected brain cells could be described by a Levy random walk., In this work, by considering a Levy distribution for the displacements, we propose a space fractional-order diffusion equation for the T cell density in the infected brain tissue. Furthermore, we derive a mathematical model representing CD8 + T cell response to infected brain cells. By solving the model equations numerically, we perform a comparison between Levy and Brownian search strategies.


Introduction
After the discovery of the Toxoplasma gondii (T. gondii) parasite in 1907, it has been extensively studied in various fields, such as morphology, immunology, and identification of diseases caused by T. gondii parasite [48,80]. From a morphological point of view, there are different forms of T. gondii parasite that lead to a complex life cycle for T. gondii and widespread infection in almost all warmblooded vertebrate hosts [20]. Initially, T. gondii oocysts are created in the body of cats, as definitive hosts and then transferred to their feces, which can contaminate the environment, such as soil, water, and food. Subsequently, through the transmission of oocysts into the bodies of secondary hosts, other types of T. gondii parasites, known as tachyzoites and bradyzoites, are formed. In humans, T. gondii can be perilous for the fetus during pregnancy and cause congenital diseases, known as toxoplasmosis, such as hydrocephalus, microcephaly, and blindness [25]. T. gondii parasites can maintain their survival in the brain tissue through the slow proliferation of bradyzoites and the creation of tissue cysts in different types of brain cells [85]. In adults, following Toxoplasma infection in the brain, in addition to the possibility of developing psychiatric disorders, like schizophrenia [75], the transformation of tissue cysts into free tachyzoites can jeopardize the life of patients suffering from immunodeficiency, particularly AIDS patients [22,47]. This issue highlights the importance of the immune system in controlling the proliferation of tachyzoites, which are capable of rapid division within the brain cells [18].
The human body is equipped with a complicated and ingenious system, known as the immune system that is composed of molecules, such as cytokines and chemokines, and numerous cells, like B cells and T cells. The goal of the immune system is to combat pathogenic organisms, like bacteria, viruses, and parasites in the body leading to the control of the pathogen growth and the elimination of the infected cells from the host body. There are two main groups of immune responses: innate immunity and adaptive immunity. Innate immunity includes the first reactions of the immune system against pathogens that enter the body, like changes in the temperature of the host body and nonsepecific mechanisms of the innate immune cells, such as the secretion of cytokines that promote the proliferation of other immune cells. These initial mechanisms are not sufficient to clear an infection. However, adaptive immunity comprises specific mechanisms that lead to the neutralization of pathogens through antibodies produced by B cells and the destruction of infected cells by CD8 + T cells. Although innate and adaptive cells have different biological functions, they are closely linked to each other. For instance, in the case of T. gondii infection, initially, the innate immune cells, such as dendritic cells respond to T. gondii tachyzoites. One of these responses is to secrete molecules, like IL−12. Then, CD8 + T cells, as the adaptive immune cells, are stimulated by IL−12 leading to the secretion of other molecules, such as cytokine IFN−γ that can control the growth of intracellular parasites and help T cells to penetrate into the brain, as protective and regulatory mechanisms, respectively [23,78]. Thus, CD8 + T cell deficiency causes inadequate production of the cytokine IFN−γ, which triggers the multiplication of a large number of tachyzoites in the brain tissue. For more details on basic concepts in immunology, refer to the book [21] by E. Paul. Over the last decade, as a fundamental characteristic of the immune system, the movement of immune cells in the body and search patterns in different tissues, including lymph nodes [79], skin [2], brain [31], and lungs [54] have been considerably studied. As soon as pathogens enter the body, the immune cells need to move continuously and search for different types of cells, like infected cells in various tissues. For instance, upon infection, naïve T cells, i.e., nonactivated T cells search for antigen-presenting cells, especially dendritic cells (DCs) in the lymph nodes. Following interactions with DCs, they can convert into effector CD8 + T cells. Effector T cells then move toward infected tissues where they search for infected cells so that they can perform their protective mechanisms. Therefore, the movement of T cells and their search for pathogens during an infection in the body are essential for immune system activities. T cell movement could have two components: directional movement (chemotaxis) and stochastic motion. By chemotaxis, T cells are directed to the infected organs by chemical signals, such as chemokines [55,83]. For example, Landrith et al. [43] have examined how effector CD8 + T cells are recruited into the T. gondii-infected brain tissue by adhesion molecules, VCAM-1, VLA-4, and chemokines CXCL9/10. However, studies of T cell motility in tissues have proposed that similar to foraging animals, T cells do a random search for their targets [50,79]. Analysis of these random motions suggests two types of search strategy: Brownian motion and Lévy movement. In complex diffusion processes, occasional large jumps between many short displacements can be considered as a distinctive characteristic of a Lévy random walk (RW) versus a Brownian motion. In that case, for a Lévy RW, the probability distribution function (pdf) for the random movements is obtained by a heavy tail distribution with an asymptotic behaviour as ∼ l −(α+1) , where l is the length of the displacements and 0 < α < 2 [49]. Figure 1 shows three sample random displacements of a particle following a Lévy RW. In all cases, initially, the walker is placed at the origin (x, y) = (0, 0) and randomly makes 10 3 steps with different values of α. As the value of α increases, the probability of large displacements decreases so that with α = 2, a classical Brownian motion is obtained. For more details on the simulation of Lévy distributions refer to [13,81]. So the existence of large jumps for Lévy RWs enables T cells to travel large distances within the tissues to find their targets. Studies have shown that due to finding rare DCs and also the lack of information on the exact location of DCs in the lymph nodes, naïve T cells perform a random search that can be described by a Brownian RW [59,65]. However, the search pattern of effector CD8 + T cells in a peripheral tissue where the tissue is injured by pathogens could be represented by a Lévy RW or a Brownian motion, depending on factors, such as the type of infection and the features of infected tissues. For instance, Effector T cells in inflamed lungs can have a similar foraging pattern to naïve T cells [53]. However, Harris et al. [31] observed that CD8 + T cells perform a Lévy RW search strategy in the brain of T. gondii-infected mouse. The difference between T cell search patterns in lymph nods and peripheral tissues may result  from factors, such as the properties of tissue and intrinsic differences between nonactivated and effector CD8 + T cells [41].
In addition to CD8 + T cells in the T. gondii-infected brains, there are many species that perform a Lévy RW search for their targets. Some of these species are listed as follows: jackals [3], spider monkeys [9], honeybees [63,64], bumblebees [45], fruit flies [62], and large marine predators [33,69]. In the field of modelling a large number of random walkers, such as living organisms following a specific random motion in an environment, diffusion equations for the density of the random walkers can be obtained. These equations consist of terms with time and space derivatives. For a Brownian motion, a Gaussian probability distribution function (pdf) for the displacements leads to a second-order diffusion equation [56]. However, for a Lévy RW, a power-law asymptotic behaviour of the pdf for the displacements results in a space-fractional-order diffusion equation with order 1 < α < 2. In some anomalous diffusion processes, the random walkers tend to pause between displacements, and hence they diffuse slower than predicted by a Brownian motion. Therefore, the pdf for waiting times between jumps exhibits a powerlaw asymptotic behaviour. In that case, a time-fractional-order equation with order 0 < γ < 1 is obtained. However, in the case of a Brownian motion, considering a Poissonian pdf for the pausing times leads to a time integer-derivative with order γ = 1 [49]. In recent years, such modelling using fractional derivatives has been applied to a broad range of problems in finance [12,67], plasma turbulence [15,16], epidemiology [28,30], ecology, and biology [61,76]. As well as the applications of the fractional calculus, designing efficient numerical methods for solving the fractional models is one of the key issues that has been widely studied [4,5,35,36,51,66].
In immunology, mathematical models have been widely applied to study the dynamics of host immune responses to tumour cells and viral infections, like hepatitis, human immunodeficiency virus (HIV), and influenza in the form of an ordinary differential system of equations with integer or fractional orders (see for instance [19,42]). In these models, the interactions between viruses or infected cells and immune responses are considered as predator-prey dynamics, i.e., viruses and infected cells are the prey and CD8 + T cells are the predators. In recent years, researchers have tried to examine the interplay between viruses or tumour cells and immune system cells based on their spread or movements in the organs of the body [8,14]. The resulting models are in the form of a system of reaction and diffusion equations. In these models, it is assumed that T cells search for viruses or tumour cells based on a Brownian motion, and viruses and tumour cells proliferate logistically and spread in the injured tissues relying on a Brownian dispersion.
Observation of the Lévy RW pattern performed by CD8 + T cells in their search for the T. gondii infected brain cells has led to more studies on such a pattern of the movement in the body [10,44]. Therefore, following such studies, it would be essential to derive appropriate mathematical models based on the Lévy RW that can shed some light on differences with the classical Brownian motion, which has been also observed in the immune system. Such modelling based on the Lévy random motion leads to space-fractional-order diffusion equations. However, recently, in the field of immunology, fractionalorder models of infections, such as HIV have been just designed relying on the time-fractional operators [1,6,34]. In these models, CD4 + T cells, as the immune cells, are the target of HIV infection and the dynamics of free HIV particles, healthy and infected CD4 + T cells have been investigated by numerical and analytical methods. For instance, Baleanu et al. [6] performed a comparison of their results by applying ordinary, Caputo and Caputo-Fabrizio time-fractional operators to their models. However, Khalid K et al. [1] introduced a new fractional model by using the conformable fractional derivative and examined the analytical and numerical solutions of their model. However, in these models, the effect of the CD8 + T cells on HIV infection has not been considered. Hence, to our knowledge, applying spacefractional-order operators to modelling the immune system response to infection has been ignored. This is obviously a serious knowledge gap since the function of the immune system is necessarily dependent on the movement of T cells. Thus, the aim of this study is to address this knowledge gap by applying the Lévy foraging pattern of T cells to the spatio-temporal modelling of the interaction between CD8 + T cells and bradyzoite-infected cells during T. gondii infection in the brain. Here, the main idea will be to understand how the Lévy RW of CD8 + T cells can lead to an efficient immune response to the infected brain cells. Therefore, the fractional-order model results will be compared with the results obtained with the classical model in order to highlight the effect of Lévy RW on killing infected cells distributed broadly within the brain tissue.
We organize the remainder of the paper in the following structure. In the next section, we derive a space-fractional-order diffusion equation for the CD8 + T cells density and also a space-time model describing the interactions between CD8 + T cells, healthy and infected brain cells. We discuss the system parameters and the estimation of their values in Section 3. A numerical method for solving space-fractional-order diffusion equations is applied in Section 4. Our model is then used to estimate the efficiency of effector CD8 + T cells for a Lévy RW versus a Brownian motion in Section 5, and in the last section, we finish the paper with a review of the results and also some suggestions for future work.

Model description
In this paper, our goal is to investigate the efficacy of a Lévy search pattern performed by effector CD8 + T cells in the T. gondii-infected brain tissue. When effector T cells enter the brain, they perform random searches to find the infected brain cells . Therefore, we can study the impact of a Lévy random search on the density of infected cells distributed throughout the brain tissue. To do so, we first derive a model with non-local diffusion representing the density of CD8 + T cells, and then according to the T. gondii life cycle in the brain, we shall formulate a model composed of three coupled equations representing the density of healthy and infected brain cells in space and time with no diffusion term .

Modelling the superdiffusion of CD8
Like viral pathogens, as a result of T. gondii parasites entering the body, naïve T cells located in the lymph nodes are activated and proliferate due to interactions with dendritic cells. Activated or effector CD8 + T cells can then do their protective mechanisms in the infected tissues of the body, including the brain tissue [37].
It is worth noting that CD8 + T cells can infiltrate into the brain parenchyma (gray and white matter) only during infection [32]. In the pathway of CD8 + T cells toward the brain tissue, there are barriers, known as blood-brain barriers (BBBs) that control the penetration of the immune cells into the brain parenchyma [82]. However, CD8 + T cells can overcome such barriers through a series of steps leading to a directed migration of CD8 + T cells into the brain tissue. In [43], for T cell infiltration into the brain parenchyma during T. gondii infection, a process consisting of three steps has been considered. Initially, the speed of CD8 + T cells moving in the blood vessels is slowed down by the adhesion molecules. At the next step, they can infiltrate into the perivascular space, and eventually chemokines lead them to the brain parenchyma. When CD8 + T cells reach the brain, they do random searches for their targets, i.e., infected brain cells. The properties of such random searches have been studied by Harris et al. [31].
From a statistical point of view, there are basic properties for Brownian and Lévy RWs that can be used to understand their differences. One of these properties is the mean-square displacement (MSD). If the MSD as a function of time, grows linearly, the random motion can be described by a Brownian RW, and the nonlinear growth of the MSD represents a Lévy RW. By considering this property, the random search strategy of living organisms in an environment can be determined. For instance, Harris et al. [31] examined the movement of CD8 + T cells in the brain of T. gondii-infected mouse and suggested that the MSD grows superlinearly, i.e., < x 2 >∼ t 1.4 . It means that a Lévy RW search pattern is adopted by CD8 + T cells in response to Toxoplasma infection in the brain. In that case, CD8 + T cells can diffuse within the brain tissue faster than expected from a Brownian motion leading to superdiffusion. Another property is the probability distribution function (pdf) for the jumps of the random walker. For a Brownian RW, the pdf is a Gaussian function. However, in the case of a Lévy RW, the pdf exhibits a power-law asymptotic behaviour. For instance, based on the statistical analysis of T cell movement in the infected brain [31], the pdf for the displacements of CD8 + T cell has a power-law asymptotic behaviour as ∼ |x| −(α+1) , where |x| is the length of displacements and α = 1.15. Thus, when CD8 + T cells enter the brain tissue from a source, with non-zero probability, they can make large displacements in the brain tissue. For more details on the properties of Lévy and Brownian RWs, refer to [39].
To derive a fractional-order diffusion equation representing the density of CD8 + T cells in the T. gondiiinfected brain, we shall consider the continuous-time random walk (CTRW) processes consisting of a large number of walkers (here T cells) that make random jumps, x 1 , x 2 , ..., x i , ... at times t 1 , t 2 , ..., t i , .... As the displacements x i and the waiting times τ i = t i − t i−1 are random variables, we define the probability distribution functions (pdf's) η(x) and ψ(τ ) for the jumps and waiting times, respectively. By considering a separable CTRW and given τ and ψ, we can get the probability of finding a walker P (x, t) at place x and time t by the Montroll-Weiss equation [52] in Fourier-Laplace form as follows: whereψ(s) andη(k) are the Laplace and Fourier transforms of the pdf's ψ and η, respectively defined as follows: Different CTRW processes can be obtained by different asymptotic behaviours of the jump and waiting time pdf's, which lead to whether integer-or fractional-order diffusion equations. For example, we suppose that the pdf's for the waiting times and jumps exhibit heavy-tailed (power-law) asymptotic behaviours as follows: where 0 < γ < 1 and 1 < α < 2. If we use the fact that where f (0) is the initial condition, −∞ D α x and x D α +∞ are the positive (left) and negative (right) spacefractional-order Riemann-Liouville operators, and c 0 D γ t is the time-fractional-order Caputo operator defined as: where Γ(.) denotes Euler's gamma function, by taking the inversion of the Laplace and Fourier transforms of Equation (1), one can then find the following space-time fractional-order diffusion equation: where coefficient K α,γ represents a generalised diffusivity whose dimension is [K α,β ] = m α s −γ and the parameter β ∈ [−1, 1] is a skewness parameter that shows a preferred direction of displacements that can be seen in heterogeneous systems. When β = 0, the distribution is symmetric and space derivative represents a symmetric Riesz derivative. For more details on the derivation of Equation (3), refer to [11,49]. Here, we note that by defining the concentration of walkers C(x, t) in one dimension as the number of walkers per unit length at position x and time t and assuming that there are N random walkers in total with the same probability P (x, t), we can easily see that C(x, t) = N P (x, t) also satisfies Equation (3).
As an application of Equation (3) to obtain the concentration of walkers following a CTRW defined by (2), by considering a large number of CD8 + T cells following a superdiffusion in the T. gondii-infected brain tissue and the fact that the pdf for the T cell displacements has a power-law asymptotic behaviour with order α, where α = 1.15, in one dimension, the concentration (density) of CD8 + T cells E(x, t) at position x and time t with the dimensions of the number of T cells per unit length is the solution of the following space-fractional-order diffusion equation: where L is the domain length, the parameter K α is a fractional diffusion coefficient whose dimension is [K α ] = m α s −1 , and 0 D α x and x D α L are the left and right space-fractional Riemann-Liouville derivatives on [0, L], respectively defined as follows: where n = 1 + [α] such that [α] = max{m ∈ Z|m ≤ α} and Z is the set of integers. Here, we get n = 2 as in this study α = 1.15. Since the Riemann-Liouville derivatives are singular on the boundaries of the domain, instead, Caputo derivatives can be used. The left and right Caputo derivatives of order α, respectively are defined as follows: For more details on the fractional derivatives, see [57,58].
It should be noted that in addition to the power-law asymptotic behaviour of the pdf for the CD8 + T cell jumps in the brain, based on the same experiments in [31], the pdf for the waiting times exhibits a power-law asymptotic behaviour as ∼ τ −(γ+1) , where γ = 0.7. Here, for the sake of simplicity, we have only considered the space-fractional-order operator in our model.

cells response to infection
Just as CD8 + T cells cross BBBs to enter the brain, T. gondii parasites can go through such barriers by certain mechanisms [40] and then reach the brain parenchyma. When T. gondii parasites (tachyzoites) enter the brain, they begin to infect brain cells, such as astrocytes, neurons, and microglia. The infection occurs as soon as free parasites attack brain cells and multiply inside them. Tachyzoites are then placed in a covering called parasitophorous vacuole (PV). This stage of infection is called the early stage of the T. gondii life cycle. In order to maintain their survival, tachyzoites need to enter into the second stage of their life cycle. Therefore, there is another form of T. gondii parasites, known as bradyzoites. Similar to tachyzoites, they multiply in a PV. But, the growth and proliferation of bradyzoites are slower than tachyzoites. PVs containing bradyzoites create intracellular tissue cysts over time that can be stable and persistent in the brain tissue for a long period of time [70].
As a result of parasites entering the brain tissue, effector CD8 + T cells move from lymph nodes to the brain parenchyma. They then exert two protective responses against infected brain cells through (IFN)-γ and perforin molecules. The type of such mechanisms depends on the type of T. gondii parasite. It means that in the case of infection with tachyzoites, effector CD8 + T cells are able to stop the proliferation of intracellular tachyzoites through the secretion of (IFN)-γ [17]. It is important to note that neurons provide favourable conditions for the creation and stability of tissue cysts in the brain tissue because neurons live longer than other brain cells, and according to in vitro studies, (IFN)-γ mechanism has little effect on the infected neurons compared to other infected brain cells [68]. Effector CD8 + T cells can also identify all bradyzoite-infected brain cells and also those infected cells that carry tissue cysts of various sizes. Following the identification, CD8 + T cells bind themselves to the surface of infected cells, enter the tissue cysts, and eventually, through the secretion of perforin, tissue cysts are eliminated from the brain tissue [73,74].
Based on the above biological information, to derive a mathematical model representing effector CD8 + T cells response against infected cells, we shall consider the density of four cell populations at position x and time t, the susceptibles (healthy brain cells) S(x, t), tachyzoite-infected brain cells I 1 (x, t), bradyzoite-infected brain cells I 2 (x, t) and effector CD8 + T cells E(x, t). We assume that susceptibles convert into tachyzoite-infected cells at a rate θ and that the transmission from tachyzoite-to bradyzoiteinfected cells occurs at a rate β. In this model, we only consider the perforin-mediated response against infected cells. Hence, it is assumed that bradyzoite-infected cells could be killed at a rate p by effector CD8 + T cells. The parameters d 1 , d 2 , and d 3 represent the mortality rate of infected cells with tachyzoite, infected cells with bradyzoite, and CD8 + T cells, respectively. Since effector CD8 + T cells present in the infected brain are recruited from lymph nodes toward the brain tissue, the model considers the function F (x, t) representing the source of CD8 + T cells. This model is schematically represented in Figure 2. Figure 2. Schematic representation of transmission from susceptibles S to infected brain cells with tachyzoite I1 at a rate θ, tachyzoite-infected cells conversion into bradyzoite-infected brain cells I2 at a rate β and effector CD8 + T cells E response to I2 at a rate p. The function F represents the source of effector CD8 + T cells. The parameters d1, d2, and d3 are the mortality rates of S, I1 and I2, respectively.
By taking the schematic representation of the model and Equation (4) into account, we obtain the following model representing the density of healthy, T. gondii-infected brain cells and CD8 + T cells in the brain tissue: If it is assumed that effector CD8 + T cells perform a Brownian RW search for infected brain cells, then the model reduces to a classical model with α = 2. Therefore, the advantage of such a model is that we can compare the results of a Lévy RW (α < 2) with those obtained by the Brownian RW assumption.

Model parametrization
In this section, we shall discuss the value of the parameters defined in Equations (9)- (12). Since there have been very little attempts to model T. gondii infection [71,72] and necessary experimental data are scant, here we try to estimate the value of the model parameters from direct measurements or by considering the parameter value in other studies.
One of the problems in relation to brain infection with T. gondii is the infection rates of different brain cells from different host species. Extensive studies have been conducted on the infection rate of brain cells by free tachyzoites (see for instance [24,46]). When we investigate the results of the experiments, we see that there are many differences in the value of rates that are mainly due to factors, such as the type of brain cells, the brain regions where the cells come from, the host species (rat, human, and mouse), and the type of Toxoplasma parasite. For instance, the rate of infection by the first type of parasites (RH strains) is higher than the second type (ME49 and Prugniaud). Neurons of the hippocampus region of the brain from rats are infected by T. gondii at a lower rate than astrocytes, while it is more likely that cortical neurons from rats are infected. According to studies in vitro [24,46], if we consider a range of infection rate 5% to 30% based on the type of brain cells, one day after infection by adding 10 5 tachyzoites to the cultures and assume that each infected cell egresses 8 to 16 tachyzoites, then we estimate a range of the infection rate given by θ = 4 × 10 −6 day −1 cells −1 to θ = 5 × 10 −5 day −1 cells −1 . For instance, if we consider the infection rate 5% by 10 5 parasites and also assume that each infected cell egresses 8 tachyzoites, we can then compute θ = 0.05 × 8 × 10 −5 = 4 × 10 −6 day −1 cells −1 . In order to estimate the transmission rate β and the death rates of infected cells d 1 and d 2 , we refer to the experiments in [60]. When free tachyzoites attack brain cells, they multiply inside infected cells by dividing. After 5 divisions, they go out of the infected cells, which leads to infecting the brain cells in their neighbourhood. Tachyzoites have a half-life of 6 hours. Thus, the estimated daily death rate of tachyzoite-infected cells is given by d 1 = ln 2/(5 × 6) h −1 ≈ 5 × 10 −1 day −1 . These tachyzoites convert into bradyzoites approximately after 20 divisions, which results in the estimated transmission rate β = 1/(20 × 6) h −1 = 2 × 10 −1 day −1 . Bradyzoites have a half-life of 15 hours. If we suppose that bradyzoites go out of the infected cells after 10 divisions, the estimated daily death rate of bradyzoiteinfected cells is given by d 2 = ln 2/(10 × 15) h −1 ≈ 1 × 10 −1 day −1 . The average per capita killing rate p is the number of targets (bradyzoite-infected cells at the chronic stage) destroyed per effector CD8 + T cell per day. There are very few studies discussing perforin-mediated response against bradyzoite-infected cells or tissue cysts at the chronic stage of infection [73,74]. Hence, experimental data are scant. However, based on the experiments during 7 days in [73], if we assume that on average 50 brain cysts per day are killed by immune CD8 + T cells and also that there are about 5 × 10 4 T cells per day, then by the definition of p, we can measure directly an average value that is given by p = 10 −3 day −1 cells −1 . Although this value may be better estimated in future works, this estimated value is sufficient for our numerical simulations because in this work, we aim to perform a comparison between the results of Brownian and Lévy RWs with the same conditions. In order to obtain an approximation of diffusivity K α , we follow the observations of Harris et al. [31]. To do so, we know that K α has SI units of m α s −1 . Thus, the mean displacement length of diffusion L can be defined as L = (K α T ) 1/α , where T is the corresponding time interval. Hence, based on the findings of Harris et al. [31] on the diffusion coefficient of CD8 + T cells in the T. gondii-infected brain, if we assume that L = 8 µm and T = 1 min, the fractional diffusion coefficient can then be estimated as K α = (8×10 −6 ) α /60 m α s −1 . It is important to note that this value could be reduced when no chemokine signals and CXCL10 are present in the infected brain. In that case, CD8 + T cells nevertheless maintain the behaviour of the Lévy search pattern.

Numerical solution of space-fractional order diffusion equations
Over the past decade, in order to solve time-space fractional-order diffusion equations, analytical methods are replaced by various numerical methods due to the complexities of fractional derivatives. Most of these numerical methods are based on the finite element (FE), finite-difference (FD), and Chebyshev pseudo-spectral schemes. In this work, we have considered the same FE method used by Hanert [26] to solve fractional-order transport models. More precisely, we have used a piecewise linear finite-element method based on a Galerkin formulation for solving Equation (12). In this method, the exact solution (unknown variable) E(x, t) is expressed as a sum of the unknown coefficients e j (t) and basis functions φ j (x) as follows: By introducing a partition of the domain [0, L] into N − 1 subintervals [x j , x j+1 ] with a constant length h, i.e., x 1 = 0, x N = L, and x j+1 − x j = h for j = 1, . . . , N − 1, we can consider the piecewise linear basis functions φ j (x) for j = 1, 2, . . . , N as follows: for j = 2, ..., N − 1, and for j = N , Now we approximate the source term F (x, t) by the same discretization as the model solution as follows: Since the basis functions φ j (x) for j = 1, 2, ..., N have a value of 1 at node x j and 0 at other nodes, f j (t) = F (x j , t). By using a Galerkin formulation, we replaceẼ(x, t) andF (x, t) in the model equation and then by orthogonalizing the discrete equation with respect to all φ j , we get the following equation: for i, j = 1, ..., N and By performing integration by parts (see [26] for details), we get a semi-discrete equation in a matrix form as follows: where e(t) = [e 1 (t) . . . e N (t)] T is the vector of unknown coefficients at time t and the elements of the T are the known values of the source term at a node j at time t. The elements of the mass matrix M and the diffusion matrix D read:

is a combination of the left and right
Riemann-Liouville derivatives of order α − 1 of φ j . Finally, we can discretize Equation (13) by using standard time integration methods. Here, a third-order Adams-Bashforth scheme has been considered. To solve Equation (12), we also need to consider initial and boundary conditions. Here, we present an example to illustrate how different values of the parameter α affect the behaviour of the solution. We have solved Equation (12) with F (x, t) = 0, and d 3 = 0. Figure 3, black curve shows a Gaussian function, as the initial condition. In the case of a Brownian motion (α = 2), the solution tails decay exponentially (Figure 3, red curve). But, for a Lévy RW (α < 2), the solution has a power-law asymptotic behaviour (Figure 3, blue curve). Since with a Lévy RW, particles can make large displacements, as we see, the solution decays more slowly compared to the Brownian motion assumption.

Effect of Lévy RW of CD8 + + + T cells on the infected cells
In this section, by using the model Equations (9)- (12) with the parameter values estimated in Section 3, we investigate the influence of a Lévy RW search performed by T cells on the density of infected cells in the brain. To do so, we can consider a comparison between the classical model based on the Brownian RW assumption (α = 2) and the fractional model relying on the Lévy RW assumption (α = 1.15).
In order to solve the proposed model, we need to define a computational domain L as a part of the brain tissue and the source function F (x, t) in Equation (12). The Brain tissue is composed of two kinds of tissue: grey matter and white matter. Based on the microscopic studies, tissue cysts can be formed in all regions of the brain tissue, but most of them are distributed in the cerebral cortex and hippocampus regions [7]. The cerebral cortex is the most outer layer of grey matter that covers the brain. Hence, we assume that the computational domain is a part of the cerebral cortex. A detailed description of the cerebral cortex is provided in [77]. We define the source term by F (x, t) := f (t)δ(x − L/2), where the time-dependent function f (t) gives the number of CD8 + T cells entering the domain in the middle per unit of time at time t and δ(x − L/2) is the Dirac delta function centred in the middle of the domain. We define the input function f (t) as a decreasing function given by f (t) := rN e −rt , where N and r represent the number of CD8 + T cells and the rate at which they enter the domain, respectively. To estimate the input rate r, we followed the observations of John et al. [37] who studied the kinetics of CD8 + T cells in the T. gondii infected brains. If we assume that the death rate of CD8 + T cells is equal to 10 −1 day −1 and consider a maximum number of CD8 + T cells in the brain by day 14, we can then estimate the value of the input rate to be given by r = 8 × 10 −2 day −1 . Here, it is assumed that CD8 + T cells enter the brain with a delay on day 3 [37]. It should be noted that in the numerical simulations, we consider the values of N according to the estimations of Harris et al. [31] (see supplementary information in [31]).
In order to compare the efficiency of CD8 + T cells for Lévy RW versus Brownian RW, we define the average number of bradyzoite-infected cells n(t) during simulation time as follows: where I 2 (x, t) is the density of bradyzoite-infected cells (tissue cysts) with units of the number of cells per mm. We discretize Equations (9)-(12) with a piecewise-linear FE method based on a Galerkin formulation. For more details on the discretization of the model equations, see Section 4 and appendix B. Equation (12) needs to be solved with appropriate boundary conditions that guarantee mass conservation in the numerical simulations (see appendix A for details).
In the first numerical simulation, we consider constant initial conditions S(x, 0) = 2.5 × 10 4 healthy brain cells and I 1 (x, 0) = 150 tachyzoite-infected cells per mm over the whole domain. We compare the average number of infected cells n(t) for Lévy and Brownian RWs of CD8 + T cells over the domain [0, L] of the brain tissue. Here, we suppose that L = 14 mm. According to our model, infected cells with tachyzoites I 1 and bradyzoites I 2 will first increase in a period [0, T 0 ] and then tend to 0 after T 0 . The time evolution of the average number of infected cells n(t) over the domain [0, 14] with different foraging patterns of CD8 + T cells is shown in Figure 4. Figure 4 shows that in the case of a Lévy RW, n(t) reaches its maximum value at time t = 8 and after that tends to 0 at time t ≈ 35. However, in the case of a Brownian motion, n(t) grows quickly and reaches its maximum value at time t = 11 and then decreases such that n(50) = 60. Here, we define the average of n(t) in a period [0, T ] byN = T 0 n(t)dt/T . That quantity is computed for Brownian RW (N B ) and Lévy RW (N L ). By taking the ratio between the values ofN B andN L , we getN L = 0.26N B in the period [0,30]. This result shows that a Lévy foraging pattern of CD8 + T cells is about four times more efficient than a Brownian motion in detecting and killing infected cells distributed throughout the brain tissue. To illustrate this result, we consider the diffusion profiles for CD8 + T cells on the domain [0, 14] with Lévy and Brownian RWs, which are shown in Figure 5.  30] for Lévy RW (NL) is equal to 26% of that value for Brownian RW (NB), highlighting the fact that a Lévy RW search strategy is more efficient than a Brownian motion in killing infected cells distributed in a large distance of the brain tissue. Figure 5 shows that the difference between Brownian and Lévy CD8 + T cell density models is most obvious in the tails of T cell distribution. The Lévy RW search strategy increases the number of T cells further away from the source more than a Brownian motion and thus leads to an increase in the thickness of T cell distribution tails. However, for a Brownian motion, CD8 + T cells are present more around the source, and hence they are not able to successfully eradicate infected cells located further away from the source (see Figure 5A and C). T cells with a Lévy search strategy can travel larger distances and hence scan the entire tissue much more rapidly (see Figure 5B and D). The objective of CD8 + T cell search is to remove all or nearly all infected cells from infected tissues. Since a Lévy RW search allows CD8 + T cells to scan further distances from the source more efficiently than a Brownian motion, our results suggest that a Lévy RW search can lead to a complete detection of T. gondii-infected cells in the brain tissue. For instance, the average number of infected cells n(t) over the domain [0, 14] with a Lévy RW vanishes after 35 days, while with the Brownian motion assumption, n(50) = 60.
We here note that from a biological point of view, when initially the infection load (burden) in the brain tissue is growing, activated CD8 + T cells migrate from lymph nodes to the brain tissue in order to eliminate the infected brain cells. In that case, the number of immune cells throughout the tissue increases so that the infection load decreases. At this stage, the spatial propagation of CD8 + T cells throughout the tissue in the case of Brownian and Lévy random walks are shown in Figure 5A and B, respectively. Once the infection is resolved, the number of CD8 + T cells declines. This is called the "contraction phase". Here, the spatial decline of T cells in the case of Brownian and Lévy random walks are compared in Figure 5C and D, respectively, and finally, T cells could persist at a level for long periods of time in the tissue. Figure 4 shows the time evolution of the infection load during the immune response to the infected brain cells.
In order to highlight the efficiency of Lévy RW versus Brownian RW on the domain [0, 14], we can also compute a time such that the average number of infected cells at that time is equal to 1% of the maximum value of n(t) for α = 2 (t B ) and α = 1.15 (t L ). By taking the ratio between the values of t B and t L for the domain [0, 14], we get t L = 0.58t B . However, that ratio for the domain [0, 6] is approximately 1. Here, if we increase the length of the domain from 6 mm to 20 mm, then t L and t B increase from 24 to 25 days and from 24 to 50 days, respectively. In other words, the ratio t L /t B decreases from 1 to 0.5 (see Figure 6). The numerical result shown in Figure 6 has been obtained by considering the time dynamics of the infection load throughout the tissue with different lengths. As the domain of infection gets larger, the average number of infected cells n(t) after reaching its maximum value with a Brownian motion decreases more slowly compared to a Lévy RW search. Therefore, this result highlights the importance of a Lévy RW search pattern for accelerating the elimination of infected cells distributed broadly within the brain tissue. The final example illustrates the differences between the elimination dynamics of infected cells with Lévy and Brownian RWs. Here, we consider the following model for the density of CD8 + T cells and bradyzoite-infected cells: with the initial condition I 2 (x, 0) = I 0 . We suppose that an equal number of CD8 + T cells per    14] with Lévy RW (α = 1.15) and Brownian RW (α = 2) obtained by solving Equations (15) and (16). The equal number of CD8 + T cells enter the right-hand side of the domain with a Lévy RW and the left-hand side with a Brownian RW, leading to a rapid elimination of infected cells over the whole domain on the right with a Lévy RW (blue curves), and a slower decline of infected cells faraway from the source with a Brownian RW (red curves). The arrows show the direction of the model solutions over time beginning from the initial condition of infected cells displayed by a black horizontal line. Here, the simulation has a duration of 5 days and the time interval between profiles is set to 1 day. and then observed a significant decrease in the number of tissue cysts. They suggested that a local expansion of CD8 + T cells in the brain tissue could lead to the high destruction of tissue cysts. In our model, we considered the impact of the protective mechanism of CD8 + T cells by perforin on the bradyzoite-infected cells. Our results suggest that the Lévy RW search enables T cells to accelerate the the elimination of tissue cysts. Thus, we predict that in addition to possible reasons for high reduction of cysts, such as a local proliferation of CD8 + T cells, the Lévy RW search strategy performed by CD8 + T cells can be a possible factor for marked decreases in tissue cyst numbers in the T. gondii-infected brain.

Conclusion
In their search for targets, many species in almost any biological environment do random searches when their intellectual capacities are constrained by the distribution of their targets. For instance, T cells can perform different search patterns for their targets across different tissues. Mathematical models are useful tools to study how the search strategies performed by T cells can lead to an efficient immune response against pathogens. Among all foraging patterns in tissues, random motion of CD8 + T cells in the T. gondii-infected brain is considerable due to a switch in the random search strategy of T cells from a Brownian motion observed in other tissues to a Lévy RW search in the brain tissue. Hence, studying the effective role of a Lévy RW search indeed requires specific tools. One of them is to use fractional diffusion equations to obtain the density of T cells in the infected brain.
In this work, we have derived a space fractional-order equation consisting of fractional-order diffusion and reaction terms to get the density of T cells in the infected brain. The fractional diffusion with order α < 2 represents a complex dispersion of T cells called superdiffusion, and the reaction terms describe the recruitment of T cells from lymph nodes into the brain and their natural death. When α = 2, the resulting equation represents the Brownian motion. Hence, we were able to perform a comparison between the results of Lévy and Brownian RWs. Based on the T. gondii life cycle in the brain, we have also derived space-time reaction equations representing the density of tachyzoite-and bradyzoite-infected cells in the brain tissue. This allowed us to study the effect of the Lévy RW search of CD8 + T cells on the bradyzoite-infected cell (tissue cysts) density by the perforin-mediated immune response. Our simulation results have shown that the average number of infected cells (N ) in the period [0, 30] with a Lévy RW on the domain [0, 14] is 26% ofN with the Brownian motion. This means that CD8 + T cells with a Lévy RW search are more successful than a Brownian motion in response to infection because T cells can travel faraway from the source, leading to more eradication of the cysts in the brain. Furthermore, as T. gondii parasites are capable of infecting large areas of the brain tissue, a Lévy RW search allows CD8 + T cells to scan large distances, which results in a complete detection of infected cells and hence all or nearly all infected cells are eliminated from the brain tissue. However, with the Brownian motion assumption, infected cells located further away from the source could be out of reach of CD8 + T cells.
We here note that Harris et al. [31] have also shown the efficiency of effector CD8 + T cells for Lévy RW versus Brownian RW. In their model, the efficiency has been defined based on the sum of T cells displacements at the moment when the first T cell reaches a target at the origin of a sphere. They showed that in the case of a Lévy RW, the sum of the displacements to reach a target is less than a Brownian motion, and hence a Lévy RW search is more efficient. However, our model provides a much more realistic picture of the efficient role of the Lévy search strategy on the protective mechanism of CD8 + T cells.
It is noteworthy to point out that the model could be further improved by considering some issues that have been ignored. One of them is the effect of the anti-parasitic mechanism of effector CD8 + T cells by IFN-γ on controlling the rapid division of tachyzoites. In [84], the two protective mechanisms of CD8 + T cells by perforin and IFN-γ are referred to as the lytic and nonlytic responses of the immune system against pathogens, respectively. This model may therefore be improved by considering both the lytic and non-lytic immune responses using the same approach that Wodarz et al. [84] applied to develop viral infection models.
Another issue is the study of non-Markovian property in anomalous dispersions. According to observations of Harris et al. [31], when CD8 + T cells make a displacement, the probability distribution function for the waiting or pausing times exhibits a power-law asymptotic behaviour as ∼ τ −(γ+1) with γ = 0.7. In that case, the anomalous diffusion of CD8 + T cells has a non-Markovian effect. By considering a large number of CD8 + T cells, for the T cell density in the brain, a time fractional-order diffusion equation with order 0 < γ < 1 is obtained. However, for a Brownian motion, the time derivative order is equal to 1 and the diffusion has a Markovian property. For solving space-time fractional-order equations, there are some efficient numerical methods based on Chebyshev pseudo spectral method. For instance, see [27,29].

Appendix A. Fractional Neumann boundary conditions
One of the delicate issues in the field of modelling using fractional diffusion equations is to apply appropriate boundary conditions to the models on bounded domains. One of the boundary conditions is Neumann (no-flux) that contains the first-order derivatives with respect to space on the boundaries of the domain. They are also known as classical Neumann boundary conditions. In this study, since our goal is to perform a comparison between the results of Lévy and Brownian RWs, the total number of CD8 + T cells entering the domain at time t should be the same in both cases. In other words, the issue of mass conservation should be guaranteed. In the case of Brownian RW, it can be done by using classical Neumann boundary conditions. However, they do not guarantee mass conservation for diffusion equations related to a Lévy RW, i.e., space-fractional-order diffusion equations. Thus, we have used the following fractional Neumann boundary conditions proposed by Kelly et al. [38]: for all t ≥ 0 at positions x = 0 and x = L and the fractional derivatives of order α − 1 are in the Caputo sense defined in Equations (7) and (8). In the numerical examples, the following boundary conditions in the symmetric case, i.e., when β = 0 have been considered: for all t ≥ 0 at positions x = 0 and x = L.

Appendix B. Numerical discretization of the model
In this section, we shall illustrate the discretization of the model equations by finite-element scheme. To do this, we approximate the exact solutions of the model equations by the following expansions in terms of basis functions φ j (x): where s j , a j , c j and e j are unknown nodal values. With the same method in Section 4, the resulting discrete Equations (9)-(11) then read: such that for any x ∈ J j , we have φ i (x)φ j (x) = 0. Now, for these values of j, we determine the possible values of index k such that φ i (x)φ j (x)φ k (x) = 0 for any x ∈ J j . Here, we have three cases. (1) for j = i − 1, by considering the two possible values of index k = i − 1, i, we have φ i (x)φ j (x)φ k (x) = 0 for any x ∈ J i−1 . (2) for j = i, by considering the three possible values of index k = i − 1, i, i + 1, we have φ i (x)φ j (x)φ k (x) = 0 for any x ∈ J i , and (3) for j = i + 1, the only possible values of index k are i and i + 1 such that for any x ∈ J i+1 , we get φ i (x)φ j (x)φ k (x) = 0. Therefore, we get the following nonzero values of the coefficients R ijk : For i = 1, the coefficients R 1,1,1 , R 1,1,2 , R 1,2,1 and R 1,2,2 , and finally for i = N , the coefficients