Atomistic Characterization of Gramicidin Channel Formation.

We investigated gramicidin A (gA) subunit dimerization in lipid bilayers using microsecond-long replica-exchange umbrella sampling simulations, millisecond-long unbiased molecular dynamics simulations, and machine learning. Our simulations led to a dimer structure that is indistinguishable from the experimentally determined gA channel structures, with the two gA subunits joined by six hydrogen bonds (6HB). The simulations also uncovered two additional dimer structures, with different gA-gA stacking orientations that were stabilized by four or two hydrogen bonds (4HB or 2HB). When examining the temporal evolution of the dimerization, we found that two bilayer-inserted gA subunits can form the 6HB dimer directly, with no discernible intermediate states, as well as through paths that involve the 2HB and 4HB dimers.

The linear gramicidins are antimicrobial peptides produced by the soil bacterium Brevibacilus brevis, which effectively kills Gram-positive bacteria. The antimicrobial mechanism has been linked to the formation of transmembrane ion-conducting channels that increase the membrane permeability to ions and water,[1] and gramicidin channels have been studied extensively as prototypical ion channels.[2] The major component of the natural mixture of gramicidins is gramicidin A (gA), with the amino acid sequence of formyl-l-Val1-d-Gly2-l-Ala3-d-Leu4-l-Ala5-d-Val6-l-Val7-d-Val8-l-Trp9-d-Leu10-l-Trp11-d-Leu12-l-Trp13-d-Leu14-l-Trp15-ethanolamine. The alternating l- and d-amino acids and the peptide sequence caps, a formyl group at the N-terminus and an ethanolamine group at the C-terminus, allow gA to fold into a β6.3-helical conformation in lipid bilayers, with the Trp-rich C-terminus localized at the lipid bilayer/water interface and the formyl-N-terminus localized at the bilayer center. The monovalent cation-selective channel is formed by the transmembrane dimerization of two gA subunits. The dimer is stabilized by a maximum of six hydrogen bonds between the 1Val1-2Ala5, 1Ala3-2Ala3, and 1Ala5-2Val1 amino acid pairs in the two subunits.[3−5] Atomic resolution gA dimer structures (PDB: 1GRM,[6]1JNO,[5]1MAG[7]) show six hydrogen bonds between the two gA monomers. This structure will be denoted 6HB dimer, but little is known about the dimerization process per se. Many simulations studies have been reported on the 6HB channel, examining structural stability, water wire, ion permeation, and lipid packing.[3,8−15] Monte Carlo simulations by Miloshevsky and Jordan[16] suggested the existence of intermediate gA dimer structures with four or two hydrogen bonds, denoted 4HB and 2HB dimers, respectively. This computational prediction was supported by temperature-jump infrared spectroscopy studies by Stevenson and Tokmakoff.[17] The molecular structures of these intermediate gA dimers, as well as their role in 6HB channel formation, however, remain poorly understood. In this work, we have run an aggregate of 1356 μs all-atom molecular dynamics (MD) simulations to investigate the mechanism and process of gA subunit dimerization in lipid bilayers of different thickness.

Results and Discussion

