Transient Superdiffusion and Long-Range Correlations in the Motility Patterns of Trypanosomatid Flagellate Protozoa

We report on a diffusive analysis of the motion of flagellate protozoa species. These parasites are the etiological agents of neglected tropical diseases: leishmaniasis caused by Leishmania amazonensis and Leishmania braziliensis, African sleeping sickness caused by Trypanosoma brucei, and Chagas disease caused by Trypanosoma cruzi. By tracking the positions of these parasites and evaluating the variance related to the radial positions, we find that their motions are characterized by a short-time transient superdiffusive behavior. Also, the probability distributions of the radial positions are self-similar and can be approximated by a stretched Gaussian distribution. We further investigate the probability distributions of the radial velocities of individual trajectories. Among several candidates, we find that the generalized gamma distribution shows a good agreement with these distributions. The velocity time series have long-range correlations, displaying a strong persistent behavior (Hurst exponents close to one). The prevalence of “universal” patterns across all analyzed species indicates that similar mechanisms may be ruling the motion of these parasites, despite their differences in morphological traits. In addition, further analysis of these patterns could become a useful tool for investigating the activity of new candidate drugs against these and others neglected tropical diseases.


Introduction
Kinetoplastids protists are responsible for numerous diseases in humans and animals. Many of these protozoa are the etiological agents of neglected tropical diseases [1]. These diseases affect the lives of approximately one billion people around the world [2] and are considered a serious public health problem in several countries. The main regions affected are developing countries located in tropical areas, where the parasites have appropriate natural conditions for life-cycle and insects vectors are abound [3]. Leishmaniasis, Chagas disease, and African sleeping sickness are examples of neglected tropical diseases caused by kinetoplastids parasites-Leishmania spp, Trypanosoma cruzi, and Trypanosoma brucei, respectively [4]. These protozoa have a complex life-cycle, alternating between invertebrate (vector) and vertebrate hosts. These parasites have a flagellum at least during one of the evolutionary forms of their life-cycle [5][6][7]. The flagellum is a multifunctional organelle, responsible for cell propulsion and associated with control cell morphogenesis, chemotaxis, and cytokinesis process during last stage of the cell division cycle [8]. The parasites motility is a key to host-cell attachment invasion and colonization of host tissues [8][9][10][11][12][13].
Diffusive motion is ubiquitous in nature and plays a fundamental role in the motility of swimming microorganisms [14][15][16][17][18][19][20]. Researchers in statistical mechanics have focused on phenomena where anomalous diffusion are present [19]. Also, the field of non-extensive statistical mechanics has made an advance in the understanding of several systems where Boltzmann thermodynamic fails to explain the results [21][22][23]. For instance, anomalous diffusion and non-Gaussian velocity distribution have been observed in the context of Hydra cells [23], where maximum entropy densities associated with non-standard entropic measures were used to describe the motion of these cells. These densities are related to nonlinear diffusive process such as the generalized Fokker-Planck equation proposed in the context of the non-extensive statistical mechanics [22]. Motility and diffusive patterns have also been investigated in protozoa [9,10,12,[24][25][26][27][28]. For instance, T. brucei studies have focused on the movement of propulsion [9], flexibility and directionality [24], and body adaptations to the environment [25]. Cellhost interaction [12] for Leishmania spp and the flagellar beating [26] for T. cruzi have also been studied. Motility is strongly related to cell viability in all flagellate kinetoplastid species and it is widely used as a proxy measurement for viability [10,27,28].
Each protozoan specie has unique adaptations depending on the different living conditions. The investigation of the dynamics, diffusion and motion behavior of these microorganisms is an advance in the understanding of the microbial pathogenesis mechanism and in the field of diffusive patterns. However, a more complete and general understanding of motility patterns of these parasites is still lacking. To overcome this gap, we study the diffusive dynamics of causative agents of the neglected tropical diseases. By tracking the positions of these parasites, we present a complete characterization of their motility patterns. Specifically, we show that the spread of the trajectories is superdiffusive for short-times and that probability distributions related to the radial positions differs from the predictions of the usual diffusion equation. We further verify that the velocity time series of individual trajectories have long-range correlations and are well approximated by a generalized gamma distribution. Our results reveal some "universal" parasite motility patterns that could facilitate the identification of novel targets for therapeutic intervention. Furthermore, it could be expanded to screen aspects of cell viability.

Parasites maintenance
Leishmania amazonensis (MHOM/BR/75/Josefa strain) and Leishmania braziliensis promastigote forms were cultivated inside cell culture flasks containing Warren's medium (brain heart infusion plus bovine hemin and folic acid, pH 7.2) supplemented with 10% of fetal bovine serum (FBS). The parasites were incubated at 25°C for 48 h. Trypanosoma cruzi epimastigote forms (Y strain) were maintained in LIT medium (Liver Infusion Tryptose, pH 7.2), supplemented with 10% of FBS and incubated at 28°C during 96 h. Trypomastigote forms of Trypanosoma brucei brucei (EATRO-427 strain) were cultivated in HMI-9 medium, supplemented with 10% of FBS, incubated at 37°C and 5% CO 2 tension for 24 h. These incubation periods are essential to harvest the protozoa in the exponential growth phase.

Experimental setup, image acquisition and tracking
After the incubation periods, we have prepared a suspension containing about 6 × 10 6 parasites for Leishmania species, T. cruzi, and T. brucei in Warren, LIT, and DMEM (Dulbecco's Modified Eagle Medium, supplemented with 2 mM L-glutamine, pH 7.4), respectively. The mediums were not supplemented with FBS. Next, 10 μL of the protozoa suspensions were placed between a glass slide and coverslip to start the image acquisition. The thickness between the glass slide and coverslip is comparable to the size of the protozoa (*5 μm) [29][30][31], reducing it to a two-dimensional problem. We used the Motic's BA410E microscope equipped with a 5.0 Megapixel CMOS camera at a resolution of 800 × 600 pixels, acquisition rate of 10 frames/second and magnification of 20 ×. The area covered by the microscope at this configuration is 285.12 × 213.84 μm 2 . The length of the acquired videos was 10 minutes for each sample. We repeated this procedure three times for each protozoan. In order to extract the trajectories from the videos, we used a Matlab algorithm of motion tracking in image sequences [32]. After, we excluded the trajectories with less than 500 time-steps to ensure that we have long enough trajectories for statistical analysis. For these trajectories, we removed the first and last 50 steps due to imprecision of the algorithm in track the protozoa positions. By visual inspection, we further removed the trajectories for which the algorithm mistook two or more microorganisms at some step along the path. All trajectories were smoothed by applying a moving average filter of length 10 and are available in S1 Dataset. The viscosity of the culture medium was measured in a viscometer (Visco Star Plus) at 25°C and 50 rpm. The viscosity values and the number of trajectories analyzed for each protozoan are shown in Table 1.
In Fig 1A, we show typical trajectories of the L. amazonensis, that is,r i ðtÞ ¼ ½x i ðtÞ; y i ðtÞ, where x i (t) and y i (t) are the horizontal and vertical components of the position vectorr i ðtÞ in the time t for the i-th track. In Fig 1B, we plot the corresponding radial velocity time series,  [33]. Additional experiments considering different boundaries and thickness between the glass slide and the coverslip could further evaluate the effects of the walls on the motility of these protozoa. We have also calculated the average ( Fig 1C) and standard deviation (Fig 1D) of the velocities for each protozoan specie. It is worth noting that the velocities of the swimming microorganisms depend on the viscosity of the culture medium [25] and other values can be found in the literature due to different viscosity [9,[24][25][26]. For instance, it was found that the wildtype bloodstream form of the T. brucei can reach much higher velocities in the blood (around 30 μm/s) [25]. Trypanosoma spp use the motility as a tool to evade the immune cells and remove the bind between their surface and antibodies molecules. In the blood vessels, where the protozoa is adapted to survive, red and white blood cells behave as support for the flagellum to propel the cell body. The same behavior is observed in more viscous liquids, justifying the increase in speed of the parasite in these environments [25]. Other factors that affect the velocities include chemical cues, oxygen content, pressure, flow and confinement [25].
Our results about the motility of Leishmania spp promastigote and T. cruzi epimastigotes suggest a very similar behavior profile. Further details about these similarities are given in the results and discussion section.

Results and Discussion
We start by characterizing the spreading of the trajectories of all protozoa. In order to do so, we have considered the time series of the magnitude of the position vectorr i ðtÞ after subtracting its initial positionr i ð0Þ ¼ ½x i ð0Þ; y i ð0Þ. We thus measure the spreading by evaluating the  time dependence of the variance of the centered radial position where hr i ðtÞi ¼ 1 N k ðtÞ P N k ðtÞ i¼1r i ðtÞ is the average radial position and N k (t) is the number of available trajectories for the k-th specie longer than t steps (see Table 1). For usual diffusive processes the variance is expected to increase linearly on time, that is, σ 2 (t)*t. This usual behavior (or the Brownian motion) is related to absence of memory along the particle trajectory as well as indicates a finite characteristic scale for the position increments and for the waiting times between flights. However, more complex diffusive processes often display deviations of this linear behavior. When this happens, a typical behavior for the variance is a power-law dependence on time [34], where 0 < λ < 1 corresponds to subdiffusion and λ > 1 to superdiffusion. In our case, the evolution of the variances are shown in Fig 2A, where is evident that the spreading of the protozoan trajectories occurs much faster than the expected by a Brownian motion. We further observe an approximate power-law behavior over two decades of the temporal scale (t < 10 seconds). By least square fitting the log-log relationships (log σ 2 (t) versus log t), we find that values of λ are actually much larger than one. As shown in Fig 2B, λ ranges from 1.69 for the L. braziliensis to 1.93 for the T. brucei; therefore, the four flagellate protozoa studied here display a strongly superdiffusive behavior for short-times. A similar exponent was observed for the motion of intracellular particles with λ % 1.8 [35]. It is worth note that, biological swimmers can present some persistence related to their drive mechanism for short-times. We expect that the diffusion may eventually approach the usual regime for longer trajectories. Thus, the values of λ obtained here represent a initial short-time behavior, and other regimes could be observed for longer trajectories. In fact, we note from Fig 2A that the curves start to bend downwards for t > 10s. In Fig 2C, we show the exponent λ calculated within a window of size 30 seconds centralized in t w as a function of t w , where we note that the values of λ decrease with t w and approach the value expected by the usual diffusion. Although all studied species showed superdiffusive spreading for short-times, we noticed a intriguing aspect of movement related to velocities and diffusion analysis for T. brucei. This protozoan shows the smallest velocity and the greatest diffusion exponent in culture medium, suggesting that the T. brucei motility is more directional than the other protozoa species considered in this work. Directional motility probably occurs because T. brucei is a free extracellular parasite and spreads across several tissues. Trypomastigote forms of T. brucei are highly adapted to live in intercellular spaces [36][37][38]. During spreading in the host T. brucei parasites penetrate between cells in the tissue where there are collagen fibers that can facilitate or hinder motion. For trypanosomes to reach the maximum forward velocity a specific density of cells is required. If the density of obstacles resembles collagen networks the protozoa swim backwards in order to avoid getting trapped [25]. These parasites have an entire antigen surface called variant surface glycoprotein (VSG) [39]. The VSG is different among individuals in a population of T. brucei and prevents specific binding with antibodies that could kill the parasites [39]. We hypothesize that the presence of surface molecules such as VSG's, which reduces the difficulty of the cell motion in its surroundings media, allows T. brucei parasite to diffuse faster. Epimastigote forms of T. cruzi and promastigote forms of Leishmania spp exhibit similar velocities and diffusion exponent. These forms are faster than the T. brucei trypomastigote forms, have smaller diffusion exponent and have different strategies in this stage of the lifecycle. In the natural environment, T. cruzi epimastigote and Leishmania spp promastigote do not need to travel long distances. Specifically, epimastigote forms of T. cruzi must migrate to the midgut of the insect vector during a specific period of the life-cycle for proliferation [40]. In the case of promastigote forms of Leishmania spp, they must be phagocytized by mammalian cells to complete their life-cycle. At this stage, the Leishmania parasites secrete a substance called promastigote secretory gel (PSG). The PSG is a mucin-like gel and has been shown to be an important factor for the amastigote growth in the intracellular environment [41,42]. PSG production seems to be responsible for the formation of agglomerates (rosettes) and could play a role as constraint for the diffusive motion [43,44]. Overall, we suggest that changes in motility parameters and surface molecules that affect mechanical or physical constraints restricting the ability of cells to freely diffuse could significantly contribute to the virulence of these parasites.
Another striking feature of diffusive processes is related to the probability distributions of the positions. For the usual diffusion with radial symmetry, this distribution can be obtained by solving the following differential equation @Pðr; tÞ @t ¼ Dr 2 Pðr; tÞ ð 3Þ where D is the diffusion coefficient and r 2 is the two-dimensional Laplacian. By considering P(r ! 1, t) ! 0 and Pðr; 0Þ ¼ d 2 ðr Þ, we can show that probability distribution of the radial position is It is worth noting that this solution leads to a linear behavior for the variance over time, that is, σ 2 (t) = 4Dt. Furthermore, the distributions given by the Eq 4 are self-similar in time and collapse into a single curve, for the rescaled position ξ(t) = r(t)/σ(t). The usual diffusion equation (Eq 3) and its solution (Eq 4) can be understood as a null model for the distributions of the radial positions of the protozoa if we assume that they behave as Brownian particles, and deviations from this prediction is another indication of anomalous diffusion. We have calculated the time evolution of the empirical distributions of the radial positions. As shown in Fig 3A, we note that these distributions shift toward positive values of r while also become broaden over time t. After considering the rescaled position ξ(t) = r(t)/σ(t), we observe that all distributions collapse into a single curve (Fig 3B). This result demonstrates the empirical self-similar nature of the protozoa trajectories, but also reveals remarkable deviations between the empirical distributions and the expected by the usual diffusion equation (gray dashed lines in Fig 3B).
In an attempt to find a better description for these empirical distributions, we propose to replace the Gaussian term in Eq 4 by a stretched Gaussian (characterized by another parameter δ > 0). Specifically, we have considered the following probability distribution for the radial where GðxÞ ¼ R 1 0 y xÀ1 e Ày dy is the gamma function. We observe that δ = 2 recovers the usual Gaussian term; however, for 0 < δ < 2 the tail of this distribution goes to zero slower than the usual case, whereas for δ > 2 it decays faster. A similar generalization was proposed by Richardson [45] in the context of atmospheric diffusion by considering a spatial-dependent diffusion coefficient, and is considered one of the first anomalous diffusion equation. The distribution given by the Eq 6 is also self-similar and for the rescaled radial position ξ(t) = r(t)/σ(t) it can be written as We have thus adjusted Eq 7 to the window average values of the empirical distributions via least squares fitting. The continuous lines in Fig 3B show that agreement is far from being perfect, but this generalization is a better description when compared with the distribution emerging from the usual diffusion equation (Eq 5). Therefore, the diffusive motion of the protozoa studied here displays simultaneously an anomalous (enhanced) scaling of variance σ 2 (t) as well as radial position distributions with much longer tails than the expected by Brownian swimmers. The values of the best fitting parameters δ are shown in Fig 3C, where we observe that the L. amazonensis and the T. brucei display larger tails (δ % 1) when compared with the L. braziliensis and the T. cruzi (δ % 1.25). It is also worth noting that the values of δ should be related to values of λ. In fact, the time dependence of the variance evaluated from Eq 6 is σ 2 (t) / t 2/δ , leading to δ = 2/λ, a relationship that is roughly valid for the values δ obtained via the fitting procedure ( Fig 3C).
In order to further characterize the diffusive motion of the protozoa, we investigate the radial velocity time series v i (t) related to the individual trajectories (Fig 1B). We first ask whether these time series have short or long-range correlations. To answer this question, we have applied the detrended fluctuation analysis (DFA) [46,47] to these time series. This technique consists of four steps: i) first, we define the integrated profile YðtÞ t v i ðtÞ with t max being the length of the i-th time series; ii) next, we split Y(t) into N n = t max /n non-overlapping segments of size n; iii) for each segment, a local polynomial trend (here we have used a linear function, but higher orders do not change our results) is calculated and subtracted from Y(t), defining Y n (t) = Y(t) − p ν (t), where p ν (t) represents the local trend in the ν-th segment; iv) finally, we calculate the root-mean-square fluctuation func- where hY n (t) 2 i ν is the mean square value of Y n (t) over the data in the ν-th segment. If the velocity time series is self-similar, the fluctuation function F(n) blue (t = 2s) to red (t = 10s). B) The same distributions after considering the rescaled position ξ(t) = r(t)/σ(t).
The black dots represent window averages of the CDFs and the error bars stand for one standard deviation. The gray dashed lines represent the distribution given by the normalized Gaussian distribution (Eq 5). The continuous lines are the cumulative version of the distribution of Eq 7 and the best values for the fitting parameter δ are shown in the plots. C) Comparison between the values of δ obtained via least squares fitting the window averages of the CDFs and the ones expected by the time dependence of the variances, that is, δ = 2/λ. presents a power-law dependence on the time scale n, that is, F(n)*n h , where h is the Hurst exponent. If h = 1/2 the velocities are either uncorrelated or short-range correlated, whereas h 6 ¼ 1/2 indicates that the time series is longe-range correlated.
We have applied the above procedure to all velocity time series and a typical behavior for the fluctuation function F(n) is depicted in Fig 4A. In this log-log plot, we adjust a linear model for obtaining the Hurst exponent h, which is h = 1.12 for the original time series and is close to 1/2 for random shuffled versions of the time series. We calcule the Hurst exponent for all time series and the average values for each protozoan is shown in Fig 4B. These averages are practically indistinguishable from each other and indicate that the velocities of the protozoa have long-range correlations. Furthermore, these velocities display a strong persistent behavior (since h % 1), that is, positive increments in the velocities are followed by positive increments and negative increments are followed by negative increments much more frequently than by chance. We have also evaluated the distributions of h for each protozoan specie for all individuals (Fig 4C), where we note that these distributions peak around h % 1 and that they are quite overlapped.
Another intriguing question is whether the velocity distributions of the protozoa exhibit a particular functional form (Fig 5A). In this case, the two-dimensional Maxwell-Boltzmann distribution (or Rayleigh distribution) is a natural null model for the protozoan velocity (T is a parameter). This distribution represents the velocity of two-dimensional gas particles in thermodynamic equilibrium (at a temperature T) and also emerges when evaluating the distribution of the magnitude of velocity vectors whose the components are uncorrelated and normally distributed (with zero mean) [48]. This functional form has been found to describe quite well the velocities of humans in a very peculiar situation (a mosh pit) [49] and it has also been used in the attempt of modeling the velocity distributions of Hydra cells in a two-dimensional setup [23]. Deviations from this model give clues about weather the velocities are correlated or not. In our case, we have tested the (twodimensional) Maxwell-Boltzmann hypothesis by adjusting this distribution for each velocity time series and verifying the goodness of the fit via Kolmogorov-Smirnov test. This hypothesis was rejected for almost all velocity time series (% 99%), a result that somehow agrees with the more complex behavior observed in the correlation analysis of these time series. Aiming to find a better description for these empirical distributions, we have considered the generalized gamma distribution where α and γ are the shape parameters and β is a scale parameter. Despite of being an ad hoc generalization, it is worth noting that this distribution recovers the two-dimensional Maxwell-Boltzmann (for α = 1, γ = 2 and β 2 = T), and it has been employed by several authors as a wind speed model [50][51][52][53]. For our data, the Kolmogorov-Smirnov test cannot reject the generalized gamma hypothesis in about 50% of trajectories (see S3 Fig), an improved description when compared with the two-dimensional Maxwell-Boltzmann distribution. Furthermore, we have tested for the usual gamma (Eq 9 with γ = 1), log-normal, Weibull and q-exponential distributions, finding that they do not outperform the generalized gamma description. In Fig 5B, we show the average values of the best fitting parameters α, γ and β. These values are practically indistinguishable among the four protozoa. An exception occurs for the T. brucei, which is characterized by a significantly smaller value of β. This fact is a direct consequence of the small standard deviation observed for the velocities of T. brucei (see Fig 1D), since the standard deviation of v calculated from Eq 9 is proportional to β.

Conclusions
We presented a description of the motility patterns of four trypanosomatid flagellate protozoa. By analyzing the time evolution of the positions of these protozoa, we identified that the spreading of their trajectories are strongly superdiffusive for short-times and characterized by self-similar probability distributions with longer tails than the expected by the usual Brownian motion. We also investigated the velocities of these protozoa, finding out that they have long- range correlations and present a strong persistent behavior. We further observed that the velocity distributions cannot be described by a two-dimensional Maxwell-Boltzmann distribution (the natural candidate for a random process). Instead, a generalized gamma distribution showed to be a better description. Thus, our results show that the motility patterns of these protozoa are anomalous in several ways and also reveal that these four protozoa exhibit similar behaviors, despite their morphological differences. Deviations from the behaviors reported here can be employed as an indicator of drug activity in drug tests.