Sequential Water and Headgroup Merger: Membrane Poration Paths and Energetics from MD Simulations.

Membrane topology changes such as poration, stalk formation, and hemifusion rupture are essential to cellular function, but their molecular details, energetics, and kinetics are still not fully understood. Here, we present a unified energetic and mechanistic picture of metastable pore defects in tensionless lipid membranes. We used an exhaustive committor analysis to test and select optimal reaction coordinates and also to determine the nucleation mechanism. These reaction coordinates were used to calculate free-energy landscapes that capture the full process and end states. The identified barriers agree with the committor analysis. To enable sufficient sampling of the complete transition path for our molecular dynamics simulations, we developed a "gizmo" potential biasing scheme. The simulations suggest that the essential step in the nucleation is the initial merger of lipid headgroups at the nascent pore center. To facilitate this event, an indentation pathway is energetically preferred to a hydrophobic defect. Continuous water columns that span the indentation were determined to be on-path transients that precede the nucleation barrier. This study gives a quantitative description of the nucleation mechanism and energetics of small metastable pores and illustrates a systematic approach to uncover the mechanisms of diverse cellular membrane remodeling processes.

Significance

The primary steps and nucleation of lipid membrane pore formation are key to membrane fusion, viral infection, and vesicular cellular transport. Despite decades of experimental and theoretical studies, the underlying mechanisms are still not fully understood at the atomic level. Using a committor-based reaction coordinate and molecular dynamics simulations, we report structural and energetics insight into the full poration process. We find that the pore nucleates via an elastic indentation rather than by forming a hydrophobic defect. Subsequently, water pierces the thinned slab as a prerequisite for the following axial merger of the first lipid headgroups from opposite monolayers, which best characterizes the transition state. We also identify a metastable pore basin, thereby explaining previous indirect experimental evidence.

Introduction

Lipid membranes can undergo a number of topological remodeling processes—such as endo- and exocytosis, vesiculation, viral entry, and fertilization—that involve the formation or closure of aqueous pore defects. Pore nucleation is key to transmembrane transport (1), ion permeation (2,3), antimicrobial peptide function (4,5), and bilayer equilibration (6,7) and is the essential step of synaptic transmission (8, 9, 10). Despite the evidence for small, metastable aqueous pores from conductivity (11) and tension (12) experiments, as well as molecular simulations (13, 14, 15), the atomistic details, driving forces, and kinetics of their formation and closure are not entirely understood (16,17).

In the long-established continuum description for pores, namely the Litster model (18), a pore’s energy and radius are governed by the balance of membrane surface tension and the pore rim line tension. This model is suited for already-formed pores but does not explain the large energies required to create the topological defect, nor can it describe metastable pore defects (sometimes called “prepores”) in the absence of applied tension (11, 12, 13,15). In a more recent continuum description, a nascent pore is treated explicitly as a hydrophobic cylinder of solvent-exposed lipid tails, with both a height and radius (19). Here, we focus on the formation of metastable pores that arise in a tensionless bilayer, using molecular dynamics (MD) simulations to fully resolve the process. We use the Berger lipid force field (20), which models lipids atomistically except for an implicit description of hydrogens on nonpolar heavy atoms.

Specifically, we will address the following questions. Is the poration pathway hydrophobic, with a penetrating water column surrounded by lipid tails, or hydrophilic, involving an indentation where lipid headgroups submerge to shield water-tail interactions? In fact, “water wires” that span the hydrophobic slab have been observed in several recent simulation studies (13,21,22), but whether or not this step is rate limiting is unclear. In one study, the energetics of pore formation were found to be insensitive to “bundling” four waters together, suggesting that the precise organization of water matters little (23). We will therefore address the question of whether or not these water wires are the energetic transition state and whether single-file water columns are sufficient to nucleate the pore. To address these questions, we aim to identify optimal reaction coordinates (RCs) in terms of collective variables (CVs) to properly describe the progress and energetics of the nucleation mechanism. (The terms CV and RC are both used throughout this study, with the intention that CV indicates a collective coordinate that may or may not be a useful RC, whereas an RC is tasked to measure reaction progress).

Whereas spontaneous pores can form on microsecond timescales for simulations of short-tailed lipids at elevated temperatures (e.g., (21) using simulations at 50°C), biasing potentials acting along a chosen RC are usually required to access relevant structural intermediates, in particular the transition state (TS). Typically, an ad hoc RC ϕ(x)—defined as a function of the combined coordinate vector x of a suitable set of relevant atoms—is used as a biasing coordinate to compute a free-energy profile (or potential of mean force, PMF) G(ϕ) via umbrella sampling or related methods (24,25). The proper choice of RCs is therefore critical for obtaining well-converged energetics and also serves to characterize the reaction mechanism in structural terms. Indeed, for diverse membrane remodeling processes (stalk nucleation, vesiculation, hemifusion), the fluid disorder of the membranes makes it far from obvious which collective motions (i.e., RCs) are best suited as mechanistic descriptors and for biasing simulations toward transient intermediates.

A recent comparison of poration RCs documented that many established and intuitively plausible RCs for poration were poor biasing coordinates because of prohibitively slow convergence and resulting hysteresis effects (17). Further, and quite generally, even for well-converged simulations, suboptimal choice of RCs may artificially lower or hide important energy barriers if different intermediates of the reaction progress are projected onto similar values along the RC (3,26). In other cases, the chosen RC could introduce spurious barriers. In general, the TS region may not be properly described by the RCs, and the position of the obtained PMF barrier may differ markedly from the true TS.

