Protein crowding in lipid bilayers gives rise to non-Gaussian anomalous lateral diffusion of phospholipids and proteins

Biomembranes are exceptionally crowded with proteins with typical protein-to-lipid ratios being around 1 ∶ 50 − 1 ∶ 100 . Protein crowding has a decisive role in lateral membrane dynamics as shown by recent experimental and computational studies that have reported anomalous lateral diffusion of phospholipids and membrane proteins in crowded lipid membranes. Based on extensive simulations and stochastic modeling of the simulated trajectories, we here investigate in detail how increasing crowding by membrane proteins reshapes the stochastic characteristics of the anomalous lateral diffusion in lipid membranes. We observe that correlated Gaussian processes of the fractional Langevin equation type, identified as the stochastic mechanism behind lipid motion in noncrowded bilayer, no longer adequately describe the lipid and protein motion in crowded but otherwise identical membranes. It turns out that protein crowding gives rise to a multifractal, non-Gaussian, and spatiotemporally heterogeneous anomalous lateral diffusion on time scales from nanoseconds to, at least, tens of microseconds. Our investigation strongly suggests that the macromolecular complexity and spatiotemporal membrane heterogeneity in cellular membranes play critical roles in determining the stochastic nature of the lateral diffusion and, consequently, the associated dynamic phenomena within membranes. Clarifying the exact stochastic mechanism for various kinds of biological membranes is an important step towards a quantitative understanding of numerous intramembrane dynamic phenomena.


