The interplay between membrane viscosity and ligand-binding receptor kinetics in lipid bilayers

Plasma membranes appear as deformable systems wherein molecules are free to move and diffuse giving rise to condensed microdomains (composed of ordered lipids, transmembrane proteins and cholesterol) surrounded by disordered lipid molecules. Such denser and thicker regions, namely lipid rafts, are important communication hubs for cells. Indeed, recent experiments revealed how the most of active signaling proteins co-localize on such domains, thereby intensifying the biochemical trafficking of substances. From a material standpoint, it is reasonable to assume the bilayer as a visco-elastic body accounting for both in-plane fluidity and elasticity. Consequently, lipid rafts contribute to membrane heterogeneity by typically exhibiting higher stiffness and viscosity and by locally altering the bilayer dynamics and proteins activity. A chemo-mechanical model of lipid bilayer coupled with inter-specific dynamics among the resident species (typically transmembrane receptors and trasporters) has been recently formulated to explain and predict how proteins regulate the dynamic heterogeneity of membrane. However, the explicit inclusion of the membrane viscosity in the model was not considered. To this aim, the present work enriches the constitutive description of the bilayer by modeling its visco-elastic behavior. This is done through a strain-level dependent viscosity able to theoretically


List of symbols and definitions
Applied membrane pressure