More recently, a poration RC based on the fraction of cylindrical slabs (spanning the membrane) occupied by hydrophilic particles (water, lipid headgroups) has been proposed, which enabled converged free-energy calculations suggesting a free-energy barrier and TS leading to the metastable (pre)pore state (14). Despite this advance, the observed barrier was sensitive to several tuning parameters, such that it remained unclear which combination provides an accurate result, and the subsequent pore could not be well resolved.

For a more systematic approach, we exploit the optimal RC for any conformational transition from a conformational state A to a state B, which is—by definition—the committor c(x) P(B|x), i.e., the probability that an unrestrained trajectory seeded from x reaches state B before state A (cf. Fig. 2 b). This concept traces back to Onsager (27) and has been more recently reviewed (28, 29, 30). In this framework, c(x) optimally reflects the reaction progress, and all configurations x with c(x) = 1/2 comprise the TS between the basins of attraction of A and B. Although G(c) could be considered a “true” free-energy profile, it is practically impossible to compute. Computing c(x) is computationally prohibitive because it involves starting sufficiently many trajectories from every point x along the transition path ensemble. In addition, c(x) is not an intuitive structural descriptor and thus offers no structural insight by itself. The real utility of c(x) stems from the fact that a good RC must be a strong committor correlate in the TS region. It follows, then, that putative RCs can be quantitatively compared by how well they correlate to c(x).

Figure 2

Schematics showing (a) a commitor regression model f, which consists of a linear combination of input features with a logistic activation. (b) A 2D TS model for the CVs ϕ1(x) and ϕ2(x). The two free-energy basins A and B are separated by the f = 0.5 isoline. The three-point MEP estimate ϕMEP = [ϕ0.2, ϕ0.5, ϕ0.8] is drawn as cyan, red, and blue dots. The solid black line is a linear fit to ϕMEP, constrained to pass through ϕ0.5, and the closed and open circles at the ends indicate the directions of the closed and open pore states, respectively. Because a is a linear function of input CVs, all f = const isolines are parallel. To see this figure in color, go online.

Indeed, this concept has been applied to several other systems (31, 32, 33, 34) as well to as a highly coarse-grained lipid system (35). For these systems, RCs were successfully identified that correlated strongly with c(x) and hence were considered close to optimal. Notably, once the expensive computation of c(x) is achieved for a sufficient number of points, any number of candidate RCs can be evaluated and ranked via postprocessing at only a little computational cost and without further simulations. The appeal of this strategy is that rather than relying on intuition, the selection of an optimal RC becomes an optimization problem with a defined cost function. This framework also allows for building more powerful descriptors via linear (32) or nonlinear (31) combinations of the initial CV pool, although overly complex hybrid coordinates (e.g., derived via a deep neural network) bear the risks of overfitting and of becoming uninterpretable. Yet, for a very broad range of systems, the approach provides a principled method to identify RCs that are sufficiently close to optimal and thus provide quantitative insight into reaction paths and energetics.

Here, we apply such a committor-driven scheme to study pore nucleation in a solvated 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipid bilayer, which is known to form metastable pores (17,23). These pore defects are not necessarily stable for all possible membrane compositions. To efficiently sample states in which c(x) ≈ 1/2, we will first develop a biasing scheme using a membrane-embedded “gizmo” that energetically biases the system toward high-energy states—including on and off-path intermediates—and functionally resembles a lipid scramblase (36,37). Using gizmo biasing, multimicrosecond simulations will be used to sample TS crossing events as input for subsequent (unbiased) c(x) calculations. The obtained c(x) estimates will then be used as a regression target to score and rank putative RCs from a diverse combinatoric pool including CVs for water, lipid headgroups, and lipid tails, as well as pairwise combinations thereof. A second set of gizmo-biased simulations will serve to sample the full poration pathway and recover unbiased PMFs projected onto optimal RCs determined from the committor analysis. Taken together, the PMF, along with the optimal RC, will serve to address the above mechanistic questions.

Theory

Biasing potential

To thoroughly sample the poration process, and in particular the TS, we have developed a pseudomolecular “gizmo” embedded within the membrane, much like (but not intended to resemble) a membrane protein. This gizmo uses repulsive particles to impose a potential energy bias onto adjacent lipid and water molecules so as to direct the membrane toward open, transition, or closed pore states. Because this bias potential is known and well-defined, the structural ensemble obtained by using this bias potential can subsequently be reweighted such that an unperturbed ensemble and resulting free energy is obtained. This approach is similar in spirit to umbrella sampling (24) but is more general because it includes additional (gizmo) degrees of freedom, which need proper treatment (see the Supporting Materials and Methods).

To this aim, the gizmo (sketched in Fig. 1 and fully described in the Supporting Materials and Methods and Fig. S1) was designed with two special-purpose structural elements for driving pore opening and closure, respectively. To drive pore opening, we used a “chain” of lipid-tail-repellent (CX) beads (red) that is aligned to the bilayer normal to create a packing defect in the hydrophobic region (light gray) of the membrane (Fig. 1 a). To drive closure, we used a circular “belt” of water-repellent (WX) beads (green) lying in the membrane midplane to constrict and pinch off the water column (blue) (Fig. 1 b). The chain and belt are flexible and are bound to a stiff frame composed of ghost particles (GG), which interact via bonded interactions, but not with any other particles of the simulation system (i.e., no van der Waals or electrostatic interactions).

Figure 1

