Liver Environment–Imposed Constraints Diversify Movement Strategies of Liver-Localized CD8 T Cells

This article is featured in Top Reads, p.

Searching for resources is a fundamental property of many biological agents, starting from individual cells to large animals. Although different agents involve various search strategies, two types of search are of particular importance due to their fundamental nature: Brownian walks and Lévy walks/flights (1). Brownian walkers display short displacements, typically consistent with a thin-tailed distribution that has a finite mean and variance, with random turning angles; in Brownian walks, mean squared displacement (MSD) increases linearly with time (1–4). If a walker exhibits preference in movement direction similar to the previous movement direction, this results in correlated random walks during which MSD increases transiently faster than linearly with time and can be characterized by so-called persistence time (5–7). In contrast, Lévy walkers display both short and long displacements consistent with a fat-tailed distribution that has either finite mean and infinite variance or both infinite mean and variance; in Lévy walks/flights, MSD increases faster than linear with time (1, 3). Various combinations of these two extremes, such as composite correlated random walks and generalized Lévy walks, have also been studied (e.g., [8]). Understanding walk types of various biological agents has been a subject of many investigations, and several techniques, such as analysis of movement length distribution and MSD dynamics, have emerged as relatively standard (7–17).

For adequate protection, T cells of vertebrates survey tissues for infection by moving within tissues. Yet, specific movement strategies and impact of tissue environments on the movement patterns remain incompletely understood for most tissues in vertebrates such as mice. In lymph nodes (LNs), naive T cells seem to search randomly for APCs, acting as Brownian walkers undergoing short displacements in correlated random walks, although some role of attracting chemokines toward dendritic cells (DCs) expressing cognate Ag has been proposed (3, 6, 18–20). In the skin, activated or memory CD8 T cells tend to perform Brownian walks (21), and in explanted murine lungs, T cells display intermittent migration resulting in normal diffusion and a large variability in average velocities between individual cells, although the specific type of the walk (Brownian, Lévy, or a combination) has not been determined (22). Interestingly, Toxoplasma gondii–specific activated CD8 T cells in murine brains undergo Lévy walks, although specific mechanisms (whether cell intrinsic or environmentally driven) allowing such T cells to display rare long displacements have not been determined (8).

The liver is a large organ in mice and humans and is the target of multiple pathogens (23). Apicomplexan parasite Plasmodium, the causative agents of malaria, infects the liver following inoculation of parasites by probing mosquitoes (24–26). At this stage of the Plasmodium life cycle, the intracellular liver stage parasites are susceptible to clearance by vaccine-induced CD8 T cells. Remarkably, liver-localized CD8 T cells alone are capable of eliminating all liver stages and preventing clinical malaria (27, 28). Given that a murine liver has ∼10⁸ hepatocytes, how exactly CD8 T cells locate all parasites within the 48-h span of the parasite liver stages is unknown (29). Furthermore, although we and others have documented movement of liver-localized CD8 T cells, the walk types employed by these cells and their potential underlying mechanisms remain incompletely understood (30–32).

Typical intravital imaging experiments allow researchers to collect three-dimensional (3D) coordinates of T cells in tissues over time, and, with such data, one can construct several movement characteristics to determine the type of walk cells exhibit (5–8, 12, 33). Specifically, the rate of MSD increase with time, the distribution of movement lengths, and the distribution of turning angles often allow researchers to determine the walk type (7, 12, 34). Scaling the distribution of movement displacements generated for different imaging frequencies also allows researchers to determine the type of cell walk, because Lévy walk–based movement distributions are thought to be relatively invariant to such rescaling (6, 8, 33). Such methods, however, have not been applied to quantitative characterization of the movement program of T cells in the liver.

Although several studies provided strong evidence that T cells undergo a specific movement strategy, mechanisms explaining such a movement strategy in general remain unknown. The constrained environment of the skin epithelium is likely responsible for the movement patterns of skin-resident CD8 T cells (21). Activated CD8 T cells in murine brains undergo Lévy walks (8); however, whether such a walk type is driven by a cell-intrinsic program (as suggested by Harris et al. [8]) or arises due to other factors has not been determined. Indeed, previous ecological studies suggested that turbulence of fluids or air may explain patterns of movement of individual cells or birds, respectively, but the actual features of turbulence that allow agents to undergo Lévy walks were not directly measured (35, 36). Furthermore, in a theoretical study, Volpe and Volpe (37) also showed that environmental constraints impact the movement patterns of agents and the efficiency of the search for targets, suggesting that in addition to cell-intrinsic programs, other factors are likely to impact movement patterns of T cells (20). The importance of physiological factors such as blood flow in highly vascularized tissues such as brain or liver in determining movement patterns of T cells has not been evaluated, as far as we know.

In this study, with the use of intravital microscopy of live mice, we provide (to our knowledge) the first accurate quantification of the movement pattern of liver-localized CD8 T cells (and openly share these data with the community). We show that in most circumstances, CD8 T cells undergo short, Brownian-like displacements characterized by a thin-tailed distribution that yet result in transiently superdiffusive displacements; that is, most T cells in the liver undergo correlated random walks. Importantly, however, the duration of the transient superdiffusion is longer than that observed in other systems (e.g., by T cells in LNs or murine brains). By accurately quantifying the structure of liver sinusoids, we show that transiently superdiffusive behavior of liver-localized T cells naturally emerges for Brownian-like, persistent walkers when simulations of movements are performed in the liver. Interestingly, our analysis of the structure of fibroblastic reticular cells (FRCs) of LNs suggested that LN-localized T cells should be superdiffusive for lengths of time similar to those for liver-localized T cells; however, LN-localized T cells are superdiffusive for only a few (∼5) minutes, suggesting that other factors such as the width of the FRC paths or chemical cues may also reduce the duration of movement persistence of T cells in the LNs. Interestingly, we found that in some circumstances, T cells in the liver may also undergo Lévy flights—large displacements occurring nearly instantaneously—due to detaching from the sinusoidal epithelium and flowing with the blood flow. Taking these data together, we have carefully characterized movement patterns of liver-localized CD8 T cells and provided examples of how environmental constraints (structure of liver sinusoids) and physiology (blood flow) may regulate the observed T cell displacements.

