Molecular Mechanism for Gramicidin Dimerization and Dissociation in Bilayers of Different Thickness.

Membrane protein functions can be altered by subtle changes in the host lipid bilayer physical properties. Gramicidin channels have emerged as a powerful system for elucidating the underlying mechanisms of membrane protein function regulation through changes in bilayer properties, which are reflected in the thermodynamic equilibrium distribution between nonconducting gramicidin monomers and conducting bilayer-spanning dimers. To improve our understanding of how subtle changes in bilayer thickness alter the gramicidin monomer and dimer distributions, we performed extensive atomistic molecular dynamics simulations and fluorescence-quenching experiments on gramicidin A (gA). The free-energy calculations predicted a nonlinear coupling between the bilayer thickness and channel formation. The energetic barrier inhibiting gA channel formation was sharply increased in the thickest bilayer (1,2-dierucoyl-sn-glycero-3-phosphocholine). This prediction was corroborated by experimental results on gramicidin channel activity in bilayers of different thickness. To further explore the mechanism of channel formation, we performed extensive unbiased molecular dynamics simulations, which allowed us to observe spontaneous gA dimer formation in lipid bilayers. The simulations revealed structural rearrangements in the gA subunits and changes in lipid packing, as well as water reorganization, that occur during the dimerization process. Together, the simulations and experiments provide new, to our knowledge, insights into the process and mechanism of gramicidin channel formation, as a prototypical example of the bilayer regulation of membrane protein function.

Significance

Membrane protein function is affected by the lipid bilayer environment. Gramicidin channels form by a transmembrane monomer ↔ dimer association, which makes them suitable to study the energetic coupling between bilayer-spanning inclusions and their host bilayer. We used extensive atomistic molecular dynamics simulations to characterize the bilayer dependence of gramicidin channel formation/dissociation to explore the channel-bilayer energetic coupling, which allowed us to capture—for the first time, to our knowledge—the monomer → dimer transition in lipid bilayers at molecular-level resolution. These simulations, supported by fluorescence experiments, show that small changes in bilayer thickness can have large effects on the channel function. The simulations also revealed the importance of lipid packing, protein structure deformation, and water organization in channel association and disassociation.

Introduction

The interactions between membrane proteins and the surrounding lipid bilayer remain a fundamental problem in biology and biophysics. The interactions can span from specific lipids directly binding to membrane protein binding sites to nonspecific interactions mediated through the bulk properties of the lipid membrane such as thickness and curvature. Sorting out the relative importance of the bulk bilayer properties compared to specific interactions is difficult because of the multitude of different lipids, the soft, fluid lipid environment, multiple overlapping timescales, and cooperative effects. Recent work combining experiments and molecular dynamics (MD) simulations are shedding light on the nonspecific lipid-protein interactions. At the molecular level, how do membrane properties such as bilayer thickness influence membrane protein activity? In this study, we approach this problem using the gramicidin channels as a prototypical membrane protein system.

The bilayer-spanning gramicidin channel forms by transmembrane dimerization of two nonconducting monomers (1,2), and the antiparallel dimer structure is stabilized by six hydrogen bonds at the formyl-N-termini (3, 4, 5). Each hydrogen bond contributes ∼4 kT to the channel stability (6, 7, 8, 9), indicating that gramicidin channel dissociation occurs on a timescale far beyond microseconds. The measured gramicidin channel lifetimes are sensitive to the lipid bilayer thickness: dependent on the lipid composition, the channel lifetimes may vary by an order of magnitude in response to subnanometer changes of lipid bilayer thickness (10,11). This variation of gramicidin channel lifetime as a function of lipid bilayer thickness has been explored by continuum elastic models and molecular simulations on gramicidin A (gA) (12,13). The hydrophobic length of the gA channel, estimated to be ∼2.2 nm based on both single-channel (11) and x-ray diffraction (14) experiments, is shorter than the lipid bilayer hydrophobic thickness. The prevailing model proposes that the energetic penalty associated with a thickness mismatch between the gA channel and the lipid bilayer causes the lipid molecules adjacent to the gA channel to deform to minimize the exposure of hydrophobic lipid acyl chains to water. This local bilayer deformation and the associated deformation energy impose a disjoining force on the gA channel that increases with increasing hydrophobic mismatch (15,16).

The gA channel has been the subject of MD simulation studies for more than five decades, beginning with the work by Mackay et al. (17). Early MD simulations focused mainly on structural refinement and ion permeation (3,18, 19, 20, 21), which are accessible with picosecond to nanosecond atomistic MD simulations. The scale of atomistic MD simulations has advanced both temporally and spatially over the last decade. It is now possible to employ atomistic MD simulations to investigate both thermodynamic properties and kinetic pathways of complex protein associations and dissociations (22,23). In their recent work, Im and colleagues performed a single 3.5-μs-long atomistic MD simulation on gA channel in mixed lipid bilayers to investigate the lipid redistribution surrounding the channel (24). With respect to the effect of lipid bilayer thickness on the gA channel stability, 100-ns-long atomistic MD simulations have been reported (25,26). These MD simulation studies, however, did not provide quantitative information about how lipid bilayer thickness influences the transmembrane gA channel formation or dissociation.

To systematically explore how lipid bilayer thickness influences gA dimerization and dissociation, we performed over-300-μs unbiased and biased MD simulations on gA channels in three bilayers with different acyl chain lengths: 1,2-dioleoyl-sn-glycero-3-phosphocholine (DC18:1PC), 1,2-dieicosenoyl-sn-glycero-3-phosphocholine (DC20:1PC), and 1,2-dierucoyl-sn-glycero-3-phosphocholine (DC22:1PC). To quantify the channel stability in the three different bilayers, we performed umbrella sampling simulations, totaling ∼120 μs, to determine the potential of mean force (PMF) for gA channel dissociation. It was found that a 0.4 nm increase in bilayer thickness can decrease the free-energy change for the gA dimer → monomer transition by 3–5 kcal/mol, and the results were supported by the fluorescence-quenching experimental results performed in this work. To elucidate the process of gA association (or dimerization) in the three different bilayers, we performed a total of 200 μs unbiased MD simulations with the distributed computing approach (27,28). These unbiased MD simulations showed that the kinetic rate of gA dimerization sharply decreased in the thick DC22:1PC bilayer, consistent with the PMF results. The unbiased MD simulations further revealed that gA dimerization process involves the cooperative induction of membrane deformation, unfolding of gA N-terminus, and reordering of lipid tails and water molecules at the subunit interface. The results from our atomistic MD simulations agree with the past and current experimental results and provide new, to our knowledge, mechanistic insights into the molecular mechanisms underlying gramicidin dimerization and dissociation in lipid bilayers.

