A New Mixed All-Atom/Coarse-Grained Model: Application to Melittin Aggregation in Aqueous Solution.

We introduce a new mixed resolution, all-atom/coarse-grained approach (AACG), for modeling peptides in aqueous solution and apply it to characterizing the aggregation of melittin. All of the atoms in peptidic components are represented, while a single site is used for each water molecule. With the full flexibility of the peptide retained, our AACG method achieves speedups by a factor of 3-4 for CPU time reduction and another factor of roughly 7 for diffusion. An Ewald treatment permits the inclusion of long-range electrostatic interactions. These characteristics fit well with the requirements for studying peptide association and aggregation, where the system sizes and time scales require considerable computational resources with all-atom models. In particular, AACG is well suited for biologics since changes in peptide shape and long-range electrostatics may play an important role. The application of AACG to melittin, a 26-residue peptide with a well-known propensity to aggregate in solution, serves as an initial demonstration of this technology for studying peptide aggregation. We observed the formation of melittin aggregates during our simulations and characterized the time-evolution of aggregate size distribution, buried surface areas, and residue contacts. Key interactions including π-cation and π-stacking involving TRP19 were also examined. Our AACG simulations demonstrated a clear salt effect and a moderate temperature effect on aggregation and support the molten globule model of melittin aggregates. As a showcase, this work illustrates the useful role for AACG in investigations of peptide aggregation and its potential to guide formulation and design of biologics.

Introduction

Peptide association is key for the structural organization within cells and for the proper functioning of biological processes[1−7] with associating partners ranging from oligo peptides, globular proteins, antibodies, and membrane proteins such as GPCRs. Misassociation of proteins, particularly proteinaceous aggregates,[8] can cause diseases such as mad cow disease (bovine spongiform encephalopathy)[9] and Alzheimer’s disease.[10,11] Control or disruption of peptide association is exploited by parasites,[12−14] venoms,[15−17] viruses,[18−20] and drugs.[21−27] Peptide–peptide interactions also need to be considered when formulating biologics, particularly with the recent trend of producing injectable biologics formulations at high concentrations which can suffer from high viscosity, aggregation, and precipitation.[28−31] The underlying problem may be that at high concentrations conventionally folded proteins can stick together or, more problematically, peptides can unfold leading to aggregation. For the latter, the aggregation has a greater propensity to be irreversible.

A range of experimental methods including size exclusion chromatography,[32] analytical ultracentrifugation,[33,34] static light scattering,[35] dynamic light scattering,[36] nanoparticle tracking analysis,[37] resonant mass measurement,[38] circular dichroism,[39] light obscuration,[40] electrophoresis,[41] NMR,[42] and mass spectrometry[43] can be applied to study association and small-scale aggregation of peptides in solution.[44,45] Usually, specific structural information on the nature of the aggregates is hard to obtain or not available making a rational approach to understanding and controlling aggregation difficult, if not impossible. Additionally, issues such as evaluating shelf life present significant challenges since problematic aggregation can develop over periods of months. Instead of a rational approach, experience-based trial and error or indirect methods such as chemical or thermal denaturation are used to detect and assist in resolving problems.

In principle, computer modeling can provide structural information on peptide aggregates. For example, computationally identified aggregation hot spots can be useful in identifying potential risks.[46] Computer simulations of molecular systems, including molecular dynamics and related techniques, can provide structural information on the specific nature of the molecular interactions involved in aggregation or increased viscosity. This information can complement experimental information in solving problems for formulating biologics. A number of simulations of protein aggregation have been performed ranging from fully atomistic simulations[47−49] to more simplified representations.[50−53] Atomistic simulations can represent the full flexibility of the peptides; however, these calculations tend to be limited to smaller systems containing only one or a few peptides and for durations much shorter than the association lifetimes. Simplified representations for proteins in solution, typically at a coarse-grained level, permit the simulation of larger systems for much longer times and have been used to understand the increased viscosity in concentrated peptidic solutions.[50−52] However, these approaches are limited in their ability to represent protein flexibility and unfolding.

Mixed resolution simulations in which all of the atoms (AA) or, in some cases, united atoms (UA, typically a single site representing a heavy atom and the hydrogen atoms bonded to it) for the peptide are embedded in a coarse-grained (CG) environment may have a useful role to play by modeling peptide flexibility while gaining increased efficiency, permitting larger system sizes and longer simulations as compared to purely atomistic models.[54−66] Models of solvated peptides using spherically symmetric water sites representing multiple molecules have been published using Gromacs CG water[61] or Martini CG water[57] in conjunction with established AA force fields. Alternatively PACE uses a new effective UA force field[60] for the peptides designed to function well with Martini CG water. All of these approaches truncate the potentials at a fixed distance. These approaches offer large decreases in computer time as compared to atomistic simulations. To improve the quality of the results, both the Gromacs and Martini approaches have been extended to include among other things dipoles[57] and explicit layers of all-atom waters between the proteins and CG water[62,65,67] to better mimic the general electrostatics of the systems and the short-range solvent–peptide interactions. These refinements do lead to improved results; however, they increase the complexity and decrease the efficiency of the models. A related approach using the WatFour model for water,[68] which represents 11 water molecules by a tetrahedron in conjunction with atomistic models for water, has also been explored.[64] Models in which a CG water site represents a single water molecule can be accurate; however, the reductions in computer time as compared to AA water are more modest. The MS-CG approach[69] uses force matching to develop such potentials from atomic simulations, typically, on a system by system basis. The AA/ELBA force field[56,63] uses an AA force field for the peptide combined with a CG model in which each water molecule is represented by a Stockmayer particle (Lennard-Jones potential + fixed dipole) which has been parametrized using solvation thermodynamics and potentials of mean force. The inclusion of dipoles in individual water molecules results in a model that has a speed that is comparable to an all-atom model.

