Kinetics of hole nucleation in biomembrane rupture

The core component of a biological membrane is a fluid–lipid bilayer held together by interfacial–hydrophobic and van der Waals interactions, which are balanced for the most part by acyl chain entropy confinement. If biomembranes are subjected to persistent tensions, an unstable (nanoscale) hole will emerge at some time to cause rupture. Because of the large energy required to create a hole, thermal activation appears to be requisite for initiating a hole and the activation energy is expected to depend significantly on mechanical tension. Although models exist for the kinetic process of hole nucleation in tense membranes, studies of membrane survival have failed to cover the ranges of tension and lifetime needed to critically examine nucleation theory. Hence, rupturing giant (∼20 μm) membrane vesicles ultra-slowly to ultra-quickly with slow to fast ramps of tension, we demonstrate a method to directly quantify kinetic rates at which unstable holes form in fluid membranes, at the same time providing a range of kinetic rates from <0.01 to >100 s−1. Measuring lifetimes of many hundreds of vesicles, each tensed by precision control of micropipette suction, we have determined the rates of failure for vesicles made from several synthetic phospholipids plus 1:1 mixtures of phospho- and sphingo-lipids with cholesterol, all of which represent prominent constituents of eukaryotic cell membranes. Plotted on a logarithmic scale, the failure rates for vesicles are found to rise dramatically with an increase in tension. Converting the experimental profiles of kinetic rates into changes of activation energy versus tension, we show that the results closely match expressions for thermal activation derived from a combination of meso-scale theory and molecular-scale simulations of hole formation. Moreover, we demonstrate a generic approach to transform analytical fits of activation energies obtained from rupture experiments into energy landscapes characterizing the process of hole nucleation along the reaction coordinate defined by hole size.


Introduction
Nearly a century ago, physiologists found that cells rupture quickly (lyse) when subjected to large osmotic pressure differences and yet remain sealed much longer before rupturing when subjected to smaller osmotic stresses (cf the classic book on red blood cell hemolysis by Ponder [1]). Similarly, 50 years ago, Deryagin and Gutop [2] observed that supported surfactant membranes (films) survive for variable periods of time and that the times diminish with increases in the film tension. Building on earlier concepts in phase nucleation [3], Deryagin and Gutop proposed that membrane rupture originates from thermal activation of an unstable hole (pore) and introduced what can now be called the paradigmatic theory for hole nucleation in tense membranes [2]. Shortly after this seminal proposition, experimenters shifted their focus to measuring permeability and charge conductance through supported membrane films, exploring the impact of voltage up to the limit of capacitive breakdown (see examples in [4][5][6][7][8][9]). Particularly relevant to thermal activation, fluctuations in conductance across membranes under low voltages reveal frequent transient-permeating defects that quickly vanish without becoming unstable and forming open nanopores [9]. Yet, in the same tests at high voltage, catastrophic holes open in the membrane that cause capacitive breakdown. In most studies of membrane conductance and capacitive breakdown, the events have been too fast and too small to capture by visible light microscopy. However, innovative approaches have been devised to briefly stabilize macroscopic holes in weakly tensed membranes, enabling quantitative studies of the fluid 3 dynamics of holes and providing an assay of the 'edge energies' ε (i.e. line tensions) for large holes [10][11][12].
In contrast to transient permeation and conductance in supported membranes, studies of cell and vesicle membrane poration have for the most part been focused on measuring the stresses (lateral tension or electrocompression) needed to break the membranes, yet neglecting the dependence of membrane survival time on the history of stress. For this reason, a few years ago we developed an approach called dynamic tension spectroscopy (DTS) [13] to quantify the most frequent tensions σ * at which giant (∼20 µm) membrane vesicles rupture when subjected to ramps of tension, σ (t) = r σ t, and thus the most frequent membrane lifetimes under increasing tension. While DTS exposes the inherent tradeoff between membrane strength and survival, it does not directly assay the kinetic rates of failure, which must be derived from the ramp speed dependence of the data using parametric or model-dependent fits. Postulating kinetic rates to be continuous functions of tension, the criterion for most frequent rupture ( * ) is specified by the crossover where reduction in the 'time' needed to open an unstable hole (1/k hole ) under rising tension surpasses the reciprocal ramp speed [13], −[∂(1/k hole )/∂σ ] * = 1/r σ . To circumvent this approach, yet utilizing the four orders in magnitude range of vesicle lifetimes obtained from DTS, we will demonstrate how the methodology can be modified to directly assay kinetic rates of vesicle failure at several times t k over the course of a repeated tension ramp experiment. Invoking Arrhenius phenomenology for thermally activated kinetics, ∼ exp(−E b /k B T ), we convert the measurements of kinetic rates k hole (t k ) on a logarithmic scale into relative changes in activation energy E b (scaled by thermal energy, k B T ≈ 4.08 × 10 −21 J), which become related to membrane tension through the product of ramp speed and time, i.e. σ (t k ) = r σ t k . The measured dependencies of energy barriers on tension are then correlated with expressions for thermal activation derived from meso-scale theory and molecular-scale simulations of hole formation. Using these correlations as examples, we end by demonstrating how analytical fits to the measured dependences of barrier energy on tension can be transformed into quantitative images of the energy landscapes governing thermal activation along a scalar reaction coordinate defined by hole size.

Stochastic process of unstable hole formation in tense membranes
When external stresses creating membrane tensions persist, as in pipette aspiration of vesicles (figure 1) or in mechanical stretch of cells, a mesoscopic hole will eventually arise that has no possibility of closing. Hence, the random process describing the formation of unstable holes can be formulated as a two-state kinetic transition between an intact and a ruptured membrane. Assuming that molecular-scale dynamics in nucleation are much faster than frequencies for appearance of unstable holes, k hole (t), we use a single Markov equation to characterize the diminishing probability S(t) of membrane survival, i.e. dS(t)/dt = −k hole (t)S(t). (1) The number of failure events per increment of time defines the probability density of rupture, p(t) = −dS(t)/dt. Hence, by establishing experimental estimators for the probability density and probability of survival and then taking their ratios at common sample times t k , we get a direct method to assay the kinetic rates for forming unstable holes in membranes, i.e.
k hole (t k ) = p(t k )/S(t k ). Of particular significance, this method of assay yields several values of the kinetic rate throughout the course of any deterministic history of tension when applied to a large population of vesicles. Dependence of the kinetic rate on tension is then established independently using the tension history, i.e. σ (t k ) = 0→tk r σ (t )dt , set by the programmed changes in tension r σ (t).