I. INTRODUCTION
The dynamics of biomembranes plays a crucial role in the regulation of numerous cellular processes. This is largely due to the fact that membrane proteins carry out a substantial fraction of cellular functions. For instance, they are involved in cellular signaling, in which the functional complex can be a single protein as well as an oligomer [1][2][3] together with an appropriate pool of lipids modulating the protein function by protein-lipid interactions [4][5][6][7]. Since the formation of protein-lipid complexes is reversible, proteins and lipids are repeatedly probing the membrane for functionally favorable surroundings, thereby allowing the diffusive motion to largely control their reaction kinetics.
Cell membranes are known to be extremely complex fluids characterized by heterogeneity [6,8] and compartmentalization [9,10]. Similar to the crowded cytoplasm of biological cells [11], membranes are crowded with proteins and other macromolecules [12]. These effects substantially complicate the relationship between the dynamics and function of cell membranes. Notably, crowding has been a neglected feature of the intracellular environment [11] until recently, when numerous studies have identified the role of crowding in multiple phenomena including, inter alia, protein stability [13], signaling [14,15], and gene transcription [16]. Most notably, crowding shapes the reaction kinetics [11,[17][18][19] by modifying the mobility [20] as well as association rates [21][22][23][24][25]. This multitude of examples together with the shift in dimensionality from three to two hints that, besides the cytoplasm, crowding plays a substantial role also in cellular membranes. Indeed, similarly to the cytoplasm [26][27][28], one can argue that reaction kinetics in two dimensions are optimized by the Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. degree of crowding as the slower mobility is compensated by the higher reactant concentration [29,30].
Cytoplasmic crowding in living biological cells is known to cause anomalous diffusion with a mean-squared displacement (MSD) of the form [31,32] with the generalized diffusivity K α of physical dimension m 2 =s α and the anomalous diffusion exponent α [33][34][35]. Normal Brownian motion corresponds to the limit α ¼ 1.
In cellular membranes crowding is known to induce heterogeneous and confined motion. While some studies argue that crowding leads to the slowing-down of diffusion in terms of an effective diffusivity in an otherwise Brownian diffusion scenario [48], others suggest that diffusion becomes anomalous instead [48][49][50][51][52][53][54][55][56][57][58]. Anomalous diffusion of the form Eq. (1) was observed in various cellular membranes [32,55,57,58], emphasizing its general nature in membranes. Intriguingly, while anomalous diffusion would be expected to hinder the efficiency of diffusionlimited processes, in suitably crowded conditions it can surprisingly also result in favorable effects, such as to enhance the search of nearby reaction partners, which further leads to increased protein complex formation and a subsequent boost in reaction rates [59,60]. Confinement and corralling effects can also lead to enhanced protein oligomerization [61] and bursts in reaction rates [62]. Another important effect is the observation of aging [63], the explicit dependence of the dynamics of different membrane proteins on the length of the observation time [55,58].
While a lot of work has been conducted on the development of reaction-diffusion theories in the subdiffusive regime [64,65], such theories have not been widely applied to biological systems [15,32,66,67]. Therefore, resolving the diffusion mechanisms of molecules in crowded and compartmentalized conditions of cellular membranes could greatly improve our understanding of diffusion-controlled reactions in the cells. Moreover, membrane protein complexes [2,68] are targeted by a major fraction of current drugs. Therefore, understanding and thereafter altering the dynamics that govern the formation of functional proteinprotein [69,70], lipid-protein [71,72], or domain-protein [73] units is an intriguing approach to the treatment of numerous diseases.
In this work, we show through extensive molecular simulations and theoretical analysis that protein crowding changes the membrane dynamics drastically. We find that the dynamics of lipids in crowded conditions is no longer described by the mechanism consistent with the fractional Langevin equation (FLE) typically associated with the lipid motion in noncrowded membranes [52,53], or by any single known mechanism. Instead, the motion becomes non-Gaussian and heterogeneous, yet maintains its ergodic properties. In particular, while the time-averaged MSD scales sublinearly, no aging is observed. Concurrently, a strong dynamic heterogeneity is observed among different lipids as well as membrane-embedded proteins. Our findings are central to resolving the mechanisms governing how molecules move along crowded membranes and, therefore, to understanding how the multitude of processes occurring in cellular membranes are controlled by anomalous diffusion-reaction dynamics.
This paper is structured as follows. After a brief introduction to the model and simulations in Sec. II, we present our results in Sec. III. In particular, we demonstrate that the motion of the membrane constituent molecules is multifractal and anomalous. In contrast to noncrowded membranes, significant non-Gaussian shapes for the particle diffusion are observed and the inhomogeneity of the motion of individual particles is shown to be slowly varying with time. In Sec. IV, we discuss our results and conclude the paper. Additional results to support those presented in the main text are given in Supplemental Material [80].

II. MODELS AND SIMULATION
For our investigation of the stochastic characteristics of the lateral molecular diffusion in protein-crowded lipid membranes, we simulate three lipid bilayer systems similar to those employed in our study of the emergence of anomalous diffusion in the presence of membrane-embedded proteins [49]. Two protein-crowded membrane systems are modeled by embedding a total of 16 NaK channels (Protein data bank (PDB) entry 2AHY), respectively, in bilayers of 1600 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and 1,2-dilauroyl-sn-glycero-3phosphocholine (DLPC) lipids. The lipid-to-protein ratio is 1∶50 per leaflet, corresponding to a protein surface coverage of approximately 34%, which is known as a typical value in cellular membranes [81]. Both systems span a lateral membrane area of approximately 25×25nm 2 with periodic boundary conditions. The main difference between the two systems, as shown in Fig. 1, is that the NaK proteins tend to strongly aggregate with the others to form a larger complex in the DPPC system while they do so to a much lesser extent in the DLPC membrane. This is most likely due to a different hydrophobic mismatch of the proteins with the two lipid environments. In this spirit, the two systems chosen for this study are used to gauge the dynamics in two physically different settings often found in cell membranes-protein aggregating and protein nonaggregating. As a reference, we also simulate a proteinpoor membrane composed of 2045 DPPC lipids and a single NaK channel.
In order to study large systems consisting of multiple proteins over a time scale of 100 μs, the computationally efficient coarse-grained Martini-based [82,83] force field [84] is employed. The protein-poor case is simulated for 25 μs, which corresponds to 100 μs of efficient time when the commonly used coarse-grained-to-atomistic comparison factor of 4 is taken into account. Compared to the protein-poor case, the protein-crowded systems turn out to have slower lateral dynamics [49]. Because of this, in this work, we simulate them for 100 μs (corresponding to 400 μs of effective time), which is fourfold longer than our previous simulation (25 μs) for similar systems. All simulation parameters are identical to our previous simulation study [49]. The analyses on the simulation results below are carried out in terms of the real simulation time.

III. RESULTS
In this section, we showcase the various aspects of the protein and lipid motion in the simulated lipid membranes. As we see, major changes are effected at higher protein crowding fractions.

A. Protein crowding induces multifractal anomalous lateral diffusion in membranes
We start our analysis by examining the diffusion properties of single phospholipid and protein molecules in our model membranes. From an individual two-dimensional trajectory r i ðtÞ encoding the motion of such a single molecule, the time-averaged MSD, is evaluated as a function of the lag time Δ and the overall observation time T. Here and throughout the study, r i ðtÞ refers to the in-plane coordinate of the center of mass of a given ith lipid or protein molecule with respect to the center of mass of the entire membrane, thus removing the effect of membrane drift during the simulations [53].
In protein-rich membranes the membrane drift does not need to be removed separately with respect to each leaflet because the drift of the monolayers with respect to each other is negligible due to the anchoring effect of the proteins; see Figs. S1 and S2 in Supplemental Material [80]. The advantage of the time-averaged MSD Eq. (2) over the ensemble average Eq. (1) is that for sufficiently long trajectories r i ðtÞ resolves differences between the motion of individual particles, as we show below. We note that in so-called weakly nonergodic systems the long-time scaling of the time-averaged MSD δ 2 i ðΔÞ may differ fundamentally from the corresponding ensemble average hr 2 ðtÞi, reflecting the nonstationarity of the underlying motion [31,32,35,85,86]. Experimentally and from simulations, such weak nonergodicity was indeed observed in the cytoplasm [40,41] and the plasma membrane [55,58] of living cells, as well as quite different systems such as blinking quantum dots [87] or granular gases [88]. (ii) protein-rich membrane system composed of 1600 DPPC lipids and 16 NaK proteins: aggregating system; (iii) protein-rich membrane system composed of 1600 DLPC lipids and 16 NaK proteins: nonaggregating system; (iv) schematic structures of a DPPC phospholipid and a NaK channel protein employed in our coarsegrained simulations. For both DPPC and the NaK channel, the transparent coarse-grained structure is shown on top of the atomistic representation.
noncrowded and crowded membranes. Along with the individual time-averaged MSDs, the mean as the blue solid line is plotted for all lipids. From this averaged curve, the anomalous scaling exponents α are estimated as a function of the lag time Δ, shown in the upper panels in Fig. 2. The results show that protein crowding significantly affects the character of the lateral diffusion in the membranes. While the lipids and the single protein in the noncrowded membrane exhibit trivial Brownian diffusion at lag times Δ longer than some 10 ns, consistent with previous studies [52][53][54]89,90], those in the crowded membrane exhibit significant anomalous diffusion. The plot of α as a function of Δ tells us that the lipid carries out multifractal subdiffusion until tens of microseconds, where the scaling exponent α is temporally varied with Δ in a broad spectrum of values. Additionally, the subdiffusive dynamics extends for around 3 orders of magnitude longer than for noncrowded lipid membranes. The variation in α indicates that the lipid motion is strongly affected by the proteins until around 10 μs. The behavior at time scales of > 10 μs is not conclusive for our analysis due to the limited simulation time. It appears that the strong fluctuation in the scaling exponent αðΔÞ for the proteins stems from the comparatively small numbers of proteins and will be different for independent runs. Interestingly, the diffusion dynamics of the proteins is affected by their own crowded state. In contrast to the normal diffusion of a single protein in a noncrowded membrane, 16 proteins in the crowded membranes exhibit strong subdiffusion. This effect is also clearly seen in the corresponding trajectories (see Fig. S3 in Supplemental Material [80]). Without any signature of transient behavior, it persists until the end of the simulations at 100 μs with an anomalous diffusion exponent that is smaller than the value for the lipids. The magnitude of α becomes smaller when the proteins are aggregated in the DPPC membrane compared to the case where they do not aggregate in the DLPC membrane. These results are consistent with our earlier study [49], even though the exact lag time dependence of α is affected by the confinement effects which strongly depend on the incidental aggregation geometry of the proteins. Note that the lateral dynamics of the proteins is about 10 times slower than the lipid dynamics. It is also noteworthy that, owing to the protein crowding, individual lipids have various diffusion patterns which are distinct from the averaged behavior hδ 2 ðΔÞi; in both protein-rich DPPC and DLPC membranes there are some lipids whose time-averaged MSD curves almost completely follow those of the proteins, which suggests that these lipids move together with the proteins during the entire simulation time, possibly in nonannular binding sites.
In line with this view, it is found that approximately 30 lipids in the vicinity of a protein are considerably slowed down compared to lipids far from the protein (see Fig. S13 in Supplemental Material [80]). For the case of the In this plot, the offset Δ ¼ 10 ns is due to the finite interval in fitting the instantaneous slope of the time-averaged MSD. Note the different scales of the ordinates in the different panels. We check that the change of the fit range does not seriously alter the αðΔÞ curves. Note that, in general, the scaling exponents α of the proteins exhibit larger fluctuation than those of the lipids. In particular, its sudden decrease for Δ > 5 μs shown in the noncrowded membrane reflects the statistically insignificant fluctuation of a single time-averaged MSD curve.
protein-aggregating DPPC membrane, some other lipids, while unbound from the proteins, undergo restricted subdiffusion in confined space in between aggregated proteins; see the sample trajectory in Fig. S3 in Supplemental Material [80]. Such a confined diffusive pattern is not observed if the lipids easily move through the nonaggregating proteins in the DLPC membrane. The majority of lipids appear to carry out unrestricted anomalous diffusion in the bulk. We quantify the effect of protein crowding on the mobility of individual lipids by comparison of the amplitude scatter in the time-averaged MSD for the protein-poor and protein-rich membranes in Fig. S4 in Supplemental Material [80]. The amplitude scatter has a markedly broader distribution in the protein-crowded membrane. Notably, in this case, the shape of the distribution becomes asymmetric with respect to the average because of the existence of a fraction of less mobile lipids caused by the presence of the proteins.

