A joint reaction coordinate for computing the free energy landscape of pore nucleation and pore expansion in lipid membranes

Topological transitions of membranes, such as pore formation or membrane fusion, play key roles in biology, biotechnology, and in medical applications. Calculating the related free energy landscapes has been complicated by the fact that such processes involve a sequence of transitions along highly distinct directions in conformational space, making it difficult to define good reaction coordinates (RCs) for the overall process. In this study, we present a new RC capable of driving both pore nucleation and pore expansion in lipid membranes. The potential of mean force (PMF) along the RC computed with molecular dynamics (MD) simulations provides a comprehensive view on the free-energy landscape of pore formation, including a barrier for pore nucleation, the size, free energy, and metastability of the open pore, and the energetic cost for further pore expansion against the line tension of the pore rim. We illustrate the RC by quantifying the effects (i) of simulation system size and (ii) of the addition of dimethyl sulfoxide (DMSO) on the free energy landscape of pore formation. PMF calculations along the RC provide mechanistic and energetic understanding of pore formation, hence they will be useful to rationalize the effects of membrane-active peptides, electric fields, and membrane composition on transmembrane pores. Graphical TOC Entry


Introduction
The opening and closing of pores over lipid membranes play critical roles in various processes of biophysical, biotechnological, or pharmaceutical relevance. Many antimicrobial peptides function via the formation of transmembrane pores, thereby causing the loss of ionic concentration gradients over bacterial membranes. 1 Pore-forming toxins are produced as virulence factors by many pathogenic bacteria. 2 Pore formation is a critical step during membrane fusion, as required for the infection of cells by enveloped viruses, for exocytosis, and for other forms of active transport. 3,4 In turn, pores are required to close during endocytosis, that is, by transport of cargo into the cell via the formation of an intracellular vesicle. Pores are also formed on purpose with applications in medicine, for instance with the aim of drug delivery, gene therapy, or killing of malignant cells. 5 For biotechnological applications, the opening of membranes pores is often induced by applying electric fields, a process termed electroporation. 6,7 In biophysical experiments, pores were also induced by the application of osmotic stress or membrane tension. [8][9][10] A quantitative understanding of such processes requires understanding of the free-energy landscape of pore formation.
Computing the free energy landscape of pore formation, for instance from molecular dynamics (MD) simulations with a biasing method such as umbrella sampling or constrained simulations, is complicated by the fact that pore formation involves two highly distinct conformational transition, here referred to as (a) pore nucleation and (b) pore expansion. 11 Because pore nucleation and pore expansion occur along different directions in conformational space, it has been difficult to define a unified reaction coordinate (RC) for the overall process. Namely, during nucleation, the membrane is locally thinning, followed by the penetration of water and head groups into the hydrophobic core of the membrane. [12][13][14] Hence, nucleation is primarily characterized by a transition of polar atoms perpendicular to the membrane plane. In contrast, expansion of an established pore involves the increase of the pore radius, that is a transition of the (already established) pore rim parallel to the membrane plane. Notably, the expansion of large pores was successfully modeled by classical nucleation theory (CNT). 15,16 CNT describes the pore free energy as function of the pore radius R as where σ denotes the line tension of the pore rim, and γ is the surface tension between membrane and water. However, because CNT breaks down at small pore radii, it is not suitable for describing the energetics during pore nucleation. 11 Early simulation studies of pore formation were carried out under non-equilibrium conditions, where the pore was induced by applying an external perturbation, such as an electric field or mechanical stress. [17][18][19][20][21] From such non-equilibrium simulations, it is difficult to obtain the free energy cost of pore formation. To obtain the potential of mean force (PMF) of pore formation, some studies used the pore radius R as RC. 22,23 Because R is not suitable for following the early steps of pore formation during the nucleation process, several alternative RCs have been proposed, including a RC that steers the lipids collectively in lateral direction away from the pore center, 24 a lipid flip-flop coordinate, steering a single lipid head group towards the membrane core, 13,[25][26][27] or the mean solvent density within a membrane-spanning cylinder. 28 In a systematic comparison of RCs of pore formation, we found that pulling along these established RCs for pore nucleation efficiently triggered pore formation. However, if the open pore is metastable (long-living), pulling along these RCs in reverse direction did not reliably close an established pore, primarily because the continuous hydrogen bond network along the transmembrane pore was not reliably disrupted by pulling along these RCs. 29 Such different pathways during forward and reverse pulling may lead to severe hysteresis problems during PMF calculations. The problems prompted us to develop a RC that probes the degree of connectivity of a polar transmembrane defect, here referred to as "chain coordinate" ξ ch , which is briefly described in the Methods. We found that pulling along ξ ch allows for reliable pore opening and closing, thereby avoiding any hysteresis problems during PMF calculations. 11,30 Using the chain coordinate, we resolved for the first time a nucleation barrier and, thereby, the metastability of the open pore in certain membranes, which had been proposed forty years ago to rationalize memory effects in electrophysiology experiments. 6 An alternative biasing potential for efficient pore opening and closing, coined "gizmo", was recently proposed based on a set of repellant beads. 14 Here, a linear, membrane-spanning arrangement of beads was used to drive pore opening, and a ring-shaped arrangement of beads was used to squeeze a transmembrane polar defect to drive pore closing. In the same study, a careful committor analysis showed that the penetration of head groups into the membrane core, and not merely of water, is critical for crossing the transition state of pore nucleation.
While the chain coordinate ξ ch has been designed to follow pore nucleation and to resolve a nucleation barrier (if present), it is not suitable to follow pore expansion because all pores with larger radii are projected onto ξ ch = 1, characterizing a fully formed transmembrane defect. Consequently, PMFs along ξ ch do not clearly resolve the free energy minimum and the size of an open, metastable pore, and they do not provide the free energy cost for further pore expansion against the line tension of the pore rim. In this study, a new RC is introduced that combines nucleation along ξ ch with expansion along the pore radius. We show that a single PMF along the new RC reveals the complete free energy landscape reflecting the unperturbed membrane, pore nucleation, a metastable pore, and further pore expansion.
We illustrate the RC by analyzing the effects of membrane size and of dimethyl sulfoxide (DMSO) on the free energy landscape of pores. The chain coordinate ξ ch quantifies the connectivity of a polar defect as the fraction of polar-occupied slices of a transmembrane cylinder. Polar-occupied slices are highlighted by blue shading. Only polar heavy atoms within the cylinder contribute to ξ ch (solid blue circles), but not polar atoms outside the cylinder (empty blue circles). We considered oxygen atoms of water and phosphate groups as polar atoms contributing to the reaction coordinate. (B) After a continuous polar defect has formed (ξ ch ≈ 1), the pore expands along the radius of the pore, quantified via the number of polar heavy atoms within a horizontal slab at the membrane center (blue solid circles inside the orange shading).