Experimental estimators for statistics of biomembrane survival
For the tests reported here (described in appendix A), the experimental approach has been to pressurize 50-100 individual membrane vesicles using micropipette aspiration to create a steady ramp of tension ( σ/ t = constant r σ ). This test was repeated with 4-6 ramps for each membrane composition selected in a range from 0.01 to 100 mN m −1 s −1 . During the course of pressurization, each vesicle was monitored by ultrafast video image processing to accurately determine its increase in surface area and its lifetime up to the moment of rupture. As the estimator for probability S(t) of survival, the vesicle lifetimes t i in a particular test were arranged into an indexed array N (t i ) of ascending order in time beginning with N total at the shortest time and ending with N = 1 at the longest time (see examples in figures 2(a) and (b)). At the same time, the statistics of failure events were organized into a discrete histogram of values N k per time interval t, thereby providing the estimator for probability density of rupture p(t k ) and the (negative) coarse-grained derivative of the survival probability at the sampling times t k (see figures 2(c) and (d). The 'raw' data of lifetimes N (t i ) from rupture experiments characterize the survival of all vesicles tested with a ramp. However, especially when made from neutral-uncharged lipids, a small percentage (<10%) of the vesicles were likely to be double-(or multi-) walled structures. Although nearly impossible to identify optically when filled with a higher-refractive-index medium, we found many years ago that multi-walled vesicles are more resistant to stretch than single-walled vesicles. Yet, elastic moduli of multi-walled vesicles do not partition into precisely quantized groups since the outermost lamellae support more stress-and displace further into a suction pipette-than the inner lamellae. Providing better discrimination, multi-walled vesicles live much longer than single-walled vesicles before rupture, creating a marked deviation in the lifetime data at its terminus (see the black crosses superposed on open black circles in figures 2(a) and (b)). Thus, when clearly evident, we truncated these short deviant tails to better represent the actual population of single bilayer vesicles (N uni ) experiencing a common welldefined tension history (see blue open circles in figures 2(a) and (b)). The few events associated with the deviant branches in figures 2(a) and (b) are also shown, highlighted by white patterned bins in the corresponding histograms of figures 2(c) and (d).

Experimental assay for kinetic rates of forming unstable holes in biomembranes
After constructing the histogram estimator for probability density of failure, p(t k ) ≈ [ N k / t]/N uni , we then interpolated amongst the array of vesicle lifetimes to acquire estimators, N (t k )/N uni , for survival probabilities at the histogram sampling times t k . As developed above, the ratios of probability densities to probabilities provided direct assays of the kinetic rates, i.e. k hole (t k ) ≈ [ N k /N (t k )]/ t, at the sampling times. Formulated in this way, the assay is independent of sample size and remains valid when transition rates vary over time as when a tension history involves increases, decreases or rapid steps to constant levels of tension. In the paragraphs that follow, we demonstrate the assay with results from tests of giant vesicles composed of single-component phospholipids and 1 : 1 mixtures of phospho-or sphingo-lipids with cholesterol (Ch). We begin by describing the results for vesicle compositions where the kinetic rates depend only on levels of tension. Then, we present results for vesicle compositions where the kinetic rates depend not only on levels of tension but also on speed of tension application. Although yet to be explained, we will show that the unexpected behaviour can be attributed to differences in the frequencies for nucleating initial symmetry breaks (i.e. 'zerosize' pores).

Kinetic rates of unstable hole formation in strong biomembranes
Recognized to be important constituents of eukaryotic cell membranes, we begin with the results obtained from experiments on mixtures of phospho-and-sphingo-lipids plus Ch, which are known to create very strong membranes [14][15][16][17]. The experimental data will show that these bilayers survive for long periods (>hours) under low tensions (<10 mN m −1 ) and yet withstand large tensions (< 20-30 mN m −1 or more) for short periods (0.1 s or less). As the first demonstration of the assay for kinetic rates (k hole ), we will focus on the tests of vesicles made from 1 : 1 mixtures of C18:0/1 phosphocholine (SOPC) plus cholesterol (1 : 1 SOPC:Ch). Compiled from rupture of ∼60 vesicles using a 10 mN m −1 s −1 tension ramp, a sample of the lifetime data is shown above in figure 2(a) along with its corresponding histogram in figure 2(c). In order to best summarize the full span of data obtained from tests with the various membrane a b c compositions, we will present cumulative log plots of the unnormalized probabilities for bilayer survival versus log time as shown for 1 : 1 SOPC:Ch vesicles in figure 3(a). Likewise, instead of plotting histograms of the data from each ramp test, we have superposed solid yellow circles on each lifetime array at the interpolation points corresponding to the histogram bin centres (see figure 3(a)); note that the bin content N k local to each yellow point is the difference between numbers of lifetimes starting from the previous time midpoint (t k − 1/2 t) to the following time midpoint (t k + 1/2 t). Determining the ratios of these bin contents to interpolated values for lifetimes, we obtain the kinetic rates, k hole (t k ) = [ N k /N (t k )]/ t. Establishing the levels of tension at each sampling time, σ (t k ) = r σ t k , we have then plotted the rates from all ramp tests on a common log scale versus tension as shown in figure 3(b) for 1 : 1 SOPC:Ch vesicles.
With each ramp test labelled by a particular solid-coloured symbol, we see in figure 3(b) that the kinetic rates of membrane failure rise continuously upward with each faster ramp speed. To confirm that the continuous progression in kinetic rates represents an explicit dependence on tension, independent of ramp speed, we have also plotted the kinetic frequencies obtained from tests using tension clamps between 17 and 21 mN m −1 (open stars) in figure 3(b). The tension clamps were created by brief periods of rapid vesicle pressurization, stopped abruptly to set a specific level of tension and held fixed until rupture. For clamps, the kinetic rates were derived  figure 3(c). Unfortunately tension clamps have a limited range of usefulness, on the one hand because vesicles can survive too long at even modest tensions and on the other hand because vesicles cannot be pressurized fast enough to reach large tensions before rupture. Even stronger and longer lived than 1 : 1 SOPC:Ch vesicles, the data for survival probabilities of vesicles made from 1 : 1 mixtures of brain sphingomyelin (SM) and Ch in figure 4(A) are seen to be shifted to much longer times when compared at similar ramp speeds. Again determined at the sample times in histograms and plotted on a log scale, the kinetic rates of failure in figure 4(b) exhibit a nonlinear dependence of failure rate on tension similar in form to that for 1 : 1 SOPC:Ch vesicles, yet shifted to significantly higher tension. Interestingly, although not as strong or long lived at any level of tension as the Ch-containing vesicles, the survival statistics (figure 5(a)) for single-component diC22:1 phosphocholine (PC) vesicles yield kinetic rates of failure that also seem to exhibit a similar functional dependence on tension (figure 5(b)), albeit shifted to lower tensions. Certainly relevant to the high strengths and longevities under stress, all three of these lipid compositions form thick membranes as indicated by their average peak-peak headgroup distances obtained from x-ray scattering using multibilayer preparations at room temperature, i.e. ∼4.24 nm for SOPC:Ch [17], ∼4.37 nm for diC22:1 PC [18] and ∼4.55 nm for 1 : 1 SM:Ch [17].

Kinetic rates of unstable hole formation in weak biomembranes
In marked contrast to the thick-strong membranes just described, bilayers made from diC13:0 PC are thin (∼3.41 nm [18]) and weak, failing in less than a minute at 1 mN m −1 and surviving no more than ∼10 ms at 4 mN m −1 . Figure 6(a) shows the cumulative log plot of survival statistics versus log time for the tension ramps applied to diC13:0 PC vesicles. Best demonstrating the fragility of diC13:0 PC bilayers, the kinetic rates of failure plotted in figure 6(b) are seen to be very fast even at low levels of tension. However, the rapid initial rise in kinetic rates follows a weakly curved form somewhat like that exhibited by strong biomembranes even though at tenfold lower tensions. Most striking, in contrast to the data for strong membranes, the kinetic rates for diC13:0 PC bilayers exhibit a marked crossover to a single exponential dependence on tension above ∼1.6 mN m −1 (highlighted by the red solid line in figure 6(b)). As discussed later, the crossover in dependence of kinetic rate on tension implies a switch in dominance between two activation barriers that exist along the thermodynamic pathway governing kinetics.

Kinetics of unstable hole formation in intermediate-strength biomembranes
Intermediate in strength and longevity between thin diC13:0 PC bilayers and thick diC22:1 PC bilayers, single-component vesicles formed from C18:0/1 PC lipid (SOPC) were found a b . Fitted by variation of pre-exponential frequency scale, the continuous blue dashed curve underlying the kinetic rates follows the dependence of activation energy on tension predicted by equation (3) in the text using the parameters appearing in table 1.
to survive for long times (∼hours) at tensions below 5 mN m −1 and less than ∼0.1 s above 15 mN m −1 . Figure 7(a) shows the cumulative log plot of survival statistics versus log time for tension ramps used in the experiments. Plotting the statistics on a log-log scale again shows decade reductions in vesicle lifetimes with tenfold increases in tension ramps as seen for weak and strong biomembranes (see figures 3(a), 4(a), 5(a) and 6(a)). Surprisingly, however, the kinetic rates of failure (figure 7(b)) obtained from the lifetime statistics are dramatically different. We see that the rates shift upward with increases in the ramp speed until the speeds exceed 10 mN m −1 s −1 . Completely puzzling, the kinetics depend on the speed of tension application-as well as tension level-suggesting that some dissipative mechanism plays a role in the process of opening a hole at the molecular scale. Moreover, we have found the same effect in all of our tests of single-component diacyl 18-carbon PC vesicles (C18:0/1, diC18:1, diC18:2 and diC18:3), although the results for each lipid differ in quantitative details (see figure 8 for diC18:2). By comparison, the data plotted in figures 5(b) and 6(b) for other single-component PC vesicles (diC22:1 and diC13:0) seem to correlate explicitly with tension-not ramp speedas do vesicles made from 1 : 1 mixtures of PC or SM with Ch (see figures 3(b) and 4(b)). Even though we have no explanation for the unusual ramp speed dependence at present, we show next that all of the kinetic rates for these PC bilayers correlate with a generic dependence on level of tension merely by rescaling each set of data at the different ramp speeds by a pre-exponential constant. a b Figure 5. (a) Data for log N (t i ) versus log t i (lifetime) encompassing the full span of tension ramps applied to vesicles formed from diC22:1 PC. Yellow solid circles identify histogram sample times at which kinetic rates of failure were determined during each ramp. (b) Kinetic rates of failure plotted on a log scale versus the tensions corresponding to sample times in ramp-test histograms (solid colored circles). Fitted by variation of pre-exponential frequency scale, the continuous dashed curve underlying the kinetic rates follows the dependence of activation energy on tension predicted by equation (3) in the text using the parameters appearing in table 2. Evidence for precursor dynamics is highlighted by data points correlating with the short red solid line (defined by r δ = 0.67 nm) vis-à-vis deviation of the continuous blue dashed curve.

Correlations of activation energy to mesoscopic models and molecular simulations
Employing Arrhenius phenomenology and treating pre-exponential factors as constants (see appendix B.3), we have converted all of the kinetic rates to increases of activation energy E b measured relative to the fastest kinetic rate k max , i.e. ln(k max /k hole ) ≡ E b /k B T . (Even though meso-scale models for rates involve pre-exponential dependences on tension (cf appendix B), the impact on correlations with the data is imperceptible in accuracy of fit; the pre-exponential scales for frequency are simply altered by a numerical factor in the range of 1-3.) Treated in this way, the kinetic measurements presented in figures 3-8 provide the experimental profiles for changes in activation energy appearing in figures 9. Illustrated by the continuous curves underlying the measures of activation energy in figure 9 and the kinetic rates in figures 3-8, we will show that the experimental data agree well with the expressions for thermal activation derived from a combination of meso-scale theory and molecular-scale simulations of hole formation. Although highlighted in the following paragraphs, formal details of the developments are to be found in appendix B.

Mesoscopic models and molecular-scale simulations
While transient permeation and opening of holes in fluid-lipid membranes have been studied with many careful experiments [4][5][6][7][8][9][10][11][12], the collective model usually employed to explain Table 2. Parameters for energy landscapes that best match the data for kinetic rates.  . Fitted separately to data for kinetic rates at each ramp speed by variation of the pre-exponential frequency scale, the continuous dashed curves underlying the kinetic rates follow the dependence of activation energy on tension predicted by equation (4) in the text using the parameters appearing in table 2. For the sake of comparison to our previous work [13], we have plotted the most frequent kinetic rates (open stars) predicted for each ramp speed by the relation the results is based on the mesoscopic theory first described by Deryagin and Gutop [2] and then reintroduced in more contemporary forms by Helfrich [19] and Litster [20]. In the D-G model, membrane holes are treated as circular defects (radius r) in a two-dimensional (2D) incompressible fluid (see the bottom sketch in figure 10(a)), contributing an excess energy proportional to the hole perimeter (i.e. 2πr ε). The material property characterizing hole energy is the edge energy/length or line tension, ε. When subjected to lateral tension, the membrane free energy is lowered by mechanical work of dilation as the hole opens [i.e. E(r ) = 2πrε − πr 2 σ ]. Well known from idealized models of cavity formation in 3D fluids, excess energy of the cavity dominates at small radii and yet is overwhelmed by work of mechanical tension at large radii. In the cavitation model, the transition to an unstable hole is marked by passing the barrier, E c = πε 2 /σ , positioned at a critical radius defined by r c = ε/σ . Thus, mediated by thermal activation, the likelihood of forming an unstable hole is characterized by an activation barrier, E b ∼ σ β /σ , that increases dramatically under tension on a scale set by edge energy and temperature [13], σ β = (πε 2 /k B T ). However, completely absent from this model are the dynamics nucleating the initial symmetry break that initiates a 'zero-size' pore. Extending the mesoscopic paradigm, recent studies of fluctuations in membrane conductance [9] indicate that collective defects of some type exist prior to opening a hole (see the top panel in figure 10(a)). Labelled 'prepores' [9], these ephemeral precursors most often vanish and, based on dynamic studies of vesicle rupture [13], seem to only rarely initiate open pores (see the bottom panel in figure 10(a)). Supporting these experimental perceptions, recent coarse-grained molecular dynamics simulations of hole formation [21,22] have revealed energy landscapes that commence from a deep energy well prior to symmetry breaks opening nascent pores (see panels in figure 10(b) taken from [21] with the permission of the authors). Once a pore opens (lower left panels in figure 10(b)), the free energies computed in the simulations were found to increase linearly with the mean-hole radius as idealized in the cavitation theory (see top right panel in figure 10(b)). Although values we present for edge energies and the length scales describing precursor structures will differ significantly from the simulation results, the models used to correlate the experimental data for activation energies in figure 9 will mirror the shapes of energy landscapes obtained in the coarse-grained molecular dynamics and illustrated by the right panels in figure 10 In contrast, collective fluctuations that thicken the membrane draw material centripetally inward as if created by a 'negative' radius hole without changing thickness.) (b) Data and images replotted from [21] with permission. Left panels: cross sections of amphiphilic bilayers computed in coarse-grained molecular dynamics simulations. For clarity, only the particles (size r 0 ) within a distance of 1r 0 normal to the planes of each cross section appear in the plots. Head-group particles in grey and tail particles in black together depict the final configurations obtained for the particular values of reaction coordinate defined by ξ = 0, 0.37, Figure 10. (Contd.) 0.80 and 0.92 (from top to bottom). Based on a particle insertion criterion, no pore was detected in the top two cross sections, whereas pores of radii R = 0.95r 0 and 1.48r 0 , respectively, were found for the bottom two panels. The circles represent the molecular-size scale; solvent particles were omitted for clarity. Top right: free energies computed in simulations of a tensionless bilayer plotted as a function of the pore radii R. The solid line denotes states with a detectable pore (ξ > 0.5). The dotted portion for R(ξ ) in the range ξ = 0 → 0.5 was derived using an extrapolation based on the requirement that R effectively measures the excess particle density. The dotted line represents a fit to the hole perimeter energy as idealized in the mesoscopic theory for cavitation. Bottom right: differences F(R) between the free-energy functions computed in the simulations and the linear fits to hole energies for zero tension (solid and dashed) and for a tension producing 3.5% membrane stretch (dotted and dashdotted). (Note, in the simulations [20], the size scale r 0 = 1/3 nm, the energy scale ∈≈ 3/4k B T (3.3 × 10 −21 J at 325 • K). Thus, the linear regime of free energy (R/r 0 > 1) divided by 2π R yields an edge energy, ε = 3.7 ∈ /r 0 ≈ 37 pN).
of the membrane at the site of pore nucleation, which is then followed by opening of the pore lumen.)

