Monte Carlo study of the frame, fluctuation and internal tensions of fluctuating membranes with fixed area

Three types of surface tensions can be defined for lipid membranes: the internal tension, $\sigma$, conjugated to the real membrane area in the Hamiltonian, the mechanical frame tension, $\tau$, conjugated to the projected area, and the"fluctuation tension", $r$, obtained from the fluctuation spectrum of the membrane height. We investigate these surface tensions by means of a Monge gauge lattice Monte Carlo simulation involving the exact, nonlinear, Helfrich Hamiltonian and a measure correction for excess entropy of the Monge gauge. Our results for the relation between $\sigma$ and $\tau$ agrees well with the theoretical prediction of [J.-B. Fournier and C. Barbetta, Phys. Rev. Lett., 2008, 100, 078103] based on a Gaussian approximation. This provides a valuable knowledge of~$\tau$ in the standard Gaussian models where the tension is controlled by $\sigma$. However, contrary to the conjecture in the above paper, we find that $r$ exhibits no significant difference from $\tau$ over more than five decades of tension. Our results appear to be valid in the thermodynamic limit and are robust to changing the ensemble in which the membrane area is controlled.


Introduction
A few decades of research has clarified that lipid bilayers can be treated as sheets obeying the continuum curvature/tension elasticity formulated by Helfrich. 1 Characterising membrane elastic constants in numerical simulations has become an important issue due to the increase in computational power and development of simulation methods. Starting with the seminal work of ref. 2, a growing number of research papers are dedicated to the determination of bending rigidity, [3][4][5][6][7][8][9][10][11][12] Gaussian modulus, 13 and surface tension, 7,8,11,14,15 from molecular simulation of all-atom systems to coarse-grained models.
A number of questions have been raised regarding the definition of membrane surface tension. When a thermally fluctuating flexible membrane is stretched by a moderate lateral tension, the membrane extends its projected area by suppressing its transverse fluctuation. 16 This leads to multiple definitions of membrane tension: (i) ''internal tension'', s, conjugated to the real area A in the membrane Hamiltonian, (ii) mechanical ''frame tension'', t, conjugated to the projected area A p , 17 and (iii) ''fluctuation tension'', r, associated to the lowest-order wavevector dependence of the inverse fluctuation spectrum. The latter two are directly experimentally measurable. 18 Some discussions on the difference between these surface tensions have been lasting for more than a decade. First of all, it is well accepted that s and t are intrinsically different, because s is a microscopic parameter belonging to the effective Hamiltonian while t is a macroscopic observable integrating the local stress tensor over all the membrane fluctuations. Based on free-energy calculations and stress tensor averages, the difference between s and t has been estimated within the Gaussian approximation. 14,19,20 Note that the validity of this estimation was questioned by Schmid. 21 It is also well accepted that s and r are different, because r is actually the renormalised version of s (at the Gaussian level they are equal). The relationship between r and t has been much debated. Several authors have argued theoretically, 5,15,22,23 or observed in numerical simulations, [2][3][4]7,8,15,21,[24][25][26][27][28] that r matches the mechanical frame tension t. Other authors have observed a difference between r and t, either in numerical simulations 11,14,20,29 or experiments. 18 These questions should be carefully revisited with the help of extensive, large-scale numerical simulations. A complication has to do with the equivalence of ensembles in the thermodynamic sense, as in the simulations found in the literature, either the areas A and A p or their conjugated tensions s and t are fixed. Different ensembles are only equivalent in the thermodynamic limit of infinite membranes, A -N at fixed A/N, with N the total number of degrees of freedom. It is therefore very important to check the convergence of the numerical results with system size. In another recent simulation, 30 r was claimed to deviate from t for t o 0 and to coincide with t for t 4 0. It is thus interesting to study how r behaves when t assumes positive values very close to zero, for which the system undergoes large fluctuations. Let us stress that investigating several orders of magnitude of t is important, because membrane tension can experimentally vary over more than five decades. 31 In this paper, we investigate the differences between the three surface tensions, s, t and r, over several decades of t, by means of Monte Carlo simulations of a lattice membrane system in the Monge gauge. In the lattice model, the internal tension s and the frame tension t can be controlled and short range fluctuations generated by molecular or particle protrusions, which make accurate estimations of r difficult, are not present. However, a measure correction is required to cure the excess entropy due to the unidirectional lattice-site motions implied by the Monge gauge. We also employ exact nonlinear expressions to calculate the membrane's curvature and area. We find that s differs from t at small positive frame tensions, while they match at moderate and high frame tensions. Our results for s À t agree well with the prediction of ref. 20 based on a Gaussian approximation, although our simulation investigates regimes where the Gaussian approximation is far from being justified. We find also that r and t match, within the accuracy of our simulation, over more than five decades, even in the limit of very small frame tensions. These relationships remain in the thermodynamic limit and seem not to depend on the ensemble in which the real area A is controlled (see Appendix A).
Our paper is organised as follows. In Section 3, we detail our lattice membrane model and the Monte Carlo simulation framework that realises an equilibrium ensemble for fixed t and s. In this framework, the average of the real membrane area, hAi, is kept constant by adjusting s while t is varied. We also perform our simulation in the equilibrium ensemble for fixed t and A. This is discussed in Appendix A. In the thermodynamic limit of large systems, we thus simulate a membrane with fixed real area subjected to a variable frame tension. In Section 3, we numerically determine and we compare the three tensions s, t and r. We give our conclusions in Section 4.

