Resolving mixed mechanisms of protein subdiffusion at the T cell plasma membrane

The plasma membrane is a complex medium where transmembrane proteins diffuse and interact to facilitate cell function. Membrane protein mobility is affected by multiple mechanisms, including crowding, trapping, medium elasticity and structure, thus limiting our ability to distinguish them in intact cells. Here we characterize the mobility and organization of a short transmembrane protein at the plasma membrane of live T cells, using single particle tracking and photoactivated-localization microscopy. Protein mobility is highly heterogeneous, subdiffusive and ergodic-like. Using mobility characteristics, we segment individual trajectories into subpopulations with distinct Gaussian step-size distributions. Particles of low-to-medium mobility consist of clusters, diffusing in a viscoelastic and fractal-like medium and are enriched at the centre of the cell footprint. Particles of high mobility undergo weak confinement and are more evenly distributed. This study presents a methodological approach to resolve simultaneous mixed subdiffusion mechanisms acting on polydispersed samples and complex media such as cell membranes.

T he plasma membrane (PM) of cells is a diverse, multicomponent complex medium through which the cell interacts with its surroundings. Proteins at the PM diffuse and interact to facilitate a wide range of cellular functions, including sensing and signalling 1 . Specifically, T cells probe the surface of antigen-presenting cells (APCs) for cognate antigens to trigger an adaptive immune response. Antigen recognition is achieved by highly specific T-cell antigen receptors (TCRs) and leads to the rapid development of a complex interface between the cells-the immune synapse. T-cell activation results in dramatic macroscopic rearrangement of protein distribution at the immune synapse 2,3 . However, much remains to be learned about the microscopic properties of the PM and protein mobility within 4 .
According to the Nicholson-Singer model, the PM can be regarded as a complex fluid, in which transmembrane proteins diffuse laterally 5 . The random motion of particles in a purely viscous and homogeneous fluid is known as Brownian motion and is characterized in two dimensions by: where the left-hand side is the mean squared displacement (MSD) of the particle from its origin, P(r,t) is the probability distribution function (PDF; or propagator) of the diffusion process, D is the diffusion constant and t is time. The MSD is typically measured in two ways. The first, as an average across an ensemble of particles (eMSD): wherer i t ð Þ is the location of particle i at time t and N is the number of particles in the ensemble. The second method, as a function of gap time for a single particle (tMSD): where T is the number of frames in the trajectory, t is the time gap measured in frames and Dt is the measurement time step. One can further take the mean of the tMSD functions of multiple trajectories to obtain an average tMSD of the ensemble: The ensemble and mean time averages converge to the same value for large N and t for an ergodic system. Complex media may lead to sublinearity of the MSD as a function of time [6][7][8][9] , that is, where 0rao1. Such a motion is said to be subdiffusive. K a is the generalized diffusion coefficient with dimensions of [m 2 s À a ]. Note should be taken that the mean of the tMSD functions in equation (4) is an arithmetic mean. In contrast, a geometric mean provides an accurate estimator for the mean value of the power a (ref. 10): The geometric mean is mathematically smaller or equal to the arithmetic mean. Because of the differences in averaging, the resultant K a values for the time averaged measurements in equation (6) are biased to be smaller than the K a values of the ensemble averages in equation (2). Several mechanisms may give rise to subdiffusive motion. These are often described by related mathematical models, including [11][12][13] : (a) Diffusion of tracer particles in a viscoelastic medium statistically results in anti-persistent motion and can be described using the fractional Brownian motion (fBM) model 14 ; (b) Tracer particles may experience trapping by specific interactions with other particles or objects in the medium. The particles may exhibit trapping events with a heavy-tailed waiting time distribution. Such motion can be described using the continuous time random walk (CTRW) model 15 ; (c) Tracer particles diffusing in obstructed or labyrinth-like environments demonstrate movement in a fractal-like space with a dimension d f smaller than the real space dimension. Such movement is modelled by a random walk on a fractal (RWF) 16 ; (d) Tracers diffusing in a confined environment due to non-permeable physical boundaries demonstrate normal diffusion within the boundaries at short timescales, appear to be subdiffusive in intermediate timescales, and will saturate to a flat MSD at long timescales. In the case of permeable boundaries, the MSD will regain normal diffusion at long timescales 17 .
A set of tests has been proposed to determine the dominating mechanism underlying an observed subdiffusive process 11,18 . Naturally, multiple mechanisms of subdiffusion may act simultaneously, which may dramatically complicate the analysis. For instance, cases of mixed mechanisms have been demonstrated and modelled by a combination of CTRW and RWF processes 19,20 . Moreover, particle mobility may be also complicated by static heterogeneity (for example, particle polydispersity) or spatial inhomogeneity of the medium. Effective or unified models have been developed to address heterogeneity in single particle tracking (SPT) results [21][22][23][24][25][26] . However, cases of mixed underlying mechanisms of subdiffusion or of static heterogeneity remain poorly understood as it is often unclear how to distinguish between the underlying mechanisms.
Here we characterize the mobility and organization of a short transmembrane protein at the PM of live T cells, using SPT 27 and photoactivated-localization microscopy (PALM). This transmembrane segment shows highly heterogeneous subdiffusive mobility consistent with ergodic underlying mechanisms. We use a segmentation approach combined with multiple statistical tests 11,[28][29][30] to discriminate and characterize subdiffusive behaviour, including fractal-like and weak confinement. This method is able to resolve simultaneous mixed subdiffusion mechanisms acting on cell membranes.