Experimental data on movement of wild-type (WT), GFP-expressing, TCR transgenic CD8 T cells (OT-I) have been generated similarly to data in our previous study (31). OT-I cells were activated in vitro as follows. Splenocytes were incubated with 1 μg/ml SIINFEKL peptide (Biomatik) at 37°C for 30 min, washed, and cultured in complete RPMI 1640 (10% FCS, 2 mM l-glutamine, 1 mM Na-pyruvate, 100 U/ml penicillin/streptomycin, 5 mM HEPES, 20 μg/ml gentamicin, and 50 μM 2-ME) at 37°C, 5% CO_2. Two days later, cells were passaged and cultured in fresh media supplemented with (12.5 U/ml) recombinant IL-2 (PeproTech). Cells were passaged again after 24 h and given fresh media with IL-2. Four days after activation, cells were ready for collection and use. Cells were collected from tissue culture flasks, washed, layered over Histopaque 1083 (Sigma-Aldrich), and centrifuged at 400 relative centrifugal force for 30 min at 21°C to remove dead cells. Live lymphocytes were collected from the interface between the Histopaque and cell suspension media. Cells were washed a further three times before transfer to mice.

To label the sinusoids, Evans blue or rhodamine-dextran was injected i.v. into mice (50 μg/mouse) immediately before imaging of the livers (Fig. 1A, 1B). In some cases, static images of the livers were acquired using a two-photon microscope. In total, we generated images of the liver sinusoids from five mice (with or without T cells).

Mice were prepared for multiphoton microscopy essentially as previously described (31). Briefly, mice were anesthetized with a mix of ketamine (100 mg/kg) and xylazine (10 mg/kg). Throughout the surgery and imaging procedure, the mouse temperature was maintained at 37°C using a heating mat attached to a feedback probe inserted in the mouse rectum. A lateral incision was made over the left lobe of the liver, and any vessels were cauterized by applying light pressure to the vessel until clotting occurred naturally. The mouse was then placed in a custom-made holder. The liver was then exposed and directly adhered to a coverslip that was secured in the holder. Once stable, the preparation was transferred to a FLUOVIEW FVMPE-RS multiphoton microscope system (Olympus) equipped with an XLPLN25XWMP2 objective (25×, numerical aperture 1.05, water immersion, 2-mm working distance). For analysis of the motility of liver-localized CD8 T cells, a 50-μm z-stack (2 μm/slice) was typically acquired using a standard galvanometric scanner at a frame rate of two to four frames per minute.

Other datasets for movement of WT and LFA-1–deficient CD8 T cells or WT CD8 T cells 1 or 4 wk after immunization with radiation-attenuated Plasmodium berghei sporozoites $Pb-CS5M$ were generated previously (31).

Raw imaging data were analyzed with Imaris (64 bit) software version 9.1.2 (Bitplane). Tracking of individual cells in a z-stack was performed using the Spots function in the surpass mode and in a single slice by the Surfaces function in surpass mode. Detection of individual cells relied upon their relative fluorescence intensity and size (diameter ≥9 μm), and the remaining default function settings were used for the calculation of the motion tracks using an autoregressive motion algorithm. Tracks were also manually inspected and curated. Basic cell characteristics such as average track speed, meandering index, and arrest coefficient were calculated using Imaris. To calculate the percentage of T-cell movements that were occurring in the sinusoids, we generated a colocalization channel that calculated colocalized signal of both red and green fluorescence within each frame. To do this, we used the colocalization function in Imaris, with the threshold defined in Imaris as 16 for uGFP OT-1 cells, which subtracted background fluorescence of the liver, and 3000 for rhodamine dextran. We then calculated the total number of frames in which a T cell was tracked using the unique cell identifier and compared this with the total frames in which each individual T cell was detected in the designed colocalization channel. Rhodamine dextran fluorescence was dependent on the depth of imaging, and in the lower z-stacks, the frequency of T cells with colocalized signal was poor (results not shown).

In our analysis, we used 3D coordinates of Plasmodium-specific CD8 T cells in the liver over time from the following experiments:

1. Data on movement of GFP-expressing activated CD8 T cells called “New13” and “New26”: In the New13 dataset, data were generated from three mice, and frames were sampled at 14.63 (14.63–14.64) s apart. In the New26 dataset, data were generated from two mice, and frames were sampled on average at 25 (24.29, 25.51) s apart. To increase the power in the analysis, we generated a combined dataset New1326 by merging data from datasets New13 and New26. We ignored the middle time points in time frames in the New13 dataset, thus computing time frames at ∼26-s frequency.

2. Previously published data on movement of WT and LFA-1–deficient CD8 T cells in the liver (31) resulted in Lfa1 WT and Lfa1 knockout datasets. Data were from seven mice, and frames were sampled, on average, at 27 (24.65–29.56) s apart.

3. Previously published data on movement of Plasmodium-specific CD8 T cells following vaccination of mice with radiation-attenuated sporozoites (RASs) (31). In the week 1 dataset, data were generated from two mice immunized with RAS 1 wk before. Frames were sampled, on average, at 6.51 (6.51–6.51) s apart. In the week 4 dataset, five mice were RAS immunized 4 wk before imaging, and frames were sampled, on average, at 3.25 (3.25–3.26) s apart.

4. Data on movement of Plasmodium-specific (OT-1) and lymphocytic choriomeningitis virus–specific (P14) CD8 T cells moving in livers of P. berghei–infected mice (38).