Stress-controlled ensemble
Let us consider an incompressible, fluctuating lipid membrane with a fixed area A = A 0 , corresponding to a fixed number of lipids. We assume that it is attached to a deformable frame, of variable area A p , that exerts onto the membrane a fixed tension t. The Hamiltonian of the system is thus given by 1 where H = c 1 + c 2 represents the sum of the two principal curvatures of the membrane. The Gaussian bending modulus 1 k has been neglected, according to the Gauss-Bonnet theorem, which is correct for fixed angular boundary conditions or periodic boundary conditions along the frame. 32 We use the Monge representation, in which the membrane shape is describe by its height z = h(x,y) above the plane of the frame. Note that this parametrisation does not require the membrane deformation to be small. The only restriction is that overhangs are forbidden, since the function h is single valued. In this representation, the total curvature H is exactly given by 33 where subscripts represent spatial derivatives, e.g., h x = qh/qx. We work here in the non-linear regime, allowing for large deformations of the membrane, i.e., we do not use the custom linear approximation H C h xx + h yy (employed in the Gaussian framework). Note that in our simulations (see the details below) the maximum slope, max(|h x |), reaches 1.3 at vanishing frame tension, but does not exceed this value. This justifies both using the exact expression for H and neglecting overhangs. Close to the buckling transition, however, i.e., for large negative frame tensions, these assumptions may not be valid.
Because it is quite difficult to work with a fixed area A = A 0 , we change ensemble in order to control the conjugated internal tension s, instead of the membrane area A. In the thermodynamic limit, the two ensembles are expected to be equivalent. Indeed, as shown in Appendix A, we find quantitatively similar results in the ensemble in which the real area A is almost prescribed by means of a quadratic potential. In this stresscontrolled ensemble, the modified Hamiltonian weighting the membrane fluctuations is given by with hAi = A 0 systematically enforced by adjusting s.