Correlations of models for activated kinetics to changes in barrier energies with tension
Applying the Kramers Brownian-kinetic theory [23,24] to minimal representations of the energy landscapes just described, we have developed analytical expressions for the dependence of activation energy on tension and thereby kinetic rates when given a pre-exponential scale for frequency. As shown in appendix B, full treatment of even a minimal landscape like that found in the simulations ( figure 10(b)) can involve many parameters, six or more in this case, creating a challenge to correlations with experiments. Fortunately, several of the parameters appear as preexponential factors scaling a paramount Arrhenius-like exponential dependence on activation energy. Hence, over the range of experimental data, the analytical expressions for kinetic rates are found to be well approximated by the product of a pre-exponential constant (frequency scale) and the exponential weight describing activation, i.e. ∼ k 0 exp(−E b /k B T ). Most important is that the pre-exponential frequency reflects the large energy barrier governing precursor dynamics and symmetry breaks (see illustrations in figure 11). Since creating a symmetry break, as well as opening a pore, causes radial expansion of the membrane, application of membrane tension acts to lower both (precursor and mesoscopic) barriers and yet impacts the exponential weight exp(−E b /k B T ) in very different ways. Correlating analytical models for kinetics to our experimental results, we will emphasize how these two barriers can be distinguished and how their properties depend on membrane composition. three parameters: (i) a frequency scale k 0 δ for escape from the deep well to initiate symmetry breaks; (ii) a thermal tension scale, σ δ = πr 2 δ /k B T , for reduction of the precursor barrier set by mean radius r δ of symmetry breaks; and (iii) a thermal tension scale, σ β = π ε 2 /k B T , for the reduction of the cavitation barrier set by pore edge energy ε, i.e.
In this model, tensions greater than σ ⊗ = ε/r δ cause the cavitation barrier E c to 'sunset' below the energy cusp E δ that characterizes the precursor barrier, thereby defining a crossover between kinetics limited by the cavitation barrier and kinetics limited by nucleation of a symmetry break. The signature for this crossover is a switch in form of the barrier energy from a curved regime ∼σ −1 at low tensions to a linear regime ∼σ at high tensions. An important experimental caveat is that the crossover tension σ ⊗ may be inaccessible due to a very small defect size r δ and large edge energies ε as discovered when correlating equation (3) with the activation energies plotted in figures 9(a) and (b) for 1:1 SM:Ch and 1:1 SOPC:Ch vesicles.
Failing to detect a crossover to linear dependence on tension for these very strong membrane materials, we can only establish upper bounds for the sizes of symmetry breaks as listed in table 1. Completing the analysis, activation energy correlations are then matched to the kinetic rates (as shown in figures 3(b) and 4(b)) to obtain the pre-exponential scales for frequency.
Again, lacking evidence for symmetry breaks, these frequency scales represent attempt rates for forming unstable holes starting from the deep-ground state and are thus very fast (see k 0 * in table 1).
By comparison, the highest rupture tensions attained for strong diC22:1 PC vesicles appear to be just sufficient to reveal a crossover in the activation energies as demonstrated by a short deviation of the cavitation regime from the few red data points for high tensions in figure 9(c). Even though short, the linear regime in equation (3) (highlighted in red) at high tensions provides a good estimate of the mean size for symmetry breaks (r δ in table 2). Matched to the curved regime ∼σ −1 over most of the tension range, the barrier energies yield an accurate measure of the hole-edge energy (also in table 2). At the same time, even though much weaker materials, the activation energies obtained from kinetic rates for diC13:0 PC vesicles are seen in figure 9(d) to also match the tension dependence given by equation (3), revealing a pronounced crossover to an extended linear regime. Here spanning more than half of the accessible experimental range of tensions, the linear regime yields a much larger size for symmetry breaks (see table 2), which seem to dominate the kinetics of failure except at very low tensions. Final fits of the activation energy correlations to both sets of kinetic rates (shown in figures 5(b) and 6(b)) have provided the pre-exponential frequency scales characterizing frequencies of symmetry breaks defined by k 0 δ (table 2). Turning to the intriguing experiments in which kinetic rates were found to depend on ramp speed as well as tension, we have found that the activation energies obtained from logarithms of kinetic rates for rupture of diC18:2 PC and C18:0/1 PC vesicles at the various ramp speeds can be collapsed to generic dependences on tension (see figures 9(e) and (f)) simply by shifting slower ramp speed data upward to match the data for the fastest ramp speeds. The shifts were accomplished by fitting each pre-exponential frequency scale (values listed in table 2). Treated in this way, the data for each type of membrane can be seen to follow the same linear regime characterizing symmetry breaks and the same cavitation regime, albeit with composition-specific parameters such as edge energy and radius scale of the symmetry break. Most simply, we find that the activation energies for diC18:2 PC bilayers correlate well with the relation given by equation (3) as evidenced by the continuous blue dashed curve and red solid line appearing in figure 9(e). While yielding a comparable low value for edge energy like diC13:0 PC, the radius scale for symmetry breaks in diC18:2 PC bilayers is found to be nearly threefold shorter (see table 2).
In contrast to the other PC bilayers, the activation energies in figure 9(f) can only be fitted by postulating a weak depression in the free energy landscape, E app * , local to the symmetry break after the precursor barrier at r δ (see figure 11(c) and the lower right panel in figure 10(b)). As developed in appendix B, the relation describing dependence of activation energy on tension (see equation (B.3.6)) involves the same three parameters in equation (3) plus the possibility of weak confinement represented by the parameter, E app * , i.e.
(4) Unlike the simpler landscape with no overshoot in energy at the precursor barrier, the energy depression, E app * , following the barrier causes the cavitation barrier to 'sunset' below the precursor barrier at a tension σ ⊗ greater than σ δ as predicted by E * = (πε 2 /σ ⊗ )(1 − σ ⊗ r δ /ε) 2 . Moreover, lacking explicit local details of the energy landscape, this condition only provides bounds on the location r * of a symmetry break, i.e. r δ < r * < r ⊗ = ε/σ ⊗ (see equation (B.3.5) for the analytical expression specifying r ⊗ ).
Correlated to the activation energies for C18:0/1 PC bilayers, equation (4) is found to closely match the superposed ramp data given a common energy depression of E * ∼ 3.6k B T . Similar to the bump in energy appearing in the simulation results, this value for E * would only produce very brief weak confinement-if at all-near the symmetry break. Still, a small energy depression is needed to achieve the quality of fit shown by the continuous blue dashed curve and solid red line underlying the data in figure 9(f). Here, the red line yields nearly the same radius scale for the location of the precursor barrier as found for diC18:2 PC bilayers. However, the energy depression E * attributed to C18:0/1 PC bilayers leads to a crossover from the cavitation regime at a lower tension σ ⊗ (i.e. large radius r ⊗ ; table 2). The correlation in figure 9(f) for tensions less than σ ⊗ yields a value for edge energy of C18:0/1 PC bilayers intermediate between the values for diC18:2 PC and diC22:1 PC. Using the correlations to activation energies defined by parameters in table 2, the apparent frequency scales characterizing symmetry breaks in C18:0/1 PC and diC18:2 PC bilayers have been derived from fits to each ramp speed data set for kinetic rates as shown by the curves underlying data in figures 7(b) and 8(b). Listed in table 2, the apparent frequency scales increase with ramp speeds for both bilayer systems and appear to reach a plateau when speeds exceed 10 mN m −1 s −1 , suggesting that fast ramps prevent relaxation to deeper ground states. While the 10-50-fold changes in frequency scale appear striking, these changes would likely only represent about 10-15% variation in depth of the ground states as sketched in figure 11(c).