Results
SPT of a short transmembrane protein. The tracers we use are part of the human immunodeficiency virus envelope glycoprotein gp41, whose assembly at the PM is critical for viral budding 31 . Our construct consists of the transmembrane domain of gp41 with a mutation that inhibits its binding to cholesterol 32 and thus, to cholesterol-enriched domains at the PM. The gp41 constructs are tagged with a fluorescent protein (PAGFP), which is specifically labelled with a primary antibody stained with the fluorescent dye Alexa Fluor 594 (AF594). Using total internal reflection (TIRF) microscopy, we visualize and localize proteincontaining particles with a resolution of B25 nm at a frame rate of 85 fps ( Supplementary Fig. 1a). Importantly, our constructs could be imaged for relatively long durations ( Supplementary  Fig. 1b) at a high signal-to-background ratio using either AF594 fluorescence for SPT or PAGFP photoactivated emission in PALM mode 33 (see 'Methods' section). Figure 1a shows a representative frame of an imaged live T cell as it adheres to a functional TCR-stimulating coverslip. The coverslip is coated with an aCD3e antibody that crosslinks the TCRs and leads to robust cell activation and spreading. The bright spots are particles highlighted by AF594 at the PM of the cell. The inset shows the bright field image of the cell. Imaging of cells over 59 s (5,000 frames) yields movies in which particles move across the cell footprint. Next, we use a tracking algorithm (modified from refs 34,35) to track single particles for durations of up to 5,000 frames before fluorescence photobleaching of the sample. Trajectories shorter than 50 frames are discarded in the analysis to increase the robustness of the data. Figure 1b,c shows a representative cell with trajectories coloured by the values of their individual subdiffusive power a. A wide range of particle mobility is detected. Figure 1d,e shows a PALM image of a representative cell.
Protein mobility is heterogeneous and ergodic-like. We next aim to study the mobility of individual particles and their underlying subdiffusion mechanisms. However, we notice a large fraction of immobile trajectories that may originate from nonspecific antibody staining, background or completely immobile particles. These trajectories may complicate our understanding of the system due to their dominating abundance relative to the rest of the mobile trajectories. Thus, we start our analysis by removing completely immobile trajectories. Previous studies achieved this by placing a threshold on the diffusion coefficient 22 , or on the radius of gyration R g of the trajectory 21 , which is calculated as: wherer i is the coordinate at time step i out of N total time steps. We consider the case of an immobile particle in which the course of the apparent trajectory is governed solely by the localization error of the system 36,37 . In two dimensions, a particle located at the origin would result in localizations distributed normally with mean 0 and variance s 2 . In such a case, R 2 g and the mean step size h|Dr|i are 2s 2 and s ffiffiffi p p , respectively (see 'Methods' section). Thus, the ratio R g = Dr is a constant independent of the localization error. Therefore, we use here both R g and h|Dr|i of the trajectory to robustly remove immobile particles. We validate our approach by comparing two controls (immobile gold beads and gp41 proteins in a fixed sample) and compare their results to data acquired in live cells (Fig. 2a). The normalized ratio ffiffiffiffiffiffiffi ffi p=2 p R g = Dr j j h i of the controls serves to set a threshold for robustly excluding immobile particles (66% of all detected particles in the live cells; see 'Methods' section).
We now turn to the characterization of the mobile particles through their time and ensemble MSDs (Fig. 2b). From the MSDs we calculate a, which is the slope of the curve in logarithmic scale, and K a , which is its value at t ¼ 1. We find that on average, the mobility of particles is subdiffusive (a t ¼ 0.72 ± 0.04; a e ¼ 0.75 ± 0.01). As expected from the differences between the arithmetic and geometric averaging, we obtain K a,t ¼ 1.09(1) Â 10 À 13 m 2 s À a , which is smaller than K a,e ¼ 3.4(1) Â 10 À 13 m 2 s À a . To account for this difference, we simulate MSDs with comparable characteristics to the experimental data and find similar differences in K a values between the arithmetic and geometric means ( Supplementary Fig. 2). The small difference in the experimental a values of the time and ensemble averages (being significantly different from 1.0) are consistent with ergodic underlying mechanisms of subdiffusion. We suggest later in the discussion that this small difference is due to weak confinement of the particles. Furthermore, our data shows high heterogeneity of a values from 0 to B1.5, and of K a values, which span over B4 orders of magnitude (Fig. 2c,d).
Resolving heterogeneity by trajectory segmentation. Closer examination of single particle trajectories reveals heterogeneity on the individual trajectory level (dynamic heterogeneity; Fig. 3) and on the ensemble level (static heterogeneity) (Fig. 2c,d). Figure 3a shows representative trajectories that display dynamic alterations in the diffusion state of the tracer particle (highlighted by arrows). Previous studies have shown or suggested such dynamic heterogeneity in the motion of transmembrane proteins in which particles can alternate their mobility state during their motion 21,38,39 .
To better understand the underlying mechanisms of subdiffusion and to resolve possible dynamic heterogeneity, we perform segmentation of individual trajectories into separate parts that correspond to different mobility states {s i }. For convenience, we order the states in this set in an ascending order according to their mobility characteristics (that is, state i þ 1 is more mobile than state i). Note that this ordering assumes no restriction on the possible transition of a particle from one state to any other state.
Multiple approaches have been suggested for segmentation of particle trajectories, often by the identification of transitions between states 26,40 . Trajectory segments with similar characteristics are then assigned to the same state. These segments can be grouped together into subpopulations for further study of their underlying mechanism. Since there are generally no 'closed-form' guidelines for performing optimal segmentation, we use two complementary approaches that highlight different aspects of the data, as explained below.
The first method is a variational Bayesian treatment of a hidden Markov model (vbSPT) 25 . Briefly, this method uses a maximum-evidence criterion to determine the underlying parameters and the number of diffusive states from the observed data. This method can robustly resolve multiple underlying diffusion states from a large number of short trajectories. However, it assumes only pure Brownian motion for each state.
As a second segmentation approach, we adapt concepts used in 'First passage time' or 'escape radius' statistics 19,41,42 , and introduce a modified method that we term 'consecutive escape radii' (CER). Through this method, we measure the particle's ability to travel a given distance in a given time. For each time point in the trajectory, we test if the particle exceeds a distance R s i ;s i þ 1 th during its next n s i ;s i þ 1 th time steps. The writing s i ,s i þ 1 above a variable of interest denotes the threshold value for discriminating between states i and i þ 1 in the ordered set of states {s i }. Figure 3b demonstrates this segmentation method for three states {s 1 ,s 2 ,s 3 } and thus two sets of thresholds R 1;2 th ; R 2;3 th È É and n 1;2 th ; n 2;3 th È É . We show the three different cases for a given time point in a trajectory, classified by their mobility from left to right as s 1 , s 2 and s 3 .
For correct segmentation, one needs a measure to assess the performance of specific segmentation parameters (for example, the set of thresholds {R th } and {n th }) according to their predictive power. For this, we turn to the property of Gaussianity of step-size distributions. Random walks originating from a distinct process effectively yield a Gaussian distribution of step sizes. The Gaussian distribution appears after many steps of the random walk by virtue of the central limit theorem 43 , regardless of the original distribution of step sizes of the underlying process. However, the step-size distribution of our measured trajectories results in a highly non-Gaussian distribution, as evident from the deviation of the step-size distribution (blue circles) from the best Gaussian fit (red line) in Fig. 3c. Step-size distributions of SPT data are often non-Gaussian, which may arise from tracer or environmental heterogeneity 21,24,[44][45][46] .
Step sizes may be inherently non-Gaussian in cases of fluctuating environments on similar timescales as the random walk itself 46 . Still, at long timescales even such cases would appear as Gaussian 44,46 . Thus, regardless of the method by which we segment trajectories, we expect that the detected distribution of step sizes should be divided into multiple Gaussian components, each representing a separate underlying random-walk process. Therefore, an efficient segmentation process can be evaluated based on the extent of data it is able to account for using a minimal number of such components. Note that in vbSPT, the maximum evidence criterion does not consider the statistics of the step-size distribution for the robustness of the state. With this consideration, we scan a large range of threshold sets {R th } and {n th } for the CER method. We then test for each set the extent of step-size Gaussianity and mixing of the resulting subpopulations recovered by the segmentation process. Specifically, we fit a finite mixture of k Gaussian components to each subpopulation, assumed by a model with an arbitrary number of subpopulations k: where P j DrÞ ð is the step-size distribution for subpopulation j, a j i is the weight of each Gaussian component i in subpopulation j and s 2 i is the variance of each Gaussian component. Thus, a i j is a matrix element describing the weight of the ith Gaussian component of the jth subpopulation. According to our model, we require that P k i¼1 a j i ¼ 1 and that 0 a j i 1. We also constrain each s i to be constant over all k subpopulation fits since we assume that the underlying Gaussian components are constant and only their weights in each subpopulation (a j 0 i s) are changing. This assumption is justified since we identify each Gaussian component with a corresponding, distinguishable physical mechanism. The segmentation efficiency is measured by minimization of the off-diagonal a It is a priori unclear how many states are present in the data (or moreover, if states can be regarded as discrete in the first place). Yet, any segmentation effort needs a model with a certain number of states to account for. Thus, we consider models that grow in their number of states, and thus in their complexity. These competing models should be compared according to their ability to account for the data. However, simpler models with the same (b) A schematic representation of the CER segmentation process for a model of three mobility states. Blue and green dashed circles and blue and green dots are the escape radii thresholds R s1;s2 th and R s2;s3 th and the time steps thresholds n s1;s2 th and n s2;s3 th , which distinguish between states 1,2 and 2,3, respectively. In this example, we chose n s1;s2 th ¼ 9 and n s2;s3 th ¼ 11. The leftmost example shows a case where the particle did not escape from R s1;s2 th in n s1;s2 th time steps and thus the tested time point is classified as s 1 , and accordingly, the middle and rightmost examples are classified as s 2 and s 3 . (c) Step-size distribution of 6,410 mobile trajectories 50-5,000 time steps long from 30 cells. The best Gaussian fit is presented as a red line. (d) Step-size distribution of segmented populations of three states, separated by their high, medium and low mobility (red, green and blue circles, respectively). Solid lines are best Gaussian fits of the segmented subpopulations. (e) The representative trajectories from a after segmentation. Segments are coloured blue, green and red according to their mobility states s 1 ,s 2 ,s 3 . Scale bar, 1 mm.
predictive power are superior to more complex ones due to Occam's razor. States that converge to similar underlying mechanisms and mobility characteristics may be assumed redundant and thus, make the suggested model more robust. Hence, a useful approach is to start with a low number of states, and then increase their number until redundancy of one of the states is achieved. Note that this redundancy can be inferred only after a thorough characterization of the mobility states and their underlying mechanisms.
Next, we consider models with an increasing number of states from two to four. Our analyses show that the three states model accounts best for the data and its underlying subdiffusion mechanisms. We first describe the complete set of statistics for the three states model and, only then, fully discuss the inadequacy of the two and four states models at the end of the text. Briefly, using a three states model, our analysis identifies three distinct subpopulations, which we term as 'low', 'medium' and 'high' mobilities. In contrast, the two states model misses a subpopulation with the distinct mobility properties of the 'high' mobility state and its unique underlying mechanism. With the four states model, we find that the 'high' mobility subpopulation has the same characteristics of the 'high' mobility subpopulation found for the three states model. However, the 'medium' and 'low' mobility states become redundant and separate into three subpopulations with essentially the same underlying mechanisms and overlapping mobility characteristics (see 'Discussion' section).
Returning to the three states model, Fig. 3d shows the step-size distributions for each of the subpopulations detected by the CER method. Each single Gaussian accounts for the majority (495%) of its fitted data. Figure 3e presents again the representative trajectories shown in Fig. 3a after segmentation into three mobility states.
At this point, we would like to clarify the relation between different terms used in the text that are invoked through our segmentation and classification processes. Through these processes, segment subpopulations are uniquely assigned to specific mobility states. Thus, we will now interchangeably use the terms 'mobility states', 'mobility subpopulations' or simply 'subpopulations' to describe the segments of common mobility characteristics (being either 'low', 'medium' or 'high'). Moreover, we note that through the segmentation process, a trajectory of a single particle may be divided into several segments that match different mobility states and may thus be assigned to several segment subpopulations.
Characteristics of the different subpopulations. Once we segment all trajectories into three subpopulations, we begin our analyses of each subpopulation to determine its underlying subdiffusion mechanism. Starting with general characteristics of the subpopulations, Fig. 4a shows the normalized distributions of fluorescence intensity of the three subpopulations. This is not a direct characterization of single particles, but rather of their segmented trajectories that might belong to different subpopulations. However, for physically stable particles (that is, they do not disintegrate or assemble over the typical trajectory duration ( Supplementary Fig. 1c)), we can infer on average the particle fluorescence intensity distributions from this analysis. Also, although the fluorescence intensity is not a precise measure of the number of proteins in individual particles, it can give a rough estimate of this number. This lack of precision is due to the possible existence of multiple emitters per tracer. For instance, single antibodies may carry multiple fluorophores (B3 in our case; see 'Methods' section), multiple antibodies may target an individual protein or gp41 proteins may self-cluster 32,47,48 .
In cases of multiple targeting of antibodies and self-clustering, the size of the tracer increases and thus its mobility decreases. We find that segments belonging to brighter particles are of a lower mobility state on average (Fig. 4a).
To better characterize the size distribution of particles, we turn to PALM, which is enabled by the unique construct we use for SPT. In contrast to immunofluorescence, PALM can visualize individually all of the genetically encoded gp41 proteins in each cluster with minimal artifacts and with a resolution down to B25 nm (ref. 49). The detected molecules can then be clustered by assigning proximal molecules beneath a specific distance threshold to the same cluster. We chose a distance threshold of 40 nm that accounts for the resolution limit in the localization of two proximal molecules (see 'Methods' section for further details). Figure 4b shows the cluster size distribution of molecules from a sample cell in logarithmic scale. We compare this size distribution to the one from a random scatter (that is, a Poisson process) of points in the given area. We find that the number of molecules per cluster is skewed towards clusters as small as dimers and trimers, but has a long tail compared to the Poisson distribution. This distribution indicates the existence of larger clusters than expected by a random scatter and high heterogeneity in their size. Figure 4c,d shows the spatial distribution of trajectory segments associated with the three mobility states. We find that the higher the mobility of a segment, the farther its corresponding particle can distribute from the centre of the cell footprint on average. Thus, lower mobility states tend to be correlated with aggregation of proteins and tend to be enriched at the centre of the cell footprint. These observations strengthen the validity of our segmentation process since they capture independent characteristics of the tracer particles without a priori assumptions.
We now return to the MSD statistics for further characterization of the subpopulations. Figure 5 shows the ensemble and time averages of the trajectory segments of subpopulations associated with each of the mobility states. The intersection of the MSD curves with the y axis shows that the three mobility states are highly separable and have variable diffusion coefficient values, K a . Note that our segmentation process excludes segments shorter than 50 frames. Thus trajectories with non-distinct mobility states may be depleted during this process and result in more separable K a values. The subdiffusive powers a and the diffusion coefficients K a obtained for each of the three mobility states are shown in Table 1. We find that all three mobility states are subdiffusive with non-negligible differences between the time and ensemble values. We account for these differences in the discussion. The low and medium mobility states exhibit a similar range of a values (with somewhat higher a values for the low mobility state), but differ strongly in their diffusion coefficients, K a . The high mobility state is only slightly subdiffusive and is substantially different from the other subpopulations in its mobility characteristics a and K a . Importantly, the vbSPT segmentation method yields three subpopulations with similar, yet less separable, K a and a values compared with the CER method (Supplementary Table 1).
Assigning mechanisms. To further study each subpopulation, we calculate their velocity autocorrelation (VAC) functions: is the particle velocity at time t calculated over a time period d. The VAC is computed separately for each trajectory and then averaged over the subpopulation ensemble. The shape of the VAC can separate well CTRW from confined motion and ergodic processes 11,29,30,50 , but cannot distinguish between fBM and RWF 18 . Figure 6 shows the VAC functions of the three subpopulations for a range of t values and changing time periods d. For both fBM and RWF, we expect to find anti-persistent correlations at initial t values that decay to 0 after some time. These negative correlations should appear also for longer timescales d. For CTRW, we expect to find no correlations at all, since the random walker can alter its direction of movement between trapping events. For a confined random walk, we expect to find no correlations at short timescales d, but negative correlations at longer timescales.
The low and medium mobility states are similar and both exhibit negative correlations that are insensitive to the time periods d over a wide range of values. This feature confirms our initial ruling out of CTRW as a possible underlying mechanism and points to either fBM or RWF as the dominating mechanisms of subdiffusion. In contrast, the high mobility state shows markedly different features than the low and medium mobility subpopulations. At the initial time period d, on the order of our time resolution of 11 ms, it shows no (or negligible) correlations. It then shows growing negative correlations for longer d values. This behaviour, together with the a value being close to 1,  Table 1. 3.0(8) Â 10 À 15 1.7(2) Â 10 À 14 0.9(1) Â 10 À 12 K a,t 2.1(5) Â 10 À 15 1.0(1) Â 10 À 14 0.6(1) Â 10 À 12 K a,MME 6(1) Â 10 À 15 3.4(5) Â 10 À 14 1.2(2) Â 10 À 12 Ka Þ 1=a % 3 ms, which is shorter than our time resolution. Thus, we conclude that for small d values only a small fraction of the particles is affected by the boundary, whereas at growing d values, an increasing fraction of the particles is affected by it. We refer to this as weak confinement of the high mobility subpopulation. Conversely, strong confinement would result in significant negative autocorrelations at the initial time period d. Furthermore, we note that our analyses are limited by the trajectory lifetimes to B0.6 s. This time frame is likely too short for the trajectory segments to regain apparent Brownian MSDs (Fig. 5) and null velocity autocorrelations at long time periods (Fig. 6c).
We now focus on the study of the underlying mechanisms for the low and medium subpopulations, since for these subpopulations, the VAC alone can not distinguish between fBM and RWF 18 . We first turn to PALM imaging, as it enables the spatial mapping of available sites for the particles and, consequently, the characterization of the fractal dimension of the cell footprint (related to RWF). We assume the particles are homogeneously spread over the structure imaged by PALM, and thus is effectively captured by it. We calculate the fractal Hausdorff dimension of the gp41 spatial distribution across the PM using the box-counting method 51 (Fig. 7a,b). We find a fractal dimension d f ¼ 1.84±0.07. Notably, this calculation excludes the lowest scales, where the resolution limit is bounding, and the largest scales, where the footprint of the cell becomes significant. The observed fractal-like structure may not be the only mechanism that affects the particles mobility. Thus, we turn below to two complementary statistical tests of the PDF Gaussianity of the SPT trajectories, namely moment ratios analysis 11,28,52 and the growing sphere analysis 11,28 .
The moment ratio analysis tests the extent of Gaussianity of the regular (MSD) and mean maximal excursion (MME) by calculating the regular and MME moment ratios. Note that earlier we tested the step-size distribution for Gaussianity while here, we test the Gaussianity of the PDF of the diffusion process. We compute the time-averaged MME for each trajectory segment as: where k is the moment and T,t and Dt are defined as in equation (3). We provide the MME values for a and K a in Table 1. Note that the MME values for a are higher than the ensemble and time-averaged values. This is expected for an fBM process, but cannot be regarded as conclusive evidence for it 28 .
The regular and MME moment ratios are defined as: and are used as measures for the PDF Gaussianity. These moment ratios should have distinct asymptotic values (that is, for tc1) for fBM and RWF processes. These asymptotic values depend on the subdiffusive power a and, in the case of RWF, also on the fractal dimension d f (ref. 28). For an fBM process, the regular asymptotic moment ratio converges to the value 2.0. There is no exact analytical solution for the MME moment ratio, but an approximation 28 : Also, in the case of RWF, there is an approximated solution, which can be numerically computed for any choice of d f and a (ref. 28). Figure 7c shows the regular (crosses) and MME (triangles) moment ratios for the low (blue) and medium (green) mobility subpopulations. We mark the expected values for fBM and RWF processes for the low and medium subpopulations for the regular and MME moment ratios (overall, seven lines in the figure). Our calculated moment ratios show asymptotic values, which for the MME moment ratios fit slightly better the RWF model, yet hard to distinguish from fBM. More distinctly, the asymptotic values of the regular moment ratios converge to the RWF expected values. This supports the possibility of an RWF process as an underlying mechanism.
The second statistical test is the growing sphere analysis 11,28 . In this method, we calculate the probability of a particle to be within a growing sphere with a radius increasing as pr 0 t a/2 : where Y(r) is the Heaviside function (that is, 1 for rZ0 and 0 for ro0) and N(t) is the number of particles present at time t. For an fBM process, this ratio should be constant, while for an RWF it should scale as / t aðd À d f Þ=2 (ref. 28), thus showing a power law growth with t. Figure 7d shows the growing sphere analysis of the low and medium subpopulations. We find that the values for the growing sphere analysis are constant, with an exception of the shortest timescales, which is found to be a property of the growing sphere analysis as demonstrated by simulated particles (black dash-dotted line) in Fig. 7d. This indicates fBM as an added underlying mechanism to RWF for the subdiffusive mobility of both subpopulations.
Cell activation effects on subdiffusion mechanisms. T-cell activation results in significant changes in its morphology as it adheres to APCs. It is unclear how such changes affect protein mobility at the cell PM 53 . To study the effect of cell activation, we now turn to results from imaging of cells under non-activating conditions, as the cells adhere to coverslips that do not stimulate the TCR (see 'Methods' section). Figure 8 shows the summary of the different tests for these cells using the three states model. As for the activated cells, we find both dynamic and static heterogeneity in the MSDs, thus we follow the same segmentation procedure. After segmentation, using the CER approach, we observe similar overall characteristics: (1) the segments can be separated into subpopulations of low, medium and high mobility (a e ¼ 0.56 ± 0.08, 0.49 ± 0.04 and 0.92 ± 0.02, respectively).
(2) These subpopulations are distributed differentially in the cell footprint as the high mobility subpopulation is evenly distributed and extend to the periphery of the cell, while the lower mobility subpopulations localize to the centre (Fig. 8a,b). (3) The time and ensemble MSDs are consistent with ergodicity of all underlying processes (Fig. 8c). (4) The VAC indicates fBM or RWF processes as the underlying mechanisms of the low and medium mobility subpopulations. The VAC of the high mobility subpopulation shows again confinement, yet a larger negative correlation at short time periods (Fig. 8d-f). (5) The moment ratios for the low and medium mobility subpopulations correlate better with fBM as their regular moments converge to a value of E2 and the MME ratios converge to their expected values (Fig. 8g). (6) The growing sphere analysis also points to fBM, in agreement with the moment ratio analysis (Fig. 8h). Thus, we conclude that without cell activation, the PM is more viscoelastic and less fractal-like as indicated by the mobility properties of the low and medium subpopulations. Still, the three distinct mobility subpopulations are common to activated and non-activated cells, under our measurement conditions. The fractal dimension is the slope of the fit. (c) Regular (MSD) moment ratios (crosses) and MME moment ratios (triangles) for low (blue) and medium (green) mobility states. The expected values are marked for both fBM (short marks on the right axis) and RWF (long horizontal lines) processes.
The expected values come in pairs for the low (blue) and medium (green) mobility subpopulations according to their respective a e values from