Pore gizmo designs. (a) A vertical chain of CX particles (red) drives pore opening by creating a packing defect in the hydrophobic slab. (b) A belt of WX particles (green) drives pore closure by pinching off the water column. The belt is drawn with six WX particles for clarity, whereas the actual belt has 12. The faces in (a) and (b) indicate in which compartments the CX/WX particles are noninteracting (smiley face) and purely repulsive (frowning face) such that these particles can occupy one bulk region of the system without introducing any perturbation. (c and d) When combined to form a gizmo, the CX chain and WX belt allow reversible control of pore formation. (c) The “hydrophobic” H-gizmo uses a chain of 11 CX beads to create a hydrophobic packing defect spanning the bilayer. (d) The “indentation” I-gizmo has its three centermost beads switched off so that two partial chains penetrate the membrane from both sides, creating an indentation defect. To see this figure in color, go online.

The gizmo was implemented via the potential energy of the combined system,where V(x) is the force field for the lipids and solvent system, V(y) that for the gizmo, V(x, y) = λCXVCX(x, y) + λWXVWX(x, y) couples the system and gizmo via repulsive CX-WX interactions, and a harmonic restraint V(ξ(y)) = controls the WX belt radius of the gizmo with a collective radial breathing coordinate ξ (ξ is defined in the Supporting Materials and Methods and illustrated in Document S1. Supporting Materials and Methods, Figs. S1–S7, and Tables S1–S3, Data S1. Example Files for Including a Gizmo in a Gromacs Simulation, Document S2. Article plus Supporting Material e). In the above equations, x denotes the (combined) atomic coordinates vector of the lipids and solvent and y the (combined) atomic coordinates vector of the gizmo. The tilde over indicates that the lipid-solvent system is coupled to, and biased by, the gizmo.

The potentials VCX and VWX are the sums of repulsive interactions for the CX and WX atoms, respectively, and were implemented using Lennard-Jones potentials (σ = 1.1 nm and ε = 0.01 kJ/mol). This choice of Lennard-Jones parameters, in combination with the van der Waals cutoff of 1.2 nm (see the Methods), effectively implements a (repulsive) Weeks-Chandler-Andersen potential (38).

The relative strengths of the gizmo’s chain and belt biases were tuned via λCX, which scales the CX-tail repulsion that drives opening, and ξ0, which controls the WX belt radius via a harmonic restraint that is constricted to drive closure (λWX was fixed to 1). Crucially, except for these biasing potential contributions (i.e., V), the beads do not interact with the physical system of lipids and solvent (see Fig. 1). Accordingly, the potential energy remains otherwise unperturbed, and the gizmo potential can be instantly switched off without creating packing artifacts, as is required for the committor simulations.

This setup also enabled us to derive PMFs for pore opening and closure for the lipid-solvent system from the biased MD ensemble. The PMF recovery scheme, as fully described in the Supporting Materials and Methods, treats the membrane (+ solvent) and gizmo as two coupled systems, recovers their joint probability distribution via the weighted histogram analysis method (39), marginalizes over nuisance and coupling terms, and then uses maximal likelihood to recover the independent, unbiased distribution for the poration CVs of interest.

The gizmo has a few additional features worth noting. First, in the open state, it functions as a lipid scramblase, stabilizing a pore indefinitely and allowing lipids to flip-flop and equilibrate the leaflet asymmetry. Also, because of its shape and hydrophobicity, the WX belt keeps the gizmo embedded and properly oriented in the membrane interior. Lastly, the gizmo’s central atom (a GG atom) provides a convenient reference position for the pore center, which is useful for deriving localized CVs for the lipids and solvent.

We used two different gizmo potentials to efficiently sample two different pore formation paths. The “hydrophobic” H-gizmo (Fig. 1 c) promotes a hydrophobic defect—in which a narrow water column pierces and hydrates the hydrophobic slab and then headgroups pivot inwards—by a chain of 11 CX beads, which creates a tail packing defect fully spanning the membrane. In contrast, the “indentation” I-gizmo (Fig. 1 d) drives indentation and thinning of the bilayer without directly perturbing lipid chain packing at the bilayer midplane. This effect is achieved by a chain that has its three centermost CX beads replaced by GG beads (acting as spacers) such that two CX chain segments, with four beads each, penetrate from opposite sides. The two gizmos also served to assess how much the resulting RCs and PMFs depend on the particular properties of the gizmos.

We used independent ensembles of gizmo-biased simulations to collect input structures for committor analysis and to sample the full poration pathway for PMF calculations. As summarized in Table 1, these ensembles used different initial structures to test convergence, and different gizmo types to allow us to compare hydrophobic and indentation mechanisms. For all of the computed ensembles, we ensured sufficient TS sampling by using extended simulations (up to 2 μs when necessary) targeted to the TS region. This ensured that the PMF calculations are well converged, irrespective of starting structure, and that a sufficiently large number of statistically uncorrelated starting structures is available in the vicinity of the TS for the subsequent committor analysis. Full details of these ensembles are provided in Figs. S2–S4.

Table 1

Simulation Ensembles Used in this Study

Ensemble Type	Gizmo	Starting State
PMF	H	A (closed)
PMF	H	B (open)
PMF	I	A (closed)
PMF	I	B (open)
Committor	H	A (closed)
Committor	I	A (closed)

CVs to describe pore formation