Lattice model
In order to implement the membrane fluctuations numerically, we introduce a N x Â N y lattice on the plane (x,y) of the frame, with lattice spacing a. The height of the membrane surface is defined by h ij on each lattice site, with 1 r i r N x and 1 r j r N y . In the reminder of this paper, we take . Furthermore, we adopt periodic boundary condition. Thus the projected area is given by A p = Na 2 .
Importantly, because the number of lipids is fixed while the frame area A p is variable, we take N constant but a variable. In other words, in our simulation the lattice spacing is not fixed. We stress, however, that a always remains uniform over the lattice: it is the overall lattice spacing that changes. Therefore the projection over the reference plane of the vertices always maps onto a square grid of finite spacing and the vertices can never overlap. Note that a similar lattice membrane model (with fixed spacing and the Gaussian approximation) has been employed for Monte Carlo simulation in order to investigate inclusions effects. 34 Monte Carlo simulations are performed with the lattice heights {h ij } and the lattice spacing a as dynamical variables. Special care is required, however, for the sampling of the height variables. In the Monge gauge, each site can only fluctuate in the z direction, whereas in the physical gauge, because the membrane is fluid, the relevant fluctuations occur by local displacements along the membrane's normal. In order to avoid an artificial entropy increase of the states where the membrane is tilted, the naive measure Here, y ij is the angle between the normal vector of the site (i, j) and the z axis, and cos y ij ¼ 1 In this way, the quantum height displacement in the normal direction is always d 0 . Equivalently, we can keep the naive measure d N [h] and add to the Hamiltonian the following correction term: where T is the temperature and k B Boltzmann's constant. Note that this correction, that would be exact if the membrane was uniformly tilted, is only approximate if the tilt changes from site to site. In practice, it works well when the membrane slope is small, as in the present simulations. For highly tilted membranes, such as buckled membranes, the correction term, eqn (5), is not sufficient, as discussed in Appendix B. Note that the vertical motion of the sites is enhanced in the absence of this entropy correction, leading to formation of steep slopes in the membrane surface at small tensions. As stated before, when this correction is present, the maximum slope does not overcome 1.3. which also assures that there are no tendencies of overhangs or formation of steep slopes. The partition function in the stress-controlled ensemble introduced above is given by with H 0 = H* + H corr . Metropolis sampling is employed. 35 In each Monte Carlo step, N trial moves of the height of a random site are attempted. The change of a (hence A p ) is attempted every 5 Monte Carlo steps. The amplitudes of the h ij and A p changes are adjusted so that the rejection rate lies in the range 40-60%. Both A and A p change during the simulations. However, for a given Monte Carlo run at a specified value of t, the internal tension s is adjusted in order for the average of the membrane area hAi to have the specific value A 0 , fixed once for all. In our simulation, the real area A of the membrane is calculated as follows. The local area associated to a site P(i, j) is obtained as one half of the total area of the four triangles The area of each triangle can be calculated by using height differences between the lattice sites. Then, the total real area can be expressed as Simulation are carried out at least for 10 7 steps after the initial relaxation. The statistical errors are calculated from four or more independent runs.

Dimensionless units and orders of magnitudes
Since our simulated membrane has a prescribed average area A 0 and a fixed number N of degrees of freedom, we can picture each degree of freedom as a patch of lipids of fixed average area a 0 2 = A 0 /N. In the following we are going to nondimensionalise the lengths by a 0 (the size of a degree of freedom in the membrane) and the energies by k B T = 4 Â 10 À21 J (room temperature energy). In other words, we shall set a 0 = k B T = 1. Hence, all tensions will be given in units of k B T/a 0 2 . As previously discussed, we shall thus adjust s for each t in order to have systematically hAi = N, in dimensionless units. Let us estimate tension scale k B T/a 0 2 . In principle, we are free to choose the linear size a 0 of the lipid patches that constitute our degrees of freedom. It is convenient, however, to choose the smallest possible size, in order to allow for all possible fluctuations and to avoid using renormalised elastic constants. Owing to the level of coarse-graining of our simulation, in which the bilayer membrane of real thickness E5 nm is treated as a mathematical surface, it is natural to think of the lipid patches as regions of typical size, e.g. a 0 E 20 nm. Then, the natural unit of tension in our simulation is k B T/a 0 2 E 10 À5 J m À2 . The tension mechanically imposed on ordinary vesicle membranes, in experiments, are in the range of 10 À7 to 10 À2 J m À2 , where the smallest value corresponds to floppy vesicles and the largest value corresponds to the lytic tension. 36 We shall therefore have t span the range 10 À2 -10 3 in dimensionless units. Let us stress that this corresponds to five orders of magnitude. In our simulation, we choose a specific value of the bending rigidity corresponding to k = 10 k B T, i.e., k = 10 in dimensionless units. We shall consider membrane sizes corresponding to N = 400, 1600, 6400 and 25 600. With the value of a 0 given above, this corresponds to membrane of linear dimension ffiffiffiffiffi ffi A 0 p ¼ ffiffiffiffi N p a 0 ¼ 400 nm, 800 nm, 1.6 mm and 3.2 mm, respectively. These values are quite small if one thinks of macroscopic experiments, but they are reasonably large on the biological scale. Because of the change of ensemble that we have performed, we shall check the convergence of our results with increasing N.
In Fig. 1, a simulation snapshot for N = 6400 is shown in the tensionless state t = 0.

