The free energy of nanopores in tense membranes

Membrane nanopores are central players for a range of important cellular membrane remodeling processes as well as membrane rupture. Understanding pore formation in tense membranes requires comprehension of the molecular mechanism of pore formation and the associated free energy change as a function of the membrane tension. Here we propose a scheme to calculate the free energy change associated with the formation of a nanometer sized pore in molecular dynamics simulations as a function of membrane tension, which requires the calculation of only one computationally expensive potential of mean force. We show that membrane elastic theory can be used to estimate the pore formation free energy at different tension values from the free energy change in a relaxed membrane and the area expansion curves of the membranes. We have computed the pore formation free energy for a dipalmitoyl-phosphatidylcholine (DPPC) membrane at two different lateral pressure values, 1 bar and 40 bar, by calculating the potential of mean force acting on the head group of a single lipid molecule. Unrestrained simulations of the closing process confirm that the intermediate states along this reaction coordinate are reasonable and show that hydrophilic indentations spanning half the bilayer connected by a hydrophobic pore segment represent the corresponding high energy transition state. A comparison of the stability of simulated membranes to experiment at high loading rates show that, contrary to expectation, pores form too easily in small simulated membrane patches. This discrepancy originates from a combination of the absence of ions in the simulations and the small membrane size.


Introduction
Lipid bilayer membranes are essential for maintaining the structural integrity of cells, for the spatial organization of the cell interior and to serve as a scaffold for many proteins. 1To fulfill these roles, lipid membranes combine structural stability with flexibility allowing them to adjust to external forces.Nonetheless, lipid bilayers can form transient pores or rupture under mechanical, electrical or chemical stress.
Such pore formation plays an important role in many cellular processes that require membrane remodeling, ranging from protein insertion 2,3 to membrane fusion [4][5][6][7] and fission. 8,9eyond these cellular processes, controlled pore formation is also important for biomedical applications such as drug delivery.For these reasons, vast experimental as well as theoretical and computational efforts have been invested in understanding pore formation in lipid membranes.
In experiments, unilamellar lipid vesicles often serve as model systems to help elucidate the mechanical properties of membranes.][23] Many theories of pore formation are based on classical nucleation theory, in which the Gibbs free energy change DG pore of opening a circular pore with radius r in a membrane at constant tension s is expressed as 24,25 DG pore = 2pLr À psr 2 , ( where L is the edge energy of the pore. In this model, a pore is never stable but either disappears or expands.Experiments, however, show well-defined levels for transmembrane conductance indicating the existence of metastable nanopores. 26Stabilized nanopores could also explain the dependence of the rupture tensions of lipid membranes under tension on the rate at which the tension was built up 27 by introducing a stabilized nanopore state to the energy function (1) close to r = 1 nm.
It appears reasonable that the continuum description does not hold for pore radii close to r = 1 nm, as a general assumption is that the length-scale of the described features is large compared to the molecular scale.For 1 nm size pores, this no longer applies and the closing of a hydrophilic pore will require considerable rearrangement of the lipid molecules at the pore edge, where the hydrated headgroups covering the pore edge have to be expelled from the membrane interior.From that consideration it has been proposed to add an additional term describing the nucleation free energy to eqn (1), [28][29][30] and its application has been extended to the cases of constant area 29,30 and fluctuating membranes. 31o observe spontaneous pore formation in simulation studies, much larger membrane tensions are required, 32 at which vesicles are highly unstable in experiments.This is due to the much smaller length and time scales accessible to the simulations, but makes a direct comparison of membrane stability in experiments and simulation difficult.3][44] Both for high tension 32 and electroporation, 36,40,45 a hydrophobic defect is observed before the formation of a hydrophilic pore.In addition, some studies have described the formation of hydrophilic 'dimples', before the appearance of the hydrophobic pore. 39,40Such 'dimples' may facilitate insertion of water molecules, 46 and similar defects are likely to be involved in the translocation of ions across the membrane. 479][50][51] Although in good agreement with the continuum description, many CG models have difficulties in reproducing the behavior observed at the atomistic scale for pores with a radius close to r = 1 nm. 52lthough the lifetimes of simulated membrane nanopores suggest that they are stabilized against closure, 53 various simulation studies with both atomistic MD 54,55 and CG approaches 29,30 that used different reaction coordinates at both resolutions have been unable to measure the corresponding free energy barrier stabilizing the nanopores.The two reaction coordinates chosen for these studies were the local density of the membrane 29,30,54 or the distance of one lipid phosphate group from the membrane center. 55In both cases, the reaction coordinate is not directly linked to the existence of a pore but spontaneous pore opening is observed from a certain value of the reaction coordinate.Presumably, a different reaction coordinate, directly linked to the pore opening/closing process, may be required to observe a barrier stabilizing the pore.Pores created via a density fluctuation in the membrane required approximately 140 kJ mol À1 to form, 54 whereas the free energy change for pores created by manipulating a lipid headgroup was approximately 80 kJ mol À1 . 55With a similar approach the effects of different lipids, 56 cholesterol 57 and the presence of peptides 44 on lipid flip-flop and pore formation were investigated.
Due to the large computational expense involved, methods to calculate the potential of mean force are applied to only one or a few fixed tension values.However, for a reasonable comparison with experiments, such as the tension spectroscopy, 27 the free energy of pore formation is required for a whole range of membrane tensions.To achieve this by calculating a free energy profile for each tension would be very computationally expensive.For a better comparison at a reasonable computational cost, we introduce an approach to estimate the free energy of a nanometer size pore in a simulated membrane patch as a function of the lateral tension, which combines the use of both all-atom MD simulations and elastic theory.Atomistic simulations are applied to calculate the free energy of pore formation at two different tension values, and to measure the elastic properties of the membrane patch, both with and without a pore.These are then used together with the continuum model equations to estimate the pore formation free energy at different tensions, as described in the next section.Further relevant properties observed in the atomistic simulations include the pathway, dynamics, and transition state for closing of the nanopores.The results of our simulations have then been applied together with the three state model used by Evans et al. 27 for a comparison with experimental results.