In the Imaris-processed videos, we found that there were random missing time frames in tracks in nearly all data. This could happen when Imaris would associate a cell that went out of focus for one time step with other time steps. Because having a consistent time step to determine cell movement characteristics such as MSD and movement length distribution was critical, we cleaned the data in the following ways. We split tracks with missing time frames and converted them to individual cells for the analyses. For example, if one time frame was missing for a cell, the track was split into two with the second part being assigned as a new “cell.” This created many more cells with shorter tracks. Taken together, cleaning up the data yielded the following. For the New13 dataset, we originally had 56 T cells out of which we generated 14 new split cells with 1 or 2 splits resulting in total of 70 T cells with 2565 time frames (movements). For the New26 dataset, we had 213 T cells with 51 new split cells with 1 or 2 splits from the original tracks totaling 264 T cells, with 7736 time frames. For the combined New1326 dataset, we had 334 T cells with 9018 time frames (smaller number of time frames in this dataset than in New13 + New26 is due to the different in time step in the two datasets). For the Lfa1 WT dataset, we originally had 198 T cells out of which we generated 19 new split cells with 1 or 2 splits from the original tracks totaling to 217 T cells with 6210 time frames. For the Lfa1-knockout dataset, we originally had 154 T cells out of which we generated 16 new split cells with 1 or 2 splits from the original tracks totaling 170 T cells, with 5097 time frames. For the Wk1 dataset, we had 99 T cells and generated 3 new split cells with 1 or 2 splits from the original tracks, totaling 102 T cells with 2525 time frames. For the Wk4 dataset, we had 44 T cells with 4169 time frames that required no cleaning. Most of the analyses have been repeated on “clean” data in which tracks with missing time frames were removed with similar conclusions reached.

Data have been provided by Mario Novkovic and colleagues from their previous publication (39). In short, 3 × 10⁶ naive P14 TCR transgenic T cells were transferred i.v. into naive mice. Three to twenty-four hours after T cell injection, intravital imaging was performed on the right popliteal LNs. To detect FRCs, anti–yellow fluorescent protein mAbs were used, and samples were imaged using a confocal microscope (LSM-710; Carl Zeiss) (39). In total, we had 3 sets of images of LNs with FRCs labeled (∼30 images per set with 2 μm depth). In addition, we analyzed 3D movement data of T cells in the LNs (six videos in total) produced in the same publication (39). The data were also cleaned, and tracks for which there were not always consecutive time frames available were split as described above. Analyses were repeated with clean data in which cells with missing time frames were removed; similar results were found. These data were generated at a frequency of 20 s per frame.

For each of the datasets, we calculated several basic characteristics. To calculate MSD, we used all tracks in the dataset in which the time frequency of measurements was approximately the same:

where d_i,t is the displacement of an ith cell from its initial position to the position at time t and N(t) is the number of cells for which coordinate data were available at time t. n_i is the number of steps in each ith cell. Note that this is different from movement lengths, which are movements T cells make between sequential time frames. To characterize the MSD, we used the following relation $MSD(t)=ctγ$⁠, where t is lag time or delay and γ was estimated by log-log transforming MSD and t. Estimate of γ < 1 suggests subdiffusion of cells, whose walks are constrained by boundaries, regardless of their step-length distribution; γ = 1 suggests Brownian (or normal) diffusion, whose walks are affected by neither boundaries nor transportation; and γ > 1 suggests superdiffusion, whose walks are indicative of transportation or directed movements. Root MSD, $ζ=RMSD$ was calculated as

where $MSD$ is defined in Eq. 1. To calculate the turning angle for cell movement between time steps (t, t + 1) and (t + 1, t + 2), we calculated the vectors of movements between two time points and then calculated the angle between these two movement vectors. To estimate the decay in the correlated movements, which is another alternative method to check if the walks are Brownian, of T cells over time, we calculated the velocity autocorrelation function given by the time-averaged dot product $<v¯1·v¯t>$⁠, where $v¯1$ is the vector for the initial T cell velocity (movement vector for a cell from t = 0 to t = 1) and $v¯t$ is the vector for later movement (movement vector for a cell between times t − 1 and t), with ϕ_1,t being the angle between the two vectors $v¯1$ and $v¯t$⁠, such that we get the normalized velocity autocorrelation as:

where, by definition, c_V varies between −1 and 1.

To investigate which model is most consistent with the distribution of movement/step lengths of activated CD8 T cells in the liver or distribution of branch lengths of liver sinusoids, we fit several alternative distributions to the data. Cauchy distribution is defined as

where $r≥rmin$ is the movement length, $rmin$ is the location parameters, and the mean of the distribution is undefined. Half-Gaussian distribution with r ≥ 0 is defined as

where the mean of the half-normal distribution is $r¯=σ2/π$⁠.

Stable Lévy distribution is defined as

where γ and $rmin$ are the shape parameter and location parameter, respectively. The stable Lévy distribution has infinite mean and variance.

Generalized Pareto (GP) distribution is defined as

where k is the shape parameter, σ is the scale parameter, and θ is the location parameter. In this distribution, the mean is $r¯=θ+σ/(1−k)$⁠.When k = 0 and θ = 0, GP becomes exponential distribution

where λ = 1/σ and mean is equal to $r¯=σ$⁠. When k > 0, θ → 0 GP becomes the Pareto (or power law) distribution

where $rmin=σ/k$ (scale parameter), α = 1/k, μ = α + 1 (shape parameter), and $r¯=(μ−1)rmin/μ$⁠. Note that in the Pareto distribution $r≥rmin$⁠. The important property of the Pareto distribution, describing either Brownian or Lévy walks, is dependence of the mean and variance on the value of μ (and thus α). For μ > 3 (α > 2), both mean and variance of the Pareto distribution are finite, corresponding to Brownian walk. In contrast, when μ ≤ 3 (α ≤ 2), mean is finite and variance is infinite, corresponding to Lévy walks. Finally, when μ ≤ 2 (α ≤ 1), both mean and variance are infinite, corresponding to bullet motion (1).

Generalized extreme value distribution is defined as