While great progress has been made, a number questions on how the nature of the mixed resolution models impact the speed and accuracy remain. The following questions are examples:

What is a suitable representation for the peptides (AA vs UA)?

Is it necessary to develop new potentials for the peptides in mixed all-atom/coarse-grained models or can well-established, purely AA, potentials be adapted for such studies?

Can models with peptides solvated in simplified CG water molecules approach the accuracy of purely AA models, and if so, what detail in the solvent, including the number of water molecules per site, is appropriate?

Are explicit features, such as dipoles, needed to adequately represent the electrostatic influence of the solvent or are simpler dielectric treatments sufficient?

Are truncated electrostatics adequate or are Ewald treatments needed?

Melittin, a small peptide composed of 26 amino acids, comprises about half of the peptide content in bee venom by weight[70,71] and is one of the most extensively studied peptides of this size. In aqueous solution, melittin exists as a random-coil monomer or a predominantly α-helical tetramer depending on concentration, pH, and ionic strength. Wilcox and Eisenberg proposed that secondary, tertiary, and quaternary structures of melittin form simultaneously.[72] The transition between a helical tetramer and a random-coil monomer has been suggested as a model for protein folding.[72−74] Melittin is also known to associate in membranes[75,76] and is regarded as a model system for understanding peptide interactions within membranes[77] in addition to peptide association in solution.[78] Given its size, melittin lies at the transition between oligopeptide and small protein. The large amount of experimental information available makes melittin an excellent prototypical system to study, yielding insights relevant to issues encountered for small peptidic drugs and for protein aggregation using a relatively small protein especially given that melittin is highly charged (+6), a characteristic of many peptides in biologics. Melittin has been subject to a number of simulation studies.[79−91] Among these is a long time scale atomistic simulation of melittin association[91] in aqueous solution, which demonstrated variations in aggregation behavior with temperature and salt concentration. The stepwise aggregation process and the tendency to form tetramers in solution, at least under some conditions, were rationalized.

Herein, we present a new mixed all-atom/coarse-grained model, which we refer to as AACG, and apply it to studying melittin aggregation in pure water. Our goal is to closely approach the accuracy attained by atomistic models for peptides to model the detailed nature of peptide aggregation including conformational variations and preferences while gaining the efficiencies provided by coarse-graining to make these calculations more practical. In this model, we use an AA representation for the peptides with many valence terms from OPLS,[92−94] while the nonbonded potentials were parametrized anew to function well in a CG environment. While a UA model for the peptide, such as that used by PACE, does reduce computational costs, for biologics the great majority of the system is solvent, which even when it is coarse-grained greatly reduces that speedup. The general trend for single resolution detailed modeling of proteins has been to use purely AA models to attain the most accurate results. We use a relatively fine-grained representation of water for a CG model, one molecule per CG site, in part to enable detailed interactions between water and the peptide including within environments in which specific water molecules can play a role (e.g., when bridging peptide side chains or within binding pockets). The combination of spherical water sites in an inhomogeneous environment with highly charged species poses challenges for calculating electrostatics. We use an inhomogeneous dielectric treatment in which the charges for simple ions are scaled using the dielectric constant for water, while the charges within the peptides are scaled using a smaller dielectric constant (35% of that for water). Additionally, in contrast to most CG models, we use particle mesh Ewald[95] as opposed to truncation to calculate electrostatic interactions, a feature that may be important given the highly charged nature of many biologics. This combination of features is unique to our AACG model and provides a useful perspective on the design considerations for such mixed resolution models.

In the following, we present our simulation procedures followed by an overall description of our AACG model for melittin including its rationalization and functional form. The speedup relative to AA systems will be characterized. We then apply this model to the study of melittin aggregation with variations in salt concentration and temperature, characterizing the nature of the aggregation in both time and structure. Following that, we summarize our results for melittin and suggest further improvements that can be made to our AACG model while highlighting its strengths.

Methods

Simulation Protocols

For studying the behavior of biophysical systems, constructing CG potentials that give the right structure and volume at the same time can be challenging. As a result, CG potentials are sometimes first constructed exclusively for constant particle number, volume, and temperature (NVT) ensemble simulations, where an appropriate volume is predetermined[96,97] so long as the CG model maintains a small positive pressure. We are using this approach for the current generation of AACG models. To obtain the appropriate volume and to take advantage of existing AA tools, we first set up and equilibrate an AA representation of the system using the constant particle number, pressure, and temperature (NPT) ensemble. This system is then converted into an AACG system. The AA simulations are also used to provide some of the information used in parametrizing the AACG potentials.

All of our AACG simulations were performed using the Desmond/Maestro simulation package[98−100] from the suite 2015-3 release. The AA model systems were constructed using the system builder utility included with Desmond. We used cubic simulation boxes with user-defined rather than automatically determined sizes. Protocols for adding water and ions within AA melittin systems, and the subsequent relaxation and simulation, are available in the Supporting Information.

Conversion from an AA representation to an AACG representation involves recognizing the components of the system present, mapping the nonpeptide atoms onto coarse-grained sites, determining the AACG site types for both the peptidic atoms and nonpeptidic sites, and assigning interaction parameters. Figure S1 depicts the atom types used, and Table S1 lists the masses used for the peptidic atoms and CG site types.