Conclusions and discussion
In our previous DTS tests of diacyl PCs [13], we concluded that the most frequent tensions for membrane failure (strengths) under slow and fast tension ramps revealed a sequence of two kinetic events. When tested with slow tension ramps, the most frequent rupture tensions increased weakly with ramp speed as predicted by the mesoscopic theory for cavitation and the relation connecting kinetic rates to the most frequent rupture, −[∂(1/k hole )/∂σ ] * = 1/r σ . On the other hand, when tested with fast tension ramps, a second regime emerged at high tensions where the most frequent rupture tensions increased exponentially with ramp speed, indicating a switch of activation barrier to some type of defect nucleation process. However, fitting the crossover between the two ramp-speed regimes was often difficult and inconclusive. Hence, we set out in this work to establish a direct assay for kinetic rates of failure. Taking advantage of the 10 000-fold span in vesicle lifetimes attainable with tension ramps, we have demonstrated the efficacy of this assay with results from tests of the same membrane compositions examined previously with DTS. Now, the measurements of kinetic rates plotted on a log scale directly reveal the changes in activation energies corresponding to the levels of applied tension, enabling critical comparisons to energy landscapes for thermal activation. Following the formal approach to Brownian kinetics introduced by Kramers [23], we have developed relations for activation energies and rates of unstable hole formation utilizing minimal representations of energy landscapes that capture the essential features expected from a long history of mesosopic theory [2,19,20] and membrane experimentation [4][5][6][7][8][9][10][11][12][13] as well as recent molecular dynamics simulations [21,22]. The close fits between these expressions and our measurements provide direct support for the postulated sequence of kinetic transitions, the first leading to symmetry breaks that initiate 'zero-size' pores and the second producing unstable cavities that cause membrane failure. Moreover, the correlations yield key nanoscale parameters that govern mechanical strength and survival of biomembranes.
When first applying the experimental assay, we expected to obtain kinetic rates consistent with the values derived indirectly in the original DTS experiments. This expectation has turned out to be correct for the data obtained in tests of SOPC:Ch, SM:Ch, diC22:1 PC and diC13:0 PC vesicles. On the other hand, completely unexpected, the kinetic rates measured in tests of C18:0/1 PC (SOPC) and diC18:2 PC bilayers (as well as results not shown for other unsaturated diacyl PC bilayers) have exposed a ramp speed dependence that was not perceived in DTS experiments with the same compositions. Interestingly, even though dependent on ramp speeds, we have shown that logarithms of these kinetic rates collapse to common dependences on tension simply by shifting slower ramp speed data upward to match the data for the fastest ramp speeds. In this way, the activation energies governing cavitation after symmetry breaks, as well as the apparent size of the nucleating defects, are found to be independent of ramp speed. Yet, for some reason still to be determined, the frequency scales for precursor nucleation increase almost exponentially with slow ramp speeds as if nucleated from shallower ground states and then become independent of ramp speed above a few mN m −1 s −1 .
In correlating the data for activation energies, we have used the dependences of barrier energies on tension derived from models of free energy landscapes formulated in terms of apparent hole radius. Alternatively, we could have performed generic fits to the data for activation energies to establish the dependences on tension, E b (σ ), even choosing piecewise continuous relations when needed to match a crossover from a steep-curved regime at low tensions to a linear-like regime at high tension (see the note below). Given the thermodynamic principles used to characterize activation, the critical hole size is defined by the scalar area a b that couples with tension to lower the energy barrier as prescribed by the derivative, a b = −(∂ E b /∂σ ) T . Transformed in this way to functions of hole area, E b (a b ), the changes in barrier energy obtained with any analytical fit can be used to map the corresponding segment of the energy landscape along a scalar reaction coordinate defined by area, i.e. E(a b ) = E b + σ a b . As a simple demonstration, we have illustrated the generic approach in figure 12 using the model fits obtained with activation energies measured in tests of diC13:0 PC, C18:0/1 PC and 1:1 SM:Ch vesicles. To distinguish the thermodynamic measure of hole size from the cylindrical radius characterizing physical models and simulations, we have used an apparent diameter D hole = (4a b /π ) 1/2 to represent the hole area [= −(∂ E b /∂σ )], which is plotted as a function of tension in figure 12(a). The size scale D hole provides easy comparison to the average cross sections of phospholipids in bilayers (∼1 nm). Figure 12(b) shows the barrier energy fits transformed to functions of the apparent hole diameter. Finally, adding the mechanical potential, σ a b , the corresponding energy landscape segments appear in figure 12(c). As a simple confirmation of validity, application of the approach to data for 1 : 1 SM:Ch bilayers reproduces the expected theoretical behaviour for idealized cavitation: i.e. D hole (= 2ε/σ ) is inversely proportional to tension in figure 12(a); barrier energy, E b /k B T [= (πε/2k B T ) (D hole )] rises linearly in proportion to D hole in figure 12(b); and E [= 2 E b ] also rises linearly in proportion to D hole at twice the slope in figure 12(c). Nontrivial by comparison and unique to this treatment, the energy landscape for C18:0/1 PC bilayers reveals a nonlinear dependence on hole diameter at the onset of the cavitation regime following the symmetry break (red stars in figure 12(c)).
(Note: the series expression, C 0 − σ/σ δ + n 2 C n (σ ) n , involving powers of an orderparameter like function, (σ ) ≡ (1 − σ/σ ⊗ ) for σ σ ⊗ and (σ ) ≡ 0 for σ > σ ⊗ , provides equally as good fits to the data for activation energies as found with the model-based equations (3) and (4)  Quantified by the approach just described, our experiments show that the sizes of rupture pores in biomembranes appear to be small, remaining molecular in scale even for the lowest rupture tensions. As such, the results can be compared directly to the molecular-scale simulations of hole formation [21,22]. First, although small, the length scales derived from the precursor regimes appearing in figures 9(d)-(f) (see r δ in table 1) for single-component PC bilayers seem to be significantly larger than those revealed in coarse-grained molecular dynamics simulations of diC16:0 PC lipids [22], which in turn are consistent with length scales obtained in simulations of structureless particles [21]. However, of greater concern, much larger edge energies (30-40 pN) are obtained from simulations [21,22,25] than derived here from the fits to the data for activation energies. Given a range from 30 to 40 pN, the thermal tension scales predicted to characterize activated hole formation would be so large, σ β ∼ 700-1200 mN m −1 , that membrane rates of rupture would be expected to rise by 3-4 orders of magnitude over a mere 10% increase in applied tension. The dramatic increase in barrier energy for small change in tension would make membranes appear 'brittle' and would imply that thermal excitations involve ultrafast 'sonic' transport (as in fracture of solids). By comparison, the span of kinetic rates obtained over a large range of tension in experiments reveals the role of overdamped-Brownian excitations in failure of fluid biomembranes. Further supporting this view, it is important to note that edge energies derived from hydrodynamic responses of micron-sized holes in similar bilayers [11,12] lie in the mid-range of ε values listed in table 1. A possible reason for the significant difference between edge energies of molecular simulations and vesicle experiments could be that membrane configurations in simulations only have enough time to equilibrate through local chain extension and density relaxation. Perhaps to fully equilibrate, collective excitations on much longer time scales may be needed that move (shear) molecules relative to one another with little change in density yet create significant viscous dissipation (see [11,12]).
suction P(t) was implemented by a motorized ground-glass syringe pump connected to the micropipette assembly. The value of tension at each video time point was calculated from the programmed suction pressure P(t) using the simple relation, σ (t) = P(t)R p /2(1 − R p /R v ), for mechanical equilibrium of a fluid membrane vesicle in which the pipette radius R p and radius R v of the vesicle segment exterior to the pipette describe the two unsupported membrane geometries. Produced either by steady ramps or steps in suction, the tension rate r σ (t) describing the stress history experienced by a vesicle was calculated directly from the change of tension with time. Measured relative to the first application of stress, each vesicle lifetime was determined by the time span up to its moment of disappearance. Compiled from ∼ 50-100 vesicle ruptures, the data for lifetimes at a particular tension history were arranged into an ascending array of lifetimes beginning with N total vesicles surviving to the shortest lifetime and ending with one vesicle surviving at the longest lifetime (see examples in figures 2(a) and (c)).
applied to the 1D energy landscapes modelling hole formation provides the same outcome as an analysis of 'mean first passage times' [24].)

