Weak ergodicity breaking of receptor motion in living cells stemming from random diffusivity

Molecular transport in living systems regulates numerous processes underlying biological function. Although many cellular components exhibit anomalous diffusion, only recently has the subdiffusive motion been associated with nonergodic behavior. These findings have stimulated new questions for their implications in statistical mechanics and cell biology. Is nonergodicity a common strategy shared by living systems? Which physical mechanisms generate it? What are its implications for biological function? Here, we use single particle tracking to demonstrate that the motion of DC-SIGN, a receptor with unique pathogen recognition capabilities, reveals nonergodic subdiffusion on living cell membranes. In contrast to previous studies, this behavior is incompatible with transient immobilization and therefore it can not be interpreted according to continuous time random walk theory. We show that the receptor undergoes changes of diffusivity, consistent with the current view of the cell membrane as a highly dynamic and diverse environment. Simulations based on a model of ordinary random walk in complex media quantitatively reproduce all our observations, pointing toward diffusion heterogeneity as the cause of DC-SIGN behavior. By studying different receptor mutants, we further correlate receptor motion to its molecular structure, thus establishing a strong link between nonergodicity and biological function. These results underscore the role of disorder in cell membranes and its connection with function regulation. Due to its generality, our approach offers a framework to interpret anomalous transport in other complex media where dynamic heterogeneity might play a major role, such as those found, e.g., in soft condensed matter, geology and ecology.