Normally, a newly converted AACG system is subjected to an energy minimization before a simulation is started. For AACG simulations, the initial velocities were selected by first randomly sampling from the Boltzmann distribution followed by the removal for any net center of mass motion of the system as a whole. Our NVT ensemble AACG simulations used Nosé Hoover thermostats[101,102] with time constant τ = 1 ps. Production simulations used time steps of 2 fs for valence, Lennard-Jones, and real space Ewald terms, while 6 fs was used for the reciprocal space Ewald calculations. These are the same values used for AA simulations because the fastest motions are related to the valence potentials of the protein which are generally similar.

The simulations reported in this article contained 1, 4, or 20 melittin molecules. Since the folding time for melittin monomers from a random conformation may be longer than what we can easily simulate, unless otherwise noted, simulations were started from monomer conformations from the crystal structure, 2MLT,[103,104] which is predominantly α-helical.

Functional Form for AACG Potential

The AACG potential has a form that is, in part, an adaptation of the AA potential functional form used by the OPLS force field[92] and encompasses valence and nonbonded terms. As with the OPLS force field, the valence terms include the following:for stretch, bend, and torsion contributions, respectively. Here, r is the distance between atoms i and j; θ is the angle formed by atoms i, j, and k; φ is the dihedral angle formed by atoms i, j, k, and l; K, K, K, r0, and θ0 are constants that depend on the types of atoms involved in the interaction.

AACG potentials retain these valence terms among peptidic atoms. In the current study, the species in the environment, water, and simple ions are represented by single sites, so no valence terms are needed for these components. The parameters used for these valence potentials are largely the same as in the OPLS2005 force field,[93] although adjustments have been made to many terms especially for torsions. Tables S2, S3, and S4 list the stretch, bend, and torsional parameters used in the AACG model for melittin.

To permit flexibility during AACG model development, we decided to use a nonbonded potential of the formwhere the c terms are constants. This potential can encompass the normal OPLS Coulombic and Lennard-Jones terms as well as more complex multiwell potentials sometimes used in CG force fields. AACG potentials do not use combining rules for nonbonded interactions, so all nonbonded interactions need to be parametrized. These potentials are truncated at 12 Å. No correction for average dispersion is used in common with some other CG approaches,[96,97,105] and hence, the potentials are fit consistent with the particular cutoff distance used. Other cutoff distances should not be used with these potentials. The current approach is focused on fairly fine-grained CG models for the solvent. For coarser-grained models, it may be necessary to adjust the cutoff distance. The parameters in eq for the nonbonded potentials are listed in Table S5.

The AACG potential retains the pairwise Coulombic and 12–6 Lennard-Jones potentials[92] used in OPLS force fields for atoms separated by three bonds, sometimes referred to as 1–4 interactionswhere A and C are constants that depend on the types of the two interacting atoms i and j; q and q are the charges on these atoms; and ϵ0 is the vacuum permittivity. The AACG 1–4 potential parameters are listed in Table S6.

In AACG simulations, as for AA simulations, hydrogen atoms within the protein were constrained using SHAKE.[106,107] See Table S7 for a list of the types of bonds restrained and the bond lengths used.

AACG Potential Parametrization

The overall parametrization process is described in more detail in the Supporting Information. Here, we focus on the main points that may influence other parametrization efforts. Although the focus of this article is on melittin, the development of the AACG model was based upon many systems including pure POPC membranes, globular proteins, a transmembrane protein, and small peptides. The performance of the AACG force field on this broader class of systems will be documented in a separate article.

We elected to use the force-matching method to generate initial AACG potentials.[69,108,109] Simulations using these potentials rapidly displayed distortions in the protein structures. As a result, we elected to make incremental changes to the potentials. Several hundred cycles of adjustment, each of which typically involved more than one type of interaction within the CG region and between the AA and CG regions, eliminated the most extreme distortions; however, it became clear that to obtain better quality results the potentials between the AA portions of the system also needed to be adjusted. The process of improving the results for the AA portions of the system showed that adequate matching of radial distribution functions for topologically distant interactions did not necessarily lead to good local conformations as evidenced by poor Ramachandran plots for some residues, poor side chain dihedral distributions, and in some cases, excessive, topologically local, hydrogen bonding. As a result, the potential parameters affecting dihedral angle distributions, the coefficients in dihedral angle potentials, and those for 1–4 interactions were adjusted. In some cases, contacts, typically involving hydrogen bonding between polar side chain groups and topologically local backbone N–H or C=O groups, were still too frequent as compared to both the crystal structures and the corresponding AA simulations. In these cases, extra terms were added to the list of 1–4 interactions even though the atoms were more than three bonds apart. Please see Table S8 for a complete list of these interactions. These extra terms are local correction terms and are applied in addition to the general nonbonded interactions between these atoms.

During the process of adjusting the potentials for all of these systems, it was clear that truncation of the Coulomb potentials even with shifting the net potential to 0 at the cutoff was problematic for two reasons: The slope of the Coulombic potential is large out to long distances, and interactions beyond the cutoff distance matter. This problem was resolved in three stages. In the first, we essentially replaced the direct Coulomb contributions in the nonbonded potentials with just the real space part of the Ewald sum, something quite similar to the Wolf approximation which has proven to be fairly effective at reproducing the energetics and dynamics of a number of systems.[110,111] In the second stage, we also included longer range effects by turning on the Fourier space terms via a PME treatment.[95] These terms use the atomistic charges from the OPLS2005 force field for the protein and formal charges on ions in solution. However, since our water model does not inherently include the dielectric constant of water, we reduced the magnitude of the charges provided to the PME terms by the temperature-dependent bulk dielectric constant of water as given byandwhere T and T are the temperature in °C and K, respectively. Here, ϵ(T) is the relative dielectric constant of water as a function of temperature[112] resulting in a temperature dependence in the effective charges, q(T). Similar introductions of Ewald (including PME) electrostatics have been used in previous CG studies.[96,97,113] While this scaling of the reciprocal space terms may appear inconsistent with not explicitly screening the shorter range real space potentials, we note that many of these had been dramatically reduced as part of our iterative adjustments so that they inherently include effective short-range screening due to water. In the third stage, given that there is evidence that the effective dielectric constant within the peptide is smaller than that for water,[114] we scaled the charges in the peptides used in the Fourier space portion of the PME calculation by a different, dielectric constant, 35% of that for water. Interestingly, even though significant adjustments had been made to the potentials before and between these three stages, there was a near universal reduction in the RMSD of the Cα atom coordinates relative to the crystal structures for the peptides used in parametrization with each of these stages. Additional iterative adjustments were made to the potentials after these changes to the electrostatic potential.

