An effective model for the cross-over from anomalous to Brownian dynamics of proteins based on a viscoelastic membrane description
AAn effective model for the cross-over from anomalous to Brownian dynamics of proteins based on aviscoelastic membrane description
Loris Di Cairano ∗ Department of Physics, Faculty of Mathematics, Computer Science and Natural Sciences, Aachen University, 52062 Aachen, Germany andComputational Biomedicine, Institute of Neuroscience and Medicine INM-9 and Institutefor Advanced Simulations IAS-5, Forschungszentrum Jülich, 52428 Jülich, Germany
Benjamin Stamm † Applied and Computational Mathematics, Department of Mathematics, RWTH Aachen University, Aachen, Germany
Vania Calandrini ‡ Computational Biomedicine, Institute of Neuroscience and Medicine INM-9 and Institutefor Advanced Simulations IAS-5, Forschungszentrum Jülich, 52428 Jülich, Germany (Dated: November 20, 2020)In this paper we derive a model to describe the effective motion of a protein laterally diffusing in a lipidmembrane. By postulating that the lipid membrane is a linear viscoelastic fluid and the protein a rigid body, wederive a continuum description of the system. Within this framework, the lipid membrane is modeled through thelinearized Navier-Stokes equations complemented by a non Newtonian constitutive equation, and the protein bya (modified) Euler’s law, where a noise term is introduced to account for the stochasticity of the system. Then,by eliminating the lipid membrane degrees of freedom, we obtain a Generalized Langevin Equation (GLE)describing the diffusive motion of the protein with a memory kernel that can be related to the response ofthe membrane encoded in the solution of the constitutive equation. By representing the viscoelastic behaviorthrough the Prabhakar fractional derivative, one generates a memory kernel containing a three-parameter Mittag-Leffler function. An additional Dirac delta-function term in the memory kernel accounts for the instantaneouscomponent of the system response.A comparison between the Mean Squared Displacement (MSD) derived in the framework of this model andthe MSD of a protein diffusing in a membrane calculated through molecular dynamics (MD) simulations showsthat the proposed model accurately describes all regimes of the diffusive process, namely, the ballistic, thesubdiffusive (anomalous) and the Brownian one, as well as the transitions from one regime to another. Themodel further provides an estimation of the membrane viscosity, which is of the order of magnitude of thevalues found by rheology experiments on analogous systems.
INTRODUCTION
Lateral diffusion in membrane is key for cellular infor-mation processing [1]. Cell membrane fluidity determineslipid and protein mixing, thus regulating diffusion-limitedbiochemical interaction rates responsible for signal transduc-tion from the extracellular to the intracellular environment.In the last years, it has been studied how biological struc-tures and features of living cells, such as membrane com-position, concentration of proteins in membrane, compart-mentalization and crowding, affect the diffusion of protein orlipids in membrane [2, 3]. A plethora of experimental studiesin cellular membranes [4–7], and in vitro lipid bilayers [8–10], as well as computer simulations of minimalistic modelmembranes in crowding conditions [11–15] have shown a de-viation from the simple Brownian diffusion, where the ran-dom displacements are described by a Gaussian probabilitydistribution and the Mean Squared Displacement (MSD) in-creases linearly in time when time-lags much longer than thetypical collision time are considered (
MSD ∝ Dt , where D ∗ [email protected] † [email protected] ‡ [email protected] is the diffusion constant, D = K B T / µ , with K B , T and µ being respectively the Boltzmann constant, the temperature,and a constant accounting for the geometry of the particleand the cinematic viscosity of the environment). One of themost familiar phenomena is indeed a sub-linear, power-lawincrease of the mean squared displacement, MSD ∝ D α t α ,with 0 < α < D α representing a diffusion constant withthe physical dimensions of [ L / t α ] . Such subdiffusive dynam-ics, also referred to as anomalous dynamics, is commonly at-tributed to the densely packed and heterogeneous structuresresulting from the crowding of the biological membranes. No-tably, atomistic and coarse-grained simulations, as well as thisstudy, have shown subdiffusive behavior also for lipids andproteins in simple lipid bilayers, in the limit of infinite proteindilution [16–18]. Moreover, as observed in many viscoelas-tic systems, subdiffusivity is a transient dynamical behavior[19]. After the short-time ballistic regime, where the taggedparticle freely diffuses ( MSD ∝ ( K B T / M ) t , with M beingthe the particle’s mass), the interactions with the mediumand its specific fluidic/mechanical properties can lead to non-trivial persistent correlations responsible for the subdiffusivebehavior [20]. Yet, for time-lags exceeding some characteris-tic time, the standard Brownian dynamics is recovered. Thecrossover from the subdiffusive to the standard Brownian dy-namics can take place over a quite large time window and thetransition onset strongly depends on packing and crowding, a r X i v : . [ phy s i c s . b i o - ph ] N ov ranging from tenth to hundreds of nanoseconds for lipids inprotein-free membranes or proteins in membranes at infiniteprotein dilution, up to arbitrarily long time scales for crowdedreal systems [17, 19, 21]. As a consequence, diffusive proper-ties can no longer be characterized by a single diffusion con-stant nor by a single exponent α . Several theoretical frame-works have been introduced to reproduce the subdiffusive be-havior and describe the physical mechanism behind it (seefor instance Ref. [19] and references therein for a completeoverview on their theoretical foundations and applications tobiological systems). Among them, we recall the so-calledGaussian models, such as the Fractional Brownian motion and the
Generalized Langevin Equation (GLE) [19], wherethe statistics of the noise in the relevant stochastic equationsis still assumed to be Gaussian like in the Langevin approachfor normal Brownian diffusion, but, differently from this, thenoise displays persistent power-law correlations (power-lawGaussian noise) leading to a sub-linear increase of the MSD.In the Generalized Langevin Equation, due to the fluctuation-dissipation theorem, this amounts to generalize the Stokesdrag term to a convolution of the velocity of the particle witha power-law memory kernel representing a generalized time-dependent friction [22]. Another important class is the one ofthe
Continuous-time random walks [19], where particles un-dergo a series of displacements in an homogeneous mediumaccording with a waiting-time distribution. Anomalous trans-port is connected with a power-law tail in the waiting-time dis-tribution, such that even the mean-waiting time is infinite. Thecentral-limit theorem does not apply in this case, since longerand longer waiting times are sampled. Finally, another impor-tant category is represented by the
Lorentz models [19] thatconsider spatially disordered environments where the particleexplores fractal-like structures, leading to anomalous dynam-ics. Recently, in the framework of the generalized Langevinequation model, hard exponential and soft power-law trunca-tion (tempering) of power-law memory kernels [23], as well asexponentially truncated three parameter Mittag-Leffler mem-ory kernels [24] have been used in order to quantitativelyreproduce the transition from subdiffusive to normal diffu-sion. Both exponential and soft power-law truncation implythe introduction of a characteristic cross-over time, related to amaximum-correlation time in the driving noise. Applicationsof these models to the analysis of the MSD of lipids in sim-ple lipid bilayers have shown that the cross-over time is hererelated to the time scale of mutual exchange between lipids[23, 25]. Compared to simply combining an anomalous anda normal diffusion law for the MSD, such phenomenologicalmodels naturally yield the emergence of subdiffusive to nor-mal crossover dynamics into a unique model, where the typeof truncation governs the crossover shape. The open challengeis to understand what are the relevant degrees of freedom andenvironment’s properties responsible for such phenomenolog-ical kernels.In our approach, we rely on a continuum description wherewe model the protein as a rigid body obeying Euler’s lawof motion, embedded in a linear non Newtonian viscoelas-tic fluid [19, 26–29] representing the lipid membrane. In or-der to account for the stochastic nature of the diffusive mo- tion of the protein, a noise term is added to the Euler’s law,which we call modified Euler’s law in the following. Themembrane is described by the Navier-Stokes hydrodynamicequations and constitutive equations, whose solutions deter-mine, respectively, the velocity vector field associated to thelipid membrane and the temporal relaxation to the equilibriumstate as a consequence of a deformation. Specifically, in theconstitutive equation, we model the viscous (instantaneous)and elastic (non instantaneous) behaviors of the membranethrough, respectively, a time-dependent Dirac delta functionand the Prabhakar fractional derivative [30, 31]. A solution ofsuch an equation can be provided in terms of a three-parameterMittag-Leffler function [32, 33] so that one obtains a non triv-ial time dependence in the response function of the mem-brane. By eliminating the lipid membrane degrees of freedom,one obtains a Generalized Langevin Equation (GLE) describ-ing the diffusive motion of the protein, whose memory ker-nel coincides with the solution of the constitutive equationabove mentioned. With this memory kernel, imposing spe-cific constraints on the three-parameter Mittag-Leffler func-tion, we derive a model for the MSD, the Velocity Auto-Correlation Function (VACF) and time-dependent diffusioncoefficient (TDC) that naturally accounts for the different dif-fusive dynamical regimes (ballistic, subdiffusive and Brow-nian) and the crossover between them, in terms of few phe-nomenological parameters that can be interpreted on a phys-ical base. The reliability of the model is tested versus theMSD data of a protein (the muscarinic M2 receptor) embed-ded in a mixed membrane, produced by MD simulations. Themodel further provides an estimation for the membrane vis-cosity, which is of the order of magnitude of the values foundby rheology experiments on mixed membrane [34].
OUTLINE
The article is organized as follows. In Section I, we firstpresent the geometry of the system (Subsection I A) and thenthe equations describing the protein (Subsection I B) and thelipid membrane (Subsection I C) dynamics. In Subsection I Dwe derive the Stokes’s equations for the membrane domain.The coupling conditions between the protein and membranedynamics are discussed in Subsection I E.In Section II, we derive the relation between the Cauchy’sstress tensor and the drag term appearing in the generalizedLangevin equation (Subsection II A), we provide a specificconstitutive equation for the lipid membrane using the Prab-hakar fractional derivative (Subsection II B) and we presentan explicit form of the GLE for the protein motion (Subsec-tion II C). The Mean Squared Displacement (MSD), the Ve-locity Auto-Correlation Function (VACF) and time-dependentdiffusion coefficient are discussed in Sub-subsection II C 1,and their asymptotic behavior in Subsection II D. In Sec-tion III, we compare the model to in-silico MSD data by MDsimulations.
I. EFFECTIVE CONTINUUM MODELA. Geometric Framework
In the following, we present the geometric assumptions ofthe model:• The system is embedded in a three dimensional Eu-clidean space, R , filled by two fluids (water and lipidmembrane) and by a rigid body (protein).• The lipid membrane, F , is modeled as a non Newto-nian viscoelastic fluid film with a finite thickness, h ,comparable with the protein’s height and infinitely ex-tended along the x - y directions.• The water domain, W = R \ ( F ∪ B ) , is modeled as aNewtonian fluid divided into two domains which com-pletely fill the physical space above and below the lipidmembrane.• The protein is modeled as a rigid cylinder, B ⊂ R ,with radius R and height equal to the thickness of thelipid membrane, h . We denote the mantle as ∂ B side and the two bases as ∂ B base . W
B. Continuum Description of the Brownian Particle: ModifiedEuler’s Laws.
Following Ref. [35], a rigid body is characterized by a do-main, B , in an oriented 3-dimensional Euclidean space andby a finite positive mass measure µ on B , such that the totalmass of B is: M : = µ ( B ) = (cid:90) B d µ . Being interested in the lateral diffusion, the only relevant de-gree of freedom is the rigid body’s center of mass velocity, U c ( t ) , which is described by the (modified) first Euler’s law[36, 37]: M d U c dt ( t ) = (cid:90) ∂ B T ( x , t ) n ( x ) d σ + Ξ ( t ) , (1)where ∂ B is the boundary of the rigid body and n the out-ward pointing normal on ∂ B ; T is the Cauchy’s stress ten-sor relative to the fluid, and Ξ is a noise term. We note thatwithout the noise term, Eq. (1) reduces to the Euler’s law,which describes the deterministic trajectory of the center ofmass of a rigid body in a fluid. Thanks to the integration overthe rigid body’s surface, we can introduce two mathematicalobjects written in two different descriptions within the sameequation: the space-time dependent field, T , in the Euleriandescription, and the time-dependent velocity vector, U c , inthe Lagrangian description. The noise term is introduced tomake the system stochastic. It is worth stressing that, at thislevel, such term is just an ansatz since we cannot formulate afluctuation-dissipation theorem. We will show later that thismodified Euler’s law, under suitable conditions, reduces to ageneralized Langevin equation. C. Continuum description of the lipid membrane and thewater domain
We adopt for both, the lipid membrane and the water do-main a continuum description. Although the size of a lipidmolecule is of the same order of magnitude of that of a pro-tein, which would make a priori questionable the continuumapproach, we are interested here in the mesoscale collectiveresponse of the lipid membrane to the protein diffusion as awhole.With this premise, the dynamical description of both, the lipidmembrane and water, is thus formulated in terms of the fluiddynamics equations [38, 39], namely, the continuity equation for the mass density, the conservation of momentum vector and the constitutive equation for the stress tensor defined inthe regions F (membrane) and W (water). Since the math-ematical derivation we are going to outline is valid for anyfluid, we introduce a generic fluid, E , and we will specializethe equations for the lipid membrane and water only at theend.The continuity equation involves the scalar density mass ρ : E × R −→ R + such that, given any portion P ⊂ E , the inte-gral m ( P ) = (cid:90) P ρ ( x , t ) d x , gives the mass of P . Imposing the conservation of mass, wededuce the continuity equation for the mass density function ∂ ρ∂ t ( x , t ) + ∇ · ( ρ ( x , t ) u ( x , t )) = , where u : P × R −→ R is the velocity vector field associatedto the portion of the considered fluid.The conservation of momentum is essentially Newton’s equa-tion for any portion of fluid. We associate to any portion P ⊂ E a linear momentum vector P ( P , t ) : = (cid:90) P u ( x , t ) ρ ( x , t ) d x . Each portion is subject to a body force , f b ( P , t ) , such as grav-ity, for example, and a contact force , f c ( P , t ) , due to thestress exerted by other fluid portions on the boundary of P .Therefore, we can combine them f ( P , t ) : = f b ( P , t ) + f c ( P , t ) , and write the conservation of the linear momentum d P dt ( P , t ) = f ( P , t ) . (2)Note that in our model f b ( P , t ) =
0, since the system is as-sumed to be rigid and gravity is neglected. Contact forces areinstead modeled using the Cauchy’s hypothesis, i.e.: f c ( P , t ) = (cid:90) ∂ P T ( x , t ) · n ( x ) d σ ( x ) , where n ( x ) is a pointing outward unit vector orthogonal tothe surface of P and T is the Cauchy’s stress tensor.In order to characterize the features of a fluid and thus fix thetensor, T , one has to introduce the constitutive equation . Thisis, in principle, a relation which could involve derivatives orintegrals of the strain tensor, D [ u ] , namely, the symmetric part of the gradient of the velocity vector field, i.e.: D [ u ] = (cid:0) ∇ u + ( ∇ u ) T (cid:1) . The explicit form of the constitutive equation depends on thespecific fluid. For a
Newtonian fluids , like water [40], the con-stitutive equation reads T W ( x , t ) = − p ( x , t ) R + µ W D [ u ( x , t )] , where the relation between stress and strain tensors is lin-ear. For a non-Newtonian fluid as we model the lipidmembrane[19, 26–29], the Cauchy’s stress tensor reads: T F ( x , t ) = − p ( x , t ) R + S [ D [ u ( x , t )]] , where the dependency on the strain tensor is in general non-linear. Specifically, for a viscoelastic fluid[41], it might belinear or non-linear but, in general, it depends on the past de-formations through the operator S , which is non-local in time.We model the constitutive equation for the lipid membrane asa linear viscoelastic fluid whose mathematical structure is T F ( x , t ) = − p ( x , t ) R + µ s D [ u ( x , t )]+ (cid:90) t G ( t − a ) D [ u ( x , a )] da , (3)where µ s is a coefficient with the physical dimensions of aviscosity representing the instantaneous response of the fluidto a deformation, whereas G is the so-called elastic-responsefunction . It strictly depends on the specific fluid we are con-sidering and it will be determined solving the constitutive dif-ferential equation that we will introduce later.We now summarize the equations governing the effective con-tinuum model, where (cid:63) stands for F or W (note that G W = ∂ ρ (cid:63) ∂ t ( x , t ) + u (cid:63) ( x , t ) · ∇ ρ (cid:63) ( x , t ) = ∇ · u (cid:63) ( x , t ) = ρ (cid:63) ( x , t ) (cid:18) ∂ u (cid:63) ∂ t ( x , t ) + u (cid:63) ( x , t ) · ∇ u (cid:63) ( x , t ) (cid:19) = ∇ · T (cid:63) ( x , t ) + f (cid:63) ( x , t ) f (cid:63) ( x , t ) Constitutive equation Deviatoric term T (cid:63) = − p (cid:63) R + µ s D (cid:63) [ u ] + S (cid:63) [ D [ u ]] S (cid:63) = (cid:90) t G (cid:63) ( t − a ) D (cid:63) [ u ( a )] da (4) D. Formal solution of the effective continuum fluid equations.
The goal here is to reduce the system of equations (4) intoa single partial differential equation for each domain. Wenote that the only difference between the water and membraneequations lies in the deviatoric term of the Cauchy’s stress ten-sor. Therefore we will work on the lipid membrane domain, knowing that the corresponding equation for the water domaincan be obtained by substituting G → µ s → µ W .The equation for the conservation of momentum in (4) is non-linear in u F . This makes harder to find the solution, whichin general may not exist. However, under the hypothesis oflow Reynolds number ( Re (cid:28) u F · ∇ u F [42] and time derivative of the velocity vector field ∂ u F / ∂ t can be disregarded. The resulting dynamic regimeand equation reduces to, respectively, the Stokes’ regime and
Stokes’ equation , and the system of equations (4) in F ⊂ R is given by ∂ ρ F ∂ t + u F · ∇ ρ F = , ∇ · u F = , − ∇ p F + µ s ∆ u F + ∇ · S F = , S F ( x , t ) = (cid:90) t G ( t − u ) D [ u F ( x , u )] du . (5)Once the linearity in u F is enforced, we can apply the Fouriertransform with respect to time and space to the scalar andvector fields above. Then, after some algebraic operation,one gets, by applying the inverse Fourier transform, a uniqueequation of the form − ∇ p F ( x , t ) + µ F ( t ) (cid:63) ∆ u F ( x , t ) = , x ∈ F , (6)which is similar to the Stokes’ equation for an incompressibleNewtonian fluid but where the viscosity is a time-dependentoperator µ F ( t − u ) : = µ s δ ( t − u ) + G ( t − u ) . (7)acting on ∆ u F as µ F ( t ) (cid:63) ∆ u F ( x , t ) = (cid:90) t µ F ( t − a ) ∆ u F ( x , a ) da . For the water domain, we change in Eq. (6) the subscript F with W and we set G →
0. The time-dependent viscosity µ W reduces to µ W ( t − u ) : = µ s δ ( t − u ) , (8)so that the convolution becomes the product of functions µ W ( t ) (cid:63) ∆ u W ( x , t ) = (cid:90) t µ W ( t − a ) ∆ u W ( x , a ) da = µ s ∆ u W ( x , t ) . Hence, the governing equation of water corresponds exactlyto the Stokes’ equation for a Newtonian fluid, i.e.: − ∇ p W ( x , t ) + µ s ∆ u W ( x , t ) = , x ∈ W . (9)Eqs. (6) and (9) describe separately the internal dynamics ofthe lipid membrane and the water, respectively.We also note that a similar derivation has already been used in[43] in the particular case of a G -function of the form˜ G ( ω ) : = µ p ( − i ωλ + ) , (10)starting from the Stokes-Oldroyd B model for the constitutiveequation in the Stokes’ regime. E. Coupling Conditions
In contrast to Eq. (9), Eq. (6) is non-local in time. There-fore, the Oseen’s tensor cannot be directly used to describeits homogeneous solution. However, by Fourier transformingwith respect to the time variable, the equation becomes local(in frequency ω ) and it can be formally solved in absence ofboundary conditions using the hydrodynamics Green’s func-tion, ˜ I F ( x , ω ) [44–46]. Indeed, the Fourier transformed Os-een’s tensor for a three dimensional fluid is given by (usingagain the notation (cid:63) = F , W )˜ I (cid:63) ( x , ω ) = π ˜ µ (cid:63) ( ω ) (cid:107) x (cid:107) (cid:18) R + x ⊗ x (cid:107) x (cid:107) (cid:19) , (11)where ˜ µ (cid:63) ( ω ) is, choosing the suitable subscript, the Fouriertransform of the time-dependent viscosity (7) or (8).Note that, from now on, we will work in the frequency domainand we will denote the Fourier transform of a function or fieldwith the overscript ˜.Let us introduce the following notation: Γ F : = ∂ F , Γ W : = ∂ W and Γ B = ∂ B base ∪ ∂ B side to be respectively the lipidmembrane, water and protein boundaries.For x ∈ Γ B ∩ Γ F , we choose the following interface condi-tions: ˜ u x F ( x , ω ) = ˜ U xc ( ω ) ˜ u y F ( x , ω ) = ˜ U yc ( ω ) , ˜ u z F ( ω ) = . (12)whereas, for x ∈ Γ B ∩ Γ W , we impose: ˜ u x W ( x , ω ) = , ˜ u y W ( x , ω ) = , ˜ u z W ( x , ω ) = ˜ U zc ( ω ) . (13)Note that the superscripts x , y , z denote the three componentsof the corresponding vector.Having specified the interface conditions, one can solveEq. (6), imposing the following ansatz for ˜ u F , ∀ x ∈ F :˜ u F ( x , ω ) = (cid:90) Γ B ∩ Γ F ˜ I F ( x − y , ω ) ˜ F F ( y , ω ) d σ ( y ) , (14)where ˜ F F : R × R −→ R is an unknown vector which hasto be determined inverting the relation between ˜ F F and ˜ u F .In order to describe the confinement of the protein within thelipid membrane, we also impose the following conditions forthe unknown forces:˜ F F ( y , ω ) = ( ˜ F x F ( y , ω ) , ˜ F y F ( y , ω ) , ) , ˜ F W ( y , ω ) = ( , , ˜ F z W ( y , ω )) . (15)Now, let us compute the solution far from the rigid body. Weconsider y ∈ Γ B ∩ Γ F fixed and we pick a point x ∈ F suchthat (cid:107) x − y (cid:107) (cid:29) (cid:107) y (cid:107) . In this way, we introduce a coordinatessystem (for the sake of simplicity, we adopt the cylindricalcoordinates) on the plane spanned by the x and y vectors.This allows us to employ the polar coordinates: x − y = ρ x cos θ x − ρ y cos θ y ρ x sin θ x − ρ y sin θ y z x − z y . In this coordinates, the underlying assumption translates to ρ x (cid:29) ρ y , so that we obtain: x − y = ρ x cos θ x − ρ y ρ x cos θ y sin θ x − ρ y ρ x sin θ yz x ρ x − z y ρ x ≈ ρ x cos θ x sin θ x , where we also assumed that z x and z y and ρ y are small com-pared to ρ x .Therefore, the tensor ˜ I F ( x − y , ω ) reduces to ˜ I F ( x , ω ) ,which corresponds to the first order Taylor approximation ofthe Oseen’s tensor in the frequency domain, and the ansatz(14) becomes: ˜ u F ( x , ω ) ≈ ˜ I F ( x , ω ) ˜ f F ( ω ) . (16)where ˜ f F ( ω ) : = (cid:90) Γ B ∩ Γ F ˜ F F ( y , ω ) d σ ( y ) Thereby, we define the leading contribution to the fluid veloc-ity vector field to be:˜ u L F ( x , ω ) = ˜ I F ( x , ω ) ˜ f F ( ω ) . (17)We can now invert the above relation by exploiting the inter-face conditions. We thus impose (12) in Eq. (14) and thus weintegrate over the protein surface so to obtain (cid:90) Γ B ∩ Γ F ˜ u L F ( x , ω ) d σ ( x ) = σ ( Γ B ∩ Γ F ) ˜ U c ( ω ) , where σ ( Γ B ∩ Γ F ) is the measure of the surface Γ B ∩ Γ F .We then obtain: ˜ U c ( ω ) = ˆ E ( R , h , ω ) ˜ f F ( ω ) , (18)where the constant (in space) tensor ˆ E is given by:ˆ E ( R , h , ω ) = σ ( Γ B ∩ Γ F ) − ˜ E ( R , h , ω ) , ˜ E ( R , h , ω ) = (cid:90) Γ B ∩ Γ F ˜ I F ( x , ω ) d σ ( x ) . (19) Therefore, the velocity of the rigid body’s center of mass justdepends on the constant (in space) force vector ˜ f F and ageometry-dependent matrix, ˜ E ( R , h , ω ) , so that the relation(18) can be easily inverted obtaining:˜ f F ( ω ) = ˆ E ( R , h , ω ) − ˜ U c ( ω ) . (20)In Appendix A, we develop explicit expressions for ˜ E andthus ˜ f F ( ω ) as a function of ˜ U c ( ω ) .Accepting the approximation introduced above, the solutionsfor the pressure scalar field and the velocity vector field read˜ u F ( x , ω ) = π ˜ Π ( x , ω ) (cid:107) x (cid:107) (cid:18) ˜ f F ( ω ) + ( x · ˜ f F ( ω )) (cid:107) x (cid:107) x (cid:19) , ˜ p F ( x , ω ) = x · ˜ f F ( ω ) π (cid:107) x (cid:107) . (21)The Fourier transform of the strain tensor, D [ ˜ u F ] , can thenexplicitly be given by: D [ ˜ u F ( x , ω )] = ( ∇ ˜ u F ( x , ω ) + ( ∇ ˜ u F ( x , ω )) T )= x · ˜ f F ( ω ) π ˜ Π F ( x , ω ) (cid:107) x (cid:107) (cid:18) R − x ⊗ x (cid:107) x (cid:107) (cid:19) . (22)Furthermore, in order to make Eq. (1) explicit, we have tocompute the Cauchy’s stress tensor: T F ( x , t ) = − p F ( x , t ) R + (cid:90) t Π F ( x , t − a ) D [ u F ( x , a )] da . (23)By applying the Fourier transform to the Cauchy’s stress ten-sor one gets˜ T F ( x , ω ) = − ˜ p F ( x , ω ) R + Π F ( x , ω ) D [ ˜ u F ( x , ω )] , (24)which yields, upon combining (21) with (24),˜ T F ( x , ω ) = − ( ˜ f F ( ω ) · x ) π (cid:18) x ⊗ x (cid:107) x (cid:107) (cid:19) . (25)The resulting rigid body dynamics turns out to be decoupledinto two contributions; the first lies on the membrane proteinplane and involves the U xc , U yc components of the body veloc-ity, whereas the second just involves the U zc component. Thisis a consequence of the interface conditions we have chosentogether with the fact that the tensor ˜ E does not combine the x , y components of the velocity with the z one (see Eq. (A9)).Therefore, being the first pair of velocity components decou-pled from the third component, we will just focus on the dy-namics throughout the lipid membrane plane neglecting thedynamics along the z -axis from now on. II. GENERALIZED LANGEVIN EQUATION FROM THEEFFECTIVE CONTINUUM MODEL
In this Section, we explain the procedure to eliminate thedegrees of freedom of the fluid (membrane) and derive a gen-eralized Langevin equation for the protein, which is modeledas a Brownian particle diffusing in a linear viscoelastic fluid.We also consider a realistic constitutive equation, which es-tablishes the response of the fluid to a deformation caused bythe protein motion.
A. Elimination of the degrees of freedom of the fluid.
Combining Eq. (20) and (25) we obtain an explicit map-ping between U c ( t ) → T F [ U c ( t )] so that we can now plugthis relation into the modified Euler’s law, Eq. (1), in order toobtain M d U c dt ( t ) = (cid:90) ∂ B T [ U c ( t )] n d σ + Ξ ( t ) , (26)where Ξ is the translational stochastic force and it is con-nected to the mechanical properties of the fluid. As a con-sequence, a fluctuation-dissipation relation can be formulatedjustifying a posteriori the introduction of the noise term. Byplugging the cylindrical coordinates x = ( R cos θ , R sin θ , z ) inEq. (25), the unit normal vector is n = ( cos θ , sin θ , ) and theintegral in Eq. (26) on the lateral surface of the cylinder canbe computed as (cid:90) h / − h / (cid:90) π ˜ T F (( θ , z ) , ω ) n ( θ ) R d θ dz = − R π (cid:90) h / − h / (cid:90) π ( x ( θ , z ) ⊗ x ( θ , z )) (cid:107) x ( θ , z ) (cid:107) ˜ f F ( ω ) R d θ dz = − ( h + R ) / ( h + hR ) ˜ f x F ( ω )( h + hR ) ˜ f y F ( ω ) h ˜ f z F ( ω ) = − π h (cid:0) h + R (cid:1) Λ (cid:16) h + Λ log (cid:16) h ( Λ + h ) R + (cid:17)(cid:17) ˜ µ F ( ω ) ˜ U xc ( ω ) ˜ µ F ( ω ) ˜ U yc ( ω ) . (27)where Λ = √ h + R . Note that, in order to pass from thesecond to the last equality, we have replaced the force com-ponents with the explicit relation given in Eq. (A11) togetherwith the interface conditions (15).Finally, by defining the geometrical coefficient: P F ( R , h ) : = π h (cid:0) h + R (cid:1) Λ (cid:16) h + Λ log (cid:16) h ( Λ + h ) R + (cid:17)(cid:17) , (28) and applying the inverse Fourier transform to the last equalityin Eq. (27), we get the final result: (cid:90) ∂ B side T F (( θ , z ) , t ) n ( θ ) R d θ dz = − P F ( R , h ) (cid:82) t µ F ( t − u ) U xc ( u ) du (cid:82) t µ F ( t − u ) U yc ( u ) du . (29)Equation (29) proves that the integral appearing in Eq. (26)is equivalent to the convolution of the time-dependent fluidviscosity with the center of mass velocity of the protein, mul-tiplied by a geometrical factor. In other words, it takes theform of the drag term appearing in the generalized Langevinequation.We have thus demonstrated that Eq. (26) can be written in ageneralized Langevin-like equation for the x , y protein veloc-ity degrees of freedom as follows M dU ic dt ( t ) = − (cid:90) t ζ F ( t − u ) U ic ( u ) du + Ξ i ( t ) , i = x , y , (30)where the friction term is given by ζ F ( t ) : = P F ( R , h ) µ F ( t ) , (31)The { Ξ i ( t ) } i = x , y are the components of a Gaussian dis-tributed correlated thermal noise such that, by means of thefluctuation-dissipation theorem which can now be formulated,verifies the relation: (cid:104) Ξ i ( t ) Ξ j ( u ) (cid:105) = K B T µ F ( t − u ) δ i j , i = x , y . (32) B. Effective response of the lipid membrane to thedeformation induced by the protein
In order to make the effective continuum model descriptionoperative, we should make an assumption about the form ofthe elastic-response G -function, which depends on the spe-cific fluid under consideration (see for instance [47] for a lipidmonolayer). In our approach, the Brownian particle can ac-tually be viewed as a probe. This probe plays the same roleas a mechanical deformation in rheology but which is now lo-calized in a small region of the fluid around the particle. Thisis the reason why one is legitimated to assume the Stokes’regime (see Section I D). The response of the fluid is encodedin the friction term of the generalized Langevin equation.The simplest model one can consider is the Jeffrey’s fluididentified by the Stokes-Oldroyd B constitutive equations inthe Stokes’ regime. However, a time-dependent exponentialresponse is not sufficient to reproduce the different dynamicalregimes featured by the MSD of a protein laterally diffusingin a membrane. Specifically, the first regime is the ballistic one, occurring at short time-lags (sub-ps), where the proteindoes not interact with the fluid yet and the MSD just dependson the protein mass and temperature. Then, after a character-istic time corresponding to the onset of the interactions withthe next neighbor lipid molecules, a transition towards a sub-diffusive regime is observed. During the transition, the fluidsurrounding the protein reacts to the deformations imposed bythe protein to reach a new equilibrium state (relaxation phase)in a non trivial way. Several lipids molecules have to rearrangeto allow the protein to move. The power law behavior suggestsa distribution of relaxation times with a long time-tail. Finally,for time-lags larger than some characteristic correlation time,the standard Brownian dynamics regime is restored.In order to correctly reproduce the response function of thelipid membrane, we model the membrane as a linear vis-coelastic material. Within this hypothesis, we adopt the
Maxwell–Prabhakar fractional constitutive differential equa-tion [48, 49] S ( t ) = P λ , ν , δ C D − δλ , − ν , A λ D [ u ( t )] , (33)where C D αβ , γ , A λ is the Prabhakar fractional timederivative[30, 31], λ , ν , δ are dimensionless parame-ters, and A λ is a parameter, depending on λ , related to acharacteristic time. Finally, P λ , ν , δ is a parameter, dependingon λ , ν and δ , which is introduced to get the right physicaldimensions for S ; it plays the role of an elastic parameter.
1. A Realistic Constitutive Equation
The Prabhakar derivative in Eq. (33) is defined by D − δλ , − ν , A λ f ( t ) : = E δλ , m + ν , A λ f ( m ) ( t ) , (34)where the index m is an integer representing the order of thestandard time-derivative and f ∈ AC m ( , t ) , i.e. the spaceof absolutely continuous functions[30]. The non-local oper-ator on the right-hand side denotes the Prabhakar fractionalintegral[32], defined as E δλ , ν , A λ f ( t ) : = (cid:90) t ( t − a ) ν − E δλ , ν (cid:16) A λ ( t − a ) λ (cid:17) f ( a ) da . (35)Here E δλ , ν is the 3-parameter Mittag-Leffler function, alsocalled Prabhakar function[32, 33] E δλ , ν ( z ) = ∞ ∑ k = ( δ ) k Γ ( λ k + ν ) z k k ! , (36) where ( δ ) k is the Pochhammer symbol and λ , ν , δ ∈ C with Re [ λ ] > Prabhakar kernele δλ , ν ( A λ , t ) : = t ν − E δλ , ν (cid:16) A λ t λ (cid:17) , (37)one can rewrite the constitutive equation (33) as S ( t ) = P λ , ν , δ (cid:90) t e δλ , ν ( A λ , t − a ) D [ u ( a )] da . (38)Note that the constitutive equation written in this form canprovide a solution for the S tensor as long as we know the timebehavior of the strain tensor D [ u ( t )] and the initial condition, D [ u ( t = )] = D [ u ] .By substituting Eq. (38) into Eq. (3), we get: T ( x , t ) = − p R + µ s D [ u ]+ P λ , ν , δ (cid:90) t e δλ , ν ( A λ , t − a ) D [ u ( a )] da , (39)thus A λ has the dimension of an inverse fractional time [t − λ ],in particular: A λ : = − τ − λ , (40)with τ being a characteristic time, so that the 3-parameterMittag-Leffler function in the Prabhakar kernel is dimension-less, and the Prabhakar kernel reads e δλ , ν ( − τ − λ , t ) : = t ν − E δλ , ν (cid:20) − (cid:16) t τ (cid:17) λ (cid:21) . (41)The latter has the dimension of a fractional power of the time( [ t ] ν − ). Therefore, in order to have the correct physical di-mension for the Cauchy stress tensor in Eq. (39) and to havea consistent notation, one has to impose P λ , ν , δ : = µ p τ ν , (42)where µ p has got the same dimension of µ s .The Cauchy stress tensor (39) thus takes the form T ( x , t ) = − p R + (cid:90) t (cid:32) µ s δ ( t − a ) + µ p ( t − a ) ν − τ ν E δλ , ν (cid:34) − (cid:18) t − a τ (cid:19) λ (cid:35)(cid:33) D [ u ( a )] da . (43) C. The generalized Langevin Equation for Protein Diffusion
By exploiting (29), which relates the Cauchy’s stress tensorto the drag term of the generalized Langevin equation, and using the result of equation (43), we can now write equation(30) as:
M d U F c dt ( t ) = − (cid:90) t Σ F ( t − a ) U F c ( a ) da + Ξ F ( t ) , (44)where the kernel tensor Σ F is given by Σ F ( t − a ) : = ζ F ( t − a ) R , ζ F ( t − a ) : = ξ s δ ( t − a ) + ξ p τ ν ( t − a ) ν − E δλ , ν (cid:34) − (cid:18) t − a τ (cid:19) λ (cid:35) , ξ s : = P F ( R , h ) µ s , ξ p : = P F ( R , h ) µ p . (45)The random force Ξ F ( t ) is a zero-mean ( (cid:104) Ξ ( t ) (cid:105) =
0) noise,related to the friction term through the fluctuation-dissipationtheorem (cid:104) Ξ F ( t ) ⊗ Ξ F ( a ) (cid:105) = K B T Σ F ( t − a ) , (cid:104) Ξ F ( t ) · Ξ F ( a ) (cid:105) = Tr [ (cid:104) Ξ F ( t ) ⊗ Ξ F ( a ) (cid:105) ]= K B T Σ F ( t − a ) . (46)One can also introduce the memory matrix, K F , defined as K F = Σ F / M , which reads K F ( t − a ) = κ F ( t − a ) R , κ F ( t − a ) = ξ s M δ ( t − a ) + ξ s M ( t − a ) ν − τ ν E δλ , ν (cid:34) − (cid:18) t − a τ (cid:19) λ (cid:35) . (47)We note that Σ F (and K F ) are diagonal matrices. This is aconsequence of the Stokes’ regime hypothesis, that is, no hy-drodynamic interactions coupling different velocity compo-nents are included, in combination with the geometrical set-ting and the first order approximation (16).Equation (44) thus reduces to a system of two decoupledstochastic differential equations, accounting for the lateral dif-fusion in the x , y plane of the membrane. This descriptionholds as long as the considered lipid membrane presents anhomogeneous distribution of lipids, regardless of whether itschemical composition is heterogeneous. We note that setting µ p =
0, Eq. (44) reduces to the standard Langevin equation,so that the model is self-consistent.
1. Solution of the Model
In order to solve equation (44), we follow the approach usedin Ref. [33]. By applying the Laplace transform, Eq. (44)reduces to: ˜ U c ( s ) = M U c ( ) Ms + ˜ ζ F ( s ) + Ξ ( s ) Ms + ˜ ζ F ( s ) , (48)where the Laplace transform of ζ F is[33]:˜ ζ F ( s ) : = L { ξ s δ ( t − u ) } + L (cid:40) ξ p τ ν ( t − u ) ν − E δλ , ν (cid:34) − (cid:18) t − u τ (cid:19) λ (cid:35)(cid:41) = ξ s + ξ p ( τ s ) δλ − ν ( + ( τ s ) λ ) δ . (49) By introducing the so-called ˜ g ( s ) - relaxation function , g -RF:˜ g ( s ) : = Ms + ˜ ζ ( s ) = M s + ξ s + ξ p ( τ s ) δλ − ν ( +( τ s ) λ ) δ , (50)equation (48) reads:˜ U c ( s ) = M U c ( ) ˜ g ( s ) + Ξ ( s ) ˜ g ( s ) . (51)By applying the inverse Laplace transform, one obtains theformal solution for the velocity vector U c ( t ) = (cid:104) U c ( t ) (cid:105) + (cid:90) t g ( t − u ) Ξ ( u ) du , (cid:104) U c ( t ) (cid:105) = M U c ( ) g ( t ) . (52)As for the position vector, by substituting U c = ˙ X c in (44)and applying the Laplace transform, one gets:˜ X c ( s ) = X c ( ) s + M U c ( ) Ms + s ˜ ζ F ( s ) + Ξ ( s ) Ms + s ˜ ζ F ( s ) . (53)By introducing the ˜ H ( s ) - relaxation function , H -RF:˜ H ( s ) : = Ms + s ˜ ζ F ( s ) = M s + s (cid:16) ξ s + ξ p ( τ s ) δλ − ν ( +( τ s ) λ ) δ (cid:17) , (54)equation (53) can be recast in the form:˜ X c ( s ) = X c ( ) s + M U c ( ) ˜ H ( s ) + Ξ ( s ) ˜ H ( s ) . (55)Note that ˜ H is related to ˜ g through the relation˜ H ( s ) = ˜ g ( s ) s . (56)Once again, by performing the inverse Laplace transform ofequation (55), one gets X c ( t ) = (cid:104) X c ( t ) (cid:105) + (cid:90) t H ( t − u ) Ξ ( u ) du , (cid:104) X c ( t ) (cid:105) = X c ( ) + M U c ( ) H ( t ) . (57)We finally introduce the so-called I - relaxation function , I -RF: I ( t ) : = (cid:90) t H ( u ) du (58)which is related to the other relaxation functions through:˜ I ( s ) = ˜ H ( s ) s = M s + s (cid:16) ξ s + ξ p ( τ s ) δλ − ν ( +( τ s ) λ ) δ (cid:17) . (59)Knowing g , H and I , the stochastic process can be fully char-acterized by the statistical observables MSD, VACF and time-dependent diffusion coefficient, D ( t ) , through the relations(see Appendix B for the derivation): C v ( t ) = K B T g ( t ) D ( t ) = K B T H ( t ) (cid:104) ∆ X ( t ) (cid:105) = K B T I ( t ) . (60)0 D. Asymptotic behaviors of Statistical Observables.
1. Short time behavior
Let us rewrite the kernel function of the GLE ζ F ( t ) = ξ s δ ( t ) + ξ p τ ν t ν − E δλ , ν (cid:20) − (cid:16) t τ (cid:17) λ (cid:21) . (61)The short time behavior of the kernel function can be obtainedby taking the first terms in the series representation of theMittag-Leffler function E ρα , β ( − z ) ≈ Γ ( β ) − ρ z Γ ( α + β ) , ( | z | (cid:28) ) , (62)obtaining the following relation ζ F ( t ) ≈ ξ s δ ( t ) + ξ p τ ν t ν − (cid:32) Γ ( ν ) − δ t λ τ − λ Γ ( λ + ν ) (cid:33) = ξ s δ ( t ) + ξ p τ ν (cid:32) t ν − Γ ( ν ) − δ t λ + ν − τ − λ Γ ( λ + ν ) (cid:33) ≈ ξ s δ ( t ) + ξ p τ ν t ν − Γ ( ν ) , ( t (cid:28) τ ) . (63)We thus define the short time kernel function ζ S to be: ζ S ( t ) = ξ s δ ( t ) + ξ p τ ν t ν − Γ ( ν ) , (64)whose Laplace transform reads˜ ζ S ( s ) = ξ s + ξ p τ ν s − ν . (65)By substituting the asymptotic kernel function in the relax-ation functions, we define the short time asymptotic relaxationfunctions: ˜ g S ( s ) : = Ms + ξ s + ξ p τ ν s − ν , ˜ H S ( s ) : = Ms + ξ s s + ξ p τ ν s − ν + , ˜ I S ( s ) : = Ms + ξ s s + ξ p τ ν s − ν + . (66)Their inverse Laplace transform yields the analytic behaviorof the RFs at short time (see Appendix B 1 for the derivation) g S ( t ) ≈ M (cid:34) − ω s t (cid:35) , H S ( t ) ≈ M (cid:34) t − ω s t (cid:35) , I S ( t ) ≈ M (cid:34) t − ω s t (cid:35) , (67) where ω s : = ξ s M , (68)is a characteristic frequency setting the time scale of the tran-sition from the ballistic to subdiffusive regime.
2. Long time behavior
Here we focus on the relaxation function ˜ g ( s ) , since theasymptotic behavior of all the other RFs can be obtainedthrough the relations (54) and (58). By replacing Eq. (49)in Eq. (50), one gets˜ g ( s ) = M s + ξ s + ξ p ( τ s ) δλ − ν ( +( τ s ) λ ) δ . In the long time limit, one can assume M ξ s s (cid:28) . This legitimates the following approximation in the denomi-nator of ˜ g ( s ) M (cid:16) s + ω s + ω p ( τ s ) δλ − ν ( + ( τ s ) λ ) δ (cid:17) = M ω s (cid:32) s ω s + + ξ p ξ s ( τ s ) δλ − ν ( + ( τ s ) λ ) δ (cid:33) ≈ ξ s (cid:32) + ξ p ξ s ( τ s ) δλ − ν ( + ( τ s ) λ ) δ (cid:33) , (69)where we have introduced the characteristic frequency ω p ω p : = ξ p M . (70)By introducing the dimensionless parameter φ = ξ p / ξ s , wedefine the long term asymptotic relaxation functions by˜ g B ( s ) : = ξ s + φτ ν s δλ − ν ( τ − λ + s λ ) δ , ˜ H B ( s ) : = ξ s s + φτ ν s δλ − ν + ( τ − λ + s λ ) δ , ˜ I B ( s ) : = ξ s s + φτ ν s δλ − ν + ( τ − λ + s λ ) δ . (71)1 Parameters Physical Role R hydrodynamic radius µ s is the viscosity component representing the instantaneous response of the lipid membrane µ p is the elastic component representing the retarded (elastic) response of the lipid membrane ξ s friction component coming from the instantaneous response ξ s / M characteristic frequency setting the time scale of the transition from ballistic to subdiffusive regime ξ p friction component coming from the retarded response. It is the leading term in the asymptotic diffusion coefficient ξ p + ξ s total macroscopic friction felt by the protein τ time scale of the retarded (elastic) response of the lipid membrane δ , λ , ν δ λ = ν to asymptotically get the Brownian regime. They shape the transition from subdiffusive to Brownian regime TABLE I: Summary of the physical role of the parameters.The inverse Laplace transforms of the above relaxation func-tions yield (see Appendix B 2 for the derivation) g B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν − Γ ( δ λ − ν ) − τ λ φ ξ s t − λ − Γ ( − λ + ( δ λ − ν )) , H B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν Γ (( δ λ − ν ) + ) , I B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν + Γ (( δ λ − ν ) + ) . (72)Let us first analyze the H -RF which is related to the time-dependent diffusion coefficient through the relation D ( t ) = K B T H ( t ) .Since the Brownian regime should be recovered in the longtime limit, D ( t ) has to approach a constant D ∞ (cid:54) = δ λ − ν = , (73)so that D ∞ ≈ K B T H B ( t ) = K B T ξ p . (74)Moreover, this derivation highlights that the leading contribu-tion to the asymptotic diffusion coefficient is provided by thecoefficient ξ p .Nevertheless, the exact mathematical expression for D ∞ canbe obtained considering the non approximated expression for˜ H ( s ) that, under the constraint (73), reads˜ D ( s ) = K B Ts M s + ξ s + ξ p ( +( τ s ) λ ) δ . By taking the limit s →
0, we get1
M s + ξ s + ξ p ( +( τ s ) λ ) δ ≈ ξ s + ξ p , which finally produces˜ D ( s ) ≈ K B Ts ξ s + ξ p . The inverse Laplace transform gives then the exact expressionfor the asymptotic diffusion coefficient D ∞ = K B T ξ s + ξ p . (75)To recast the previous result, we now observe that if ξ s (cid:28) ξ p ,then K B T ξ p + ξ s = K B T ξ p + ξ s ξ p ≈ K B T ξ p , which coincides with the expression in Eq. (74).We conclude the section by discussing the asymptotic behav-ior of g -RF and I -RF which are related to the VACF and MSD,respectively.The first term of the g -RF expansion vanishes since wehave Γ ( ) at the denominator which, by definition, is zero ( / Γ ( ) = ) . The second term is well-defined provided that λ / ∈ Z and it produces: g B ( t ) ≈ − τ λ ξ p t − λ − Γ ( − λ ) , (76)where we used φ ξ s = ξ p .Finally, we observe that the MSD, asymptotically, linearly de-pends on time: (cid:104) ∆ X ( t ) (cid:105) = K B T ξ p t . (77)so to recover the Brownian regime.We summarize the physical role of each parameter in Table I. III. SIMULATION
As a proof of concept, we test the developed effectivemodel against in-silico MSD data of the center of mass of a protein laterally diffusing in a fully hydrated membrane. The2in-silico MSD data are extracted from quasi-atomistic Mar-tini [50–52] MD simulations. The model is thus fitted tothe numerical data in order to estimate the best values of themodel parameters.
A. System and MD simulations setup
We used as model system the M2 muscarinic acetylcholinereceptor (M2 receptor), an inhibitory class A G-protein-coupled receptor expressed both in the central and parasym-pathetic nervous systems [53]. The inactive form of M2 [53](pdb 3UON) was embedded in a mixed lipid bilayer contain-ing some of the most abundant species in neuronal cell mem-branes [54]. The antagonist and the fusioned T4 lysozyme,which in the original pdb structure replaces the third intra-cellular loop of the receptor, was removed from the system.About 14000 water molecules were added as well as Na + andCl − ions in order to neutralize the system and reach a physio-logical concentration of 0.15 M. A snapshot of the simulatedsystem is shown in Fig. C.1. The membrane compositionand the characteristic size of the system are reported in Ap-pendix C.The whole system was prepared by using CHARMM-GUImembrane builder web server [55]. The simulations were per-formed through the GROMACS program suite [56–59], ver-sion 2016.3, in double precision, using a Martini force field[50–52] with the leap-frog integration algorithm [60] and anintegration time step of 20 fs. Electrostatic interactions werecalculated using the reaction-field scheme and a cutoff of 1 . τ = τ =
12 ps, T =
310 K and reference pressure of P = analyze with the flag -msd . The maximum time lag consid-ered corresponds to 10% of the trajectory length. B. Comparison of the Effective Continuum model with thenumerical data
The generalized Langevin Equation based on the effectivecontinuum model introduced in Section II B 1 has been testedagainst the lateral MSD data extracted from the MD simula-tions. Two different flavors of the model have been initiallyimplemented, corresponding to different choices on the freeparameters θ :(M1) θ = ( ω s , ω p , τ , λ , δ ) , with ν = δ λ ;(M2) θ = ( ω s , ω p , τ , λ ) , with δ =
1. Together with the con-dition δ λ = ν , this leads to set ν = λ . Model (M1) corresponds to a three parameter Mittag-Lefflerfunction, it is the more generic model, whilst model (M2) re-duces to a two parameter Mittag-Leffler function, which hasbeen already used in the literature in the context of viscoelas-ticity [49]. The mass of the protein and the temperature arefixed to M = T =
310 K in both implemen-tations. The expected values of the parameters θ and their un-certainties have been evaluated by using a Bayesian approachwith a flat (uninformative) prior distribution for the param-eters [63], i.e. all the parameters values are assumed to beequally probable before analyzing the data (see Appendix Dfor details about the data analysis). Due to the complexityof the MSD model in the time domain, Eq. (60), the modelfor a given set of parameters is first calculated in the Laplacespace. Combining Eqs. (59), (54) and (49), under the condi-tion δ λ = ν , it reads (cid:104) (cid:103) ∆ X (cid:105) ( s ) = K B TM s s + ω s + ω p ( +( τ s ) λ ) δ . (78)Then, in order to carry out the Bayesian analysis for the bestfit parameters estimation, the latter is inverse Laplace trans-formed using the numerical Mathematica routine by Horváthet al. [64–66] for each set of sampled parameters (see Ap-pendix D for details).Both models provide a reasonable fit in terms of reduced ˜ χ .In order to compare them in a fair way and avoid overfitting(the number of free parameters is different in the two mod-els), we used the Bayesian information criterion [67] (BIC),which provides a posterior estimation of the probability forthe two models, assuming that the models under considera-tion are equally plausible a priori, and taking into considera-tion the number of free parameters. The posterior probabilityfor model (M1) and (M2) were respectively 0.30 and 0.70,thus indicating a week evidence for model (M2). The modelsare not really distinguishable. We thus report here only thedata and the analysis relevant to model (M2). The parametersand the 95% confidence intervals are summarized in Table IItogether with the parameters ξ s and ξ p derived thereof andthe Brownian diffusion coefficient D ∞ . See Appendix D fordetails about the parameters and confidence intervals estima-tion.The comparison in log-log scale between the model and thenumerical MSD data is shown in Fig. (2). In order to high-light the physical role of the model parameters, we report inthe same figure the asymptotic MSD obtained from the model (cid:104) ∆ r ( t ) (cid:105) t → ∞ = D ∞ t , with D ∞ given by Eq. (75), the tran-sition time scale τ from subdiffusive to Brownian diffusion,and the transition time scale 1 / ω s from the ballistic to subdif-fusive regime.The evolution of the time dependent diffusion coefficient asobtained from Eq. (60) upon normalization to D ∞ is shown inFigure 3. It has been obtained computing the inverse Laplacetransform of the H -RF provided by (54).Both Figs. 2 and 3 show that the subdiffusive to Browniandynamics transition lasts from tens to hundreds of nanosec-onds and that the Brownian dynamics is fully recovered onlyat time lags of the order of 500 ns. Within this time-frame,the protein moves about ≈ Parameters 3UON § δ ∗ λ = ν . . . ∗ τ (ns) 13 . . . ∗ ω s (ps − ) . . . ∗ ω p (ps − ) † ξ s (Pa · s · µ m) 0 . . . † ξ p (Pa · s · µ m) 0 . . . † D ∞ ( µ m · s − ) 0 . . . TABLE II: Best fit parameters ( ∗ ), fixed parameters ( § ), andparameters derived from the best fit parameters ( † ). The sub-and superscripts represent the 95% confidence intervals.radius of an average lipid, which is about ≈ < Δ � � > �� < Δ � � > ��� ω - � � τ < Δ � � > � →∞ ���� � ��� ����� �� × �� � �� × �� - � ������������� � ( �� ) < Δ � � ( � ) > ( � � � ) FIG. 2: Comparison between the MSD data from MDsimulations (red) and the Effective Continuum Model (blue).The long-time asymptotic MSD (black), and thecharacteristic times τ (green) and ω − s (orange) are shownadditionally. ��� ��� �� � ������������ ��� ��� �� � ������������� ( �� ) � ( � ) � ( ∞ ) FIG. 3: Time-dependent diffusion coefficient normalized tothe value D ∞ = K B T / ( ξ s + ξ p ) .timation of the viscosity components µ s and µ p , arising fromthe instantaneous and retarded response respectively. Plug-ging ξ s , ξ p (see Table II) and the radius and height of the pro-tein into Eqs. (45) and (28) (that we recall here for simplicity) ξ s : = P F ( R , h ) µ s , ξ p : = P F ( R , h ) µ p , P F ( R , h ) : = π h (cid:0) h + R (cid:1) Λ (cid:16) h + Λ log (cid:16) h ( Λ + h ) R + (cid:17)(cid:17) , one obtains µ s = . . . Pa · s , µ p = . . . Pa · s , with R = . h = . P F ( R , h ) = . R appearing in the geomet-ric factor P F is, in principle, the hydrodynamics radius, weused here the geometric radius. Both R and h have been de-termined through the tools furnished in the Protein Data Bankweb-site [68]. The order of magnitude of the total viscositycalculated through the model is comparable with those foundby rheology experiments for mixed membranes [34]. IV. CONCLUSIONS
The derived continuum effective model based on theviscoelastic assumption of the lipid membrane proved to beadequate for describing all the dynamical regimes appearingin a diffusive process. Specifically, a viscoelastic responsemade of an instantaneous component δ ( t ) , and a retardedcomponent, modeled by t δλ − E δλ , δλ ( − ( t / τ ) λ ) , naturallydescribes the emergence of the transition from the ballisticto the subdiffusive regime and from the subdiffusive to theBrownian regime. The retarded component underlies theexistence of a distribution of relaxation times to the stationarystate, lasting till ≈
500 ns, which corresponds to the time scalewhen the Brownian dynamics is fully recovered. Within thistime frame, the protein undergoes a displacement of the orderof 0.6 nm, which is comparable to the radius of an averagelipid molecule, 0.4 nm. The subdiffusive behavior can thus be4ascribed to the local interactions between the protein and thefirst lipid shell, leading to a “transient trapping” effect. Theparameters τ , δ , and λ shape the underlying distribution ofrelaxation times, while ω − s indicates the characteristic timeof the onset of the interactions between the protein and theneighbour lipids. The model further provides the viscosity ofthe membrane (arising from the instantaneous and retardedresponse), µ p + µ s ≈ .
787 Pa · s, which is of the same orderof magnitude of the reported experimental data on mixedmembranes.[34]As a next step we plan to study the underlying distribution ofrelaxation times as a function of the membrane compositionand protein size. The long term perspective is to derive amodel able to fully predict the different dynamical regimesof a protein in terms of few physical parameters such as theviscosity of the medium, the mass and the geometry of theprotein. ACKNOWLEDGMENTS
The authors would like to offer their special thanks to TrifceSandev for his advises about the theoretical analysis of ourstochastic equation. The authors also would like to thank Ric-cardo Capelli and Luca Pesce for fruitful discussions and helpin MD simulations setting.This work was supported by the Jülich-Aachen Research Al-liance Center for Simulation and Data Science (JARA-CSD)School for Simulation and Data Science (SSD).
Appendix A: Computation of the long-range part of the Oseen’stensor.
We are going to compute the integral (19) introduced inSection I E that we recall below˜ E ( R , h , ω ) = (cid:90) B side ˜ I F ( x , ω ) d σ ( x ) , (A1)being Γ B ∩ Γ F = ∂ B side . We will develop an explicit expres-sion of ˜ E in (A1) depending on the geometrical parameters R , h of the rigid body.We remind here that the protein is modeled as cylinder of ra-dius R and height h , centered around the origin and alignedwith the z -axis. The cylinder’s bases are in contact with thewater domain located above and below the lipid membrane,while the mantle is in contact with the lipid membrane. Wethus write ∂ B side = S R × [ − h / , h / ] , where S R denotes the circle of radius R .We then introduce the following parametrization: ∂ B side : = { s sR ( θ , z ) : = ( R cos θ , R sin θ , z ) | ( θ , z ) ∈ [ , π ) × [ − h / , h / ] } , (A2)The gradients g sz , g s θ of the parametrization s sR with respect tothe parameters z , θ then read g s θ : = − R sin θ R cos θ , g sz : = . (A3)The normal vector to the surface, oriented as outward point-ing, is given by: g s ⊥ = g s θ × g sz , with (cid:107) g s ⊥ (cid:107) = R . (A4)and the differential element on the mantle reads: d σ ( x ) = R d θ dz , x = s sR ( θ , z ) ∈ ∂ B side . (A5)Therefore, the integral (A1) becomes:˜ E ( R , h , ω ) = R (cid:90) h / − h / (cid:90) π ˜ I F ( s sR ( θ , z ) , ω ) d θ dz (A6)Now, we consider the specific form of the Oseen’s tensorgiven in (11). In the new coordinates (A2), Eq. (A6) canbe split into two terms:˜ I F ( s sR ( θ , z ) , ω ) = ˜ I F ( z , θ , ω ) + ˜ I F ( z , θ , ω ) , ˜ I F ( z , θ , ω ) = π ˜ µ F ( ω ) √ z + R , ˜ I F ( z , θ , ω ) = R (cid:0) z + R (cid:1) − / π ˜ µ F ( ω ) R cos ( θ ) R sin ( θ ) cos ( θ ) z cos ( θ ) R sin ( θ ) cos ( θ ) R sin ( θ ) z sin ( θ ) z cos ( θ ) z sin ( θ ) z . (A7)The splitting above automatically divides the tensor ˜ E into two terms. Indeed, the integration of both ˜ I F and ˜ I F pro-duces:5˜ E ( R , ω , h ) = R µ F ( ω ) log h (cid:16) √ h + R + h (cid:17) R + , ˜ E ( R , ω , h ) = R µ F ( ω ) h √ h + R h √ h + R
00 0 2 csch − (cid:0) Rh (cid:1) − h √ h + R . Introducing the auxiliary terms H x , y F : = π h (cid:16) h Λ + log (cid:16) h ( Λ + h ) R + (cid:17)(cid:17) H z F : = π h (cid:16) − h Λ + log (cid:16) h ( Λ + h ) R + (cid:17) + − (cid:0) Rh (cid:1)(cid:17) (A8)with Λ = √ h + R , we obtain the diagonal tensor ˜ E as:˜ E ( R , ω , h ) = π RhH x , y F ˜ µ F ( ω ) (cid:40) H x , y F ( e x ⊗ e x + e y ⊗ e y )+ H z F e z ⊗ e z (cid:41) (A9)By plugging the explicit expressions (A9) into Eq. (20), thatis ˜ f i F ( ω ) = σ ( ∂ B side )[ ˜ E − ( R , ω , h )] ii ˜ U ic ( ω ) , (A10)together with the conditions (15), we obtain the following re-lation making explicit the dependency on the viscosity:˜ f i F ( ω ) = H i F ˜ µ F ( ω ) ˜ U ic ( ω ) , i = x , y , z . (A11) Appendix B: Series representation of the relaxation functions
Let us first analyze the correlation functions originatingfrom the GLE in (44).As already shown in Refs. [33, 69], the two-time correla-tion functions (cid:104) U c ( t ) · U c ( t (cid:48) ) (cid:105) and (cid:104) X c ( t ) · X c ( t (cid:48) ) (cid:105) , can beobtained starting from Eqs. (52) and (57). Since both X c and U c belong to a two dimensional vector space, we have: (cid:104) U c ( t ) · U c ( t (cid:48) ) (cid:105) = g ( t ) g ( t (cid:48) ) (cid:0) M (cid:104)(cid:107) U c ( ) (cid:107) (cid:105) − M K B T (cid:1) + K B T g ( | t − t (cid:48) | ) , (B1) (cid:104) X c ( t ) · X c ( t (cid:48) ) (cid:105) = (cid:104) X c ( ) · U c ( ) (cid:105) ( G ( t ) + G ( t (cid:48) ))+ (cid:0) M (cid:104)(cid:107) U c ( ) (cid:107) (cid:105) − M k B T (cid:1) g ( t ) g ( t (cid:48) )+ K B T ( I ( t ) + I ( t (cid:48) ) − I ( | t − t (cid:48) | ))+ (cid:104)(cid:107) X c ( ) (cid:107) (cid:105) . (B2) We can now provide the expression for the MSD and theVACF. Setting t (cid:48) = t (cid:48) = t in Eq. (B2)and using the thermal initial conditions X c ( ) = (cid:104)(cid:107) U c ( ) (cid:107) (cid:105) = K B T / M one obtains: C v ( t ) = K B T g ( t ) , (cid:104) ∆ X c (cid:105) = K B T I ( t ) , (B3)where we used I ( ) = g ( ) = / M that will be provedlater by Eq. ((B13)). The time-dependent diffusion coefficientis defined by: D ( t ) : = d (cid:104) ∆ X c (cid:105) dt ( t ) , (B4)then, explicitly, we have: D ( t ) = K B T H ( t ) . (B5)
1. Computation of short time behavior of the relaxationfunctions
We consider the short time asymptotic g -relaxation functiongiven by (66):˜ g S ( s ) = Ms + ξ s + ξ p τ ν s − ν = M − s ν s ν + + ω s s ν + ω p τ − ν . (B6)This s -dependent function can be immediately transformedusing the formula [70] L − (cid:104) s ρ − s α + As ν + B (cid:105) = t α − ρ ∞ ∑ r = ( − A ) r t ( α − ν ) r E r + α , α + − ρ +( α − ν ) r ( − Bt α ) . (B7)Indeed, by setting A = ω s , B = ω p τ ν , α = ν + , ρ = ν + , in formula (B7), we obtain g S ( t ) = M ∞ ∑ r = ( − A ) r t r E r + ν + , r + ( − Bt ν + ) . (B8)6We apply once more relation (62) to the ML function obtain-ing: g S ( t ) = M ∞ ∑ r = ( − At ) r E r + ν + , r + ( − Bt ν + ) ≈ M ∞ ∑ r = ( − At ) r (cid:34) Γ ( r + ) − ( r + ) Bt ν + Γ ( r + ( ν + )) (cid:35) = M ∞ ∑ r = ( − At ) r Γ ( r + ) − Bt ν + M ∞ ∑ r = ( r + )( − At ) r Γ ( r + ( ν + ))= M (cid:34) e − At − Bt ν + ∞ ∑ r = ( r + )( − At ) r Γ ( r + ( ν + )) (cid:35) = : g S ( t ) . (B9)Noting that ( r + )( − At ) r = − A ddt ( − At ) r + , we can rewrite the series in the last equality in Eq. (B9) asfollows: g S ( t ) = M (cid:34) e − At − Bt ν + ∞ ∑ r = ( r + )( − At ) r Γ ( r + ( ν + )) (cid:35) = M (cid:34) e − At − Bt ν + ddt (cid:32) t ∞ ∑ r = ( − At ) r Γ ( r + ( ν + )) (cid:33) (cid:35) = M (cid:34) e − At − Bt ν + ddt ( t E , ν + ( − At )) (cid:35) . (B10)The term involving the derivative can be further simplified(see Eq. (B1) in [71]) yielding g S ( t ) = M (cid:34) e − At − Bt ν + ddt ( t E , ν + ( − At )) (cid:35) = M (cid:34) e − At − Bt ν + E , ν + ( − At ) (cid:35) . (B11)The form of other RFs are obtained just integrating the previ-ous one (see Eq. (B2) in Ref. [71]) H S ( t ) = M (cid:34) − e − At A − Bt ν + E , ν + ( − At ) (cid:35) I S ( t ) = M (cid:34) At + e − At − A − Bt ν + E , ν + ( − At ) (cid:35) . (B12)Before to continue and make a Taylor expansion of the RFs,it is worth noting that the only parameter appearing in the ar-guments of the exponential and in M-L functions is A = ω s .This means that it is the only characteristic frequency govern-ing the transition to subdiffusive regime.Let us go back to the Taylor expansion of the RFs, by using the representation formula for the 2-parameter ML functionand keeping just the first two terms, i.e.: E , β ( − At ) = ∞ ∑ r = ( − At ) r Γ ( r + β ) ≈ Γ ( β ) − At Γ ( β + ) , we obtain, neglecting again higher order terms, g S ( t ) ≈ M (cid:34) − At (cid:35) , H S ( t ) ≈ M (cid:34) t − At (cid:35) , I S ( t ) ≈ M (cid:34) t − At (cid:35) . (B13)
2. Computation of long time behavior of the relaxationfunctions
Let us now compute the inverse Laplace transform of thelong time asymptotic relaxation functions defined by Eqs.(71) that we recall below:˜ g B ( s ) : = ξ s + φτ ν s δλ − ν ( τ − λ + s λ ) δ , ˜ H B ( s ) : = ξ s s + φτ ν s δλ − ν + ( τ − λ + s λ ) δ , ˜ I B ( s ) : = ξ s s + φτ ν s δλ − ν + ( τ − λ + s λ ) δ . (B14)By focusing on the relaxation function ˜ g B , we can computethe inverse Laplace transform by first rewriting ˜ g B in terms ofits series representation and then employing the relation[33]: L − (cid:34) s λγ − α ( s λ + A ) γ (cid:35) = t α − E γλ , α ( − At λ ) . (B15)Indeed, setting: α = ν n , γ = δ n , A = − τ − λ , (B16)we obtain: g B ( t ) = ξ s ∞ ∑ n = (cid:18) − φτ ν (cid:19) n t ν n − E δ n λ , ν n (cid:20) − (cid:16) t τ (cid:17) λ (cid:21) . (B17)Now, using again the relation (B2) given in Ref. [71], we canwrite the remaining relaxation functions: H B ( t ) = ξ s ∞ ∑ n = (cid:18) − φτ ν (cid:19) n t ν n E δ n λ , ν n + (cid:20) − (cid:16) t τ (cid:17) λ (cid:21) , I B ( t ) = ξ s ∞ ∑ n = (cid:18) − φτ ν (cid:19) n t ν n + E δ n λ , ν n + (cid:20) − (cid:16) t τ (cid:17) λ (cid:21) . (B18)7Then, the asymptotic behavior of ML function for large timeis given by [71]: E ρα , β ( − z ) ≈ z − ρ Γ ( β − ρα ) − z − ( ρ + ) Γ ( β − ( ρ + ) α ) , ( | z | (cid:29) ) . (B19)Then, substituting such an expression into (B18), the relax-ation functions read: g B ( t ) ≈ t − ξ s E ν − δλ , (cid:16) − φ τ δλ − ν t ν − δλ (cid:17) − t − λ − τ − λ ξ p E ν − δλ , − λ (cid:16) − φ τ δλ − ν t ν − λδ (cid:17) , H B ( t ) ≈ ξ s E ν − δλ , (cid:16) − φ τ δλ − ν t ν − δλ (cid:17) , I B ( t ) ≈ t ξ s E ν − δλ , (cid:16) − φ τ δλ − ν t ν − δλ (cid:17) , (B20)where we kept also the second order of Eq. (B19) in the ex-pansion of g -RF (see Eq. (B17)).By being interested in the asymptotic behavior at very largetime, we can apply again the expansion (B19) to the 2-parameters Mittag-Leffler function in (B20). In fact, setting ρ =
1, the representations in (B20) can be rewritten as fol-lows: g B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν − Γ ( δ λ − ν ) − τ λ φ ξ s t − λ − Γ ( − λ + ( δ λ − ν )) , H B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν Γ (( δ λ − ν ) + ) , I B ( t ) ≈ τ − δλ + ν φ ξ s t δλ − ν + Γ (( δ λ − ν ) + ) . (B21) Appendix C: Numerical details
In this appendix, we report further details of the numericalmethods presented in Section III A. A snapshot of the sim-ulated system is shown in Fig. C.1, while details about theconstituents of the system, the membrane composition, the ge-ometrical parameters of the system, and the trajectories lengthare reported in Tables C.2, C.1, C.3, and C.4, respectively.FIG. C.1: Snapshot of the simulated system. Water and ionsare not illustrated for the sake of clarity.
Lipid Type Upperleaflet Nº Lowerleaflet Nº
CHOL 300 300DOPC 12 12DPPC 30 30POPC 54 54DGPE 12 12POPS 18 18PAPI 18 18DPSM 54 54PNSM 18 18POSM 6 6
TABLE C.1: The number of lipids for each species.
Nº of water molecules 14077Nº of Na + ions 825Nº of Cl − ions 161Nº of lipids 1044Nº of proteins 1 TABLE C.2: The number of each constituent of the system.
Appendix D: Data Analysis
In order to evaluate best fit values of the model parametersand their uncertainties, we used a Bayesian approach with aflat (uninformative) prior distribution for the parameters [63],i.e. all the parameters values are assumed to be equally prob-able before analyzing the data.According with the Bayes’s theorem, the posterior distribu-tion of a model described by the parameters θ , given the mea-sured data d , reads: P ( θ | d ) = L ( d | θ ) Π ( θ ) E ( d ) , (D1)where L ( d | θ ) = Z e − χ / = Z exp (cid:34) − N ∑ i = (cid:18) m ( θ , i ) − d i σ i (cid:19) (cid:35) , (D2)is the likelihood of the MSD data d = { d i } , given the MSDmodel m ( θ , i ) calculated for the parameters θ , with i being theindex running over the points of the data set of size N , and σ i the error of the data. Z is a normalization constant that guar-antees that (cid:82) L ( d | θ ) d θ =
1. The evidence of the data, E ( d ) ,amounts to a normalization constant whose value is fixed and Π ( θ ) is the prior distribution of parameters, which is assumedto be uniform [63].To obtain an accurate estimation of the subdiffusive to Brow-nian dynamics crossover, many microscopic fluctuations mustbe observed. Thus, the relative error on the observed MSDhas been estimated as the reciprocal of the square root of thenumber of observed fluctuations and so as ( τ ∗ / T ) / , where τ ∗ is the characteristic time of the fully Brownian dynamicsrecovering, ≈
500 ns, and T is the observation time (i.e. thelargest trajectory length).The posterior distribution of the model is sampled over a8FIG. D.1: One-dimensional and 2-dimensional marginal distributions. The highlighted area below the one-dimensionalmarginal distributions corresponds to 95%. The ellipses in the 2-dimensional marginal distributions correspond to 68% and95% confidence regions. The reduced ˜ χ is 1.03. Upperleaflet Lowerleaflet
Protein Area 1152.12249 Å Lipid Area 26331.6 Å Total Area 27483.72249 Å Protein X Extent 26.06 ÅProtein Y Extent 24.69 ÅBox Y Extent 165.74 ÅBox X Extent 165.74 ÅBox Z Extent 109.31 Å
TABLE C.3: Geometrical parameters of the simulatedsystem.
Regimes Sampling time step Time lengthBallistic 20 fs 1000 psSubdiffusive 20 ps 3 µ sDiffusive 5000 ps 100 µ s TABLE C.4: Sampling time step and total length of thenumerical simulations. regular grid of the parameter space, ω s ∈ [ . , . ] ps − , ω p ∈ [ . , . ] ps − , τ ∈ [ , ] ps, and λ (= ν ) ∈ [ . , . ] . The expectation value of each parameter isestimated by the maximum of the one-dimensional marginal-ized posterior distribution, i.e. after marginalizing P overall but one parameter in turn (Fig. D.1). The parametersconfidence intervals correspond to the 95% of the area be-low the one-dimensional marginalized posterior distributionstarting from the maximum. In order to estimate the correla-tions between the parameters of the model, the 2-dimensionalmarginal distribution for each pair of parameters (obtained bymarginalizing over n − n is the number offree parameters) is calculated as well (Fig. D.1). [1] T Springer and RE Lechner. Diffusion in condensed matter ,volume 1. Springer New York, 2005. [2] Paolo Mereghetti, Razif R Gabdoulline, and Rebecca C Wade. Brownian dynamics simulation of protein solutions: structuraland dynamical properties.
Biophysical journal , 99(11):3782–3791, 2010.[3] Ann E Cowan, Ion I Moraru, James C Schaff, Boris MSlepchenko, and Leslie M Loew. Spatial modeling of cell sig-naling networks. In
Methods in cell biology , volume 110, pages195–221. Elsevier, 2012.[4] Takahiro Fujiwara, Ken Ritchie, Hideji Murakoshi, Ken Jacob-son, and Akihiro Kusumi. Phospholipids undergo hop diffusionin compartmentalized cell membrane.
The Journal of cell biol-ogy , 157(6):1071–1082, 2002.[5] Kotono Murase, Takahiro Fujiwara, Yasuhiro Umemura,Kenichi Suzuki, Ryota Iino, Hidetoshi Yamashita, MihokoSaito, Hideji Murakoshi, Ken Ritchie, and Akihiro Kusumi.Ultrafine membrane compartments for molecular diffusion asrevealed by single molecule techniques.
Biophysical journal ,86(6):4075–4093, 2004.[6] Aubrey V Weigel, Blair Simon, Michael M Tamkun, and DiegoKrapf. Ergodic and nonergodic processes coexist in the plasmamembrane as observed by single-molecule tracking.
Proceed-ings of the National Academy of Sciences , 108(16):6438–6443,2011.[7] Ellen Gielen, Jo Vercammen, Jan S`ykora, Jana Humpolickova,Alés Benda, Niels Hellings, Martin Hof, Yves Engelborghs,Paul Steels, Marcel Ameloot, et al. Diffusion of sphingomyelinand myelin oligodendrocyte glycoprotein in the membrane ofOLN-93 oligodendroglial cells studied by fluorescence corre-lation spectroscopy.
Comptes rendus biologies , 328(12):1057–1064, 2005.[8] Margaret R Horton, Felix Höfling, Joachim O Rädler, andThomas Franosch. Development of anomalous diffusion amongcrowding proteins.
Soft Matter , 6(12):2648–2656, 2010.[9] Sivaramakrishnan Ramadurai, Andrea Holt, Victor Krasnikov,Geert van den Bogaart, J Antoinette Killian, and Bert Poolman.Lateral diffusion of membrane proteins.
Journal of the Ameri-can Chemical Society , 131(35):12650–12656, 2009.[10] MA Deverall, E Gindl, E-K Sinner, H Besir, J Ruehe,Michael J Saxton, and CA Naumann. Membrane lateral mobil-ity obstructed by polymer-tethered lipids studied at the singlemolecule level.
Biophysical journal , 88(3):1875–1886, 2005.[11] Jae-Hyung Jeon, Matti Javanainen, Hector Martinez-Seara,Ralf Metzler, and Ilpo Vattulainen. Protein crowding in lipidbilayers gives rise to non-Gaussian anomalous lateral diffusionof phospholipids and proteins.
Physical Review X , 6(2):021006,2016.[12] István P Sugár and Rodney L Biltonen. Lateral diffusionof molecules in two-component lipid bilayer: a Monte Carlosimulation study.
The Journal of Physical Chemistry B ,109(15):7373–7386, 2005.[13] Michael J Skaug, Roland Faller, and Marjorie L Longo. Corre-lating anomalous diffusion with lipid bilayer membrane struc-ture using single molecule tracking and atomic force mi-croscopy.
The Journal of chemical physics , 134(21):06B602,2011.[14] Axel Kammerer, Felix Höfling, and Thomas Franosch. Cluster-resolved dynamic scaling theory and universal corrections fortransport on percolating systems.
EPL (Europhysics Letters) ,84(6):66002, 2008.[15] Michael J Saxton. Anomalous diffusion due to obstacles: aMonte Carlo study.
Biophysical journal , 66(2):394–401, 1994.[16] Jae-Hyung Jeon, Hector Martinez-Seara Monne, Matti Ja-vanainen, and Ralf Metzler. Anomalous diffusion of phospho-lipids and cholesterols in a lipid bilayer and its origins.
Physicalreview letters , 109(18):188103, 2012. [17] Matti Javanainen, Henrik Hammaren, Luca Monticelli, Jae-Hyung Jeon, Markus S Miettinen, Hector Martinez-Seara, RalfMetzler, and Ilpo Vattulainen. Anomalous and normal diffusionof proteins and lipids in crowded lipid membranes.
Faradaydiscussions , 161:397–417, 2013.[18] Gerald R Kneller. Communication: A scaling approach toanomalous diffusion, 2014.[19] Felix Höfling and Thomas Franosch. Anomalous transport inthe crowded world of biological cells.
Reports on Progress inPhysics , 76(4):046602, 2013.[20] Michael J Saxton. A biological interpretation of transientanomalous subdiffusion. I. Qualitative model.
Biophysical jour-nal , 92(4):1178–1191, 2007.[21] Daniel Molina-Garcia, Trifce Sandev, Hadiseh Safdari, GianniPagnini, Aleksei Chechkin, and Ralf Metzler. Crossover fromanomalous to normal diffusion: truncated power-law noise cor-relations and applications to dynamics in lipid bilayers.
NewJournal of Physics , 20(10):103027, 2018.[22] Robert Zwanzig.
Nonequilibrium statistical mechanics . OxfordUniversity Press, 2001.[23] R Metzler, J-H Jeon, and AG Cherstvy. Non-Browniandiffusion in lipid membranes: Experiments and simula-tions.
Biochimica et Biophysica Acta (BBA)-Biomembranes ,1858(10):2451–2467, 2016.[24] André Liemert, Trifce Sandev, and Holger Kantz. GeneralizedLangevin equation with tempered memory kernel.
Physica A:Statistical Mechanics and its Applications , 466:356–369, 2017.[25] Jae-Hyung Jeon, Hector Martinez-Seara Monne, Matti Ja-vanainen, and Ralf Metzler. Anomalous diffusion of phospho-lipids and cholesterols in a lipid bilayer and its origins.
Physicalreview letters , 109(18):188103, 2012.[26] Shao-Hua Wu, Shalene Sankhagowit, Roshni Biswas, ShuyangWu, Michelle L Povinelli, and Noah Malmstadt. Viscoelasticdeformation of lipid bilayer vesicles.
Soft Matter , 11(37):7385–7391, 2015.[27] GE Crawford and JC Earnshaw. Viscoelastic relaxation of bi-layer lipid membranes. Frequency-dependent tension and mem-brane viscosity.
Biophysical journal , 52(1):87, 1987.[28] DS Dimitrov. Electric field-induced breakdown of lipid bilay-ers and cell membranes: a thin viscoelastic film model.
TheJournal of membrane biology , 78(1):53–60, 1984.[29] Mohammad Rahimi and Marino Arroyo. Shape dynamics,lipid hydrodynamics, and the complex viscoelasticity of bilayermembranes.
Physical review E , 86(1):011932, 2012.[30] Mirko D’Ovidio and Federico Polito. Fractional diffusion-telegraph equations and their associated stochastic solutions. arXiv preprint arXiv:1307.1696 , 2013.[31] M Caputo. Elasticitá e dissipazione (Elasticity and anelasticdissipation).
Zanichelli, Bologna , 1969.[32] Tilak Raj Prabhakar et al. A singular integral equation with ageneralized Mittag Leffler function in the kernel. 1971.[33] Trifce Sandev, Živorad Tomovski, and Johan LA Dubbeldam.Generalized Langevin equation with a three parameter Mittag-Leffler noise.
Physica A: Statistical Mechanics and its Appli-cations , 390(21-22):3627–3636, 2011.[34] Yilei Wu, Martin Štefl, Agnieszka Olzy´nska, Martin Hof,Gokhan Yahioglu, Philip Yip, Duncan R Casey, Oscar Ces, JanaHumpolíˇcková, and Marina K Kuimova. Molecular rheometry:direct determination of viscosity in L o and L d lipid phases viafluorescence lifetime imaging.
Physical Chemistry ChemicalPhysics , 15(36):14986–14993, 2013.[35] Chao-Cheng Wang.
Mathematical principles of mechanics andelectromagnetism: Part A: Analytical and continuum mechan-ics , volume 16. Springer Science & Business Media, 2013. [36] Carlos Conca, H Jorge San Martín, and Marius Tucsnak. Mo-tion of a rigid body in a viscous fluid. Comptes Rendus del’Académie des Sciences-Series I-Mathematics , 328(6):473–478, 1999.[37] Nikolai V Chemetov and Šárka Neˇcasová. The motion of therigid body in the viscous fluid including collisions. Global solv-ability result.
Nonlinear Analysis: Real World Applications ,34:416–445, 2017.[38] Morton E Gurtin.
An introduction to continuum mechanics .Academic press, 1982.[39] Jerrold E Marsden and Alexandre J Chorin.
A mathematicalintroduction to fluid mechanics . Springer-Verlag, 1993.[40] Ray M Bowen.
Introduction to continuum mechanics for engi-neers . Plenum Press, 1989.[41] Daniel D Joseph.
Fluid dynamics of viscoelastic liquids , vol-ume 84. Springer Science & Business Media, 2013.[42] The physical reason of this assumption will be justified later inSection II.[43] Shuvojit Paul, Basudev Roy, and Ayan Banerjee. Free and con-fined Brownian motion in viscoelastic Stokes–Oldroyd B fluids.
Journal of Physics: Condensed Matter , 30(34):345101, 2018.[44] William Bailey Russel, WB Russel, Dudley A Saville, andWilliam Raymond Schowalter.
Colloidal dispersions . Cam-bridge university press, 1991.[45] Roberto Mauri.
Non-Equilibrium Thermodynamics in Multi-phase Flows . Springer Science & Business Media, 2012.[46] Giovanni Galdi.
An introduction to the mathematical theory ofthe Navier-Stokes equations: Steady-state problems . SpringerScience & Business Media, 2011.[47] Gabriel Espinosa, Iván López-Montero, Francisco Monroy, andDominique Langevin. Shear rheology of lipid monolayers andinsights on membrane fluidity.
Proceedings of the NationalAcademy of Sciences , 108(15):6008–6013, 2011.[48] Francesco Mainardi. Fractional calculus: theory and applica-tions, 2018.[49] Andrea Giusti and Ivano Colombaro. Prabhakar-like fractionalviscoelasticity.
Communications in Nonlinear Science and Nu-merical Simulation , 56:138–143, 2018.[50] Yifei Qi, Helgi I Ingólfsson, Xi Cheng, Jumin Lee, Siewert JMarrink, and Wonpil Im. CHARMM-GUI Martini maker forcoarse-grained simulations with the Martini force field.
Journalof chemical theory and computation , 11(9):4486–4494, 2015.[51] Siewert J Marrink, H Jelger Risselada, Serge Yefimov, D PeterTieleman, and Alex H De Vries. The Martini force field: coarsegrained model for biomolecular simulations.
The journal ofphysical chemistry B , 111(27):7812–7824, 2007.[52] Siewert J Marrink, Alex H De Vries, and Alan E Mark. Coarsegrained model for semiquantitative lipid simulations.
The Jour-nal of Physical Chemistry B , 108(2):750–760, 2004.[53] Kazuko Haga, Andrew C Kruse, Hidetsugu Asada, TakamiYurugi-Kobayashi, Mitsunori Shiroishi, Cheng Zhang,William I Weis, Tetsuji Okada, Brian K Kobilka, Tatsuya Haga,et al. Structure of the human M2 muscarinic acetylcholinereceptor bound to an antagonist.
Nature , 482(7386):547–551,2012.[54] Robin B Chan, Tiago G Oliveira, Etty P Cortes, Lawrence SHonig, Karen E Duff, Scott A Small, Markus R Wenk, Guanghou Shui, and Gilbert Di Paolo. Comparative lipidomicanalysis of mouse and human brain with Alzheimer disease.
Journal of Biological Chemistry , 287(4):2678–2688, 2012.[55] CHARMM-GUI web-site. .[Online].[56] Herman JC Berendsen, David van der Spoel, and Rudi vanDrunen. GROMACS: a message-passing parallel moleculardynamics implementation.
Computer physics communications ,91(1-3):43–56, 1995.[57] Erik Lindahl, Berk Hess, and David Van Der Spoel. GRO-MACS 3.0: a package for molecular simulation and trajectoryanalysis.
Molecular modeling annual , 7(8):306–317, 2001.[58] David Van Der Spoel, Erik Lindahl, Berk Hess, Gerrit Groen-hof, Alan E Mark, and Herman JC Berendsen. GROMACS:fast, flexible, and free.
Journal of computational chemistry ,26(16):1701–1718, 2005.[59] Sander Pronk, Szilárd Páll, Roland Schulz, Per Larsson, PärBjelkmar, Rossen Apostolov, Michael R Shirts, Jeremy CSmith, Peter M Kasson, David van der Spoel, et al. GROMACS4.5: a high-throughput and highly parallel open source molec-ular simulation toolkit.
Bioinformatics , 29(7):845–854, 2013.[60] Szilárd Páll and Berk Hess. A flexible algorithm for calculat-ing pair interactions on SIMD architectures.
Computer PhysicsCommunications , 184(12):2641–2650, 2013.[61] Alfredo Di Nola, Danilo Roccatano, and Herman JC Berend-sen. Molecular dynamics simulation of the docking of sub-strates to proteins.
Proteins: Structure, Function, and Bioin-formatics , 19(3):174–182, 1994.[62] Michele Parrinello and Aneesur Rahman. Polymorphic tran-sitions in single crystals: A new molecular dynamics method.
Journal of Applied physics , 52(12):7182–7190, 1981.[63] Edwin T Jaynes.
Probability theory: The logic of science . Cam-bridge university press, 2003.[64] Gábor Horváth, Illés Horváth, Salah Al-Deen Almousa, andMiklós Telek. Numerical inverse Laplace transformation us-ing concentrated matrix exponential distributions.
PerformanceEvaluation , 137:102067, 2020.[65] Inverse Laplace transform package. http://inverselaplace.org . [Online].[66] GitHub Repository with inverse Laplace transform codes. https://github.com/ghorvath78/iltcme . [Online].[67] Eric-Jan Wagenmakers. A practical solution to the perva-sive problems ofp values.
Psychonomic bulletin & review ,14(5):779–804, 2007.[68] RCSB PDB. . [On-line].[69] Noëlle Pottier. Aging properties of an anomalously diffusingparticule.
Physica A: Statistical Mechanics and its Applica-tions , 317(3-4):371–382, 2003.[70] R Figueiredo Camargo, Ary O Chiacchio, R Charnet, andE Capelas de Oliveira. Solution of the fractional Langevin equa-tion and the Mittag–Leffler functions.
Journal of mathematicalphysics , 50(6):063507, 2009.[71] R Figueiredo Camargo, E Capelas de Oliveira, and J Vaz Jr. Onanomalous diffusion and the fractional generalized Langevinequation for a harmonic oscillator.