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, a new RC capable of driving both pore nucleation and pore expansion in lipid membranes is presented. The potential of mean force (PMF) along the RC computed with molecular dynamics 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. The RC is illustrated by quantifying the effects of (i) simulation system size and (ii) the addition of dimethyl sulfoxide 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.


■ 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 membrane pores is often induced by applying electric fields, a process termed electroporation. 6,7 In biophysical experiments, pores were also induced by applying osmotic stress or membrane tension. 8−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 transitions, 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−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 a function of the pore radius R as follows (1) 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 nonequilibrium conditions, where the pore was induced by applying an external perturbation, such as an electric field or mechanical stress. 17−21 From such nonequilibrium 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 the lateral direction away from the pore center, 24 a lipid flip−flop coordinate, steering a single lipid head group toward the membrane core, 13,25−27 or the mean solvent density within a membrane-spanning cylinder. 28 Complementarily, a theory for pore formation based on continuum elasticity was recently proposed, which quantified the pore free energy in terms of the pore radius together with the thickness of a hydrophobic pore rim. 29,30 In a systematic comparison of RCs of pore formation, it was found that pulling along established RCs for pore nucleation efficiently triggered pore formation. However, if the open pore is metastable (long-living), pulling along these RCs in the 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. 31 Such different pathways during forward and reverse pulling may lead to severe hysteresis problems during PMF calculations. These 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. It was found that pulling along ξ ch allows for reliable pore opening and closing, thereby avoiding any hysteresis problems during PMF calculations. 11,32 Using the chain coordinate, a nucleation barrier, and therefore, the metastability of the open pore in certain membranes, were resolved for the first time, which had been proposed 40 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 arrange-ment 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 closure. 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. It is shown 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. The RC is illustrated by analyzing the effects of membrane size and dimethyl sulfoxide (DMSO) on the free-energy landscape of pores.