Overall roughly 1000 cycles of adjustment were employed to obtain the potentials used in the current study. Aside from the scaling and including the electrostatic real-space terms in the nonbonded potentials, the majority of the nonbonded polynomial potentials have been modified. By comparison, changes to the valence terms involved less than 10% of the torsional potentials and less than 5% of the 1–4 interactions along with the addition of another approximately 2.5% of the terms, representing the more topologically distant pair interactions treated in the same manner as 1–4 interactions. Most adjustments involved reducing or correcting the interactions of atoms with larger charges with the solvent or with each other. The final potentials needed for the melittin simulation are provided in the Supporting Information.

Results

Computational Speed

The main motivation for coarse-graining in molecular simulations is to increase the amount of time simulated per unit of computer or wall clock time. AACG simulations attain a reduction in CPU time for a model involving CG sites due to the relatively small number of atoms mapped onto each CG site for only a portion of the system. The retention of the normal masses and fast local motions in the portions of the system represented atomistically means that with the current approach we must use the same time steps for AA and AACG simulations. Table gives the relative speeds for AA and AACG simulations of the same melittin system. Overall the speedups are 3.6–4 on CPUs and 3.15 on GPGPUs for melittin, values which are similar to those obtained for globular proteins and for a GPCR embedded in a membrane.

Table 1

System Size and Simulation Speed (ns/day) for a 10 mM Melittin System Containing 4 Melittin Molecules, 24 Cl– Ions, and 20,762 Water Molecules Simulated at 310 K for both AA and AACG Representations

	AA	AACG	ratio AACG/AA
number of atoms or sites	64,054	22,530	0.352
1 × 2.2 GHz AMD Opteron processor	0.335 (ns/day)	1.35 (ns/day)	4.03
8 × 2.2 GHz AMD Opteron processors	2.85 (ns/day)	10.2 (ns/day)	3.58
1 NVIDIA GeForce GTX 780 GPGPU	27.5 (ns/day)	86.7 (ns/day)	3.15

CG potentials often yield artificially faster dynamics.[67,97,115] For the AACG model, the diffusion rate for a melittin monomer alone in water with 100 mM NaCl is 2.5 ± 0.4 × 10–9 m2/s, which is 7.2 times faster than in the corresponding AA system, 3.4 ± 0.6 × 10–10 m2/s. This artificial speedup can be problematic if one is interested in detailed kinetics; however, if one is interested in attaining equilibrium, then it can be helpful since the simulation is effectively proceeding faster than it would otherwise. The effective overall speedup can be regarded as the product of the increase in the amount of time formally simulated per unit of CPU time and the increase in the speed of the relevant process for the system being studied. For instance, if one simulates an isolated protein, then the overall speedup will be close to that from the reduction in computer time needed (a factor of 3–4) with little benefit from faster diffusion. However, if the phenomena of interest depend on diffusion, such as the association of melittin molecules leading to aggregation in the current study, then the speedup is closer to 22.

Melittin Simulations Performed

We performed simulations of melittin for the systems listed in Table . All simulations are 200 ns in duration. We do not adjust time scale reported in this paper for the 7 times faster diffusion rate of melittin in water as is sometimes done for coarse-grained models. For simulations with multiple melittin molecules, this factor suggests an upper bound of ≈1.5 μs on the actual time simulated. The techniques used for analyzing these simulations are described in the Supporting Information. The two simulations with a monomer were primarily to explore the solvent-accessible surface area for individually solvated melittin molecules starting from two very different conformations, one with a high initial α-helical content from the melittin crystal structure (2MLT) and the other with a low initial α-helical content (prepared by a short simulation at a high temperature) to provide independent checks on monomer characteristics. Simulations starting from a tetramer extracted from the crystal structure (2MLT) were performed using the AACG, OPLS2005, and OPLS3 force fields to provide an assessment of the quality of the AACG force field for melittin. Given that melittin is believed to form aggregates consisting of approximately 4 monomers, we elected to use simulations containing 20 melittin molecules permitting the study of a number of separate aggregates. As the concentration of melittin in the simulation box is lowered, more solvent needs to be included resulting in a larger simulation box. We chose to use a melittin concentration of 10 mM, which is on the higher end of the range of concentrations used in experimental studies, while being near the lower bound of melittin concentrations that we can practically simulate. The melittin monomers (10 from chain A and 10 from chain B of 2MLT) were manually placed at well-separated positions in a simulation box and with orientations that are not parallel to nearby monomers, as depicted in Figure at 0 ns. We followed the time evolution of the size distribution of melittin aggregates for three different temperatures each with and without 100 mM salt. The methods used to analyze these simulations are described in the Supporting Information.

Table 2

Systems Simulateda

	number of molecules (concentration)