Theory
To circumvent the high computational cost of estimating the free energy required to open a nanopore as a function of the membrane tension S by performing many PMF calculations, we propose the scheme illustrated in Fig. 1 and described below.The process of creating a pore in a bilayer at a particular tension S is divided into the three consecutive processes of (a) relaxing the membrane tension s from s = S to s = 0, (b) creating a pore in the relaxed membrane, and (c) stretching the membrane containing a pore from s = 0 back to s = S.The corresponding free energy changes DG n (S), DG*, and DG p (S) for the three processes are estimated separately as follows.
The free energy changes DG n and DG p corresponding to the two processes (a) and (c), respectively, can be estimated using the area A(s) as a function of the membrane tension.The dependence of A(s) on s can be measured directly in the simulations.The free energy change DG* corresponding to part (b) can be determined with the help of a PMF calculation restraining the lipid in the center of the bilayer as described in ref. 36 and 56 and in the following section.
In general, the change in the system's free energy dG is given by the expression where, as usual, V, P, S, T, m i and N i are the volume, pressure, entropy, temperature, chemical potential, and number of molecules, respectively.Similarly, the membrane excess free energy dG m , defined in analogy with Gibbs surface thermodynamics, 58 is given by where S m and N m i are the surface excess entropy and molecule number, respectively.DG m is the free energy change typically considered in elastic theory calculations.T and N i in the simulation box are constant.Furthermore, the pressure is constant at P = 1 bar and isotropic in the solvent phase, which was verified numerically using the local pressure version of gromacs. 59he change in the system's free energy upon reduction of the lateral tension in the bilayer patch from S to 0 is therefore given by For the membrane excess free energy change, the values of the molecules attributed to the membrane N m i = N i À c i V, where c i is the number density of species i in the bulk solvent phase, have to be taken into consideration.The concentration of lipids in the bulk solvent phase is zero, so that for the lipids N m i is the total number of lipids in the system and therefore constant.For water, on the other hand, N m i is not necessarily constant.Therefore, N m i = N i À c i V has been monitored numerically and both c i and V are found to remain constant for the range of lateral pressures simulated here, so that N m i can be considered constant and the change in the membrane excess free energy is given by DG n;m ðSÞ ¼ À ð 0 S AðsÞds: (5) While the pore is metastable and the membrane does not rupture, the same relations apply to stretching a patch of membrane containing a small nanopore from s = 0 to s =S: and Using eqn ( 4) to ( 7) the free energy changes corresponding to steps (a) and (c) can be estimated if A(s) and A p (s) are known.The total free energy of a nanopore in a tense membrane with s = S can be estimated from or In practice, A(s) and A p (s) are determined by a series of simulations at different values of s.This protocol therefore requires only two simulations for each tension value compared to the 50 umbrella windows required for each PMF profile.Therefore our protocol reduces the computational cost by about 25 times.The computational details are described in the next section.

