Membrane microdomains emergence through non-homogeneous diffusion

  • Hédi A Soula1, 3Email author,

    Affiliated with

    • Antoine Coulon4 and

      Affiliated with

      • Guillaume Beslon2, 3

        Affiliated with

        BMC Biophysics20125:6

        DOI: 10.1186/2046-1682-5-6

        Received: 7 November 2011

        Accepted: 29 February 2012

        Published: 30 April 2012

        Abstract

        Background

        In the classical view, cell membrane proteins undergo isotropic random motion, that is a 2D Brownian diffusion that should result in an homogeneous distribution of concentration. It is, however, far from the reality: Membrane proteins can assemble into so-called microdomains (sometimes called lipid rafts) which also display a specific lipid composition. We propose a simple mechanism that is able to explain the colocalization of protein and lipid rafts.

        Results

        Using very simple mathematical models and particle simulations, we show that a variation of membrane viscosity directly leads to variation of the local concentration of diffusive particles. Since specific lipid phases in the membrane can account for diffusion variation, we show that, in such a situation, the freely diffusing proteins (or any other component) still undergo a Brownian motion but concentrate in areas of lower diffusion. The amount of this so-called overconcentration at equilibrium issimply related to the ratio of diffusion coefficients between zones of high and low diffusion. Expanding the model to include particle interaction, we show that inhomogeneous diffusion can impact particles clusterization as well. The clusters of particles were more numerous and appear for a lower value of interaction strength in the zones of low diffusion compared to zones of high diffusion.

        Conclusion

        Provided we assume stable viscosity heterogeneity in the membrane, our model propose a simple mechanism to explain particle concentration heterogeneity. It has also a non-trivial impact on density of particles when interaction is added. This could potentially have an impact on membrane chemical reactions and oligomerization.

        Keywords

        Membrane domains Non-homogeneous diffusion Individual-based model

        Background

        In the classical fluid-mosaic model of membrane, membrane components undergo isotropic random motion akin to Brownian motion [1, 2]. In this model, since specific interactions between individual molecules are not considered, resulting equilibrium distribution of components is homogeneous. Thus, the so-called Singer-Nicolson model implicitly assumes that no lateral heterogeneities within the membrane plane exist.

        Recently, this picture has considerably evolved towards a non-homogeneous distribution of cell-membranes components [35]. Experimental evidence has accumulated that cell membranes, particularly cell surface membranes, are indeed laterally heterogeneous on scales that range from tens of nanometers to a few microns. These heterogeneities are commonly referred to as “microdomains” to contrast them with the membrane macrodomains, the functionally differentiated surfaces of epithelial and other morphologically polarized cells. More evidence points towards the fact that some of these domains are enriched in various lipids such as cholesterol. These are sometimes called lipid rafts [6].

        The organization of membranes into microdomains can biologically be interesting because microdomains could strongly affect membrane functions as interacting species are likely to be in higher concentrations in the domains. Indeed a wide collection ofmembrane proteins involved in such processes have been shown to colocalize with rafts. They are thought to play an important role in various cellular processes such as trafficking and signaling to cite but a few [7, 8]. In addition, perturbations (like disruption by cholesterol removal) have important effects on cell responses compared to control [9, 10]. Surprisingly some recent works even suggest that the non-raftparts of the membrane are protein free [11]. It has been reported that membrane proteins with at least one transmembrane domain or with a hydrophobic modification are enriched in lipid rafts [12].

        It is not a surprise that various models have emerged to explain the existence and stability of membrane heterogeneity [13]. In particular, membrane protein diffusion alteration has been extensively studied in the context e.g oftrapping systems [1416], specific interactions [1719] or crowding [20, 21]. Membrane-skeleton fence or membrane-skeleton corralling models have also been proposed. Transmembrane proteins protrude into the cytoplasm and, in this model, their cytoplasmic domains collide with the membrane skeleton, inducing temporary confinement of the transmembrane proteins within the membrane-skeleton mesh [2224]. The transmembrane proteins then hop to an adjacent compartment. This model was supported by Monte Carlo simulation results [25].

        The latter model therefore rests on the assumption that rafts are well defined stable structure with mixed effects on diffusion. Indeed constrained diffusion within bounds imply an increase in concentration only temporarily since outbound protein have low probabilities to come back. In addition several other works suggest a more complicated picture [26]. Rafts (domains) are not well defined entities. Several measurements seem to indicate that some proteins are stably associated over long but finite period of time with discrete domains. These domains can either diffuse across the cell surface [27] or are immobile, such as cell surface caveolae [28]. Contrary to this view is the dynamic assembly. In this model proteins are transiently occupying raft domains and undergo diffusion both inside and outside the raft [2931].

        Being bound to a domain affects the protein diffusivity [32]. The same applies for lipids as their lateral mobility in a liquid-ordered (raft) phase is slower than in a liquid-disordered phase (non-raft) [29, 33, 34]. Indeed, several studies suggest that proteins and lipids undergo constrained and/or slowed diffusion within rafts [30, 31]. However, individual raft proteins do not appear to undergo correlated diffusion with one another [35].

        We propose in this article a simple mechanism which leads to the formation of protein/lipid enriched microdomains. This mechanism is based on non-homogeneous diffusion (NHD). Indeed, many structural constituents of the membrane can alter its viscous (or diffusion) properties. As such, a membrane cannot be approximated anymore by a simple 2D manifold with a constant diffusion but rather should be described by a diffusion profile: a space (position) dependent diffusion function. The membrane composition is therefore by itself a source of heterogeneity in the displacement of trafficking proteins.

        We will assume throughout this paper that the cell membrane has a non-spatially constant diffusion tensor that is temporally stable within the timeframe of the diffusion. These assumptions allow a greater generality since one never needs to assume any structurally stable component to define a domain. However we will not address the underlying mechanism behind such diffusion variation and simply assume it exists. Recently, several works have provided plausible mechanisms for viscosity alterations [3639].

        By solving the corresponding equation of motion, we show that the diffusion profile relates simply to the equilibrium concentration profile. In the case of punctual particles without interaction we show that this relation is extremely simple and we give a closed form of the equilibrium. In some cases we are able to provide a closed form for the whole solution and derive a FRAP (Fluorescence Recovery After Photobleaching) estimation of the diffusion in such conditions.

        We expand the model by performing simulations with non-punctual interacting particles. We show that the classical phase transition observed for this kind of system is altered by non-homogeneous diffusion – namely resulting in a shift in the transition diagram. This shift allows for extremely high concentration in the slow zones – much higher than for non-homogeneous diffusion or interaction separately – as well as higher clustering properties.

        Methods

        Continuous model

        In this first model, membrane particles are independent, punctual and subjected to Brownian motion. The membrane is either a one or two dimensional manifold with periodic boundary conditions. This imposes periodic boundary conditions on the equations of motion.

        For simplicity’s sake, we first derive the model and equations for the 1D case (equations for the 2D case are provided afterward using another method). Therefore, we consider the membrane as a segment of length 2 L centered on zero ([–L; L]); when a particle crosses one of the boundary it will appear at the other side). We assume that these particles undergo a classical Brownian motion in the position space (i.e. overdamped regime). To describe membrane heterogeneity, the diffusion coefficient depends on the position of the particle on the membrane. This so-called multiplicative noise can be efficiently modeled, in the 1D case, with the Stratonovich formalism yielding the following stochastic differential equation
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ1_HTML.gif
        (1)

        Here X t describes the position of the molecule at time t in the segment [–L; L], dZ/dt is the classical Brownian noise (with zero mean and unit variance) and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq1_HTML.gif a nonnegative 2L-periodic function. Here http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq2_HTML.gif is thespace-dependent standard deviation of the noise. Then, the actual diffusion (denoted as http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq3_HTML.gif ) is equal to http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq4_HTML.gif . We will suppose that http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq5_HTML.gif is continuously differentiable and > 0 everywhere.

        By setting ρ(x,t) the probability density function of a particle, we obtain the associated Fokker-Planck equation [40] derived from Eq. 1:
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ2_HTML.gif
        (2)

        Looking for continuous and differentiable solutions, boundary conditions emerge naturally from the 2 L periodicity. Namely, ρ(L,t) = ρ(–L,t) and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq6_HTML.gif for all t ≥ 0. Moreover, fluxes of particles are opposed on each side of the border. Writing this condition leads to http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq7_HTML.gif .

        The solution of Eq. 2 can be found at equilibrium, yielding the constant flux J:
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ3_HTML.gif
        (3)
        The boundary conditions impose J(–L) = –J(L) = J = 0 and consequently
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ4_HTML.gif
        (4)

        with http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq8_HTML.gif So, up to a normalization constant, equilibrium density is the inverse of the square root of the diffusion coefficient (or equal to the square root of viscosity).

        This result remains the same for 2D (or any higher dimension): Assuming a non-homogeneous brownian motion in higher dimension via
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ5_HTML.gif
        (5)
        in the square [–L; L]2 and dZ i /dt (i = 1,2) being the classical Brownian noise (with zero mean and unit variance). The Fokker-Planck equation for the density function ρ(x,y) will be
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ6_HTML.gif
        (6)
        with the boundary conditions
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ7_HTML.gif
        (7)
        for all t ≥ 0. To solve Eq. 6, we multiply by http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq9_HTML.gif and integrate with respect to both x and y to obtain http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq10_HTML.gif since http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq11_HTML.gif An integration by parts yields http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq12_HTML.gif . The first integral on the right hand side of the previous equation is taken over the square [–L; L]2 and vanishes due to toric boundary conditions (Eq. 7), yielding
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ8_HTML.gif
        (8)

        Since http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq13_HTML.gif is > 0 everywhere the left hand side of Eq. 8 is a decreasing function and its derivative must reach zero at equilibrium.

        This equilibrium therefore verifies
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ9_HTML.gif
        meaning the gradient of http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq14_HTML.gif is zero everywhere. That is
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ10_HTML.gif
        (9)

        The same result holds for any dimension: any diffusion over a non-constant diffusion profile yields, at equilibrium, to a non constant concentration namely its inverse. Note that as it should be, in the case of homogeneous diffusion ( http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq15_HTML.gif ) the equilibrium density (ρ(x) ≡ (2 L)–1) is independent of the diffusion coefficient D 0.

        As a consequence zones of slow diffusion (or high viscosity) will tend to gather more particles at equilibrium than any faster zones. Note that we used the Stratonovich formalism to describe non-homogeneous brownian motion of particles. This formalism make assumptions on the diffusion at the microscopic scale [41]. With other assumptions, we could have derived the same equations using the Ito formalism to obtain a slightly different result.

        Briefly: the associated Fokker-Planck would have been:
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ11_HTML.gif
        (10)
        yielding an equilibrium distribution
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ12_HTML.gif
        (11)

        The main results hold in a stronger way (due to the square) in any dimensions. In the continuous model, all simulations described below have also been done using the Ito formalism for the random motion scheme. All the results at equilibrium described below hold qualitatively in the Ito formalism but in a stronger way quantitatively. Thus, to stress our point, we will focus on the ‘weaker’ form of our results and remain in Stratonovich formalism throughout the rem`aining of the paper.

        The previous results hold for equilibrium distribution only. General transient solutions can be obtained for the one dimensional case by solving Eq. 2 entirely. Let http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq16_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq17_HTML.gif Eq. 2 becomes
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ13_HTML.gif
        (12)

        which is exactly the classical diffusion equation. One can immediately see that the same boundary conditions apply to p by replacing x by http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq18_HTML.gif . Equilibrium solution for classical diffusion on a circle is a constant density. So http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq19_HTML.gif . In addition, we obtain the transient solutions (on a circle) by replacing x by http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq20_HTML.gif in the general solutions of the classical diffusion. Note that depending on the initialconditions, particles can take a longer time to reach equilibrium in the slow diffusing zones. Therefore, transiently it may happen that zones of slow diffusion will gather less particles than in faster zones.

        On biological membranes, transient solutions are impossible to measure. The classical solution is to estimate them indirectly via imaging techniques such as Fluorescence Recovery After Photobleaching (FRAP). Bleaching particles amounts to create a new initial distribution ρ(x,0) = g(x). Particles inside the bleaching are removed (actually bleached) while the others are untouched: g(x) = 0 for http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq21_HTML.gif (for a bleaching beam of radius r centered on zero). In that case, the classical approximation of the half-time to recovery [42, 43] is
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ14_HTML.gif
        (13)
        where γ is a constant. Replacing r by http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq22_HTML.gif and using http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq23_HTML.gif , we can account for an heterogeneous diffusion dependent half-time recovery.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ15_HTML.gif
        (14)

        Equation 14 indicates that the time constant of the recovery is proportional to the square of the spatial average of the inverse of the square root of diffusivity. Using this equation, one can, theoretically, extracts the diffusion profile from FRAP data [44]. Unfortunately, for a general diffusion profile, no such transient solution can be found for higher dimensions.

        At this point, we can already make some useful predictions on membrane particle distribution. Indeed, some membrane components are known to alter membrane viscosity. Our model suggests that, around components that increase viscosity, diffusing particles will be in higher concentration – overconcentration – compared to faster zones. A good candidate could be cholesterol: Indeed, cholesterol is found in high amount on virtually all mammal cells and is known to rigidify membranes and cholesterol-enriched membranes display lower diffusion coefficient (by around one order of magnitude) [45, 46]. In this case our model predicts that membrane cholesterol islands (patches) are originating membrane domains. In other words and without any more assumptions, stable cholesterol domains should imply by means of slow diffusion a protein-riched domain.

        In addition, our model estimates this domains to concentrate more particles in proportion to the slowing. Indeed this overconcentration is simply proportional to the square root of the ratio of diffusion. For example, diffusion an order of magnitude slower should provide slow domains with around three times more particles.

        Conversely, when cholesterol concentration is low enough or high enough for the membrane to exhibit an approximatively constant diffusion profile, our model predicts no domain.

        Force field interactions

        This first model is neglecting some crucial features of membrane diffusing particles. As such, our model assumes independent particles that never collide nor hinder each other and that cannot interact. But both crowding and particle interactions are known to yield inhomogeneous particle distribution. We expand the previous model to take into account steric hindrance and possible particle interactions (e.g. protein-protein interactions) [47].

        To study the combined effect of non-homogeneous diffusion and particle interaction, we used 2D particle simulations. As such, particle interactions must be modeled at the same scale as the diffusion itself: the mesoscopic scale. In order to do so, particles were subjected to non-homogeneous Brownian motion and short-distance interactions [48, 49]. Assuming there are N particles, we introduce a new equation of motion for a particle j as:
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ16_HTML.gif
        (15)
        where u ij is the unitary vector of the semi-line starting from particle i to particle j, r ij being the distance between the two particles. The interaction force f is derived from the 6 – 12 Lennard-Jones empirical potential accounting for van der Waals interactions:
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ17_HTML.gif
        (16)

        the force f itself being the opposite of the gradient of the potential http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq24_HTML.gif . The two parameters and σ describe the shape of the interaction. For a given σ increasing increases the strength of the interaction (as well as its range) and clusters of particles appear. The average size ofclusters with increasing yields a phase transition from N clusters of size 1 to one cluster of size N[50, 51].

        However, in the case of non-homogeneous diffusion, the phase transition will be shifted to lower interaction strength in the slower zones. As such, fast zones will act as a particle provider (well) and slow zones as particles sink. We expect therefore to obtain in the slow zones a higher over-concentration compared to the fast zones but also a higher degree of clustering. This over-concentration should be higher than both effects (interaction and non-homogeneous diffusion) taken separately.

        Results and discussion

        Continuous model

        We present in this section the results obtained via simulation of particles that undergo equation of motion such as Eq. 1. We used the Milstein numerical schema for the Stratonovich formalism [52],
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Equ18_HTML.gif
        (17)

        Where Xi (t) is the position of particles i at time t (t k  = (k + 1)Δt) and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq25_HTML.gif and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq26_HTML.gif are the diffusion profile and its gradient respectively. Various diffusion profiles were tested as well as different number of particles. For the 1 dimensional experiment 100,000 particles were used whereas 1,000,000 were used for the 2 dimensional case. We tested various time step and took the highest value (Δt = 0.01)) that yielded no change in equilibrium.

        Equilibrium distribution behaved as expected. The first set of experiments is for the one dimensional case where we test various diffusion profiles (continuously differentiable 2 L-periodic). As an example, we present two of them: In one case, Figure 1, the diffusion profile is a Gaussian function http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq27_HTML.gif with D 0 = 1, d = 0.9 and ω = 0.1 (displayed in inset). The theoretical distribution http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq28_HTML.gif and the simulated particle distribution are depicted by a plain line and dots respectively.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig1_HTML.jpg
        Figure 1

        Results for 1 d experiment. Simulation (dots) versus theoretical (plain line) distribution of particles that undergo non-homogeneous diffusion. The diffusion profile is a simple inverse gaussian function (in inset) with equation: http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq29_HTML.gif and σ = 0.1 (theoretical:plain line, simulations: triangles).

        Transients for this simulation are displayed on Figure 2 as ρ(t,x) for 0 ≤ t ≤ 10 and x∈[–1:1]. The initial condition is ρ(0,x) = δ (x – 0.5) a Dirac function centered on 0.5.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig2_HTML.jpg
        Figure 2

        Transient distribution. Simulation of ρ(t,x) for t ∈ [0:5] and x ∈ [–1:1]. Setting ρ(0,x) = δ(x–0.5) as the initial condition. The diffusion profile is as in Figure 1 and therefore the distribution obtained near t = 5 is the same as the one at equilibrium described in Figure 1.

        In this case, due to the slow diffusivity the concentration at the peak of viscosity ρ(t,0.0) is the lowest value for a long period: transiently the slow zone is under populated.

        In addition, we simulated FRAP experiments and compared the diffusion coefficients computed through simulation to the theoretical ones as a function of the FRAP radius r (the FRAP experiment was centered on the slowest diffusion point). Figure 3 shows the results for three experiments, one classical homogeneous diffusion with http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq30_HTML.gif and two NHD with a centered gaussian http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq31_HTML.gif with d = 0.9 with ω = 0.1 and ω = 0.01 (shown in inset). Triangles, squares and circles are diffusion coefficient computed from the particle simulation for the cases http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq32_HTML.gif , http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq33_HTML.gif with ω = 0.1 and ω = 0.01 respectively. The lines are their theoretical counterparts computed using Eq. (14).
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig3_HTML.jpg
        Figure 3

        FRAP experiment. Results for FRAP experiments compared to theoretical prediction. Three experiments with different diffusion profiles are depicted. D = 1 (theoretical:plain line, simulations: triangles) and two gaussian profile http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq34_HTML.gif with σ = 0.01 (middle, theoretical: dotted line, simulations: squares) and σ = 0.1 (bottom, theoretical: dashes line, simulation: circles).

        Next, we present the results of NHD on a two-dimensional torus. The diffusion profile consists of 4 randomly positioned Gaussian patches centered on (x i ,y i ) ∈ [–L; L]2 with 1 ≤ i ≤ 4 and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq35_HTML.gif with D 0 = 1, d = 0.9 and ω = 0.1. The overall diffusion profile is the average http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq36_HTML.gif . Figure 4 shows the diffusion profile (light spectrum heightmap, Figure 4A) and the 2D histogram of the resulting equilibrium distribution of particles (Figure 4B). The last panel (4C) compares a cross-section (x = 0.2 is constant) of theoretical and simulated distributions. This cross-section can be visualized as a black line on the heightmap of Figure 4A and as a white line on Figure 4B.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig4_HTML.jpg
        Figure 4

        Results for 2 d experiment. A) The diffusion profile in light spectrum (toward red is slower) for 4 randomly “cholesterol” islands. B) 2 d histogram of simulated particles positions at equilibrium. C) Comparison with theoretical prediction for a cross section (x = 0.25). This section is represented by a black line on the heightmap and a white line on the 2 d histogram.

        In both one and two dimensional cases, the equilibrium distribution of particles that undergo non-homogeneous diffusion matches the theoretical prediction. Slower zones concentrate more particles than the others. In addition this over-concentration isproportional to the amount of slowing.

        Force field interactions

        Simulations were performed for the particles with interaction. To take into account steric hindrance, the numerical schema (Eq. 17) was modified by adding an interaction strength vector Δtf(X i (t)).

        Collision detection was included to allow a simple crowding experiment (i.e  = 0 in Eq. 16). We present the following experiment: On a 2D square (with toric boundary conditions) we insert a squared slow patch (on the upper leftcorner) where diffusion http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq37_HTML.gif is http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq38_HTML.gif times smaller than the diffusion ( http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq39_HTML.gif ) in the remaining of the space: a “cholesterol” patch. Another square of the same surface away from the cholesterol patch is used as control patch. Numbers of particles are counted in both patches at equilibrium and normalized by the total number of particles and the size of the patch.

        Assuming cholesterol diffusion is reduced by around one order of magnitude [5356], those numbers are displayed on Figure 5 for two conservative (i.e higher) values of the diffusion ratio http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq40_HTML.gif (black square) and r = 1.4 (black circles). The dashed lines are the theoretical predictions for the concentration in the patch in case of punctual particles (that is y = r). For low interaction strength, the punctual particles prediction remains a good approximation. However when the interaction increases, concentration of particles in the cholesterol patch increases dramatically while there is no such increase in the control patch. On Figure 6 screenshots of particles positions are displayed (r = 0.5) for  = 0 (A),  = 500(B),  = 4000(C) and  = 10000(D) – see also supporting videos.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig5_HTML.jpg
        Figure 5

        Particles interactions results. Concentration of particles in the “cholesterol” patch for two diffusion ratios http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq41_HTML.gif with r = 2 (black square) and r = 1.4 (black circle). Plain line is the control concentration with http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq42_HTML.gif . Dashed line is the theoretical NHD concentration.

        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig6_HTML.jpg
        Figure 6

        Particles maps. Position map of particles for A) ε = 0 and r = 2, B) ε = 500 and r = 2, C) ε = 4000 and D) ε = 10000. Note that the actual values of ε are scaled down by 10–12. See also supporting Additional file 1: Video S1, Additional file 2: Video S2, Additional file 3: Video S3 and Additional file 4: Video S4 for r∈{2,10} and ∈{0,1e5}.

        As expected, the overconcentration effect due to non-homogeneous diffusion is dramatically increased via particles interaction. As Figure 5 shows it, a great portion of the particles is concentrated in the cholesterol patch. In this extreme case, the patch is virtually composed of one big cluster of particles in interaction whereas there are no clusters in the remaining of the space.

        To numerically assess this situation, we compute the clustering coefficient as a function of the interaction strength: c(f) = N c /N i where N + N c is the number of clusters and N i is the number of particles involved. Clusters were computed by creating a graph of particles whose nodes are the particles themselves and link between two nodes exists if the distance between their centers is below ar 0 (a = 2.7 was chosen a bit higher than σ = 2.5). A value close to 1 indicates no clustering while a more pronounced clustering yields smaller values. Results are displayed in Figure 7. The values are computed in three situations: in the first two, diffusion is homogeneous with diffusion http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq43_HTML.gif (squares) and http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq44_HTML.gif (circles). The remaining case (triangles) is the one with a cholesterol patch with ratio r = 2 ( http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq45_HTML.gif as in the previous experiments). In all cases, the clustering coefficient decreases with the interaction strength and thereforeclusters appear (one needs higher interaction strength for the case http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq46_HTML.gif ). Moreover, the lower the viscosity the deeper the clustering. However, for high enough interaction strength, the cholesterol patch does not behave as an homogeneous slow zone: whereas cluster sizes seem to reach a plateau in the slow diffusion experiment, clusters were bigger in the cholesterol patch.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig7_HTML.jpg
        Figure 7

        Effects on clustering. Clustering coefficient (see text for definition) as a function of the interaction strength ε for control http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq47_HTML.gif (squares), slow experiment http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq48_HTML.gif (circles) and cholesterol square patch withration r = 2 (triangles).

        To mimic cholesterol effect, we deliberately chose diffusion ratios http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_IEq49_HTML.gif to be conservative i.e lower than one order magnitude of the diffusion ratio (r 2 ranges from 2 to 4 instead of 10) in almost all experiments [5356]. However, for completeness’ sake, we finally test several values of diffusion ratios. Results are summed up in Figure 8 for three values of the interaction strength (∈) and for ratios going from r = 1 (no slowing) to r = 10. Typical extremal dynamics are available on the supporting Additional file 1: Video S1, Additional file 2: Video S2, Additional file 3: Video S3 and Additional file 4: Video S4.
        http://static-content.springer.com/image/art%3A10.1186%2F2046-1682-5-6/MediaObjects/13628_2011_38_Fig8_HTML.jpg
        Figure 8

        Diffusion value dependence. Concentration ratio as a function of the ratio of diffusion for three interactions strength ε = 0 (triangles), ε = 1000 (squares), and, ε = 10000 (circles). Bold dashed line is the theoretical prediction of concentration for the case of no interaction (punctual particles) y = x. To guide the eye the line y = 1 is displayed.

        These results first confirm the predictions for punctual particles. With the parameters tested and as shown on Figure 8, the crowding only experiment did not differ much from the predictions. In addition, we show that particle interactionsin the form of Lennard-Jones potential can strongly modify both the magnitude and the nature of the resulting overconcentration. Around conservative values for diffusion, the so-called cholesterol patch can gather up to 9 times more particles than expected for an homogeneous medium. In addition, these particles tend to be more in clusters than any other situation. Our results can only be obtained and explained by combining both NHD and particle interactions. The cholesterol patch – the domain – spontaneously gathers particles and increases the stability of their interactions.

        Conclusion

        Cell membrane can display a wide variety of heterogeneity in its physical properties, viscosity and constituents composition. Our model describes the possible emergence of what we call a minimal membrane domain based on only one feature of the membrane, namely its local viscosity. We show that, everything else being equal, domain with high viscosity tends to gather statistically a larger proportion of the diffusive particles. Particles are not trapped within these domains but merely tend to spend more time in them. In addition, we show that the ratio of concentration of particles in and out domain is proportional to the square root of the inverse ratio of their local diffusion. As such, protein free membrane zones are seen as a product of high diffusivity and not from forbidding constraints – contrary to the hypothesis developed in [11] – and in accordance with the results found in [26]. Provided we assume that viscosity is indeed non constant for relevant time scales, this provides an extremely simple and parsimonious explanation for protein-enriched membrane domains.

        Moreover, while this model is fairly general, we emphasize the particular role of cholesterol. Cholesterol is ubiquitous in cell membranes and its influence on membrane diffusion has been well documented. Therefore, it is not surprising, from our point of view, that cholesterol has been found in higher quantity in virtually all microdomains exhibited so far. For many authors, its presence even characteries membrane domains [3, 11]. Our model shows that gradient in cholesterol concentration leads naturally to protein domains. This property depends only on the relative diffusion between cholesterol-rich zones and the remaining of the membrane.

        Looking for experimental results, authors have reported the impact of temperature and cholesterol depletion on membrane heterogeneities. Cholesterol removal should decrease the colocalization of membrane components as in [57, 58]. In addition, cooling decrease the diffusion ratio between membrane with and without cholesterol. For low temperature it even flips: our model predicts no colocalization or domains in the case of low temperature as it has been reportedin [58, 59].

        Important diffusion variation leads not only to a quantitative particles concentration variation but also to a qualitative variation. Protein interactions alter significantly the particles concentration landscape when combined with NHD. We show that cluster of particles are more stable in slower zones and appear for lower values of interactions. Our model predicts that domains of slow diffusion will alter affinity interactions. That is, for example oligomerization will happen preferentially inside these domains. This latter prediction has been observed on model membrane where modifications of the cholesterol content trigger oligomerization [60]. Thus, chemical reactions will have their equilibrium shifted in these zones, revealing a potential function for membrane microdomains.

        Declarations

        Acknowledgment

        This work has been partly supported by the Institut of Complex System of Rhones-Alphes (IXXI).

        Authors’ Affiliations

        (1)
        Université de Lyon Inserm UMR1060
        (2)
        Université de Lyon, CNRS INSA-Lyon, LIRIS
        (3)
        EPI BEAGLE, INRIA
        (4)
        Laboratory of Biological Modeling, NIDDK, NIH

        References

        1. Singer SJ, Nicolson GL: The fluid mosaic model of the structure of cell membranes. Science 1972,175(23):720–731.ADSView Article
        2. Saffman PG, Delbrück M: Brownian motion in biological membranes. Proc Natl Acad Sci USA 1975,72(8):3111–3113.ADSView Article
        3. Lingwood D, Simons K: Lipid rafts as a membrane-organizing principle. Science 2010,327(5961):46–50.ADSView Article
        4. Simons K, Ikonen E: Functional rafts in cell membranes. Nature 1997,387(6633):569–572.ADSView Article
        5. Simons K, Vaz WLC: Model systems, lipid rafts, and cell membranes. Annu Rev Biophys Biomol Struct 2004, 33:269–295.View Article
        6. Hancock JF: Lipid rafts: contentious only from simplistic standpoints. Nat Rev Mol Cell Biol 2006,7(6):456–462.View Article
        7. Kannan KB, Barlos D, Hauser CJ: Free cholesterol alters lipid raft structure and function regulating neutrophil Ca2+ entry and respiratory burst: correlations with calcium channel raft trafficking. J Immunol 2007,178(8):5253–5261.
        8. Hanzal-Bayer M, Hancock JF: Lipid rafts and membrane traffic. FEBS Lett 2007,581(11):2098–2104.View Article
        9. Nishimura SY, Vrljic M, Klein LO, McConnell HM, Moerner WE: Cholesterol depletion induces solid-like regions in the plasma membrane. Biophys J 2006,90(3):927–938.View Article
        10. Niemelä P, Hyvönen MT, Vattulainen I: Structure and dynamics of sphingomyelin bilayer: insight gained through systematic comparison to phosphatidylcholine. Biophys J 2004,87(5):2976–2989.View Article
        11. Lillemeier BF, Pfeiffer JR, Surviladze Z, Wilson BS, Davis MM: Plasma membrane-associated proteins are clustered into islands attached to the cytoskeleton. Proc Natl Acad Sci USA 2006,103(50):18992–18997.ADSView Article
        12. Foster LJ, Hoog CLD, Mann M: Unbiased quantitative proteomics of lipid rafts reveals high specificity for signaling factors. Proc Natl Acad Sci USA 2003,100(10):5813–5818.ADSView Article
        13. Burrage K, Hancock J, Leier A, Nicolau DV Jr: Modelling and simulation techniques for membrane biology. Brief Bioinform 2007,8(4):234–244.View Article
        14. Saxton MJ: Anomalous diffusion due to obstacles: a Monte Carlo study. Biophys J 1994,66(2 Pt 1):394–401.View Article
        15. Saxton MJ: Single-particle tracking: the distribution of diffusion coefficients. Biophys J 1997,72(4):1744–1753.View Article
        16. Saxton MJ: Anomalous subdiffusion in fluorescence photobleaching recovery: a Monte Carlo study. Biophys J 2001,81(4):2226–2240.View Article
        17. Nicolau DV Jr, Hancock JF, Burrage K: Sources of anomalous diffusion on cell membranes: a monte carlo study. Biophys J 2007, 92:1975–1987.View Article
        18. Nicolau DV Jr, Burrage K, Parton RG, Hancock JF: Identifying optimal lipid raft characteristics required to promote nanoscale protein-protein interactions on the plasma membrane. Mol Cell Biol 2006, 26:313–323.View Article
        19. Saxton MJ: Anomalous diffusion due to binding: a Monte Carlo study. Biophys J 1996,70(3):1250–1262.View Article
        20. Saxton MJ: Modeling 2D and 3D diffusion. Methods Mol Biol 2007, 400:295–321.View Article
        21. Berry H: Monte Carlo simulations of enzyme reactions in two dimensions: fractal kinetics and spatial. Biophys J 2002,83(4):1891–1901.View Article
        22. Kusumi A, Sako Y, Yamamoto M: Confined lateral diffusion of membrane receptors as studied by single particle tracking (nanovid microscopy)Effects of calcium-induced differentiation in cultured epithelial cells. Biophys J 1993,65(5):2021–2040.View Article
        23. Kusumi A, Nakada C, Ritchie K, Murase K, Suzuki K, Murakoshi H, Kasai RS, Kondo J, Fujiwara T: Paradigm shift of the plasma membrane concept from the two-dimensional continuum fluid to the partitioned fluid: high-speed single-molecule tracking of membrane molecules. Annu Rev Biophys Biomol Struct 2005, 34:351–378.View Article
        24. Kusumi A, Suzuki K: Toward understanding the dynamics of membrane-raft-based molecular interactions. Biochim Biophys Acta 2005,1746(3):234–251.View Article
        25. Niehaus AMS, Vlachos DG, Edwards JS, Plechac P, Tribe R: Microscopic simulation of membrane molecule diffusion on corralled membrane surfaces. Biophys J 2008,94(5):1551–1564.View Article
        26. Kenworthy AK, Nichols BJ, Remmert CL, Hendrix GM, Kumar M, Zimmerberg J, Lippincott-Schwartz J: Dynamics of putative raft-associated proteins at the cell surface. J Cell Biol 2004,165(5):735–746.View Article
        27. Pralle A, Keller P, Florin EL, Simons K, Hörber JK: Sphingolipid-cholesterol rafts diffuse as small entities in the plasma membrane of mammalian cells. J Cell Biol 2000,148(5):997–1008.View Article
        28. Thomsen P, Roepstorff K, Stahlhut M, van Deurs B: Caveolae are highly immobile plasma membrane microdomains, which are not involved in constitutive endocytic trafficking. Mol Biol Cell 2002, 13:238–250.View Article
        29. Dietrich C, Yang B, Fujiwara T, Kusumi A, Jacobson K: Relationship of lipid rafts to transient confinement zones detected by single particle tracking. Biophys J 2002,82(1 Pt 1):274–284.View Article
        30. Niv H, Gutman O, Kloog Y, Henis YI: Activated K-Ras and H-Ras display different interactions with saturable nonraft sites at the surface of live cells. J Cell Biol 2002,157(5):865–872.View Article
        31. Shvartsman DE, Kotler M, Tall RD, Roth MG, Henis YI: Differently anchored influenza hemagglutinin mutants display distinct interaction dynamics with mutual rafts. J Cell Biol 2003,163(4):879–888.View Article
        32. Wang Q, Zhang X, Zhang L, He F, Zhang G, Jamrich M, Wensel TG: Activation-dependent hindrance of photoreceptor G protein diffusion by lipid microdomains. J Biol Chem 2008,283(44):30015–30024.View Article
        33. Dietrich C, Bagatolli LA, Volovyk ZN, Thompson NL, Levi M, Jacobson K, Gratton E: Lipid rafts reconstituted in model membranes. Biophys J 2001,80(3):1417–1428.View Article
        34. Kahya N, Scherfeld D, Bacia K, Poolman B, Schwille P: Probing lipid mobility of raft-exhibiting model membranes by fluorescence correlation spectroscopy. J Biol Chem 2003,278(30):28109–28115.View Article
        35. Vrljic M, Nishimura SY, Brasselet S, Moerner WE, McConnell HM: Translational diffusion of individual class II MHC membrane proteins in cells. Biophys J 2002,83(5):2681–2692.View Article
        36. Khelashvili G, Weinstein H, Harries D: Protein diffusion on charged membranes: a dynamic mean-field model describes time evolution and lipid reorganization. Biophys J 2008,94(7):2580–2597.View Article
        37. Lee DW, Min Y, Dhar P, Ramachandran A, Israelachvili JN, Zasadzinski JA: Relating domain size distribution to line tension and molecular dipole density in model cytoplasmic myelin lipid monolayers. Proc Natl Acad Sci USA 2011,108(23):9425–9430.ADSView Article
        38. Brewster R, Safran SA: Line active hybrid lipids determine domain size in phase separation of saturated and unsaturated lipids. Biophys J 2010,98(6):L21-L23.View Article
        39. Machta BB, Papanikolaou S, Sethna JP, Veatch SL: Minimal model of plasma membrane heterogeneity requires coupling cortical actin to criticality. Biophys J 2011,100(7):1668–1677.View Article
        40. Risken H: The Fokker-Planck Equation: Methods of Solution and Applications. Berlin: Springer; 1989.MATHView Article
        41. Schnitzer M: Theory of continuum random walks and application to chemotaxis. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 1993,48(4):2553–2568.MathSciNetView Article
        42. Reits EA, Neefjes JJ: From fixed to FRAP: measuring protein mobility and activity in living cells. Nat Cell Biol 2001,3(6):E145-E147.View Article
        43. Lopez A, Dupou L, Altibelli A, Trotard J, Tocanne JF: Fluorescence recovery after photobleaching (FRAP) experiments under conditions of uniform disk illumination. Critical comparison of analytical solutions, and a new mathematical method for calculation of diffusion coefficient D. Biophys J 1988,53(6):963–970.View Article
        44. Siggia ED, Lippincott-Schwartz J, Bekiranov S: Diffusion in inhomogeneous media: theory and simulations applied to whole cell photobleach recovery. Biophys J 2000,79(4):1761–1770.View Article
        45. Lindblom G: Johansson LB. Arvidson G: Effect of cholesterol in membranes. Pulsed nuclear magnetic resonance measurements of lipid lateral diffusion. Biochemistry 1981,20(8):2204–2207.
        46. Owicki JC, Mcconnell HM: Lateral diffusion in inhomogeneous membranes. Model membranes containing cholesterol. Biophys J 1980,30(3):383–397.View Article
        47. Mayawala K, Vlachos DG, Edwards JS: Spatial modeling of dimerization reaction dynamics in the plasma membrane: Monte Carlo vs. continuum differential equations. Biophys Chem 2006,121(3):194–208.View Article
        48. Tucci K, Kapral R: Mesoscopic model for diffusion-influenced reaction dynamics. J Chem Phys 2004,120(17):8262–8270.ADSView Article
        49. Tucci K, Kapral R: Mesoscopic multiparticle collision dynamics of reaction–diffusion fronts. J Phys Chem B 2005,109(45):21300–21304.View Article
        50. Soula HA, Robardet C, Perrin F, Gripon S, Beslon G, Gandrillon O: Modeling the emergence of multi-protein dynamic structures by principles of self-organization through the use of 3DSpi, a multi-agent-based software. BMC Bioinforma 2005, 6:228.View Article
        51. Zeng X, Oxtoby D: Gas–liquid nucleation in Lennard-Jones fluids. J Chem Phys 1991,94(6):4442.ADSView Article
        52. Kloeden PE, Platen E: Numerical solution of stochastic differential equations. Berlin: Springer; 1999.
        53. Niemelä PS, Ollila S, Hyvönen MT, Karttunen M, Vattulainen I: Assessing the nature of lipid raft membranes. PLoS Comput Biol 2007,3(2):e34.ADSView Article
        54. Meer GV, Voelker DR, Feigenson GW: Membrane lipids: where they are and how they behave. Nat Rev Mol Cell Biol 2008,9(2):112–124.View Article
        55. Lenaz G: Lipid fluidity and membrane protein dynamics. Biosci Rep 1987,7(11):823–837.View Article
        56. Miersch S, Espey MG, Chaube R, Akarca A, Tweten R, Ananvoranich S, Mutus B: Plasma membrane cholesterol content affects nitric oxide diffusion dynamics and signaling. J Biol Chem 2008,283(27):18513–18521.View Article
        57. Parpal S, Karlsson M, Thorn H, Strålfors P: Cholesterol depletion disrupts caveolae and insulin receptor signaling for metabolic control via insulin receptor substrate-1, but not for mitogen-activated protein kinase control. J Biol Chem 2001,276(13):9670–9678.View Article
        58. Fujita A, Cheng J, Hirakawa M, Furukawa K, Kusunoki S, Fujimoto T: Gangliosides GM1 and GM3 in the living cell membrane form clusters susceptible to cholesterol depletion and chilling. Mol Biol Cell 2007,18(6):2112–2122.View Article
        59. Fretz MM, Penning NA, Al-Taei S, Futaki S, Takeuchi T, Nakase I, Storm G, Jones AT: Temperature-, concentration- and cholesterol-dependent translocation of L- and D-octa-arginine across the plasma and nuclear membrane of CD34+ leukaemia cells. Biochem J 2007,403(2):335–342.View Article
        60. Ishitsuka R, Kobayashi T: Cholesterol and lipid/protein ratio control the oligomerization of a sphingomyelin-specific toxin, lysenin. Biochemistry 2007,46(6):1495–1502.View Article

        Copyright

        © Soula et al.; licensee BioMed Central Ltd. 2012

        This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://​creativecommons.​org/​licenses/​by/​2.​0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.