To elucidate the molecular details for the 6HB gA dimer (channel) formation, we first performed 156 μs all-atom replica-exchange umbrella sampling (REUS) simulations. The configurational space of two trans-bilayer gA subunits was explored in three lipid bilayers with linearly 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). Figure (a) shows two simulation snapshots of gA monomers and the dimer in the DC18:1PC bilayer. From the REUS simulations, we obtained the one-dimensional potential of mean force (PMF) profiles for the 6HB dimer → monomer transition in the three lipid bilayers using a reaction coordinate of the gA–gA center-of-mass (COM) distance (dgA–gA), see Figure S1 for PMF convergence analysis. Of note, the one-dimensional PMF profiles can be used to derive two-dimensional PMF maps with the dgA–gA and the Z-direction (the direction normal to the bilayer surface) component of dgA–gA as two reaction coordinates (see Figure S2). All PMF profiles shown in Figure (b) exhibit two V-shaped basins centered at ∼1.3 and ∼1.5 nm along the dgA–gA reaction coordinate. Inspection of the REUS simulation trajectories shows that the REUS umbrella windows, which centered around the two free energy minima, were dominated by the 6HB and 4HB dimer configurations. The 6HB dimers are stabilized by six hydrogen bonds between the 1Val1-2Ala5, 1Ala3-2Ala3, and 1Ala5-2Val1 amino acid pairs at the dimer interface (Figure (c)) with the interformyl distance being 0.45 ± 0.04 nm and dgA–gA being 1.33 ± 0.01 nm. The channel pore supports a continuous water wire owing to the proximity of two formyl groups, which seals the channel gate at the dimer interface (Figure (c)). The simulated 6HB channels are comparable to the experimentally determined structure (PDB: 1JNO, see Figure S3). In the 4HB dimer, the two gA monomers are joined by a maximum of four hydrogen bonds between the 1Val1-2Ala3 and 1Ala3-2Val1 amino acid pairs and the average interformyl distance, and dgA–gA’s have increased to 0.81 ± 0.05 and 1.48 ± 0.01 nm, respectively. The increased interformyl distance causes disruption of the order of water wire at the dimer interface (Figure (c)). No 2HB dimer structure was identified in the REUS simulations.

Figure 1

(a) Snapshots of two gA monomers, gA1 and gA2 (cyan and orange, respectively), and the dimeric channel in a DC18:1PC lipid bilayer. (b) PMF profiles for the gA dimer → monomer transition in DC18:1PC, DC20:1PC, and DC22:1PC bilayers obtained from REUS simulations. The reaction coordinate is the COM distance between two gA monomers, and the COM distance is defined using only the backbone Cα atoms. The three PMF profiles were manually shifted to make the PMF to be 0 at dgA–gA = 3.9 nm, and no Jacobian correction[18] was made to the PMF profiles. (c) The REUS simulations demonstrated the presence of 6HB and 4HB dimers, evident by the two local free energy minima in the PMF profiles. The unbiased MD simulations demonstrated not only the presence of 6HB and 4HB dimers but also the less stable 2HB dimers. The formyl groups were drawn using the licorice in the VMD software.[19] (d) Schematic illustration of MD simulations with different monomer–monomer orientations: the face-to-face (FF) and face-to-back (FB) orientations (lipids and water outside the monomers are not shown).

The PMF profiles show that the 6HB dimer is thermodynamically more favorable than the 4HB dimer in all three bilayers. In the thin DC18:1PC bilayer, the 6HB dimer is 3.5 kcal/mol more favorable than the 4HB dimer, indicating that the equilibrium (time-averaged) number of 6HB dimers is about 300-fold larger than the number of 4HB dimers, which compares well with the experimental results of Sigworth and Shenkel[20] that the channel is 60–200 times more likely to dissociate from a sublevel state than from the major (6HB) state.

The PMF results thus are consistent with the experimental findings that the dominant gA channels are 6HB dimers.[5,7,21] To form the channel, the two gA monomers have to produce a local bilayer thinning and therefore overcome a bilayer-imposed transition barrier (ΔG‡). In DC18:1PC and DC20:1PC bilayers, the ΔG‡ estimates are very close (∼0.8 kcal/mol), whereas ΔG‡ is increased to ∼2.2 kcal/mol in the DC22:1PC bilayer. ΔG‡ thus increases nonlinearly with the linear increase of acyl chain length, in agreement with the predictions of continuum elastic models of protein-induced bilayer deformations.[22,23] The rate constant for gA dimerization, k, is associated with ΔG‡ via the relationship of k = A × e–Δ, with A being a system-dependent frequency prefactor. The PMF results suggest that the rate constant for gA dimerization in DC22:1PC bilayers is reduced by a factor of ∼10 compared to DC20:1PC bilayers, assuming the frequency factor A is similar in the three lipid bilayers, which is, in general, in agreement with the observation that an increase in the acyl chain length of two CH2 groups means that one needs to add 10 times more gA in order to observe the same channel appearance rate.[24]