Computational methods
All atom MD simulations were used to simulate a number of different DPPC bilayer patches.These simulated systems were: (i) unperturbed bilayer patches at different lateral pressures in the range of 0 to À45 bar, for which a small pore was found to be stable against rupture; 53 (ii) bilayers containing a nanometer sized pore, stabilized by one lipid headgroup restrained in the bilayer center, at the same lateral pressure values; (iii) tensionfree membrane patches, with a lipid phosphate group restrained at the bilayer center to evaluate pore opening and closing rates in the biased system; (iv) membrane patches at different lateral pressures initially containing a nanopore, but without lipid restraints for stabilization; and (v) umbrella sampling simulations for both a relaxed bilayer and a bilayer at a lateral pressure P lat = À40 bar, to compute the potential of mean force (PMF) corresponding to the displacements of a DPPC phosphate group from its equilibrium position to the membrane center.
All simulations were performed in the NPT ensemble using the GROMACS molecular dynamics package, 60 version 3.3.1.The temperature was maintained at 323 K by separately coupling lipids and water to a heat bath via the Berendsen scheme 61 with a relaxation time of 0.1 ps.The average pressure was controlled using the Berendsen scheme 61 with a relaxation time of 1.0 ps and the box dimensions parallel and normal to the bilayer were scaled independently.Periodic boundary conditions were employed and the overall center of mass motion of the system was removed.
Lipids and water molecules were described using parameters from Berger et al. 62 and the simple point charge (SPC) model [63][64][65] respectively.A 1 nm cutoff for Lenard-Jones parameters was employed, and full electrostatic interactions were evaluated using the particle mesh Ewald technique 66 with a cutoff distance of 1 nm in direct space, a 0.12 nm grid spacing and a 4th order polynomial.Covalent bonds in the lipids were constrained using LINCS 67 and water molecules were kept rigid with SETTLE 68 allowing the use of a 4 fs timestep.
The details of the simulated systems are as follows.Each simulated system contained 64 DPPC lipids surrounded by 3846 SPC water molecules.The initial configurations of systems (ii) and (iii) containing a nanopore were generated by restraining the phosphate group of one of the lipids at the bilayer center using a harmonic potential with a force constant of 5000 kJ mol À1 nm À2 until a small pore formed.The pores were stabilized during the entire simulation time by the same restraining potential.For the tense bilayers in (ii) this ensured that the pores remained open.For (iii), on the other hand, pores can open and close in the presence of the restraining potential.Therefore 20 simulations were started both from initial states with an open and with a closed pore.The simulations were run for up to 400 ns, until a transition to the closed and open pore, respectively, was observed.Starting structures for the closing simulations (iv) were obtained from 10 equally spaced snapshots from the simulations (ii) at each lateral pressure.Pore closing rates and pathways at different tensions were obtained from simulations without restraining potential.
In the PMF calculations (v), the position z relative to the center of the membrane of the restrained phosphate group was chosen as the reaction coordinate.For each PMF, 51 umbrella windows between z = 0 nm and z = 2.5 nm were used.Starting configurations for each window were created with 1 ns simulations using an umbrella potential with a lower force constant of 500 kJ mol À1 nm À2 and subsequently equilibrated for 10 ns with a final force constant of 5000 kJ mol À1 nm À2 .Data were collected from 100 ns simulations.The weighted histogram analysis method 69 was used to construct the PMF from the biased distributions.Standard errors for the PMFs were obtained by dividing the trajectories into four blocks and calculating the standard deviation between the corresponding profiles, which were aligned at the lipid's equilibrium position.

