pH-Responsive peptide nanopores are stabilized by lipid and water-mediated hydrogen bonding networks

Abstract

Membrane-spanning nanopores that allow controlled passage of macromolecular cargo across cell membranes can empower many biomedical applications. Such nanopores are formed, in a pH-responsive manner, by the synthetically evolved “pHD peptide” family. pHD peptides fold into amphipathic α-helices, but have many charged and polar residues and are thus not predicted by classical hydropathy analyses to fold into membrane-spanning structures. Yet, when the pH is below ∼6, pHD peptides readily self-assemble into nanopores, even at low concentration. Knowledge of the molecular structure of the pHD peptide pore is needed for further rational design and optimization of nanopore-forming activity targeted to specific membranes and pH conditions. To this end, we have carried out extensive atomistic molecular dynamics simulations to explore the protonation-dependent structure and dynamics of nanopores created by the peptide pHD108. Simulations and graph-based analyses of hydrogen bonding reveal that, in the nanopore, the numerous carboxylate and carboxyamide sidechains form a dense, water-bridged H-bond network across the bilayer. In this network, direct H-bonds between neighboring peptides are few. Instead, the network is dominated by water-bridged intrapeptide interactions and by water-bridged interactions with the headgroups of many lipid molecules with unusual conformations and orientations. The lipids in the H-bond network make critical contributions to nanopore stabilization. These studies reveal a non-classical means of stabilizing nanopores in bilayers formed by highly charged peptides, creating an avenue towards engineering of membrane-embedded structures.

---

Introduction

Many peptides that permeabilize membranes via a wide variety of activities, mechanisms, and inferred structures have been described in the literature.¹ Of these, peptides that form macromolecule-sized, aqueous nanopores are very rare, likely because the free energy penalty for a water-filled pore scales with the area of the pore-bilayer boundary.^2,3 A large macromolecular pore can thus form only if created by peptides that minimize the free energy penalty for such a structure. The fundamental principles of engineering peptides that self-assemble into nanopores are currently not known. However, we previously used synthetic molecular evolution, multiple generations of iterative library design and high-throughput screening, to select unique peptides that enable the passage of macromolecules of hydrodynamic radii up to 10 nm, which requires a nm-scale aqueous pathway through the membrane, created by the peptides.^4–7 These nanopore-forming peptides, which we call pHD peptides, are active at very low concentrations, but only when pH < 6.⁸ By delivering uptaken protein cargoes to the cell cytosol, we have shown that pHD peptides also self-assemble into nanopores in the endosomal membranes of living cells.⁹ As published functional experiments lack direct information on how the pHD peptides stabilize nanopore-bilayer interfaces, here we explore the mechanism with extensive atomistic molecular dynamics (MD) simulations and graph-based analyses of dynamic interaction networks sampled by the peptides.

The pHD peptides are highly amphipathic, with a polar and highly charged helical face that has five acidic glutamate or aspartate residues, two basic histidine or lysine residues, and three polar glutamines. An opposing nonpolar face is comprised mostly of aliphatic residues. Even when the multiple acidic sidechains are deprotonated, and thus charged, the pHD peptides can insert into membrane-spanning orientations and form nanopores^4,5,7 which is surprising because they are highly polar, overall.¹⁰ With unmodified pHD peptides, the self-assembly of the pHD peptides into the unique nanopores occurs upon a slight acidification to pH ≤ 6, likely due to protonation of just one or two glutamic acid sidechains.^4,5,7 pHD peptide activity, defined as the ability to form pores as a function of pH, scales only with membrane binding,⁵ which shows that pores can form when all carboxylic groups are negatively charged as long as some peptide is bound. This conclusion is confirmed by the observation that acylation of the amino terminus of a pHD peptide with chains as short as 6 carbons promotes membrane binding at pH 7 and gives the peptide potent nanoporation activity at all pH values.⁸ While peptide amphipathicity is necessary for nanoporation, alone it is not sufficient to explain the propensity of the pHD peptides to insert across the membrane and form macromolecule-sized nanopores.

Critically important questions remain about the sequence–structure–function relationships of the pHD peptide nanopores. First, why are pHD peptides stable when inserted into membranes with membrane-spanning orientations? There are few, if any, peptide sequences, synthetic or natural, that insert across membranes while carrying so many charged and polar groups. Many other amphipathic helices have been studied, and nearly all of them, including the parent sequence melittin that was evolved into the pHD peptides, are interfacially bound with equilibrium orientations that are mostly parallel to the membrane surface. Second, which interactions stabilize the nanopores? Among those few amphipathic peptides that span membranes, most form compact structures that release small molecules preferentially.¹¹ Third, how does nanopore stability depend on the charged state of the ionizable sidechains?

Nanopores formed by pHD peptides are heterogeneous in size⁵ and likely dynamic, properties which preclude structure determination using experimental structural biology methods. Therefore, to address these questions, we have used MD simulations to probe the structure and dynamics of pHD108 peptides inserted into membranes. Based on the size distribution of the nanopores measured by atomic force microscopy (AFM),⁴ we anticipate a nanopore can form from as few as 8 peptides. Thus, here we studied 8-mers with distinct protonation states and initial configurations. Taken together, the ensemble of the simulations we report suggest that the long Glu and Gln sidechains selected during the evolution of the pHD peptides stabilize the pore via a collective network of direct and water-mediated H-bonds that includes many lipid molecules. This H-bond network enables the pHD peptides to overcome their high net hydrophilicity and insert across the membrane. This work outlines a path to circumvent the classical requirement of high net hydrophobicity for membrane insertion.

Results

Rationale for the choice of the protonation states of titratable sidechains

The synthetically evolved nanopore-forming peptide pHD108 contains seven ionizable residues, 5 glutamates (E) and 2 histidines (H), along with several polar glutamine (Q) residues. These ionizable and polar residues share ten helical i to i + 3 or i to i + 4 spacings all along the length of the helix from residue 4 to residue 25, as shown in Fig. 1A. This pattern places the polar residues along one continuous face of an α-helix, the secondary structure that they are known to attain in membranes.^4,5,7 The amphipathicity of pHD108 is shown by a helical wheel diagram, Fig. 1B and idealized helical structures in Fig. 1C & D. By studying the pH-dependence of acylated and non-acylated pHD peptides, we have shown that the ionizable glutamate residues are mostly in their charged states in nanopores comprised of pHD108,^5,7,8,12 In fact, these charged residues are likely critical for the nanopore-forming activity of the pHD peptides as they were strongly selected for during the original high throughput screening.⁴