system label	melittin	water	NaCl	temperature (K)	time (ns)
1melXtal	1	11,092	0	310	200
1melRandom	1	11,071	0	310	200
4melXtal	4	21,050	0	310	200
20mel283	20 (10 mM)	105,031	0	283	200
20mel310	20 (10 mM)	105,031	0	310	200
20mel330	20 (10 mM)	105,031	0	330	200
20mel283NaCl	20 (10 mM)	102,530	189 (100 mM)	283	200
20mel310NaCl	20 (10 mM)	102,530	189 (100 mM)	310	200
20mel330NaCl	20 (10 mM)	102,530	189 (100 mM)	330	200

In all cases, six Cl– ions for each melittin molecule were also included to maintain overall charge neutrality.

Figure 1

Aggregation of melittin during the 200 ns course of the “20mel310NaCl” simulation at the top. Sequential snapshots of the entire system, labeled with the simulation time in ns, are depicted. The outline corresponds to the simulation box which is cubic with sides 149 Å in length. NaCl and most of the water present in the simulation are omitted from these images to focus on the behavior of melittin. Most monomers are depicted with blue to mauve colors; however, four monomers that eventually form one aggregate are drawn with brighter colors: yellow, green, red, and cyan. These four monomers are also depicted in close-up at the bottom with all other monomers hidden. Despite starting out widely separated, first contact between any of these monomers, the green and yellow molecules, occurs at 2 ns. By 3 ns, the other pair of molecules, red and cyan, has come into contact, and soon thereafter contact between the red and yellow molecules results in the initial formation of the tetramer. The structure of this aggregate progressively becomes more compact as the simulation proceeds.

Melittin Tetramer Characterization

To compare the AACG force field results with those produced by all-atom force fields, we also ran equivalent simulations of “4melXtal” (a solvated melittin tetramer starting from the crystal structure) using the all-atom force fields, OPLS3, and OPLS2005. Figure shows the secondary structure profile per residue as well as the secondary structure content averaged over all residues. The starting structure of melittin (taken from 2MLT) is highly ordered with 90% α-helical content. Previous NMR studies have shown that melittin tetramers in aqueous solution consist of two helical segments (residues 2–11 and 13–23) and have an unstructured C-terminal region.[116] Melittin tetramers in our AACG simulations have two helical segments and follow a per-residue secondary structure profile similar to the one we observed in the OPLS3 (all-atom force field) simulation. Melittin tetramers in OPLS2005 (another all-atom force field) simulations are less ordered. AACG gives an overall % helicity (57%) which compares favorably with the experimentally measured % helicity of melittin tetramers in aqueous solution, about 60%.[73] The buried surface area (Figure ) and contact (Figure ) profiles from the AACG simulation are similar to those observed in OPLS3 and OPLS2005 simulations, with the exception of the C-terminal residues (23–26) which are known to be unstructured in aqueous solution.[116]

Figure 2

Percent helicity by residue number in melittin for simulations starting from the tetramer from the crystal “4melXtal” using (a) AACG, (b) OPLS3, and (c) OPLS2005 force fields averaged over the second half of the simulations (100–200 ns). For comparison, the % helicity from the tetramer in the crystal structure is given in panel (d). The overall helicity (averaged over all residues) is given at the top right corner of each plot. The analysis of the simulations gave 0% β-strand content in all cases.

Figure 3

Buried surface area (in Å2) per residue from the “4melXtal” simulations for the AACG, OPLS3, and OPLS2005 force fields averaged over the second half of the simulations (100–200 ns).

Figure 4

Number of inter- and intrachain contacts for simulations starting from the “4melXal” structure for the AACG, OPLS3, and OPLS2005 force fields. These results are averages over the second half of the simulation (100–200 ns).

Riniker et al.[61] found a large increase (41%) in the total number of intraprotein hydrogen bonds when a supramolecular coarse grain solvent (four water molecules are represented by one CG site) was used to solvate an all-atom protein as opposed to when an AA solvent is used. We integrated the radial distribution function between polar hydrogens and hydrogen bond acceptors involving only melittin atoms (PHA) for the four systems corresponding to (a) through (d) in Figure to obtain comparable information on polar hydrogen contacts from our simulations. The results are described in detail in the Supporting Information. Our PHA counts are also higher for the AACG model relative to the OPLS3 (5.6%) and OPLS2005 (28%) simulations, although these increases are smaller than found by Riniker. Interestingly, the PHA count for our AACG simulation is lower (−3.9%) than that for the tetramer in the crystal.

Aggregation Events Examined Using Snapshots

The evolution of the entire “20mel310NaCl” system as a function of time is depicted in Figure (top). The simulation started from a configuration with the monomers dispersed throughout the cell without intermelittin contacts. The formation of aggregates proceeds quickly with the initial intermolecular contacts being established within 1 ns. Within 8 ns, nearly all of association events observed within the 200 ns simulation have occurred. In Figure (bottom), we follow the evolution of four monomers that eventually form one aggregate. The tetramers observed in our simulations are dynamic, internally disordered, and on average more bent than in crystalline tetramers. They show amphiphilic folded structures but are less compact than crystalline tetramers (Figure S2). After monomers come into contact to become a member of an aggregate, they rarely separate.

Salt and Temperature Effects on Aggregation

Time evolution of size distribution of melittin aggregates, under four different conditions, is shown in Figure . Nearly all aggregation events happen in the first 10–30 ns for the systems containing 20 melittin molecules as compared to a recently published AA simulation study of the aggregation of 4 melittin molecules with events occurring on time scales longer than 200 ns,[91] illustrating the effect of the higher diffusion coefficient in AACG models. Figure shows the distribution of melittin counts among the different sizes of aggregates at the end of the simulations (200 ns) for all six simulations with 20 monomers. The presence of salt clearly enhances aggregation, leading to the formation of larger aggregates at all three temperatures studied. Since the addition of salt should screen the net electrostatic repulsions between the highly positively charged melittin monomers, this shift to larger aggregates with the addition of salt is expected. Figure also displays a weak tendency to form larger aggregates as the temperature is lowered, both with and without added salt. This trend makes sense as higher temperatures should favor states with more disorder and thus smaller aggregates. Similar trends were noted in the aggregation of four melittin molecules in a recently published study using an AA model,[91] although in the current study, we observe the formation of some aggregates containing five melittin molecules.