Methods

Simulation systems and force field

The starting configurations for the MD simulations, with the [Val1]gA channel (Protein Data Bank (PDB): 1NRM) (29) embedded in the DC18:1PC, DC20:1PC, and DC22:1PC bilayers were generated using CHARMM-GUI (30). The CHARMM36 force field (31) with the cross-term energy correction map was used to model the L- and D-amino acids. The CHARMM-compatible force field parameters for the formyl and ethanolamine groups were from previous work (3,32). Lipids were modeled using the CHARMM36 lipid force field (33). Each simulation system contained 200 lipids (100 lipids in each leaflet) and 1 gA channel. The different systems were fully hydrated, containing 12,000–14,000 TIP3P water molecules and 0.15 M concentration of KCl.

MD simulations of bilayer-incorporated gA channels

We initially ran 2-μs-long MD simulations of gA in each bilayer; this was sufficiently long to equilibrate the gA channel structure in the lipid bilayers, as inferred from previous work by Allen et al. (3). All simulations were done in the semi-isotropic ensemble at 310.15 K and 1 bar using the GPU version of GROMACS (version 2018.3) (34). The Nosé-Hoover thermostat (35,36) and the Parrinello-Rahman barostat (37) were used to maintain the temperature and pressure. The LINCS algorithm (38) was used to constrain water geometry and covalent bonds involving a hydrogen atom. Lennard-Jones interactions were switched off smoothly at 1–1.2 nm, and the particle mesh Ewald method (39) was used to treat long-range electrostatics with a real space cutoff distance of 1.2 nm. Long-range dispersion corrections to the energy and pressure were not applied. Snapshots of each simulation were saved every 20 ps.

Umbrella sampling simulations of gA channel dissociation

Two different collective variables (CVs) were used for the umbrella sampling simulations. The first CV was the center-of-mass (COM) distance between the two gA monomers, where the COM is defined using all atoms in each monomer. Starting from the configuration at the end of the 2 μs equilibration, we used steered MD simulations to generate 50 configurations (windows) to gradually separate the two gA monomers in 0.05 nm steps: first, in the direction perpendicular to the bilayer plane (Z direction) until the COM distance in the Z direction reached predetermined values (1.8, 2.0, and 2.5 nm for DC18:1PC, DC20:1PC, and DC22:1PC bilayers, respectively. These distances correspond to the averaged monomer-monomer COM distance in the Z direction (see Fig. S1 a)), then in the direction parallel to the bilayer surface. The weighted histogram analysis method algorithm (40), as implemented in GROMACS, was used to derive the PMF profile.