Formation and closing of membrane pores
Previous simulation studies 55 have shown that when a lipid head group is restrained in the hydrophobic center of a bilayer, a small hydrophilic pore may form.The corresponding free energy change for DPPC has been estimated to be 80 kJ mol À1 from the corresponding potential of mean force. 55Similarly, in the simulations described here, a small hydrophilic pore shown in Fig. 2(a) forms when the phosphate group is restrained at z = 0 nm from the bilayer center.In the tension-free membrane, this pore is found to coexist with a half-pore state, i.e. a pore spanning only one of the monolayers as depicted in Fig. 2(b).The additional simulations of pore opening and closing with the potential at z = 0 show that both states are approximately equally populated, with mean lifetimes of approximately 100 ns.
The PMF profile shown in Fig. 2(d) steeply increases when the restrained lipid headgroup is displaced from its equilibrium position at z E 1.70 nm to the center of the bilayer at z = 0.A close-up view of the PMF in the vicinity of z = 0, however, shows that large differences of the order of 6 kJ mol À1 exist between the profiles from different trajectory segments.These large error bars are directly related to the two states observed for the system, i.e. the full-pore and half-pore state.That these will have a considerable effect on the PMF value is apparent from Fig. 2(c), where the transition between the full pore and a half pore is clearly visible in the distinct jump of the position of the restrained head group relative to the umbrella potential from 0 AE 0.001 nm to 0.016 AE 0.0005 nm when the pore closes.When the pore is present, there is no opposing force from the bilayer on the lipid and it can remain in the center of the harmonic potential.Similar large errors close to z = 0 are also shown for PMF profiles in other simulation studies where a nanopore may form close to z = 0. 44,52 A consequence of the two states is that, close to z = 0, the probability of forming a membrane spanning pore in the simulations becomes larger than zero and the PMF represents a mixture of the two states.Their relative population at each value of z depends on the relative biased free energies, G p biassed (z) and G hp biassed (z) of the biased system with a pore and a half-pore, respectively.As transition times between the two states are on the order of 100 ns at z = 0 and even longer for slightly larger values of z, the equilibrium distribution between the two states will not be sampled in simulations.Furthermore, the PMF at z = 0 will not correspond to the pore free energy, but to that of the mixture of states.
Instead, for a more accurate estimate of the pore free energy, here we calculate the PMF in the presence of a half-pore reaching z = 0 by including only trajectory parts which do not contain a full pore (Fig. 2(e)).We then estimate the free energy difference between a pore and a half-pore in the presence of the harmonic potential U(z = 0) by sampling many transitions between the two states.The corresponding average transition times for both states and their distributions are very similar.Taking into account also the mean position of the harmonic potential of the two states, this suggests that G p biassed (z = 0) E G hp biassed (z = 0.016 nm) + U(z = 0.016 nm).From the PMF (Fig. 2(e)), G hp biassed (0.016 nm) = 81.4AE 1.2 kJ mol À1 .This leads to an estimate of the free energy of the nanopore of 82.1 AE 1.2 kJ mol À1 for the bilayer.This value agrees with previous simulations 55 which only observed a full pore state.As our reaction coordinate is not directly linked to pore formation, but depends on the reaction coordinate z and the probability of overcoming an energy barrier close to z = 0, it is important to ensure that our reaction coordinate z leads to a realistic pore formation pathway.It is not possible to observe the formation of nanopores in unrestrained simulations of bilayers at reasonably low tensions, but simulations of the opposite process, i.e. nanopore closing, can give insight into intermediate steps on the pathway.Such simulations can reveal whether the monolayer spanning half-pore represents a reasonable high energy intermediate state or whether under unconstrained conditions, pore formation would proceed via a different pathway.
Trajectory snapshots from the two very similar closing pathways observed in such pore closing simulations are shown in Fig. 3.These include either a half-pore intermediate state similar to that in the restrained simulations (Fig. 3a), or a similar transition state (Fig. 3b) with smaller hydrophilic indentations present in both monolayers instead of the half-pore in one monolayer, similar to the 'dimples' observed in some simulations of electroporation. 39,40In both cases the hydrophilic transition states are preceded by a similar conformation, in which the hydrophilic pore segments are connected by a narrow water column also shown in Fig. 3.This water column dissolves quickly and is likely to represent the energy barrier state.A full hydrophobic pore spanning the entire bilayer is not observed in the 96 closing simulations.The similarity between the intermediate states in the unrestrained closing simulations and the opening pathway in the umbrella simulations increases the confidence that z is a reasonable reaction coordinate for nanopore formation.It would be interesting to see if a reaction coordinate controlling the presence and length of the hydrophobic pore segment can reveal an energy barrier for nanopore formation and closure.