Fig. 1 Primary and secondary structure of pHD108. (A) The amino acid sequence of pHD108. Acidic sidechains are red and basic sidechains are blue. Polar, uncharged glutamines are orange. The proline at position 14 is green and the tryptophan at position 19 is purple. Helical i-to-i + 3 and i-to-i + 4 spacings between polar groups are shown with black and red lines. B. Helical wheel diagram of pHD108 showing the amphipathic segregation of polar and nonpolar sidechains, assuming a completely helical structure. (C and D) pHD108 shown as an α-helix with the sidechains colored in the same way as in panel A. The helix in panel D is rotated 180° around the helical axis relative to panel C.

Initial conditions of the simulations

Here we used MD simulations to explore the stability and interactions in systems composed of 8 pHD108 peptides initially inserted symmetrically in a circular arrangement in the bilayer. We chose pHD108 as the representative pHD peptide as this peptide is the best studied. We rely on two sets of observations to develop the simulated system. First, AFM data on the pore size and shape indicates that stable, circular pores of 2–3 nm diameter can be formed.⁵ About eight pHD peptides arranged approximately symmetrically relative to a central axis are sufficient to create a pore of this observed size. Second, simulations of the very similar pHD15 peptide inserted as single monomeric peptides into phosphatidylcholine (PC) membranes show that the peptide samples transmembrane orientations that are tilted relative to the membrane normal, with a bend near the central proline residue so that the peptide termini are in contact with the two bilayer interfaces.⁷ In the simulations, we placed the peptides in hydrated bilayers made from 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) as pHD peptides are highly active in this bilayer.⁷ For the starting peptide configuration we used the orientation as sampled by the pHD15⁷ and created 8 copies of the pHD108 peptide which we arranged symmetrically without direct contact between peptides, without water in the bilayer, and with lipids placed in the inner region delineated by the peptides. This state mimics a bilayer with inserted peptides that are not yet hydrated or in contact with other peptides (SI Fig. S1). With this initial peptide placement, we performed simulations with distinct protonation states as discussed below. To explore the role of the insertion pathway on nanopore formation, we have additionally performed three independent simulations initiated with water in the area delineated by the peptides (SI Fig. S2). Here, initial lipids or initial water in the inner area delineated by the peptides are designated by a “/L” or “/W” at the end of the simulation name (Table 1).

Table 1 Twelve atomistic MD simulations, amounting to 7.8 µs, were carried out in this work. Simulation names and the protonation states of the glutamate and histidine residues are shown in the first two columns. Whether the area delineated by the eight peptide initially contained lipids or water is indicated in the third column. Main and repeat simulations are indicated in the fourth columns, and the simulation time for each individual simulation is indicated in the last column

Name	Protonation	Lipids or water initially in pore?	Main or repeat?	Simulation time (ns)
Sim 5-/L	E4, E8, E11, E15, E18 charged (−)	Lipids	Main	603
	H7, H21 neutral, Nd1		Repeat	474
Sim 3-/L	E4, E8, E11, E15, E18 charged (−)	Lipids	Main	622
	H7, H21 charged (+)		Repeat	621
Sim 1-/L	E4, E8, E15 charged (−)	Lipids	Main	851
	E11, E18 neutral		Repeat	526
	H7, H21 charged (+)
Sim 2+/L	E4, E8, E15, E11, E18 neutral	Lipids	Main	674
	H7, H21 charged (+)		Repeat	604
Sim 5-/L*	7 peptides as in Sim 5-/L	Lipids	Main	1300
	1 peptide all neutral
Sim 5-/W	E4, E8, E11, E15, E18 charged (−)	Water	Main	509
	H7, H21 neutral, Nd1
Sim 3-/W	E4, E8, E11, E15, E18 charged (−)	Water	Main	533
	H7, H21 charged (+)
Sim 2+/W	E4, E8, E11, E15, E18 neutral	Water	Main	506
	H7, H21 neutral, Nd1
			Total: 7.8 µs

---