B.2. Rates of forming unstable holes predicted by idealized cavitation
In the cavitation-based theory, thermal fluctuations δr/δt in hole radius space are confined relative to a steric barrier at zero radius by the linear rise in hole energy, 2πε r, as set by the material 'edge energy' ε. At low-to-modest tensions, the thermal spread in states near r = 0 is little affected by the applied tension and can be approximated by l 0 = l ε ≈ k B T /(2π ε). However, the added mechanical potential, −πr 2 σ , overcomes the hole energy at large distance to produce a finite-energy barrier at a critical radius r c = ε/σ . When normalized by thermal energy, the energy barrier, E c = π ε 2 /σ , leads to a generic Arrhenius factor, exp(−σ β /σ ), that describes the mechanical susceptibility of kinetics to tension as set by the thermal scale for tension, σ β = π ε 2 /k B T . Since the curvature of the energy landscape |∂ 2 E/∂r 2 | is proportional to tension, the asymptotic (Gaussian) integral for thermal width local to the transition state is approximated to first order by l ts = l c ≈ (k B T /σ ) 1/2 when σ β σ . Hence, ignoring the nontrivial dynamics creating incipient holes and symmetry breaks local to the origin, we obtain a compact-generic expression for kinetic rate of hole formation predicted by the idealized energy landscape for cavitation [13], In this case, the attempt frequency k 0 * scaling the rate is simply given by the ratio of thermal tension to the viscous coefficient that damps the thermal fluctuations, i.e. k 0 * = (σ β /ζ ). (Note: for 2D incompressible membranes, fluctuations in hole radius r drive radially divergent collective flows impeded by a surface-shear viscosity η and locally characterized by the energy dissipation per area, ∂e/∂t = 4πη[∂ln(r )/∂t] 2 . When integrated over the entire membrane, it follows [13] that the coefficient of damping fluctuations in hole radius becomes ζ = 4πη.)