where $t(r)={(1+ξ(r−μ)σ)−1/ξ,if ξ≠0,e−(r−μ)/σ,otherwise.$

We fit the models (Eqs. 4–9) to either actual movement lengths r of cells between two time points or scaled movements ρ. Scaled movements are defined as $ρ=r/r¯$⁠, where $r¯$ is the average movement in the particular dataset (⁠$r¯=1n∑i=1nri$ for n movements). Using scaling of movement lengths by the square root of the mean of squared movement lengths (⁠$r2¯$⁠) gave similar results. It is important to emphasize that scaled movement ρ allow one to convert cell tracks to be independent of the effect of different mean speeds of cells and independent of the effect of different sampling frequency (8). The likelihood of the model parameters given the data, is given as

where x_i is either the cell movement length data x_i = r_i or the scaled movement data x_i = ρ_i(t) for i = 1…n, for n data points, for fixed time interval t. The probability density functions f(x_i|m₁,m₂,…) are given in Eqs. 4–9. Parameters m₁, m₂, … of the models were estimated by minimizing the log of negative likelihood, $L=−lnL$⁠.

In addition to fitting different distributions to the movement length/displacement data, we performed an analysis to determine the shape parameter of the powerlaw/Pareto distribution μ that is most consistent with the tail of the displacement distribution data (34). Given that Pareto distribution has two parameters (⁠$rmin$ and μ; Eq. 9) one needs to define $rmin$ that determines the cutoff in the data and that allows the best description of the tail of the data with Pareto distribution. The steps to define $rmin$ are as follows (34): 1) for each possible choice of $rmin$ from the smallest to the largest unique values of r, we computed μ by maximum likelihood estimate given by the analytic expression $μ^=1+n(∑i=1nlnrirmin)−1$⁠, where r_i, i = 1 … n, are the observed values of r and $ri≥rmin$⁠. In doing so, we also set a limit to the number of data points in the analysis by setting a criterion by which the analysis is terminated if the SE on $μ^$⁠, approximately given by $σ=μ^−1n$⁠, was >0.1, because higher-order corrections of the SE are generally positive. This condition was set because the estimates of μ can be erroneously large when $rmin$ is so large that the number of observed $r≥rmin$ is very small. It was also assumed that μ > 1, because the distributions with μ ≤ 1 are not normalizable and hence not realistic in nature. 2) We selected $rmin$⁠, at which the Kolmogorov-Smirnov goodness-of-fit statistic D was the minimum of all D values given by the fits over all possible unique $rmin$ values of the data. Here, we also made a finite correction to the estimates of μ for small sample sizes (n < 50) in addition to the condition set by the SE of μ.

From the above analyses, we carried forward the estimated $rmin$⁠, μ, and D to the next stage of testing the probability of the data, given that the estimated model is the plausible one, by the Kolmogorov-Smirnov goodness-of-fit test. To do that, we generated data from the hypothesized (or estimated) Pareto model, sampling using a semiparametric bootstrap method, computing D for each sample, which is the synthetic D, and thus computing a fraction of those was greater than the empirical D of the given data, which yields the p value of the hypothesis. We then computed the uncertainty in the estimates of $rmin$⁠, μ, and n, given by SDs, by resampling the data with replacements by Markov chain Monte Carlo simulations.

We analyzed 23 z-stacks of Evans blue–stained 512 × 512–μm images of the liver (mouse 1, 171207; mouse 2, 190624), as well as 3 sets of images/videos with T cells in which liver sinusoids were labeled with rhodamine-dextran. Resulting structures were similar based on calculated characteristics such as branching angle distribution and straightness index (results not shown). Light-adjusted 2D images were converted from RGB into grayscale using rgb2gray, filtered through MATLAB filters (1) imfilter with Gaussian N(10,2) to remove noise, (2) medfilt2 with “symmetric” setting to remove all the lines from the image, (3) uint8 to subtract image without lines from image with lines, (4) setting threshold at >130 of the subtracted image, and (5) removing more noise by imopen, and saving all 2D images in a stack converting to a digitized 3D image, with “1” giving the sinusoids and “0” giving the hepatocytes and other unstained tissues (Fig. 3). The digitized 3D image was then skeletonized, with 1’s giving the mean path of sinusoids and 0’s elsewhere, converting to a voxel image, using parallel medial axis thinning method as per the MATLAB vectorized implementation of the algorithm (40). The 3D binary voxel skeleton was converted to a network graph, computing nodes (junctions) and edges (branches) using previously published methods (41) and mapping the respective skeleton to a 3D graphics function. The nodes and the branches were used to compute branch angles, branch lengths, node-to-node direct distances, and number of branches per node. The straightness of branches metric was computed as the node-to-node direct distance over the branch length. The binary 3D liver image was used to run simulations of cells walking on the 1’s on the lattice, the scale of which can be decided, for example, as 1 U equals 0.01 μm (Fig. 4).

We simulated cell walks on (1) digitized 3D real sinusoidal branch structures (five imaging sets) and (2) no structures (five simulations), imposing no physical constraints and giving freedom for cells to turn in any direction after every step. For each generated structure, we selected a random point on a mean sinusoidal path (i.e., on the digital skeleton given by 1’s, on the matrices generated from digitized 3D sinusoidal structures of 1’s and 0’s in the real case, and on a branch position coordinate in the simulated case) as the starting point of a T cell. Then we selected a cell step length from the Pareto distribution with a given parameter μ. Note that we considered all movements as forward movements, because reverse movements of cells along the sinusoids were the minority in the data (∼20%; Supplemental Fig. 3). The ability of cells to move in one preferred direction in constant external conditions can be viewed similarly to the physics postulate of momentum conservation (or “biological inertia” [42]). We disregarded deviation of tracking cell centers from the mean sinusoidal path because the cell diameter, on average, 10 μm, was greater than the cylindrical diameter of the sinusoids, 7 μm, and thus T cells likely occupied most of the cross-sectional area of the sinusoids. Note also that step lengths and branch lengths were thus simulated disregarding their autocorrelations, assuming they are independent and identically distributed. For cell walk simulations on no structure, the cell turning angles were random from a uniformly spherical distribution. Note that the walk simulations of cells and sampling coordinate were terminated once a cell hit the boundary of the structure before the simulation time was elapsed in the case of real sinusoidal structures. We computed MSD of all simulation scenarios. To simulate T cell walks in digitized 3D network of FDCs in LNs, we proceeded similarly to the method for liver sinusoids, assuming that T cells can only walk in the middle of the network path and turns can be taken only when a T cell reaches a branching point.