To further probe the 6HB gA dimer formation, we ran a total of 1.2 ms unbiased all-atom MD simulations. The starting configurations for the unbiased MD simulations are two bilayer-incorporated gA monomers with different relative orientations, Figure (d): a face-to-face (FF) orientation and a face-to-back (FB) orientation. We ran 200 independent 1 μs-long MD simulations in each bilayer system with both monomer–monomer orientations and extracted 12 million simulation snapshots of the gA monomers. To identify channel-like structures, we used the k-means clustering approach to classify the 12 million snapshots into 100 groups using monomer–monomer Cα atom contact maps as features (see Simulation Methods in the Supporting Information for details). From these 100 groups, we were able to identify the 6HB and 4HB dimer structures, as well as the 2HB dimer structure (Figure (c)). In the 2HB dimer, the gA monomers are joined by hydrogen bonds between the 1Val1-2Val1 pair, with the average interformyl distance of 0.72 ± 0.06 nm and dgA–gA of 1.52 ± 0.03 nm. The small difference in the COM distance between 4HB and 2HB dimers (less than the distance between two neighboring umbrella windows of 0.05 nm) explains why the REUS simulations using the three-dimensional COM distance as the single reaction coordinate failed to capture the conversion from the relatively stable 4HB dimer to the much less stable 2HB dimer within the simulation time. Multidimensional REUS simulations, i.e., using the lateral and vertical components of the COM distance as two distinct reaction coordinates, may overcome this sampling problem, though the use of two reaction coordinates is computationally demanding in REUS simulations.

To identify all gA dimerization events during the unbiased MD simulations, we implemented a supported vector machine (SVM) model. The SVM model works by assigning a given simulation snapshot to one out of five different states: monomer (M), monomer → initial dimer transition (T), 2HB, 4HB and 6HB, such that the high-dimensional MD simulation trajectories were mapped onto the five-state trajectories (see Simulation Methods in the Supporting Information for details). The resultant 1200 five-state simulation trajectories are shown in Figure S4. Briefly, Figure S4(a–f) shows each of the five-state simulation trajectories initiated with different starting gA configurations (face-to-face vs face-to-back) in the DC18:1PC, DC20:1PC, and DC22:1PC bilayers, respectively. All 200 simulations for each condition are plotted together in Figure S4(g). It can be seen from Figure S4 that the MD simulations produced a few numbers of 6HB, 4HB, and 2HB channels. We observed frequent transitions between the T and 2HB states (see Figure S4(a-146) for example), confirming that the 2HB dimer is short-lived. We also observed 4HB channel dissociation (see Figure S4(a-19, a-46, a-73, a-119, a-196, and d-190)) within the 1 μs simulation time, whereas all 6HB channels that formed remained stable through the remainder of that simulation. To summarize the dimer formation during the 1200 MD simulations, we plotted in Figure S5 the time-dependent fractions of gA in the different states. The number of 4HB and 6HB dimers formed at the end of the 1200 × 1 μs-long unbiased MD simulations is summarized in Table . Almost all 6HB dimers were formed directly by the dimerization of two gA monomers via the M → T → 6HB pathway; in seven trajectories, the 6HB dimer was formed via intermediate states (the 2HB or 4HB dimers). To unveil details of 6HB dimerization with and without the intermediate states, we used the five-state simulation trajectories to identify the dominant M → T → 6HB (Figure (a)) and the M → T → 4HB → 6HB (Figure (b)) pathways, shown along with selected snapshots. In the M → T → 4HB → 6HB pathway, the 4HB dimer sometimes dissociated (Figure (b)), and the two gA monomers could further rotate to form the 6HB dimer. In both pathways, gA dimerization involved an orientational/rotational movement of the two gA monomers, as inferred from the rotational root-mean-squared displacement (RRMSD) of the monomers. We also analyzed the four tryptophan residues’ χ1 and χ2 dihedral angle change during the two distinct dimerization pathways, and the results are shown in Figure S6. We observed an abrupt change of the dihedral angles in the 4HB and 6HB channel states; however, no such change was found during the T → 6HB, T → 4HB, and 4HB → 6HB transition process, suggesting that the tryptophan’s orientation change is not a determining factor for the gA dimerization. Figure S7 shows the χ1 and χ2 tryptophan dihedral angle distributions for gA in three lipid bilayers separated by state, and no clear trend can be seen between the monomer vs dimer states. When the starting orientation of the two gA monomers was changed from FF to FB, the number of 6HB dimer formed at the end of the MD simulations decreased from 74 to 26 in DC18:1PC bilayers and from 55 to 11 in DC20:1PC bilayers, whereas the number of 4HB dimers increased from 21 to 29 and from 26 to 43, respectively. Although most 4HB dimers remained stable throughout the 1 μs-long simulations, they were likely to dissociate over longer time scales, as determined from the number ratio of 4HB and 6HB dimers obtained in the unbiased MD simulations, which was much higher than that estimated from the PMF simulations.