Figure 5

Characterization of aggregation during the simulations. Panels (a) and (b) are for 310 K without and with NaCl, respectively. Panel (c) is for the least aggregating condition, 330 K without salt, and panel (d) is for the most aggregating condition, 283 K with NaCl. The left-hand plot for each condition tracks the number of aggregates of each size as a function of time, while the right-hand bar graphs give the breakdown of the monomers present in aggregates of each size at times: 10 ps, 10 ns, 20 ns, and 200 ns. Blue is used for melittin in unaggregated monomers, red for dimers, green for trimers, purple for tetramers and light blue for pentamers.

Figure 6

Effect of added salt and different temperatures on aggregation. These bar graphs depict the number of monomers in aggregates for each size at the end (200 ns) of the six simulations containing 20 melittin monomers. Blue is used for melittin in unaggregated monomers, red for dimers, green for trimers, purple for tetramers and light blue for pentamers.

Buried Surface Area

To better understand the nature of melittin aggregates, we examine residue-based properties including buried surface area (Figure ), solvent-exposed surface area (Figure S3), and number of inter- and intrachain residue contacts (Figure ) for the least aggregating condition (330 K without salt, “20mel330”) and for the most aggregating condition (283 K with NaCl, “20mel283NaCl”). For comparison, the same analysis is performed on a melittin tetramer extracted from the crystal structure. Since TRP19 has been identified as a key residue for characterizing aggregation, we examine a breakdown of its contacts with other residues (Figure ) as well. For the averaged buried surface area per residue in Figure , we plot time averages over all monomers for 0–10, 50–100, and 100–200 ns. Here, 0–10 ns represents the early-to-midstages of aggregate formation during which 80% or more of the monomers come into contact and aggregate growth proceeds (Figure ). The buried surface areas for 0–10 ns are distinctly lower than those for 50–100 ns, which in turn is quite similar to those for 100–200 ns. So the melittin aggregates seem to be mature throughout the second half of our simulations, and we will omit the 50–100 ns in subsequent discussions.

Figure 7

Average buried surface area (in Å2) per residue. The top panel gives the results from the least aggregating simulation, “20mel330”, while the middle panel contains the results for the most aggregating simulation, “20mel283NaCl”. The bottom panel gives the results obtained using a melittin tetramer extracted from the crystal structure (2MLT).

Figure 8

Number of inter- and intrachain contacts. Interchain contacts based upon whole residues (including backbone atoms) for “20mel330”, “20mel283NaCl”, and the tetramer from the crystal structure are given in panels (a), (b), and (c), respectively, while the corresponding intrachain contacts for these systems are given in panels (d), (e), and (f), respectively. LEU9 and TRP19 undergo the greatest increase in interchain contacts upon increased aggregation.

Figure 9

Number of residues in contact with TRP19 averaged over the second half of the simulation (100–200 ns). The distance cutoff for contacts was 4.0 Å. The top three panels are for intermolecular contacts, and the bottom three panels are for intramolecular contacts.

In the crystal structure of a melittin tetramer (Figure ), ILE2 and TRP19 have the highest buried surface areas. Hydrophobic residues VAL5, LEU6, VAL8, LEU9, LEU13, LEU16, and ILE20 are also highly buried. All of these are among the most buried residues in our simulations. Among charged residues, LYS23 and ARG24 are the most buried in the crystal tetramer. With a few exceptions, the buried surface area of the residues (Figure ) is noticeably higher for the most aggregating condition (“20mel283NaCl”) as compared to the least aggregating condition (“20mel330”) with the largest increases for most of the hydrophobic residues ILE2, LEU6, VAL8, LEU9, LEU13, and LEU16. The polar residues are generally more buried in “20mel283NaCl” than in “20mel330” with the exception of SER18, while the charged residues are relatively unaffected with the exception of ARG24, which is significantly more buried in “20mel283NaCl”. TRP19 is highly buried in the melittin aggregates formed during our simulations. The buried surface area of TRP19 roughly doubles to 83 Å2, about 81% of that for the tetramer for the crystal structure, under the most aggregating condition (“20mel283NaCl”) becoming the residue with the highest buried surface area. LYS23 is relatively unburied in our simulations, suggesting that packing effects may be responsible for its large buried surface area in the crystal.

The general pattern of the exposed surface area, broken down by residue, is similar in both of the simulations of individual monomers with higher levels of hydration near the ends of the chain for residues 7–11 and for residues 21–24 (Figure S3). The most exposed residues in the crystal also tend to be the most exposed in solution. Interestingly, LYS21 and more distinctly LYS23 have smaller surface areas in the crystal than in our simulations, perhaps reflecting packing effects specific to the crystal structure. Overall, these results suggest that solvent exposure of residues in melittin depends mostly on the nature of the residue (i.e., hydrophobic vs hydrophilic) and proximity to either end of the chain with aggregation significantly reducing the exposure of the more hydrophobic residues.

Residue Contacts

Interchain contacts follow a similar overall pattern under the least aggregating conditions (“20mel330”, Figure a) and the most aggregating conditions (“20mel283NaCl”, Figure b) with generally higher counts for the latter except for residue VAL5 which drops slightly. Most residues with high interchain contacts in the crystal (Figure c) also have high levels of contact within the aggregates with the exception of LYS21 which loses essentially all contacts in solution. The level of intramolecular contacts for the two aggregating conditions (Figure d and e) closely correspond with a small overall drop under conditions of increased aggregation (Figure e). Consistent with the shifts in buried and exposed surface areas for LYS23, there is a dramatic drop in the number of intrachain contacts for LYS23 from about 1.5 to less than 0.25 between the crystal structure and aggregates in solution.