Definition of the reaction coordinate (RC)
The new RC is defined as a combination of two coordinates. First, we quantify the state of pore nucleation with the chain coordinate ξ ch (r) defined in Ref. 30, in which it was simply denoted by "ξ". The symbol r denotes the coordinates of all atoms of the system. ξ ch (r) was designed to quantify the connectivity of a polar transmembrane defect. In brief, ξ ch (r) is defined via a transmembrane cylinder that is decomposed into N s slices with a thickness of typically 1Å, as illustrated in Fig. 1A. ξ ch (r) is approximately given by the fraction of slices that are filled by polar heavy atoms (Fig. 1A, blue shaded slices). As polar atoms contributing to ξ ch , we typically used oxygen atoms of water and lipid phosphate groups. Hence, both water and head groups contribute to the degree of connectivity that is quantified by ξ ch . Upon pulling the system along ξ ch , the slices are filled one-by-one, thereby enforcing the formation of a chain of hydrogen-bonded polar atoms -hence the notion "chain coordinate". The length of the cylinder is chosen such that ξ ch takes a value of ∼0.2 for the flat, unperturbed membrane, indicating that 20% of the cylinder slices located at the two head group regions are filled by polar atoms. A value ξ ch ≈ 1 indicates a fully formed pore. A thin continuous water needle, flanked by indented head groups, is typically given by ξ ch (r) ≈ 0.85, corresponding to the transition state of pore formation in membranes that exhibit a pore nucleation barrier. Because all larger pores with increased radii are projected onto ξ ch = 1, the chain coordinate is inappropriate for studying pore expansion. Instead, pulling along ξ ch may efficiently open or close transmembrane defects.
Second, we quantify the state of pore expansion via the approximate radius of the pore R(r), which we define via the number of polar heavy atoms n P (r) within a horizontal layer of thickness D at the center of the membrane (parallel to the membrane plane, see Fig. 1B).
Because the defect takes approximately the shape of a cylinder, we define R(r) as We used D = 1 nm in this work. The symbol v 0 denotes the volume per polar atom, here taken from the molecular volume of water, v 0 = 29.96 Å 3 . Hence, pore expansion is quantified via the number of polar atom n P inside the hydrophobic membrane core; we merely translate the latter quantity into an approximate pore radius to obtain a more intuitive interpretation of the RC. Similar to previous studies, 28,30 n P was computed using a differentiable indicator function that takes unity inside the horizontal layer, and zero outside, as follows: Here, θ(x; h) is a differentiable step function that is zero for |x| > 1 + h, and unity for For the definition of θ(x; h) and its derivative we refer to previous work. 30 Z mem denotes the center of mass of a reference group, here defined by all heavy atoms of the lipid tails.
We combine these two coordinates for pore nucleation (ξ ch (r)) and pore expansion (R(r)) into a single RC as follows: Here, H (x) denotes a smoothed, differentiable Heaviside step function switching between 0 and 1 in the interval [− , ], implemented using 3 rd -order polynomials as follows: Hence, H ξ ch (r) − ξ s ch switches gradually between zero and one as the chain coordinate ξ ch (r) passes a switching value ξ s ch . We used ξ s ch = 0.925, characterizing a continuous polar defect, and we used = 0.05. The parameter R 0 serves as a normalization constant with two purposes: first, it renders the second expression on the right hand side of Eq. 4 unitless; second, by setting R 0 to the value that R(r) takes where ξ ch (r) ≈ ξ s ch , the term [R(r)−R 0 ]/R 0 takes zero in the switching region between pore nucleation and pore expansion, and then grows to finite values as the pore expands to larger radii.
Taken together, the proposed RC ξ p (r) quantifies the radius of the pore in units R 0 . The range 0 < ξ p (r) ξ s ch follows the nucleation of the pore, up to the formation of a thin pore with a radius of approximately R 0 . Larger values of ξ p (r) follow the expansion a fully formed defect with a radius of approximately ξ p (r)R 0 .
Choosing the parameter R 0 The parameter R 0 does not change the overall shape of the free energy landscape, as different values of R 0 merely shift and stretch the RC. However, to make sure that ξ p (r) quantifies the pore radius in units of R 0 , we choose R 0 as the average value of R(r) in the switch region, i.e., where Z is the partition function, H the Hamiltonian, β the inverse temperature, and δ the Dirac delta function. Specifying R 0 via Eq. 6 would require a simulation of pore formation prior to the production simulations.
However, we found that a reasonable estimate for R 0 can be obtained as follows, thereby respectively. Then, we estimated the number of polar atoms in the central layer (and thereby the value of R(r)) that leads to ξ ch ≈ ξ s ch . To this end, we assume that the N s,o slices outside of the central layer are fully occupied by polar atoms before a continuous polar column is formed inside the layer, yielding ξ ch = N s,o /N s (the fraction of fully occupied slices). Next, we randomly add polar atoms to the N s,i slices inside the central layer, until ξ ch ≈ ξ s ch . This procedure is repeated 10000 times the beginning of the simulations. Finally, the average number of polar atomsn P added to central layer was used to estimate R 0 via Eq. 2.
We found that this estimate is in good agreement with the average R(r) in simulations with ξ ch ≈ ξ s ch . For instance, during umbrella sampling simulations with 128 DMPC lipids (N s = 26, D = 1 nm, d s = 0.1 nm, see below), in which the average value of the switch function was H ξ ch (r) − ξ s ch = 0.50 (indicating ξ ch (r) = ξ s ch ), we found R(r) = 0.448, close to the value 0.443 for R 0 obtained with the simple procedure presented above. As such, R 0 is not a free parameter but is specified at the beginning of each simulation based on the parameters N s , D, d s , and ξ s ch .