■ METHODS
Definition of the RC. The new RC is defined as a combination of two coordinates. First, the state of pore nucleation was quantified with the chain coordinate ξ ch (r) defined in ref 32, 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 Figure 1A. ξ ch (r) is approximately given by the fraction of slices that are filled by polar heavy atoms ( Figure 1A, blue shaded slices). As polar atoms contribute to ξ ch , typically, oxygen atoms of water and lipid phosphate groups were used. Hence, both water and head groups contribute to the degree of connectivity that is quantified by ξ ch . It was observed previously that using additional polar atoms for defining ξ ch , such as oxygen atoms of ester groups, had only a small effect on the PMFs. 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". Whether a slice was occupied or 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). Oxygen atoms of water and phosphate groups are considered as polar atoms contributing to the RC. (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). not was defined with two differentiable functions: (i) the number of polar atoms within a slice was defined by a differentiable indicator function that takes unity inside the slice and zero outside, yet with a differentiable boundary; (ii) the number of polar atoms was translated into the occupancy ∈[0, 1] with a smooth switch function, which strongly increases (e.g., by 0.75) upon the addition of the first polar atom and only slightly increases upon the addition of the second or third atom. 32 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 the pore expansion. Instead, pulling along ξ ch may efficiently open or close transmembrane defects.
Second, the state of pore expansion was quantified via the approximate radius of the pore R(r), which is defined via the number of polar heavy atoms n P (r) within a horizontal layer of thickness D in the center of the membrane (parallel to the membrane plane, see Figure 1B). Because the defect takes approximately the shape of a cylinder, R(r) is defined as (2) D = 1 nm is used 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 atoms n P inside the hydrophobic membrane core; the latter quantity is merely translated into an approximate pore radius to obtain a more intuitive interpretation of the RC. Similar to previous studies, 28,32 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 |x| < 1 − h, where h = 0.1. For the definition of θ(x; h) and its derivative, one shall refer to previous work. 32 Z mem denotes the center of mass of a reference group, here defined by all heavy atoms of the lipid tails.
These two coordinates for pore nucleation (ξ ch (r)) and pore expansion (R(r)) are combined 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 third-order polynomials as follows Hence, H ε [ξ ch (r) − ξ ch s ] switches gradually between zero and one as the chain coordinate ξ ch (r) passes a switching value ξ ch s .
ξ ch s = 0.925 is used to characterize a continuous polar defect and ε = 0.05 is used. 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) ≈ ξ ch s , 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) ≲ ξ ch s 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 , R 0 is chosen as the average value of R(r) in the switch region, that is 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, a reasonable estimate for R 0 can be obtained as follows, thereby avoiding the need for such a prior simulation. Let the chain coordinate ξ ch be defined by a cylinder with height N s d s , composed of N s slices with thickness d s . The thickness D of the central layer is smaller than the height of the cylinder, D < N s d s ; hence N s,i = D/d s and N s,o = N s − D/d s cylinder slices are located inside and outside of the central layer, respectively. Then, the number of polar atoms in the central layer (and thereby the value of R(r)) that leads to ξ ch ≈ ξ ch s is estimated. To this end, it is assumed 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, polar atoms are randomly added to the N s,i slices inside the central layer, until ξ ch ≈ ξ ch s . This procedure is repeated 10,000 times at the beginning of the simulations. Finally, the average number of polar atoms n ̅ P added to the central layer was used to estimate R 0 via eq 2.
This estimate is found to be in good agreement with the average R(r) in simulations with ξ ch ≈ ξ ch s . 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) − ξ ch s ]⟩ = 0.50 (indicating ⟨ξ ch (r)⟩ = ξ ch s ), it is found that ⟨R(r)⟩ = 0.448, close to the value of 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 ξ ch s . Implementation. The implementation of the RC ξ ch was taken from our previous study, 32 but merged with Gromacs 2018.8. 33 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 approximately 7% slower as compared to a free simulation on a server equipped with a six-core Intel Xeon E-2136 and an Nvidia GeForce GTX 1070Ti graphics card.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Simulation Setup and Parameters. The membrane systems were set up with the MemGen webserver (http:// memgen.uni-goettingen.de). 34 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 modeled with the force field by Berger et al. 35 and the SPC water model was used. 36 Parameters for DMSO were taken from ref 37. Bonds and angles of water were constrained with the SETTLE algorithm. 38 All other bonds were constrained with LINCS. 39 The temperature was controlled at 300 K using velocity rescaling and by coupling water and lipid to separate heat baths. 40 During equilibration, the pressure was controlled at 1 bar using a semi-isotropic weak coupling scheme (τ = 1 ps). 41 Although weak coupling does not yield a well-defined ensemble, it is used here owing to its numerical stability. Electrostatic interactions were calculated using the particle-mesh Ewald method. 42,43 Dispersion interactions and short-range repulsions were described by a Lennard-Jones 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. 44 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 or 120 ns, respectively. If not stated otherwise, pulling simulations were carried out in the opening direction, that is from the state with a flat membrane to the open pore. To test for undesired hysteresis effects, starting frames for the POPC membrane were alternatively taken from a pulling simulation in the closing direction.
For PMFs along ξ ch , 24 umbrella positions were selected between 0.065 and 0.625 in steps of 0.08, and between 0.7 and 1.0 in steps of 0.02. For windows restrained to ξ ch ≤ 0.7, an umbrella force constant of k = 3000 kJ mol −1 was used, and k = 5000 kJ mol −1 otherwise. For PMFs along ξ p , 63 umbrella positions were selected between 0.065 and 0.625 in steps of 0.08, between 0.68 and 1.10 in steps of 0.03, and between 1.12 and 6.97 in steps of 0.15. 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, 45 as implemented in the gmx wham module of Gromacs. 46 Statistical errors were estimated using 50 rounds of Bayesian bootstrapping of histograms. 46 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. 32 The number of slices N s was chosen such that ξ 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, a PMF for pore nucleation along the ξ ch coordinate was computed using two setups: (a) a 4 fs time step with a stochastic dynamic integration scheme 47 and (b) a 5 fs time step and the temperature coupling indicated above. Because the PMFs were virtually identical, the computationally more efficient setup [option (b)] was used.
■ RESULTS AND DISCUSSION Fixed Versus Mobile Transmembrane Cylinder. Before turning to the new RC, the influence of the lateral mobility of the transmembrane cylinder was investigated by comparing two setups: first, with a cylinder at a fixed lateral position and second, with a cylinder allowed to follow the defect as the defect travels laterally in the membrane plane. As outlined in the Methods, the chain coordinate ξ ch is defined by a transmembrane cylinder spanning the two polar head group regions and the hydrophobic membrane core. Using such a cylinder is critical because it avoids that two laterally displaced partial defects, one connected with the lower and one connected with the upper head groups, are misinterpreted as a continuous membrane-spanning defect. In other words, the cylinder enforces that the formed defect is localized in the membrane plane.
Previously, it was found that it is critical to define the lateral position of the cylinder dynamically, allowing the cylinder to follow the defect as it travels laterally in the membrane plane. Otherwise, if the cylinder would be fixed, the system could evade the unfavorable transition state by shifting the defect laterally out of the cylinder, whereupon the nucleation barrier would be integrated out (Figure 2, flat black curve at ξ ch ≈ 0.85). Using a mobile cylinder avoids these problems, thereby clearly resolving the nucleation barrier (Figure 2, magenta).
With a small simulation system of 128 lipids, simulations with a mobile cylinder are numerically stable. However, upon Figure 2. Importance of using a mobile transmembrane cylinder: PMFs of pore nucleation over a membrane with 128 DMPC lipids, computed along the chain coordinate ξ ch . The PMF was computed with the transmembrane cylinder fixed in all umbrella windows, where the nucleation barrier was integrated out (black), with the cylinder fixed in none of the umbrella windows (magenta), or fixed only in windows whose umbrella center is <0.7 (blue). The difference (blue vs magenta) is rationalized by increased lateral entropy of the defect with a laterally mobile cylinder, as indicated by the label −TΔS lateral ; see the text for details. Shaded areas indicated standard errors computed by bootstrapping.
simulating pore formation in large membrane patches, it was found that using a mobile cylinder in umbrella windows with a nearly flat membrane may lead to integration problems, merely because water fluctuations in the glycerol region cause rapid lateral movements of the cylinder. To solve such numerical problems, the cylinder was fixed in umbrella windows for which the reference position was ξ ch < 0.7, thereby avoiding integration problems, while leaving the cylinder mobile in windows for which the reference position was ξ ch ≥ 0.7, thereby properly resolving the nucleation barrier. As shown in Figure 2 (blue and magenta curves), fixing the cylinder at small ξ ch leads to a slightly increased nucleation free energy by approximately 8 kJ/mol, rationalized by the entropic cost for forming the initial defect at a lateral position near the fixed cylinder center, as compared to forming the defect at an arbitrary lateral position.
In addition, the position of the PMF minimum is shifted to a larger ξ ch when using a mobile cylinder as compared to a cylinder that is fixed at small ξ ch (Figure 2, magenta and blue, ξ ch ≈ 0.25). This effect is merely a manifestation of small water fluctuations at the water/hydrophobic interface in a flat unperturbed membrane. Namely, as the lateral position of the cylinder is defined to follow any water penetration into the membrane, 32 the cylinder also jumps onto such small water fluctuations into the hydrophobic core. In consequence, more cylinder slices are on average water-filled, leading to a larger ξ ch for the flat membrane. However, this effect has no physical significance and could have been removed by choosing a smaller number of slices N s when using a mobile cylinder. Taken together, the overall shape of the PMF remains similar. Hence, fixing the cylinder at small ξ ch allowed numerically stable simulations of pore nucleation in both small and large membranes.
Single PMF for Nucleation and Pore Expansion. Figure  3 presents PMFs covering both pore nucleation and pore expansion across membranes of 512 DMPC or 512 POPC lipids, computed along the new RC ξ p . During umbrella sampling, the lateral position of the transmembrane cylinder in windows with a reference position of ξ p < 0.7 was fixed, and the cylinder was mobile in all other windows, equivalent to the setup used for Figure 2 (blue curve). Typical simulation snapshots, taken from the final frames of various umbrella windows, reveal the progression from a flat membrane, via a partial defect with a locally thinned membrane due to indented head groups, a thin water needle, up to a fully formed large pore ( Figure 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 ξ p positions, Movie S2 presents the average densities of head groups, tails, and water during all windows. Evidently, the PMFs reveal an approximately quadratic regime (ξ p < 0.85), corresponding to the indentation of the head groups and penetration of a water needle into the membrane. Broad maxima in the PMFs around ξ p ≈ 1 indicate the transition states of pore formation, followed by shallow, broad minima around ξ p ≈ 2, indicating metastable open pores. It was found previously that, despite the only shallow minima, these metastable pores can be stable for 10 μs or longer. 31 During further expansion of the pores, the PMFs reveal a linear regime at ξ p > 4, which is compatible with a constant line tension of the pore rim.
To exclude that the PMFs are compromised by undesired hysteresis effects, the PMF was recomputed for POPC from a set of umbrella sampling simulations with starting conformations taken from a pulling simulation in the reverse direction, that is, from the open pore toward the flat membrane (Figure 3, green dashed line). Evidently, the differences between the PMFs from opening and closing simulations are relatively small, implying reasonably well converged PMFs and the absence of major  hysteresis effects. However, for truly qualitative convergence, slightly longer umbrella simulations as compared to the 100 ns windows used here may be required.
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 its 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 that of POPC. The radius is given by ξ p R 0 at the shallow PMF minimum, where the parameter R 0 was 0.40 and 0.38 for DMPC and POPC, respectively (see Methods), implying the radii of the metastable pores of 0.80 and 0.72 nm for DMPC and POPC, respectively.
Third, the slope of the PMF in the linear regime (ξ p > 4) reveals the line tension of the pore edge. Classical nucleation theory models the free energy of the open pore by an unfavorable contribution from the line tension of the pore edge plus a relief of surface tension between membrane and water (eq 1). The term owing to surface tension vanished in the simulations because simulations were performed (i) without an externally applied surface tension in the membrane plane but instead at a constant pressure and (ii) with a constant number of lipids, hence the growth of the pore leaves the area of the membrane−water interface nearly constant. Fitting a line to the PMFs in the interval ξ p ∈ [4, 6.6] yields line tensions of (21.83 ± 0.04) and (27.69 ± 0.04) pN for DMPC and POPC, respectively. These lines tensions are slightly smaller than the values (23.3 ± 0.6) and (35.5 ± 0.05) pN for DMPC and POPC, respectively, reported by West et al. 48 for pores in membranes modeled with the Berger lipids, as also used here. The slightly increased line tensions by West et al. might originate (i) from a different definition for the pore radius or (ii) from simulating a smaller membrane patch of 256 lipids, as compared to the patch of 512 lipids used here. In a previous study, likewise, it was observed that the line tension in a smaller DMPC system with 200 lipids steadily increased up to ∼30 pN upon growing the pore. 11 A detailed analysis of the effect of system size on the PMFs is presented below.
Additional support for the PMFs and for the presence of metastable long-living pores if given by free simulations. Figure  5A,B presents five independent replica simulations of the DMPC and POPC membranes, starting from a wide open pore (ξ p = 7). As expected from the PMFs in Figure 3, the pores rapidly shrink in free simulations, until the pore becomes metastable, corresponding to the local free-energy minimum in the PMFs at ξ p ≈ 2 (cf. Figure 3). Evidently, among the simulations of 250 ns, none of the pores close, demonstrating that even the shallow free-energy minimum leads to pore lifetimes on the order of at least several 100 ns.
Taken together, with a single PMF along the new RC, the full free-energy landscape of pore nucleation and pore expansion was obtained. The PMF displays the free-energy cost for nucleating the pore, the metastability, size, and free energy of the open pore, and the line tension of the pore edge.
Switching the RC from Nucleation to Expansion. To illustrate how the RC ξ p accounts for the switch from pore nucleation to pore expansion, Figure 6 presents the progression of the most relevant parameters during nucleation and expansion. The values in Figure 6 report the respective mean . Five replicas of 250 ns were carried out for each lipid. As expected from the PMFs, the pores rapidly shrink until the simulation reaches the metastable pore size (ξ p ≈ 2), corresponding to the local free-energy minimum (compared with PMFs in Figure 3). (C,D) Histograms of ξ p taken from the last 100 ns of all replicas. and standard deviation taken from umbrella windows of the system with 512 POPC lipids. Evidently, during pore nucleation (ξ p ≲ 0.9), the new RC ξ p (red symbols) is simply given by the chain coordinate ξ ch (black symbols) because the switch function H ε (green symbols) still takes the value of zero (see also eq 4). As the chain coordinate approaches the switch value ξ ch s set to 0.925 in this simulation, the smoothed Heaviside step function H ε (ξ ch − ξ ch s ) switches from zero to unity and, hence, the new RC ξ p starts to deviate from the chain coordinate ξ ch (see eq 4, and Figure 6, black/red/green symbols). Henceforth, ξ ch saturates at unity (indicating a continuous water defect), whereas ξ p continues to increase as (R − R 0 )/R 0 (orange symbols).
In addition, Figure 6 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, Figure 6 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.
Simulation Box Size Influences the Size and the Free Energy of the Open Pore. Previous simulations investigated the pore formation in the membrane of various sizes, spanning systems between 64 and 512 lipids. 17,20,24−26,28,31,32,49−52 However, it was previously found that the system size may strongly influence the free-energy landscape of pore formation. In small membrane systems, the open pore may be destabilized because there is no flat membrane between the pore and its periodic image. 31 To probe the influence of system size on the free-energy landscape of pore formation, the PMF was computed along the newly proposed RC using systems of 128, 200, or 512 DMPC lipids. As shown in Figure 7, 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 because the membrane between the pore and its periodic image becomes overly bent.
Specifically, with 128 DMPC lipids, the free energy of the metastable pore is increased as compared to that of larger membranes, suggesting a reduced lifetime of the pore (Figure 7, black). In addition, the local minimum of the PMF shifted to smaller ξ p , demonstrating that the metastable pore with 128 DMPC lipids exhibits 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 (Figure 7, 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 suggests 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.
Case Study: Effect of DMSO on the PMF of Pore Nucleation and Expansion. 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. 53−61 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. 62 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 the presence of DMSO. 62−65 However, the effects of DMSO on the free energies of pore formation and expansion are not well understood. Here, DMSO is revisited to illustrate how the proposed RC may be used to rationalize the effects of DMSO on membrane stability and pore formation in energetic terms. Figure 8 presents PMFs of pore nucleation and expansion in a membrane with 512 POPC lipids plus increasing amount of DMSO in POPC/DMSO ratios of 2:1, 1:1, 1:1.5, and 1:2. Evidently, the addition of DMSO has multiple, drastic effects on  the PMF. (i) The pore nucleation free energy decreases by up to 20 kJ/mol, indicating that the rate of pore opening is increased by a factor of ∼3000; (ii) the free energy of the open, metastable pore decreases by up to 40 kJ/mol, demonstrating a greatly increased probability for an open pore. (iii) As a consequence of points (i,ii), the free energy barrier that must be overcome to close the pore increases with DMSO, implying greatly decreased rate of pore closing (or increased pore lifetime) by a factor of ∼3000. Hence, the pore becomes more metastable. (iv) The local PMF minimum corresponding to the open, metastable pore is shifted to larger ξ p , implying larger pore radii. However, because the minima are shallow, the pores likely adopt a wide distribution of pore radii. (v) The slope of the PMFs at large ξ p decreases, reflecting a decreased line tension of the pore rim ( Figure 8, dashed lines). The line tensions obtained from the slope of the PMFs are summarized in Table 1. The larger radii of the metastable pores at increased DMSO content (point iv) are likely a consequence of the smaller line tension. In addition, the reduced line tension affects the shape of large open pores. Whereas pores adopt approximately circular shapes at low DMSO content ( Figure 9A−C), as favored by a high line tension of the rim, pores at high DMSO content adopt increasingly irregular shapes ( Figure 9D,E). This trend was confirmed by visual inspection of many more MD frames. (vi) Finally, the minima at ξ p ≈ 0.3 corresponding to the flat, unperturbed membrane are shifted to slightly larger ξ p . This shift reflects that, in the flat membrane, more slices of the transmembrane cylinder ( Figure 1A) are filled by polar atoms; hence, the flat membrane is thinner in the presence of DMSO. 62 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 a comprehensive understanding of the influence of DMSO on the energetics of membrane stability and pore formation.

■ DISCUSSION
A new RC ξ p is presented 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 different system sizes. For DMPC membranes, a system of 128 lipids was sufficient for obtaining the correct barrier for pore nucleation, 200 lipids 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, DMSO is shown to have drastic effects on pore formation, affecting the free-energy barrier of pore nucleation, the free energy, metastability, and radius 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, for the following reason: after pore nucleation, as the radius of the pore starts to increase, the metastability ensures that the transmembrane defect remains intact, the chain coordinate ξ ch (r) remains close to unity and, consequently, the switch function H ε (ξ ch (r) − ξ ch s ) also remains close to unity. Therefore, while restraining the open pore along ξ p , the switch function H ε does not fall back to small values, and a stable restraint is obtained along the radius R. With a highly unstable pore, in contrast, while restraining purely the number of polar atoms in the membrane plane, the membrane might evade the high free-energy state of a continuous transmembrane pore and prefer to fall back into a state with two disconnected partial defects. Indeed, such problems, where the membrane evades states of high free energy by moving orthogonal to the restraint, render PMF calculations of topological transitions generally challenging. As such, the design of the new RC makes use of the specific freeenergy landscape of pore formation and thereby avoids such problems. It can be anticipated that, for PMF calculations with highly unstable pores, alternative strategies might be advantageous. For instance, expanding the pore with a repulsive cylindrical potential acting on the lipid tail atoms might be less prone to undesired effects as compared to the attractive potential acting on polar atoms used here. 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 the simulations follow approximately the minimum free-energy pathway. As visualized in Movies S1 and S2 and by Figure 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 the 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 headgroups 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 headgroups is lower in the toroidal pore as compared to that of the flat membrane, reflecting a suboptimal 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. 66 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 carriers and to model the controlled drug release from such carriers. 67 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 the membranes adopt nonplanar geometries. In the context of membrane fusion, PMFs could rationalize the role of the lipid composition and transmembrane helices in fusion pore opening and expansion. 68 In addition, it can be anticipated that the RC may be used to study stalk formation and stalk widening, during which the hydrophobic membrane cores carry out a similar topological transition as the water phase during pore formation and expansion. 3,69 ■ APPENDIX Gradients of ξ p (r) Restraining the simulation system along ξ p (r), as carried out during umbrella sampling simulations, requires the calculation of the gradient ∇ i ξ p of the RC with respect to the Cartesian coordinates of atom i. From eq 4, one obtains Here, H ε ′ is the derivative of the smoothed Heaviside step function. The gradient of the chain coordinate ∇ i ξ ch was computed as described previously. 32 The last term of eq 7 was set to zero if R(r) = 0. The gradient of the number of polar atoms is i k j j j j y { z z z z i k j j j j j j j j j y  (8) where θ′ is the derivative of θ.

Modification of the RC and Simulation Setup for Highly Floppy Membranes
By default, the radius of the pore was defined based on the number of polar atoms within a horizontal slab at the membrane center (see Methods and Figure 1B, orange shading). For large, floppy membranes, such as membranes with 512 lipids in the presence of DMSO, this definition leads to problems because polar atoms may enter the slab not only by expanding the pore but also by membrane undulations. In such cases, the horizontal slab was replaced with a cylinder. The axis of this cylinder coincided with the axis of the membrane-spanning cylinder used for defining the chain coordinate ξ ch . The radius was chosen to be 0.8 nm larger than the current approximate radius of the pore given by ξ p ref R 0 , where ξ p ref is the reference position of the umbrella potential along ξ p . This definition reduced the undesired effects from membrane undulations, while at the same time minimized undesired effects of the cylinder on shape fluctuations of the pore.
For the floppiest membranes (512 POPC plus 768/1024 DMSO), an additional measure was needed to avoid problems from heavy membrane undulations. Here, flat-bottomed position restraints are used on z-coordinates for all lipid tail atoms, V fb ( Here, z i is the z-coordinate (membrane normal) of tail atom i, Z mem is the center of the membrane, z fb = 2.2 nm is the half-width of the flat region, k = 100 kJ/mol/nm 2 is the force constant, and H is the Heaviside step function. This potential allowed natural local fluctuations of head groups but excluded large-scale membrane undulations, thereby rendering the simulations more stable.

* sı Supporting Information
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.0c01134. 62 simulation snapshots of a membrane of 512 POPC lipids during pore formation, taken from the final MD frames of umbrella sampling windows (MP4) Density of head groups, water, and lipid tails, averaged over each of 62 umbrella windows of a membrane of 512 POPC lipids, which is plotted in color code as functions of the lateral distance r and vertical distance z from the pore center (MP4)

Notes
The author declares no competing financial interest.

■ ACKNOWLEDGMENTS
The author thanks Greg Bubnis for stimulating discussions and for critical reading of the manuscript. This study was supported by the Deutsche Forschungsgemeinschaft via the grant SFB 1027/B7.