Table 1

Numbers of 6HB and 4HB Dimer Formed via Different Pathways in Three Bilayer Systems

	DC18:1PC		DC20:1PC		DC22:1PC
system	FF	FB	FF	FB	FF	FB
M → T → 6HB	71	23	55	10	2	0
M → T → 2HB → 6HB	1	3	0	1	0	0
M → T → 4HB → 6HB	2	0	0	0	0	0
M → T → 2HB → 4HB	0	0	0	0	0	0
→ 6HB
M → T → 4HB	18	25	26	38	1	2
M → T → 2HB → 4HB	3	4	0	5	0	0
M → T → 6HB → 4HB	0	0	0	0	0	0

Figure 2

Five-state trajectories for the (a) M → T→ 6HB and the (b) M → T→ 4HB → 6HB pathways. The five different macrostates are monomer (M), monomer →initial dimer transition (T), 2HB, 4HB, and 6HB dimers. The gA dimerization process involves orientational/rotational movements of the two gA monomers relative to each other, which can be characterized by the rotational root-mean-squared displacement (RRMSD) of the two gA monomers with the 6HB dimer as the reference structure. (c) PMF profiles for a single K+ ion permeating the 6HB and 4HB channels. Snapshots of each channel structure are shown with the ion in the middle of the channel.

The PMF for a single K+ translocating along the dimer pore (Figure (c)) suggests that both 6HB and 4HB dimers can conduct ions, with similar free energy barriers relative to the bulk aqueous phase (8.2 kcal/mol in the 4HB channel and 8.5 kcal/mol in the 6HB channel). Nevertheless, considering that the equilibrium number of 6HB dimers in DC18:1PC bilayers is around 300-fold higher than that of 4HB dimers, as estimated from the PMFs in Figure (b), the 6HB dimers will be the dominant ion-conducting gA channels. Sigworth and Shenkel’s experimental work[20] reported a low-conductance sublevel channel state, which was characterized as the incomplete closure of the channel and rapid channel dissociation. The sublevel state may reflect the formation of 4HB dimers. Our simulated PMFs would suggest the two types of channels have similar ion conductance (Figure (c)) but should only be used as a qualitative comparison. The work by Paulino et al.[15] suggested that the hydrogen bond interactions between water and the carbonyl groups lining the gA channel are underestimated using classical force fields, hinting that the quantum treatment of the water wire may lead to improved results including the PMF for ion permeating the gA channel. Additionally, the computational work by Allen et al.[14] and Peng et al.[25] has shown that inclusion of electronic polarizability would dramatically lower the energy barrier for K+ permeating the gA channel by as much as 3 kcal/mol, and it remains to be investigated whether a polarizable force field would predict a big difference in the PMFs for K+ permeating across the 6HB and 4HB channels.