Implementation
The implementation of the RC ξ ch was taken from our previous study, 30 but merged with Gromacs 2018.8. 31 The new RC ξ p was implemented directly into the pull code of Gromacs 2018.8. Because the calculations of the RC and the pull forces use OpenMP parallelization, the code ran only approx. 7% slower as compared to a free simulation on a server equipped with a 6-core Intel Xeon E-2136 and an Nvidia GeForce GTX 1070Ti graphics card.

Simulation setup and parameters
The membrane systems were set up with the MemGen webserver (http://memgen.unigoettingen.de). 32 Equilibrium simulations were carried out with Gromacs 2019.6, and simulations using the new RC with an in-house modification of Gromacs 2018.8. Lipids were modelled with the force field by Berger et al., 33 and the SPC water model was used. 34 Parameters for DMSO were taken from Ref. 35. Bonds and angles of water were constrained with the SETTLE algorithm. 36 All other bonds were constrained with LINCS. 37 The temperature of the simulations was controlled at 300 K using velocity rescaling and by coupling water and lipid to separate heat baths. 38 During equilibration, the pressure was controlled at 1 bar using a semi-isotropic weak coupling scheme (τ = 1 ps). 39 Although weak coupling does not yield a well-defined ensemble, we used it here owing to its numerical stability. Electrostatic interactions were calculated using the particle-mesh Ewald method. 40,41 Dispersion interactions and short-range repulsion were described by a Lennard-Jones (LJ) potential, which was cut off at 1 nm. An integration time step of 5 fs was applied. All systems were equilibrated until the box dimensions and potential energies were fully converged.

Simulations of pore opening and umbrella sampling simulations
PMFs were computed using umbrella sampling simulations. 42 Starting frames for umbrella sampling simulations were taken from constant-velocity pulling simulations either along ξ ch or along ξ p carried out for at least 60 ns or 120 ns, respectively. The following force constants were used: k = 3000 kJ mol −1 for ξ p ≤ 0.75, k = 5000 kJ mol −1 for 0.75 ≤ ξ p ≤ 1.1, and k = 400 kJ mol −1 otherwise. For systems containing no DMSO, each window was simulated for 100 ns, and the first 40 ns were omitted for equilibration. For systems containing DMSO, each window was simulated for 50 ns, and the first 10 ns were omitted for equilibration. During umbrella sampling, the pressure was controlled with the Parrinello-Rahman barostat (τ = 1 ps). All other parameters were chosen as described above.
The PMFs were computed with the weighted histogram analysis method, 43 as implemented in the gmx wham module of Gromacs. 44 Statistical errors were estimated using 50 rounds of Bayesian bootstrapping of histograms. 44 Here, each PMF computed from bootstrapped histograms was defined to zero at the PMF minimum, corresponding to the flat membrane, before computing the standard deviation between the bootstrapped PMFs. The thickness of the cylinder slices was set to d s = 0.1 nm. The degree to which a cylinder slice was considered as filled upon the addition of the first polar atom was taken as ζ = 0.75, as described previously. 30 The number of slices N s was chosen such that the ξ ch ≈ 0.2 referred to the flat membrane. This convention yields 26 and 32 slices for DMPC and POPC membranes, respectively.
To exclude that the integration time step or the algorithm used for temperature coupling influence the PMFs, we computed a PMF for pore nucleation along ξ ch coordinate using two  Previously, we found that it is critical to define the lateral position of the cylinder dynamically, allowing the cylinder to follow the defect. This protocol excludes that the system moves along the RC by shifting the defect laterally out of the cylinder, whereupon the nucleation barrier would be integrated out (Fig. 2, black curve). With a small simulation system of 128 lipids, simulations with such a mobile cylinder are numerically stable. However, upon simulating pore formation in large membrane patches, we found that using a mobile cylinder in umbrella windows with a nearly flat membrane may lead to integration problems, merely because water fluctuations at the glycerol region cause rapid lateral movements of the cylinder. To solve such numerical problems, we here fixed the cylinder in umbrella windows with reference position ξ ch < 0.7, thereby avoiding integration problems, while leaving the cylinder mobile in windows with ξ ch ≥ 0.7, thereby properly sampling the nucleation barrier. As shown in Fig. 2 (blue and magenta curves), fixing the cylinder at small ξ ch leads to a slightly  A single PMF for nucleation and pore expansion defect with a locally thinned membrane due to indented head groups, a thin water needle, up to a fully formed large pore (Fig. 4A-H). The final snapshots of all umbrella windows are shown in Movie S1. In addition, to reveal the average shape of the defect at various Comparing the PMFs for DMPC and POPC reveals a number of differences. First, the free energy of pore nucleation, given by the height of the barrier, is increased in POPC as compared to DMPC by 22 kJ/mol. This finding is compatible with previous studies that demonstrated that pore nucleation free energies correlate with membrane thickness. 11,26 Consequently, the free energy of the open pore is lower in DMPC than in POPC by 27 kJ/mol. Second, the depth of the shallow local minimum ∆G cl at ξ ≈ 2 (relative to the barrier) is 10 kJ/mol for DMPC, but only 5 kJ/mol for POPC. Assuming that the rate of closing the pore is proportional to exp(−β∆G cl ), these values imply that the pore in DMPC has an ∼8 times longer lifetime as compared to the pore in POPC. Structurally, this difference is rationalized by the larger spontaneous curvature of DMPC: with it larger headgroup-totail volume ratio, DMPC can better accommodate the large curvature at the pore rim. In contrast, the radius of the metastable pore is only marginally larger in DMPC as compared to POPC. The radius is given by ξ p R 0 at the shallow PMF minimum, where the parameter

Switching the reaction coordinate from nucleation to expansion
To illustrate how the RC ξ p accounts for the switch from pore nucleation to pore expansion, In addition, Fig. 5 demonstrates the relation between the pore radius R and ξ p . During pore nucleation (0 < ξ p 0.9), R increases up to the unit radius R 0 , corresponding to the radius of a continuous thin polar defect, where R 0 = 0.38 in this simulation (blue symbols).
Henceforth, during pore growth, R increases linearly as ξ p R 0 . Overall, Fig. 5 illustrates how a restraint along ξ ch is automatically turned into a restraint along R, thereby driving pore nucleation and expansion by restraining only a single degree of freedom. To probe the influence of system size on the free energy landscape of pore formation, we computed the PMF along the newly proposed RC using systems of 128, 200, or 512 DMPC lipids. As shown in Fig. 6, the free energy of pore nucleation, as indicated by the local maximum at ξ p ≈ 0.9, is not influenced by system size. In contrast, in the regime of expansion (ξ p > 1.5), major effects of system size are apparent. As expected, pores are destabilized in smaller systems due to interactions of the pore with its periodic image.
Specifically, with 128 DMPC lipids, the free energy of the metastable pore is increased as compared to larger membranes, suggesting a reduced lifetime of the pore (Fig. 6, black). In addition, the local minimum of the PMF is shifted to smaller ξ p , demonstrating that the metastable pore with 128 DMPC lipids has a reduced radius.
The PMFs for 200 or 512 DMPC lipids, in contrast, agree around the minimum of the metastable pore, suggesting that a 200 DMPC system is sufficient for capturing the energetics and radius of the metastable state (Fig. 6, blue and magenta). However, upon further expansion of the pore (ξ p > 3.5), the PMF for 512 lipids reveals a constant slope indicating a constant line tension (see above), whereas the slope of the PMF for 200 lipids continuously increases. This behavior suggest that the 200 lipid system cannot be used to compute the line tension of an expanding pore because the system passes directly from (i) the regime with reduced (or even zero) line tension around the metastable state to (ii) a regime with spuriously increased line tension due to periodic boundary artifacts. Table 1: Line tension of the rim of large pores over a POPC membrane in absence or presence of DMSO, obtained from the slope of the PMFs at large pore radii (Fig. 7, dashed lines Case study: Effect of DMSO on the PMF of pore nucleation and expansion Dimethyl sulfoxide (DMSO) is a widely used organic solvent, whose effect on lipid membranes has been studied in great detail by experiments and simulation. DMSO influences the structure, hydration properties, phase behavior, and permeability of membranes. [51][52][53][54][55][56][57][58][59] MD simulations showed that DMSO preferentially binds at the interface between the polar and apolar regions of the membrane, where it may act as a spacer, thereby rendering the membrane thinner and floppier. 60 At high DMSO concentrations, the formation of water defects has been observed in unbiased simulation studies, rationalizing the increased permeability of membranes for electrolytes in presence of DMSO. [60][61][62][63] However, the effects of the DMSO on the free energies of pore formation and expansion are not well understood. Here, we revisit DMSO to illustrate how the proposed RC may be used to rationalize the effects of DMSO on membrane stability and pore formation in energetic terms. Taken together, DMSO shows a multitude of effects on the free energy landscape of pore formation, affecting the rates of pore nucleation and closure, the probability and size of the pore size, as well as the free energy cost for further expansion. The PMFs along ξ p reflect all these effects, hence providing comprehensive understanding of the influence of DMSO on the energetics of membrane stability and pore formation.

Discussion
We have presented a new reaction coordinate (RC) ξ p for pore formation by combining two established RCs, namely (i) the "chain coordinate" ξ ch to follow pore nucleation, and (ii) the pore radius R to follow pore expansion. This was achieved by introducing a continuous switch from one coordinate to the other. Because pulling the system along the new RC ξ p allows for reversible pathways from the flat, unperturbed membrane up to a large pore with a diameter of several nanometers, the RC allows one to compute the free energy of the pore (at any stage of pore formation) relative to the flat, unperturbed membrane.
The capability of new RC was illustrated by computing the PMFs of pore formation with different lipids and with different system sizes. For DMPC membranes, a system of 128 lipids was sufficient for obtaining the correct barrier for pore nucleation, 200 lipid were required to obtain the correct free energy and size of a metastable pore, but 512 lipids were required for reliably obtaining the line tension of the pore rim. Hence, the appropriate system size needed to avoid periodic boundary artifacts depends on the quantity under study. In addition, we showed that DMSO has drastic effects on pore formation, affecting the free energy barrier of pore nucleation, the free energy and metastability of the open pore, and the line tension as given by the slope of the PMF at large pore radii.
Pore nucleation and pore expansion occur along different directions in conformation space.
Therefore, it may appear surprising that restraining only a single degree of freedom, as given by ξ p , is sufficient for steering the system along such different directions. However, the success of the RC used here critically depends on the metastability of the pore. After pore nucleation, and as the radius of the pore starts to increase, the metastability ensures that the transmembrane defect remains intact, that the chain coordinate ξ ch (r) remains close to unity and, consequently, that also the switch function H (ξ ch (r) − ξ s ch ) remains close to unity. Therefore, while restraining the open pore along ξ p , the switch function H does not fall back to small values, and we obtain a stable restraint along the radius R. As such, the design of the new RC makes use of the specific free energy landscape of the pore formation.
We anticipate that PMF calculations with highly unstable pores might require an additional restraint to avoid a depletion of the part of the pore during pore expansion. For instance, a second restraint keeping ξ ch (r) ≈ 1 could maintain a continuous transmembrane defect.
The molecular mechanism of pore formation observed here is largely compatible with previous studies, suggesting that pore formation by pulling along ξ p does not perturb the membrane more strongly than needed for pore formation. In other words, this suggests that is the simulations follow approximately the minimum free energy pathway. As visualized in Movies S1 and S2 and by Fig. 4, pore formation starts by a local thinning of the membrane, that is by an elastic indentation of the headgroup/tail interface. 14 Next, a thin water needle penetrates the hydrophobic core. 26 The water needle widens by penetration of additional water and a substantial number of lipid headgroups. Interestingly, the density analysis in Movie S2 reveals that small pores with a radius of 4-7Å contain more headgroup than water density. This observation is compatible with an important role of the headgroups for crossing the transition state of pore formation, as suggested by a recent study. 14 Larger pores with a radius of more than 12 Å reveal a nearly ideal toroidal shape, in which the head groups shield the membrane-spanning water defect from the hydrophobic lipid tails. However, the density of the headgroups is lower in the toroidal pore as compared to the flat membrane, reflecting a sub-optimal lipid packing at the pore rim (Movie S2).
The RC proposed here will be useful for a range of applications. It will be interesting to quantify the influence of antimicrobial or cell-penetrating peptides on pore formation, as required for a rational design of such membrane-active peptides. 64 PMF calculations with the new RC could reveal the free energy landscape of electroporation. 20 The RC may be further applied to quantify the stability of liposome-based drug carries, and to model the controlled drug release from such carriers. 65 Routine calculations of the pore's line tension may be useful for further refinement of lipid force fields. Indeed, agreement of lipid membrane simulations with experimental line tension data would be important for simulations that involve topological transitions of membranes, during which membranes adopt non-planar geometries. In the context of membrane fusion, PMFs could rationalize the role of the lipid composition and of transmembrane helices during fusion pore opening and expansion. 66 In addition, we anticipate that the RC may be used to study stalk formation and stalk widening, during which the hydrophobic membrane cores carry out similar topological transition as the water phase during pore formation and expansion. 3,67 zero if R(r) = 0. The gradient of the number of polar atoms is where θ is the derivative of θ.
Modification of the RC and simulation setup for highly floppy mem-

Acknowledgement
This study was supported by the Deutsche Forschungsgemeinschaft via grant SFB 803/A12.

Supporting Information Available
The following files are available free of charge.
• Movie S1: 62 simulation snapshots of a membrane of 512 POPC lipids during pore formation, taken from the final MD frames of umbrella sampling windows.
• Movie S2: density of head groups, water, and lipid tails, averaged over each of 62 umbrella windows of a membrane of 512 POPC lipids. The average ξ p and the radius R of pore are shown in the top right corner. The density is plotted in color code as function of the lateral distance r and vertical distance z from the pore center.