Discussion
In this study, we detect heterogeneity of protein mobility at the PM of activated and non-activated T cells. This heterogeneity appears to be both static and dynamic. We observe polydispersity of particles containing a transmembrane protein derived from gp41 and dynamic heterogeneity in which particles may change their diffusion state during their motion. After segmentation of trajectories into subpopulations with distinct mobility states, we notice a third heterogeneity in the spatial distribution of the segment subpopulations, where more mobile trajectories can distribute farther from the centre of the cell footprint. This complexity, in the form of multiple levels of heterogeneity, is likely a common feature of biological samples and is generally difficult to capture and characterize on the whole ensemble level or with small data sets of particle trajectories. Such a complexity is often simplified by using, artificial and in vitro systems including vesicles and lipid bilayers 38 , model fluids, computer simulations and so on. 7,40,54 . However, multiple studies have characterized macromolecular diffusion within the PM of live cells 7,40,54 . A few recent studies have detected underlying heterogeneity in SPT results and suggested unifying models to explain their findings [21][22][23][24][25][26] . However, these studies typically do not attempt to separate the data into subpopulations with distinct mobility states. Thus, they cannot resolve potentially distinct underlying mechanism acting differentially on different particle subpopulations and during the diffusion of individual particles. The dynamic heterogeneity of single trajectories and the non-Gaussianity of their step-size distribution led us to apply segmentation of trajectories according to local variations in particle mobility. We use two complementary segmentation approaches: vbSPT 25 and our modified CER approach. Alternative segmentation approaches may be applied to SPT data, including temporal sliding window, Bayesian algorithms, HMM or supervised classification 40 . Each approach incorporates dynamic variables that serve for classification. However, the target of classification may be subjective. Here we offer the property of Gaussianity of the step-size distribution as a universal classifying target for effective segmentation of distinct subpopulations. This approach is applicable to a wide range of mixed subdiffusive or Brownian random walks, however, excluding systems of fluctuating environments on similar timescales as the random walk itself 46 . Thus, an optimal classifier for mixed random walks may include a robust classifier (for example, vbSPT) with a consideration of the possibility of anomalous diffusion of the subpopulations and with a requirement for the Gaussianity of their individual step-size distributions.
Notably, our analyses rely on a first robust step of rejecting immobile species based on their radius of gyration R g and mean step size h|Dr|i. Such rejection is important since this immobile fraction often dominates the population of single trajectories 21,22 , and may complicate attempts to distinguish the step-size Gaussianity of distinct subpopulations.
In our segmentation process, we consider models with an increasing number of mobility states, from two to four. A two states model fails to separate the medium and high mobility subpopulations as it fails to capture the VAC properties of the high mobility state, while showing dominating features of the medium mobility state ( Supplementary Fig. 3, compare panels (d,e) to Fig. 6). Next, we consider a three states model. This model accounts for 495% of the data of each subpopulation. As expected, the subpopulations that correspond to the three mobility states differ in their diffusion coefficient K a . The subpopulation corresponding to the high mobility state undergoes weak confinement as its dominating underlying mechanism. Surprisingly, the two subpopulations with the lower mobility states are found to share a common underlying mechanism as they have similar mobility characteristics, including comparable a values, and similar VAC curves, moment ratios, growing sphere analysis and spatial distribution. These statistical tests indicate either RWF, fBM or their combination as the underlying mechanism for subdiffusion. Having two subpopulations of low and medium mobilities that differ predominantly in their K a values, but otherwise share multiple statistical characteristics, strongly suggests that they belong to a common subpopulation of particles. The particles in this subpopulation are affected by the environment through the same mechanism and only differ in their size-dependent mobility. To confirm this conclusion, we test the data using a four states model (Supplementary Fig. 4) and find, yet again, a highly mobile subpopulation with the same characteristics indicating weak confinement. The remaining three lower mobility states have differing K a values, but otherwise, a common underlying mechanism (Supplementary Table 2).
For all subpopulations, we notice non-negligible differences between the mobility values a and K a of the time and ensemble averages (for example, Table 1). This opens the possibility of the involvement of non-ergodic mechanisms (for example, CTRW). A previous study 20 predicts specific dependencies of the ensemble and time MSDs on the subordination of a CTRW on a fractal lattice. This study suggests a relation of a e ¼ a CTRW Â a RWF and a t ¼ 1 À a CTRW þ a CTRW Â a RWF . This implies that a CTRW ¼ 1 þ a e À a t . Since a t oa e in all of our results, this yields erroneous values of a CTRW 41.0. An alternative explanation is the general effect of the weak confinement, which we find for the high mobility subpopulation, on all subpopulations. Indeed, weak confinement can generate inequivalent time and ensemble averages in ergodic systems 55 . To test this explanation, we simulate combinations of weak confinement and either Brownian motion or fBM (see 'Methods' section). We find that in both cases, the confinement leads consistently to a t oa e with a difference of B5-20%, depending on the relative radius of confinement, simulated a, localization uncertainty and MSD fitting accuracy (for example, Supplementary Fig. 6). This explanation also holds for the small observed difference in the a values of the ensemble and time averages before the segmentation process (Fig. 2). Thus, we conclude that weak confinement can sufficiently account for the inequivalence in the a e and a t values in all presented cases. We also observe that consistently, K ða;tÞ oK ða;eÞ . This is understandable since the MSD of the time average is computed as the geometric average of all trajectories, while the MSD of the ensemble is an arithmetic average. Therefore, K ða;tÞ consistently underestimates the true K a .
It is difficult to distinguish between RWF and fBM as the underlying mechanisms of the 'low-medium' mobility subpopulation. The moment ratios test and the PM fractal dimension found by PALM indicate RWF as the dominating mechanism. However, the asymptotic behaviour of the growing sphere analysis (Fig. 7d) indicates a possible fBM process, which cannot be excluded by VAC 11 . Thus, our data suggest the existence of both mechanisms acting simultaneously. Statistical tests have been suggested and applied for resolving the subordination of subdiffusion mechanisms 11,19 . However, these tests have been applied to the separation of ergodic (for example, fBM or RWF) from non-ergodic (for example, CTRW) processes, while clear guidelines and concrete examples for the separation of mixed ergodic processes seem missing.
Both fBM and RWF have been associated with molecular crowding 56,57 . Structural patterning of the PM may additionally cause RWF, while a viscoelastic mesh may give rise to fBM 58,59 . We further show that these mechanisms depend on cell activation as RWF is dominant for activated cells, while fBM is dominant for non-activated cells. Indeed, cell activation has been shown to cause global changes in the morphology of the PM and the cytoskeleton 2,3 . These changes include coalescence of signalling proteins into micro-clusters, polymerization of cytoskeletal proteins and tighter adhesion with the APC and spreading. These changes might alter the diffusive motion of proteins at the PM. For instance, signalling proteins might get trapped in clusters 39 or in actin corrals 58,60 , experience diffusion barriers 61 , while the PM might change its three-dimensional patterning 62 .
Here we use a relatively inert transmembrane protein, namely the gp41 transmembrane domain with a mutation that suppresses its interactions with specific lipids 32 and cellular proteins at the PM. However, specific interactions of this protein with the PM environment cannot be completely ruled out. Naturally, other transmembrane proteins in the cell may undergo significant and specific interactions to facilitate their action, and thus may exhibit different underlying subdiffusive mechanisms than are found here (for example, non-ergodic processes) 40 .
Taken together, our results and statistical tests demonstrate the simultaneous existence of multiple subdiffusion processes, namely confinement, RWF and fBM at the PM of cells. These mechanisms act differentially on particles depending on their size and may change their dominance under varying cell activation conditions. Thus, our results shed new light on the multiple mechanisms that govern protein mobility at the PM of activated and non-activated T cells. More generally, our approach and results demonstrate how one can methodologically resolve mixed mechanisms of subdiffusion that act simultaneously on a polydispersed sample within the complex medium of the cell PM. This approach can also serve as a framework for resolving mixed mechanisms in a wide range of random walks of arbitrary nature.