Frame tension and internal tension
From now, unless otherwise specified, all quantities will be given in dimensionless units (see Section 3.3). We investigate here the relationship between s and t. As explained above, we place ourselves in the (s,t) ensemble; however, for every imposed frame tension t, we determine the internal tension s that achieves hAi = N. Our simulated membrane thus effectively has a constant real area (up to fluctuations that become irrelevant in the thermodynamic limit). Simulations in the (A,t) ensemble, where A is nearly fixed by means of a steep quadratic potential, are discussed in Appendix A.
The determined values of s are plotted in Fig. 2a against the frame tension for different system sizes: N = 400, 1600, 6400 and 25 600. For any value of t we find s 4 t. Data for both positive and negative frame tensions are displayed together as a function of |t|: the upper branch corresponds to t 4 0 and the lower branch corresponds to t o 0. Note that for t 4 0 the difference between s and t is sizeable only for small values of the frame tension. For t = 0 we find that the residual internal tension is equal to s 0 C 0.500 AE 0.001, where the error bar takes into account the determinations using different values of N. In other words, the residual tension turns out to be fairly independent of system size.
When t decreases below some negative threshold t b o 0, the membrane buckles either into the x or the y direction: the average membrane shape undergoes a transition from flat to non-flat through a symmetry breaking (like a rod under compression). Note that the membrane height function remains single valued although it acquires a bimodal distribution. The absolute value |t b | of the buckling threshold decreases as the system size N gets larger. In the lower branch, the data points for |t| 4 |t b |, in the buckling state, are eliminated from the plots in Fig. 2a. It is notable that the condition s o 0 is not required for the buckling transition. We also simulated our lattice model in the (A p ,s) ensemble, i.e. with constant lattice spacing a and freely changing membrane area A. We found that the membrane buckles then at s = 0. These results suggest that s is not the mechanical force that drives membrane buckling.
In ref. 20, the difference between the tensions s and t was calculated within the Gaussian approximation, yielding . The lengths are given in units of a 0 , which corresponds to the size of a coarse-grained degree of freedom in the membrane. In this simulation, the equilibrium lattice spacing is a C 0.983a 0 , which corresponds to (A À A p )/A C 0.034. With typically a 0 E 20 nm and k B T = 4 Â 10 À21 J, the residual internal tension is about s C 5 Â 10 À6 J m À2 .  where L is the ultraviolet wavevector cutoff and all quantities have their normal dimensions. Within our simulation, this ultraviolet cutoff corresponds to the lattice spacing. We thus expect ; with a % 1: Whereas our simulation is fully nonlinear, the expression for s À t calculated within the Gaussian approximation fits the data quite well (Fig. 2b). Since the numerical value of hai changes by only 1.4% from t = 0 to 1000, we use for the fit the extreme value at high tension, a = 1. The best fit of the data gives then a = 1.12, very close to unity. Thermodynamic limit. In Fig. 3, we plot the two tensions s and t against the excess area (A À A p )/A for the various system sizes N = 400, 1600, 6400 and 25 600. Here, only positive values of t are displayed: the downwards divergence of the curves associated to the plain circles corresponds to t -0. The system size dependence is marked for small sizes, but the data appears to converge to a universal thermodynamic limit at large system sizes. Finite size effects are more or less irrelevant for N Z 6400.

Fluctuation tension
We now discuss the fluctuation tension r. It is derived from the spectrum density, assumed to take the standard form: 33 where h(q) is the Fourier transform of the height function, with q the wavevector. The fluctuation tension r and the bending rigidity k r , appearing in the fluctuation spectrum, are the renormalised counterparts of the Hamiltonian parameters s and k.
In estimating these quantities using simulation results, it is crucially important to consider large-scale fluctuations. However, because we always treat finite-size systems, the wavevector spectrum is limited, which results in systematic deviations. Previously, two of the present authors have addressed in detail the question of the estimation of the bending rigidity from simulation data. 8 For a planar membrane, the correct value of the bending rigidity can be determined uniquely by first estimating the rigidity k r via a fit of the inverse spectrum (rq 2 + k r q 4 )/(k B T) in the wavevector range 2p a 0 ffiffiffiffi N p À Á q o q cut , then by extrapolating the results to the limit q cut -0, where q cut is a varying upper fitting wavevector limit. Note that because the projected area A p fluctuates, there is some ambiguity regarding the definition of the quantified wavevectors. As discussed in Appendix C, our results are not sensitive to those details. This procedure has been carried out for the present system. We find that the value of k r coincides with the bare value k, up to our numerical precision. This result makes sense because renormalisation group calculations predict up to a numerical factor of order unity 33 (here all quantities have their normal dimensions). This prediction gives Dk/k E 0.03 for N = 6400, which is less than the error bars, and thus not detectable in our simulation.
The correct value of the fluctuation tension r can be estimated likewise in the limit q cut -0. For finite values of q cut , we find that r is overestimated but converges in the limit q cut -0 to a well defined value. The corresponding extrapolation is performed by using a quadratic function of q cut . The fitting was performed in the range 0.015 r (q cut /p) 2 o 0.15. Fig. 4 shows the values of the fluctuation tension r. We find that r is fairly close to t for all the investigated values of t, as indicated in Fig. 4a. In Fig. 4b, the difference r À t is plotted against t. For large system sizes, N = 6400 and 25 600, we observe a small deviation between r and t; however, as discussed in Appendix B, we believe that this deviation is smaller than a possible systematic error arising from the measure correction.