Free energy estimate for nanopores in tense membranes
For a comparison of membrane stability with experimental data, the corresponding free energies of nanopore formation as a function of the membrane tension are required.However, repeating the above calculations for a number of different tension values would be extremely time consuming.Therefore here we propose instead to divide the poration of a membrane at tension S into the three subprocesses of relaxing the intact membrane, creating a nanopore and then stretching the porated membrane back to S, as described in detail above (the Theory section).
To apply this scheme, we need the membrane area as a function of the tension both for a membrane without pore, A(s), and with pore, A p (s).These relations are obtained from a number of simulations of membrane patches -both with and without nanopores -at different lateral pressures.A(s) and A p (s), shown in Fig. 4, are fitted well by a second order polynomial function.The stretching modulus K A can be determined from the gradient of s(A) close to the tension-free membrane.The values K A = 283 mN m À1 and K A = 281 mN m À1 are found for the membrane with and without nanopore, respectively, which is in reasonable agreement with the experimental value of approximately 250 mN m À1 . 11ubstituting the fits for A(s) and A p (s) into eqn ( 4)-( 9) and using the above estimate of DG* = 82.1 kJ mol À1 for the tension free membrane, the free energy change of the simulation box, DG total (S), and of the membrane, DG m total (S), can be estimated as a function of tension as shown in Fig. 5(a).As expected, the membrane free energy change upon nanopore formation, DG m total (S), decreases with the tension.Due to the increased area  of the porated membrane and subsequent deformation of the simulation box, the free energy of the simulation box on the other hand increases.
A simple estimate for the radius of the nanopores in the simulations is obtained from the difference between the area of the membrane patches with and without a pore as pr(s) 2 = A p (s) À A(s).Similar to the membrane area, r(s) is fit well by a second order polynomial, shown together with r(s) in Fig. 4b.If the tension S is the same for the integrals ( 5) and ( 7), AðsÞds can be written as DG m total ðSÞ À DG Ã ¼ Àp Ð S 0 rðsÞ 2 ds.The expression obtained using the fit for r(s) in this integral is shown in Fig. 5(b) together with DG m total (S) À DG*.To test the accuracy of the free energy values calculated in this scheme, we have calculated a second PMF at a large lateral pressure of P lat = À40 bar.The corresponding free energy difference for pore formation DG total is estimated similarly as described above for the tension-free bilayer patch, and is found to be 65.7 AE 3 kJ mol À1 , i.e. approximately 16 kJ mol À1 lower than that for the relaxed bilayer patch.
To compare this value directly with the elastic theory estimate, it is important to note that in the simulations described here, the lateral pressure is controlled by semi-isotropic pressure coupling treating the z direction and the x-y plane independently.In the PMF calculation, the membrane area increases and the box deforms when the lipid is constrained at different z positions.As a consequence, the membrane tension changes continuously with z in the PMF from S = 33.5 mN m À1 at z = 1.6 nm to S = 31.3mN m À1 at z = 0 nm.Therefore, to compare the free energy difference found from the PMF with the integration results, these different values are used as the start and end values of S in the integrals (4) and ( 6) for the membrane without and with a pore, respectively.As the tension of the porated bilayer (corresponding to z = 0) is lower than that of the intact bilayer, the positive contribution DG p (S = 31.3mN m À1 ) becomes smaller compared to the negative contribution at DG n (S = 33.5 mN m À1 ).Using these different tension values as the integration boundaries, in fact, DG total À DG* changes sign and the system's free energy change upon pore formation is found to decrease to DG total = DG* + DG n (33.5 mN m À1 ) + DG p (31.3 mN m À1 ) = 68 AE 1.2 kJ mol À1 .This is in excellent agreement with the value 65.7 AE 3 kJ mol À1 obtained from the PMF calculations.