Methods
Sample preparation. Gp41-PAGFP 63 was generated from previously published constructs of gp41 (D1 mCRAC gp41 (ref. 32)). DNA was introduced into E6.1 Jurkat T cells (a kind gift from the Samelson lab at the NIH) using the Neon electroporation system (Invitrogen). Cells were next incubated at 37°C for 48 h. Since PAGFP is not stable or bright enough for SPT, we used it as a specific target for immunofluorescence labelling. For labelling, 0.5 mg anti-GFP, rabbit polyclonal antibody, Alexa Fluor 594 (AF594) conjugate (Life technologies Cat. A21312 Lot 1476604) was added to 5 Â 10 5 cells suspended in FACS buffer for 30 min on ice. According to vendor specifications, there are three AF594 fluorophores bound to each antibody molecule. The cells were washed three times with PBS and suspended in an imaging buffer. Cells were dropped onto a coated coverslip for live cell imaging. For fixed cell imaging, the cells were first incubated for 20 min and then PFA was added at a ratio of 3/5 for cell fixation over another 30 min (in an incubator) before washing and imaging.
For imaging, we used coverslips (Ibidi m-8-well glass bottom No. 1.5H) coated with 0.01% poly-L-lysine (Sigma) and aCD3e (UCHT1; BD Biosciences) or aCD45 (BD Biosciences) antibodies at a concentration of 10 mg ml À 1 . aCD3e coating robustly stimulates the TCR and leads to T-cell activation and spreading. In contrast, the aCD45 coating does not stimulate the TCR and serve for T-cell adherence to the coverslip without its TCR-dependent activation.
Imaging. Imaging was performed on a TIRF microscope (Nikon) with a CFI Apo TIRF Â 100 oil objective (NA 1.49, WD 0.12 mm). Imaging in TIRF mode served to visualize molecules at the PM of spreading cells in close proximity to the coverslip (up to B100 nm). Cells were first dropped on the coverslip in imaging buffer and given B10 min to spread before imaging commenced. Imaging was conducted for 20-30 cells per experiment at a rate of B3 min between cells, hence B60-90 min typically passed between the first and last cells imaged. Imaging was conducted at room temperature.
For live cell imaging, cells were illuminated by a 561 nm excitation laser at 80% power for 5,000 frames at 85 fps of an EMCCD Ixon þ camera. Fixed cells were further imaged by PALM. For that, they were illuminated using a 405 nm laser for photoactivation of the PAGFP using a changing intensity over the duration of the imaging sequence (typically, using 1-20% power). Illumination with a 488 nm laser at 100% power was used for PAGFP excitation.
Particle localization. Single molecule localization was performed using an ImageJ ThunderStorm plugin 64 . The default values were used for the analysis (B-spline wavelet filter-order 3 and scale 2.0, approximate localization by eight-neighbourhood local maximum, Subpixel localization by PSF integrated Gaussian with the weighted least squares fitting method with a 3 pixels fitting radius and 1.6 pixels initial sigma). Particle coordinates and statistical properties were exported and further analysis was conducted using Matlab 2013a (MathWorks). Image drift and vibrations were corrected by mean image displacement values using a custom code in MATLAB.
SPT. Linking particle locations into trajectories were done using a modified MATLAB version of refs 34,35. Briefly, the algorithm links particles in a given frame with the particles in the previous one. Particles are linked according to their displacement distance between frames. In cases of ambiguity, a choice has to be made. This is done by minimizing the overall distance travelled by all particles between two consecutive frames: P N i¼1 Dr 2 i where N is the number of particles and Dr i is the displacement of each one. The algorithm uses a threshold for the maximum displacement allowed for a particle in a single time step. We use a distance corresponding to three times the expected displacement: Assuming a normal diffusion with D ¼ 0.5 mm 2 s À 1 , at a frame rate of 85 Hz, this corresponds to d max E459 nm. We chose for the algorithm to have 'no memory', that is, if a particle is lost for one frame or more, the linkage is terminated and a new track would begin once it returns. After the linking process is done, only tracks longer than 50 consecutive frames were kept for further analysis. These stringent measures ensure that erroneous linkage of random localizations would be highly suppressed since the probability of persistent linking of randomly occurring points for 50 consecutive frames is negligible at the particle density we work with.
Immobile trajectories threshold. As discussed in the 'Results' section, we set a threshold on the ratio of the radius of gyration R g and the mean step size h|Dr|i. For ideal immobile particles, this ratio is a constant, while for a particle that actually propagates in space the ratio will increase. For an immobile particle, where the apparent motion is due to localization errors, the localizations would be distributed in a 2D Gaussian distribution around the true location, namely, N 0; s 2 ð Þ. The radius of gyration for a normal distribution with variance s 2 is simply R 2 g ¼ 2s 2 , which is just the second moment. For the mean step size, we want to compute the expected value: , where X 1,2 and Y 1,2 are random variables drawn from a normal distribution with mean 0 and variance s 2 . The addition or subtraction of two such values is a normal distribution with twice the variance of the first: X 1 À X 2 $ N 0; 2s 2 ð Þ . Next, the square of a normal distribution is the w 2 distribution with 1 degree of freedom (DOF): The sum of w 2 distributions gives a w 2 distribution with 2 DOF: Þ 2 $ 2s 2 w 2 2 and the square root then gives the w 2 with 2 DOF: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Then, the mean step size, which is the expected value, is Dr . Therefore, for a trajectory composed of points taken from a normal distribution with mean 0 and variance s 2 , the ratio of the radius of gyration and the mean step size is a constant: R g = Dr j j h i¼ ffiffiffiffiffiffiffi ffi 2=p p . Importantly, this ratio is independent of the localization error. Note that the localization error is due to multiple factors 36,37 , including factors that depend on the detection system (for example, dark noise, quantization noise and detector pixel size), the sample (for example, background level), as well as particle-related factors (for example, detected fluorescence intensity, its size, velocity and so on).
For fixed cells and for immobile 40 nm gold beads fixed to the coverslip, we find that 95% of the normalized ratio values ffiffiffiffiffiffiffi ffi p=2 p R g = Dr j j h i fall below the value 2.11 and 99% below 3.89. We use the threshold for 5% false positives (that is, 2.11) for the rest of the analyses of the mobile particles. Importantly, this straight-forward method allows the detection of immobile particles, where R g and h|Dr|i have large and comparable values. In contrast, using only the diffusion coefficient or R g for rejecting immobile particles would falsely classify such particles as mobile. Notably, our excluded immobile particles account for 66% of all detected particles.
Statistical tests and calculations. The values of a and K a are calculated by a linear fit to the first 50 time points in any one of the MSDs. Errors for K a and a are computed by the s.d. of values found between different threshold sets, which met the guidelines of the CER segmentation (see 'Methods' section and Supplementary  Fig. 5 for choice of the threshold set for CER segmentation). For the whole ensemble of trajectories, before segmentation, the values are calculated for all time points for a e and only the linear regime (time points 10-30) for a t . The errors, in this case, are taken from the linear fit while using weights according to the s.e.m. of the values of the MSD curves.
Segmentation is done as described in the main text. Gaps up to three time points of inconsistency in the mobility state are smoothed out. For example, a trajectory with mobility states y 1 1 1 1 1 2 2 1 1 1 1 1 1 y would be smoothed to become all 1s. This is done to minimize noise-related errors in segmentation.
Gaussian fitting for step-size distributions is done based on the least square curve fitting function in MATLAB (Optimization toolbox). The code is adjusted from https://www.mathworks.com/matlabcentral/fileexchange/40613-multiplecurve-fitting-with-common-parameters-using-nlinfit so that boundaries can be introduced. Multiple curves are given simultaneously for fitting with the same set of Gaussian parameters as discussed in the main text.
Segmentation by the vbSPT method was done by a published MATLAB code 25 . Measurement of the fractal dimension on PALM images is done by the box-counting method. MATLAB code was custom written for this.
Localization error corrections for second and fourth moments. As shown by Kepten et al. 10 , localization errors affect the apparent MSD and may result in an underestimation of a and K a by raising the MSD at short times. Here we will repeat the calculation done for the MSD by ref. 10 and extend it further to the fourth moment, which we use for the moment ratio test.
Let x and y be the true locations of the particle, and e x and e y be the corresponding localization errors. Then the observed locations arex ¼ x þ e x andŷ ¼ y þ e y . We further assume he x i ¼ he y i ¼ 0 and he 2 wherer ¼x þỹ. The error in the displacement between two locations x 1 and x 2 due to localization amounts to: Dx ¼x 2 Àx 1 ¼ Dx þ e x;2 À e x;1 À Á . Therefore, the observed mean squared displacement is: x;2 D E À 2e x;1 e x;2 zfflfflfflfflffl ffl}|fflfflfflfflffl ffl{ In two dimensions: . Thus, the correction needed for the MSD is: We now turn to the correction needed for the fourth moment: The last argument on the RHS of equation (17) can be expanded to: and the second argument is: Thus, we arrive at: As noted earlier, e 2 x ¼ e 2 y ¼ 1 2 e 2 r and e 4 x can be realized by: Thus: And equation (20) can be written as: We now note that: and so: And so finally, the correction for the fourth moment is: The values for e 2 r and e 4 r are taken from the uncertainty estimates given by the ThunderSTORM localization algorithm 64 .
Choosing the threshold set for CER segmentation. In the case of a three states model, we use a set of four thresholds: R s1;s2 th ; n s1;s2 th È É ; R s2 ;s3 th ; n s2 ;s3 th È É . This means that we are looking for an optimal set in a four-dimentional space. To simplify this process, we fix two thresholds related to the transition between (for example, s 1 ,s 2 ) and scan a range of relevant values for the other two thresholds (for the same example, between s 2 ,s 3 ). For each set of thresholds, we get three subpopulations and a matrix of weights a j i as defined Sec. IV in the main text. As a reminder, the weight matrix reports on the weight of each of the three Gaussian components fitted for each of the three subpopulations. For clarity, we provide a detailed example of a representative case: This matrix reports that the first subpopulation consists of 95% of the low mobility state s 1 and 5% of the medium mobility state s 2 . The second subpopulation consists of a mixture of 3%, 92% and 5% of mobility states s 1 ,s 2 ,s 3 , respectively. The third subpopulation consists of a mixture of 7% and 93% of s 2 and s 3 , respectively. Note that the elements in each row must sum to 1. On the basis of such a matrix, we compute a score of the extent of mixing between states as follows: In the case of our example, we would have a score of 2.603 out of the possible maximum of 3.0. For a given range of thresholds, we construct a 2D heat map showing the score for each set of thresholds of the variable thresholds, while the two other values are fixed. Importantly, the total number of segments that are accounted for by each set of thresholds may differ. This is due to the fact that only trajectory segments longer than 50 frames are kept. Thus, long trajectories may be excluded from the analysis if they break down into several short segments that are shorter than 50 frames each. We take this restrictive approach to be able to classify these segments properly in the next stages of the analysis where we perform multiple statistical tests for the identification of the underlying mechanisms. Supplementary Fig. 5a,b shows two representative heat maps. Supplementary Fig. 5a  . For each set, we specify the on-diagonal values a j i¼j and the number of trajectory segments in each subpopulation. Having the scores and the number of segments, we now provide general guidelines for the choice of a preferable threshold set: 1. R th should not be too close to the resolution limit. Here in Supplementary  Fig. 5b, we prefer threshold sets where R th Z40 nm. 2. Sets that result in all three values of a j i¼j ! 90 % are preferable. 3. Sets that result in any of the a j i¼j ¼ 100 % should be treated with caution since the fitting algorithm puts a bound on the value which cannot exceed 1.0. Thus, such values may indicate a corrupted fit for which the value would have exceeded 1.0 without the boundary. Hence, we do not consider such sets in our choice. 4. We prefer threshold sets that result in numbers of segments that represent well the distinct subpopulations. For instance, diluting a subpopulation to only its most representative trajectory segments would result in a high fitting score, but with too few segments for further analysis.
On the basis of these guidelines, we now turn to discuss the choice of threshold sets in our data. We observe a tendency of high scores to appear along with a ridge in the heat maps with a negative diagonal trend. Along this ridge, there are several sets that meet the above guidelines. We tested these sets and found that they overall provide very similar results and underlying mechanisms for the subpopulations. We conclude that our results are insensitive to the exact choice of a threshold set. Supplementary Figure 5c shows a and K a values for the various threshold sets tested (grey points) together with the threshold sets that met the above guidelines (squares). As can be seen, a and K a vary in their resultant values. This variability is reduced by a factor of B3-4 when using the guidelines. The s.d. of a is typically B1-2% of the mean value and B17% for K a . This high variability for K a is understandable since it varies over several orders of magnitude. The values we find for a and K a are summarized in Table 1. Furthermore, we note that, as discussed in the main text, the low and medium subpopulations are in fact part of a continuum of mobility states and the effort for distinct classification is somewhat subjective. To conclude, we recommend testing the robustness of threshold sets in the resultant a and K a across multiple sets when employing this method.
Clustering of molecules detected by PALM. Individual PAGFP-tagged proteins were clustered using a custom code in MATLAB for non-hierarchical clustering. Two molecules were assigned to the same cluster if they were closer than a threshold distance. We chose a distance threshold of 40 nm that accounts for the resolution limit in the localization of two proximal molecules in our study (see mode at 28 nm for the localization uncertainty in Supplementary Fig. 1a). We also tested thresholds between 10 and 60 nm and verified the resemblance of their resultant cluster size distributions, and thus, the insensitivity of this distribution to the chosen threshold.
Confinement simulations. Simulations were written in custom MATLAB code. Five hundred particles were simulated for 500 time steps. Confinement was introduced by an un-crossable circular boundary with a radius of 10 units (a.u). Particles were distributed randomly inside the boundary at the initiation of the simulation. For an fBM random walk, the 'wfbm' function was used with an a value of 0.5. The mean step size is B0.9 units. Resulting a values are a e ¼ 0.56±0.01; a t ¼ 0.50 ± 0.01; a MME ¼ 0.63 ± 0.02. For pure Brownian motion, the mean step size is B0.4 and the resulting a values are a e ¼ 0.96 ± 0.08; a t ¼ 0.90 ± 0.03; a MME ¼ 0.98±0.03.
These results demonstrate the effects of confinement on random walks with or without an additional subdiffusive mechanism. Importantly, notice the difference between the time averages and the ensemble or MME average values. The choice of simulating fBM rather than RWF is due to the relative simplicity of simulation and demonstration of results.
Code availability. All custom codes used in this work are freely available at https:// github.com/ShermanLab/SubDiffusion. These codes use MATLAB 2013a (MathWorks).
Data availability. The authors declare that all data that supports the findings of this study are available within the article and its Supplementary Information files and from the corresponding author on reasonable request.