Comparison with the Gaussian approximation
In the present work, all the nonlinearities of the problem have been taken into account: (i) the nonlinear expression of the membrane elementary area dA (ii) the nonlinear expression of the membrane bending energy density H 2 dA, and (iii) the nonlinear measure correction H corr dealing with the excess entropy of the Monge gauge. The Gaussian approximation, which is frequently used in basic calculations, consists in replacing the bending energy density H 2 dA by H L 2 dA p , where H L = h xx + h yy is the Laplacian approximation of the mean-curvature and dA p the projected elementary area, and in neglecting the entropy measure correction H corr . Indeed, at quadratic order, the measure correction Fig. 4 (a) Plot of t, r, and s as a function of excess area (A À A p )/A for N = 6400, which roughly corresponds to the thermodynamic limit. The data for t and s are the same as in Fig. 3, and r is plotted with the error bars. As for the statistical errors in (A À A p )/A, they are small, within the symbol marks. (b) The difference between r and t, normalised by s, is plotted against t for N = 1600, 6400, and 25 600. The convergence of the data with increasing N indicates that the thermodynamic limit is reached for N = 6400.  (1), and we keep the entropy correction term H corr . The corresponding results are denoted by the subscript EL in Fig. 5. Then we also remove the entropy correction term. The corresponding calculations are denoted by the subscript L.
We find that the internal tension s is not changed by these modifications. However, r EL and r L are significantly different from r, as shown in Fig. 5. Using the Laplacian approximation while keeping the entropy correction term results in a fluctuation tension r EL about three times larger than s. Further removing the entropy correction yields a tension r L very close to s, in agreement with the well-known coincidence of r and s at the Gaussian level. Therefore, the effects of the measure correction and of the nonlinear curvature are opposite to each other: the former raises the value of r, the latter lowers it. Both of these effects need to be taken into account for a faithful description of the membrane fluctuations beyond the level of the Gaussian approximation. The estimated r deviates from the true value if either of them is missing.

Conclusions
In this paper, we have examined the various surface tensions of membranes: frame tension t, internal tension s, and renormalised ''fluctuation'' tension r. We have compared them quantitatively using a Monte Carlo simulation with control over both t and s, the latter being slaved to the former in order to keep the average membrane area hAi constant. We have also validated our results in the conjugated ensemble where the real membrane area A is fixed (Appendix A). We have investigated large systems in order to reach the thermodynamic limit. Our model being a lattice model instead of a particle-based model, protrusion effects and other artificial elements are eliminated. Gaussian approximation artefacts are also excluded because curvatures and areas are computed with exact formulas (although we employ the Monge gauge) and because we have corrected the excess entropy associated with the measure of the Monge gauge. Moreover, because our simulated membrane does not exhibit rupture at extremely large t, all tensions from vanishingly small to very high can be investigated.
There have been two findings in our results. The first is that the residual internal tension s 0 = s À t remains finite at large N (large systems). While it is well known that s and t are different, we have confirmed that the theoretical (albeit Gaussian) estimation of ref. 20 predicts correctly and quantitatively the difference s À t from vanishingly small frame tensions to very large frame tensions.
The renormalised ''fluctuation'' tension r has also been investigated, and compared with the other two tensions. We conclude r C t within an accuracy of 0.1k B T/a 2 (estimated C 10 À6 J m À2 ), which is consistent with the proposition of ref. 22.
Note finally that the shape of the buckled membrane, described by elliptic functions, breaks the symmetry between the x and y directions and yields an anisotropy in the frame tension. 9 How exactly the thermal fluctuations of the membrane modify the mechanical buckling transition should be investigated in further studies.