Comparison with the continuum membrane model
The good agreement between the PMF and the integration results despite the very large tension difference demonstrates that the integration method produces valid results over a large range of tensions.We can therefore use the results to obtain estimates for the density of nanopores per lipid, r p (s), as a function of membrane tension s using r p (s) = exp(ÀDG m total (s)/k B T).The corresponding probability of finding a nanopore in a membrane of a certain size is found by scaling with the membrane area.The tension values for which the predicted probabilities of finding a pore reach values of approximately 10% are 60 mN m À1 for a bilayer with the area of simulation box and 15 mN m À1 for a vesicle with a diameter of 20 mm.These values correlate fairly well with the order of magnitude of rupture tensions observed for the two systems, despite the very different timescales.

Comparison with tension spectroscopy experiments
The hypothesis that membrane nanopores represent a stabilized intermediate state could successfully explain the dependence of the membrane rupture tension on the tension loading rate R = ds/dt observed by Evans et al. 27 In their three state Markov model, membrane rupture proceeds in two steps.First, a membrane nanopore forms, which can either reseal or expand to a large membrane pore.The master equations describing the population of the three states, (i) an intact membrane, S 0 , (ii) a membrane containing a nanopore S * , and (iii) a ruptured membrane S h , are : S 0 (t) = Àn 0-* (s)S 0 (t) + n * -0 (s)S * (t) (10) : S * (t) = À(n * -0 (s) + n * -h (s))S * (t) + n 0-* (s)S 0 (t) (11)  : S h (t) = n * -h (s)S * (t) (12)   where, in principle, the transition rates n 0-* , n * -0 and n * -h could depend on the membrane tension.Experimental data showed two different regimes for pore formation at different loading rates.For fast loading, rupture is limited by the formation of a nanopore (a defect limited regime), whereas for low loading rates rupture is limited by the expansion of the pore (a caviation limited regime).
From our simulations, a number of properties defining the system can be estimated.The nanopore closing rates n * -0 and their possible dependence on the membrane tension can be obtained from unrestrained simulations starting from an initial state containing a nanopore.The pore opening rate per lipid, n 0-* , can then be estimated from a detailed balance condition.These two rates are sufficient to describe the defect limited regime.In addition, an approximate crossover tension between the two regimes can be estimated from the point at which expansion of the nanopore begins to compete with closure, indicating that the barrier for expansion E c = pL 2 /s is reduced to a similar magnitude as that for pore closure.This journal is © the Owner Societies 2014 The pore closure times in unrestrained simulations for a range of different lateral pressures between 1 and À45 bar were studied using 10 independent simulations at each tension value.Surprisingly we found that pores typically closed within 10-20 ns, independent of the membrane tension.The tension independent timescale agrees with the hypothesis used in the model in ref. 27 and simulations of bilayer self-assembly. 33owever, longer pore lifetimes were observed in a systematic simulation study for tense membranes. 53A possible explanation for this discrepancy is the small bilayer size used here.Test simulations of a larger membrane patch with 128 lipids under otherwise unchanged conditions showed mean closing times of 47 ns for a tension-free bilayer as also described in ref. 70.In fact, simulations of electroporation have suggested that up to 140 lipids can be affected by a nanopore. 40Note that for the pore closure simulations, the type of pressure coupling in the x-y plane again comes into play.If semi-isotropic pressure coupling is used for tense membranes, the membrane tension increases slightly upon pore closure, making the latter less favorable.At lower tensions, this effect is small; nonetheless, surface tension coupling was used in the closing simulations.
To compare the stability of the simulated membranes to the experimental results 27 for the defect limited regime, we calculate the maximum in the probability distribution for the rupture tensions, P(s), derived in ref. 27.In this regime, the membrane will rupture as soon as a nanopore defect is formed, so that the rupture time can be approximated by the average time for nanopore formation.Using the constant pore closing rate of n * -0 = 1/10 ns À1 and the fit for DG m total (s) À DG*, the nanopore formation rate per lipid n 0-* (s) is obtained from . This opening rate depends on the membrane tension and has to be scaled by the number of lipids A/A lipid for a given membrane size A and area per lipid A lipid .The condition dP(s)/ds = 0 for the maximum in the probability distribution leads to the expression (A/A lipid )n 0-* (s*) = R(qln(n 0-* )/qs) s* (13)   which we solved numerically for the rupture tension s* with different loading rates.The rupture tensions as a function of the loading rate R for a simulation box and for a 20 mm vesicle predicted this way are shown in Fig. 6(a) and (b), respectively.In both the plots, the range of loading rates shown was adjusted to obtain the order of magnitude for the rupture tensions observed experimentally for various PC lipids. 27The loading rates which achieve this for the small membrane patch in the simulation box are approximately an order of magnitude lower than the experimental ones, whereas those predicted for a 20 mm vesicle are five orders of magnitude larger.This unexpectedly indicates that, unlike often discussed, the simulated membrane patches are not too stable, at least with respect to nanopore formation, but they rather form pore defects several orders of magnitude too easily.
For lower loading rates, the caviation barrier remains high after a pore defect has formed, and pore expansion becomes the limiting step.The crossover between the two regimes is reached when rupture becomes more probable than closure of pores.
Here this is the case at 35 to 40 mN m À1 .This is in good agreement with the critical tension found by previous simulations of DPPC bilayers 53 and the predictions for the critical tension in the experimental study, 27 which range from 30 to 130 mN m À1 .
There are a number of different explanations that may contribute to the discrepancy between the loading rates predicted from the simulations and experiment.First, the physico-chemical conditions of the system, such as temperature, buffers and ion concentrations, can affect the membrane stability and edge tension significantly.AFM Force spectroscopy 71 and micropipette aspiration 72 experiments as well as previous MD simulations 73 have shown that glucose or salts like sodium chloride stabilize PC bilayers against indentation and pore formation.The MD simulations of small pores in the presence of sodium and chloride ions indicate that the pore edge energy increases by a factor of two relative to the pure bilayer system at a salt concentration of 0.2 M. 73 A second factor is the small size of the simulated membrane patch, and the associated faster pore closing times.6][57] The longer pore lifetimes observed for the larger simulation patches 53,70 improve the agreement of the results with experimental values already by approximately one order of magnitude.These considerations indicate that the pore opening and rupture properties of lipid bilayers are highly sensitive to the physicochemical conditions and system setup.In order to elucidate the dependence on these different factors and to evaluate the agreement between the stability of simulated and experimental membranes, further studies using larger system sizes and different physiochemical conditions will be required.Such an evaluation will become much more feasible using the protocol to estimate the free energy of nanopore formation introduced here.