B. Protein crowding induces non-Gaussian anomalous diffusion
We now proceed to investigate whether the fractional Langevin equation, identified previously as a governing process for the lateral diffusion in pure lipid bilayers [52][53][54], is still able to reproduce the behavior of the lateral motion in protein-crowded membranes. The FLE, a special case of the generalized Langevin equation [91] with power-law friction kernel compensating the power-law correlations of the driving fractional Gaussian noise, typically describes the motion of tracer particles in viscoelastic media [92].
We first examine the (non-)Gaussianity of the diffusion process by looking at the cumulative distribution of the squared displacements Πðr 2 ; ΔÞ ¼ R r 0 Pðr 0 ; t ¼ ΔÞ2πr 0 dr 0 for the two-dimensional motion [55,93,94]. Here, P is the propagator, i.e., the probability density that the spatial increment of the particle is found to be at [r, r þ dr] over the lag time interval Δ. Gaussian anomalous diffusive processes such as the FLE are described by the propagator in free space Pðr; Thus, the plot of − log½1 − Πðr 2 ; ΔÞ versus r 2 displays a power-law scaling with the exponent 2 [94]. As we show in Fig. 3, the lipids in the noncrowded membrane are in very good accord with the Gaussian scaling law.
However, both protein-crowded systems in Fig. 3 clearly show that the lateral diffusion is not Gaussian. It appears that the cumulative distributions for the crowded membranes follow a power-law relation − log½1 − Πðr 2 ; ΔÞ ∼ r δ , with a single or multiple scaling exponent δ < 2 over the entire range of r 2 plotted. It turns out that these non-Gaussian scaling curves are found to be consistently preserved at different overall observation times T (see Fig. S5 in Supplemental Material [80]). Such a power-law scaling with an exponent δ ≠ 2 may be explained with a non-Gaussian propagator P of the form which leads to a cumulative distribution − log½1 − Πðr 2 ; ΔÞ ∼ ðr=cΔ α=2 Þ δ in the large-displacement limit.
Here, P r ðr; ΔÞ is the radial part of the propagator given by P r ðr; ΔÞ ¼ R 2π 0 dθPðr; θ; ΔÞ, with r ≥ 0, and c is a scaling constant of dimension [m=s α=2 ]. We confirm the non-Gaussianity of the propagators in Fig. 4 by examining the profiles of log P r ðr; ΔÞ at several lag times with the fit curves based on the above non-Gaussian (solid lines) and the Gaussian (dashed lines) propagators. The lipids in the protein-poor membrane (Fig. 4, top left) have the anticipated Gaussian propagator at all given times. In the protein-crowded DPPC membrane (Fig. 4, top right and bottom), the two fit curves unambiguously suggest that the propagators are non-Gaussian with exponents δ < 2. Moreover, we find that the propagators require a composite fit function involving two non-Gaussian components (see the caption of Fig. 4 and also Fig. S7 in Supplemental Material [80]) to fully explain the behavior at both small and large displacements. For the two-component non-Gaussian propagator (shown in the caption of Fig. 4), its cumulative distribution can be exactly obtained as where γðs; xÞ ¼ R x 0 t s−1 expð−tÞdt is an incomplete gamma function. The solid lines in Fig. 3 depict the theoretical cumulative distribution Eq. (5) expected from the non-Gaussian propagator, which fits the simulation data in Fig. 4. For all cases, the theoretical curves display excellent agreement with the cumulative distribution, which in turn supports the validity of the proposed non-Gaussian propagator. We notice that the theoretical curves do not show a tip at the end. This suggests that the sharp bending at the tip in the simulation likely results from insufficient statistics. As we discuss further below, the character of the above non-Gaussian propagator is different from the non-Gaussian propagators in other complex soft matter systems found in recent experimental studies [95][96][97][98][99]. There is no single master curve describing all propagators obtained at different times, since the estimated δ values vary with time. The variation trend of δ in time appears to be system specific depending on the extent of the transient confinement effect due to protein crowding; see Fig. S6 in Supplemental Material [80]. For the protein-crowded DPPC system (where many lipids experience strong transient confinement), the values of effective δ tend to decrease with increasing lag time Δ. This indicates that the non-Gaussianity is enhanced until, at least, Δ ¼ 5 μs as time elapses. In contrast, the tendency is opposite in the nonaggregating DLPC system, in which δ increases and thus P becomes more Gaussian over time.
As our analysis shows, the non-Gaussian lateral anomalous diffusion revealed by Π and P cannot be explained by the model of multiple-component Gaussian processes that has often been introduced for modeling complicated anomalous diffusion processes found in experiment [93,100]. For our protein-crowded membranes, one can invoke a twocomponent (anomalous) diffusion model composed of slowly diffusing lipids in the proximity of the proteins and fast lipids in bulk [101], with the propagator Pðr;ΔÞ¼ mexp½−r 2 =ð2σ slow Þ=ð2πσ slow Þþð1−mÞexp½−r 2 =ð2σ fast Þ= ð2πσ fast Þ, with 0 < m < 1. As demonstrated in Fig. S7 in Supplemental Material [80], this model does not properly fit the distribution of displacements at various lag times. Additionally, we examine in Fig. S8 [80] whether the propagators are explained by the two-component model composed of a Gaussian center and a non-Gaussian tail as employed in some recent studies [95][96][97]. While this model explains the propagators better than the two-component Gaussian model, it also shows a deviation from the propagators, especially when the center is sharply peaked.
The protein diffusion exhibits similar characteristics. The corresponding cumulative distributions and propagators are plotted in Fig. 5. This analysis demonstrates that an isolated single protein in the noncrowded membrane (top panels) shares the property of Gaussian diffusion. The plot of − log½1 − Πðr 2 ; ΔÞ scales as r 2 , and P r ðr; ΔÞ follows the Gaussian propagator (dashed lines). The inconsistent cumulative distribution curve for Δ ¼ 5 μs is statistically inconclusive because it was obtained from a single trajectory of length T ¼ 25 μs. The proteins in the crowded DPPC (middle panels) and DLPC (bottom panels) membranes display non-Gaussian scaling with δ being significantly below 2 at any analyzed displacements. The depicted solid lines are the expected cumulative distributions [Eq. (5)] from the fitted non-Gaussian propagators. The non-Gaussian protein diffusion appears to be more complex than the non-Gaussian behavior of the lipids, presumably due to the slow lateral dynamics of the proteins and nontrivial, effective lipid-mediated interactions among proteins. The cumulative distributions are typically characterized by multiple scaling curves.
We perform additional analyses to cross-check the failure of the FLE model for the protein-crowded membranes. Figure 6 presents the results of the moment ratios hr 4 i=hr 2 i 2 (regular) and hr 4 MME i=hr 2 MME i 2 (mean maximal excursion [102]) for lipids in the protein-poor and proteinrich membranes [53,102]. Note that in the literature the  Fig. S7 in Supplemental Material [80]. regular moment ratio is alternatively redefined to a non-Gaussian parameter α 2 ðtÞ¼ð1=2Þhr 4 ðtÞi=hr 2 ðtÞi 2 −1 [32,103]. It was theoretically shown that for a twodimensional Gaussian process including FLE motion the moment ratio should be 2 for the regular moment and < 1.49 for the MME [53,102]. The lipids in the noncrowded membrane are in good agreement with these criteria. In contrast, in Fig. 6, the lipids in the crowded membranes are shown to have regular moments > 2 and the MME > 1.49 (dashed line), thus disobeying the Gaussianity criteria. Again, the obtained moment ratios are inconsistent with the FLE Gaussian model. We also note that the obtained moment ratios do not satisfy the criteria for anomalous diffusion in fractal media, which is characterized by moment ratios < 2 for the regular and < 1.49 for the mean maximal excursion [102]. This rules out the possibility that the protein-induced non-Gaussian diffusion in our model membranes is caused by a fractallike obstacle structure formed by proteins.
The velocity autocorrelation function is also found to be inconsistent with the FLE model. From the trajectory, we obtain CðΔÞ for average velocities v δt ðtÞ ¼ ½rðt þ δtÞ − rðtÞ=δt, with a varying time interval δt [53,104], comparing them to the corresponding theoretical curve of the FLE model. Previously it was shown  Here, the MME distance r MME ðtÞ refers to the maximum distance reached from the origin until time t [102]. For two-dimensional FLE motion with 0 < α < 1, the moment ratio is 2 for regular and < 1.49 (dashed line) for MME. (Bottom right) Velocity autocorrelation function CðΔÞ ¼ hv δt ðt þ ΔÞ · v δt ðtÞi for the average velocity v δt ðtÞ ¼ ½rðt þ δtÞ − rðtÞ=δt with the chosen time interval δt ¼ 100 ðnsÞ. The result for DPPCs in the protein-rich membrane is compared to the corresponding theoretical curve of FLE [53]. that the CðΔÞ is excellently explained by the FLE model in the noncrowded case [53,54]. Figure 6 (bottom right) shows the result for DPPC lipids in the protein-rich membrane. While the antipersistent tendency in the lipid motion is observed, the relaxation profile does not follow the FLE model.

C. Protein crowding induces spatiotemporal heterogeneity in lateral diffusion
In order to obtain more insight into the physical nature of the protein-induced non-Gaussian diffusion, we visualize the spatiotemporal character of the lateral diffusion through diffusion maps at different times in Fig. 7. Since in our membrane systems the diffusion is anomalous with a timedependent scaling exponent, neither the map of the generalized diffusion coefficient K α ðx; y; tÞ nor the scaling exponent αðx; y; tÞ properly quantify the degree of the local diffusion preference. We thus generate a 2D contour map of the local mean-squared displacements ½rðt þ t 0 Þ − rðtÞ 2 at time t over a short time interval t 0 set to be 1 ns in our analysis. Usually there are only 1-2 particles (or even none) at each site of the square grid in space, so the averaged local mean-squared displacements can be highly fluctuating in space regardless of the true local diffusional preference. To obtain a statistically meaningful contour map averaged over these fluctuations, we evaluate the diffusion map at time t by averaging 10 consecutive contour maps evaluated at times t, t þ 1; …; t þ 9 ns, based on the observation that during 10 ns the protein position remains almost unchanged compared to that of the lipids. By this procedure, the diffusion maps of lipids are evaluated at t ¼ 1, 10, 100, and 1000 ns. Figure 7 compares the diffusion maps of DPPC molecules in protein-poor and protein-rich cases. Therein, the blue regions represent the unoccupied space by the lipids, which almost overlaps with the dashed circles denoting the positions of the NaK proteins at that moment, as estimated from the trajectory. In the case of the protein-poor membrane, the profile of the diffusion map changes with time. Except for the tendency that the lipid diffusion becomes slowed down in the proximity of a protein, the hot regions in which large lipid diffusion occurs emerge randomly in space. The rapidly varying spatiotemporal heterogeneity effectively produces homogeneous lipid diffusion over time, as seen above in the time-averaged MSD curves and their distribution. In sharp contrast to this, the diffusion map for the protein-crowded membrane exhibits a heterogeneous profile that strongly depends on the protein configuration in space. The diffusion tends to be slower in the protein-crowded regions and faster in between two distant protein complexes. It is surprising that the proteininduced diffusion heterogeneity has a very long life span. The given pattern of the diffusion map is maintained for more than 1 μs. Thus, a lipid molecule will experience a space-dependent diffusivity while diffusing across the membrane, as quantified below. The persistent pattern of the diffusion map originates from the slow diffusion dynamics of the protein complexes compared to the lipid motion in the investigated time range.
The spatiotemporal heterogeneity in the lateral diffusion of the lipids is further corroborated in Fig. 8 the increase of T, we evaluate the ergodicity breaking parameter [35,104,105] EBðTÞ ¼ In the noncrowded membrane the heterogeneity of the time-averaged MSD in Fig. 8 is consistent with that for typical homogeneous ergodic diffusion, such as FLE motion. While individual time-averaged MSD traces show an erratic profile, their mean depicted by the thick gray line is independent of the observation time T, hδ 2 ðTÞi ≃ T 0 , and the time-averaged MSD heterogeneity decreases with increasing T in accordance with the expected convergence law EBðTÞ ∼ T −1 as for FLE [106][107][108][109]. Meanwhile, the lipids in the protein-rich membranes suffer strong heterogeneous diffusion. There exist lipids showing continuously decreasing, aginglike time-averaged MSDs upon increasing T. Additionally, it is observed that a fraction of lipids experience a sudden increase in the time-averaged MSD at a certain T 0 while having been independent of T for shorter observation times T < T 0 . From the diffusion map it can be inferred that those lipids diffusing successively through slow and fast diffusivity regions in space can have such temporal mobility heterogeneities. This new picture on the lipid diffusion in protein-crowded membranes is corroborated in Fig. 9, where we plot the temporal fluctuation of the diffusivity K α for some individual lipid molecules and the probability density PðK α Þ for the obtained diffusivity from all trajectories. The traces of K α ðtÞ indeed show that the individual lipids display fluctuating diffusivity. The bimodal profile of P tells us that the diffusivity fluctuates within a finite range of values with two favorable diffusivity states. For some lipids (e.g., Fig. 9, bottom left), the transitions between the high and low diffusivity states are clearly seen. In this case, the lipids explore the low diffusivity region with a surprisingly long sojourn time of ∼10 μs. Such patterns of diffusivity fluctuation lead to the decrease or the increase of time-averaged MSD δ 2 ðTÞ with increasing T, as observed in Fig. 8. A supporting simulation described in Sec. S7 in Supplemental Material [80] suggests that the temporal fluctuation in δ 2 ðTÞ may be attributed to a temporal change of the diffusivity whether or not the particle is under transient confinement. Geometrical effects like confinement are unable to change the amplitude of δ 2 ðTÞ, which can be decreased or increased if the diffusivity is changed with time. In spite of the strong heterogeneity on the single molecule level, on average the lipid diffusion does not age in the sense of hδ 2 ðTÞi ≃ T 0 . This ergodic nature of the lipid diffusion is consistent with the diffusivity fluctuation in Fig. 9 in that a typical sojourn time for a given diffusivity state is always bounded. In other words, no diverging sojourn time exists for a low diffusivity state, which may be necessary for inducing nonergodicity, as shown in recent studies on the fluctuating diffusivity model [58,110,111]. The EB parameter indicates that the space-dependent diffusivity only leads to a slower ergodic convergence than the conventional Brownian case EBðTÞ ∼ T −1 . This observation indicates the importance of the proper averaging procedure in the analysis of the crowded membrane dynamics.

IV. DISCUSSION AND CONCLUSIONS
Based on extensive computer simulations and single trajectory analyses, we here provide a systematic  PðK α Þ for the obtained K α value from all DPPC lipid molecules in the protein-crowded membrane. The dashed line represents a double-Gaussian fit to the data. Plots corresponding to the protein-poor system and nonaggregating systems can be found in Figs. S11 and S10 in Supplemental Material [80].
investigation of the anomalous lateral molecular diffusion of lipid molecules and embedded proteins in proteincrowded membranes. Recently, several experimental and computational studies reported that anomalous lipid diffusion is present in protein-crowded membranes up to time scales of microseconds to milliseconds and even up to 100 s in living cells [49,55,58,100]. Going beyond these studies, our focus here is to pin down the stochastic character and physical origins of the anomalous lateral diffusion induced by protein crowding. The three major results of this work are as follow. (1) The protein crowding significantly enhances the anomalous lateral diffusion in membranes. The subdiffusive lipid dynamics becomes much slower and spans up to several orders of magnitude longer than the typical crossover time of a few nanoseconds in a proteinpoor membrane. This observation is in line with the abovediscussed previous studies for crowded artificial or cellular membranes [49,55,58,100]. The protein channel motion observed in the plasma membranes of living cells also exhibits pronounced aging effects [55,58,63], a phenomenon not observed in our much simpler model membranes.
(2) The protein crowding effects a non-Gaussian character in the anomalous diffusion, which is incompatible with the known models describing lateral diffusion. Previously, several groups, including us, found that in protein-free lipid bilayers, at any lipid phases, the lateral diffusion is consistently governed by the Gaussian FLE process (or, equivalently, fractional Brownian motion [112][113][114]), which thus appears to be a universal platform for the description of lateral diffusion in membranes [52][53][54]. However, our current study demonstrates that the inclusion of proteins in the membranes drastically changes the diffusion character: the picture suggested earlier no longer holds under protein crowding, thereby implying that the current paradigm of anomalous diffusion in membranes has to be revised. We show in this work that in protein-rich membranes, regardless of the lipid species, the lateral diffusion of phospholipids and proteins becomes highly non-Gaussian, and can no longer be explained by the Gaussian FLE model. This is clearly seen in the propagator P r ðr; ΔÞ obtained from the simulated lipid trajectories, which obeys a compressed (or stretched) exponential form [Eq. (4)]. The analyses also show that the observed non-Gaussian character is not explained by some other lateral diffusion models, such as the anomalous diffusion in fractal space [32,115], the free-volume jump diffusion [116,117], or the nonergodic continuous time random walk model [55,89]. (3) Protein crowding creates significant spatiotemporal heterogeneity in lateral diffusion dynamics with a very long life span of > 1 μs (Fig. 7). An individual lipid, while diffusing across such an environment, undergoes a spatially varying diffusivity, stemming largely from slowed-down diffusion in the vicinity of membrane proteins [101] compared to fast diffusion in protein-free regions. This results in strong heterogeneity in the motion of individual lipids, as seen in the amplitude scatter of the time-averaged MSD curves shown in Fig. 8 and also in Fig. S4 in Supplemental Material [80], but does not give rise to nonergodic properties. Collecting these findings, it can be concluded that the protein crowding completely reshapes the lateral anomalous dynamics in membranes in its duration and stochastic character.
There are a few remarks on the non-Gaussian lateral diffusion observed in this work. Often, complicated diffusion patterns found in living cells have been modeled heuristically by a multicomponent Gaussian model [93,100]. Our computational study shows that this approach is not valid even in our crowded model membranes although they are much simpler than crowded cellular membranes. That is, the lateral diffusion in highly protein-crowded environments is not decomposed into the superposition of slow Gaussian diffusion in the proximity of proteins and fast Gaussian diffusion away from the proteins. This is mainly because a lipid molecule, whose dynamics is much faster than the protein dynamics, explores both slow and fast regions in space multiple times during its lateral diffusion and, thus, there is no distinction between slow and fast particles. While our study does not invalidate the use of multicomponent Gaussian approximation in general, it cautions the excessive interpretation of the results in terms of multiple-component Gaussian models when no further information is available.
Recently, a series of experimental and computational studies reported non-Gaussian diffusion dynamics in various soft complex systems [95][96][97][98][99]. Examples include colloidal beads on lipid bilayer tubes, particle diffusion in entangled actin networks, liposomes in nematic actin filaments, and colloidal suspensions. In these studies it was shown that the complex environments give rise to transient non-Gaussian diffusion on a certain time scale over which the propagator P typically has an exponential tail notwithstanding the diffusion dynamics maintaining the normal Einstein scaling law. On the one hand, the non-Gaussian lateral diffusion reported in this work is in line with these cases in that it is a phenomenon arising from the spatial complexity of the system. On the other hand, it is distinguished from them in that the non-Gaussian displacements are described by the compressed-exponential propagator [Eq. (4)] and lead to the anomalous diffusion dynamics δ 2 ∼ Δ α .
What physical origins govern the observed non-Gaussian anomalous diffusion present in the protein-crowded membranes? Is it a system-specific phenomenon only shown in the protein-crowded membrane or a universal character commonly valid for similar crowded systems? Through additional computational work, we find that a very similar non-Gaussian anomalous diffusion can take place in a much simpler crowded, quasi-two-dimensional Lennard-Jones (argon) system. Simulation results are summarized in Fig. 10. The inert atoms, otherwise displaying normal Brownian self-diffusion, undergo multifractal non-Gaussian anomalous diffusion when placed in an obstacle-crowded space where the immobile obstacles are aligned to give a transient confining effect to the diffusing particles. It is remarkable that an analogous profile of αðΔÞ over a long time window, similar to that observed in the protein-rich DPPC membrane (Fig. 2), is obtained when the obstacles are aligned in a similar fashion to the aggregating proteins [ Fig. 10(e)]. The propagators are also shown to have a similar non-Gaussian character as in the case of the crowded membranes studied [ Fig. 10(h)]. These results strongly suggest that the obstacle-induced hindrance and the transient confinement have a major responsibility for the non-Gaussian anomalous diffusion found in this work. We also find that the specific aggregation dynamics of proteins and the associated transient effect, as well as finite-size effects, are not significantly attributing to the non-Gaussian anomalous lipid dynamics. This is corroborated by the fact that the spatial correlation lengths and correlation times of proteins and lipids are considerably smaller than the membrane size and the simulation time (see Sec. S9 in Supplemental Material [80]). Our speculation on the origin of the observed non-Gaussian diffusion also gives an explanation of why the FLE Gaussian model, seemingly a universal dynamic model for the anomalous diffusion in protein-free membranes in all physical phases, fails to reproduce the diffusion characteristics in the protein-crowded membranes. This is because the observed non-Gaussianity is an effect of geometrical origin due to protein alignment, not from the viscoelasticity of a lipid membrane stemming from the lipid polymeric tail.
What dynamic model then replaces the Gaussian FLE model for the description of the non-Gaussian heterogeneous diffusion occurring in protein-crowded membranes? While it will be a challenging task to establish the quantitatively accurate model, we speculate that it will be based upon a hybrid model combining the obstructed diffusion with a diffusion process with space-dependent or fluctuating diffusivity. Several versions of the latter model were introduced recently where the local diffusivity is deterministically given by a specific functional form K 1 ðx; yÞ or randomly given with a probability density PðK 1 Þ [58,96,99,110,[118][119][120]. These models were shown to generate heterogeneous diffusion processes having pronounced fluctuations in the time-averaged MSDs and in some cases induce weak ergodicity breaking. It was demonstrated that the fluctuating diffusivity model describes the nonergodic subdiffusion of a DC-SIGN receptor in the living-cell membranes [58] and the reptation dynamics in entangled polymer systems [111]. As we learn from the diffusion map and the trace of the diffusivity fluctuation in Figs. 7 and 9, the heterogeneous lateral diffusion in our protein-crowded membranes may be understood within the framework of the fluctuating diffusivity model. While in the literature heavy-tailed distributions of the diffusivity and the sojourn time are usually introduced to explain a nonergodic subdiffusion [58,110,111], in our membrane systems the diffusivity fluctuation is shown to have a finite variance and bounded sojourn times. Such moderate diffusivity fluctuation induces heterogenous lateral diffusion without weak ergodicity breaking.
Expanding the idea introduced in Ref. [99], we consider a Gaussian subdiffusive process with a fluctuating diffusivity K α ðtÞ, where the evolution dynamics of K α is modeled by stochastic motion of a particle confined in the double-well potential landscape of − log PðK α Þ shown in Fig. 9. The physical scenario for this description is that the individual lipids experience an annealed fluctuating diffusivity induced by the protein crowding, while they carry out a FLE-like Gaussian subdiffusion if the presence of the obstacles (i.e., the proteins) is ignored. Then, in the long-time limit where K α ðtÞ reaches the stationary state with PðK α Þ, the average propagator for this heterogeneous process is given by It can be shown that a non-Gaussian propagator P is attained for the given form of P found in Fig. 9 (see Fig. S15 in Supplemental Material [80]). This means that the non-Gaussian, heterogeneous anomalous diffusion we find in our crowded membranes can be reasonably explained within the framework of the fluctuating diffusivity model. However, it turns out that the propagator Eq. (8) based on the Gaussian kernel does not successfully explain the obtained propagator for the lipids in the simulations (see Fig. S15 in Supplemental Material [80]). Previous studies [58,96,110,118,120] inform us that a multifractal anomalous diffusion is not seen in these space-dependent diffusion processes. It is speculated that the heterogeneous lateral diffusion attains the additional non-Gaussian, multifractal nature by the effect of obstruction in space. Thus, the proper diffusion model for the protein-crowded membranes should integrate the non-Gaussian feature associated with the obstruction in space into the above-discussed fluctuating diffusivity model. In line with this idea, we find that the propagator Eq. (8) with a non-Gaussian kernel stemming from the presence of obstacles shows reasonable agreement with the propagator for the lipids; see Fig. S16 in Supplemental Material [80]. A thorough theoretical investigation of this model remains for future work. Complementary stochastic analysis tools will be of crucial assistance in this task [35,[121][122][123]. Finally, we note in passing that our proposed diffusion mechanism for the protein-crowded membranes does not negate other possible diffusion mechanisms in other crowded artificial or living-cell membranes [124][125][126]. Various diffusion patterns markedly different from the current non-Gaussian heterogeneous diffusion can emerge due to the presence of additional complexity.
Is the non-Gaussian anomalous diffusion expected to be biologically relevant? There are many reasons to assume so. First, as cell membranes are crowded with proteins [81], their influence on membrane dynamics has to be accounted for. Second, a great fraction of cellular functions are due to proteins working in unison as protein dimers or higher oligomeric complexes, and in order to function, the proteins have to form a complex where the relative orientation and the distance of individual proteins renders the function possible. This implies that the proteins should sample their relative conformational space for sufficiently long times in order to find the functional structure for the protein complex, and this sampling process is obviously fostered by correlated motion of the proteins. In the same spirit, in agreement with the lipid raft paradigm [3], it is known that for a number of membrane proteins there are lipid binding sites [127] where a specific lipid binds to a protein in a manner where it is able to modulate protein conformation and dynamics, and thus also activation and function. Here, too, the formation of the protein-lipid complex may be a slow process where correlated motion of the lipids next to the protein would increase the formation rate of the functional protein-lipid unit.
As we learn from this study, a new component in complexity may significantly affect the diffusion dynamics in protein-crowded membranes, thereby reshaping the stochastic character of the motion. A systematic understanding of the role of various complex components in living cellular membranes for the lateral membrane dynamics will be a challenging task in the future.

ACKNOWLEDGMENTS
We thank Touko Apajalahti for the tool used to estimate the correlation of lipid motion. We acknowledge financial support from the Academy of Finland Here, we support our observations on the anomalous non-Gaussian lipid dynamics reported in the main text with additional simulations on a much simpler two-dimensional crowded system, i.e., a Lennard-Jones (LJ) fluid. The dynamics of this fluid are studied by simulating both obstacle-free [ Fig. 10(a)] and obstacle-containing [Figs. 10(b) and 10(c)] systems. In the obstacle-containing case, two different arrangements of the obstacles are considered in order to produce the effects of transient weak and strong confinement.
The obstacle-free system consists of 9000 particles, whereas in the two obstacle-containing systems 20 obstacles are placed among 1900 particles. Each of these obstacles is built of an aggregate of 55 particles that are permanently fixed in space.
The simulations are run in the NVT ensemble at the boiling point of argon (87.3 K), whose LJ parameters (ϵ ¼ 0.996 kJ=mol, σ ¼ 0.3405 nm [128]) are used for the particles. The area of the obstacle-free system is adjusted so that its 2D diffusion coefficients correspond to the value calculated for the same model in 3D at the same temperature and at a pressure of 1 bar. The areas of the two obstacle-containing systems are fixed to provide a very similar surface density and liquidlike behavior.
The obstacle-free system is simulated for 100 ns, whereas the obstacle-containing systems showing weak and strong confinement are simulated for 500 ns and 1 μs, respectively. All simulations are performed using an integration time step of 2 fs. The LJ interactions are cut off at 1 nm and a dispersion correction is applied to both energy and pressure. The temperature is kept constant at 87.3 K using the Berendsen thermostat [129]. Periodic boundary conditions are employed in the xy plane where diffusion takes place. The simulations are performed with the GROMACS simulation package version 4.5.x [130]. Figure 10 presents the highlights of the simulation results. The time-averaged MSD and the anomalous diffusion exponent α in Figs. 10(d) and 10(e) show that the obstacles in space induce anomalous diffusion. In the case of strong transient confinement [ Fig. 10(c)], mimicking the protein-crowded DPPC (aggregating) case in Fig. 2(b), the variation of α in time reproduces the behavior of α for the lipids in the crowded membrane. Further, the duration of anomalous diffusion in the argon system, as also in the DPPC membrane, is significantly elongated.
It is important to note that here the anomalous diffusion (consistent with the results we discuss in the main text) is observed for simple LJ particles with static obstacles. This suggests that the effect of protein-aggregating dynamics in Fig. 2(b) is not essential for the observed strong anomalous lipid diffusion. Once the obstacles give rise to transient but strong confinement, the diffusing particles carry out anomalous dynamics. Moreover, the analyses of the propagators for the three systems in Figs. 10(f)-10(h) show that the anomalous diffusion in the strongly confined case [ Fig. 10(c)] is indeed non-Gaussian. Analogously to the membrane system, P is described by a combination of two stretched-exponential propagators. Further, to rule out possible finite-size effects, we repeat this analysis for a LJ argon system whose size is 9 times (3 × 3) larger than the one discussed here. Consistent non-Gaussian diffusive behavior is observed (data not shown).
Summarizing, our simulation study on the simple LJ fluid systems corroborates the validity of the non-Gaussian anomalous lateral diffusion observed in protein-crowded membranes. The non-Gaussian lateral diffusion is not a system-specific, out-of-equilibrium transient dynamic property dependent upon the system preparation. Instead, it is a general dynamic equilibrium property induced by obstacle crowding.