B.3. Kinetic rates of forming unstable holes seeded by precursor nucleation
Requiring many more parameters, we now apply Kramers theory to an energy landscape for hole opening that involves two regimes of nucleation, perhaps even punctuated by a weak metastable intermediate. As evidenced in earlier experiments [9,13] and supported by recent molecularscale simulations [21,22], collective precursor excitations seem to precede symmetry breaks that initiate an incipient pore, thereby seeding the process of unstable cavity formation. Even when closed, these excitations act to locally thin/thicken the membrane introducing small area expansions/condensations that can couple with tension to lower the activation barrier governing hole opening. Although seeming to be dominated by the energetics of cavity formation, nucleation of the symmetry break can become the rate-limiting process at high levels of tension.
Snapshots from molecular-scale simulations (left panels in figure 10 taken from [21]) suggest that collective 'pinches' and 'bulges' of the membrane structure seem to occur local to the symmetry break that creates an incipient pore. Concomitant to these precursor excitations, the simulation reveals a free-energy landscape that commences from a deep well rising by ∼ 15-20 k B T to terminate in a 'cusp' at a very small radius r δ . Beyond the cusp (r > r δ ), the free energy appears to increase in proportion to the pore perimeter [∼2π r] with no perceptible effect of the applied tension, yielding an apparent edge energy ε ∼9 k B T nm −1 (upper right panel in figure 10 taken from [21]). However, subtracting out the linear increase (i.e. as if subjected to very large tension σ ≈ ε/r δ ), the authors exposed a small excess (∼ 2-3 k B T ) in the free energy just beyond the cusp (lower right panel in figure 10 taken from [21]). The actual break in symmetry creating a pore seems to follow the small energy bump at a radius r * , which is much larger than r δ . Even though small and imperceptible at low tension, the excess in energy suggests a possibility that a weak metastable intermediate could exist at the symmetry break. Interestingly, in related coarse-grained simulations using detailed chemical models of phospholipid molecules [22], the results indicate a similar excess in free energy relative to a linear open-pore regime characterized by nearly the same edge energy. Motivated by these features, we will apply Kramers theory to a modified energy landscape that commences with a precursor regime nucleating symmetry breaks and initiating 'zero-size' pores.