In the unbiased MD simulations, the total number of formed 6HB dimers decreases as the bilayer thickness is increased (see Table and Figure S5). Particularly, the number of gA channels formed in the thickest DC22:1PC bilayer is sharply decreased compared to the two thinner bilayers. From the unbiased MD simulations, the rate constant for 6HB gA formation (kA) can be approximated assuming the dimerization is dominated by a single activation barrier. Two related methods were used to estimate kA for 6HB dimer formation:[26,27] (a) kA can be estimated from p(Δt) = 1 – e–Δ,[26] where Δt = 1 μs is the simulation time for each MD simulation, τ is the ensemble-averaged time needed for dimerization, and p is the probability of 6HB channel formation within the 1 μs-long simulations. For one monomer dimerizing with another in the simulations, k = τ/C (where the gA monomer concentration C is 2.3 × 10–11 mol/cm2).[26] (b) When Δt ≪ τ, k can be estimated from the total number of 6HB dimers formed during the simulations (N6HB) divided by the total amount of elapsed MD simulation time before the 6HB dimer formation event occurs (tf) scaled by C.[27] The rate constant estimates are listed in Table . The rate constant for gA channel formation in the DC22:1PC bilayer is reduced ∼60- and ∼40-fold compared to the DC18:1PC and DC20:1PC bilayers, respectively. This steep reduction in rate fits the ∼10-fold increase in gA concentration needed to maintain a similar gA appearance rate for every two CH2 added to the lipid tails.[24] The simulation-derived k6HB estimates in DC18:1PC bilayers are about 2 orders of magnitude larger than experimental estimates in DC18:1PC/n-decane bilayers,[28,29] as n-decane increases bilayer thickness by ∼1.5 nm,[30,31] meaning that the dimerization rate is expected to be higher in the n-decane free bilayers used in the simulations. We cannot exclude, however, that the discrepancy may be due to the relatively shorter MD simulation time scale as compared to experimental conditions; though most of the 4HB intermediate dimers remain stable in the one microsecond simulations, they would not be detectable experimentally. The simulation-derived k should be considered an approximation of the 6HB dimerization rate through the M → T → 6HB pathway. Additionally, we caution that the gA dimerization process in lipid bilayers is mainly driven by the polar interactions embedded in the hydrophobic lipid bilayer center, and therefore employing a polarizable force field to treat the gA/bilayer system is expected to improve the results.

Table 2

Estimates for the Number (N), Time (tf), and Rate Constant (k) for the 6HB Dimer Formation in the Three Lipid Bilayersa

system	N6HB	tf (μs)	kA(a) (×1016)	kA(b) (×1016)
DC18:1PC	100	336.73	1.25	1.29
DC20:1PC	66	354.00	0.78	0.81
DC22:1PC	2	399.43	0.02	0.02

The k’s were predicted with two different methods, defined as (a) and (b) in the main text, and the unit for the k’s is cm2 × mol–1 × s–1.

The 400 × 1 μs unbiased MD simulations performed in the thick DC22:1PC bilayer produced only two trajectories that led to 6HB dimer formation through the M → T → 6HB pathway. Figure plots the end-to-end distance between the gA monomers in one of these trajectories, which shows that the two formyl-NH-termini are partly unwound before the monomers dimerize. This structural flexibility effectively increases their length and promotes polar interactions between the two N-termini (see snapshot in Figure ). The flexibility of the monomers is reduced once the 6HB dimer is formed.

Figure 3

Five-state trajectories for 6HB dimer formation in the DC22:1PC bilayer. The two monomers (gA1 and gA2) can partially unfold their N-terminus (shown as VDW beads) and thereby promote the monomer–monomer association as shown by the monomer end-to-end distance plots.

In summary, we have elucidated the molecular features of gA channel formation in lipid bilayers of different thickness by running extensive classical force field-based MD simulations (156 μs REUS and 1.2 ms unbiased MD simulations). The simulations unveil three different dimer substates and that the most stable 6HB dimers can be formed by two major pathways: two gA monomers can either form the 6HB dimers directly or form the 6HB dimers via intermediate 2HB or 4HB dimers as hubs. These simulations demonstrate the complexity of protein–protein and protein–lipid coupling and time scales necessary to capture them even for a small membrane protein like the gA channel.