Conclusions
We have demonstrated a simulation cost efficient way to determine free energy of nanopores in simulations as a function of membrane tension from a single PMF and the membrane's area versus tension curves.The accuracy of the method has been tested by comparison to a second PMF calculation at a high lateral tension.Unrestrained simulations of pore closing show that the same intermediate membrane conformations are found along the closing pathway as in the restrained PMF simulations of pore opening.This demonstrates that the chosen reaction coordinate leads along a reasonable pathway for pore formation.The high energy intermediate state is a hydrophilic pore spanning part of the bilayer, combined with a segment of a hydrophobic pore spanning the remaining distance.This state does not correspond directly to a specific value of the reaction coordinate, so that a free energy barrier does not show up in the PMF.
Using simple Markov equations and the free energies and rates measured in the simulations, rupture tensions for the simulated membranes have been estimated for high tension loading rates.A comparison of the stability of simulated membranes to experimental results show that, contrary to expectation, pores form too easily in small simulated membrane patches as used here.This may be partially attributed to the very small membrane size, as well as to the presence of ions in the experiments.For an exact evaluation of the stability of simulated membranes a similar analysis for different membrane sizes and under correct physico-chemical conditions will be required.Given the great sensitivity to these conditions, force-fields should not be parametrized to match membrane rupture stability.Membrane stability under accurate physio-chemical conditions on the other hand may represent a good quality check.