Introduction
Early findings assumed the eukaryotic cell membranes as a bi-dimensional assembly of lipids organized in a fluid bilayer where transmembrane proteins can laterally diffuse [1].Lipids self-assemble in a ∼ 5nm thick bilayer [2] and achieve an areal stretch of the order of 5% [3].Phospholipids can move in the planar direction and, so, plasma membranes are characterized by quasi-fluid deformable surfaces that express solidfluid-like behavior, resulting in systems wherein in-plane fluidity and elasticity may simultaneously emerge [4].Such fluidity is measured through the viscosity, whose available literature data are, however, highly experiment dependent, sometimes varying by orders of magnitude [5].A possible explanation for this huge variability could be that membrane surface viscosity is a macroscopic quantity modeled at scales where the bilayer is assumed to behave like a 2-dimensional quasiincompressible fluid.For this reason, micro-or nano-scale measurements may not be sufficient to catch the effective continuum viscosity but, rather, the so-called "microviscosity".The latter is a local quantity influenced by the environment [6].
Membrane fluidity is therefore associated with the high molecular mobility inside the lipid bilayer, enabling for a lateral diffusion of the embedded proteins [7].Hence, viscosity results to be measured through the estimation of lipid diffusion coefficient [5].It is indeed confirmed that the ligand-binding of receptors -as for example the G-Protein Coulped Receptors (GPCRs)-requires the presence of molecules that are able to move within the membrane [8].In this regard, it has been established the difference, in terms of viscosity, among the resistance to flow under an applied shear stress and the capability of molecules to move and diffuse inside the membrane [9].In the latter, it has been demonstrated that high diffusion mobility could be linked to a finite macroscopic shear viscosity, however discussing many cases of gel-phase of single saturated phospholipids or solid ceramide lipids that are able to pack themselves into a solid structure with high shear stiffness and viscosity.Quantitative stability analyses of viscoelastic lipid bilayers with properties deduced by [9], have been provided in [10].Furthermore, in complex bio-membranes gel domains may coexist with fluid ones, thus promoting regions with vastly distinct viscosities [11].Actually, evidences show that the mammalian cell membrane has a time-varying force response as nonlinear function of strain, so behaving as a visco-elastic or non-Newtonian fluid [12].Related to this phenomenology, one can recall that lipid bilayers undergo various stages at which they may experience area expansion, thereby responding with compression and shear moduli [9].Such a variation in the local mechanical properties seems to be responsible for the majority of cellular processes [13].
Several experimental strategies have been used to quantify the dynamical visco-elasticity of lipid systems [14,15].Recently, AFM measurements were performed to capture both the elastic and viscous properties of lipid systems that resulted to affect the propagation or attenuation of mechanosignaling across the cell membrane [16].Also, high frequency experiments, modeled through a continuum mechanical theory, revealed that the plasma membrane displays a visco-elastic behavior [17].In particular, it has been estimated that the cell surface responds like an elastic material on short time scales of around 1s, while exhibiting properties of a viscous body on longer time scales ∼ 10 − 100s [18].Bulk membrane viscosity and transverse stiffness are therefore correlated but also influenced by lipid packing density [19].
Modulation of membrane behavior has been demonstrated to be fundamental in various diseases [20][21][22][23][24].For instance, it is indeed confirmed that changes in membrane viscosity influence the evolution of the metastatic progression of cancerous cells [25,26].In [27] it is shown that the latter are softer than healthy cells and that they are also characterized by a more fluid membrane.For these reasons, the measure of membrane viscoelasticity leads to the possibility of discriminating between normal and cancerous cells through the application of multi-frequency vibrations [17].
Lipid rafts have been demonstrated to be involved in cardiovascular signaling as determinant regulators of vascular endothelial and smooth muscle cells, and in particular in signal transduction across the plasma membrane, of primary importance to many functional activities.At present, little is known about the specific role of lipid rafts in cardiac function and dysfunction, increasing attention focusing on their contribution to the pathogenesis of several structural and functional processes including cardiac hypertrophy and heart failure, as well as atherosclerosis, ischemic injury and different myocardial functions [28].Lipid rafts in cardiomyocyte membranes are enriched in signaling molecules and ion channel regulatory proteins, therefore contributing to calcium handling and Ca2+ entry that control excitation-contraction of heart muscle cells.Thus, they can actively participate in differential cardiomyocyte ion channel targeting and regulation [28,29].
Ordered microdomains result fundamental to stabilize signal transduction activities required for angiogenesis.In fact, it has been observed that VEGF receptor-2 (VEGFR-2), which stimulates angiogenic signaling, co-localizes with lipid rafts to regulate its activation.Also, long-term VEGFR2 relocation closely depends on lipid raft integrity, disruption of lipid rafts directly causing receptors' depletion and inefficacy.In this sense, therapeutic strategies are more and more oriented towards the possible modulation of lipid rafts to control cells' sensitivity to VEGF expression [30,31].Also, GPCRs have a primary influence in cardiac remodeling.Activation of epidermal growth factor receptors is in fact mediated by a large repertoire of GPCRs in the heart, and promotes cardiomyocyte survival, thus suggesting innovative therapeutic scenarios based on their targeting [32,33].
Despite available pure mechanical descriptions of the lipid bilayers [34,35] or purely diffusive approaches where the influence of micromechanical stimuli is neglected [36], there is still no modeling approach that takes into account the synergistic influence of membrane viscosity on transmembrane proteins activation and mobility and/or viceversa the role of proteins and lipids in membrane fluidity.Actually, it is well known that physical and chemical events act together to form the complexity of processes responsible for cell functions [37].Therefore, a multiphysics analysis becomes manifest to provide new insights into the very complex world of plasma membranes.In this regard, mathematical production provided in Carotenuto et al. [38] confirmed the common knowledge that active receptors prefer to cluster on the so-called lipid rafts -wherein high cholesterol concentration increases bilayer rigidity [39]-through a chemo-mechanical coupled model.In [38], the model was regulated by the coupling of the membrane remodeling and its energetics dependent on the active proteins involved in the system, i.e. β2−adrenergic receptors.Moreover, recent findings [40] highlighted the effects produced by the receptors and transporters on raft formation and coalescence through Cahn-Hilliardtype dynamics in a two-dimensional hyper-elastic framework.
Neverthless, as aforementioned the lipid bilayer is characterized by viscous properties and so, in order to obtain a more faithful solidliquid description of this kind of system, a viscohyperelastic model should be considered.This may provide an explicit interaction between the characteristic time evolution of the populations of transmembrane proteins and the relaxation time of the lipid bilayer.This is because, at the microscopic level, single protein re-arrangement and configurational changes are known to occur within milliseconds and are likely to locally produce elastic pressures at the membrane-protein interfaces [41,42].This can be extended at the population level through the presented continuum approaches, in which the dynamics of entire protein clusters is followed in response to the ligand time-varying precipitation stimulus.The morphoelastic reconfiguration of the membrane thus can produce maps of heterogeneous stress and deformation that could project at the continuum scale the instantaneous packing of lipids and protein activation occurring within the ordered phase.
All this considered, the aim of the present study is to enrich well-grounded hyper-elastic models [38,[43][44][45] of cell membranes by incorporating a material viscous component in the constitutive model.This provides an explicit interaction between the characteristic time evolution of the population of transmembrane proteins and the relaxation time of the lipid bilayer, by so calling into play a possible competition between the pseudo-viscous and the characteristic viscous terms.

Chemo-Mechanical characterization of the membrane behavior
It is well established that the plasma membrane undergoes a thickness change due to an ordereddisordered phase transition occurring at the lipid scale.This thickness variation is mainly caused by the lipid re-arrangement that, in assuming an ordered configuration, have straightened tails and appear tightly packed together as it occurs in functional micro-domains of the lipid membrane denoted as raft phase [46].Several approaches have been adopted to analyze the mechanical behavior of membrane systems when experience phase transition based on either molecular dynamics simulations or, at the continuum scale, phase separation and elasticity models [47][48][49][50].Recently, a nonlinear hyperelastic response of the plasma membrane has been used to build up a fullycoupled framework describing the membrane's macroscopic remodeling and functional reorganization as regulated by the leading biochemical events occurring among interacting protein species in forming lipid raft domains [38].In the subsequent work by Bernard et al. [40], this evolutionary approach has been further enriched by Cahn-Hilliard energetics and kinetics for the involved species, thereby accounting for rafts nucleation and coalescence.The time-varying nature of the involved biological species associated to configurational remodeling terms gave to the system a pseudo-visco-elastic nature (with eventual dissipation), the rate of the internal species kindling a viscous-type (chemical) stress.However, in [40] the explicit role of intrinsic visco-elasticity of the lipid membrane and the possible influence of the fluid component of the bilayer on raft development was not considered.To this purpose, we here analyze a two-dimensional system capable to experience a lipid phase separation and manifest raft coarsening within a visco-elastic environment.The whole phenomenon will be the result of the coupling between the conformational remodeling guided by the presence of the active protein species and the energetics of the membrane.In particular, the elastic part of the membrane response -in line with well-established literature [51][52][53]-is modeled by assuming a neo-Hookean type behavior [40], by neglecting for now the spontaneous trends of lipids to reorganize themselves in co-existing phases (this can be accounted for not convex energy terms [54]).At the molecular scale, the activation of a single transmembrane protein within the lipid environment provokes a re-arrangement of its sub-units, which induces a stress in the surrounding membrane in the form of an in-plane pressure.This, inevitably, calls into play the adaptation of the neighboring lipids.In the absence of any viscous component, the adaptation of the lipid membrane is entirely dictated by the dynamics of the protein populations.In this sense, at the macroscopic scale the overall deformation and morphological remodeling of the lipid membrane is seen as the averaged result of the overall behavior of protein densities.The latter will pass to their active state asynchronously by introducing delays and by exchanging (positive or negative) chemical feedbacks.These give rise to more complex spatial and temporal patterns of the membrane heterogeneity.Noteworthy, the characteristic times of the membrane evolution do not simply follow the activation times of single units (of the order of few milliseconds).Rather, instead ensue the collective dynamics of active resident proteins and their progressive recruitment.Indeed, lipid and proteins' clusters have a much larger lifespan (from seconds to several minutes [55][56][57]).In this sense, the micro-and macro-scopic scales of the ordered macro-islands could potentially describe multi-scale kinematics in a cascade manner.Through the above described mechanisms, in [40] an interspecific protein dynamics, enriched with a Cahn-Hilliard energetics and kinetics phenomena, has been adopted to successfully trace back the complex spatio-temporal adaptation of the membrane.Of course, the chemo-mechanical coupling becomes absolutely crucial to theoretically explain how protein density dynamics affects the structural remodeling of the membrane, leading to the nucleation of raft domains.The heterogeneity noticed in lipid bilayers has to be indeed addressed to the coexistence of disordered and ordered lipid phases [58].To this end, well-grounded observations show the formation of zones with different concentration levels [59].In particular, regions with high concentration of proteins have been recognized in lipid rafts [60], where the clustering phenomena give rise to the initiation of most of cellular processes [61][62][63].For this reason, the introduction of a phase-separation diffusive model able to predict coalescence of different species becomes apparent.Within this framework, the Cahn-Hilliard equation is typically used to describe two-phase separation problems [64][65][66] that are mathematically described by a diffusion equation for the species concentration [67].
In this respect, the theoretical model proposed in [40] described the evolution of protein species through Cahn-Hilliard-like energetics and kinetics wherein reaction interspecific terms account for the mutual influence among protein populations, i.e. the above mentioned GPCRs and their antagonist the Multidrug Resistance Proteins (MRPs), while non-local species momenta are enriched by strain-dependent morphotaxis terms.The latter enable the movement of protein species along the gradients of lipid order distribution, so promoting the tendency of signaling proteins to reside on raft domains by favoring spatial co-localization of such species on raft islands.When the viscous component of the membrane is introduced and a visco-elastic behavior of the membrane is considered, the above described dynamics can be altered by the direct competition between both the characteristic adaptation and the intrinsic bilayer relaxation times.Indeed, it is expected that viscosity may affect the membrane deformation triggered by proteins through creep-associated effects in raft emergence, thus so influencing its chemical stability and persistence.On the other hand, stress relaxation phenomena could occur as well by redistributing internal stresses with effect on the residual stress-induced stiffness and membrane tension.However, rough estimations of the visco-elastic and lipid raft characteristic times -respectively of microseconds and tens of seconds-would suggest that these phenomena would minimally concur together in determining the structural re-organization of the membrane.More important effects could be rather produced by the synergy of protein dynamics with nonlinear deformations and viscous response, which could lead instead to more significant changes into the material remodeling of membrane properties.This would meet some experimental evidences showing that rafts are highly viscous and stiff zones of the membrane.To do this, in what follows we present the governing equations of the coupled model within a visco-elastic framework.This will enable to investigate how membrane fluidity is influenced by the dynamical re-organization.In particular, we will initially consider the effects of a constant (i.e.linear) viscous term on raft persistence.While afterwords a strain-level dependent viscosity will be considered to explore if the increase of viscosity of heterogeneous lipid membranes plays a key influence on co-evolving with lipid rafts.

Uploading visco-elasticity in the coupled chemo-mechanical model
The lipid bilayer can be assumed as a twodimensional quasi-incompressible hyperelastic thin body, wherein areal and thickness stretches locally vary with the corresponding changes of the lipid order [51][52][53].Herein, the membrane is assumed flat in its natural configuration and its kinematics is supposed to be confined in the class of normal preserving deformations (see e.g.[34,54,68]).The natural configuration of the membrane B 0 is partitioned in a two-dimensional domain x = x e 1 + ye 2 and the thickness z.Hence, the material particles x ∈ B 0 are described as x = x + ze 3 , at time t.Accordingly, the displacement field characterizing the kinematics of the membrane can be written as follows: (1) where the function ϕ (x, y, t) represents the thickness stretch in the direction e 3 , at time t.The displacement (1) yields the deformation gradient to which the chosen strain measures, as well as strain rates, can be readily associated: By restricting the problem to the mid-plane of the membrane (see e.g.[34,54,68]) and by accounting for a volumetric incompressibility constraint restricted to such mid-plane, the determinant of F at z = 0 reads: where ϕ (x, y, t) = 1 J0 , and J 0 denotes the areal stretch in the membrane plane, i.e.J 0 = det F 0 with F 0 defined as the dimensional reduction of F on the membrane mid-plane, i.e.F 0 = 2 α,β=1 δαβ + ∂u α /∂x β e α ⊗ e β , where δαβ is the Kronecker delta.Incompressibility on the midplane also implies that tr(D) = 0, once the trace is restricted to operate on D in such a plane.
Following [40], the energetics of the system is assumed to be governed by the Helmholtz-free energy density W (F, n i , ∇n i , ϕ), where n i is the concentration of the i-th active species.Hence, by considering an additive decomposition of such energy, the contributions given by the potential associated with the hyperelastic energy of the membrane and the one related to the transmembrane proteins are introduced: Herein, the contribution W ni contains a coupling term that explicitly depends on the out-ofplane stretch ϕ, accounting for the influence that changes in species concentration have on membrane deformation and vice-versa.In fact, protein re-organization at the micro-level exerts work on the surrounding membrane, thus calling into play the bilayer deformation and stress.On this account, besides an intrinsic species-dependent energy density, Ψ ni , the potential W ni provides the coupling term due to the above mentioned interaction which reads as follows: Here w i is a coupling parameter connected to the exchange of mechanical work between activating proteins and membrane: such w i directly emerges from the sub-macroscopic scale as shown in [38].As discussed above, the energy contribution Ψ ni is actually given in terms of the Ginzburg-Landau phase separation energy [69]: defining the coefficients ϵ, γ > 0, and the gradient term ∇ n i − n 0 i so written to ensure thermodynamic consistency [40].More in detail, in relation ( 6) a double-well potential is assumed to model the energy contribution of each species in passing from the inactive to the active state.This is done by deriving conditions for chemical equilibrium that could explicitly, although phenomenologically, take into account the effect of the fundamental mechanical coupling (i.e. the second term of ( 5)), by so modifying the energetic convenience of the system.Indeed, the cell membrane undergoes shape deformations in terms of phase transition between states separated by energy barriers.
The energy landscape of lipid membranesand biphasic systems in general-is modeled by a parameterized double-well potential characterized by two fixed degenerate minima standing for the coexistence of such phases [70].In the case of the proposed model, in presence of a varying mechanical micro-environment, the membrane mechanical state directly influences the chemical activation of the protein species.More in detail, given that in a classical double-well potential the two minima uniquely identify the active/inactive state of the proteins in a completely symmetric way, the presence of the stretch-dependent coupling term here alters such symmetry.This occurs by moving the position of the minima and so determining a non-symmetric and variable convenience of certain protein species to be in their active or inactive state on the base of the surrounding conditions.This constitutes an important mechano-signaling pathway contributing to co-localization.In fact, when the transverse stretch ϕ > 1 the coupling term makes the active state more energetically favorable with respect to the inactive one.Viceversa, as the membrane is thinning (i.e.0 < ϕ < 1) the disordered state results to be more energetically convenient (see Figure 1).In this present paper, in order to characterize the elastic part of the bilayer response, a standard incompressible neo-Hookean strain energy [40,51,52] is considered: where I 1 = tr F T F is the first invariant of the Cauchy-Green strain tensor and is the tangent shear modulus with the Poisson's ratio ν approaching 0.5 due to the incompressibility constraint, and p is the associated lagrangian pressure.Consistency with linear elasticity, suggests a finite value of the elastic modulus G, as these two material constants are connected to each other through well-established Lamé relations.This is done coherent with evidence arising while observing that lipid bilayers may possess rigidity and elastic compressibility [9].In fact, as reported in Espinosa et al. [9], biological membranes -for which fluidity is associated to the high molecular mobility inside the lipid bilayer enabling for a lateral diffusion of the embedded proteinsalso can account for a nonzero shear modulus as structural intrinsic property needed for biological functions.Moreover, in the light of thermodynamics, as in [40] it is possible to introduce specific constitutive assumptions upon which one can evaluate the stresses and the chemical potentials associated to each protein species in the presence of the chemomechanical coupling.In doing this, it is assumed that the kinematics of the remodeling membrane provides a multiple configuration path, in which the membrane is first hypothesized to undergo a geometry-preserving activation step (see Fig. 2).There, part of the proteins pass to the active state by experiencing conformational switches at the sub-macroscopic scale [38].At the macro-scale, this virgin-to-active state can be attained through a jacobian remodeling term, say K r , derived in the framework of Structured Deformations [71][72][73][74][75][76].More in detail, this remodeling is due to submacroscopic re-arrangements of lipids clusters incorporating activated receptors.Obviously, the latter activates through conformational changes of some of their transmembrane domains during ligand-binding across the membrane.Thus, this depends on the amount of proteins entering the active state and it can be derived by imposing mass conservation between the virgin configuration -where material points have a virgin mass dm 0 = ρ 0 dV 0 -and the active (macroscopically Fig. 2: Active species conformational changes induce the remodeling of the lipid membrane where rafts are formed.This process is modeled through the theory of Structured Deformations [71][72][73][74][75], a multiscale geometric framework that allows for tracing back sub-macroscopic changes in combination with classical macroscopic deformation between the active reference and the current deformed state.In the model, an inactive (undeformed) configuration is first mapped onto a geometrically identical configuration in which transmembrane proteins pass to their active state, this being characterized by the conformational jacobian K r (standing for the change in volume induced by disarrangements that are here caused by the submacroscopic remodeling).Material points in the active (reference) configuration are then mapped onto the current (deformed) one by means of the pair (x, F) representing the classical motion/deformation path.Here F = ∇y (X), and x = y (X), where X is a material point in the active configuration and y represents the macroscopic deformation of the body.undeformed) state, where the active mass of the material points instead read as dm a = ρ a dV a (see Fig. 2).Conservation of mass at the local level leads to K r = dV a /dV 0 = ρ 0 /ρ a , with the densities ρ (k) in the heterogeneous medium being calculated as the sum of the true densities of lipids and proteins weighted by the respective fractions (see e.g.[38]).With this in mind, thermodynamical principles allow for expressing the chemical potential as: where, by virtue of ( 5) and ( 6), the species' chemical potentials µ i write as follows: (9) On the other hand, in deriving the mechanical stresses, the Clausius-Duhem inequality leads to: with S * denoting the second Piola-Kirchhoff stress tensor with respect to the virgin configuration.
In the present consitutively enriched model, a viscous dissipation potential W v (C, Ċ) is introduced to take explicitly into account the energy dissipation due to the inherent viscosity of the membrane medium that, in the case under exam, is a pure lipid system.In this way we exclude more complex mixtures involving other structural macro-molecules such as cholesterol, whose presence in different percentages affects the membrane properties.Under these assumptions, the non-negative condition (10) equates the internal dissipation such that [77,78]: or This can be expressed also in terms of the Cauchy stress through a standard push-forward operation from the reference (active) to the current configuration.By considering volumetric incompressibility, one obtains: where the right-hand side of (2) has been considered.Therefore, visco-elasticity of the membrane will depend on the specific choice of the dissipation potential.As aforementioned, the plasma membrane behaves as a visco-elastic material that experiences a vast variety of physical states with both liquid-like and solid-like behaviors [9].For these reasons, viscous components could be included in a straightforward manner in order to account for such a liquid-solid description [79].
Herein, the stress-strain relation ( 13) can be particularized through a Kelvin-Voigt-type nonlinear viscous term proportional to the rate of deformation, in order to account for rapid system variations.The Kelvin body does indeed return to its original configuration when the load, or more in general the source of deformation, is released, as typical of visco-elastic bodies [80].To this extent, it is possible to study the interplay between the characteristic relaxation time of the membrane and the protein activation dynamics in order to capture differences in lipid rafts behavior.Under these assumptions, the Cauchy stress tensor, with respect to the current configuration, reads as follows (see e.g.[81][82][83]): The viscous part of the stress is thus defined through the viscosity term η > 0, which can be either constant as in the case of linear viscoelasticity or can be a function of polynomial scalar invariants involving the strain and the strain rate tensors [77,78,82].In what follows, we will focus on the effects of both possible constant viscosities as well as a strain-sensitive viscosity.In the light of this, it is worth highlighting that the particular constitutive choice in (14) corresponds to considering a dissipation potential of the type: where the right-hand side of (2) has been used (the pulled-back fourth order identity tensor is defined such that [A⊗B] ijhk = A ih B jk ).In addition, by considering the free energy of the system (4) involving the coupled potential ( 6) and the neo-Hookean strain energy contribution (7) of the membrane, the Cauchy stress assumes the following expression: (16) Under the assumption of plane stress, the outof-plane stress component σ 33 = e 3 • σ • e 3 vanishes thus leading to estimate the pressure p.By restricting the deformation gradient in the mid-plane of the membrane, one has that: This allows to obtain the in-plane Cauchy stress σ 0 as follows: in which I 0 and D 0 are respectively the in-plane identity operator and the strain rate.In order to write equilibrium with respect to the reference domain, the in-plane nominal stress tensor can be obtained through a Piola transformation as P 0 = σ 0 F −T 0 , so having: where the relation φ = −ϕ( Ḟ0 : F −1 0 ) is employed because of incompressibility.Consequently, the pulled-back stress reads as follows: By neglecting body forces and inertia terms, the mechanical equilibrium of the membrane reads: with ∇ 0 representing the in-plane nabla operator in the virgin configuration.
As said, the mechanical stress terms involve the co-action of resident transmembrane protein species, whose dynamics induce the rearrangement of the membrane and, in turn, its overall deformation.Therefore, the coupled system at hand must provide the presence of species-related mass balances.The generic mass balance equations for the i-th species ṅi , given in terms of the species' reference flux Q i and the interspecific rates Γ i , are thus calculated according to the above attained chemical potential: The flux term Q i = −L i ∇µ * i refers to the driving force ∇µ * i generating species momentum in the mass balance and mediated by the scalar diffusion mobility parameter L i .While, the source term Γ i measures chemical interactions between the two protein populations, namely GPCRs and MRPs indicated with ξ and ζ respectively.Given their mutual interaction extensively explained in [40], through Volterra-Lotka-like interspecific terms, the mass conservation equations write: where such dynamics is regulated by the decay rates δ i , the interspecific terms β ij and the activation term α ξ that regulates the activity of GPCRs.More specifically, the uptake function α ξ accounts for the response of the receptor to the ligand precipitation rate whose kinetics is controlled in time by a generic Gamma distribution γ (t) and spatially by a distribution function ι(x).Therefore, one can write where k b is defined as the binding constant, and Q is the total quantity of ligand averaged over the membrane area [40].
All the values adopted for the numerical study are reported in Table 1.

Governing equations of the model
Given the well-established interplay between GPCRs structural and functional organization of the cell membrane and the bilayer thickness and stress variations [40], we now present the governing equations regulating the modeled dynamics.In this sense, the mechano-biological process turns out to be governed by the balance of linear momentum in (21) and the time-evolution laws in (23) for the two protein fractions GPCRs and MRPs involved in the ligand-binding.Indeed, these species have been selected as the main families of transmembrane proteins that participate to the regulation of the membrane microenvironment.Therefore, one has the following set of coupled equations: Numerical solutions of such system have been implemented in the software COMSOL Multiphysics ® [93], by adopting a monolithic scheme of fully coupled PDEs solved simultaneously by using a Newton nonlinear method and by discretizing the domain through a Delaunay tessellation.This by considering a circular domain Ω = {(x, y) ∈ R 2 : x 2 + y 2 ≤ R 2 } with R = 5µm, and a time span t ∈ [0, t max ], where t max = 1h [40].Provided constant initial conditions for the protein fractions ζ (x, y, 0) = ζ 0 and ξ (x, y, 0) = ξ 0 , the in-plane displacements are both set with null initial values u (x, y, 0) = 0. Also, null species fluxes imply the boundary condition ∇n i • N = 0 for the proteins and a stress-prescribed situation with a non-zero radial stress at the boundary is considered to simulate the Laplace membrane tension due to the intracellular pressure.Therefore, the nominal traction in the radial direction at the outer radius writes P * 0 • N = T R N, which can be evaluated through a prescribed outer (actual) pressure p o by imposing the equivalence p o h ds = T R h 0 dS 0 that leads to T R = p o (1 + u R /R)/J 0 , where u R stands for the magnitude of the in-plane displacement at the boundary.In the following section, we will show the influence of viscous dissipation on the solid-liquid behavior of plasma membranes under different conditions able to reproduce scenarios in which membrane's morphology and mechanical adaptation lead to various situations.

Results and discussion
Within the framework of membrane viscoelasticity, we here present numerical results that permit to observe the viscosity landscape of the  phase-separated domains, by focusing on possible differences in terms of raft lifespan and heterogeneity.To this aim, sensitivity analyses will be carried out to map the evolution of an initially (geometrically and materially) homogeneous membrane, by observing how raft domains and viscosity change.This will be mainly investigated as a function of the membrane's (elastic and viscous) tangent properties and initial protein distributions.In the light of the pivotal role of mechanics in the spatio-temporal dynamics of the raft-associated proteins, we analyze proteininduced adaption processes.Indeed, conformational changes of GPCR and MRP populations are capable to induce the overall remodeling of the bilayer at the membrane scale.With this in mind, in order to trigger the activation dynamics, we consider the realistic situation in which extracellular molecules randomly precipitate on the domain.This is done by assigning a random distributions to the ligand precipitation rate functions used in (23) and by modulating the amount of precipitating ligand to induce differential receptor responses, thus orienting the membrane dynamics towards various patterns.
In numerical analyses, we start from studying the effects of a constant viscosity on the spatiotemporal behavior of the ordered phase.To then investigate more in depth the material adaptation of the bilayer in terms of the evolution of viscous properties of the rafts through a strainsensitive viscosity term.This enrichment allows to follow the strain-induced remodeling of the lipid phase.In particular, this is done by meeting wide literature evidences demonstrating that viscosity of ordered clusters tends to increase as the phase order increases [94].Starting from the initial Newtonian hypothesis, sensitivity analyses are carried out by varying the viscosity over a range compatible with literature data.In this respect, surface shear viscosity seems to exhibit a large variability depending on the particular composition of the mixed lipid system, on the specific conditions in which tests are performed as well as on the adopted experimental methods.Typical values of tangent viscosity for the most of biological membranes result of the order of 10 −3 − 10 2 P a.s [5,9,10,89,95].Fewer cases were found to instead exhibit significantly higher tangent viscosities ranges of 10 5 − 10 6 P a.s [9,89], up to peaking to unusual values 10 9 P a.s in case of the so-called tough visco-elastic systems [91,92].However, it is worth highlighting that these experimental observations report significant differences when cholesterol is introduced in the mixed lipid systems.In particular, cholesterol highly affects the stiffening and the viscosity increase of the membranes and it has a direct impact on raft stabilization as well [89,96,97].In the present model, we limit our analyses to pure and mixed lipid systems, for now excluding the explicit modeling of cholesterol as a structural component of the membrane medium, which could be instead taken into account through the suitable determination of homogenized material properties depending on the extent of cholesterol fraction.

Insights on the influence of tangent stiffness and viscosity on membrane remodeling from a Newtonian model
First, we assume the simplest case with a constant viscosity term η, whose range of variability is reported in Table 1.This is considered as a mean shear viscosity, evaluated on the whole membrane, that does not take into account the fluidic variation in phase transitions.When η is a constant, given the wide range of viscosity values, outcomes have been organized and presented by referring to two classes of visco-elastic responses, denoted as the weak and the tough visco-elastic systems.The former case indicates Newtonian viscosities lying in the wide range 10 −3 −10 5 P a.s, which characterizes most of the biological membranes encountered throughout the literature.Their behavior varies from that one of a low viscosity fluid to that of a visco-elastic gel.In such a situation, linear visco-elasticity results to minimally interfere with the chemo-mechanical activity of the membrane and the overall dynamics almost entirely proteindominated.The most important differences are indeed appraised by varying the initial stiffness of the membrane, which really does affect the coupling.The tangent Young's modulus is assumed to vary so that the membrane can undergo different configurations in the solid-fluid transition.Indeed, the stiffness of the environment mediates the mechanical work performed by proteins on the lipid medium.By considering as representative, and most frequent, cases for the weak visco-elastic systems the values η = η 1 = 100P a.s and η = η 2 = 10 −3 P a.s, Fig. 3A shows that the thickness stretch is mostly determined by variations in the elastic part rather than the dissipative one.It indeed increases at higher Young's moduli, though it does not significantly change when different viscosity values are employed.Coherently with literature findings [98], the out-of-plane deformation results to be in a range of about 20 − 50%.It is worth to note that the coupling parameters w i vary proportionally with the elastic modulus by so influencing the overall membrane activity and deformability.In fact, as such coefficient translates the microscopic mechanical interaction at the protein subunit-membrane interface, it results to be proportional to the local surface tension.That inevitably involves the stiffness of the lipid medium [38].Moreover, for the higher viscosity η 1 = 100P a.s, the influence of the elastic part results in both the activation time of the raftassociated proteins GPCRs and the persistence of L o phase in the bilayer (see Fig. 3B).As shown, in the case of a more deformable system, the receptor-ligand biding occurs at t ≃ 430s accompanied by a faster raft duration of about 10s.Stiffer membranes instead produce a slower response of GPCRs, although a larger duration of the L o domain up to a lifespan of 100s is ensured.Noteworthy, these delays in the activation times of Fig. 3B can be produced by the competition of the viscosity with the internal protein dynamics.The latter emerges from the complex interplay of protein intrinsic rates and stiffness-associated work terms influencing their spatio-temporal evolution through the species' momentum terms.
The low influence of Newtonian viscosity de facto suggests to adopt nonlinear viscosity models.To get more insights into the influence that a constant viscosity term can have on membrane dynamics, we carried out -at least as illustrative theoretical cases-simulations that take in consideration the extreme situation of tough viscoelastic membranes.This is reported to the best of Authors' knowledge in few literature works concerning the characterization of red blood cells' membranes [91,92].By thus prescribing steep values of viscosity capable to interfere with membrane dynamics, it is possible to observe a drastic change of the bilayer's morphological response to the activation of protein populations.Indeed, as shown in Fig. 3C, GPCRs evolve in a substantially analogous manner both in the weak and tough Fig. 3: Lipid membrane response to elastic and dissipative variations.A: Thickness stretch ϕ measured at constant viscosities with varying Young's modulus.Viscosity variation does not significantly affect the out-of-plane deformation that is instead influenced by changing in membrane rigidity.B: At fixed η = 100P a.s, membrane undergoing deformability and rigidity results in changing the activity of GPCRs and the formation of rafts domains.C: Influence of weak and tough viscosities on the morphological reorganization of the membrane in response to analogous GPCRs activity.D: Thickness stretch and raft domains persistance measured for weak and tough visco-elastic systems.Highly viscous system leads to variations in membrane remodeling.visco-elastic cases, since they respond to the same imposed chemical stimulus.On the other hand, in the fluid case, after the initial contraction due to the applied tension, membrane thickening grows with strong synergy and has a reduced relaxation delay following the GPCRs' decay.Conversely, in the tough system, raft emergence forms with much slower velocity.There, the extremely viscous environment highly reduces the proteins' mobility, by preventing their capability to exert mechanical work against the membrane, and by also inducing high retardation in the morphological adaptation of the plasma medium to receptors' desensitization.This is confirmed in Fig. 3D at different viscosities.In the fluid-gel regime, dynamics leads to co-localized and almost synchronous progression with similar morphological rearrangement, this drastically decelerating in tough visco-elastic systems with a consequent decline of the out-ofplane reconfiguration.In the light of these considerations, the latter cases demonstrate that high initial viscosity contrasts the highly dynamic and heterogeneous character of plasma membranes, by compromising the co-evolution capability.That allows the bilayer to exhibit a sufficiently reactive Fig. 4: Surface plots showing the active GPCRs domains in the visco-elastic system with fixed η = 100P a.s and varying elastic moduli.Such a variation influences membrane remodeling and configuration.It is indeed evident that a more rigid surface leads the rafts islands to be more persistent in time by reducing the lateral mobility of transmembrane proteins.morphological adaptation able to favor the formation of ordered domain working as necessary sights for chemical signaling.
Then, with reference to more common viscoelastic gel-like systems (at η 1 = 100P a.s), differences in durability can be captured in terms of prolonged protein activity in stiffer environments.In fact, as reported in Fig. 4, variations in the persistence of receptor ligand-binding reflect the spatial organization of the bilayer in terms of raft emergence and membrane relaxation.Although the maximum activity of GPCRs occurs at slightly different times, as observable starting from t ≃ 400s, the thickened L o domains decay faster in the softer membranes -being they almost extincted already at 800s-while the formed GPCRs clusters are still active in membranes with a higher degree mechanical interaction.

Effects of strain-sensitive viscosity and evolution of membrane fluidity
Further information can be envisaged by introducing a more complex viscous term in the model.Indeed, nonlinear effects could occur during moderate-to-large strains.In turn, this could involve non-Newtonian responses for the shear viscosity.In this way, it is possible to capture the effective fluidity of the membrane upon large strength motions [9].For this reason, a strainlevel dependent viscosity is assumed in a purely phenomenological fashion.This allows us to investigate situations able to theoretically confirm that the viscosity depends on membrane composition, thus it varies following ordered-disordered phase transition [94].
To this aim, among the possible constitutive choices and in order to introduce an essential functional variability (see e.g.[77,78,82]), we assume that the viscosity term is a function of the right Cauchy-Green strain tensor through its first invariant.This is done here by means of the expression η m = η 0 [1 + τ 0 (tr(C) − 3)].Herein, the tangent (Newtonian) viscosity η 0 has been set equal to η 1 , being it compatible with the order of magnitude of the most of lipid systems.Furthermore, the coefficient τ 0 is a nondimensional parameter modulating the sensitivity to the strain.In order to determine a proper value of this latter coefficient, we exploited data in Kelley et al. [99], reporting experiments and associated scaling relationships for the viscosity of mixed lipid membranes as a function of the lipid area per unit molecule.In particular, as also shown in Fig. 5A the lower is the available area per lipid the higher results the viscous term.In the present continuum approach, the area per unit lipid molecule can be put in direct correlation with the in-plane areal stretch J 0 .To this end, by assuming a homogeneous deformation, one can fit experimental points to calibrate the proposed strain-dependent viscosity law, so deriving a reference value for the fitting parameter τ 0 (τ 0 = 17.35).However, in order to account for the large variability of membrane fluidic properties and investigate the influence of strain sensitivity, possible variations of the parameter τ 0 have been prescribed during the numerical simulations (three values proportional to τ 0 have been assumed).The proposed phenomenological law for the viscosity proposed above has been then uploaded in the coupled model in order to analyze the evolution of raft viscosity during membrane activity.In particular, the effective viscosity of raft domains has been evaluated as the tangent viscosity at the achieved strain level as η raf t = A −1 raf t A f (ϕ) η 0 K r [1 + τ 0 (tr(C) − 3)] dA, with the auxiliary function f defined to select raft zones as f (ϕ) = (1 + tanh(χ ϕ − ϕ ), while the raft area coverage results A raf t = A f (ϕ) dA (see the Appendix for details on tangent viscosity).As it can be noticed in Fig. 5B, the numerical simulations show that raft viscosity intensifies from four up to ten times at the moment of maximum activity, depending on the strength of strain sensitivity.These increments are consistent with many experimental works reporting that L o phases exhibit a higher viscosity than the L d domains [5,89,94,99,100].Thus, this approach suggests that the adopted nonlinear viscosity can represent a proper strategy to predict the dynamic changes of membrane fluidity during order transitions.Fig. 5: Fitting parameter τ 0 .A: Determination of the viscosity sensitivity to membrane strain.Data adopted from [99].B: Analysis of strain-induced viscosity, at maximum protein activity, for different strain sensitivity values τ .
Noteworthy, the strain-dependent membrane shear viscosity can be affected by the intra-cellular tension that acts on the bilayer in both structural and dynamical properties [101].Therefore, we performed simulations with different pressures p 0 at the stress-prescribed boundary.Outcomes are shown in Fig. 6 where, according to literature findings [102], the membrane tension ranges from Fig. 6: Membrane mechanical properties evaluated at different membrane tensions.The viscosity of the ϕ L0 domain decreases as the pressure p 0 increases in the range of 0 − 1.2M P a, as well as membrane thickening, suggesting that such mechanical properties varies with the intracellular stimuli.0.1M P a to 1.2M P a.Such values are consistent with the levels of intracellular pressures (Laplace's law implies that p 0 ∝ p cell × R cell /2h 0 ≃ 10 3 p cell , being the intracellular pressure of the order of 0.01 − 1 kPa [103]) and keep below the estimated rupture tension of 2M P a [104].From Fig. 6 one can also show that, at fixed τ = τ 0 , the effective raft viscosity η raf t /η 0 tends to decrease as the intra-cellular pressure increases.Such behavior is reasonable with the established relationship between membrane tension and bilayer mechanical response [101,105].Indeed, increasing pressure reduces membrane thickness and works for areal expansion.It competes against the morpho-taxis phenomena involving membrane thickening and contrasting the tendency of transmembrane proteins to aggregate, thereby reducing the ligandbinding effectiveness and resulting in lower L o volume fraction.
It is then apparent that membrane shear viscosity varies with lipid phase order.This is due to the fact that ordered-phase islands exhibit a higher level of lipid packing compared to L d domains, by so resulting to be less polar and more viscous [106].In particular, according to literature measurements, the L o regions seem to be characterized by a membrane viscosity higher than the one of the L d phase [5,[107][108][109].
To appraise these differences, we studied the viscosity behavior as a function of the volume fraction of the disordered phase ϕ L d .This was done numerically by varying the amount of precipitating ligand, by so influencing the activation potential of the transmembrane proteins.As analyzed in Fig. 7, the theoretical curve shows a two-fold viscosity ratio passing from a predominantly disordered phase to a domain mostly occupied by ordered clusters.These numerical outcomes have been put in direct comparison with two different sets of experimental data available in the literature.First, Sakuma et al. [94] correlated the order parameter with the measured viscosity for different lipid systems.In such a case, the relative viscosity variations obtained from theoretical predictions well fit with these literature findings in the range 0.5 ≤ ϕ L d < 1.0.Below such an interval, i.e. for 0 < ϕ L d ≤ 0.5, the here presented model is far from capturing the experimental data obtained in Sakuma et al., as the reported values refer to lipid mixtures in which ordered and disordered phases coexist with a high cholesterol percentage.It is indeed confirmed that significant cholesterol percentages increase membrane viscosity [97,110] and can impact on the change of membrane properties by chemically altering the lipid micro-environment.In the case at hand, for 0 < ϕ L d ≤ 0.5, these bilayers turn out to be rich in cholesterol content (about the 30% more than the average ones) produced a different trend.In this sense, the lack of such species in the system represents a limitation, and more faithful results could be achieved by introducing a finer description of its role in the multi-physics model.More interestingly, the increase in viscosity predicted in silico results that are remarkably compatible with additional literature findings over the entire range of phase order.In fact, the numerical curve is found to be in excellent agreement with data points derived from the experiments performed on giant unilamellar vesicles (GUVs) performed by Wu et al. [100], in which lower Chol concentrations were employed.Noteworthy, they obtained a more gradual change of viscosity variation that increases to 2.1 for ordered membrane configurations, so demonstrating the dynamic change of viscosity involved also in lipid rafts.Fig. 7: Numerical measured viscosities compared with experimental data adapted from Sakuma et al. [94] and Wu et al. [100].By assigning different spatial distributions in the ligand precipitation rate, in order to modulate the volume fraction of disordered domains, the model is capable to find consistent values with both the experimental findings in the range 0.5 ≤ ϕ L d < 1.0.Cholesterol rich membranes, 0 < ϕ L d ≤ 0.5, lead to variation in the measured viscosities that differ from the ones measured in absence of cholesterol percentages and the ones numerically found.Surface plots of disordered phase volume fractions are shown above and viscosity maps are visible on the right (adopted parameters p 0 = 0.8M P a and τ = τ 0 ).

Conclusions
Following a recent theoretical formulation describing the mechanobiology of lipid membrane remodeling and raft formation carried out in [38,40], the current study aims at investigating the dynamic visco-elastic response of plasma membranes to chemo-mechanical stimuli.Through in silico analyses accounting for viscous-associated terms in the constitutive model, the multiphysics coupling between chemical events and mechanical adaptation highlights how the solid-fluid behavior of the bilayer evolves with the activity of the membrane.The evolved processes are strongly influenced by the dynamics of the transmembrane proteins activation and their interaction with the lipid medium.By considering both the cases of a Newtonian shear viscosity and a strain-sensitive viscosity, in this present paper we investigate the relationship between the reconfiguration of an initially inactive membrane micro-environment as a function of the competition between the internal viscous dissipation and the kinetics of phase transitions governing the emergence of lipid islands.
Numerical outcomes allowed one to observe that the shear viscosity varies in phase-separated membranes resulting in higher values for orderedphase domains, i.e. lipid rafts.Hence, this provides a mechanically-based explanation of a well-known phenomenon highlighted by a large number of biophysical studies by means of various experimental methods.The synergy between active protein regions and raft emergence leads the system to re-organize itself by creating thicker and more viscous domains.Also, sensitivity analyses revealed how the visco-elastic behavior is influenced by the intra-cellular pressure applied at the boundary.That alters the mechanical properties of the membrane, and the volume fraction of the liquid-disordered phase.Hence, our visco-elastic approach enriches the existing studies regulating the mechanisms on the lipid membrane's behavior.This could help to earn some insights in characterizing the role of lipid rafts in membrane mechanics and in mediating important cellular biochemical processes.
By refining the modeling of species interspecificity, one would have the opportunity to include some other agents influencing membrane dynamics in the analysis.This may allow one to enlarge the complex multi-species environment under exam, as well as to further enrich the membrane constitutive framework.To this aim, the self-reconfiguration of lipids could be studied by considering non-convex terms in the elastic strain energy (see e.g.[34,38] and reference cited therein).Moreover, enriched coupling terms may be considered in the model in order to have deeper insights into the influence of the mechanical stress on the interspecific dynamics.In fact, through ad hoc mechanical feedback functions, it would be possible to better investigate the processes of cell mechano-sensing and mechano-trasduction, that inevitably involve the mediation of membrane selectivity during cell-environment communication.Also, as emerged from the presented analyses, one of the main components that can be included to further refine and enrich the description of membrane visco-elastic adaptation could be the cholesterol.This has a direct responsibility for lipid rafts stabilization and bilayer lateral diffusion, GPCRs re-configuration and activity, besides its participation to determine the membrane effective properties.For this significant reason, this will be object of future investigations.

where [A⊗B]
. By focusing on the response to the incremental strainrates, the tangent viscosity tensor can be evaluated as follows: ) where S = (I⊗I)/2 is the identity fourth-order tensor mapping symmetric tensors.By virtue of (4) and (A.2), and on account of constitutive expressions ( 16) and (17), after some passages one has: To measure the effective surface shear viscosity, a planar shear velocity v = v 1 e 1 + v 2 e 2 is imagined to be applied on a generic point of the upper membrane surface, by producing a shear deformation γs such that dv = γs dx 3 , or dv 1 = ( γs dx 3 ) cos θ s and dv 2 = ( γs dx 3 ) sin θ s .Then, the corresponding strain rates are linked to the shear γs throught the relations:

Fig. 1 :
Fig. 1: Qualitative influence of the membrane stretch ϕ on the equilibria of the double well coupled potential when a generic homogeneous density fractions is considered, i.e.W ni = W (n i , 0, ϕ)

Table 1 :
Summary of the numerical values for the coefficients used in the model.