To obtain structural insight into how lipids and solvent reorganize when crossing the TS, we assembled a set of CVs that capture diverse structural changes during pore nucleation. The pool of CVs was chosen to quantify the merger and depletion of specific atom groups α (see Table 2) relative to the pore center, defined by the center y0 of the gizmo (i.e., providing a local coordinate system). To study how different groups of atoms collectively reorganize, we defined atom groups that combine and isolate individual components of the system (e.g., lipid headgroups with and without water oxygen). To compare the effects of local versus nonlocal collective motions, we varied the cutoff number N of atoms closest to y0 used to compute the CVs. Finally, to distinguish different (an)-isotropic symmetries, we implemented isotropic, axial (z axis), and lateral (xy plane) collective coordinates.

Table 2

Atom Groups, α, Used for r, rxy, and rz CVs

Group Label (α)	Group Description
OW	water oxygens
P	lipid phosphorus
NP	lipid nitrogen and phosphorus
NPO	NP ∪ OW
CT	(terminal) lipid tail carbons
C	(all) lipid carbons

α, atom subsets.

The isotropic mean radius was computed aswhere is the position of atom i belonging to the atom group α.

To measure lateral merger or spreading of atoms, the “xy” variant was computed viawhere ‖‖ denotes the distance between points projected onto the xy plane.

The axial variant was computed by projecting and sorting atoms on the membrane-normal z axis to capture the penetration of the lipid slab. When all system coordinates are projected onto the z axis, thermal fluctuations can create nonlocal projection artifacts, giving the false appearance of membrane penetration. These include “hanging droplets” (described in (14)), in which two laterally distant indentations partially penetrate the bilayer, and low-wavenumber thermal undulations (26). For the current system of a 128-lipid patch, these artifacts are expected to be small. However, to avoid these effects and have a CV definition that would be robust to larger systems, we preselected the 100 atoms closest to y0 before projection and sorting.

From this set, was computed as the maximum over N-tuples (i.e., a set of N successive, z-sorted atoms) of the average axial distance between the atoms and their mean positionwhere ‖..‖ denotes the distance between points projected onto the z axis and is the moving average (of width N).

This CV measures the largest depletion in the vertical column without requiring it to be centered at any reference (y0 or bilayer midplane) and will therefore detect a depletion defect that is centered above or below the bilayer midplane. The exact nature of the depletion is not further specified; this CV can measure both a gap between merging atoms as well as a connected, but thinned, column. The fact that this CV does not saturate when the gap is closed is a potential benefit compared to similar formulations that do (14,23). A graphical example of computing this CV is presented in Fig. S7.

Committor regression models

Next, we tested the above geometric descriptors, alone and in combination, for their similarity to the committor RC. To this aim, a committor regression scheme was used to find optimal RCs and to build low-dimensional TS models for nucleation.

Given a configuration x, a committor estimate ĉ(x|N) was computed from N = 16 unrestrained (i.e., without gizmo potential), velocity-randomized trajectories of 6 ns length, each spawned from x. By counting how many trajectories reached absorbing boundaries in the A and B basins first (N and N, respectively), we obtained ĉ(x) = N/(N + N); trajectories that reached neither boundary (13–16% of cases) were discarded and not rerun. The absorbing boundaries for A and B were placed at  = 1.2 and  = 0.35, respectively. We denote a data set of configurations and their committor estimates as  = {(x, ĉ(x|N))}. (The Supporting Materials and Methods describe full details of the committor ensembles).

