Cytoskeleton mediated effective elastic properties of model red blood cell membranes
aa r X i v : . [ q - b i o . S C ] J un Cytoskeleton mediated effective elastic properties of model redblood cell membranes
Rui Zhang and Frank L. H. Brown
Department of Chemistry and Biochemistry and Department of Physics,University of California, Santa Barbara, California (Dated: November 20, 2018)
Abstract
The plasma membrane of human red blood cells consists of a lipid bilayer attached to a regularnetwork of underlying cytoskeletal polymers. We model this system at a dynamic coarse-grainedlevel, treating the bilayer as an elastic sheet and the cytoskeletal network as a series of phantomentropic springs. In contrast to prior simulation efforts, we explicitly account for dynamics of thecytoskeletal network, both via motion of the protein anchors that attach the cytoskeleton to thebilayer and through breaking and reconnection of individual cytoskeletal filaments. Simulationresults are explained in the context of a simple mean-field percolation model and comparison ismade to experimental measurements of red blood cell fluctuation amplitudes. . INTRODUCTION Membranes are essential components of all biological cells [1]. In addition to their biolog-ical importance, lipid bilayers and biomembranes have also attracted considerable attentionfrom physicists due to their fascinating and unusual properties [2, 3]. One particularly wellstudied system is the human red blood cell (RBC) membrane. The RBC membrane is acomposite structure, consisting of a lipid bilayer adhered to an underlying network of fila-mentous cytoskeletal proteins (spectrin) via integral membrane protein anchors (see Fig. 1).The spectrin network is observed to be quite regular [4], with an approximate hexagonalsymmetry extending over the entire cell surface. (It is worth emphasizing that the spectrinbased cytoskeletal network of RBCs is completely different from the three dimensional actinnetworks common to other types of animal cells [1].)Given the apparent simplicity of the RBC membrane, it is tempting to attempt modelingwith elementary elastic models. Indeed, there is a rich history of elastic RBC models to befound in the literature [5, 6, 7, 8, 9, 10]. Without attempting a full historical review here, wecomment that no single elastic model has yet been identified that is capable of reproducingthe full set of experimental data available for RBC membranes. For example, while thework of Lim, Wortis and Mukhopadhyay [8] is capable of capturing the range of observedRBC shapes (stomatocyte, discocyte, echinocyte and non-main-sequence shapes as well) seenunder various chemically induced stresses (e.g. pH, salt, ATP, etc.), this model has not beenapplied to explain the thermal fluctuation amplitudes observed in RBC membranes. And,while the models of Gov, Zilman and Safran [7] and Fournier, Lacoste and Raphael [9] appearto do a good job fitting thermal fluctuation data [11], these models do not appear capableof explaining mechanical deformation experiments on RBCs [12, 13, 14]. The failure ofsimple models to consistently explain both thermal fluctuations and mechanical deformationexperiments has been recognized for years [6]. One recent model does explain both types ofdata within a single elastic model [10], however this model treats the spectrin network as anincompressible and homogeneous viscoelastic plate coupled to the lipid bilayer. It remainsunclear why such an approximation should suffice for the sparse cytoskeletal network presentin RBCs. Additionally, this model seems too simplistic to capture the full range of shapebehaviors explained in reference [8].The models mentioned in the preceding paragraph make no mention of the role of energy2xpenditure in the behavior of RBC membranes. This despite the fact that it is known thatRBC membranes possess kinase and phosphatase activities capable of altering the proper-ties of spectrin and other network associated proteins via (de)phosphorylation [15, 16]. And,certain measurements of RBC fluctuation amplitudes show a correlation between ATP con-centration and fluctuation magnitude [17, 18]. One might argue that the difficulty in fittingall RBC behavior to a single elastic model stems from the fact that a truly comprehensivemodel must incorporate the effects of energy expenditure by the cell in a realistic fashion.Gov and Safran [19] are the first to seriously consider active energy expenditure within theRBC from a theoretical standpoint. They have proposed that ATP induced phosphoryla-tion and dephosphorylation of the RBC cytoskeletal network leads to a continual dynamicevolution of the integrity of the spectrin network. While this picture remains hypothetical,without direct proof, it is consistent with the general observations relating ATP concentra-tion to membrane fluctuation amplitudes. Gov and Safran (GS) [19] have used this pictureto motivate a simple picture for RBC fluctuations under the presence of ATP. Local breakingand reforming of the spectrin network is captured via non-thermal forces imparted on anelastic membrane model. This model has provided the first plausible explanation for theviscosity dependence of RBC fluctuations [20].While elastic models with proper accounting for non-thermal energetics may eventuallyprove adequate in describing the long wavelength physics of the RBC, it is clear that wave-lengths near or below the spectrin mesh size ( ∼ nm ) must be considered within a moremicroscopic picture. A recent model by Dubus and Fournier (DF) [21] has extended thetraditional elastic modeling of RBC membranes to explicitly include the cytoskeletal net-work at a molecular level of detail. Within this model, the spectrin network is consideredas a completely regular hexagonal network of phantom entropic springs attached to a fluidlipid bilayer. Over wavelengths significantly longer that the spectrin mesh size, the networkso modeled becomes mathematically equivalent to an imposed surface tension on the fluidbilayer. At wavelengths comparable to and shorter than network spacing, the system be-haves differently from a simple membrane with applied tension. This model was used tocompute the spectrum of thermal fluctuation amplitudes for the RBC membrane, but madeno attempt to account for non-thermal consumption of energy and only computed thermal(non-dynamic) observables.In this paper, we extend the DF entropic spring model of the cytoskeleton meshwork to3nclude dynamic evolution of the entire system. We allow the anchor points between spectrinand membrane to laterally diffuse and we allow for dynamic dissociation and associationof spectrin links as a molecular level manifestation of the non-thermal energetic pictureproposed by GS (fig. 2).One important consequence of the GS picture is that sufficiently high ATP concentra-tions lead to an appreciable fraction of dissociated spectrin links at the membrane surface.Depending upon the timescales for spectrin (re)association, the effective long-wavelengthelastic properties of the membrane interpolate between two limiting cases. If spectrin(re)association kinetics are much slower than all other timescales in the problem, the ef-fective tension imposed by the network (as inferred by out-of-plane bilayer undulations) iswell predicted by a simple percolation-theory argument (see fig. 3). In the opposite limitof fast spectrin kinetics, the effective tension is well predicted by assuming each link in thenetwork has a reduced spring constant proportional to the probability of the link being in-tact at steady state. At intermediate rates, simulations are seen to interpolate between thetwo extreme cases.This paper is organized as follows. In Sec. II, we present our mathematical model forthe RBC membrane. In Sec. III details of our simulation methods are discussed. In Sec. IVresults for a fully intact spectrin meshwork are presented, while in Sec. V we generalize to themore interesting case of a randomly broken network (both static and dynamically broken).We discuss our results in relation to experiment in Sec. VI and conclude in Sec. VII. II. MODEL
We treat the RBC membrane as a Helfrich fluid sheet [22] coupled via mobile anchorpoints to a network of springs. Our Hamiltonian is H = Z A d r κ ∇ h ( r )] + X h i,j i ξ ij ( t ) n µ r i − r j ) + µ h ( r i ) − h ( r j )] + E o . (1)The first term (integral portion) is the standard bilayer bending energy for a Monge gaugesheet assuming small deformations [22] and bending modulus κ . A is the projected area4f the membrane, r = ( x, y ) is the position vector in the xy plane and and h ( r ) is thelocal displacement of the membrane away from the flat reference configuration specified by h ( r ) = 0 (see Fig. 2). We assume that the lipid bilayer itself has a negligible (bare) surfacetension. In more general situations, eq. 1 is easily modified to account for non-vanishingtension inherent to the lipid bilayer portion of the membrane [2, 22]. The second term (sumportion) accounts for the energy of the 2D cytoskeletal meshwork, modelled as a network ofentropic springs (or in other words, ideal chains of polymers) with effective spring constant µ ( µ = 3 k B T /ℓ k ℓ c with ℓ k and ℓ c the Kuhn and contour length of the polymer, respectively).In a fully intact cytoskeletal meshwork, all spring end points (nodes) are restricted to lieon the surface of the bilayer, so their coordinates are specified by r i and h ( r i ), with indices i (or j ) labelling different nodes. The sum is over all distinct nearest neighbor node pairs(equivalently, over all polymer springs), denoted as h i, j i . The factor ξ ij ( t ) is included toaccount for the possibly incomplete connectivity of the network. It is equal to 1 when the ij link is connected and 0 otherwise. The constant E reflects all other free energy changeassociated with connecting a detached filament end to a node besides the elastic energyof the spring (e.g., the binding energy to the node); we assume this energy is negativeand significantly larger than thermal energy scales to insure stability of the network in theabsence of non-thermal energy sources.Although our starting point is very similar to the DF model, we emphasize a few keydifferences. In DF, the lateral positions of nodes are held fixed in the geometry of minimumenergy for a flat bilayer surface; i.e. the r i variables are treated as set constants in DF,not as variables capable of influencing the energetics and/or dynamics of the system. Thisapproximation renders the Hamiltonian analytically tractable, however it is not immediatelyclear that such a choice fully captures all relevant physics in this system. For example, withfixed node positions equally spaced on a regular lattice, the “spring” contribution to eq.1 amounts to a finite differenced version of the usual surface tension contribution to theHelfrich Hamiltonian. At long wavelengths, eq. 1 with fixed nodes is guaranteed to behaveas a Helfrich sheet under tension. If the nodes are mobile, as physically expected, the springnetwork represents a simple manifestation a tethered membrane. Such membranes are knownto exhibit more complicated fluctuations than expected for Helfrich fluid bilayers [23]. Wealso emphasize that much of the following work is concerned with membrane dynamics,which was not considered in DF. One of the most interesting aspects of our model is the5ynamic breaking and reformation of cytoskeletal filaments, which could not be studied withthe equilibrium approach adopted by DF.Dynamics in our system are overdamped, owing to the low Reynolds number environmentpresent at cellular length scales [24]. For bilayer height fluctuations, we have the followingLangevin type equation of motion, which accounts for hydrodynamic flow in the surroundingcytoplasm [25, 26, 27] ∂h ( r , t ) ∂t = Z ∞−∞ d r ′ Λ( r − r ′ )[ F ( r ′ , t ) + ζ ( r ′ , t )] . (2)Here, Λ( r ) is the diagonal part of the Oseen tensor [28], given byΛ( r ) = 18 πη | r | , (3)where η is the viscosity of the surrounding fluid. The above integral is taken over the entire x, y plane; it is assumed that the area of interest, A , is embedded within a periodicallyrepeating environment of identical subsystems. F ( r , t ) is the force per unit area on themembrane, F ( r , t ) = − δHδh ( r , t ) = − κ ∇ h ( r , t ) − X h j i µξ ij ( t ) δ ( r − r i )[ h ( r , t ) − h ( r j , t )] , (4)where δ ( r ) is the Dirac delta function, and the sum is over all nearest neighbors of node i ,denoted as h j i . ζ ( r , t ) is a spatially local Gaussian white noise satisfying the fluctuation-dissipation relation h ζ ( r , t ) i = 0 , (5) h ζ ( r , t ) ζ ( r ′ , t ′ ) i = 2 k B T Λ − ( r − r ′ ) δ ( t − t ′ ) , (6)where k B is Boltzmann’s constant and T temperature of the system.We have another set of Langevin equations describing lateral diffusion of the nodes withinthe bilayer ∂ r i ( t ) ∂t = Dk B T [ F i ( t ) + ζ i ( t )] , (7)with F i ( t ) = − ∂H∂ r i = − X h j i µξ ij ( t ) (cid:26) ( r i − r j ) + [ h ( r i ) − h ( r j )] ∂h ( r i ) ∂ r i (cid:27) , (8)6nd h ζ i ( t ) i = 0 , (9) h ζ i ( t ) · ζ j ( t ′ ) i = 4( k B T ) D δ ij δ ( t − t ′ ) , (10)where D is the lateral diffusion constant of the node across the membrane surface. Eqs. 2 and7 completely specify the dynamics of the lipid bilayer and the attached 2D meshwork. Noticethat the two sets of Langevin equations are coupled (and must be solved simultaneously)via the shape of the membrane surface. We note that eq. 8 neglects the purely geometriceffect of non-flat membrane geometry on the ( x, y ) motion of node points [29, 30]. Thisapproximation significantly simplifies our modeling and it has recently been demonstratedthat such geometric effects are very small for the physical parameters studied herein [30].It is convenient to recast eq. (2) in a Fourier basis [26]. ∂h k ( t ) ∂t = Λ k [ F k ( t ) + ζ k ( t )] (11)with k = (2 πm/L x , πn/L y ) for integer m and n . Here for the latter convenience of thesimulation, we assume in general a rectangular sample with size L x × L y in real space, andtherefore the two lattice constants in k space are different. The quantities h k , F k and ζ k derive from functions periodic in x and y , due to the assumed periodicity of the system.The Fourier transform pair for an arbitrary function, f , with such periodicity is f ( r ) = 1 L x L y X k f k e i k · r , (12) f k = Z A d r f ( r ) e − i k · r . (13)The Fourier transformed Oseen interaction,Λ k = Z ∞−∞ d r e − i k · r Λ( r ) = 14 ηk , (14)in contrast, reflects transformation over the full 2D plane. By construction, the dynamicsspecified by eq. 11 reflects an infinite network of periodic membrane images interacting viathe long range 1 /r Oseen hydrodynamic kernel. The random forces obey h ζ k ( t ) i = 0 , (15) h ζ k ( t ) ζ k ′ ( t ′ ) i = 2 k B T L x L y Λ − k δ k , − k ′ δ ( t − t ′ ) . (16)7or the dynamics of the breaking and reforming of spectrin springs, we consider a simpletwo state kinetic model. We define the rate to reconnect a link as k c , and the rate todisconnect a link as k d , irrespective of the instantaneous position of the endpoints of thespring. Accordingly, the value of each ξ ij jumps back and forth between 1 and 0. Defining p as the steady-state probability of a link to be connected at any moment, we have p = k c k c + k d . (17)It should be emphasized that our simple model for the breaking and reformation of spec-trin filaments does not obey detailed balance since we do not account for the variationsin energetics caused by positions of the network nodes within the kinetic scheme. Thisapproach necessarily corresponds to a non-equilibrium situation, with the dynamics of thespectrin network driving membrane fluctuations in a non-thermal manner. Qualitatively,this corresponds to the picture proposed by GS, however the detailed kinetics involved inthe spectrin (re)association process are unknown and likely differ substantially from thepicture adopted herein. Our simple two-state picture represents one possible manifestationof non-equilibrium driving.The coupled equations implied by eqs. 11 and 7 are not amenable to analytical solu-tion and will be solved via simulation as detailed below. Before proceeding, we note thatsimulations are run on a discrete square lattice. In other words, Eq. (13) is replaced by f k = X r a f ( r ) e − i k · r , (18)where a is the lattice constant and r now only take discrete values on the lattice. Cor-respondingly, in k space, the reciprocal lattice (with lattice constant 2 π/L x and 2 π/L y )contains ( L x /a ) × ( L y /a ) points and the summation in eq. 12 is finite.Two types of meshworks are considered in our simulations. First, in accordance withthe true geometry of RBC membranes [31], we consider a hexagonal meshwork (Fig. 4a).In such a meshwork, when all links are connected, the average positions of all nodes forma hexagonal lattice (this lattice of nodes should not be confused with the lattice used todiscretize the surface as discussed above). We consider a finite sample and assume periodicboundary conditions in a rectangular geometry of size L x × L y approximately commensuratewith the embedded hexagonal meshwork (see Fig. 4a). The average distance between thenearest neighbor nodes, A , is determined from the average surface density of nodes in a8BC membrane. For theoretical interest, we also consider a square meshwork with squarelattice symmetry for the average positions of nodes (Fig. 4b). Here we simply take a periodicsquare box, i.e., L x = L y = L .To connect with experiment and prior theoretical work, our simulations are used primarilyto calculate the mean square height fluctuation of the membrane surface in k space, h| h k | i (angular brackets represent non-equilibrium averages as well as equilibrium averages in thiswork). At long wavelengths ( λ ≫ A ) we observe that the relation h| h k | i = k B T L x L y κ eff k + σ eff k . (19)holds fairly well. This expression corresponds to that expected [2, 22] for a thermal fluidbilayer sheet with effective bending rigidity κ eff and effective surface tension σ eff , but westress that its use in interpreting the simulations is empirical. The composite membranesurfaces studied in this work can not be regarded as fluid-like due to the assumed connectivityof the cytoskeleton matrix. Much of our analysis is presented in terms of σ eff , which isobtained by fitting the simulated data to eq. 19 (while assuming κ eff = κ unless notedotherwise). For future reference, we define the “free fluctuation spectrum” of the sheet, h| h k | i f , as the result anticipated from eq. 1 in the absence of any cytoskeletal contributions(i.e. all ξ ij = 0) h| h k | i f = k B T L x L y κk . (20) III. SIMULATION METHODS
Two simulation methods were used to study the composite membrane system: FourierMonte Carlo (FMC) [32, 33] and Fourier space Brownian dynamics (FSBD) [26, 34, 35].In the limit of infinitely slow breaking and reformation of spectrin filaments (quencheddisorder), the ξ ij variables are static over the course of any finite simulation and the systemrelaxes to a thermal equilibrium dependent upon the connectivity of the network. In thislimit, both FMC and FSBD simulations may be used to calculate thermal averages andboth should agree with one another. When spectrin links are breaking and reforming intime following the non-equilibrium scheme presented above, we must run dynamic FSBDcalculations. Our primary interest in this work is the non-equilibrium case, however FMCcalculations were performed both to verify the accuracy of our FSBD simulations (in the9tatic network limit) and to examine the scaling of some our equilibrium results with systemsize (the FMC method is computationally more efficient than FSBD). A. Fourier Monte Carlo (FMC)
The FMC scheme [32, 33] is a standard Metropolis algorithm [36], which uses Fouriermodes of the bilayer, h k , as the bilayer degrees of freedom. There are two kinds of degreesof freedom in our system: membrane shape modes and lateral position of the nodes of thespectrin network. For the former, we attempt MC moves on the Fourier modes of the systemto enhance computational efficiency relative to naive real space schemes [32, 33]. While forthe latter, we attempt moves that displace the x, y position of the nodes to adjacent sitesof the real space lattice inherent to the simulations. Note that the shape of the membranesurface is continuously variable in our description, while the lateral position of nodes isdiscrete. The primary advantage to evolving the shape of the bilayer surface in Fourier spaceis that we may tune the maximal size of attempted MC moves for each mode separately. Inthe case of a simple fluid bilayer (without attached cytoskeleton) under finite surface tensionit is clear that a good choice for maximal move sizes is (see Eq. (1) in Ref. [37])( κk + σk ) δ k L x L y = const., (21)where δ k is the maximum attempted jump size of mode k . In our case, a similar choiceworks well only for wavelengths sufficiently long that eq. 19 is obeyed, however it is a simplematter to tune each δ k individually to optimize simulation efficiency. In practice, maximaljump sizes were tuned to insure that the acceptance ratio of all trials was approximately1 / B. Fourier space Brownian dynamics (FSBD)
The FSBD method has been fully detailed in our previous work [26, 34, 35]. Here wepresent only a minimal discussion to introduce the method. In this section we consider thesimple case when p ij in Eq. (1) is independent of time and postpone the discussion of p ij ( t )for the next subsection. 10ntegrating Eqs. (11) and (7) from t to t + ∆ t for small ∆ t , we have h k ( t + ∆ t ) = h k ( t ) + Λ k [ F k ( t )∆ t + R k (∆ t )] , r i ( t + ∆ t ) = r i ( t ) + Dk B T [ F i ( t )∆ t + R i (∆ t )] . (22)where R k (∆ t ) = Z t +∆ tt dt ′ ζ k ( t ′ ) , R i (∆ t ) = Z t +∆ tt dt ′ ζ i ( t ′ ) . (23)The statistical properties of R k (∆ t ) and R i (∆ t ) follow directly those of ζ k ( t ) and ζ i ( t ), h R k (∆ t ) i = 0 , (24) h R k (∆ t ) R k ′ (∆ t ) i = 2 k B T L ∆ t Λ k δ k , − k ′ , (25) h R i (∆ t ) i = 0 , (26) h R i (∆ t ) · R j (∆ t ) i = 4( k B T ) ∆ t/Dδ i,j . (27)In the simulation, R k (∆ t ) and R i (∆ t ) are drawn from Gaussian distributions with meansand variances specified above. Since h k = h ∗− k must be satisfied to insure the height fieldof the membrane is real valued, only about half of the R k need to be generated in eachtime step; the remaining follow via complex conjugation. The precise formulation of thisstatement is somewhat complicated by the finite number of modes in our discrete Fouriertransform. Readers are referred to ref. [38] for a detailed discussion of how the full set ofrandom forces are to be generated while preserving the real valued nature of h ( r ).At each time step of the FSBD simulation, the following calculations are performed:(1) Take h ( r , t ) and r i ( t ) from the last time step.(2) Evaluate F ( r , t ) and F i ( t ) using Eqs. (4) and (8).(3) Compute F k ( t ) by Fourier transforming F ( r , t ).(4) Draw R k (∆ t ) and R i (∆ t ) from the Gaussian distributions specified above.(5) Compute h k ( t + ∆ t ) and r i ( t + ∆ t ) using Eq. (22).(6) Get h ( r , t + ∆ t ) through inverse Fourier transformation for use in the next iteration.There is one complication in step 2 of the above procedure. While we evolve each r i ( t ) asa continuous variable, the height field of the membrane is only specified at points on the realspace lattice (more precisely, it is only readily obtainable via Fast Fourier Transformationat these points). To compute the forces for use in eqs. 4 and 8, both r i and h ( r i ) are ap-proximated by assuming the position of node i directly coincides with the nearest real space11attice site. While this approximation could potentially cause problems due to discontinuousjumps in the forces, the surfaces under study are only weakly ruffled. It was verified (inthe thermal case) that FMC simulations with node positions strictly confined to lattice sitesand FSBD simulations as outlined here were in good agreement. In practice, we choose ∆ t small enough that further reduction has no consequence for the reported results, typicallyon the order of 0 . µ s depending upon choice of lattice size a . C. Kinetic Monte Carlo (KMC) for spectrin (re)association kinetics
For the case of a dynamic network, every ξ ij ( t ) jumps back and forth between 1 and 0,according to the two rates k c and k d (see Sec. II). We use the stochastic simulation algorithmof Gillespie (kinetic Monte Carlo) [39] to pick random times for breaking and reformationevents to occur. Let us consider the event of reconnecting a link. If the link is disconnectedat time t = 0, then the waiting time distribution for link reconnection is W ( t ) = k c e − k c t [40].A random number consistent with this distribution is obtained by applying the followingtransformation to a uniform random deviate x [41] t = − k c ln x. (28)A similar transformation, replacing k c with k d is used to determine the time at which anintact filament breaks apart. The set of times so obtained provides a complete trajectoryfor the behavior of each ξ ij for use in eqs. 4 and 8. The continuously distributed times arerounded off to the nearest time point in the discrete FSBD procedure. We note that ourmodel assumes each link must always connect the same two nodes (or be broken) - there isno provision for a spectrin filament to dissociate from one node and reconnect elsewhere. IV. EFFECTIVE SURFACE TENSION IN THE PRESENCE OF AN INTACT CY-TOSKELETAL MESHWORK
In this section we assume a fully connected cytoskeleton meshwork ( ξ ij = 1 for all linksat all times). In this limit, there are no non-thermal effects and either KMC or FSBDsimulations may be performed to calculate the equilibrium spectrum, h| h k | i . The qualitativefeatures of this fluctuation spectrum have been predicted by Fournier et. al. [9]. They12rgued that the behavior of the composite membrane surface should behave simply in twolimits. In the short wavelength limit, the effect of the cytoskeleton might be expectedto play a minor role; neglecting the cytoskeletal terms in eq. 1 leads to the prediction h| h k | i = k B T L x L y /κk . In the long wavelength limit, the cytoskeleton should play animportant role, but may be regarded as a continuous medium imparting an effective surfacetension to the bilayer (and possibly modifying the bare bilayer bending rigidity) as in eq. 19.“Short” and “long” wavelengths referenced above are understood to be interpreted relativeto the size of the individual spectrin links ( A , see fig 4). The original work of Fournier et.al. assumed a sharp transition between these two regimes and fit experimental data with atransition at a wavelength of approximately 2 πA . The simulations below are in qualitativeagreement with this picture, but place the transition wavelength at A and predict a finitewidth to the crossover region.Both FMC and FSBD simulations were performed with identical results (as expected).Simulations were seeded from an initially flat membrane with the cytoskeletal anchor pointsarranged in a perfect lattice (as indicated in fig. 4). The initial configuration was allowedto fully equilibrate before collecting any data for analysis. The real-space lattice constantused in the simulations was a = A/ L x = 32 a and L y = 42 a for thehexagonal network and L x = L y = L = 32 a for the square network (except where indicatedotherwise, these values of a and L x , L y and L were used in all reported simulations). Thiscorresponds to 96 independent nodes (192 triangular corrals within the periodic box) in thecase of six-fold connected anchors and 64 independent nodes (64 square corrals within theperiodic box) in the case of the four-fold connected anchors. In preliminary runs, it wasverified that neither increasing the sample size ( L x , L y ) nor decreasing a by a factors of 2significantly altered results; the values outlined above were thus chosen to insure convergedresults with minimal computational expense. For convenience, all physical parameters usedin the simulations are listed in Table I.Our simulation results are shown in Fig. 5. In the long wavelength limit, the results arewell fit by eq. 19, assuming κ eff = κ and using σ s = µ, σ h = √ µ. (29)for the effective tension in the square and hexagonal symmetry simulations, respectively.These tension values are not fit constants, but rather may be inferred from the cytoskeletal13ontribution to eq. 1. Surface tension is an energy per unit area, so we may calculate itsvalue as the ratio of entropic spring (cytoskeleton) energy per corral to the area per corral.In the square geometry, there are effectively two springs per corral (each spring is shared bytwo adjacent corrals) and using an idealized zero temperature geometry, a single corral hasarea A and total spring energy of 2 × µ A . The reported value for σ s follows immediately.A similar calculation leads to the somewhat larger value of σ h in the hexagonal geometry,reflecting the higher density of springs in this case. The numerical value of σ h so calculatedis 1 . × k B T µ m − = 5 . × − J m − , which is close to a theoretical fit [9] of theexperimental result [11] (see Fig. 5).The fluctuation spectra in fig. 5 indicate that free membrane predictions hold reasonablywell out to wavelengths of approximately A ( k − = A/ (2 π ) ∼ . µ m), with good con-vergence to long-wavelength behavior established by 5 A ( k − ∼ . µ m). The intermediateregime between wavelengths of A to 5 A encompasses the crossover between the two limitingcases. While this fact is unfortunate in light of the experimental data for the RBC [11],which displays a crossover at longer wavelengths, the simulation predictions are unambigu-ous. The experimental data remains a mystery, but we do note that it may be accountedfor by adoption of an ad hoc harmonic confining potential [7].One of the motivations for this work was to test the performance of the analytical theorydeveloped in DF. In particular, we anticipated that the mobility of cytoskeletal anchorswould lead to some quantitative deviations from the DF theory at short wavelengths. Infact, the mobility of anchor points leads to insignificant changes in the fluctuation spectrumfor physical parameters relevant to the RBC (see fig. 6). We also note that the detailedanalysis of DF does slightly better in reproducing the simulated spectrum than does theadoption of surface tensions implied by eq. 29 (see fig. 7). At intermediate wavelengths,the small deviations of κ eff away from κ predicted by DF do improve the fit as comparedto our naive arguments, however the effect is very slight in comparison to the leading ordereffect of introducing a finite surface tension. As a final point, we note that the anisotropicnature of the square network over all wavelengths leads to some variance in fluctuationsfor a given magnitude of k depending on the direction taken. In principle, the hexagonalnetwork should not suffer from this effect [31], however the underlying square lattice takenfor our simulations does introduce some anisotropy to the spectra; these effects are mostsevere at short wavelengths. The spread of data points plotted in figs. 5 - 7 reflects this14irectional dependence in the spectra and should not be taken as evidence of statisticalnoise. As previously mentioned, statistical errors are less than the size of the symbols usedin plotting. V. EFFECTIVE SURFACE TENSION IN THE PRESENCE OF A CYTOSKELE-TAL MESHWORK WITH DYNAMICALLY EVOLVING CONNECTIVITY
We now generalize the results of the previous section to include RBC membranes withrandomly (and dynamically changing) broken cytoskeletal meshworks. As outlined in sectionII, we assume the dynamics of spectrin association and disassociation are governed by simplerate processes. The rate constants for connecting a broken cytoskeletal filament, k c , andbreaking an intact filament, k d , are assumed independent of the distance between filamentend points on the bilayer surface and independent of all other connections within the network.This immediately leads to the conclusion that the probability for a filament to be intact ata given time is p = k c / ( k c + k d ). Equivalently, p is the average percentage of intact filamentsin the cytoskeletal network. For the moment, we simply take this picture as a hypotheticalmodel for non-equilibrium dynamics in the RBC membrane. A discussion regarding theconnection between this model and experiment will be provided in sec. VI.Our analysis in this section centers around the calculation of σ ( p ), the effective tensionof the composite membrane as inferred from long wavelength fluctuations. As indicated bythe notation, this tension depends on the degree of connectivity within the network. Notapparent in the notation is the fact that this tension also depends upon the magnitude ofthe rates k d and k c (and not just the ratio of them in p ) due to the non-equilibrium natureof the dynamics. Extraction of σ ( p ) follows the same general prescription as in the previoussection; the fluctuation spectrum is collected and compared to the empirical result of eq.19. At the longest wavelengths modeled in the simulations, it is found that eq. 19 does agood job of reproducing the simulation data. Obvious theoretical estimates for σ eff = σ ( p )are only available in the limiting cases of very fast spectrin (re)association kinetics and veryslow kinetics. In general, the effective tension as a function of p (and kinetic rates) must beextracted from simulation.In the limit that spectrin (re)association rates are much faster than all other time scales inthe problem, each cyotoskletal link in the network behaves as an intact link with a diminished15pring constant. The numerical value for this weakened spring constant is simply the timeaverage of the spring as it flips between the two possible values of µ and zero. This pictureimmediately leads to σ s ( p ) = pµ, σ h ( p ) = √ pµ. (fast spectrin kinetics) (30)Similar arguments have been invoked previously to suggest a possible ATP dependence inthe shear modulus of the RBC membrane [19]. Such an ATP concentration dependence mayprovide an explanation for RBC shape changes as a function of ATP concentration [19].In the opposite regime, when network connectivity kinetics are far slower than any othertimescale in the system, the membrane may be regarded as evolving thermally under theinfluence of quenched network disorder. The behavior of the system in the quenched disorderlimit is analogous to the 2D percolation problem [42]. As such, it is expected that theeffective tension within the system must vanish at some finite critical p = p c . For values of p less than p c , no global connectivity within the network exists and, consequently, there isno restoring force possible in response to area dilations of the membrane surface. The 2Dconnectivity percolation limit is known to occur at p c = 2 /z , where z is the bond valencefor each node (i.e. z = 4 with p c = 1 / z = 6 with p c = 1 / σ as p drops from one to p c is linear. That is σ s ( p ) = 2 µ ( p − / , σ h ( p ) = 3 √ µ ( p − / p > p c and σ = 0 for p < p c . A brief justification for the percolation theory resultsquoted above may be found in the Appendix.The discussion of “fast” and “slow” kinetics in the previous paragraph was intentionallyleft ambiguous, without specification of what time scales these quantities were to be com-pared with. It seems prudent to avoid a detailed discussion of this issue, due to the manydifferent dynamic scales in this particular problem, however a crude discussion is appropri-ate. Given our focus on long wavelength elastic properties, the most obvious timescale forcomparison is the membrane relaxation time for the longest wavelength under observation.For equilibrium membranes at tension σ and bending rigidity κ , it is well known [5, 25] that16elaxation of h k is exponential with a characteristic time τ m ( k ) = 4 ηκk + σk . (32)This result follows immediately from eq. 11. Considering only the longest wavelengthmodes in our simulations ( k = 2 π/L max with L max = L for the square and L max = L y forthe hexagonal simulations) and allowing for tensions between zero and the fully connectednetwork values, we find that τ m falls in the range of 9 × − s - 8 × − s. While these valuesare only required to hold for homogeneous equilibrium membranes, they were verified to holdreasonably well for the non-equilibrium membranes studied here when σ ( p ) (determined fromthe fluctuation spectrum) is naively used in eq. 32. The relaxation rate associated with thetwo-state spectrin dynamics is ( k c + k d ), which provides a time scale τ cd = 1 / ( k c + k d ). Thefast kinetics limit discussed above would be expected to apply for τ m ≫ τ cd and the slowkinetics limit would apply for τ m ≪ τ cd .FSBD simulations were run including KMC breaking and reforming events in the spec-trin meshwork (as detailed in sec. III). Although no equilibrium can be reached in thesesimulations due to the intrinsically non-equilibrium nature of the simulation, the systems doconverge to a steady state regime. h| h k | i as a function of k were collected, at steady-state,for both the square and hexagonal systems described in the previous section. Three differ-ent sets of simulations were run in each geometry corresponding to different spectrin kineticregimes. Specifically, k d was chosen to assume the values 10 s − , 200 s − , and 10 s − . Avariety of different connectivity percentages were simulated, spanning 0 < p <
1. Givenvalues for k d and p , k c = k d p/ (1 − p ) follows immediately allowing for complete specificationof the model. The three choices for k d roughly correspond to the regimes τ cd ≪ τ m , τ cd ∼ τ m ,and τ cd ≫ τ m respectively.Figures 8 and 9 display our simulation results, plotted in the form of effective surfacetensions as a function of network connectivity. The tensions are calculated as in the lowerpanel of Fig. 5 using the largest wavelengths available in the simulations. In addition tothe FSBD/KMC simulations, FMC results are plotted for the quenched disorder case corre-sponding to ( k c + k d ) = 0. In these simulations, an ensemble of random network connectivitiesconsistent with the prescribed p values were run without allowing the network connectivityto evolve. The results for h| h k | i were averaged over the ensemble to obtain σ ( p ). FMC wasused in these cases for numerical efficiency and is valid since the evolution under conditions17f static network connectivity is purely thermal. The speed of the FMC algorithm alloweda set of simulations to be run for larger membrane sizes in order to clarify the role of finitesize effects.The figures clearly display both expected regimes. Fast network reorganization yieldsresults in good agreement with eq. 30, whereas slow reorganization yields results consistentwith eq. 31. Intermediate kinetic regimes interpolate between these two limiting cases.Finite size effects appear not to play any role, except perhaps for the case of p values veryclose to p c in the quenched disorder case, which is to be expected due to the percolation phasetransition at this point. At values of intermediate membrane connectivity, it is possible forthe systems to display a range of effective surface tension values, depending upon the relativetimescales for spectrin kinetics as compared to membrane relaxation. The implications ofthis fact will be discussed further in the next section. VI. CONNECTION TO EXPERIMENTS
Evidence that RBC membrane shapes and shape fluctuations depend upon non-thermalenergy expenditure comes from two different experimental sources. Measurement of cellshape [43] and shape fluctuations [17, 18] under a variety of MgATP concentrations suggestthis fact directly. Increases in MgATP concentration lead to enhanced membrane surfacefluctuations. Less directly, it has been shown that RBC membrane fluctuation amplitudes inthe presence of MgATP depend upon the viscosity of the surrounding solvent [20]. Thermalproperties of a physical system should not be influenced by transport coefficients, such asthe viscosity, but should depend only upon system energetics via the partition function.Biochemically, it is clear that the presence of MgATP leads to the phosphorylation ofspectrin [15], which has been implicated in the observed shape changes and fluctuationamplitudes in RBC membranes under varying MgATP concentrations [43]. Whether or notthis phosphorylation event (and/or similar events in other molecular components of thenetwork) coincides with breaking of spectrin filaments is less clear, however the hypothesisthat this is in fact the case [19] is compelling and served as one of the major motivations forthis study. In what follows, we assume that phosphorylation of spectrin induces individualfilaments to dissociate from the network. We stress that this model is hypothetical and,indeed, that some studies refute the correlation between spectrin phosphorylation and RBC18hape/elastic properties [44, 45].In a naive model for spectrin phosphorylation, we assume the following kinetic equationsfor spectrin dissociation and association with the network as a whole.S a + ATP k d → S d + ADP (33) S d + H O k c → S a + HPO − . In the above, S a and S d stand for the associated and dissociated forms of spectrin, respec-tively, and it is assumed that the reactions proceed via catalysis by kinases and phosphatases.We assume ATP concentration is held fixed, either by endogenous biochemistry within theRBC in vivo or by experimental control in vitro. Using our previously introduced notation,we write [S a ] = Cp and [S d ] = C (1 − p ), where C is the total concentration of spectrin on thecell surface. The steady state condition is formulated as dp/dt = 0, which, in conjunctionwith kinetic eqs. 33, leads to p = n c n c + 2 n . (34)Here, both kinetic constants have been wrapped up into the single constant n c . n is theconcentration of ATP in solution. n c has been defined such that the condition n = n c implies p = p c = 1 / n = n c specifies the concentration of ATP necessaryto reach the percolation phase transition in systems with quenched disorder. While thismodel neglects the observed dependence of spectrin dephosphorylation on n [15], it hasthe advantage of simplicity in the absence of quantitative kinetic data and follows the lineof reasoning previously introduced [19], allowing for comparison to earlier work. Eq. 34predicts the connectivity of the network as a function of ATP concentration and allows forcomparison between the modeling of the previous section and experimental data.Experimental data for RBC fluctuations as a function of ATP concentration is availableonly in terms of real space height amplitudes, h ( r ) [18] (see fig. 10). Following Eqs. (12)and (19) (and for simplicity assuming L x = L y = L ), h h ( r ) i = X k k B TL ( κk + σk ) , (35)where the interval between adjacent k x or k y is 2 π/L . In the limit of large L , the sum P k can be approximated by the integral L / (2 π ) R d k , and we have h h ( r ) i = k B T πσ ln(1 + σL π κ ) . (36)19he limiting expressions and simulations of the preceding section provide a route towardestimation of σ ( p ) under various rates of spectrin association/dissociation. Eq. 34 enablesus to translate these results into tensions as a function of ATP concentration. The resultingtensions may be used in eq. 36 to calculate the real space membrane fluctuation amplitudesas a function of ATP concentration. RBC’s have a diameter of 7 microns and we use L = 7 µ m [46] in our comparison to experiment. It is also important to note that experimentalobservations of RBC fluctuations were carried out on cells that were allowed to “firmly adhereto [a] glass substratum” [18] under the influence of their natural (presumably of electrostaticorigin) affinity for such surfaces. The adhesion of the bilayer to an underlying supportingmatrix will be assumed to contribute a bare surface tension, σ bare to the membrane of theRBC [47, 48]. That is, we assume σ = σ bare + σ ( p ).We take n c and σ bare as two fitting parameters to reproduce the experimental data. Theresulting fits, assuming the two extreme cases of fast spectrin kinetics (eq. 30 assumed)and slow spectrin kinetics (eq. 31 assumed) are displayed in fig. 10. The values adoptedby the fitting parameters in these two cases are σ bare = 0 . µ = 150k B T /µ m (both cases),and n c = 0 .
22 mM for the fast kinetics case and n c = 0 . n c naturally leadsto the observed plateau in the experimental data.The slowest possible time-scale associated with RBC membrane fluctuations in our model( τ m ) is on the order of a couple of seconds. This assumes a tension-free membrane and L = 7 µ m. Our data analysis suggests that spectrin breaking/reformation dynamics aresubstantially slower than this, since the quenched-disorder (percolation theory) results pro-vide the best fit to the experimental data. This finding would appear to be consistent withthe experimental observation [15] that establishment of steady-state levels of spectrin phos-phorylation occurs on the time scale of minutes to tens of minutes with changes in MgATPconcentration. Previous theoretical work exploring the consequences of spectrin dissociationon membrane fluctuations [19] has assumed the time-scale for spectrin kinetics to be on theorder of milliseconds. Within the context of the present model, rates this fast seem unlikely.Experimentally, it is observed that RBC fluctuation amplitudes depend not only uponMgATP concentration, but also upon the viscosity of the solvent surrounding the mem-20rane [20]; increasing viscosity above that of aqueous buffer solution leads to decreases inobserved fluctuations. Membrane fluctuations at thermal equilibrium should yield identicalstatistics regardless of solvent viscosity (although the dynamics will certainly differ), so theseexperiments provide additional proof of the non-thermal character of RBC fluctuations.Within the present model, Eq. (32) clearly displays the linear relationship between inverseviscosity and the relaxation timescales for individual membrane modes. In the precedingsection, it was argued that the ratio between membrane relaxation rates and spectrin kineticrates determines the magnitude of effective tension that must be used if one desires tofit membrane fluctuations to the functional form of eq. 19. At a given connectivity, p ,fluctuation amplitudes decrease ( σ eff increases) as the spectrin kinetic cycle increases inspeed. Alternatively, since it is the ratio of kinetics to membrane relaxation that is relevant,we expect fluctuation amplitudes to decrease as viscosity is raised assuming all other systemparameters are held fixed. Qualitatively, this trend is in agreement with the experimentalresults.To obtain quantitative comparison with the variable viscosity data it is necessary toevaluate h h ( r ) i explicitly from direct simulation, without assuming a single effective tensionas in eq. 35. Within our model, fluctuation amplitudes are affected by viscosity changes onlyby shifting the timescale for membrane dynamics relative to spectrin kinetics. A relativeshift in timescale (at constant p ) corresponds to a vertical motion between the limitingregimes plotted in fig. 9. In order for viscosity changes to yield any effect, the timescales forspectrin kinetics and membrane relaxation must be comparable. (e.g. if spectrin dissociationtakes an hour, increasing viscosity by a factor of 10 will not remove you from the limiting“slow kinetics regime”. Similarly, if dissociation took a nanosecond, no reasonable viscositychange could promote the system out of the “fast kinetics regime”.) And, if these timescalesare comparable, differences in the individual relaxation rates at each wavelength will lead todifferent effective tensions as a function of wavelength. In practice, it is not feasible to runsimulations for a full 7 µm × µm patch with sufficient sampling to obtain reliable results.It is possible to run a 1 µm × µm patch. The effective tensions as a function of wavevectorwere extracted from this rectangular geometry and were used in a generalized version of eq.35 with σ → σ ( k ) to obtain h h ( r ) i for the full 7 µm × µm system.To simultaneously find agreement with the ATP concentration data (fig. 10) and thevariable viscosity data it is necessary to assume spectrin rates that are slow relative to21embrane relaxation when the solvent is pure buffer solution (as in the conditions of ref.[18] and fig. 10), but are poised to move out of this limiting regime with small viscosityincreases. By trial and error, it was found that τ cd = 5 × µ s for n = 1 . L = 7 µ m, n c = 0 . σ bare = 150k B T /µ m .In fig. 11 we plot real space membrane fluctuations as a function of η − for viscosity valuesspanning an order of magnitude. We observe a reduction in fluctuation amplitudes of about20% over the entire range of viscosities studied, which falls just outside the errors of theexperimental data [20]. It is important to stress that figs. 10 and 11 simultaneously fit twovery different experiments to a single set of physical parameters. Our results agree withboth experiments quite well. VII. DISCUSSION
From a mathematical standpoint, the results of the preceding section could be viewedas a success; we fit the available experimental data with our theoretical model. However,the physical implications of the derived fit parameters are troubling. Most worrisome is thecritical concentration of ATP for onset of the percolation limit, n c = 0 . n ∼ n = 1 . n ≥ n c , since we donot allow the connectivity of the network to evolve away from its starting point (i.e. a givenspectrin tetramer is assumed to always link the same two vertices or be disassociated). Ina severely compromised network, how could a given filament ever be expected to find itsassigned association point following a dissociation event?As mentioned previously, there does not appear to be definitive experimental data relat-ing ATP consumption to spectrin network dissociation. Rather, this is a plausible hypothesisadvocated in reference [19], based upon the limited biochemical data that is available forthe various constituents of the RBC cytoskeleton. We can not reconcile this picture withthe simulations carried out in this work, primarily for the reason outlined in the previous22aragraph. It is possible that ATP consumption could lead to a significant change in theelastic properties of individual spectrin filaments without completely dissociating the fil-ament from the network. In such a picture, the two forms of spectrin introduced above, S d and S a , would correspond to filaments with weak and strong effective spring constants,respectively (or different natural lengths, etc.). In this context, the analysis we present issensible since network connectivity is always maintained, but the elastic properties of indi-vidual links are changing. The percolation limit in this picture corresponds to the point atwhich it is impossible to follow a path across the spectrin network without stepping on atleast one weak S d link. Although such a picture is completely hypothetical, we note thatit is theoretically possible to alter the effective elastic properties of spectrin filaments whilemaintaining lateral integrity of the network [49, 50].The underlying physical picture pursued in this work was motivated by and is similar inspirit to the work of Safran and Gov [19] (i.e. a dynamically breaking and reforming spectrinnetwork). However, the simulation model implemented to study this system is more similarto the work of Dubus and Fournier [21] (with the additional possibility of dynamicallybreaking bonds). We have interpreted our simulation results within the context of a simplepercolation theory. While this picture shares some similarity with the analysis of ref. [19] inthe context of global effects associated with variable ATP concentration, we find differencesbetween the present model and ref. [19] in the description of local fluctuation dynamics.For example, we find no evidence in our simulations for localized forces of the sort discussedin ref. [19]. Instead, the increased amplitudes of fluctuation at high ATP concentration areattributable to global properties of the spectrin network within our model. We also findthat the rates associated with spectrin kinetics must be orders of magnitude slower thanassumed in ref. [19] in order to reconcile our simulation model with the experiments of refs.[18, 20]. However, it should be pointed out that the particular rates assumed in ref. [19]were only rough estimates and are not essential to the qualitative predictions of that work[51].The starting point for all of our simulations, eq. 1, assumes a central-force network for thecytoskeletal matrix and a phantom entropic spring description for the spectrin filaments. Inaddition, we do not allow for the possible formation of 5-fold or 7-fold defects [19, 52] in ournetwork. Inclusion of direct interactions between spectrin and the bilayer surface, generaliz-ing the description of the network to deviate from the central-force picture, and/or allowing23or defects in network connectivity could all potentially alter the quantitative findings ofthis work. The present simulation model was adopted for ease of numerical implementationand to facilitate comparison to the existing work of Dubus and Fournier [21]. This studyprovides a detailed analysis of the behavior of a specific model for a composite membranesystem (lipid bilayer plus actively labile cytoskeletal network). This model is motivated byour (incomplete) picture for the structure and kinetics of the RBC membrane surface andmay prove useful in the development of more refined models in the future. In particular, itshould be noted that more refined descriptions of the spectrin network have been developedthat allow for interactions between the bilayer surface and the filaments and a more completedescription of polymer elasticity beyond the simple Gaussian-chain model [53, 54, 55]. Fu-ture modeling, incorporating some of these approaches for the behavior of the cytoskeleton,would be highly desirable and is currently under investigation.We note in closing that a very recent study by Evans et. al. [56] directly contradicts theexperimental results in refs. [18, 20]. This recent series of experiments finds no correlationbetween ATP concentration and RBC membrane fluctuations. Given the severe disagree-ment between experimental results on the RBC system, one has almost complete freedomin interpretation of our simulation results. The picture outlined above suggests that we canobtain quite good agreement between simulation and experiments that do measure ATPdependence of the RBC fluctuations. Within this picture, it seems necessary to assume thatour “broken” spectrin links actually correspond to intact links with diminished spring con-stant and/or a non-zero natural length. Another, equally valid, interpretation of our resultsis that it is impossible to reconcile an actively breaking cytoskeletal meshwork picture withavailable experimental data due to the high degree of network dissociation predicted underphysiological conditions within such a picture. Given the uncertain experimental landscape,it seems impossible to make a definitive statement in favor of either viewpoint. Acknowledgments
We thank Golan Bel, Nir Gov and Kim Parker for useful discussions. This work wassupported by the National Science Foundation (grant No. CHE-0321368 and grant No.CHE-0349196). F. B. is an Alfred P. Sloan research fellow and a Camille Dreyfus Teacher-Scholar. 24 ppendix
The calculation of an effective tension in the limit of quenched disorder within the spectrinnetwork is closely related to the 2D bond percolation problem [57]. A mean field treatmentfor the percolation problem was first developed in the context of electric conductivity [42]and was later extended to calculate the elastic moduli of networks of Hookean springs withfinite [58, 59] and zero natural lengths [60, 61]. The below discussion summarizes thearguments of Ref. [58] for the readers’ convenience.Let us consider a randomly connected meshwork with only a fraction ( p ) of the linksintact. Each intact link is a spring with spring constant µ and a natural extension ofzero. In the spirit of mean field theory, we assume that the global elastic properties of thisimperfect meshwork are well approximated by a complete meshwork comprised of springswith spring constant µ m at every link (see Fig. 12) even though individual links will beeither completely intact or fully broken. Our problem is to determine an expression for µ m .For this purpose, let us consider an arbitrary link AB in the meshwork. The real springconstant, µ AB , associated with this link is either µ if the link is connected, or 0 if it is not.If we apply a force f m across this link (see Fig. 12), the extension x AB is given by x AB = f m µ AB + µ eff , (37)where µ eff is an effective spring constant reflecting the contribution of all network compo-nents except the link AB . It can be shown that [58] µ eff = ( z/ − µ m , (38)where z is the connectivity of the meshwork, for example, z = 4 for the square meshworkand z = 6 for the hexagonal meshwork considered in this work.On the other hand, if we assume the mean spring constant µ m for link AB , we predict amean extension, x m = f m µ m + µ eff . (39)To achieve consistency within the mean field approximation, we must have h x AB i = x m (40) pf m µ + µ eff + (1 − p ) f m µ eff = f m µ m + µ eff µ m = p − p c − p c µ, (41)where p c = 2 /z. (42)To calculate σ ( p ), we simply replace µ in eq. 29 by µ m from eq. 41, since the latter is nowthe average spring constant of the equivalent complete meshwork. Taking z = 4 and z = 6for the square and hexagonal meshwork respectively, we arrive at eq. 31.One should note that eqs. 31 and 42 are only results of a mean field approximation. Inthe vicinity of p c , it is well known that the correlation effects are strong and that mean fieldtheory breaks down. However, the mean field results are adequate to describe the behaviorof σ outside the immediate vicinity of p c as evidenced by figs. 8 and 9, which is sufficientfor the purposes of this work. [1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology ofthe Cell (Galland, New York, 2002).[2] S. A. Safran,
Statistical Thermodynamics of Surfaces, Interfaces, and Memebranes (WestviewPress, Boulder, CO, 2003).[3]
Structure and Dynamics of Membranes , edited by R. Lipowsky and E. Sackmann (ElsevierScience, Amsterdam, 1995).[4] S. Liu, L. Derick, and J. Palek, J. Cell. Biol. , 527 (1987).[5] F. Brochard and J. F. Lennon, J. Phys. (Paris) , 1035 (1975).[6] M. A. Peterson, Phys. Rev. A , 4116 (1992).[7] N. Gov, A. G. Zilman, and S. Safran, Phys. Rev. Lett. , 228101 (2003).[8] G. H. W. Lim, M. Wortis, and R. Mukhopadhyay, Proc. Nat. Acad. Sci. USA , 16766(2002).[9] J.-B. Fournier, D. Lacoste, and E. Raphael, Phys. Rev. Lett. , 018102 (2004).[10] S. B. Rochal and V. L. Lorman, Phys. Rev. Lett. , 248102 (2006).[11] A. Zilker, H. Engelhardt, and E. Sackmann, J. Phys. (Fr.) , 2139 (1987).[12] D. Discher, N. Mohandas, and E. A. Evans, Science , 1032 (1994).
13] V. Heinrich, K. Ritchie, N. Mohandas, and E. Evans, Biophys. J. , 1452 (2001).[14] H. Engelhardt, H. Gaub, and E. Sackmann, Nature (London) , 378 (1984).[15] W. Birchmeier and S. J. Singer, J. Cell Biol. , 647 (1977).[16] V. Bennett, Biochim. Biophys. Acta , 107 (1989).[17] S. Levin and R. Korenstein, Biophys. J. , 733 (1991).[18] S. Tuvia, S. Levin, A. Bitler, and R. Korenstein, J. Cell Biol. , 1551 (1998).[19] N. Gov and S. A. Safran, Biophys. J. , 1859 (2005).[20] S. Tuvia, A. Almagor, A. Bitler, S. Levin, R. Korenstein, and S. Yedgar, Proc. Natl. Acad.Sci. USA , 5045 (1997).[21] C. Dubus and J.-B. Fournier, Europhys. Lett. , 181 (2006).[22] W. Helfrich, Z. Naturforsch. pp. 693–703 (1973).[23] D. R. Nelson and L. Peliti, J. Phys. (Paris) , 1085 (1987).[24] E. M. Purcell, American Journal of Physics , 3 (1977).[25] R. Granek, J. Phys. II (Paris) , 1761 (1997).[26] L. C.-L. Lin and F. L. H. Brown, Phys. Rev. E , 011910 (2005).[27] F. L. H. Brown, Annual Review of Physical Chemistry , 685 (2008).[28] M. Doi and S. Edwards, The Theory of Polymer Dynamics (Clarendon, Oxford, 1986).[29] E. Reister and U. Seifert, Europhys. Lett. , 859 (2005).[30] A. Naji and F. L. H. Brown, J. Chem. Phys. , 235103 (2007).[31] D. H. Boal, Mechanics of the cell (Cambridge University Press, Cambridge, UK, 2002).[32] N. Gouliaev and J. F. Nagle, Phys. Rev. E. , 881 (1998).[33] N. Gouliaev and J. F. Nagle, Phys. Rev. Lett , 2610 (1998).[34] L. C.-L. Lin and F. L. H. Brown, Biophys. J. , 764 (2004).[35] L. C.-L. Lin and F. L. H. Brown, Phys. Rev. Lett. , 256001 (2004).[36] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, San Diego,2002), 2nd ed., pp. 525-532.[37] D. Bouzida, S. Kumar, and R. H. Swendsen, Phys. Rev. A , 8894 (1992).[38] L. C.-L. Lin and F. L. H. Brown, J. Chem. Theory Comp. , 472 (2006).[39] D. T. Gillespie, J. Comput. Phys. , 403 (1976).[40] N. G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amster-dam, 1992), pp. 63,83,220-221.
41] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
Numerical Recipes in C (Cambridge University Press, Cambridge, 1994).[42] S. Kirkpatrick, Rev. of Mod. Phys. , 574 (1973), and references therein.[43] M. P. Sheetz and S. J. Singer, J. Cell Biol. , 638 (1977).[44] V. P. Patel and G. Fairbanks, J. Cell. Biol. , 430 (1981).[45] L. Bordin, G. Clari, I. Moro, F. D. Vecchia, and V. Moret, Biochem. Biophys. Res. Comm. , 249 (1995).[46] H. Lodish, D. Baltimore, A. Berk, S. L. Zipursky, P. Matsudaira, and J. Darnell, MolecularCell Biology (Scientific American Books, New York, 1995), 3rd ed.[47] F. Brochard, P. G. DeGennes, and P. Pfeuty, J. Phys. (Paris) , 1099 (1976).[48] L. C.-L. Lin, J. T. Groves, and F. L. H. Brown, Biophys. J. , 3600 (2006).[49] N. S. Gov, New Journal of Physics , 429 (2007).[50] Q. Zhu and R. J. Asaro, Biophys. J. , 2529 (2008).[51] N. S. Gov, personal communication.[52] H. S. Seung and D. R. Nelson, Phys. Rev. A , 1005 (1988).[53] D. H. Boal, Biophys. J. , 521 (1994).[54] J. Li, G. Lykotrafitis, M. Dao, and S. Suresh, Proc. Nat. Acad. Sci. USA , 4937 (2007).[55] Q. Zhu, C. Vera, R. J. Asaro, P. Sche, and L. A. Sung, Biophys. J. , 386 (2007).[56] J. Evans, W. Gratzer, N. Mohandas, K. Parker, and J. Sleep, Biophys. J (2008), in press.[57] D. Stauffer and A. Aharony (1994).[58] S. Feng, M. F. Thorpe, and E. Garboczi, Phys. Rev. B , 276 (1985).[59] S. Feng and P. N. Sen, Phys. Rev. Lett. , 216 (1984).[60] W. Tang and M. F. Thorpe, Phys. Rev. B , 3798 (1987).[61] W. Tang and M. F. Thorpe, Phys. Rev. B , 5539 (1988).[62] M. Tomishige, Y. Sako, and A. Kusumi, J. Cell Biol. , 989 (1998). ABLE I: Physical parameters in simulationsParameter Description Value Reference κ bending modulus 4.8 k B T a [5] η cytoplasm viscosity 1.5 k B T s µ m − [5] D diffusion constant of nodes 0.53 µ m s − [62] A the average distance between two nearest nodes (see Fig. 4) 0.112 µ m [31] ℓ c contour length of the spectrin tetramer 0.2 µ m [31] ℓ k Kuhn length of the spectrin tetramer 0.02 µ m [31] µ spring constant of the spectrin tetramer 750 k B T µ m − b a T ≡ b Calculated using the Gaussian chain model through two parameters above. b)(c)(a) FIG. 1: A schematic rendering of the RBC membrane. (a): Top View: Looking down on thecell with the lipid bilayer invisible for clarity. The nodes of the cytoskeletal network are proteincomplexes embedded in the lipid bilayer. The links are spectrin tetramers, with approximatelateral end-to-end distance of 100 nm between the nodes. (b): Side View: The protein anchors areconfined to the lipid bilayer, but may diffuse laterally. The coils underneath the bilayer are spectrintetramers. The contour length of the tetramers is approximately 200 nm , so there is considerableextra length coiled up below the membrane surface. (c): One end of a link is disconnected fromthe node, possibly with the help of ATP. IG. 2: (Color online). A snapshot the RBC membrane as treated in our simulations. The lipidbilayer is connected to spectrin via a series of laterally mobile anchoring proteins (white dots). Theanchors are interconnected via harmonic potentials (not shown for clarity), modeling the entirecytoskeletal network as a series of interconnected entropic springs. Our model does not account forany steric repulsion between cytoskeleton and bilayer surface, nor does it account for inextensibilityof spectrin beyond its natural contour length; the individual spectrin links are treated as simpleGaussian chains. Individual links in the network are allowed to break and reform following kineticequations that do not obey detailed balance; this provides a model for energy (ATP) consumptionat the membrane surface. A is the average mesh size ( ∼ nm ). a) (b) FIG. 3: A schematic illustration of a randomly broken spectrin meshwork. The solid lines areconnected links while the dotted lines represent disconnected ones. The percentage of intact links(equivalently, the probability for any single link to be intact) will be denoted “ p ”. All lines together(solid and dotted) correspond to a complete hexagonal meshwork. (a): Above the percolationthreshold, p > p c , the connected links have global connectivity, i.e., they form a cluster extendingacross the entire cell. The corresponding effective surface tension will always be finite. (b): Belowthe percolation threshold, p < p c , there is no global connectivity in the meshwork, i.e. connectedlinks form finite islands of connectivity that do not span the cell. The corresponding effective surfacetension becomes zero if the links are breaking and reforming slowly. (c): Even at instantaneousconnectivities below the percolation threshold, there can be a finite surface tension if spectrinbreaking/formation kinetics are sufficiently fast. In this limit, each link can be thought of as aspring with a reduced restoring force since out of plane membrane undulations evolve more slowlythan the lifetime of the link. A LLL L xy A (a) (b) FIG. 4: A schematic illustration of the connectivity of the nodes of the meshwork. (a): thehexagonal meshwork. The dotted lines are a guide of eyes for the rectangular sample used. (b):the square meshwork. IG. 5: UPPER: Spectrum of fluctuations, h| h k | i , (normalized by h| h k | i f ) versus inversewavenumber for the our RBC membrane model assuming a perfectly intact cytoskeletal network.The squares and triangles indicate simulation results for square and hexagonal meshworks respec-tively. The error bars are approximately the same size as the symbols. The solid lines are theoreticalresults using Eq. (29) for σ eff and κ eff = κ in Eq. 19. The wavevector k m = 2 π/A is indicatedon the horizontal axis to show where crossover from free membrane to long-wavelength behaviormight be expected to occur. The circles are experimental results taken from Ref. [11]. LOWER:The approach of σ to its limiting long-wavelength value is displayed. In this case, each point inthe spectrum on the left was individually fit to eq. 19 by adjusting σ eff (assuming κ eff = κ ). Theresulting σ eff as a function of k − are shown. At the longest wavelengths simulated, the inferredtensions converge to the predicted values. IG. 6: Simulation results of h| h k | i for a hexagonal lattice. The circles are the data for the caseof immobile nodes placed on a perfect lattice. While the triangles are the data for the case ofmobile nodes (identical data to fig. 5). Mobility of the cytoskeleton attachment points does notsignificantly alter the simulation results. IG. 7: Comparison of different theories for a hexagonal lattice. The circles are the simulationresults. The dashed line is the prediction introduced in fig. 5. The solid line is the DF modelprediction. At intermediate wavelengths DF does a slightly superior job in fitting the simula-tions. Both approximations fail at small wavelengths and converge to the same behavior at longwavelengths. p σ ( p ) µ FIG. 8: σ ( p ) /σ ( p = 1) versus p is plotted for a square meshwork. The circles, triangles, anddiamonds are FSBD/KMC results for k d = 10 s − ( τ cd ≪ τ m ), k d = 200 s − ( τ cd ∼ τ m ) and k d = 10 s − ( τ cd ≫ τ m ) respectively. The squares and the stars are FMC results ( k d = k c = 0). Allsimulations use box sizes identical to those introduced in the previous section, except for the stars.Stars represent FMC simulations on systems that are twice as large in both linear dimensions. Forclarity, only error bars for the circles and the stars are shown, but all simulations have similarerrors. The dashed line is the prediction of Eq. (30). The solid line is the prediction of Eq. (31). p σ ( p ) µ FIG. 9: σ ( p ) /σ ( p = 1) versus p is plotted for a hexagonal meshwork. All symbols and lines havethe same meaning as in Fig. 8, only the meshwork connectivity has been changed. (The obviousdiscrepancy between simulation and theory at p = 1 is due to the rectangular geometry of oursimulation box. It is impossible to perfectly fit a portion of a hexagonal network within a rectangle,so the theoretical results that assume perfect hexagonal symmetry are only approximately valid inthis case.) [MgATP] (mM) q h h ( r ) ih h ( r ) i FIG. 10: Real space RBC height fluctuations, p h h ( r ) i , as a function of ATP concentration(normalized by the zero ATP values, p h h ( r ) i ). The circles with error bars are the experimentaldata [18]. The solid line is a theoretical fit assuming σ h in Eq. (31), corresponding to the case of τ m ≪ τ cd (slow spectrin kinetics). The dashed line is a theoretical fit assuming σ h in Eq. (30),corresponding to the case of τ cd ≪ τ m (fast spectrin kinetics). /η (cp − ) q h h ( r ) ih h ( r ) i FIG. 11: Normalized height fluctuations, p h h ( r ) i , as a function of inverse viscosity of the solvent.The circles with error bars are the simulation results assuming L = 7 µ m, σ bare = 150k B T /µ m , n c = 0 . n = 1 . τ cd = 5 × µ s. The solid line is a fit to our simulation data.The dashed line and the two dashdotted lines are the experimental fitting line and the range ofexperimental errors taken from Fig. 2 of Ref. [20]. µ m µ m µ m µ m µ m µ m (a) (b) AB AB ff f eff mm m m f AB AB µ µ µ
FIG. 12: (a) A randomly broken network of springs is approximated by a complete meshwork ofsprings with diminished spring constant µ m . Here we imagine an attempt to extend the link AB using a force of f m . (b) The force f m is resisted both directly by link AB and the rest of themeshwork. The latter can be considered as a single spring with effective spring constant µ eff ,which acts in parallel with AB . (Adapted from ref. [58]). (Adapted from ref. [58])