Fig. 1
Fig.1Scheme for estimating the free energy of pore formation DG total (S) at a membrane tension S, with the membrane shown in orange (headgroups) and yellow (tails) and the pore indicated in blue (waters).The process is divided into the following three steps: (a) relaxing the membrane tension s from s = S to s = 0, (b) creating a pore in the relaxed membrane and (c) stretching the membrane containing a pore from s = 0 back to s = S.

Fig. 2 A
Fig. 2 A restrained phosphate group of a DPPC lipid in the hydrophobic center of a bilayer.(a) and (b) Simulation snapshots of two alternative membrane conformations at z = 0.The system forms either (a) a hydrophilic pore or (b) a half-pore spanning only one of the monolayers.Lipid head groups are represented by orange spheres, lipid tails by green sticks and water molecules by blue spheres.The restrained phosphate group is shown in red and indicated by arrows.(c) Position z of the restrained group relative to the bilayer center and thus from the center of the umbrella potential at z = 0, as a function of time.The transition from the nanopore to a half-pore after 44 ns is reflected in a jump -indicated by the arrowsof the average position from z = 0 AE 0.001 nm to z = 0.016 AE 0.0005 nm.(d) and (e) The potential of mean force (PMF) on the phosphate group as a function of its distance z from the center of mass of the membrane (green) and the standard deviation between four different trajectory segments (red), (d) including the entire trajectories, (e) including only half-pore trajectories.

Fig. 3
Fig. 3 Snapshots of pore closing in unrestrained simulations.(a) Half-pore intermediate state, similar to the half-pore state in the restrained simulations, (b) alternative intermediate state with smaller hydrophilic indentations in both monolayers, (c) and (d) short lived hydrophobic pore segments observed in the trajectory snapshots immediately preceding the configurations in (a) and (b), i.e. 20 ps beforehand.

Fig. 4
Fig.4(a) The projected membrane area A(s) and A p (s) of the intact and porated membrane, respectively, as a function of the tension s, (b) pore radius r(s) as a function of the membrane tension.

Fig. 5
Fig. 5 Difference between the free energy of a pore in the membrane, DG total (eqn 8) and DG m total (eqn 9), and that of a pore in the tension-free membrane, DG*, as a function of the membrane tension.(a) The free energy change of the simulation box is shown in green, the membrane excess free energy change in red.(b) DG m total (S) À DG* shown together with the expression obtained using a 2nd order fit for r(s) in DG m total (S) À DG*.

Fig. 6
Fig. 6 Rupture tensions predicted from the simulations as a function of the loading rate R for (a) a membrane patch with the size of the simulation box and (b) for a 20 mm vesicle.