The second CV was the intermolecular distance root mean-square deviation (I-DRMSD) of the two gA monomers’ backbone Cα atoms with respect to the reference gA channel structure (PDB: 1NRM). This CV was used in previous work (41) to study the dimerization of membrane proteins, and it is implemented in the PLUMED software package (version 2.4.2) (42). The I-DRMSD between the reference and the instantaneous configurations, X and X, is defined aswhere N is the number of Cα atoms in one gA monomer and D(x, x) is the distance between the coordinates x and x of two Cα atoms i and j belonging to two different gA monomers. A total number of 46 umbrella sampling windows, with I-DRMSD varying between 0 and 2.3 nm, were generated. The umbrella spring constant was 5500 kJ/mol/nm2 and the WHAM program from Alan Grossfield (http://membrane.urmc.rochester.edu) (43) was used to derive the PMF profile.

For each umbrella sampling window, we performed 400-ns-long MD simulations for a total umbrella sampling simulation time of 115.2 μs for the different bilayers and CVs. The last 300 ns of each window were used to construct the PMF profiles.

MD simulations of gA channel formation

At the end of umbrella sampling simulations, three configurations (one in each bilayer) with COM distances between the gA monomers close to 3 nm were selected for unbiased MD simulations of gA dimerization using the distributed computing approach (44,45). For two gA monomers embedded in the DC18:1PC and DC20:1PC bilayers, 500 independent 100-ns-long MD simulations were performed for each system. In the thicker DC22:1PC bilayer, a total of 1000 independent 100-ns-long MD simulations were performed for gA channel formation. The total unbiased MD simulation time was 200 μs.

Analysis and visualization of MD simulations

The MD simulation trajectory analysis was done using GROMACS, PLUMED, and MDAnalysis (see Supporting Materials and Methods for additional simulation analysis tools) (34,42,46). VMD was used for visualization (47).

Experimental materials

DC22:1PC in chloroform (25 mg/mL), DC20:1PC in chloroform (10 mg/mL), and DC18:1PC in chloroform (25 mg/mL) were from Avanti Polar Lipids (Alabaster, AL). [Val1]gA was a generous gift from Roger E. Koeppe II (University of Arkansas). The naturally occurring mixture of the linear gramicidins from Bacillus brevis was from Sigma-Aldrich Co. (St. Louis, MO). Historically, this mixture has been called gramicidin D (gD) after R. Dubos discovered the gramicidins (48); it contains 80–85% [Val1]gramicidin A (gA), 6–7% gramicidin B (gB, [Val1, Phe11]gA), and 5–14% gramicidin C (gC, [Val1, Tyr11]gA (49). gA, gB, and gC form structurally equivalent antiparallel, dimeric channels with very similar properties (50), meaning that ∼2/3 of the measured ion flux will be through symmetric gA/gA homodimeric channels, ∼1/5 will be through asymmetric gA/gB or gA/gC heterodimeric channels, and the remaining will be through symmetric gB/gB and gC/gC homodimeric channels and asymmetric gB/gC heterodimeric channels). Thallium nitrate (TlNO3), sodium nitrate (NaNO3), and HEPES were also from Sigma-Aldrich. 8-Aminonaphthalene-1,3,6-trisulfonic acid (ANTS) was from Invitrogen (Eugene, OR). All materials were used without further purification. Stock solutions of buffers and quenchers were prepared using Millipore Milli-Q deionized water Millipore Sigma (Burlington, MA); the pH was adjusted to pH 7 using NaOH and HNO3. The standard buffer was 140 mM NaOH plus 10 mM HEPES; the quench buffer was 50 mM TlNO3 plus 94 mM NaNO3 and 10 mM HEPES. The ANTS buffer was made with 25 mM ANTS, 100 mM NaNO3, and 10 mM HEPES and stored in the dark. The gD stock solutions were 500 μg/mL of the gramicidin mixture produced by B. brevis, which was dissolved in methanol and stored at a temperature of −40°C.

Gramicidin fluorescence assay

We quantified the gramicidin channel activity using a fluorescence quench assay (51,52), which is based on the quenching of an intravesicular fluorophore (ANTS) by the gramicidin-channel-permeant monovalent cation thallium (Tl+). ANTS-loaded large unilamellar vesicles (LUVs) were prepared with varying gD:lipid ratios by mixing the phospholipid with gD at the given molar ratio. The lipid-gramicidin mixture then was dried under nitrogen and further dried overnight in vacuum. The lipid was rehydrated with ANTS buffer in the dark; the volume of the ANTS buffer was adjusted to give a lipid concentration of 10 mM. The sample was equilibrated at room temperature for 3 h and sonicated for 1 min. After six freeze-thaw cycles, the gramicidin-lipid suspension was extruded 21 times with an Avanti miniextruder and a 0.1 μm polycarbonate filter to form LUVs. Extravesicular ANTS was removed with a 2.5 mL PD-10 desalting column (GE Healthcare, Piscataway, NJ), and the solution was covered and stored in the dark at 12.5°C. The vesicle size distribution was determined using dynamic light scattering with a Litesizer 500 (Anton Paar, Graz, Austria); the average diameter of the monodisperse vesicles was 130 ± 5 nm, with little variation among the three vesicle types. Assuming a molecular area/lipid molecules of 0.7 nm2 for all three lipids (53), each vesicle contains ∼150,000 lipid molecules.

The time course of the Tl+-induced quenching of the ANTS fluorescence was measured at 25°C using a SX-20 stopped-flow spectrofluorometer (Applied Photophysics, Leatherhead, UK) with a sampling rate of 5000 points/s and an instrument dead time of <2 ms. The excitation wavelength was 352 nm, and the emitted light was recorded above 450 nm using a high-pass filter (Applied Photophysics). Each LUV sample was incubated at 25°C in the dark for at least 10 min before the measurement. The gramicidin channel activity was quantified by first determining the fluorescence in the absence of the quencher (Tl+), and then determining the time course of fluorescence quench when mixing the LUVs with the Tl+ buffer. The inevitable variation in LUV sizes means that the time course of fluorescence quenching cannot be described by a single exponential decay. Rather, the time course is the sum of exponential decays with time constants and weights that reflect the distribution of vesicle sizes and the number of conducting channels in the vesicle membrane, which can be described as a so-called stretched exponential (54)where F(t) denotes the fluorescence intensity of ANTS at time t, F(0) and F(∞) are the initial and final fluorescence values, β (0 < β ≤ 1) accounts for the dispersity of the vesicle population, and τ0 is the time constant. F(0), F(∞), β, and τ0 were determined from a nonlinear least-squares fit of Eq. 2 to the first 200 ms of the quench curve, and the quench rate was defined as (51)evaluated at t = 2 ms.

Results

To characterize the gA channel’s sensitivity to subtle changes in bilayer thickness and the gA dimerization process in lipid bilayers, we performed extensive umbrella sampling and unbiased MD simulations in bilayers of three different thicknesses: DC18:1PC, DC20:1PC, and DC22:1PC, respectively. The starting configurations of gA monomers and gA channels embedded in these bilayers are illustrated in Fig. 1 a.

Figure 1

PMF profiles for gA channel dissociation into two monomers in the three lipid bilayers. (a) Snapshots of gA monomers and dimers in the three different bilayers tested. Water molecules (data not shown) fill the monomers and the dimers at the start of the umbrella sampling simulations. In (b), the used collective variable is the COM distance between the two gA monomers, and in (c), the collective variable is the I-DRMSD of the backbone Cα atoms. The insets illustrate configurations for the monomer, transition, and dimer states of the two gA monomers, which are colored orange and cyan. The lipid phosphorous atoms are colored tan. Water and lipids are not shown for clarity.

Free energy of gA channel dissociation

To address the thermodynamics of gA dimer ↔ monomer transition, which is closely related to the distributions of gA dimers versus monomers in the lipid bilayers, we conducted atomistic umbrella sampling calculations to map out the free-energy landscape for gA channel dissociation using two different CVs: monomer-monomer COM distance and the I-DRMSD.

The PMF profiles obtained with the two CVs (Fig. 1, b and c) are in overall agreement. First, the PMF profiles for gA channel dissociation, with the dimeric state and the monomeric state separated by a transition state, are reminiscent of a two-state model for ligand unbinding from proteins (55). Second, the PMF profiles indicate that the gA dimer is energetically favored relative to the monomer in the three lipid bilayers. The energy differences between the monomers and dimers states in DC18:1PC, DC20:1PC, and DC22:1PC lipid bilayers are −14.4 ± 0.3, −11.5 ± 0.7, and −6.2 ± 0.6 kcal/mol with COM distance as the CV; these values are −15.6 ± 0.3, −12.8 ± 0.3, and −5.0 ± 0.3 kcal/mol using the I-DRMSD as the CV. Thus, it is estimated that the thermodynamic equilibrium distribution of dimers versus monomers is sharply shifted toward the monomer state when going from DC20:1PC to DC22:1PC bilayers.

The activation barrier for gA channel dissociation , which determines the lifetime of the channel, can be derived from the PMF profiles. Using the COM distance as the CV, -values for gA channel dissociation in DC18:1PC, DC20:1PC, and DC22:1PC lipid bilayers are 14.8 ± 0.3, 11.4 ± 0.6, and 8.6 ± 0.6 kcal/mol, respectively; the corresponding values using the I-DRMSD as the CV are 15.7 ± 0.3, 13.1 ± 0.3, and 8.1 ± 0.2 kcal/mol, respectively. The channel lifetimes (τ) vary with ΔG via ; our PMF results thus provide insights into the channel lifetimes decrease with increasing lipid bilayer thickness (11).

The PMF profiles also allow for estimating the activation barrier for gA channel formation , which is equal to the peak values of the PMF profiles in Fig. 1. Using the COM distance as the CV, the predicted values in DC18:1PC, DC20:1PC, and DC22:1PC bilayers are 0.6 ± 0.1, 1.1 ± 0.1, and 2.8 ± 0.2 kcal/mol. Using the I-DRMSD as the CV, the corresponding values are 1.0 ± 0.2, 1.1 ± 0.2, and 3.2 ± 0.2 kcal/mol, respectively. -values in the DC18:1PC and DC20:1PC bilayers are comparable to the thermal energy (0.6 kcal/mol). The increase in in the thicker DC22:1PC bilayer may be attributed to the larger bilayer deformation associated with channel formation, which will be further explored below.

The PMF profiles (Fig. 1 b) indicate a threshold monomer-monomer COM distance for the two subunits reaching the transitions state. This distance is an important parameter when relating the kinetics of gA channel dissociation to continuum elastic models (12,56). The PMF in DC22:1PC bilayer suggests that this distance is ∼1.6 Å, in good agreement with previous estimates (8,57). The PMFs in DC20:1PC and DC18:1PC also show a shoulder at ∼1.6 nm monomer separation, but the definition of the transition state is not as clear.

The energetic coupling between gramicidin channels and their host bilayer can be described using continuum models of elastic lipid bilayer deformations (12,58, 59, 60), where the bilayer contribution to the free energy of dimerization can be expressed aswhere and denote the bilayer deformation energies associated with the gramicidin monomer and the dimer, respectively; HB is a phenomenological spring coefficient; and dmin is the bilayer thickness where the has its minimal value, . Fitting the energy minima in Fig. 1 to Eq. 4 allows us to estimate (= −15 kcal/mol), HB (= 6.9 kcal/mol × nm2), and dmin (= 2.5 nm), where HB and dmin are in good agreement with the estimates (HB = 8.6 kcal/mol × nm2 and dmin ≈ 2.6 nm) obtained by Sodt et al. (61), who used the three-dimensional elastic model to analyze the gramicidin-channel-induced bilayer deformations deduced from MD simulations (24).

gA channel stability

The lifetime of the gA channel in the DC18:1PC bilayer is on the order of seconds (62), meaning that one would expect that the gA channel would remain stable during a few-microseconds-long atomistic MD simulations. Indeed, throughout the 2 μs simulations, the channels remained stable in all three systems because of the strong interactions between the two monomers (the interaction energies between the two monomers were ∼−45 kcal/mol (see Fig. S2), irrespective of bilayer thickness). The monomers were tightly associated, with no detectable monomer-monomer wobbling or rotation. To quantify this further, we calculated the time evolution of the COM distance between the two monomers in the dimer. Fig. 2 a shows that the average COM distances between the two gA monomers are ∼1.45 ± 0.02 nm in three bilayers.

Figure 2

(a) Time evolution of the COM distance between the monomers in the gA dimers (channels). (b) Time evolution of RMSD for the channel’s backbone Cα atoms in the three lipid bilayers. The light-colored traces depict results for all backbone Cα atoms; the darker traces exclude Trp-15 and the ethanolamine. The inset shows gA channel backbones structures with large and small RMSD fluctuations (colored cyan) with the NMR structure (PDB: 1NRM, colored orange) superimposed for reference. To see this figure in color, go online.

Fig. 2 b shows the root mean-squared deviation (RMSD) profiles for the channel’s backbone Cα atoms. The RMSD in the DC18:1PC lipid bilayer shows relatively large oscillation between 170 and 590 ns (light-colored line). The RMSD profiles also exhibit sporadic spikes in the DC20:1PC and DC22:1PC bilayers. To understand the molecular basis for these larger RMSD fluctuations, we show two gA channel configurations with large and small RMSD fluctuations in Fig. 2 b. The gA channel configurations with large RMSD fluctuations (>0.1 nm) have an outward splay of one of the ethanolamine-C-termini, whereas the gA channel configurations with small RMSD fluctuation (∼0.05 nm) overlaps with the PDB: 1NRM structure. When we exclude the C-terminal ethanolamine and Trp-15 (darker colored traces in Fig. 2 b), the RMSD fluctuations are reduced in all three systems. We conclude that there is little gA monomer-monomer wobbling and twisting at the microsecond timescale in the lipid bilayer environment (though such movements could be important for gA channel dissociation at much longer timescales).

Kinetics of gA channel formation

To elucidate the process of gA dimerization in lipid bilayers, we did a total of 2000 independent 100-ns-long unbiased MD simulations and calculated the I-DRMSDs for the two gA peptides’ backbone Cα atoms, with I-DRMSD < 0.08 nm defined as the dimer state and the I-DRMSD > 0.4 nm as the monomer state (Fig. 3). In DC18:1PC bilayers, 56 out of the 500 systems ended up forming a channel within the 100 ns MD simulation. In DC20:1PC lipid bilayers, 44 out of the 500 systems ended up with channel formation. In DC22:1PC bilayers, none of 1000 systems ended up with a channel. These results are consistent with the PMF results, which show a large (∼3 kcal/mol) kinetic barrier for channel formation in the DC22:1PC bilayer as compared to ∼1 kcal/mol in the DC18:1PC and DC20:1PC bilayers, as well as a shallower energy minimum for the dimer in DC22:1PC bilayer.

Figure 3

Unbiased MD simulations of gA channel formation. (a) Scatter plots for the I-DRMSD of the two gA monomers’ backbone Cα atoms at the end of the 100 ns MD simulations. In all systems, the simulations start with two separated gA monomers. In the DC18:1PC and DC20:1PC bilayers, 500 independent MD simulations were performed in each system; in the DC22:1PC bilayer, 1000 MD simulations were performed. The dotted lines denote the I-DRMSD values of 0.08 and 0.4 nm, with the I-DRMSD < 0.08 nm defined as the dimer state and the I-DRMSD > 0.4 nm as the monomer state. The transition state is defined with 0.08 ≤ I-DRMSD ≤ 0.4 nm. (b) Histograms for the data in (a). (c) Snapshots for the gA monomer, transition, and dimer states. To see this figure in color, go online.

From these unbiased MD simulations, we can estimate the rate of gA dimerization assuming that the gA dimerization proceeds through a single rate-limiting activation barrier (i.e., a two-state model). The probability of dimerization within a time Δt = 100 ns iswhere k = 1/τ is the rate and τ is the ensemble-averaged time needed for dimerization in the simulations. The estimated rates for gA dimerization in DC18:1PC and DC20:1PC bilayers are kDC18:1PC = 1.2 × 106 s−1 and kDC20:1PC = 0.9 × 106 s−1 at a gramicidin:lipid ratio of 1:100. The rate for gA dimerization in the DC22:1PC bilayer can be estimated via the relationship of , assuming that the pre-exponential factor A is the same in three lipid bilayers. Using kDC20:1PC and the -values for the DC20:1PC (1 kcal/mol) and DC22:1PC (3 kcal/mol) systems, kDC22:1PC ≈ 4 × 104 s−1, almost two orders of magnitude smaller than kDC18:1PC and kDC20:1PC. These simulated rate constants can be compared with previous experimental estimates for the rate constant for gA channel formation in DC18:1PC/n-decane planar bilayers (63,64), ∼1014 cm2 × mol−1 × s−1. Assuming the molar area of DC18:1PC is ∼0.7 nm2 (24), a gA/DC18:1PC molar ratio of 1:100 corresponds to 1 gA per ∼70 nm2 (1.4 × 1013 gA monomers per cm2) or 2.3 × 10−11 mol/cm2. The simulated second-order rate constant for gA dimer formation in n-decane-free DC18:1PC bilayers is ∼1.2 × 106 s−1/(2.3 × 10−11 mol/cm2) ≈ 5 × 1016 cm2 × mol−1 × s−1. Our simulation results thus predict association rate constants that are two orders of magnitude higher than the experimental estimates for the DC18:1PC/n-decane bilayers. n-Decane-containing DC18:1PC bilayers are ∼1.5 nm thicker than n-decane-free DC18:1PC bilayers (65,66), and one would expect the rate constant in n-decane-containing bilayers to be less than the one in n-decane-free bilayers. We also note that the experimental rate constant in DC18:1PC/n-decane bilayers is close to the one we estimate in DC22:1PC bilayers (∼2 × 1015 cm2 × mol−1 × s−1).

Process of gA dimerization

The gA dimerization process in the thin DC18:1PC and DC20:1PC bilayers can be captured with our unbiased MD simulations. To explore the gA dimerization process in the thick DC22:1PC bilayer, we employed an adaptive sampling method (67): we first did 100 independent 100-ns-long MD simulations and then used the configuration with the smallest gA monomer-monomer COM distance to spawn a new epoch of 100 independent MD simulations. This process was repeated until the two subunits formed a bilayer-spanning channel. Three epochs, totaling 30 μs of MD simulations, produced nine simulation trajectories (each 300 ns long) that led to gA channel formation (see Fig. S3).

For each of the three bilayers, one simulation trajectory was selected to investigate the details of gA dimerization. We first calculated the COM distance and the I-DRMSD between the two gA monomers, and the results are shown in Fig. 4 a. The variations of the COM distance and the I-DRMSD are consistent as two gA dimerize. Fig. 4 b shows the time evolution of the numbers of lipid acyl chain carbon atoms, water oxygen atoms, and gA2’s Cα atoms within a cylindrical volume (radius = 0.8 nm and height = 1.5 nm) below the COM of one monomer (gA1) (see Scheme S1). The results show that during gA dimerization, the number of water oxygen atoms and the gA2 Cα atoms surrounding the N-terminus of gA1 begin to increase as the numbers of lipid acyl chain carbon atoms begins to decrease, suggesting that the dimerization process involves rearrangement of the lipid tail packing. Fig. 4 c shows the time evolution of the bilayer hydrophobic thicknesses in bilayer regions where the lipid C1 atoms have different lateral distances (r1 ≤ 1.5 nm, the near region, and 3.5 < r2 ≤ 4 nm, the far region) from the COM of the two gA monomers (see Scheme S2). As the two gA monomers begin to associate (highlighted by the vertical dotted lines), the bilayer hydrophobic thickness near the monomers begins to decrease (with the bilayer thickness far from the channel remaining largely unchanged). That is, gA dimerization is coincident with the induction of membrane deformation adjacent to the gA monomers. Fig. 4 c further reveals that the bilayer thickness changes near the gA monomers are more dramatic in the DC22:1PC bilayer. The induction of larger membrane deformation and the associated increase of unfavorable membrane deformation energy would reduce the probability of gA dimerization in the thick DC22:1PC lipid bilayer, as observed in our unbiased MD simulations.

Figure 4

Representative simulations that started out with two gA monomers and resulted in a dimer are shown for each of the three different bilayer thicknesses. (a) gA monomer-monomer COM and I-DRMSD. (b) The change in lipid, water, and protein density in a volume below the gA1 monomer; see Supporting Materials and Methods and Scheme S1. (c) The time evolution of the bilayer thickness next to the channel (r ≤ 1.5 nm) and away from the channel (3.5 < r ≤ 4 nm). The vertical dashed lines denote the start of the gA dimerization process. To see this figure in color, go online.

Membrane deformation around gA monomers and dimers

To further characterize the membrane deformation, we divided the lipids in each bilayer leaflet into six different regions with different lateral distances (r ≤ 1.5, 1.5 < r2 ≤ 2, 2 < r3 ≤ 2.5, 2.5 < r4 ≤ 3, 3 < r5 ≤ 3.5, and 3.5 < r6 ≤ 4 nm) from the COM of the embedded gA monomer (see Scheme S2). As shown in Fig. 5, there is a marked increase of the bilayer’s hydrophobic thickness from the first to the second lipid region. The bilayer’s hydrophobic thickness reaches the unperturbed bilayer hydrophobic thickness at a distance of ∼2.3 nm from the monomers and ∼2.3, ∼2.5 and ∼2.8 nm from the dimer in the DC18:1PC, DC20:1PC, and DC22:1PC bilayers, respectively. The magnitude of the gA-induced bilayer deformation can be estimated as the difference between the bilayer thicknesses of the first (r1) and the sixth (r6) lipid regions, r6–r1 (the deformation of each leaflet will be half that value). Consistent with previous work (24,26), both gA monomers and dimers produce bilayer deformations that arise from a combination of interactions between the phospholipid headgroups and the polar groups around the mouth of the channel pore (both monomers and dimers; Fig. S4) and hydrophobic adaptation between the phospholipid acyl chains and the bilayer-spanning dimer. The monomer-induced deformations are almost identical (∼0.13 ± 0.01 nm per leaflet) in the three bilayers. The average dimer-induced leaflet deformations were 0.18 ± 0.01, 0.28 ± 0.01, and 0.43 ± 0.01 nm for the DC18:1PC, DC20:1PC, and DC22:1PC bilayers, respectively.

Figure 5

The region-specific lipid bilayer thickness for the three lipid bilayers with two gA monomers or a gA dimer. The dotted lines denote the bilayer thicknesses for pure DC18:1PC, DC20:1PC, and DC22:1PC bilayers. To see this figure in color, go online.

The formyl-N-termini unfolds partially during gA dimerization

Fig. 6 shows the root mean-squared fluctuation profiles for the Cα atoms of one gA monomer in the monomeric, transition, and dimeric states. The formyl-N-termini of the gA monomers in the monomeric and transition states are much more disordered than they are in the dimeric state. This partial unfolding at the formyl-N-termini may increase the effective length of the monomer and promote the initial polar interactions between the two formyl-N-termini that eventually lead to formation of the hydrogen-bond-stabilized dimer.

Figure 6

Root mean-squared fluctuation (RMSF) of the backbone Cα atoms and the C-terminal ethanolamine (ETA) in the monomer, transition, and dimer states. To see this figure in color, go online.

Reordering of water molecules during gA dimerization

The simulations reveal a number of water molecules at the end of the formyl-N-termini. To further characterize the behavior of water in the channel, we calculated the number of water molecules within 0.8 nm of the formyl Cα in the DC22:1PC bilayer. Fig. 7 a shows that water fills the gA pore even when gA is a monomer and that the number fluctuates between 0 and 6 (Fig. 7 b). For the gA dimer, there are consistently five or six water molecules in the channel. At the transition state for dimer formation in DC22:1PC, there were between 5 and 15 waters in the cavity between the two monomers (Fig. 7 b). Fig. S5 shows the number of bridging waters between the two formyl-N-termini during a channel formation simulation, which increases abruptly at the transition state. There is a clear difference in the relative orientation of water molecules in the dimer channel and water molecules that are associated with the gA termini at the transition.

Figure 7

Water molecules and their orientation in gA monomers, the transition state, and gA dimers. (a) Snapshots highlighting water molecules within the monomer, at the transition state, and within the dimer. We show only water molecules within 1.5 nm of the Cα of the formyl group. (b) The number of waters within 0.8 nm of the Cα of the formyl group of a gA monomer. (c) Orientation of the gA-associated water molecules in the three different states (black, green, and red for the monomer, transition, and dimer, respectively). Orientation is based on the dipole angle that each water forms with respect to the bilayer normal. To see this figure in color, go online.

To analyze water ordering, we calculated the angle between the dipole of each water molecule and the Z axis (between 0 and 180°). Because of strong interwater hydrogen bonds within the gA channel, the single-file waters are highly ordered, with angles of ∼30 or 150° (19). The water dipoles align with the channel; they are highly correlated with each neighboring water and can be relatively stable for tens of nanoseconds but also jointly fluctuate on the picosecond timescale (Fig. S6). Water within the gA monomer also has an orientation preference, but not as strong as within the dimer (Fig. 7 c). At the transition state, the water between the gA monomers has little orientational preference but a clear peak at 90° (Fig. 7 c).

gA channel activity

The main conclusion from the MD simulations, that the gA monomer ↔ dimer equilibrium is shifted toward the nonconducting monomer state as the bilayer thickness is increased, is consistent with experimental studies on the relation between the gA single-channel lifetimes and bilayer thickness in hydrocarbon containing bilayers (9, 10, 11,62). There is less information about the changes in the equilibrium per se in experimental systems comparable to those examined in the simulations.

We explored this question using a fluorescence quench method, which allows us to quantify the gramicidin channel activity in solvent-free bilayers of defined composition using stopped-flow spectrofluorometry (51). For these experiments, we used both gA and the naturally occurring mixture of gA, gB, and gC (see Experimental Materials), the same mixture we have previously used in fluorescence quench experiments (52,68, 69, 70, 71, 72, 73).

The experimental setup is illustrated in Fig. 8 a. LUVs formed with different gramicidin:lipid mole fractions and encapsulating a water-soluble fluorophore, ANTS, were mixed with a solution containing a gramicidin-channel-permeant ion (the thallous ion, Tl+) that can quench the ANTS fluorescence in a stopped-flow spectrofluorometer, and the time course of fluorescence quench (caused by the influx of Tl+ through conducting gramicidin channels) is recorded. The rate of quenching at 2 ms (the instrumental dead time is ∼1.5 ms) was determined by fitting a stretched exponential decay to the first 200 ms of the fluorescence quench trace (see Gramicidin Fluorescence Assay).

Figure 8

Experimental studies on gramicidin channel function in bilayers of different thickness. (a) Schematic description of the stopped-flow fluorescence quench experiments. An LUV that encapsulates the reporter fluorophore ANTS, with gramicidin monomers and dimers embedded in the membrane is shown at the bottom of the panel. Dark blue LUVs represent LUVs that have not been quenched, and light blue LUVs depict LUVs that have been quenched because of Tl+ influx through bilayer-spanning gramicidin channels. (b and c) Fluorescence quench traces obtained with DC18:1PC, DC20:1PC, and DC22:1PC LUVs and the indicated gD:lipid ratios (b) or gA:lipid ratios (c). The gray envelopes show the results from a single mixing reaction, and the colored traces depict the average of the mixing reactions in that experiments. Table 1 summarizes the results. To see this figure in color, go online.

Table 1

Fluorescence Quench Rates at Different Gramicidin:Lipid Molar Ratios

Phospholipid	Gramicidin/Lipida	Gramicidin/LUVab	Quench Ratec (s−1)	ΔFluorescencec	n
gD
DC18:1PC	1:120,000	∼1	69.1 ± 1.9	0.14 ± 0.01	3
DC18:1PC	1:80,000	∼2	73.7 ± 9.4	0.14 ± 0.01	3
DC18:1PC	1:40,000	∼4	99 ± 15	0.17 ± 0.03	9
DC20:1PC	1:30,000	∼5	75 ± 13	0.11 ± 0.01	3
DC20:1PC	1:20,000	∼8	100 ± 5.9	0.14 ± 0.01	3
DC20:1PC	1:10,000	∼15	101 ± 21	0.17 ± 0.01	6
DC22:1PC	1:4,000	∼38	6.4 ± 2.2	0.15 ± 0.01	3
DC22:1PC	1:2,000	∼75	19.2 ± 5.7	0.23 ± 0.06	28
DC22:1PC	1:700	∼210	75.3 ± 7.5	0.18 ± 0.01	3
gA
DC18:1PC	1:120,000	∼1	26.8 ± 2.8	0.15 ± 0.02	3
DC18:1PC	1:80,000	∼2	54.8 ± 8.0	0.17 ± 0.01	3
DC18:1PC	1:40,000	∼4	71.5 ± 5.2	0.18 ± 0.02	3
DC20:1PC	1:30,000	∼5	105.9 ± 4.6	0.20 ± 0.02	3
DC20:1PC	1:20,000	∼8	127.9 ± 4.1	0.21 ± 0.02	3
DC20:1PC	1:10,000	∼15	133.8 ± 6.7	0.21 ± 0.01	3
DC22:1PC	1:6,000	∼25	17.0 ± 27.4	0.14 ± 0.04	6
DC22:1PC	1:4,000	∼38	28.9 ± 1.8	0.11 ± 0.00	3
DC22:1PC	1:2,000	∼75	49 ± 22.9	0.14 ± 0.02	12

gD denotes the naturally occurring mixture of gA, gB, and gC. gA denotes [Val1]gramicidin A.

Assuming 150,000 phospholipids/LUV.

Mean ± standard deviation.

Fig. 8, b and c shows results with LUVs formed from DC18:1PC, DC20:1PC, and DC22:1PC and the indicated gD:lipid (b) or gA:lipid (c) mole ratios. Table 1 summarizes the results. As expected, the molar ratio of gramicidin:lipid required for observing comparable quench rates increases with increasing bilayer thickness, and the gramicidin concentration in the thickest DC22:1PC bilayer is sharply increased compared to DC20:1PC and DC18:1PC bilayers, is in agreement with the PMF results in Fig. 1. Overall, gD and gA showed similar results, but gA has somewhat higher potency. Comparison of fluorescence quench traces at different gD:lipid ratios in each bilayer are shown in Fig. S9.

From our experimental results, we can estimate the bilayer thickness contribution to the changes in the thermodynamic equilibrium distributions of gramicidin monomers versus dimers,where M denotes the nonconducting monomers, D the conducting dimers, and KA the association constant for channel formation. The distribution between monomers and dimers is given bywhere the curly brackets denote surface concentrations (mol/unit area) and the total gramicidin concentration in the membrane is given by {T} = {M} + 2 × {D} (74).

Equation 7 reduces to {D} ≈ {T}2 × KA (and {M} ≈ {T}) in the limit {T} × KA → 0 and to {D} = {T}/2 (and {M} ≈ 0) in the limit {T} × KA → ∞. In the former case, in which there is an excess of nonconducting monomers, it is possible to estimate the bilayer contribution to the free energy of gramicidin dimerization ; in the latter case, the monomer ↔ dimer equilibrium is shifted so strongly to the conducting dimers that it precludes estimating . The results in Fig. 8, b and c were obtained under conditions in which the LUV preparations had comparable time-averaged numbers of conducting channels/LUVs, as evident by the quench rates and changes in fluorescence signal (ΔF = F(0) – F(∞)); see also Table 1.

Given the gramicidin:lipid ratio in the DC18:1PC LUVs whether the gramicidin was gD or gA, there can be at most two conducting channels per DC18:1PC LUVs (at 1:40,000 gramicidin:lipid), and only a fraction of the LUVs will have a conducting channel at 1:120,000 gramicidin:lipid, as evident also from the reduction in ΔF. These systems were in the {T} × KA → ∞ limit, which precludes estimating . Increasing the gramicidin:lipid ratio therefore reduced the fraction of LUVs with no conducting channels, with less effect on the initial quench rate.

At the other extreme, the gramicidin:DC22:1PC LUVs were in the {T} × KA → 0 limit. In the gD experiments, the quench rate at a gramicidin:lipid molar ratio of 1:700 was similar to that observed with the gramicidin:DC18:1PC LUVs at 1:120,000 gramicidin:lipid, meaning that the two populations have similar time-averaged numbers of conducting channels per LUV. In the gA experiments, the quench rate at a gramicidin:lipid molar ratio of 1:2,000 was similar to that observed with the gramicidin:DC18:1PC LUVs at 1:40,000 gramicidin:lipid. For either gramicidin, the quench rate increased three- to-fourfold when the gramicidin:lipid ratio increases threefold, less than the expected increase for a monomer ↔ dimer reaction. We do not understand the reason for this discrepancy.

The gramicidin:DC20:1PC LUVs are intermediate. As was the case for gramicidin:DC18:1PC LUVs, there was little change in quench rate when the gramicidin:lipid ratio was increased threefold, whether the gramicidin was gA or gD, and a modest decrease in the fraction of LUVs that had no conducting channels during the 1 s quench experiment (increase in ΔF).

The measured time course of fluorescence quenching reflects an average of the fluorescence quench time courses for the LUVs in the system. For a given LUV, the instantaneous quench rate depends on the Tl+ influx (actually the gramicidin channel-mediated Tl+/Na+ exchange), which varies with the fluctuating number of conducting channels in the LUV membrane (which varies with the rate constants for channel formation and dissociation and the number of gramicidin subunits in the LUV membrane). The measured quench rate thus varies with the time-averaged number of conducting gramicidin channels per LUV in the sample.

Thus, if the gramicidin:DC20:1PC LUVs were in the {T} × KA → 0 limit, we can estimate the magnitude of the changes in going from DC20:1PC to DC22:1PC from the gramicidin:lipid ratios that give similar quench rates in the two systems. When that is the case, the time-averaged number of gramicidin channels in the two systems (time-averaged number of channels/LUV) will be equal andorandor

Inserting the gramicidin:lipid ratios from Table 1, we can estimate to be between 1.9 kcal/mol (using the 1:10,000 gA:DC20:1PC and 1:2,000 gA:DC22:1PC results) and 4.5 kcal/mol (using the 1:30,000 gD:DC20:1PC and 1:700 gD:DC22:1PC results). The average value—based on the 1:10,000 gA:DC20:1PC and 1:2,000 gA:DC22:1PC results, the 1:10,000 gD:DC20:1PC and 1:700 gD:DC22:1PC results, the 1:30,000 gA:DC20:1PC and 1:2000 gA:DC22:1PC results, and the 1:30,000 gD:DC20:1PC and 1:700 gD:DC22:1PC results—is 3.2 ± 1.0 kcal/mol (mean ± standard deviation), in good agreement with the results from the PMF calculations.

The conformational preference of gramicidin varies with the channel-bilayer hydrophobic mismatch: in DC18:1PC and DC20:1PC membranes, the gramicidins are predominantly monomeric and form the conventional β6.3-helical dimers (75,76); in DC22:1PC membranes, however, ∼40% of the gramicidins form helical dimers with at least two different conformations (77), and at least some of these form conducting channels. This conformational polymorphism has consequences for our estimate of because first, the ratio of monomeric gramicidin to DC22:1PC will be only ∼60% of the total gramicidin:DC22:1PC ratio, which means that the above reasoning will overestimate by ∼0.6 kcal/mol; second, some of the Tl+ influx will be through double-stranded gramicidin channels, which means that we have overestimated the time-averaged number of β6.3-helical channels and therefore underestimated the magnitude of .

The experimental results do not allow for estimates of the parameters in the continuum model, Eq. 4, because only the gramicidin:DC22:1PC LUVs were in the {T} × KA → 0 limit. Given the similarity in the computational and experimental estimates for , however, HB is likely to be comparable to the estimate obtained from the analysis of the MD results or about twofold less than the estimates obtained from analysis of the changes in gA single-channel lifetimes as function of bilayer thickness (12,56), ∼16.5 kcal/mol × nm2. The difference between these estimates most likely reflects that the MD-based estimate of is the difference between and ( ), where is likely to differ from 0, whereas the analysis based on single-channel lifetimes focused solely on , effectively assuming that .

Discussion

Lipid bilayers are fluid and deformable at the molecular level. The lipid molecules can adjust their positions and conformations to pack around and interact with membrane proteins. Whereas the bulk properties of lipid bilayers can be characterized with many techniques, understanding the molecular-level interactions of individual lipids and single protein channels remains challenging. We have used extensive MD simulations to quantitatively probe the effect of lipid bilayer thickness on gramicidin dimerization. Our simulation results are in good agreement with our experimental results and provide an in-depth thermodynamic description for the process of gramicidin dimerization. Overall, there is exquisite coupling between lipid-protein-water, which results in a complex mechanism for dimerization, even for a channel as simple as gramicidin.

gA dimers were stable on the timescale accessible to MD simulations in all three lipid bilayers. To quantify the channel stability, we calculated the PMF for dimer → monomer transition. We found that a subnanometer change of 0.4 nm in bilayer thickness led to a free-energy change of 3–5 kcal/mol for the gA dimer → monomer transition. The PMFs also revealed a bilayer-thickness-dependent activation barrier for gA dimerization. The barriers for gA dimerization in DC18:1PC and DC20:1PC bilayers were similar (∼1 kcal/mol), whereas the barrier in DC22:1PC bilayers was sharply increased to ∼3 kcal/mol. The activation barrier arises from the energetic penalty for lipid bilayer deformation surrounding the two gA monomers as they start to associate. The rate constant for gA dimerization in the DC22:1PC bilayer is estimated to be ∼2 orders of magnitude slower than the DC18:1PC and DC20:1PC systems. However, the simulation-predicted rates for gA dimerization in the three bilayers may be overestimated because of lack of sampling of gA monomer rotations in our MD simulations. It has been suggested that the gA channel formation may involve gA monomer rotation (78). The effect of gA monomer rotation on the kinetics of gA channel formation will be investigated in our future work.

Detailed analysis of the simulations revealed a number of contributing factors to gA dimerization. The bilayer perturbations around the gA channel are larger in the thicker DC22:1PC bilayer than the thinner DC18:1PC and DC20:1PC bilayers, which accounts for the large shift in gA channel activity observed in our experiments. Significant reordering of lipid tails between two gA monomers in opposing leaflets is necessary for channel formation. In the monomeric and transition states, the bilayer-embedded N-terminus of gA can partially unfold, dipping further into the bilayer hydrophobic core. Therefore, when two gA monomers are in close proximity, polar interactions between two partially unfolded gA N-termini promote channel formation. Additionally, when two gA monomers come close, water enters the cavity between them and helps stabilize the transition state.

The PMF results were validated by the fluorescence-quenching experiments performed in this work, which showed that a sharply increased gramicidin concentration is required to have the approximately same number of dimers in the DC22:1PC bilayer as compared to the thinner DC18:1PC and DC20:1PC bilayers, and the experimental estimates for change in bilayer deformation between DC20:1PC and DC22:1PC bilayers are in good agreement with the energies deduced from the MD simulations.

Conclusion

Extensive MD simulations, together with experiments, delineate a consistent mechanism for bilayer thickness influencing gA channel stability and the gA dimerization process. The gramicidin channel has been developed as an experimental tool to probe lipid bilayer properties, including the subtle changes in bilayer thickness induced by exogenous molecules (69,71), and our characterization of gA dimerization with atomistic detail allows for deeper understanding of gramicidin channels as probe to for exploring the effects of changing membrane properties. The atomic details of dimerization showed the concerted action of lipid packing, water coordination, and protein structural rearrangement, demonstrating the intrinsic complexity of protein-bilayer coupling and the need for further exploration.

Author Contributions

D.S., W.F.D.B., and H.I.I. designed and performed the simulations. T.A.P. and O.S.A. designed and performed the experiments. D.S., W.F.D.B., T.A.P., O.S.A., F.C.L., and H.I.I. wrote the manuscript.