B.3.1. Modelling an energy landscape punctuated by weak metastable confinement.
Lacking specific details of the free-energy contour in radius space, we can only crudely characterize the energy barriers and minima that describe the sequence of two kinetic regimes. Guided by the discussion above, we begin by defining a 'cusp' energy E 0 δ at the terminus r δ of a deep-initial well, which represents the precursor dynamics seeding symmetry breaks. Next, we assume that the symmetry break arises at a position r * beyond the cusp. The energy level at the symmetry break is defined relative to the continuation of an apparent hole-perimeter energy profile mapped outward from the cusp, E 0 δ + 2π(r * − r δ )ε, yet reduced locally by a small energy, − E * , possibly leading to weak confinement, i.e. E 0 * = E 0 δ + 2π(r * − r δ )ε − E * . Although shifted downward by E * at r * , the linear hole-perimeter energy then continues to rise indefinitely upward as E 0 * + 2π(r − r * )ε, characterizing larger hole radii, r > r * . Taking mechanical tension into consideration, we lower the energy everywhere along the contour by the area-dependent potential, −πr 2 σ . Thus, eventually the linear rise in hole-perimeter energy will be overcome by the mechanical potential to establish a transition state at the radius r c beyond which holes are unstable. When described relative to the symmetry break, E * = E 0 * − πr 2 * σ , the height of the energy barrier for cavitation is given by At the same time, the mechanical potential will also lower the cusp energy in proportion to a local area (neglecting higher-order contributions), i.e. E δ ≈ E 0 δ − πr 2 δ σ . Hence, the energy difference, E δ − E * ≈ E * − 2π(r * − r δ )ε + πσ (r 2 * − r 2 δ ), characterizes the possibility of weak confinement. Notably, E δ − E * has to be greater than zero for confinement; hence, at zero tension, E * must exceed 2π(r * − r δ )ε. At the same time, this simple model predicts that increasing tension enhances weak confinement at the symmetry break.
Turning now to the Kramers stationary-flux approximation, we represent the ground state and intermediate state densities [ρ(0),ρ(r * )] as probabilities for being in the vicinity of each position divided by length scales (l 0 , l * ) that characterize local confinement, i.e. ρ(0) = P 0 /l 0 and ρ(r * ) = P * /l * , recognizing, P 0 + P * = 1. Hence, given a constant flux over the 1D landscape, we establish the two probabilities by integrating the energy weighted flux {k hole exp[E(r )/k B T ]} first from the ground state to the symmetry break, i.e.
and then integrating from the symmetry break to the adsorbing-end state, Implicit in each step, integration over an energy weight becomes the product of a Boltzmann weight for the transition-state energy multiplied by a thermal length scale, i.e. one (l δ ) local to the cusp and the other local to the hole transition (l c ). Along with the confinement lengths (l 0 , l * ), these four length scales would in principle be evaluated using detailed models of the energy landscape in the Boltzmann-weighted integrals. However, since these parameters appear in pre-exponential coefficients, we will ignore details of their computation and show where they become embedded in phenomenological parameters scaling frequency.
Rearranging the integration results, we obtain the kinetic rate for opening unstable holes in terms of activation over two barriers given the probability of being in the ground state, and determine the ratio of probabilities, In the case of large barrier energies and weak metastable confinement, the ratio of probabilities is extremely small, P * /P 0 1. Thus, we can use equation (B.3.2) with P 0 ≈ 1 to approximate the kinetic rate. Although the energy barriers impeding each transition are the primary factors governing kinetics in (B.3.2), we will see that weak metastable confinement can also have a significant effect.
Incorporating all exponential dependences on tension in equation (B.3.2), we acquire an expression that predicts the kinetic rate for forming unstable holes nucleated by precursor excitations and the subsequent transition to a symmetry break, where there are now two thermal tension scales, σ β ≡ πε 2 /k B T , σ δ ≡ k B T /(πr 2 δ ), and a frequency scale, k 0 δ ≡ [k B T /(ζ l 0 l δ )] exp(−E 0 δ /k B T ), for escape from the deep well, τ δ = 1/k 0 δ . We see that the effect of an energy depression at the symmetry break, E * , is to cause the cavitation barrier to 'sunset' below the cusp energy [E c < E δ ] when tension reaches a value σ ⊗ determined by E * = (π ε 2 /σ ⊗ )(1 − σ ⊗ r δ /ε) 2 , thereby establishing the crossover between kinetics limited by passing the cavitation barrier and kinetics limited by nucleation of a symmetry break. However, at this level of approximation, there is no explicit dependence on location r * of the symmetry break other than that it lies between r δ < r * < r ⊗ = ε/σ ⊗ , i.e. r δ < r * < r δ [1 + e * + (2e * + e 2 * ) 1/2 ], (B.3.5) as defined by the parameter e * = E * /(2πr δ ε).
(Note: the pre-exponential dependence on the ratio of lengths (l c /l δ ) in equation (B.3.4) can be incorporated as a hidden contribution to the apparent energy depression at the symmetry break, i.e. E app * /k B T = E * /k B T − ln(l c /l δ ).)