To predict ĉ(x) given CVs (ϕ1(x), ϕ2(x), …), we fit a logistic regression model f (cf. Fig. 2 a),where a(x) = w0 + wϕ(x) is a linear combination of CVs (including a constant offset w0) and ε is the residual error associated with x. The logistic function was chosen because ĉ(x) is expected to have a sharp transition between its limiting values 0 and 1. The weights for each model (w0, w1, …) were fitted by gradient descent to minimize the mean-squared error, MSE = , for a training data set (equivalent to maximizing the coefficient of determination R2 = 1 − . The R2-value for a set of (omitted) cross-validation data was used as the model’s score.

Each regression model f was used with to build a low-dimensional TS model (cf. Fig. 2 b) in the corresponding CV space. This minimal TS description comprises the CV basis (ϕ1, ϕ2, …), the regression model f, and an estimate of the minimal free-energy path (MEP) that connects the A and B basins. Whereas the f isosurfaces are directly accessible by solving f = f′, f does not specify where the separatrix is most likely crossed (saddle point), nor does it specify the orientation of the MEP (28). In fact, the MEP may not be orthogonal to the committor isosurfaces in the chosen coordinate basis. Therefore, a partial MEP segment was estimated with a three-point “string” ϕMEP = [ϕ0.2, ϕ0.5, ϕ0.8], where ϕ0.2, ϕ0.5, and ϕ0.8 are averages of points in binned by ĉ(x) into the intervals (0, 0.4), (0.4, 0.6), and (0.6, 1.0), respectively. With this definition, ϕ0.5 estimates the free-energy saddle point. Because the points in are drawn from , the MEP segment is only an estimate of the true MEP. We tested this estimate by comparing the committor-derived models with the computed PMFs for the unbiased system.

To test to what extent the H and I-gizmos bias sampling to the same TS regions, we built two separate committor data sets and using H- and I-gizmo-biased simulations. Further details on the committor data sets are given in the Supporting Materials and Methods.

Methods

MD simulations of the membrane-solvent-gizmo system were carried out essentially as if the gizmo were a transmembrane protein, using Gromacs 5.05 (40). The simulated system was a patch of 128 Berger DMPC lipids (20) that were solvated with single point charge (SPC) water (41) to a 43:1 water/lipid ratio using the MemGen web server (42). The gizmo was placed approximately at the membrane center. The simulation box used periodic boundary conditions and had dimensions 6.2 × 6.2 × 7.9 nm.

In all simulations, a temperature of 323 K was set with independent, velocity-rescale thermostats (43) for solvent, lipids, and the gizmo, all using a coupling of 2.0 ps. A pressure of 1.0 bar, in NPT simulations, was maintained via semi-isotropic weak coupling (44) using a time constant of 1.0 ps. Center-of-mass motion was removed from the lipids, gizmo, and solvent groups independently. The center atom of the gizmo was position restrained in the horizontal (x and y) dimensions with harmonic spring constants of 1000 kJ/(mol nm2). Electrostatics were computed with particle mesh Ewald (45) using a real space cutoff of 1.2 nm. A 1.2 nm cutoff was also used for van der Waals interactions.

Before NPT production runs, the system was prepared using a steepest descent minimization followed by brief NVT and NPT runs (100 ps apiece) using a 2 fs timestep. Production runs used a 4 fs timestep. All NVT and NPT simulations used a stochastic integrator (46). Water bonds and angles were constrained by the SETTLE algorithm (47), and all other bonds were constrained using LINCS (48).

The Lennard-Jones interactions between the gizmo chain (CX beads) and lipid tails were scaled by using an alchemical mutation to turn CX beads into (noninteracting) GG beads, with λCX implemented via the “vdw-lambdas” free-energy option in Gromacs. The restraint potential for the gizmo belt radial breathing mode V (see Eq. 1) was implemented using the essential dynamics options in Gromacs. Data S1 includes example files for including a gizmo in a Gromacs simulation.

Results

Committor-guided search for optimal RCs

To test how accurately a single RC describes TS crossing, we first used the above sketched scheme to train and score one-dimensional (1D) models (setting ϕ2 = 0) before building two-dimensional (2D) models. These 1D models (derived from ) yielded scores as high as R2 = 0.67. The atom group was found to be highly predictive of R2 (Fig. 3), with only NP- and P-based CVs having R2 > 0.38. The best CT CVs had weaker predictive power (R2 ≈ 0.35), and the CVs involving water (NPO and OW) had R2 ≤ 0.3. In Fig. 3, three examples of regression models are shown for the top-scoring NP-, CT-, and OW-based CVs.

Figure 3

Committor regression models derived using , the H-gizmo committor data set. (a–c) Top-scoring 1D models for CVs using the OW (a), CT (b), and NP (c) atom groups are shown. The ϕ terms are abbreviations for these CVs. Training and test data sets are drawn as purple and green dots, respectively, and a vertical Gaussian jitter was added to the (discrete) ĉ(x) values. The solid black lines show f and have closed and open caps indicating the directions of the closed and open pore states, respectively. The dashed lines are ± two standard deviations from f, assuming 16 (N) samples drawn from a binomial distribution. (d) All regression models, 1D and 2D, are shown co-sorted by R2. The labeled orange points are the best 1D models for each atom group listed in Table 1. (e) The regression fit for the 2D model using CVs ϕCT and ϕNP, with the same formatting as (a)–(c). (f) 2D view of the model in (e) is given, showing the training data set and the resulting TS model . Individual points (small markers) are colored by ĉ(x), with cyan, red, and blue dots belonging to the intervals (0, 0.4), (0.4, 0.6), and (0.6, 1.0), respectively. Gray squares and diamonds indicate ĉ(x) = 0 or 1, respectively. Large markers are averages of the respective subsets, and the ellipses are one-standard-deviation contours. To see this figure in color, go online.

This sorting suggests that TS passage is best characterized by localized headgroup (NP) merger. Lipid tail (CT) expulsion is also involved, but to a much lesser extent, and the presence or absence of water (OW) in the nascent pore is, unexpectedly, not predictive of barrier crossing (at least by itself). The top-ranked CV, with R2 = 0.67, was , which measures the axial merger of the four N or P atoms penetrating the bilayer. Presumably, this captures a pair of lipid headgroups merging from opposite sides. This CV is henceforth abbreviated as ϕNP. We also abbreviate and , the top-ranked CT and OW CVs, as ϕCT and ϕOW, respectively.

According to this ranking, the NPO CVs, which group together the charged N and P atoms as well as water oxygens, are less predictive than NP and P alone. This result suggests that water and headgroup merger are not tightly coupled at the TS; it also suggests that the mechanism might proceed in a stepwise fashion.

Next, we considered 2D models, testing whether any linear combinations of two CVs were more predictive than ϕNP alone. As summarized in Fig. 3, a few 2D models slightly outperformed ϕNP. These improved models all contain ϕNP paired with another CV; however, the increase in R2 (at most R2 = 0.69 as opposed to R2 = 0.67) was modest.

The high ranking of NP and P regression models is consistent between and data sets, and the absolute R2-values differ slightly (see Fig. S5 for a side-by-side comparison). For the other atom groups with lower R2 scores (OW, CT, C, and NPO), there are differences in the relative rankings between the two data sets that we attribute to the gizmos facilitating similar, but not identical, pathways. These pathways were probed further by computing free-energy landscapes.

Free-energy landscapes

Next, we constructed free-energy surfaces by projecting the PMF ensembles onto the optimized CVs constructed in the previous section. First, to examine headgroup merger and tail depletion, we used CVs ϕNP and ϕCT. Later, to compare headgroup and water penetration, we use ϕNP and the analogous water penetration coordinate. (PMF projection and unbiasing required a modified version of the weighted histogram analysis method (WHAM) (39), described fully in the Supporting Materials and Methods.) We also tested how computed PMFs vary with gizmo design (H versus I) and tested for convergence using ensembles that were started from closed (“A”) and open (“B”) states.

Fig. 4 shows the free-energy landscapes as a function of the two CVs ϕNP and ϕCT. In both the 2D and marginal PMFs, two clear minima can be seen, corresponding to the unperturbed membrane and the metastable pore at a relative energy of +10 kBT. The barrier heights are between 14 and 17 kBT, which is consistent with the range of previously reported values for the same system (17,21,22). The PMFs also resolve the full metastable pore that has been experimentally predicted for decades (11,12).

Figure 4

2D (a and b) and 1D (c and d) PMFs of pore nucleation. (a and b) PMFs computed using the I-gizmo (a) and H-gizmo (b), with ensembles seeded from the closed (A) state, are shown. Contours have 1 kBT spacing. The insets show snapshots of prebarrier structures that have identical positions in (ϕCT, ϕNNP) space, as indicated by the orange and yellow stars. In the snapshots, water is drawn as a stick with oxygen and hydrogen atoms in red and white, respectively. Lipid headgroup N and P atoms are drawn as blue and orange spheres, respectively, and the lipid tail carbons are illustrated as a gray slab that has been cut away. (c and d) Marginal 1D PMFs are shown. For each RC, 1D PMFs are shown for four cases: both gizmo types (I and H) and both starting structures (A and B). In all panels, the overlays are as described in Fig. 2. Severely undersampled bins (bin count ≤ 20), which make the PMF perimeters jagged, were not plotted; however, this has no effect on the energetics or barrier heights shown. To see this figure in color, go online.

In all cases considered in Fig. 4, the quantitative agreement between A/B pairs suggests that the PMFs are sufficiently converged and do not suffer from marked hysteresis effects. An autocorrelation analysis (detailed in the Supporting Materials and Methods) showed that intermediate windows in each ensemble exhibited slow, two-state switching corresponding to TS crossing, with autocorrelation times of 50–100 ns. Because this switching was clearly the sampling bottleneck, we extended these windows to 2 μs to ensure convergence. The resulting PMFs all show a net free-energy change ΔG of 9–11 kBT irrespective of the particular choice of gizmo, CV, and starting structure. Only the barrier heights vary, which is to be expected.

The PMFs in Fig. 4 exhibit some clear differences that shed light on the nucleation mechanism. The location of energy barriers, and to some extent their heights, depends on the CV choice and the gizmo design. The I-gizmo simulations give barriers that are lower and later than the H-gizmo simulations. This result suggests that the two gizmos are not both directing sampling along exactly the same reaction path; otherwise, the PMFs projected onto the same RC would be identical. The snapshots shown in Fig. 4 (Fig. 4, a and b, insets), reveal the difference in paths. In Fig. 4 a, the membrane has partially thinned (I-gizmo), but in Fig. 4 b, a hydrophobic column already connects the opposite sides (H-gizmo). Together, this suggests that hydrophobic defects are energetically unfavorable and unlikely to initiate nucleation, whereas an indentation pathway is energetically preferred by 3–4 kBT.

Fig. 4 also shows that the computed energy barriers somewhat depend on the chosen CV, which underscores the well-known fact that proper choice of RC is crucial. The marginal PMFs, shown in Fig. 4, c and d, show that ΔG‡(ϕNP) is always higher than ΔG‡(ϕCT), suggesting that ϕNP is better at distinguishing the A and B basins and that state densities partially overlap when projected onto ϕCT. This projection overlap for ϕCT is evident in the 2D PMFs (Fig. 4, a and b) as well. Given that ϕCT primarily captures in-plane, radial depletion of lipids, this corroborates that pore radius is a suboptimal RC for pore nucleation but is better suited to track pore expansion, such as would occur with applied surface tension.

Based on these findings, we take the PMF G(ϕNP), for the I-gizmo, to be the most accurate estimate of the true PMF (for a single RC in our pool). Here, ϕNP resolves the TS at values between 0.6 and 0.75, which is well before the headgroups fully reach the pore center and before the RC saturates, at ϕNP ≈ 0.2. The PMF G(ϕNP) indicates a barrier for pore closure of around 5 kBT. This value has not been reported previously because the full PMF, including the entire metastable basin, was not resolved with a single RC. Just recently, however, this result has been confirmed, independently (49), using an empirical RC that combines the slab scheme of (14) with a lateral expansion RC. A previous study also reported that for 20 independent, unbiased simulations of DMPC pores, zero closed within 500 ns (17), which suggests that the closure barrier is substantial. The energy landscape suggests that applied surface tension should increase the radius of this metastable defect, in accordance with a Litster type energy model, as has been reported in other studies (13,49).

Our previous committor analysis provides an independent estimate of the TS and an additional control to determine which PMF results are the most accurate. Accordingly, we have overlaid the corresponding TS models, onto the PMFs in Fig. 4. For the I-gizmo ensemble (Fig. 4 a and marginals Fig. 4, c and d), the TS model features align well with the projected PMF, and the f(x) = 0.5 isosurface (black dashed line) cuts the PMF nearly at its saddle point. This effect is also seen in the marginals (Fig. 4, c and d), where ϕ0.5 (red dots) are near the PMF maxima (blue and green curves). In contrast, the PMF barriers are higher and much earlier for the H-gizmo ensembles and do not align well with ϕ0.5. This result suggests that forming a hydrophobic column is energetically too costly and thus provides additional evidence for the indentation path.

The role of water

Despite our finding that water CVs correlate poorly with the committor, one would expect water to partially hydrate the lipid headgroups that submerge to create the pore. Thus, to examine this idea and to better resolve the sequence of events and the role of water, we computed a 2D PMF using a water CV paired with ϕNP. The CV was chosen because it measures the axial merger of the first (four) water molecules, analogous to how ϕNP measures the merger of headgroups.

Fig. 5 shows that as was the case for ϕCT, does not cleanly distinguish the A and B basins and gives rise to projection overlaps, thus rendering it a poor RC by itself. The marginal PMFs in Fig. 5 d do not reveal this overlap directly, but the PMF maximum and ϕ0.5 (red dot in the overlay) are in poor agreement. Indeed, in Fig. 5, a and b, the MEPs (cyan, red, and blue dots) are nearly orthogonal to , suggesting that the water penetration is largely orthogonal to TS crossing. In addition, is a poor committor correlate (R2 = 0.32 and R2 = 0.07 for and data sets, respectively).

Figure 5

2D (a and b) and 1D (c and d) PMFs of pore nucleation for headgroup and water CVs, ϕNP and . The formatting is as in Fig. 4. The overlays in (c) and (d) are from the H-gizmo committor ensemble. To see this figure in color, go online.

In contrast to the linear MEPs seen in Fig. 4, both of the 2D PMFs in Fig. 5, a and b exhibit a kinked pathway in the TS region. The kink results from saturating at its lowest values before the TS. The MEP overlay also shows this kink at the cyan dot, ϕ0.2.

Taken together, these findings suggest a “water first” mechanism, in which a small number of waters reach the pore center before the headgroups and before ϕ0.5. The subsequent barrier crossing, ϕ0.2 → ϕ0.5 → ϕ0.8, follows a straight line in the projected space, primarily along the ϕNP coordinate.

To further examine the sequential steps of nucleation and the roles of hydration and headgroups, we used the trained TS model for the I-gizmo, , to select representative structures at different stages of nucleation. Fig. 6 shows those frames (points x from ) from just before, at, and just after the TS (rows Fig. 6, c, e, and g, respectively) that agree best with f (residuals |ε| ≤ 0.1). The prebarrier frames (row Fig. 6 c) all yield ĉ(x) < 0.2 and illustrate the diversity of water defects that pierce the membrane before the TS. This finding is consistent with the PMFs and kinked MEP described above, indicating that hydration defects precede the TS. Next, as the barrier is crossed, the merger of lipid headgroups (ϕNP decreasing) is also clearly visible by comparing rows Fig. 6, c, e, and g. Within these rows, there are variations in hydration, shearing or skewing, pyramidal indentations, and axial asymmetry, suggesting that this variability is orthogonal to the MEP and thus does not correlate with the progress of the reaction. The snapshots show that the submerged headgroups are hydrated, as expected; however, the size and shape of the water clusters appear too variable for water to serve as a precise measure of reaction progress (3).

Figure 6

Configurations showing TS crossing, taken from the data set. The regression model (a) and the corresponding TS model (b) for the CVs ϕNP and are shown. In (a), f is shown as a solid line, and the dotted lines are ± two standard deviations, assuming 16 (N) samples drawn from a binomial distribution. In (b) the solid and dashed lines correspond to those in Fig. 5a. In (a) and (b), the committor data set is shown as gray dots. The colored dots indicate configurations that are illustrated in rows (c)–(g). These illustrated configurations were chosen based on ĉ(x)-values and residual errors. Rows (c), (e), and (g) show configurations from before, at, and after the TS, respectively, with small residual errors. Configurations in rows (d) and (f) have large negative and positive error residuals, respectively. To see this figure in color, go online.

For comparison, we also plotted snapshots from cases in which ĉ(x) is severely under- or overpredicted in rows Fig. 6, d and f, respectively. The structures in rows Fig. 6, d–f are visually similar, and no striking structural differences stand out to suggest other predictive CVs.

A considerable amount of the ĉ(x) variance here stems from using imprecise (N = 16) committor estimates. Whereas upwards of 100 shots are needed to converge a single ĉ(x) estimate, this is largely due to statistical uncertainty. Indeed, binomial deconvolution methods can be used to estimate committor histograms rigorously, at 10× reduced cost (50). Here, as Figs. 6 a and 3, c and e show, the observed variance in ĉ(x|N = 16) is indeed comparable to the statistical uncertainty (dotted lines above and below the regression curve f). This suggests that, for the purpose of ranking CVs, only low-resolution committor estimates are required, provided a sufficient number of uncorrelated configurations near the TS. (For more information about the comittor data sets, see the Supporting Materials and Methods).

Discussion

We have used a combination of MD simulations, free-energy calculations, and committor analysis to probe the mechanism of metastable pore formation in a lipid membrane.

Several previous simulation studies of poration have proposed and employed one or a few possible RCs, considering specific atom groups and CVs in isolation. These approaches include restraining a single lipid headgroup relative to the bilayer midplane (17,22,23,51,52), biasing hydrophilic atoms to occupy a stack of cylindrical slabs spanning the membrane (14), biasing water to occupy a cylindrical column (53), and growing a lateral depletion of lipid centers of mass (54).

Here, we adopted a more general approach and constructed a combinatoric pool of CVs to systematically vary the different factors relevant to the mechanism, including the atom group, geometric bias, and locality (via N, the number of atoms in the CV definition). These CVs were then assessed systematically as RCs, using a committor regression scheme. For comparison, hybrid RCs built from pairwise linear combinations of individual CVs were also tested but not found to be substantially superior.

Our results showed that the achievable RC quality largely depends on the choice of atom group. Headgroup (NP and P) based CVs yielded the highest correlation with the committor and, in this sense, represented the best descriptors for the RC in the vicinity of the TS. In contrast, and somewhat unexpectedly, CVs based only on the position of water molecules (the OW family) provided remarkably poor RCs. Also, the NPO CVs (grouping N, P, and (water) O atoms) described the RC less accurately than NP alone, suggesting that these hydrophilic components have distinct mechanistic roles in the TS crossing. Lipid-tail-based CVs (CT) were also less predictive than headgroups.

Among the headgroup RCs (and among all RCs tested), (ϕNP), which essentially measures the gap between two penetrating lipid headgroups, provided the best description of the commitor RC. Interestingly, none of the combinations of ϕNP and CVs involving larger numbers (10–20) of atoms, which measure larger collective motion, substantially outperformed ϕNP alone. This result suggests that barrier crossing is more localized than one might have expected.

The particular RC ϕNP somewhat resembles the hydrophobic belt height suggested by Akimov et al. (19). The fact that our systematic and unbiased search highlighted this RC corroborates its mechanistic relevance and utility in continuum modeling. Our RC ϕNP is also similar to the slab occupancy metric (14) referred to above. Both RCs essentially measure the gap between the penetrating hydrophilic material but for different groups of atoms, which turned out to be an important distinction. This RC, as well as the other RCs considered here, also provides a rather smoother measure of merger progress, which does not require spatial discretization and thereby avoids Poisson noise due to near-zero bin counts.

Overall, our committor analysis and RC ranking underscore the importance of headgroup merger, both as a mechanistic descriptor and biasing coordinate. However, regarding the precise molecular mechanism, it is crucial to systematically define and assess RCs against the, by definition, optimal committor RC. Unexpectedly, water turned out to be rather weakly correlated with the progress of pore formation; further, the degree of hydration fluctuates markedly during poration, which, taken together, renders hydration a rather poor RC. Water-based CVs (by themselves), therefore, seem to provide less suitable biasing coordinates for poration and possibly similar systems and processes such as membrane fusion. Lateral depletion coordinates, in contrast, such as ϕCT, are better for describing pore expansion than nucleation.

Future studies could use a similar methodology to probe even more localized details of nucleation, such as headgroup orientation and hydrogen bonding networks, in search of still better committor correlates.

To obtain sufficient sampling of the hard-to-reach TS region, we developed a biasing scheme, the pore gizmo. The intended effect of this gizmo bias potential, exerted via repulsive beads, is similar in spirit to other schemes using repulsive plates (55), hydrophilic beads (56), depletion coordinates (54), or constrictive “cuff potentials” (57). What distinguishes the gizmo and, as it turned out, renders it particularly suitable for this purpose, is its flexibility, which provides the combined system with sufficient flexibility to enhance sampling and thereby to also probe asymmetric states. In particular, the axial flexibility of the CX chain allows for proper sampling of asymmetric or skewed tail-packing defects.

Using our gizmo biasing potential enabled us to obtain sufficient sampling to study, by atomistic simulations, the full poration pathway and the mechanistic steps of a DMPC lipid bilayer before nucleation. Two different biasing schemes were used to direct a membrane thinning (indentation) path and a hydrophobic defect path, both leading to open pores. By computing PMFs, projected onto the optimal RCs determined before, and by comparing snapshots close to the TS, a sequential penetration of water and headgroups was resolved.

Taken together, our simulations and committor-based RCs suggest the following mechanism. A pore nucleates via an elastic indentation that is energetically preferred to forming a hydrophobic defect. After thinning but just before the TS, water pierces the thinned slab. This water defect is required but does not suffice to finally nucleate the pore. Instead, it is the axial merger of the first lipid headgroups from opposite monolayers (ϕNP) that precedes—and best characterizes—the subsequent nucleation and thus crossing of the TS.

This mechanism could be biologically relevant for the structure and rupture of the hemifusion diaphragm (HD) (58,59). In this context, the HD perimeter is a three-way membrane junction that is predicted to stabilize and potentially lower the nucleation barrier for transient pore defects (10). These rim pores are likely on-path precursors to rupture and therefore might also explain flickering before HD opening (60). A similar combination of enhanced sampling and committor analysis could enable one to study this and other biologically relevant topological membrane remodeling mechanisms, such as fusion stalk formation.

This methodology could also be used to systematically study how nucleation mechanisms vary with respect to membrane properties such as thickness and rigidity by testing different compositions and physical conditions.

Similar to lipids, intrinsically disordered proteins (IDPs) can regulate transport but are challenging to control and characterize structurally. For example, there is an ongoing discussion on how the IDPs that form the interior of nuclear pore complexes (Phe-Gly nucleoporins or FG-Nups) are organized and spatially distributed within the nuclear pore to achieve selectivity for karyopherins (61). Gizmo-based enhanced sampling may enable one to energetically distinguish between these hypotheses.

In a broader context, future studies that seek to uncover atomistic mechanisms and to compute energetics with kT precision, particularly for permutationally frustrated systems such as membranes, solvent surface layers, or IDPs, will benefit from hybrid schemes that combine committor-analysis-based RC optimization and free-energy calculations using gizmo-biased enhanced sampling.

Author Contributions

G.B. and H.G. designed research. G.B. performed research and analyzed data. G.B. and H.G. wrote the manuscript.