To calculate the efficiency of the T cell search for the malaria liver stage, we calculate the time to target because this is the most relevant for CD8 T cell–mediated protection against malaria. In simulations, the T cell started the search from a random point but with a given distance to the parasite. We assumed that the T cell found the infected cell when it was within 40 μm of the target, an approximate size of the murine hepatocyte (43, 44). We repeated the simulation with five different distances to targets (results not shown). For 3D simulations, if the searching T cell reached the boundary, we assumed that it missed the target (and thus was not included in the calculation of the average time). To simulate T cells searching for a target in the liver, we simulated T cells moving forward on one path, effectively reducing the search problem to 1D. This is reasonable, given that the difference between Brownian and Lévy walks is only in the movement length distributions and not in turning angles that cells take (1). To simulate n T cells searching for a single sample, assuming that search is done by different cells independently, we resampled times to target with replacement n times from the dataset on a single T cell searching for a single target and calculated the minimum of the time in the sample. This procedure was then repeated 1000 times to generate a distribution to target when n T cells search for the infection site.

All animal procedures were approved by the Animal Experimentation Ethics Committee of the Australian National University (protocol numbers A2016/17, 2019/36). All research involving animals was conducted in accordance with the National Health and Medical Research Council’s Australian Code for the Care and Use of Animals for Scientific Purposes and the Australian Capital Territory Animal Welfare Act 1992.

T cell movement/track data generated in this study or in our previous publication are provided in an Excel spreadsheet with different tabs for different datasets. Imaging data of liver sinusoids (TIFF images) are also provided. The data from a previous publication on movement pattern of naive CD8 T cells in murine LNs and imaging data of FRC networks in the LNs (39) are also provided. These files are stored together with the MATLAB codes used to process and analyze the data on GitHub.

All analyses have been primarily performed in MATLAB (version 2019b) and codes used to generate most of the figures in the paper are provided on GitHub: https://github.com/vganusov/Levy_walks_liver.git. All code questions should be addressed to H.R. (hrajakar@utk.edu).