A Simulation with nearly constant real area
In the body of this paper, the real area A of the membrane is not fixed, but its average value hAi is controlled by the parameter s (much as the average number of particles in a grand-canonical ensemble is controlled by the chemical potential). In principle, we have checked the validity our results by studying the thermodynamic limit of increasing system sizes. It is nonetheless interesting to test our results by working in the ensemble in which the area A is prescribed. This is the purpose of the present Appendix.
While it is very difficult to prescribe exactly the membrane area, because it is difficult to identify the sampling condition that satisfies this constraint, one can easily almost prescribe the area by using a quadratic potential of tunable strength K A . We thus worked with the Hamiltonian We have performed the corresponding Monte Carlo simulations exactly in the same manner as in the body of the paper. Fig. 6 shows the result of these additional simulations, for K A = 1, 10 and 100, with N = 6400, and compares them with the results of Fig. 4. With increasing K A , the standard deviation of the real area decreases: 0.01%, 0.001%, and 0.0003%, for K A = 1, 10, and 100, respectively.

B Entropy correction for membrane tilt
To check the reliability of the measure correction H corr , eqn (5), we performed simulations in which the entire membrane is tilted by an angle y a . The periodic boundary condition along the x direction is modified to h(N x , j) = h(0, j) + N x a sin y a . The x side Fig. 5 Excess fluctuation tensions r, r EL and r L , measured with respect to t and normalized by s, plotted against t, for the nonlinear case (r), and for the Gaussian curvature approximation with the entropy correction (r EL ) or without it (r L ). Here, N = 6400. For comparison, a plot of (s À t)/s is shown.
length of the lattice is changed to a cos y a in order to maintain a square grid in the tilted projected plane. The height spectrum is calculated for the tilted plane, of normal vector (Àsin y a , 0, cos y a ). The fluctuation tension r increases with increasing y a (see Fig. 7), while s is independent of y a . This implies that the excess entropy associated with the membrane tilt in the Monge gauge is not completely removed. However, the deviation is small for small y a , and the overestimation of r is less than 0.1 for y a o 0.1p. In our simulation, at t = 0 we obtain hcos y ij i = 0.98, i.e., the mean angle is 0.066p. We conclude that the obtained values of r are reliable within the accuracy of 0.1 (in dimensionless units).

C Definition of the spectrum wavevectors
In Section 3.2, the renormalised tension r is calculated by fitting the height spectrum h|h(q)| 2 i. We use a sequence of instantaneous height configurations, and we perform averages after the simulation finishes. Because the projected area A p changes during the simulations, the wavevectors on the reciprocal lattice are not fixed. They are given by where n and m are integers. All these data are mapped to a onedimensional wavevector q = |q|. In plotting h|h(q)| 2 i, the data of |h(q)| 2 are averaged over close values of q (the width of q bins is bp ffiffiffiffi N p . Here, b = 1.2 is taken) and over all the height configurations.
We have confirmed the robustness of our estimation of r, by using a different set of wavevectors, defined by ðn; mÞ: For each of integers n, m, the values of h|h(q 0 )| 2 i is averaged over the sequence of instantaneous height configurations. The results obtained by these two methods do not essentially differ from one another: the differences in the value of r are smaller than 0.025 k B T/a 0 2 for all t, which is typically smaller than the error bars. Fig. 6 The deviation of the fluctuation tension r from the frame tension t is plotted in the form (r À t)/s. The red points are pasted from Fig. 4b for comparison. The results for three values of K A (=1, 10, and 100) are shown by using symbols and error bars of different colours. For all the plotted data, the system size is N = 6400. Although we normalise the data by s, no estimate is available in the present simulation; therefore, we use the values of Fig. 4a as substitutes. Fig. 7 Dependence of the fluctuation tension r for a membrane whose projected plane is entirely tilted by the angle y a , for t = 0 and t = 1, at N = 6400.