We have broken out the contacts TRP19 forms by residue within the aggregates and the crystal structure (Figure ). The overall level of intermolecular contacts for TRP19 is dramatically higher under the most aggregating condition (Figure b) as compared to the least aggregating condition (Figure a) in our simulations, while the level of intramolecular contacts is similar (Figure d and e).

In the crystal structure of a melittin tetramer, TRP19 of each melittin molecule is in contact with five residues (ILE2, TRP19, ILE20, LYS23, and ARG24) from other melittin molecules (Figure c) and residues ARG22 and LYS23 within the same molecule (Figure f). In aggregates formed during our simulations, the TRP19 of each melittin molecule comes in contact with most of the residues of other melittin molecules. This indicates that melittin aggregates in solution are less ordered than in the crystalline tetramer. The most frequently observed contacts are with TRP19 and ARG24 of other melittin molecules which also occurs in the crystal structure, while contacts with ILE2, ILE20, and LYS23 are comparatively uncommon in contrast to the crystal structure. Intermolecular contacts with ARG22, ALA15, and LEU16 are moderately common in solution, while they are completely absent in the crystal structure. In the melittin aggregates formed during our simulations, the most common intramolecular contact residue for TRP19 is also ARG22, while intramolecular contact with LYS23 is rare. The drop in both intra- and intermolecular contacts between crystal and aggregated forms for LYS23 with TRP19 is consistent with the shifts in buried and exposed surface area for LYS23 noted earlier.

TRP-ARG and TRP-LYS interactions are π-cation interactions, a well-known interaction motif in structural biology, with TRP-ARG being more frequently observed.[117] So the extensive contacts between TRP and ARG in our simulations (Figure ) are expected, as is the preference for ARG-TRP contacts over LYS-TRP contacts. Intermolecular TRP-TRP interactions are π-stacking interactions, which is another type of well-recognized interaction motif in folded proteins and is observed in many crystal structures. Figure depicts such interactions which in our simulations seem to be tightly coupled to π-cation interactions.

Figure 10

π-cation interactions involving TRP19 observed during our simulations. The left panel is a close-up view showing an example. TRP19 and ARG22 of the same melittin molecule are colored dark orange, and ARG24 from a neighboring melittin molecule is colored green. The right panel shows this same example in the environment of the entire tetramer. The side chains of TRP19, ARG22, and ARG24 from all four melittin molecules are depicted in CPK rendering. The carbon atoms and ribbon color are given a distinct color (dark orange, green, red, light blue) for each melittin molecule to make it easier to see which interactions are intramolecular and which are intermolecular. Intermolecular and intramolecular π-cation interactions are observed 54% and 83% of the time, respectively, in “20mel283NaCl” (the most aggregating condition). They are observed 18% and 91% of the time, respectively, in “20mel330” (the least aggregating condition).

Figure 11

π-stacking interactions observed during our simulations. In each of the top images, a pair of TRP19 residues is shown. In the bottom images, surrounding ARG residues are also shown. As in Figure , the carbon atoms and ribbons for each molecule are given a distinct color. These interactions are observed 50% of the time in “20mel283NaCl” and 15% in “20mel330”.

Discussion

Aggregation of Melittin

In our simulations starting from well-separated melittin molecules at 10 mM, we have observed aggregation of melittin into dimers, trimers, tetramers, and pentamers. Aggregation events are readily observed, with most of them occurring within the first 10–30 ns (roughly 10% of the total time simulated) as compared to a recently published AA simulation study of the aggregation of four melittin molecules with events occurring on time scales longer than 200 ns,[91] illustrating the effect of the higher diffusion coefficient in AACG models. Since we rarely see dissociation events, we cannot conclude that we have reached the equilibrium size distribution. On the other hand, the frequency of association events has dramatically decreased by the end of our simulations suggesting that our simulation time is long enough to characterize melittin aggregation. In order to more thoroughly study equilibrium behavior, we would need to either lengthen the duration of simulation by orders of magnitude or employ a method for enhancing the frequency of dissociation events. We observe both a clear increase in aggregation with the addition of salt and a subtle but noticeable increase in aggregation as the temperature is lowered.

Residue Characterization

We observe significant increases in the amount of buried surface area after the initial period of large-scale monomer aggregation. This is consistent with the idea that the aggregates are fluid and that the initial monomer contacts evolve. All residues with little buried surface area in the tetrameric structure from the crystal structure, especially the C-terminal residues, are significantly more buried in the simulated aggregates because of the overall disordered state of the aggregates.

The residue properties for the most aggregating and the least aggregating conditions follow similar patterns except that the results are shifted toward compactness for the former, namely, higher buried surface areas, slightly lower solvent exposed surface areas, slightly higher levels of inter chain contacts, and slightly lower levels of intrachain contacts. For the buried surface area and interchain contacts, there are dramatic changes between the initial 0–10 ns and final 100–200 ns period. This indicates large-scale evolution of the systems and the need for relatively long simulations. Melittin aggregates formed during our simulations have amphiphilic folded structures, but they are less compact and less ordered than crystalline tetramers. We observed π-cation and π stacking interactions involving TRP19 as key characteristics of melittin aggregates in solution.

Conformational State of Melittin Aggregates