I. INTRODUCTION
Cell function relies heavily on the occurrence of biochemical interactions between specific molecules. Encounters between interacting species are mediated by molecular transport within the cellular environment. A fundamental mode of transport for molecules in living cells is represented by diffusion, a motion characterized by random displacements. The quantitative study of diffusion is thus essential for understanding molecular mechanisms underlying cellular function, including target search [1], kinetics of transport-limited reactions [2, 3], trafficking and signaling [4]. These processes take place in complex environments, crowded and compartmentalized by macromolecules and biopolymers. A prototypical example is the plasma membrane, where the interplay of lipids and proteins with cytosolic (e.g., the actin cytoskeleton) and extracellular (e.g., glycans) components generates a highly dynamic and heterogeneous organization [5].
The diffusion of a single molecule j, whose position x j is sampled at N discrete times t i = i∆t, is often characterized by the time-averaged mean-square displacement (T-MSD): which for a Brownian particle scales linearly in the time-lag t lag . Application of fluorescencebased techniques to living cells has evidenced striking deviations from Brownian behavior in the nucleus [6], cytoplasm [7][8][9][10] and plasma membrane [11,12]. Indeed, numerous cellular components show anomalous subdiffusion [13], characterized by a power law dependence of the MSD ∼ t β , with β < 1 [14][15][16]. Owing to the implications of molecular transport for cellular function and the widespread evidence of subdiffusion in biology, major theoretical efforts have been devoted to understand its physical origin. Subdiffusion is generally understood to be the consequence of molecular crowding [17] and several models have been developed to capture its main features. In general, subdiffusion can be obtained by models of energetic and/or geometric disorder, such as: (i) the continuous-time random walk (CTRW), i.e., a walk with waiting times between steps drawn from a power law distribution [18]; (ii) fractional Brownian motion, i.e., a process with correlated increments [19]; (iii) obstructed diffusion, i.e., a walk on a percolation cluster or a fractal [15]; (iv) diffusion in a spatially and/or temporally heterogeneous medium [20][21][22]. Some of these models have been associated with relevant biophysical mechanisms such as trapping [23], the viscoelastic properties of the environment [24] or the presence of barriers and obstacles to diffusion [25].
Advances in single particle tracking (SPT) techniques have allowed the recording of long single-molecule trajectories and have revealed very complex diffusion patterns in living cell systems [5,11]. Recently, it has been shown that some cellular components show subdiffusion associated with weak ergodicity breaking (wEB) [9,10,12], with the most obvious signatures being the non-equivalence of the T-MSD and the ensemble-averaged MSD (E-MSD). The experimentally determined ensemble-averaged MSD over a time interval m∆t is defined by: where J is the number of observed single-particle trajectories and t i is the starting time relative to first point in the trajectory.
Moreover, ergodicity breaking has been further confirmed by the presence of aging [26,27], i.e. the dependence of statistical quantities on the observation time. Based on these findings, several stochastic models presenting nonstationary (and thus nonergodic) subdiffusion have been proposed [20,[28][29][30][31]. Among these, CTRW has been used to model nonergodic subdiffusion in living cells [9,10,12] and has begun to provide theoretical insight into the physical origin of wEB in biological systems [28], associating the nonergodic behavior with the occurrence of particle immobilization with a heavy-tailed distribution of trapping times.
At the same time, these intriguing findings have generated new questions: Is nonergodic subdiffusion a strategy shared by other biological systems? Can biophysical mechanisms other than trapping lead to similar behaviors? What is its functional relevance? Elucidating these issues is crucial to unravel the role of nonergodic subdiffusion in cellular function. The main aim of the present work is to explore other forms of transport in biological systems to provide answers to these questions.
Here we used SPT to study the diffusion of a prototypical transmembrane protein, the Hamster Ovary (CHO) cells reproduces the essential features of the receptors naturally occurring on dendritic cells [34,35], thus serving as a valid model system. To characterize its dynamics, we performed video microscopy of quantum-dot labeled DC-SIGN stably transfected in CHO cells in epi-illumination configuration ( Fig. 1A-B, see Appendix A for details on cell culture and labeling procedures). In order to follow the standard biology nomenclature and to differentiate it from its mutated forms, in this manuscript we refer to the full receptor as the wild-type DC-SIGN (wtDC-SIGN).
We tracked quantum dot positions with nanometer accuracy by means of an automated algorithm [36]. We acquired more than 600 trajectories, all longer than 200 frames with some as long as 2000 frames, at a camera rate of 60 frame · s −1 to allow the evaluation and the comparison of time and ensemble averaged MSD. The T-MSD of individual trajectories displayed a linear behavior (β ∼ 1), consistent with pure Brownian diffusion (Fig. 1C). The fitting of the average T-MSD provided a value β = 0.95 ± 0.05. In addition, the distribution of the exponents β obtained by nonlinear fitting of the T-MSDs of the individual trajectories (inset of Fig. 1D) showed an average β = 0.98 ± 0.06.
Since the T-MSD values corresponding to different trajectories were broadly scattered, for each trajectory we calculated the diffusion coefficient D s by a linear fit of the T-MSD at time lags < 10% of the trajectory duration [37]. As expected, the resulting values of D s were found to have a very broad distribution, spanning more than two orders of magnitude ( Fig. 1D).
However, in marked contrast with the T-MSD, the E-MSD deviated significantly from linearity, showing subdiffusion with an exponent β = 0.84 ± 0.03 (Fig. 1E). The difference between the scalings of T-MSD and E-MSD is a clear signature of wEB [38]. To inquire whether DC-SIGN dynamics also exhibits aging, we computed the time-ensemble-averaged MSD (TE-MSD) by truncating the data at different observation times T : and extracting the corresponding diffusion coefficient D TE by linear fitting [37]. In systems with uncorrelated increments, it can be shown under rather general assumptions that D TE ∼ T β−1 [21,39]. The observed D TE indeed scaled as a power law with an exponent of −0.17 ± 0.05 (Fig. 1F), yielding a value of β in good agreement with the exponent determined from E-MSD. These results thus demonstrate that wtDC-SIGN dynamics exhibits aging.

III. FAILURE OF THE CTRW MODEL
The motion of some biological components, including the Kv2.1 potassium channel in the plasma membrane [12], lipid granules in yeast cells [9] and insulin-containing vesicles in Pancreatic β-cells [10], has been reported to exhibit subdiffusion compatible with the coexistence of an ergodic and a nonergodic process. The nonergodic part of the process has been modeled within the framework of the CTRW [28,38,39].
CTRW is a random walk in which a particle performs jumps whose lengths have a finite variance, but between jumps the walker remains trapped for random dwell times, distributed with a power-law probability density ∼ t −(1+β) , which for β ≤ 1 has an infinite mean. The duration of trapping events is independent of the previous history of the system. The energy landscape of this process is characterized by potential wells with a broad depth distribution. Such energetic disorder yields nonergodicity, since no matter how long one measures, deep traps cause dwell times on the order of the measurement time. Within the biological context, these traps generally have been associated with chemical binding to stationary cellular components (e.g. actin cytoskeleton [12] or microtubuli [10] In SPT experiments, the limited localization accuracy for determining the particle position sets a lower limit for the diffusivity value that can experimentally be measured. In our case, this lower threshold lies at D th = 6 · 10 −4 µm 2 s −1 . Therefore, a segmentation algorithm [41] was applied to the x-and y-displacements of our trajectories in order to detect events with diffusivity lower than D th . Surprisingly, transient trapping was only detected over less than 5% of the total recording time ( Fig. 2A-C). The detected trapping times displayed an average duration of 330 ± 30 ms (Fig. 2D). An alternative analysis, based on the transient confinement zone algorithm [42], gave comparable results [43].
In order to understand the nature of these trapping events, we attempted to fit their distribution by means of both an exponential and a power law distribution function ∼ t −(1+β) , as expected for CTRW [12]. The power law pdf provided a better fit to the data, yielding an exponent β = 0.83 ± 0.05 ( In addition, we constructed the distribution of escape times by identifying the duration of the events in which a trajectory remains within a given radius R TH (Fig. 2E). For a CTRW, the long-time dynamics is dominated by anomalous trapping events and, as a result, this quantity is expected to be independent of R TH [12]. In strong contrast to the CTRW model,

IV. DC-SIGN DISPLAY CHANGES OF DIFFUSIVITY
Recently, diffusion maps of the cell membrane have shown the presence of patches with strongly varying diffusivity [36,44,45]. Based on this evidence, we have recently proposed a class of models describing ordinary Brownian motion with a diffusivity that varies randomly, but is constant on time intervals or spatial patches with random size [21]. These models describe anomalous diffusion and wEB in complex and heterogeneous media, such as the cellular environment, without invoking transient trapping.
To address whether the observed nonergodic dynamics of DC-SIGN can be described with this theoretical framework, we further analyzed individual trajectories by means of a change-point algorithm to detect variations of diffusivity in time [41]. In brief, the algorithm consists in a likelihood-based approach to quantitatively recover time-dependent changes in diffusivity, based on the calculation of maximum likelihood estimators for the determination of diffusion coefficients and the application of a likelihood ratio test for the localization of the changes. Notably, DC-SIGN trajectories displayed a Brownian motion with relatively constant diffusivity over intervals of varying length, but that changed significantly between these intervals ( Fig. 3A-C). Similar features were identified in a large fraction of trajectories, with 63% showing at least one diffusivity change (Fig. 3D), in qualitative agreement with the models of random diffusivity [21].
To obtain a comprehensive understanding of our data, we considered an annealed model in which randomly diffusing particles undergo sudden changes of diffusion coefficient [21].
The distribution of diffusion coefficients D that a particle can experience is assumed to have a power-law behavior ∼ D σ−1 for small D (with σ > 0) and a fast decay for D → ∞. Given D, the transit time τ (i.e., the time τ a particle moves with a given D) is taken to have a probability distribution with mean ∼ D −γ (with −∞ < γ < ∞). Since the motion during the transit time τ is Brownian, particles explore areas with radius r ∼ √ τ D, and the radius of the region explored with such diffusion coefficient has probability distribution with mean We performed in silico experiments of 2D diffusion ( Fig. 4A-B), assuming a distribution of diffusion coefficients D given by: and a conditional distribution of transit times τ given by: where b and k are dimensional constants and Γ(x) is the Gamma function. The remarkable agreement between simulations and experimental data strongly supports heterogeneous diffusion as the origin for DC-SIGN nonergodicity.
It must be noticed that, in contrast to CTRW, our model does not assume particle immobilization, but a continuous distribution of diffusivity, with P D (D) ∼ D σ−1 for small D.
However, from the experimental point of view, it is not possible to distinguish immobilization from very slow diffusion. In fact, the limited localization accuracy of SPT experiments translates into a lower limit for the diffusivity value D th that can be detected. Therefore, in our analysis, trajectories, or portion of trajectories, with diffusivity lower than this threshold value (D th = 6 · 10 −4 µm 2 s −1 ) are identified as immobile, as shown in Fig. 2A-B. From the model described above, the distribution of the duration of these "apparent" immobilizations can be calculated as: We neglect here the possibility that the trajectory of an in-silico particle contains two consecutive segments characterized by diffusivities D i and D i+1 which are both smaller than D th , as this probability is vanishingly small for the parameter regime of our setup.
Independently of D th , the integral in Eq. (6) scales asymptotically as τ −1−β with β = σ γ , providing for the distribution of immobilization events the same behavior predicted by the CTRW [12]. Therefore, the distribution of immobilization times in Fig. 2D is fully compatible with the prediction of our model, further confirming its agreement with the experimental data.

V. DYNAMICS OF RECEPTOR MUTANTS
From the structural point of view (Fig. 5A), DC-SIGN is a tetrameric transmembrane protein, with each of the four subunits comprising: (i) an extracellular domain that allows binding of the receptor to pathogens, i.e., ligand binding domain; (ii) a long neck region; and (iii) and a transmembrane part followed by a cytoplasmic tail that allows interactions with inner cell components and facilitates the uptake and internalization of pathogens. [46,47]. Moreover, DC-SIGN contains a single N -glycosylation site mediating interactions with glycan-binding proteins [43]. To gain insight into the molecular mechanisms of DC-SIGN nonergodic diffusion, we generated three mutated forms of the receptor (Fig. 5A). These mutations have been reported to modify the interaction of DC-SIGN with other cellular components, strongly affecting DC-SIGN function [43,48,49]. The N80A mutant lacks the N -glycosylation site. This defect hinders interactions of DC-SIGN with components of the extracellular membrane that bind to sugars [43]. The ∆35 mutant lacks a significant part of the cytoplasmic tail, preventing interactions with cytosolic components such as actin [48].
Finally, the ∆Rep mutant lacks part of the neck region, abrogating interactions between different DC-SIGN molecules [49].
We found that each mutation has a very different effect on the dynamics of the receptor.
The N80A mutant (Fig. 5C-F Overall, these results demonstrate that the molecular structure of the receptor strongly influences its diffusive behavior on the cell membrane and the occurrence of weak ergodicity breaking.

VI. NONERGODICITY AND BIOLOGICAL FUNCTION
Together with our previous biophysical studies on DC-SIGN [43,49], the data and analysis presented in this paper allow us to link the dynamical behavior of DC-SIGN to its functional role in pathogen capture and uptake (known as endocytosis).
In terms of steady-state organization, wtDC-SIGN, N80A and ∆35 preferentially form nanoclusters on the cell membrane, which are crucial for regulating pathogen binding [43,49], whereas removal of the neck region (∆Rep) reduces nanoclustering and binding efficiency to small pathogens, such as viruses [49]. Our results thus show that the diffusive behavior of the receptor is strongly linked to nanoclustering, but not merely due to size-dependent diffusivity and/or time-dependent cluster formation and breakdown. In fact, dual-color SPT experiments performed at high labeling density do not reveal correlated motion between nearby DC-SIGN nanoclusters, excluding the occurrence of dynamic nanocluster coalescence [43].
Moreover, although superresolution imaging has revealed that wtDC-SIGN, N80A and ∆35 form nanoclusters with similar distributions of size and stoichiometry [43,49], our dynamical data evidence significant differences in their diffusion patterns (Figs. 1 and 5).
Our data are in fact consistent with the view of the plasma membrane as a highly dynamic and heterogeneous medium, where wEB stems from the enhanced ability of DC-SIGN nanoclusters to interact with the membrane environment, including components from the outer and inner membrane leaflet. This interaction is inhibited (or strongly reduced) in the case of the ∆Rep mutant since it does not form nanoclusters [49]. As a result, the motion of ∆Rep is Brownian and ergodic, and interestingly this dynamic behavior correlates with its impaired pathogen binding capability [34,49].
In contrast, we observed that both wtDC-SIGN and N80A, which show a similar degree of nanoclustering [43], exhibit wEB. But, the distribution of diffusivity of N80A is significantly broader than that of wtDC-SIGN, and is shifted towards lower diffusivity values ( Fig.   5C-F). This increased heterogeneity correlates with altered interactions of the N80A with extracellular components, resulting from the removal of the glycosylation site. Indeed, we have recently shown that the N80A mutant has a reduced capability to interact with extracellular sugar binding partners [43]. Thus, it appears that the extracellular milieu next to the membrane contributes to the degree of dynamical heterogeneity sensed by the receptor.
Remarkably, this correlation also extends to the functional level, as we have recently shown that interactions of DC-SIGN with extracellular sugar-binding proteins influence encounters of DC-SIGN with the main endocytic protein clathrin. In turn, this resulted in reduced clathrin-dependent endocytosis of the receptor and its pathogenic ligands [43].
Finally, the ∆35 mutant exhibits nanoclustering [49] and wEB similar to that of wtDC-SIGN, From the biological point of view, however, this mutant is not able to interact with cytosolic components in close proximity to the inner membrane leaflet, including actin [48].
Therefore, in contrast to the extracellular influence observed for the N80A mutant, the results obtained for the ∆35 mutant indicate that interactions with the actin cytoskeleton, responsible for the CTRW-like behavior of other proteins [12], do not play a major role in DC-SIGN wEB. Nevertheless, it should be mentioned that the reduced endocytic capability of the ∆35 could not be uniquely attributed to its dynamic behavior on the cell membrane but rather to its impaired interaction with downstream partners involved in internalization [48,49].

VII. CONCLUSIONS AND OUTLOOK
We have demonstrated that the DC-SIGN receptor displays subdiffusive dynamics, char- might provide a deeper comprehension of these mechanisms at the molecular level.
The model used to interpret our data provides a flexible and realistic framework to describe anomalous motion in cell membranes. Although in the present work we have focused our simulations on time-dependent changes of diffusivity, similar conclusions can be obtained assuming a spatial dependence, with constant diffusivity on membrane patches of random size [21]. The current data do not allow discrimination between the two scenarios. The application of techniques that combine dynamic and spatial mapping at high labeling conditions [45,50] would be necessary to verify the occurrence of spatial maps of diffusivity.
In addition, numerical simulations of spatial-dependent random diffusivity require the construction of 2-dimensional diffusivity maps consistent with the model's probability densities, which is a non-trivial task.
While the work presented here focuses on the cell plasma membrane, we point out that these results have much broader implications. In fact, our model and analysis are very general and can be applied to any diffusive system that shows wEB, in order to investigate the role that heterogeneous diffusivity plays in observed anomalies. Fundamental questions about the nature of anomalous and nonergodic diffusion in disordered media arise in many fields, such as life sciences [28], soft condensed matter [51,52], ultracold gases [53,54], geology [55] and ecology [56]. Our work provides an alternative conceptual framework and specific tools for answering these questions. controller (TC-324B, Warner Instrument) and a digital CO 2 controller (DGTCO2BX, Okolab) at 37 • C and in 5% CO 2 atmosphere. Trajectories were analyzed with custom Matlab code based on the algorithm described in Ref. [36]. In order to avoid artifacts in trajectory reconnection caused by quantum dots blinking dynamics and/or high local density of quantum dots, each trajectory was terminated at the first video frame not displaying a clearly identifiable bright spot in the surrounding of the quantum dot localization obtained in the previous frame. Similarly, the trajectory reconstruction was also interrupted if the presence of multiple bright spots did not allow unambiguous identification of the same quantum dots in successive video frames. As a further check for false-positive reconnection, trajectories were overlaid to raw movies and visually inspected.

Appendix C: Data analysis
Time-, ensemble-and time-ensemble-averaged mean-squared displacements were calculated as described in [12]. Exponents of the E-MSD, average T-MSD and D t,ens were obtained by linear fitting of the log-log transformed data. Errors were calculated as the 99% confidence interval of the fitting parameters. Short-time diffusion coefficients were extracted from the linear fit of the first 10% of the points of T-MSD curves [37].

Measurements of the apparent diffusion of quantum dots on fixed cells and glass coverslips
were used to estimate the smallest detectable diffusivity. Short-time diffusion coefficients were obtained as described above for trajectories of immobilized quantum dots and the corresponding probability distribution was calculated. 95% of the immobilized quantum dots trajectories showed values lower than 6 · 10 −4 µm 2 /s, which was therefore set as the threshold (D th ) for classifying a trajectory as mobile.
Dynamical changes in the motion of DC-SIGN receptors were identified by application of the change-point algorithm described in Ref. [41]. In brief, the trajectories were recursively segmented and a maximum-likelihood-ratio test was applied to the trajectory displacements (∆x, ∆y) in order to identify sudden changes of diffusivity. The critical values for Type I error rates were set to a confidence level of 99%, corresponding to 1% probability of having a false-positive identification of a change-point. For each dynamical region identified by the algorithm, the short-time diffusion coefficient was calculated from a linear fit of the first 10% of the points of the corresponding MSD curves [37]. Regions showing a short-time diffusion coefficient lower than D th were considered compatible with transient immobilization.

Appendix D: Simulations
Simulated trajectories (500 per parameter set) were obtained by generating random diffusion coefficients D according to the probability distribution given in Eq. (4). For each diffusion coefficient, the corresponding transit time τ was calculated as a random number drawn from the distribution given in Eq. (5). Particle coordinates r = {x, y} were generated as: