Size-Dependent Diffusion of Membrane Inclusions

INTRODUCTION

Diffusion is (he basic means of transport in membranes, and the mixing of membrane components via Brownian motion is often very efficient: a lipid can rotate about its own axis within <10µs(1) and explores ~ 16µm^sup 2^ within a second (2). Given that speedy movement, membrane-anchored reaction partners can meet and align via diffusion and. upon doing so, create new (active) complexes or trigger signaling events downstream.

In the spirit of the above, biological membranes for a long time have been viewed as an unstructured, two-dimensional fluid into which individual proteins are embedded ("fluid mosaic model" (3)). Yet, more recently this picture was replaced by a more structured picture of cellular membranes (see Engclman (4) and Simons and lkonen (5) for review). Several lines of evidence have been given that membrane domains ("rafts"), consisting of lipids and/or clusters of membrane proteins, compartmenlalize in particular the plasma membrane. In fact, the size of raft-like inclusions in hiomemhranes has been reported to cover a wide range, from a few nanometers to some KM) nm (o). Bearing this in mind, the si/e-dependent diffusive mobility of membrane inclusions (from single proteins/lipids to raft-like domains) becomes a topical and important issue, even more so as measurements of diffusion coeflicients are frequently used to determine complex formation via extracting the si/e of the tracked object.

Several experimental studies have given support to Eq. 1 and Eq. 2 (e.g., (9,K))), whereas very recent experiments have indicated strong deviations from Eq. I tor small and intermediate radii (11) (see also Discussion). In fact, a rigorous experimental test of the predicted size-dependences in Eq. I and Eq. 2 was and is very challenging due to a lack of appropriately sizable inclusions, systematic limitations in recording the diffusive movement (e.g.. problems with photobleaching protocols (12)). unavoidable membrane undulations, etc.. which perturb the measurement and increase the error bars of the determined diffusive mobility. This level of uncertainty in experiments underlines the importance of a comprehensive test of Eq. 1 (and Eq. 2) by alternative means. The need for a quantitative test is further highlighted by the fact that Eq. 1 spuriously predicts negative diffusion coefficients for h?^sub m^ < R?^sub c^. In fact, extensive liydrodynamic calculations have predicted that Eq. I only holds for small radii whereas for large radii a scaling D ~ 1/R should emerge (13). It is worthwhile to note that the latter prediction as well as the derivation of Eq. 1 are based on incompressible, cylindrical membrane inclusions surrounded by incompressible fluids and the assumption of no-slip boundary conditions. All of these assumptions, although valid on the macroscopic scale, may not hold true on the meso- and nanoscale: water as well as lipid bilayers have a finite compressibility; nanoscopic membrane inclusions and larger, oligomeric (raft-like) structures can be expected to have internal degrees of freedom that reduce their lateral mobility; and finally, evidence has been given that the stick boundary condition is in general not appropriate on the molecular scale (14).

Here, we have used mesoscopic simulations of lipid bilayers with embedded, transmembrane inclusions to study the validity of Eq. 1 and Eq. 2 over a wide range of radii. We find that Eq. 1 gives a good quantitative description of the lateral diffusion coefficient up to a critical radius R^sub c^ [asymptotically =] h?^sub m^/ (2?^sub c^). Beyond this radius, the numerically determined diffusion coefficients strongly deviate from the SaffmanDelbrück relation and a new scaling D ~ 1/R~ emerges. We give theoretical arguments that the latter arises due to the combination of the asymptotic hydrodynamic drag and internal degrees of freedom that are anticipated to become relevant for large, raft-like inclusions, especially on short timescales. The rotational diffusion coefficient on the contrary is well described by Eq. 2 over the entire range of tested radii.

METHODS

Simulation details

We have set the interaction cutoff of lhe ihominstal r^sub 0^, lhe head mass m (all heads were assumed to have the same mass), and the thermostat temperature k^sub B^T to unity and used these parameters as basic units. We further have chosen the dissipation and noise parameters for all beads to be s^sub ij^ = s = 3 and ?^sub ij^ = ? = 9/2. For the cxplieit solvent model, we have chosen N = 4 beads per lipid chain, and the remaining interaction constants were chosen in accordance with Laradji and Kumar (19) (indices W, H. T = water, hydrophilic. hydrophobic head): a^sub HH^ = a^sub ww^ = a^sub TT^ = a^sub HW^ = 25 k^sub B^T, a^sub HT^ = a^sub WT^ = 200 k^sub B^T, k^sub 2^ = 100 k^sub B^T, = 0.45 r^sub 0^. The interaction constants for the implicit solvent model were: a^sub HH^ = a^sub HT^ = 48 k^sub B^T. a^sub TT^ = 96 k^sub B^T, w = r^sub 0^, e = 1.4 k^sub B^T, k^sub 2^ = 120 k^sub B^T/r^sup 2^^sub 0^, = 20 k^sub B^T, l^sub 0^ = 0.6 r^sub 0^. For small inclusions, the linear size of the simulation box in the plane of the bilayer was chosen as L = 40 r^sub 0^; for large inclusions, L was chosen to be at least fourfold bigger than the inclusion's diameter. In all cases, the height of the simulation box was lixed to 16 r^sub 0^. which is about fourfold bigger than the membrane thickness.

We have integrated the equations of motion with a velocity Verlet scheme (20) (time increment ?t = 0.01) and imposed periodic boundary conditions. During the initial relaxation of the membrane, we used a harosiat that has been adapted for the use with DPD (21 ) to achieve a tensionless membrane. In all simulations, we lirst relaxed the membrane with a single inclusion for at least 1.5.10^sup 4^ time steps. During this lime, the barostat was used to achieve a tensionless bilayer. From the tinal slate, we iterated 10^sup 6^ time steps during which we tracked the center-of-mass position of the inclusion and its orientation (see arrow in Fig. 1 b).

Conversion to Sl units

For the implicit solvent approach, we related our data to SI units by choosing the length scale as r^sub 0^ [Lef-right arrow] 1 nm, which yielded a bilayer thickness h [asymptotically =] 3.3 nm similar to synthetic membranes (22). This value was obtained by averaging the distance of all hydrophilic lipid head beads in the opposing leaflets of the hilayer ("phosphate-lo-phosphate distance"). The internal limescale was determined hy comparing lhe numerically obtained diffusion coeflicient of a single lipid with experimentally measured values (2) D [asymptotically =] 3 . 10^sup -2^ [Lef-right arrow] 4µm^sup 2^/s (see also (23)). A single lime step (?t = 0.01) corresponded Io a real time of ~80 ps. i.e.. the total time simulated was in all cases 80 nm. whereas the typical membrane patch was ~80 nm X 80 nm. In lhe explicit solvent case, r^sub 0^ and dt = 0.01 corresponded to 1.1 nm and 97 ps, respectively. Using the mentioned conversion to Sl units, we determined here also the viscosity ?^sub c^ of the pure solvenl by monitoring the diffusion coefficient of differently sized cylinders with hexagonal cross section (diameter 2 k: + I beads: length 2 k beads: head-bead distance l^sub 0^, = 0.45 r^sub 0^) in a "water box" with particle density Q = 3/r^sup 3^^sub 0^. For a single solvent bead, a radius R^sub 0^ = (1/Q)^sup 1/3^ was assumed and the hydrodynamic radii of the diffusing cylinders were set Io R^sub h^ = kl^sub 0^ + R^sub 0^/2. which was ~15% larger than lhe radius of gyration. From the Einstein-Stokes equation (modified wilh the contribution by internal modes, sec Results and Discussion), we determined the viscosity of the solvent to be ?^sub c^ [asymptotically =] 0.03 Pa s. Bearing in mind lhe somewhat vaguely defined radii R^sub h^ (due to lhe use of soft-core potentials) and the uncertainly if slick or slip boundary conditions are more appropriate, the value for ?^sub c^ may be slightly higher or lower.

Data evaluation