27
Viewed in the context of Arrhenius phenomenology, −ln(k hole ) ≈ E b /k B T , we will use the kinetic rate in equation (B.3.4) to define the effective barrier energy for opening unstable holes, By comparison, in the absence of an intermediate state, we obtain a simpler relation for the effective activation barrier characterized by the dynamical crossover, σ ⊗ = ε/r δ , E b /k B T + ln(k 0 δ ) ≈ −σ/σ δ + (σ β /σ )(1−σ r δ /ε) 2 . (B.3.7) As r δ → 0, and ignoring pre-exponential factors, equation (B.3.7) approaches the classical expression in equation (B.2.1) for idealized cavitation, E b /k B T + ln(k 0 δ ) ≈ σ β /σ . [21,22] Treating the evolution of the bilayer structure and metastable prepore as states in an activated reaction process, the authors of [21,22] have computed a free-energy landscape for the process using a method called the potential of mean constraint force (PMCF). The following quote nicely summarizes the coarse-grained approach in the first paper [21]: 'The amphiphiles consist of one head particle and four tail particles, the latter representing three to four CH 2 groups each. The solvent particles, representing two water molecules, are identical to the head particles. Most non-bonded interactions are modeled by the usual Lennard-Jones potential LJ (r ) = 4 ∈[(r/r 0 ) −12 − (r/r 0 ) −6 ], with r 0 = 1/3 nm and ∈= 2 kJ mol −1 ≈ 3/4 k B T . Exceptions are the hydrophobic interactions between tail particles and head or water particles, which are mimicked by a purely repulsive soft core potential sc (r ) =∈ (r/r sc ) −9 , with r sc = 1.05r 0 .... The initial simulation boxes were created by constructing two parallel square lattice layers holding 24×24 straight amphiphilic molecules each, with their heads pointing outward and by randomly adding 10 800 solvent particles. The bilayers were oriented parallel to the square xy face of the box. In the first box the sides of this face were matched to produce a tensionless state, which for the current model amounts to an area of A 0 = 1218r 2 0 (Note: to avoid confusion with the variable σ we use for tension, we have replaced the same variable used by the authors to define molecular length scale by r 0 and the scale for figure abscissas.)

B.4. Comments on the coarse-grained molecular dynamics simulations in
When equilibrating the free energy at positions constrained along a specialized reaction coordinate ξ , the authors used a particle insertion routine to determine when a pore opened and to provide a measure of its radius R. Evaluated in this way, open holes were established for values, ξ > 0.5. However, the particle insertion routine yielded spurious-negative radii at reaction coordinates of ξ < 0.5. To circumvent the anomaly, the authors used linear extrapolation to determine the radius as ξ → 0 from the last detectable positive radius. Interestingly, the extrapolated regime (R < 0.25r 0 ) encompassed a deep minimum in the free energy. In figure 10(b) (left) taken from [21], we see examples of bilayer cross sections obtained along the reaction coordinate in the vicinity of the symmetry break opening a hole. Adjacent in figure 10(b) ( top right), we see that free energies rise rapidly from the minimum up to ∼ 15 k B T at R = 0.25r 0 . As noted by the authors, the tail particles were gradually expelled from the centre of the bilayer, reducing the thickness of the bilayer at the origin and increasing values of R even before the hole opens. When pushing the tail particles away from the centre, the molecular directors tilt in order to minimize contact between the remaining tail particles and the solvent. Then, at R > 0.75r 0 , the symmetry break occurs and the diametrically opposed head-group