Guided by experimental observations from pH-dependent pore formation,^4,5,7 we chose protonation states of selected titratable sidechains to explore key questions. Specifically, we simulated five scenarios distinguished by the protonation states of the Glu and His sidechains. The simulations are summarized in Table 1. First, we studied simulations in which all titratable sidechains have protonation state expected at neutral pH, i.e., all Glu sidechains are negatively charged, and both His sidechains are neutral and described as Nδ1 tautomers. Because this state has a net sidechain charge of minus 5, we call this group, Sim 5-/ denoted as Sim 5-/L when initiated with lipids in the area delineated by the peptides, and Sim 5-/W when initiated with water in this area. Second, we studied a scenario representing a possible pH 5 state in which all glutamates remain deprotonated, but the His residues are protonated, and thus charged. This is called Sim 3-/. Third, we studied an acidic pH state we have proposed previously that has three charged glutamates, two neutral glutamates and two charged histidines to test our hypothesis that E11 and E18 are important for pH-dependent membrane binding, while E4, E8, and E15 are important for the structure of the pore.¹² This group is called Sim 1-/. Fourth, we test our non-classical hypothesis that the charged residues are critical for stabilizing the membrane inserted nanopore state, by studying pHD peptides with all Glu sidechains uncharged and the His residues either uncharged or charged. We hypothesized that these simulations, called Sim 0-/ or Sim 2+/, will not result in stable nanopores. Fifth, we studied an 8-mer in which 7 peptides have all their Glu sidechains negatively charged as in Sim 5-/, but the 8^th peptide has all Glu sidechains neutral. This is a variant we call Sim 5-/*.

Nanopore structures initiated with lipids inside the pores

We simulated these eight scenarios for a total time of 7.8 µs. Most simulations were run for ∼500 ns or more, see Table 1 and the four most important simulations were repeated in independent trajectories from the same starting coordinates. Interactions and dynamics of peptide pores in membranes occur on timescales from nanoseconds, for sidechain motions, to seconds or longer, for large-scale collective structure changes. All-atom molecular dynamics (MD) simulations are performed on timescales from nanoseconds to about 1 microsecond, providing atomic level insights into the fundamental interactions and structural fluctuations within a peptide pore. MD simulations thus offer a mechanistic foundation for interpreting slower experimental phenomena.

As we show in SI Fig. S3 & S4, all simulations reach steady state within 50–100 ns. In the simulations starting with lipids inside the pore, peptides with a net negative charge contributed by the Glu sidechains, including Sim 5-/L, Sim 3-/L and Sim 1-/L, were observed to be compatible with nanopore formation, Fig. 2A–C. In all cases, these peptides remain inserted across the membrane, with the N-terminal part of the peptides farther apart and the C-terminal ends closer together, creating a cone-like pore shape. Extensive amounts of water enter the pore along the polar face of the peptides, especially on the more open N-terminal side, and some lipids leave the central pore area, see also SI Fig. S5. Multiple lipid headgroups are found unusually deep in the bilayer in these simulations, Fig. 2A–C. Sim 3-/L, especially, quickly evolved into a complete water-filled, membrane-spanning nanopore of 2 × 3 nm diameter.

Fig. 2 Illustration of the pore architecture and interactions sampled in the main simulations initiated with lipids in the area delineated by the peptides. Peptides are shown as green ribbons and lipid phosphate groups are shown as brown van der Waals spheres. Water molecules are shown as van der Waals spheres with oxygen atoms colored red and H atoms, white. (A–D) Coordinate snapshots from the end of simulations that each have distinct protonation states. Lipid–peptide interactions are shown on the left and water-peptide interactions are shown on the right side. The corresponding number density profiles are shown in Fig. 3.

In Sim 3-/L and Sim 1-/L the arrangement of the peptides remains roughly circular, while in Sim 5-/L, the initial circular arrangement of the peptides is distorted into oval or “U” shaped arrangements. In both cases the pores are stable on simulation timescales, with no evidence of collapse or peptide de-insertion from the membrane.

Many of these simulations were repeated over independent trajectories starting from the same initial conditions. These repeat simulations produced similar equilibrated structures, SI Fig. S6.

Peptides with neutral glutamates do not form nanopores

In the simulation with neutral glutamates, Sim 2+/L, nanopore formation does not occur. Almost no water, and only few lipid phosphate groups, enter the core of the bilayer, despite the continued presence of the 8 highly amphipathic pHD108 molecules with many polar sidechains in the membrane, Fig. 2D. The eight peptides collapse into a compact, but dry structure, suggesting that charges are needed for nanopore formation.

Sim 5-/L*, containing a single peptide with neutral glutamates, plus 7 charged peptides, also formed an apparently stable pore under these conditions. However, the single uncharged peptide has a final orientation that is nearly parallel to the membrane surface indicating that it is not stable in the membrane-spanning pore state, and is de-inserting into an interfacially bound state.

Water density in the simulated membranes

To explore the pore, we determined normalized number densities of peptides (green), water (blue), and lipid phosphate groups along the bilayer normal, Fig. 3. Lipid phosphates are divided into those that are within 15 Å of a peptide (orange), and those that are farther away (black). In a bilayer without peptides, both water and lipid phosphates are known to be essentially excluded from the bilayer interior.¹³ The pore-compatible simulations demonstrate dramatic water penetration into the membrane core region that extends from −10 to +10 Å from the bilayer midplane. This is especially true at the N-terminal opening of the pore structure, Fig. 3A–C. These simulations also show significant penetration of lipid phosphates into the core of the bilayer. In the 2+/L structure that has uncharged glutamates, and does not form nanopores, neither water nor lipid phosphates are present in the bilayer core, Fig. 3D. Peptide number densities show peptides all across the core of the bilayer. The most highly charged pore structures, Sim 5-/L (Fig. 3A) and Sim 3-/L (Fig. 3B) show asymmetry in the peptide distributions and in the lipid distributions, resulting from the conical pore shape with the peptide N-termini farther apart than the C-termini.

Fig. 3 Number density calculations for simulations initiated with lipids initially in the area delineated by the peptides. For each simulation we show the number densities for peptides (green), for lipid phosphate groups within 15 Å of the peptides (orange), and for all other lipid phosphate atoms (gray). Separately in the right column we show the water distribution (blue). The blue rectangles show the hydrocarbon core of the lipid bilayers. For clarity, the water number density is shown only for ±35 Å along the membrane normal. (A–D) Density profiles for the four main simulations indicated. Qualitatively similar profiles were obtained from the corresponding repeat simulations.

Water-bridged H-bond network

To evaluate efficiently the dynamics of peptide-water H-bond networks that could stabilize the periphery of the pores, we used the graph-based algorithm Bridge^14,15 to dissect the dynamic H-bond networks of each system studied. For each simulation performed, we examined the molecular interactions of the peptides, water, and lipids using H-bond graphs computed separately for direct side chain H-bonds, for one-water bridges between sidechains, and for bridges with up to three H-bonded waters which can bridge distances of about 10–12 Å.¹⁶ We determined the average % occupancy of each H-bond to derive a coarse estimate of the stability of the bridges. Most H-bond graphs are created with a cutoff of 30% occupancy as we have done previously.¹² For the graphs computed only for direct H-bonds between sidechains, which are uncommon, we found it necessary to lower the minimum H-bond occupancy to 10% in order to identify even transient direct H-bonds between sidechains. To illustrate the participation of water molecules to the H-bond networks, we determined the average number of waters per each bridge as shown in SI Fig. S7.

In Fig. 4 we show an H-bond analysis for Sim 3-/L as an example. H-bond analyses for Sim 5-/L, Sim 2+/L and Sim 5-/L* are shown in SI Fig. S8–11. Each H-bond in Fig. 4 is shown with its sampling frequency, or % occupancy. The 16 direct H-bonds identified in Sim 3-/L are shown in Fig. 4A. Direct H-bonds between adjacent peptides are uncommon, with only 3 of 16 in this simulation. Most direct H-bonds occur between sidechains in the same peptide, either between adjacent residues or between residues separated by one helical turn. The 21 H-bonds bridged by one water are shown in Fig. 4B. Here, also, H-bonds between adjacent peptide are rare with all but one of these one-water H-bonds occurring between residues in the same peptide. All of these intrapeptide, one-water interactions occur between residues separated by one helical turn of 3 or 4 residues. All H-bonds with three or fewer waters bridging them are shown in Fig. 4C. At this level, the extensive H-bond network across the membranes, supporting the pore structure, is apparent.

Fig. 4 H-bond graphs computed from Sim 3-/L. The red asterisks in panels A and C indicate sites where peptides inter-connect. Values on the connecting lines are % occupancy for that interaction over the last 200 ns of the simulation. (A) Graph of the direct H-bonds between peptide sidechains. (B) Graph of the one-water-mediated bridges between peptide sidechains. (C) Graph of the direct and water-mediated H-bonds with three or fewer waters. The minimum H-bond occupancy shown is 10% in panel A, and 30% in panels B–C. (D) View, from the membrane surface, of the eight peptides color coded to match the nodes in panels A, B and C. Glu, Gln, His, and Trp sidechains of peptide G are shown as van der Waals spheres. (E—G) Close view of individual peptides with labels for selected sidechains.

Interesting patterns are revealed when we map the specific H-bonds from Fig. 4C onto the sequence of pHD108, Fig. 5. Here, using Sim 3-/L as an example, it is clear that all polar groups have H-bonds to other polar groups, and that most polar-polar helical spacings are represented in the map. By far the most abundant H-bonds in the network utilize the glutamate residues, and many are water bridged glutamate-glutamate interactions. As we have discussed above, H-bonds between different peptides are much less common than H-bonds within each peptide. While several intrapeptide H-bonds are found in all eight peptides in Sim 3-/L, no individual inter-peptide H-bond is represented more than twice.

Fig. 5 H-bond network from Sim 3-/L. (A) The intra peptide H-bonds with 3 or fewer waters and at least 30% occupancy within the eight peptides of the pore are shown as solid lines. The thickness of the line indicates how many peptides had the interaction indicated (see scale). (B) The interpeptide, or non-self H-bond interactions drawn symmetrically. The thickness of the lines indicates how many peptides had the interaction indicated, however nearly all of these were observed in only one peptide. (C) The H-bonds to lipids observed in Sim 3-/L. The thickness of the lines indicates how many peptides had the interaction indicated.

We have also mapped the number of H-bonds to lipids, which are shown in Fig. 5C and discussed in detail below. Here again, glutamate residues are the most important contributors to the H-bond network. Both intramolecular H-bonds within the peptides, Fig. 5A, and H-bonds to lipids, Fig. 5C, are more abundant on the N-terminal half of the peptide where the peptides are farther apart from each other. H-bonds between peptides are more common in the C-terminal half of the sequence, where the peptides are close together.

The H-bond network depends on protonation state

In the left column of Fig. 6 we show H-bond networks for Sim 5-/L, Sim 3-/L and Sim 2+/L with the lipids that also participate in the networks. In the right column, we show the numbers of intrapeptide, interpeptide and peptide-lipid H-bonds for the three simulations. In the pore-compatible simulations, Sim 5-/L and Sim 3-/L, we observe a rich H-bond network involving all the polar and charged residues along with many lipids. In these pore-compatible simulations, more than 140 direct and water-bridged H-bonds that have at least 30% occupancy are observed. The most abundant H-bonds are the intramolecular bonds within each peptide (blue hatch). The five glutamate residues have the greatest number of these intramolecular H-bonds. H-bonds between different peptides are shown in red. These intermolecular H-bonds are less abundant, but constitute a critical crosslinking element of the network. In the Sim 5-/L and Sim 3-/L octamers, there are at least 60 direct and water bridged H-bonds between peptides and lipids. These are shown in black. H-bonds to lipids that bridge two different peptides are shown in black hatching. We note that lipid-bridged interactions between different peptides (black hatching) are similar in abundance to direct and water-bridged interactions between different peptides (red) and are thus similarly important contributors to the stabilizing H-bond network.

Fig. 6 Lipid-water-peptide H-bond networks in pHD108. We computed H-bond graphs for peptide sidechains, lipid phosphate groups, and up to three H-bonded water bridges; the minimum H-bond occupancy shown is 30%. (A) H-bond graph of Sim 5-/L. Nodes representing lipid phosphate groups are shown as pink circles, which are filled when between two distinct peptides; lipid nodes and edges that lack connections to peptide nodes, and only indicate lipid-lipid interactions, are excluded. Graph nodes colored other than pink represent peptide sidechains. Lines (edges) between two graph nodes represent H-bonds, which can be direct or mediated by bridges of up to three H-bonding water molecules. (B) Density of H-bond interactions for each sidechain in Sim 5-/L, presented as the sum of interactions per residue, for all eight peptides of the pore. (C) H-bond graph for Sim 3-/L. (D) H-bond density for Sim 3-/L. (E) H-bond graph for Sim 2+/L. (F) H-bond density for Sim 2+/L.

The H-bond network is significantly diminished in the 2+/L simulation with glutamates uncharged, Fig. 6E & F. This is despite the fact that protonation of glutamates creates H-bond donors in the membrane that are not present when glutamates are deprotonated and charged. Further, the protonation of glutamates will greatly diminish the polarity of the sidechain¹⁰ which would normally be expected to increase the stability of a membrane-spanning structure. Yet, we observe a significant decrease in the number of H-bond interactions. Intrapeptide interactions are especially diminished in this simulation. Peptide-lipid interactions are similar in number but occur mostly in the bilayer interface in Sim 2+/L, and there are no lipid phosphates deep in the bilayer.

In Fig. 7, we show a schematic diagram of the entire interaction network identified in Sim 3-/L. Only the polar and charged sidechains of each peptide are shown, and they are arranged in columns. Interactions that do not bridge two peptides are shown as thin lines and interactions that bridge two peptides are shown as thick lines. The many interactions with lipids that participate in the network are also shown, and those lipids that participate in crosslinks between two peptides are outlined in black. More than 150 interactions are sampled in the octameric pore simulations, and ∼40 interactions to lipid molecules were also sampled.

Fig. 7 Schematic illustration of the entire interaction network observed in Sim 3-/L. All interactions shown were bridged by three or fewer waters and had occupancies of at least 30%. Thin lines represent interactions that do not connect two different peptides. Thick lines show interactions that connect two peptides. Lipids that bridge two peptide are outlined in black.

The pore structure includes lipids with unusual conformations

One of the dominant features of the pHD108 nanopores in simulations 5-/L, 3-/L and 1-/L is the extensive involvement of lipids in the pore structure. These structures have many lipids in the pore, with phosphates deep in the membrane where they are not normally found. These lipid phosphates contribute extensively to the H-bond network, Fig. 7. This kind of lipid-rich pore architecture has long been called a “toroidal pore”,¹⁷ although actual evidence for this type of structure has previously been scarce. By significantly involving lipids in the pore structure, a toroidal pore contrasts with a “barrel stave” pore in which the pore is stabilized by extensive lateral interactions between close-packed peptides. To demonstrate the unusual conformations of the lipids in the pore, we show, in Fig. 8, example lipid molecules that are near the peptides in Sim 5-/L. For visualization, we have separated them loosely into three broad, non-exclusive conformations that we call “splayed”, for lipids with headgroups in the pore that encompass a peptide between their acyl chains; “between” for lipids with headgroups that are in the aqueous pore between two peptides, and “perpendicular” for lipids with any chains that are orientated roughly perpendicular to lipids in an unperturbed bilayer. For reference we also show lipids from a “normal” unperturbed part of the bilayer in the same simulation, but away from the pHD peptides. As we discuss below, these many acyl chains with unusual orientations are likely to be an essential mechanism for stabilizing the boundary between the pore and the lipid bilayer. Repeat simulations starting from the same initial structures produced very similar equilibrated structures, SI Fig. S6, and H-bond networks, SI Fig. S12–S14.

Fig. 8 Lipids with unusual conformations in the pore. Endpoint snapshots of Sim 3-/L are shown here to demonstrate the variety of lipids in contact with the peptides which have conformations and orientations that are not frequently observed in bilayers away from the pore. (A–H) Side and top views of four lipid classes. These non-exclusive structural categories are defined as follows: Splayed lipids have a peptide in the area between the two acyl chains. Sideways lipids have acyl chains that are roughly perpendicular to the plane of the bilayer. Between lipids have headgroups in the inner region of the pore, between two peptides. Normal lipids are taken from a portion of the Sim 3-/L bilayer away from the area of the pore.

Nanopore structures with initial water-filled pores

Three simulations, in which the area delineated by the membrane inserted peptides was initially filled with water instead of lipids, are denoted as Sim 5-/W, Sim 3-/W and Sim 0-/W, Table 1. This starting scenario enables us to further explore potential pathways to nanopore formation, which are unknown, as well as the role of lipids and water in nanopore stability.

The initial structure for simulations with water initially in the area delineated by the peptides is shown in SI Fig. S2, and end point structures are shown in Fig. 9. Time series for these simulations are shown in SI Fig. S15. Other structural views and H-bond numbers for the initial water-filled pore are shown in SI Fig. S16–S22. The pores remain filled with water during the course of the simulations, but all three structures show evidence of instability. For example, in Sim 5-/W and Sim 3-/W several of the peptides de-insert almost completely from the bilayer and have interfacial orientations by the end of the simulation. For Sim 5-/W, especially, all eight peptides move towards interfacial locations on the N-terminal side of the pore, leaving the C-terminal side empty of peptides, and filled by lipids. Sim 3-/W and Sim 0-/W both have pore structures that are much more compact than the starting configuration, suggesting that they are collapsing into a dense oligomer in the membrane, rather than an open nanopore.

Fig. 9 Illustration of the pore interactions sampled in the simulations that were initiated with water in area delineated by peptides. We used coordinate snapshots from the end of each simulation. Peptides are shown as green ribbons and lipid phosphate groups are shown as brown van der Waals spheres. Water molecules are shown as van der Waals spheres with oxygen atoms red and H atoms white. (A–C) Images from the three main simulations indicated. The corresponding number density profiles are shown in Fig. 10.

Number densities along the bilayer normal, for the simulations that were initiated with water in the pore, are shown in Fig. 10. De-insertion of peptides in the Sim 3-/W pore can be seen in the peptide peak which has shifted to +12 Å from the bilayer center, compared to a peak at +10 Å for Sim 3-/L. For Sim 5-/W, as noted above, the peptide density is low on the C-terminal side of the bilayer, in the region of −10 to 0 Å. As the peptides are moving toward de-insertion on the N-terminal side lipid phosphates appear to be filling in the space. In the Sim 0-/W structure, the collapse of the peptides into a compact non-pore state is indicated by a peptide peak number density at +2–8 Å and lack of water penetration. It is likely that the initial water-filled nanopore structures are not stable when initiated with water in the pore.

Fig. 10 Number density calculations for simulations initiated with water initially in the area delineated by the peptides. For each simulation we show the number densities for peptides (green), lipid phosphate atoms within 15 Å of the peptides (orange) and all other lipid phosphate atoms (gray) and, separately for water (blue). The blue rectangles show the hydrocarbon core of the bilayer. For clarity, the water number density is shown only for ±35 Å along the membrane normal. (A–C) Number density profiles computed for the three simulations in this group.

In SI Fig. S23 we show the radial distributions of water molecules and lipid phosphates around the peptides in the various simulations. The radial distribution of water around the peptides has distinct peaks at about one, two, and three water shell distances from the peptides. Qualitatively, more waters than lipids are found in simulations that have stable pores, panels A, B, D, E, G and H. In simulations without stable pores, panels J & K, more lipids than water are found within close distance due to the more compact structures and lack of bulk water in the pore. Simulations that were initiated with water in the region delineated by peptide, panels C, F and I, have intermediate ratios of waters to lipids.

Discussion

pHD peptides defy classical hydrophobicity

By classical hydropathy analysis, and by comparison to known membrane-spanning helices, the pHD peptides appear to be far too polar to be stable in a membrane-spanning state, especially when the glutamates are charged. For example, pHD108 with charged glutamates and neutral histidines has an unfavorable Wimley-White octanol scale hydrophobicity¹⁰ of +14.6 kcal mol⁻¹ (4.18 J cal⁻¹). This is very unfavorable compared to the values for known TM helices. Even when all ionizable amino acids are neutral, pHD108 has an unfavorable hydrophobicity of +1.5 kcal mol⁻¹. In contrast, the well-studied TM helices of human glycophorin A and the human EGF receptor have favorable values of −9.5 kcal mol⁻¹ and −9.6 kcal mol⁻¹, respectively. Even the notably polar S4 sensor helix of voltage gated potassium channels,^18,19 with four arginine residues, has a favorable Wimley-White hydrophobicity of −8.9 kcal mol⁻¹. It is thus remarkable that the pHD peptides exist in membrane-spanning configurations^4,5,7 in the form of large nanopores, despite the high net polarity of the pHD108 sequence. These findings defy classical physical principles of hydrophobicity-driven membrane insertion.

Mechanism of pH dependence

At pH 7, the pHD peptides are poised at the brink of being hydrophobic enough to bind to membranes and fold into α-helical nanopores.^4,5 Experimentally, a slight shift in that equilibrium is sufficient to drive membrane binding and nanopore formation. This shift can be accomplished by a slight acidification, which may protonate one or two glutamate residues^5,7 or by addition of acyl groups as short as 6 carbons.⁸ This idea is supported by our observations made with the macrolittins, peptides that were selected for potent, pH-independent nanoporation at pH 7 using the same screen and library that was used for the pHD peptides. The macrolittins are nearly identical to the pHD peptides, but they have three glutamate residues, compared to five for the pHD peptides. This strongly supports the idea that the protonation of one or two glutamates can drive membrane binding and pore formation in the pHD peptides.

Here we selected various protonation states to test specific ideas about nanopore stability. All protonation states that had at least some of the glutamates in the charged state formed stable nanopore-consistent structures. Interestingly the most realistic protonation state for pH 5–6, Sim 3-/L, with both glutamates and histidines charged had the most ideal nanopore in the simulation. Sim 1-/L, which mimics a macrolittin and a partially protonated pHD peptide by having three charged glutamates and two protonated glutamates, is shown in these MD simulations to form a stable pore with a rich H-bond network. On the other hand, simulations with the glutamates uncharged did not result in stable inserted nanopore structures, strengthening our paradoxical conclusion that charged glutamates are needed for maximum pore stability.

H-bonds stabilize the nanopore

In a pure bilayer, a very large free energy penalty is associated with forming an aqueous pore,^2,3 due to the loss of cohesive interactions between lipids and due to the necessary exposure of the bilayer hydrocarbon core to water. The cost of pore formation scales with the perimeter of the pore, so that larger pores are more costly. Consequently, peptides that form large nm-scale pores are very rare. To form such a pore, the large free energy penalty associated with the pore perimeter must be reduced by the pHD108 helices. Here, we used MD simulations to investigate how this is acheived.

By definition, a pore-forming peptide must interact with the bilayer-water interface on the pore perimeter to eliminate the free energy penalty associated with exposure of lipid hydrocarbon to water. A peptide that forms a large, macromolecule-sized nanopore must stabilize a large pore better than a small pore. The atomistic MD simulations described here have begun to explain the basis for the “non classical” stability of the pHD peptide nanopores. Specifically, we show that the pHD peptide nanopores are stabilized by a collective membrane-spanning, water-mediated H-bond network involving anionic, basic, and polar sidechains, as well as many lipid head groups. Direct H-bonds are rare, but water-bridged H-bonds within peptides, between peptides, and to lipid headgroups are abundant. The abundance of i–i + 3 and i–i + 4 helical spacings between ionizable and polar sidechains, which were selected in the high throughput screen used to identify the pHD peptides, allows for many H-bonding interactions to be sampled simultaneously. In fact, more than 150 interactions that have occupancies over 30% are identified for each octameric nanopore.

Lateral interactions between different peptides, which create the network, include a few direct H-bonds, multiple water-bridged H-bonds between peptides, and multiple lipid-bridged interactions between peptides. Taken together, these many intra- and inter-peptide interactions create a dense collective interaction network that spans the bilayer.

An interesting observation in this work is that the most abundant participants in the collective H-bond network are deprotonated glutamate sidechains and lipid phosphates. These moieties have no H-bond donor hydrogens and cannot H-bond directly to each other. Instead, the network is comprised of these many H-bond accepting oxygens on peptides and lipids, bridged to each other by water molecules. When the glutamate sidechains are protonated, placing many new H-bond donor groups in the bilayer, the nanopore is less stable, not more stable.

Contributions of the lipids to nanopore formation

The simulations suggest that the many H-bonds to lipids are a critical feature of stable peptide nanopores. For example, simulations with neutral glutamates, such as 2+/L, that do not show stable pores, have almost no lipid phosphates deep in the bilayer, and have few stabilizing H-bonds to lipid polar groups. This is despite the fact that protonated glutamate sidechains have H-bond donors that one might expect to improve H-bonding in the pore. We conclude that anionic sidechains interacting with lipid headgroups through water-bridged H-bonds are a critical stabilizing feature of pHD peptide nanopores. In fact, the interaction network, overall, contains more lipids than peptides. Further, many of the lipids that participate in the interaction network have unusual conformations in the pore, including many with acyl chains that interact with the pore periphery. These lipids are likely essential for stabilizing the pore periphery by helping to complete the interface between the bilayer hydrocarbon and the aqueous milieu of the pore. This observation is consistent with what many researchers have previously called “toroidal” pores. The simulations initiated with water in the pore provide additional strong evidence for the critical role of the lipids in pore stabilization, as these pores are unstable, with peptides de-inserting or the octameric pores collapsing into more compact structures suggesting that the insertion pathway may begin with peptide insertion, followed by water entry into the membrane along the carboxyl and other polar groups.

In conclusion, the surprising stability of macromolecule-sized nanopores formed by the pHD peptides is the result of a dense membrane-spanning H-bond network comprised of peptide polar and charged groups, lipids, and waters. Ongoing simulations are being used to explore the processes of insertion and nanopore self-assembly.

Methods

Starting coordinates for pHD108 simulation systems

To generate starting coordinates for pHD108 we started from a coordinate snapshot of pHD15 from a previous simulation.⁷ Using Chemistry at Harvard Molecular Mechanics (CHARMM), we mutated the sidechains of pHD15 groups D11, D12, and D15, into the corresponding E11, G12 and E15 of pHD108. We multiplied the coordinates of the peptide and used rigid-body rotations and translations to generate the coordinates of a system composed of eight peptides symmetrically distributed relative to the central axis, with the Glu sidechains towards the inner region delineated by the peptides (SI Fig. S1).

Protonation of titratable sidechains

The pHD108 peptide contains 4 Glu and 2 His amino acid residues (Fig. 1A). To explore how protonation of these six titratable amino acid residues impacts the dynamics of pHD108 systems, we performed five independent simulations distinguished by the protonation of the Glu and His sidechains, as described below and summarized in Table 1.

In Sim 5-/L, we consider all eight peptides with standard protonations: all Glu sidechains are negatively charged, and all His sidechains are neutral and described as Nδ1 tautomers. We anticipate that a pore might be observed during the simulations. Sim 3-/L, with all Glu sidechains negatively charged, and with all His sidechains protonated, approximates pH ≤ 6.5, and could be compatible with pore formation. Sim 1-/L has all peptides with negatively charged E4, E8, and E15, neutral E11 and E18, and positively charged H7 and H21; this protonation is used to test our hypothesis that E11 and E18 are important for proton binding, and E4 and E8, for the structure of the pore. Sim 2+/L has all Glu sidechains neutral, and all His sidechains positively charged, approximating protonation at pH < 3. Sim 5-/L* has 7 peptides negatively charged, with the same protonation as Sim 5-/L, and the 8^th peptide with all Glu sidechains neutral. We use this simulation to test our hypothesis that a neutral peptide is incompatible with pore formation. We list all simulations in Table 1, where we also show the lengths of all main and repeat simulations.

Protocol for MD simulations of the 8-pHD108 peptide systems

We used CHARMM-GUI^20,21 to place the 8 pHD108 peptides in hydrated POPC lipid membranes with 0.1M NaCl neutralizing salt concentration. Each simulation systems contains ∼532 lipid molecules and ∼26 000 waters, for a total number of ∼152 660–152 890 atoms. We used CHARMM36m force-field parameters^22–26 to describe protein groups, lipid molecules, and ions, and the TIP3P water model.²⁷

MD simulations were performed with NAMD.^28,29 Following heating and brief equilibration with weak harmonic constraints as set in the standard CHARMM-GUI equilibration protocol, all harmonic constraints were switched off and production runs were performed in the NPT ensemble (constant number of atoms N, constant pressure P = 1 bar, constant temperature T = 303.15 K) using a Langevin dynamics scheme with a Nosé–Hoover piston.^30,31 Covalent bonds to H atoms were fixed.³² We used an integration step of 1 fs during equilibration and the first 1 ns of the production run; for the remaining of the production runs we used a multiple timestep integration scheme^33,34 with 1 fs for the bonded forces, 2 fs for short-range non-bonded forces, and 4 fs for long-range electrostatics. We used a switch function between 10 and 12 Å for the real space interactions, and a smooth particle mesh Ewald summation for Coulomb interactions.^35,36 We saved coordinates each 10 ps. Trajectories are available upon request.

H-bond criteria and graphs of dynamic protein-water H-bond networks

We consider that two groups are H-bonded when the distance between the donor and acceptor heavy atoms is ≤3.5 Å, and the H-bond angle is ≤60°. The occupancy of an H-bond is given by the percentage of coordinate sets in which the H-bond is present according to the criteria we chose.

A water-mediated bridge between two protein groups consists of one or more H-bonded water molecules that inter-connect the H-bond donor and acceptor protein groups. Here we performed separate H-bond graph computations for one-water bridges, and for bridges of maximum 3 waters between protein sidechains. All H-bond computations were performed with the graph-based algorithm Bridge2.^14,15

Lipid, water, and peptide H-bond networks

For each simulation system we used the last 200 ns to perform separate H-bond graph computations for (i) direct H-bonds between sidechains; (ii) short one-water bridges between protein sidechains; (iii) three-water bridges between protein sidechains, and direct sidechain-sidechain H-bonds; (iv) lipid phosphate groups, water bridges of up to three H-bonded waters, and direct sidechain H-bonds. For simplicity, as graph computations with lipid headgroups give extended H-bond networks, we considered only lipids that in the last segment of the simulation trajectory have their phosphate atoms within 15 Å of the peptides. To monitor the number of H-bonds of the pore along the entire simulation, for each of the simulations performed we used the entire production run to perform an additional graph computation of all direct sidechain-sidechain H-bonds and three-water bridges.

Number density profiles

To evaluate the hydration of the pores we computed, separately, the normalized number density profiles for the peptides, water molecules, lipid phosphate atoms within 15 Å of the peptides, and all other lipid phosphate atoms. Average values, H-bond graphs, and hydrophobic contact graphs were computed from the last ∼200 (20 000 equally spaced coordinate snapshots) of each simulation. In time series, the origin of time is the start of the production run.

Conflicts of interest

The authors have no conflicts of interest to declare.

Data availability

The authors will make all data presented in this manuscript available upon request. The Supplementary Information contains additional data and analyses from the simulations reported here, including information from repeat simulations.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5nr03276h.

Acknowledgements

Supported by NIGMS 1R01GM151326 to W.W. A-NB gratefully acknowledges computing time on the supercomputer JURECA³⁷ at Forschungszentrum Jülich under grant no. PHDPORES, and Dr Emiliano Ippoliti for excellent technical support. An initial set of simulations was performed on the North-German Supercomputing Cluster, HLRN. S. L. acknowledges support from the Slovenian Research and Innovation Agency (ARIS) program and project grants P2-0438, J1-50034, and J1-60001.

References

1. S. Guha, J. Ghimire, E. Wu and W. C. Wimley, Mechanistic Landscape of Membrane-Permeabilizing Peptides, Chem. Rev., 2019, 119, 6040–6085.
2. M. T. Lee, F. Y. Chen and H. W. Huang, Energetics of pore formation induced by membrane active peptides, Biochemistry, 2004, 43(12), 3590–3599.
3. P. F. P. Almeida and A. P. Almeida, Interactions of Antimicrobial Peptides with Lipid Bilayers, Comp. Biophys., 2012, 5, 189–222.
4. G. Wiedman, S. Y. Kim, E. Zapata-Mercado, W. C. Wimley and K. Hristova, PH-Triggered, Macromolecule-Sized Poration of Lipid Bilayers by Synthetically Evolved Peptides, J. Am. Chem. Soc., 2016, 139, 937–945.
5. S. Y. Kim, A. E. Pittman, E. Zapata-Mercado, G. M. King, W. C. Wimley and K. Hristova, Mechanism of Action of Peptides That Cause the pH-Triggered Macromolecular Poration of Lipid Bilayers, J. Am. Chem. Soc., 2019, 141(16), 6706–6718.
6. S. Li, S. Y. Kim, A. E. Pittman, G. M. King, W. C. Wimley and K. Hristova, Potent Macromolecule-Sized Poration of Lipid Bilayers by the Macrolittins, A Synthetically Evolved Family of Pore-Forming Peptides, J. Am. Chem. Soc., 2018, 140(20), 6441–6447.
7. S. Y. Kim, A. N. Bondar, W. C. Wimley and K. Hristova, pH-triggered pore-forming peptides with strong composition-dependent membrane selectivity, Biophys. J., 2021, 120(4), 618–630.
8. E. Wu, R. M. Jenschke, K. Hristova and W. C. Wimley, Rational Modulation of pH-Triggered Macromolecular Poration by Peptide Acylation and Dimerization, J. Phys. Chem. B, 2020, 124, 8835–8843.
9. E. Wu, A. Ellis, K. Bell, D. L. Moss, S. J. Landry, K. Hristova and W. C. Wimley, pH-Responsive Peptide Nanoparticles Deliver Macromolecules to Cells via Endosomal Membrane Nanoporation, ACS Nano, 2024, 18(50), 33922–33936.
10. S. H. White and W. C. Wimley, Membrane protein folding and stability: physical principles, Annu. Rev. Biophys. Biomol. Struct., 1999, 28, 319–365.
11. G. Wiedman, T. Fuselier, J. He, P. C. Searson, K. Hristova and W. C. Wimley, Highly efficient macromolecule-sized poration of lipid bilayers by a synthetically evolved peptide, J. Am. Chem. Soc., 2014, 136(12), 4724–4731.
12. L. Sun, K. Hristova, A.-N. Bondar and W. C. Wimley, Structural Determinants of Peptide Nanopore Formation, ACS Nano, 2024, 18, 15831–15854.
13. M. C. Wiener and S. H. White, Structure of a fluid dioleoylphosphatidylcholine bilayer determined by joint refinement of X-ray and neutron diffraction data. III. Complete structure, Biophys. J., 1992, 61, 434–447.
14. M. Siemers, M. Lazaratos, K. Karathanou, F. Guerra, L. S. Brown and A.-N. Bondar, Bridge: A graph-based algorithm to analyze dynamic H-bond networks in membrane proteins, J. Chem. Theory Comput., 2019, 15, 6781–6798.
15. M. Siemers and A.-N. Bondar, Interactive interface for graph-based analyses of dynamic H-bond networks: application to spike protein S, J. Chem. Inf. Model., 2021, 61, 2998–3014.
16. A.-N. Bondar, Mechanisms of long-distance allosteric couplings in proton-binding membrane transporters, Adv. Protein Chem. Struct. Biol., 2022, 128, 199–239.
17. L. Yang, T. A. Harroun, T. M. Weiss, L. Ding and H. W. Huang, Barrel-stave model or toroidal model? A case study on melittin pores, Biophys. J., 2001, 81(3), 1475–1485.
18. J. A. Freites, D. J. Tobias and S. H. White, A voltage-sensor water pore, Biophys. J., 2006, 91(11), L90–L92.
19. T. Hessa, S. H. White and H. G. von, Membrane insertion of a potassium-channel voltage sensor, Science, 2005, 307(5714), 1427.
20. S. Jo, T. Kim, V. G. Iyer and W. Im, CHARMM-GUI: a web-based graphical user interface for CHARMM, J. Comput. Chem., 2008, 29, 1859–1865.
21. E. L. Wu, X. Cheng, S. Jo, H. Rui, K. C. Song, E. M. Dávila-Contreras, Y. Qi, J. Lee, V. Monje-Galvan and R. M. Venable, et al., CHARMM-GUI Membrane Builder toward realistic biological membrane simulations, J. Comput. Chem., 2014, 35, 1997–2004.
22. B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus, CHARMM: a program for macromolecular energy, minimization, and dynamics calculations, J. Comput. Chem., 1983, 4, 187–217.
23. A. D. MacKerell Jr., M. Feig and C. L. I. Brooks, Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations, J. Comput. Chem., 2004, 25, 1400–1415.
24. B. R. Brooks, C. L. I. Brooks, A. D. MacKerell Jr., L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels and S. Boresch, et al., CHARMM: the biomolecular simulation program, J. Comput. Chem., 2009, 30, 1545–1614.
25. S. E. Feller and A. D. MacKerell Jr., An improved empirical potential energy function for molecular simulations of phospholipids, J. Phys. Chem. B, 2000, 104, 7510–7515.
26. J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Votrobyov, A. D. MacKerell Jr. and R. W. Pastor, Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types, J. Phys. Chem. B, 2010, 114, 7830–7843.
27. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926–935.
28. L. Kalé, R. Skeel, M. Bhandarkar, R. Brunner, A. Gursoy, N. Krawetz, J. Phillips, A. Shinozaki, K. Varadarajan and K. Schulten, NAMD2: greater scalability for parallel molecular dynamics, J. Comput. Phys., 1999, 151, 283–312.
29. J. C. Phillips, B. Braun, W. Wang, J. Gumbart, E. Takjkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale and K. Schulten, Scalable molecular dynamics with NAMD, J. Comput. Chem., 2005, 26, 1781–1802.
30. S. E. Feller, Y. Zhang, R. W. Pastor and B. Brooks, Constant pressure molecular dynamics simulation: The Langevin piston method, J. Chem. Phys., 1995, 103, 4613–4621.
31. G. J. Martyna, D. J. Tobias and M. L. Klein, Constant-pressure molecular-dynamics algorithms, J. Chem. Phys., 1994, 101, 4177–4189.
32. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, Numerical integration of the Cartesian equations of motion of a system with constraints. Molecular dynamics of n-alkanes, J. Comput. Phys., 1977, 23, 327–341.
33. H. Grubmüller, H. Heller, A. Windemuth and K. Schulten, Generalized Verlet algorithm for efficient molecular dynamics simulations with long-range interactions, Mol. Simul., 1991, 6, 121–142.
34. M. Tuckermann, B. J. Berne and G. J. Martyna, Reversible multiple time scale molecular dynamics, J. Chem. Phys., 1992, 97, 1990–2001.
35. T. Darden, D. York and L. Pedersen, Particle mesh Ewald: an N x log(N) method for Ewald sums in large systems, J. Chem. Phys., 1993, 98, 10089–10092.
36. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103, 8577–8593.
37. P. Thörnig, Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre, J. Large-Scale Res. Facil., 2021, 7, A182,  DOI:10.17815/jlsrf-7-182.

---