We briefly derive Eq. 3. Starting with a panicle at r = 0 at lime t = 0, one obtain from the two-dimensional diffusion equation lhe probability of finding the particle in an infinitesimal area element dA = 2 prdr around the locus r at time t as p(r)dA = 2 exp(~|r|^sup 2^/(4/Dt)/(4Dt)rdr. Changing variables to the quadratic distance ? = r^sup 2^, one obtains the differential distribution of squared increments as p(?) = exp(-?/(4Dt))/(4Dt)d?. For the purpose of fitting, the integrated distribution P(?x^sup 2^) = ?^sup ?x^sup 2^^^sub 0^ p(?)d?(Eq. .1) is more convenient as it does not suffer from the choice of the bin size. This approach has also been applied successfully in single-particle tracking studies (24).

RESULTS AND DISCUSSION

To study the diffusion of large inclusions in a self-assembling lipid bilayer, we have used efficient mesoscopic simulations. We used two related coarse-grained molecular dynamics simulation methods that belong to the class of DPD schemes (15-17) (see Methods for simulation details): an explicit solvent model ("standard DPD") and an implicit solvent approach. The latter has recently been studied in some detail ( 18,25,26) and yields a very efficient way to simulate large membranes that arc virtually untractablc by standard DPD. In the remainder, we will thus concentrate on the implicit-solvent model as this approach allowed us to investigate much larger inclusions than with an explicit solvent approach. We will however compare the results to those obtained with the explicit solvent whenever possible.

In the simulations, individual lipids were considered as chains of N = 3 beads connected by Hookean springs with each head representing a number of atoms (e.g., several methyl groups). The first bead represented the hydrophilic headgroup. the N - 1 consecutive beads represented the hydrophobic tail (cf. Fig. 1 a). The chain was given a bending rigidity, i.e., a straight chain was energetically preferred. Inclusions were modeled by cylinders with a length of 2 N heads and a hexagonal cross section (cf. Fig. 1 b). The first and last bead in each chain were taken to be hydrophilic. the remaining N - 2 were taken to be hydrophobic. In correspondence with the lipids. the beads in each chain were connected by Hookean springs and the chain was given a bending rigidity. These 2 N head chains were positioned on all inner vertices of a plane hexagon with edge length K + 1 and were connected laycrwise by Hookean springs. In total, the inclusion consisted of 2 N{3 K(K + 1) + 1} beads.

All beads interacted via a pairwise soft-repulsive potential, where the strength of repulsion was tuned to he stronger between hydrophobic-hydmphilic pairs. All beads were further subject to a Galilean-invariant, momentum-conserving DPD thermostat that included dissipative and random forces (see Methods). The solvent-induced attraction of the lipids was mimicked by an attractive pairwise potential among the hydrophobic heads in agreement with Cookc et al. (1K). For comparison with Eq. I and Eq. 2, we have transferred the simulation units to SI units (see Methods).

In Fig. 2, the numerically obtained data for two representative inclusions (K = 2. 19) are shown together with the host lit according to Eq. 3. A clear shift to smaller quadratic distances ?x^sup 2^ is visible for the larger inclusion, highlighting the reduced diffusive mobility. We next determined systematically the lateral diffusion coefficient D for inclusions of various sizes. For comparison with Eq. 1, we assigned each hexagon a radius R = l^sup 0^(K + 1) = (K + 1).0.6 nm (cf. Methods). The dependence of D on R is shown in Fig. 3 together with the best fit according to Eq. 1. from which we obtain the viscosity ?^sub m^ [approximate] 0.25 Pa s of the bilayer (via h [approximate] 3.3 nm, Methods). This value is in good agreement with typical data from the literature (22). For the simulations with explicit solvent, we also find a very good agreement with Eq. 1 for small radii (Fig. 3. inset) and the determined parameters (h [asymptotically =] 3.5 ?^sub m^ [asymptotically =] 0.19 Pa s) agree well with those obtained for the implicit solvent approach. From the fit, we further obtain via h?^sub m^/?^sub c^ [asymptotically =] 14.7 nm the effective viscosity ?^sub c^ [asymptotically =] 0.056 Pa s of the surrounding (explicit solvent: ?^sub c^ [asymptotically =] 0.039 Pa s). Both values for ?^sub c^ correspond well to the independently determined viscosity of the solvent (see Methods) albeit ?^sub c^ is not a well-delined quantity in the implicit-solvent approach (see Discussion).

Although for small radii the Suffman-Delbrück relation yields a very good description, strong deviations are visible beyond a critical radius R^sub c^ [asymptotically =] h?^sub m^/(2?^sub c^) [asymptotically =] 7.4 nm. (Fitting the entire numerical data with Eq. 1 alone results in a very bad description (data not shown.)) We would like to note that the critical radius R^sub c^ emerges naturally here when comparing the flux of energy dissipated by the bilayer and the solvent, respectively: whereas the former is. Jm ? 2 pRh?^sub m^, the latter is given by J^sub c^ ? 2 pR^sup 2^?^sub c^, i.e., a crossover is expected at K = 2 R^sub c^ beyond which the friction due to the solvent-facing area dominates. For R [much greater than] R^sub c^, the problem can thus be reduced to the edgewise motion of a thin disk in a fluid of viscosity ?^sub c^ for which one finds (with appropriate prefactors) D = k^sub B^T/(16 R?^sub c^) (13).

We finally monitored the size dependence of the rotational diffusion coefficient D^sub r^. For all tested inclusion sizes, the decrease of D^sub r^, is well described by Eq. 2 (Fig. 4). The prefactor of the best fit to the data yields k^sub B^T/(4p?^sub m^h) [asymptotically =] 0.45 µm^sup 2^/s, which is in excellent agreement with the value 0.41 µm^sup 2^/s obtained from h and ?^sub m^ as determined by fitting the lateral diffusion coefficient for small radii. Similarly, we found good agreement with Eq. 2 for the data obtained with the explicit-solvent model (Fig. 4, inset).

DISCUSSION

In conclusion, we have shown by means of extensive simulations that the rotational diffusion coefficient of membrane inclusions follows indeed the predicted form. Eq. 2. The lateral diffusion coefficient, however, does only follow the Salfman-Delbrück relation. Eq. I. for small radii, whereas for large radii we find substantial deviations and an asymptotic scaling according to Eq. 4. The latter takes into account the internal degrees of freedom of the inclusion, which should arise naturally as larger inclusions are most likely loosely associated protein oligomcrs and/or rail-like entities. The proposed scaling Eq. 4 is thus expected to he more realistic than the results derived for large incompressible cylinders. Although the internal modes can be expected to subside when following the diffusion trajectory over asymptotically long times, they clearly contribute significantly to the diffusion on small and intermediate timescales and may thus be accessible experimentally, e.g.. by fluorescence correlation spectroseopy. Small deviations from Eq. 1 are also expected (and observed, cf. Fig. 3) for radii R < 1 nm, since in this regime the discrete composition of the membrane from individual lipids must be taken into account ("free-volume model" (27,28)).

It is worthwhile to note that we found similar results by two simulation approaches (i.e.. using an implicit and an explicit solvent) that also differed in the type of lipids (implicit solvent N = 3 beads per lipid, explicit solvent N = 4). We furthermore have used the approach of Shillcock and Lipowsky (29) (i.e., a lipid with N = 7) and did find similar results for the diffusion coefficient (D = 1 µm^sup 2^/s, ?^sub m^ [asymptotically =] 0.2 Pa s). We are therefore confident that the lipid model inlluences the presented results only weakly, e.g., by slightly altering the value of ?^sub m^.

At first glance it is surprising that the implicit-solvent approach, where ?^sub c^ = O by definition, can reproduce hydrodynamic relations like Eq. 1 in which a finite value for ?c^sub ^ is needed. The reason for this can be traced back again to the internal modes of the inclusion and the dissipative forces imposed by the thermostat. The erratic impact of the surrounding lipids excites shear modes within the inclusion with a polarization perpendicular to the bilayer normal, and these modes are dissipated by the action of the thermostat. A "neutral layer" of the inclusion located roughly in the midplane of the bilayer therefore feels a friction with respect to the layers that lie above and below in-plane with the hydrophilic headgroups of the lipids. Hence, these shear modes mimic an apparent solvent viscosity that should change, when the dissipation strength ? in the thermostat is altered for the beads within the inclusion. Indeed, we have observed that a reduction of ? within the inclusions leads to an effective reduction of ?^sub c^ (data not shown). Efforts to thoroughly quantify the nontrivial connection between the emergence of an apparent solvent viscosity and the inner modes beyond this qualitative argument as well as an investigation that hydrodynamie quantities can be faithfully reproduced by an implicit-solvent scheme are currently under way.

In a recent report, strong deviations from Eq. 1 that indicate a scaling D ~ 1/R have been found experimentally (11). In fact, the authors claim deviations already for radii R [asymptotically =] I nm, which is in contrast to earlier studies (9,10). We would like to point out here that substantial deviations from Eq. 1 should be possible only for R = 1 nm (due to the free-volume model, which also predicts D ~ 1/R for sufficiently small R (28)) or for R > R^sub c^ [asymptotically =] 10 nm. Bearing in mind the unavoidable systematic limitations when experimentally assessing the diffusive mobility and given that D ~ 1/R can provide a reasonable lit to Eq. 1 when error bars are big enough, it is likely that the signatures of Eq. 1 have been masked in the data presented in Gumbin et al. (11). Nevertheless, it will be interesting to revisit the approach of Gambin et al., i.e.. i), to also use a complementary experimental techniques (e.g., fluorescence correlation spectroscopy) and ii), to extend the study to larger radii and higher temporal resolution where Eq. 4 can be expected to become visible.

We thank R. Bminsma. T. Liverpool. U. Schwarz. and U. Seiferl for helpful discussions. This work was supported by the Institute for Modeling and Simulation in the Bioseiences (BIOMS) in Heidelberg.