To understand what strategies Plasmodium-specific CD8 T cells may be using to survey the liver for infection, we have established a protocol to generate stable numbers of liver-localized CD8 T cells (31). In these experiments, naive CD8 T cells are activated in vitro and then transferred i.v. into naive mice that either are left uninfected or can be later infected with Plasmodium sporozoites (Fig. 1A) (31). By injecting fluorescent label (rhodamine dextran) i.v. immediately before surgery and imaging, we can efficiently label the liver sinusoids (Fig. 1B). Intravital imaging showed that CD8 T cells in the liver are highly motile cells with the average speed of 6.6 μm/min, a high meandering index of 0.52, and a low arrest coefficient of 0.25 (Fig. 1C and Supplemental Video 1). Importantly, these CD8 T cells were observed to be located inside the liver sinusoids most of the time (89–100%; (Fig. 1D). Due to incomplete labeling of the sinusoids by rhodamine-dextran and reduced red signal from deeper liver areas, it is likely that a lower bound is an underestimate, and, in this system, nearly all T cell movements (and, by inference, most T cells) were localized in the sinusoids and thus were intravascular. As compared with the blood flow rate in the liver, a relatively low average speed of these liver-localized T cells and T cells in our previous experiments with high-frequency intravital imaging suggests that the cells are crawling on the endothelial wall of the sinusoids and not flowing with the blood (31). Thus, our robust experimental setup allows for the first time, to our knowledge, accurate quantification of movement patterns of liver-localized CD8 T cells along with the well-defined tissue structures (i.e., hepatic sinusoids).

In our experiments, movements of liver-localized Plasmodium-specific CD8 T cells were recorded with two-photon microscopy. Interestingly, T cells in the liver displayed transiently superdiffusive overall displacement characterized by MSD increasing faster than linearly with time (with slope γ = 1.55) for up to 20 min (Fig. 2A). Saturation in MSD with time in these and other experiments typically arises in part due to rapidly moving T cells leaving the imaging area (Supplemental Fig. 1). Superdiffusion is a hallmark feature of Lévy walks (3, 45); yet, liver-localized T cells only displayed short, Brownian-like movements (Fig. 2B). Tail analysis (34) of the distribution of actual or scaled movement lengths (movement lengths normalized by mean movement length) suggested the shape parameter of the power law distribution μ = 5.6, an indication of Brownian (μ > 3) and not Lévy (μ ≤ 3) walks. Interestingly, scaling the movement length distribution for subsampled data (i.e., for data sampled at a lower frequency than was used at imaging) resulted in relatively similar distributions (Supplemental Fig. 2), suggesting that even thin-tailed movement length distributions may be similar for subsampled scaled data. Thus, large overall cell displacements with small movement lengths arose due to small turning angles and a strong correlation between initial and subsequent velocity vectors for at least 10–15 min, resulting in cells moving in one direction for substantial periods of time (Supplemental Fig. 3).

To complement the movement length distribution tail analysis, which is the standard in the field (7, 8, 12), we fitted a series of alternative mathematical models characterizing the distribution of movement lengths of T cells. We found that GP distribution fitted the data with the best quality (Fig. 2C, 2D, and Supplemental Table II). Best fit parameters of the GP distribution predicted a finite mean and variance for the movement length distribution, which again are the chief features of the Brownian and not Lévy walkers (1). Thus, in our experiments, liver-localized CD8 T cells performed correlated (persistent) random walks with transient superdiffusive movements that were longer than those observed for CD8 T cells in the brain (8).

To verify that this result was not due to a potential impact of the injected rhodamine-dextran dye, we analyzed four other datasets from our previous study on movement of liver-localized CD8 T cells (31). These experiments involved transfer of in vitro generated WT and LFA-1–deficient activated CD8 T cells or in vivo generation of liver-localized T cells by transfer of naive OT-I cells to mice followed by immunization of these mice with malaria RASs expressing the SIINFEKL epitope recognized by the T cells (29, 31). Both in vitro WT activated cells and in vivo activated cells imaged at a memory time point (4 wk) after immunization also displayed transiently superdiffusive behavior with γ > 1 for 10–20 min (Fig. 2E, Supplemental Fig. 4). Interestingly, LFA-1–deficient T cells, which cannot strongly attach to endothelial cells and roll under the blood pressure (31), displayed displacement closest to Brownian and yet were transiently superdiffusive (γ = 1.18 in (Fig. 2E). We also observed high heterogeneity in MSD in mice imaged 1 wk after RAS immunization (Fig. 2E). At this time point after in vivo priming, many CD8 T cells do not undertake the crawling behavior characteristic of liver-localized T cells, most likely due to intermediate expression of LFA-1 shortly after immunization (31).

Likewise, all but one cell population had relatively small movement lengths. The exception was for CD8 T cells in mice 1 wk after RAS immunization (Fig. 2F), when cells do not attach firmly to the liver endothelium and often flow in the bloodstream, thus performing Lévy flights (μ < 3; (Fig. 2F and [31]). To detect such floating events, we imaged the livers at higher than usual frequency (3–6 frames/s). Removing displacements that occur at high speeds (e.g., ≥40 μm/min) resulted in movement length distribution being consistent with thin-tailed distribution (i.e., Brownian-like movements) and yet transiently superdiffusive MSD. Interestingly, in our more recent experiments, some i.v. transferred naive CD8 T cells were also able to localize in the liver; such cells displayed superdiffusive displacement primarily due to floating with the blood flow (O’Connor, J. H., H. A. McNamara, Y. Cai, L. A. Coupland, E. E. Gardiner, C. R. Parish, B. J. McMorran, V.V. Ganusov, and I. A. Cockburn. 2021. Platelets are dispensable for the ability of CD8+ T cells to accumulate, patrol, kill, and reside in the liver, manuscript posted on bioRxiv, DOI: 10.1101/2021.11.09.467964). Transiently superdiffusive behavior with small, Brownian-like movements was also observed for Plasmodium-specific activated CD8 T cells and activated CD8 T cells with irrelevant specificity in P. berghei–infected mice (Supplemental Fig. 5). Taken together, our results establish that most crawling liver-localized CD8 T cells display correlated random walks, resulting in transiently superdiffusive behavior and small movements (46). T cells that poorly attach to the sinusoidal endothelium undergo a combination of correlated random walk with Lévy flights also resulting in superdiffusion and fat-tailed distribution of movement lengths.

Given the network structure of liver sinusoids (Fig. 1B) and the universal pattern of transient but relatively long superdiffusive behavior of liver-localized CD8 T cells, we hypothesized that a linear structure of the sinusoids may be forcing T cells to exhibit persistent walks. Therefore, we analyzed the 3D structure of liver sinusoids. For these intravital imaging experiments, we injected a contrast dye (Evans blue or rhodamine-dextran) and imaged murine livers (Fig. 3A). Z-stacks of the images were processed in MATLAB to generate highlighted sinusoidal maps (Fig. 3B). A subset of the full image was further processed to generate nodes and edges of the sample (Fig. 3C, 3D). Analysis of the processed digital structures suggested that the branch length distribution was best described by the half-normal distribution, with GP distribution being second best (Fig. 3E, Supplemental Tables I, II). The longest branch detected did not exceed 100 μm, and the tail of the branch length distribution was characterized by a relatively high shape parameter of the power law distribution (μ = 7.65; (Fig. 3F) (34). Importantly, sinusoids tend to branch at a small angle, which is likely to improve blood flow (Fig. 3G) (47). Finally, analysis of the neighboring nodes in the network (Fig. 3D) suggests that sinusoids connecting two adjacent nodes tend to be straight (Fig. 3G). Taken together, properties of the liver sinusoids such as small branching angles, straightness, and shortness appear to match the observed pattern of movement of liver-localized T cells (e.g., Figs. 2F, 3F).

To investigate how the structure of the liver sinusoids may impact the movement patterns of cells with an intrinsic movement program, we simulated movement of cells assuming that they are either Brownian or Lévy walkers (Fig. 4). Specifically, we assumed that the distribution of movement lengths follows a power law distribution, with the shape parameter μ being different between Brownian (μ = 4.5) and Lévy (μ = 2.5) walkers while fixing the average movement length for two walk types, which is possible when μ > 2 (1). As expected in free 3D space, Brownian walkers displayed normal diffusion (γ = 1), whereas Lévy walkers superdiffused (γ = 1.5; (Fig. 4A–4C). In contrast, when cells were run on the digitized network of liver sinusoids (Fig. 3A–3D, Supplemental Fig. 6) with cells allowed to make a random choice of turning when reaching a branching point, the movement pattern changed dramatically (compare (Fig. 4A with (Fig. 4D and (Fig. 4B with (Fig. 4E). Both Brownian and Lévy walkers displayed superdiffusion in the liver structure for a substantial period of time (∼55 time steps; (Fig. 4F), with the diffusion becoming normal at longer time scales. This pattern is similar to that of small particles displaying transient superdiffusion early and normal diffusion later due to motion of self-propelling bacteria in media (48). Such a long transient superdiffusion was not due to the assumption of forward movement by Brownian-like walkers, because simulations in free 3D space assuming that Brownian-like walkers can only move forward resulted in superdiffusion for only approximately seven time steps (results not shown). Thus, the structure of the liver sinusoids along with the requirement for T cells to move forward was sufficient to allow transient yet long-term superdiffusion of Brownian-like walkers such as crawling liver-localized CD8 T cells.

By finding that the structure of liver sinusoids may be sufficient to allow transient superdiffusion of T cells, we wondered if the same phenomenon may apply to T cells in other tissues. Because previous studies have suggested that T cells in LNs are close to Brownian walkers (3, 6, 49, 50), we sought to determine if this is because of the structural constraints of the LN environment that does not allow superdiffusion. It has been suggested that T cells use the network of FRCs/DCs as a guide for movement (51, 52). Recently, the FRC network of the LNs has been imaged using confocal microscopy (39). Using a similar computational pipeline, we digitized the FRC network using previously published data (39) and calculated its basic characteristics (Supplemental Fig. 7A–C). Interestingly, similarly to the liver sinusoids, the FRC network had relatively short branch lengths, with the tail of the distribution being well described by the power law distribution with a shape parameter μ = 6.7 (Supplemental Fig. 7D). A half-normal distribution fitted the data on branch length distribution with best quality (Supplemental Fig. 7E, Supplemental Tables I, II). Similarly to the liver sinusoids, the FRC network has relatively small branching angles at the nodes, and connections between adjacent nodes are mostly straight (Supplemental Fig. 7F, 7G). By assuming that T cells can only walk in the middle of the digitized FRC networks with preferred movement, forward simulations of T cell walks suggested that, independently of whether cells do Brownian- (μ = 4.5) or Lévy-like (μ = 2.5) walks, they display transient superdiffusion similar to that of T cells simulated to walk in the liver (Supplemental Fig. 7H). This is perhaps unsurprising, given how similar the derived characteristics of FRC network is to the characteristics of the liver sinusoids. In contrast, analysis of the experimental data on movement of naive CD8 T cells in murine LNs (39) suggested that such T cells are superdiffusive for a very short time period (Supplemental Fig. 8). The duration of a persistent walk of CD8 T cells in LNs is ∼4–5 min, whereas CD8 T cells in the liver superdiffuse for nearly 20 min (Supplemental Fig. 9). Analysis of the actual moving tracks of T cells in LNs suggests that many T cells tend to turn after a short directed movement (Supplemental Fig. 8A, 8E). Indeed, it is well known that naive T cells in LNs tend to be restricted to the T cell areas, most likely by the gradient of chemokines. Also, it is possible that FRC/DC-based scaffolds are wide, thus allowing more freedom of T cells to turn. In contrast, liver-localized T cells are likely to be more constrained by the size of the sinusoids, and, in the absence of any specific liver infection, there are likely not any chemoattractant cues that would force T cells to turn more. Our results thus suggest that although structural constraints of the liver sinusoids or FRC network in LNs may allow long-term superdiffusion of T cells, other factors, such as details of the physical scaffolds used for crawling and/or attraction cues in the tissues, may impair a cell’s ability to superdiffuse.

Determining the movement strategy of T cells in tissues may be important because different walk strategies may have different efficacies (8–11, 13–17, 33, 37, 50, 53–57). Search efficiency may depend on details of a specific biological system, however. Plasmodium-specific, liver-localized CD8 T cells need to locate and eliminate few Plasmodium-infected hepatocytes before these parasites mature and escape into the circulation. We therefore simulated a search by a single liver-localized CD8 T cell for a single Plasmodium liver stage, assuming that the parasite induces no attraction to the infection site and restricting the average movement speed to be the same for Brownian or Lévy walkers. We used Pareto (power law) distribution to describe the movement lengths with scale parameters μ = 4.5 for Brownian walk and μ = 2.5 for Lévy walk (see Materials and Methods for more details). When a single CD8 T cell searches for the parasite in constraint-free 3D space, on average, Lévy walkers require less time (⁠$T¯$⁠) to locate the target (Fig. 5A, 5C), which is consistent with previous results (8, 9). Lévy walks also result in smaller variability (characterized by the SD of the time to target σ) in the time to find the target. Yet, because of the large variability in the time to target, the difference in the time to target between Brownian and Lévy walkers in these simulations may seem biologically nonsignificant.

Interestingly, when the search is performed in the digitized liver structures, Brownian-like walkers are, on average, more efficient at finding targets than Lévy walkers (Fig. 5B, 5D), with smaller variance in the time to target. Furthermore, in 55% of simulations, Lévy walkers would take longer than the slowest Brownian walker to find the target (Fig. 5D). This may be particularly important in the case of malaria, where a failure to find all the targets in the liver within a fixed time period (e.g., 48 h) will result in blood stage infection and fulminant disease. Alternatively, in the liver, Brownian walkers may be unable to locate the infection within a given time period if they start their search too far away from the parasite, whereas Lévy walkers have a nonzero chance (Fig. 5D).

The difference in search efficiency between Lévy and Brownian walkers in free 3D space and in the liver is due to differences in the paths to the target. In the digitized liver, there are essentially few short pathways between a T cell and the target, and Brownian walkers make consistently small movements in moving closer to the target. In contrast, due to small average movement length, Lévy walkers spend long periods of time pausing before performing a long run, sometimes resulting in them reaching the target rapidly. However, long pauses also result in much longer overall time to the target, which increases the average search time substantially.

In these simulations so far, we assumed that a single T cell searches for the single infection site. However, in reality, multiple T cells may be searching for the same parasite. Indeed, our recent analysis suggests that for ∼10⁵ liver-localized CD8 T cells in a mouse, at least 1 T cell should be within 100 μm of the parasite with 99% certainty (38). In our experiments, we rarely see parasites surrounded by more than 5–8 T cells, suggesting that, in many situations, <10 T cells are close enough to the parasite to find it within several hours (42, 43). As expected, the average time to find the target declines with an increasing number of searchers, independently of whether the search is done in free 3D space or in a constrained liver environment (Fig. 5E, 5F). Intuitively, the reduction in time to find the target arises due to variability in search time because we assume that search by individual T cells is done independently; therefore, with many T cells searching for infection, the one T cell having shortest search time will define the success of the search, and Lévy walk allows a larger fraction of faster searches. Yet, for physiologically relevant numbers of T cells searching for the infection site in the liver (≤10/parasite), we find a relatively small difference in the search efficiency between Brownian and Lévy walkers (Fig. 5E, 5F).

T cells use various movement strategies, depending on the tissue they are in. By combining experiments, data analysis, and computational modeling, we have shown that T cells in the liver undergo two types of walk: correlated random walk (by T cells that crawl on sinusoidal endothelium) and a combination of correlated random walk and Lévy flights (by T cells, e.g., 1 wk after RAS immunization, that poorly attach to sinusoidal endothelium and float with the blood flow; (Fig. 2F). Interestingly, even crawling T cells display transient yet relatively long-term superdiffusion, which we attribute to the special properties of liver sinusoids that constrain movement of liver-localized intravascular CD8 T cells (Fig. 4). Our results thus suggest that observed patterns of superdiffusion or rare long-range displacements observed in some experimental data need not arise from cell-intrinsic movement programs but can be explained by the physical and/or physiological properties of the environment in which T cells move (e.g., by floating with the blood flow). Interestingly, although we found that LFA-1–deficient liver-localized CD8 T cells do not firmly attach to sinusoidal endothelium and often flow with the blood flow (31), we did not detect Lévy flights for such data (μ > 3 in (Fig. 2F). This was due to a lower frequency of imaging in these experiments than the imaging frequency in experiments in week 1 after RAS vaccination (27 s/frame versus 6.5 s/frame, respectively; see Materials and Methods for more details). Therefore, the frequency at which imaging is performed can influence the inference of the type of movements cells seem to exhibit.

Specifics of the physical structures in tissues and other factors, such as chemoattractant cues, are likely to contribute to the movement patterns of T cells in tissues. In particular, we showed that the FRC network of the LNs also allows transient superdiffusion of agents in simulations; however, actual T cells in LNs superdiffuse for very short time periods (Supplemental Figs. 8, 9). Thus, other factors such as chemical cues (chemokines), finer details of the physical paths on which T cells walk such as FRCs/DCs, infection-induced inflammation, and preference for forward movement in the absence of external disturbances are likely also to be important in determining overall patterns of T cell movement in tissues. In particular, the role of environmental cues in determining agent movement patterns has been suggested in several systems (3, 22, 58). In the liver, other cells, such as platelets and APCs, have been shown to contribute to movement patterns of liver-localized CD8 T cells (30, 59).

Our results have broad implications for other systems in which superdiffusion and rare long-range displacements have been observed or suspected (11, 60). It will be important to reexamine the evidence supporting intrinsic search programs by matching movement patterns of agents to the experimentally measured environmental details (e.g., airflows for albatrosses or currents for marine predators) and evaluating the role of cues (chemical or otherwise) in defining overall movement patterns of agents (13, 14, 37). Ultimately, future studies will need to quantify the relative contribution of intrinsic programs, environmental physical constraints, and cues in determining the movement pattern of agents (20).

Given our results, it is tempting to speculate about potential mechanisms resulting in Lévy walks of brain-localized CD8 T cells observed previously (8). First, one possibility is that Lévy-like displacements arose due to special structures in the brain used by T cells for crawling. Such structures will then have to be poised at the threshold of percolation (61, 62). Second, another possibility is that some brain-localized CD8 T cells are in fact located in brain vasculature and thus undergo Lévy flights with the blood flow, similarly to what we have observed for liver-localized CD8 T cells under some circumstances. Third, T cells may be moving faster in places where other T cells had left pathways (e.g., as has been proposed for T cells moving in collagen-fibrin gels) (5). Finally, Lévy walks may arise in the analysis due to data generation or processing artifacts. We have found that when many T cells are present in the tissue, it is sometimes difficult to segment experimental data to identify individual cells and their movements. We have also found in our and other publicly available datasets that cell coordinate data often have missing time frames, and, if such tracks are not cleaned, data analysis may generate an impression of a long displacement between two assumed to be sequential time frames. Whether such artifacts contributed to the overall detection of Lévy walks of T cells in the brain remains to be determined, and in general it can be beneficial when raw data from imaging data analysis are made publicly available following publication (e.g., the data of Harris et al. [8] are not publicly available).

Our study has several limitations. Only a small area of the liver (∼500 × 500 × 50 μm) is generally imaged with intravital microscopy, which restricts the number of cells recorded and creates a bias for retention of slowly moving cells in longer videos. Although we show that the digitized structure of liver sinusoids can override the intrinsic movement program of T cells, we did not estimate how much of the experimentally observed movement pattern of liver-localized T cells is due to a cell-intrinsic program, such as a rare ability of cells to turn back, physical constraints of the liver sinusoids, or other environmental cues. (We did, however, find that even if Brownian-like walkers are not able to turn back in free 3D space, they superdiffuse for a relatively short time.) Methods for discriminating between these alternatives will have to be developed. Finally, we have focused on explaining our data with the simplest possible movement models involving movement length distribution, MSD, and turning angles. Other, more complicated Lévy walk models have been suggested, such as generalized Lévy walks (GLWs) or Zigzag-GLW, that allow runs and pauses and assume or generate heterogeneity in cell velocities (8, 33, 45). However, such models in general have more parameters, and fitting GLW-like models to our data are likely to result in overfitting. We could explain our data with relatively simple movement models such as correlated random walks or a combination of correlated random walks with Lévy flights; Lévy flights due to T cell floating in the blood were responsible for the fat-tailed distribution of movement lengths in our week 1 post-RAS immunization data. Our results suggest that future studies focusing on benefits of evolutionary optimality of Lévy walks would gain from more thorough investigation of the structural constraints and physiological tissue details that may be responsible for observed Lévy-like movement patterns.

Our analyses raise many additional questions for future research. In particular, the impact of the Plasmodium infection on movement patterns of liver-localized CD8 T cells remains to be explored further. Our analysis of one unpublished dataset suggests similar long-term superdiffusion and short displacements for CD8 T cells following the infection (Supplemental Fig. 5); however, this remains to be more fully investigated in the future. In particular, our recent work suggests that although, initially, Plasmodium-specific CD8 T cells search for the infection site randomly, following formation of a T cell cluster near the parasite, other T cells become moderately attracted to the cluster (38, 44). Whether movement patterns of liver-localized T cells change following formation of T cell clusters near the parasite remains to be investigated. Furthermore, infection may also influence blood flow, allowing T cells in circulation to be recruited to the infection site. Therefore, it may be advantageous for T cells that are crawling on endothelial cells of the liver sinusoids to detach from the endothelium and rapidly float with the blood, thus undergoing Lévy flights. Whether a combination of local search by crawling and then large displacements with the blood flow allows more efficient localization of all parasites in the liver also remains to be determined.

We thank Viktor Zenkov for contributing to the start of this work by manually digitizing liver sinusoidal maps and Marco Novkovic for providing imaging data of lymph nodes and data on movements of naive CD8 T cells in murine lymph nodes.

The authors have no conflicts of interest to disclose.