Hartings et al.[118] obtained fluorescence and circular dichroism data on single-residue mutants of melittin, V5F, and V5Y(NO2) and observed that helicity increases as the concentration of melittin is raised from 2 to 20 μM but with no blue shift in TRP19 fluorescence indicating that TRP19 side chain remains largely solvent-exposed. From these observations, they suggested that these peptides assemble into a dimer rather than a tetramer. An alternate explanation may be that they assemble into less ordered, more dynamic tetramers, in which TRP19 is more solvent exposed than in the crystalline tetramer as we are observing in our simulations of melittin. Hagihara et al.[119] found that melittin aggregates are in the molten globule state based on circular dichroism and differential scanning calorimetry results as well as the radius of gyration. Our results are consistent with the molten globule state.

The crystal structure results in Figure –9, show a strong contrast with many of the residues having values of 0, while the largest values are comparable to or larger than the largest ones obtained from simulation. Generally, the simulation results for surface areas and contacts are distinctly nonzero for all residues where the measure applies in contrast to the crystal structure, again consistent with a dynamic molten globule model for melittin aggregates. However, in most cases, residues with large values of buried surface area in the crystal also have large values in the simulations.

Conclusions

We have developed a new AACG (all-atom/coarse-grained) model. The structural results obtained using the AACG force field for a melittin tetramer in solution compare favorably with those for the OPLS2005 and OPLS3 AA force fields. We subsequently applied the AACG force field to study melittin aggregation behavior in aqueous solution obtaining results that are in good agreement with the experiment. Starting from a configuration consisting of 20 well-separated melittin molecules in a periodic box, we are able to observe the formation of melittin dimers, trimers, tetramers, and pentamers during the time scale of our simulations. The extent of aggregation, monitored by the size distribution over time, the buried surface area, and the counts of residues that come in contact with each residue is found to increase when salt is present and to a lesser extent as temperature decreases. We found that TRP19 and hydrophobic residues are important for the formation of aggregates and noted the prominent role of π-cation and π-stacking interactions. Overall our results are consistent with the molten globule model for melittin aggregates.[119] Our results, including the influence of salt and temperature on melittin aggregation, suggest that the AACG model described in this paper performs reasonably well at mimicking the behavior of melittin in water.

The 3–4 fold reduction in CPU time, as compared to AA simulations, can be useful in some cases when studying the behavior of individual peptides; however, one should weigh this benefit against the more approximate results obtained. When the diffusion of molecules is important for the phenomena being studied, such as for peptide aggregation, the 7 times higher diffusion rates for peptides in our AACG model are relevant, resulting in simulations that are at least 20 times more efficient than a simulation of the same system using an AA model. Simulating large systems, particularly for larger aggregating peptides or proteins (e.g., globular proteins or antibodies) at the concentrations of interest, would require systems with about 107 atoms. Systems with this many particles are hard to fit into the memory on current GPGPU cards. AACG models require roughly 3 times fewer particles than AA models to mimic the same amount of aqueous peptidic solution, permitting substantially larger systems to be simulated. The rarity of monomer separations from aggregates in our AACG studies of melittin suggests that care is needed when interpreting simulations of peptide aggregation.

In common with the PACE model,[60] we have elected to derive anew interactions or modify interactions, particularly the nonbonded interactions, within the finer grain portions of our mixed resolution models as compared to existing UA or AA models. This approach seems reasonable given that the interactions between coarse-grained solvent particles are typically significantly weaker than those between atoms in purely AA models. Other mixed resolution models may benefit from similar adjustments.

While the various mixed resolution models for proteins have their own strengths, our AACG model has advantages for studying the aggregation of peptides. The most significant of these advantages is the built-in compatibility with particle mesh Ewald treatments, while most other models truncate electrostatic interactions sometimes in combination with a reaction field method treatment.[57,60−62,67,69] Given that experimentally controlling the net charge on peptides is a key tool for controlling aggregation issues for biologics, a careful treatment of electrostatics on long-length scales would be critical. Evidence for this is available from a few sources including the general shift to Ewald methods in atomistic simulations to avoid artifacts noticed when cutoffs are used[120] and our own observations that as we progressively introduced particle mesh Ewald into our calculations our results systematically improved. Other simulation approaches with explicit solvent particles that do not rely on explicit representation of the solvent’s electrostatics (such as dipoles or partial charges) might benefit from our two-level dielectric treatment of electrostatics which shares some of the spirit of how electrostatics is handled in implicit solvation methods such as the Generalized Born[121] and Poisson–Boltzmann[122] approaches.

The coverage of the AACG force field is currently limited to peptides, proteins, POPC lipids, water, a small set of ions and cofactors, and a few types of ligands and cosolvents (see the Supporting Information for a more detailed list). Given that the model was developed using peptides, globular proteins, and GPCR, it should be reasonably transferable within these types of systems. To explicitly test and verify transferability of the AACG model beyond peptides in aqueous solution, we plan to characterize AACG simulations of globular and membrane-bound proteins in future publications. The AACG approach has required and will continue to require time-consuming parametrization. Following the experience of other force fields after their introduction, it will require further refinement. We plan to extend the AACG modeling approach by further refining the peptidic potentials and developing a mechanism to make fitting parameters for new atom types easier so that coverage of chemical space can be more readily extended. In addition, we plan to relax the requirement for NVT simulations to permit the direct construction and simulation of AACG systems without the need to prepare and relax an AA system prior to creating an AACG representation. It should be possible to achieve further speedup of AACG simulations by using hydrogen mass redistribution.[123,124]

We believe that the AACG model fills a useful niche in the modeling field between fully atomistic models and fully coarse-grained models where electrostatics and molecular flexibility are combined with longer time frame phenomena. We anticipate that the AACG model will prove useful for studying the association of proteins in many biological processes and also in biologics formulations.