Transferring the PRIMO Coarse-Grained Force Field to the Membrane Environment: Simulations of Membrane Proteins and Helix-Helix Association.

An extension of the recently developed PRIMO coarse-grained force field to membrane environments, PRIMO-M, is described. The membrane environment is modeled with the heterogeneous dielectric generalized Born (HDGB) methodology that simply replaces the standard generalized Born model in PRIMO without further parametrization. The resulting model was validated by comparing amino acid insertion free energy profiles and application in molecular dynamics simulations of membrane proteins and membrane-interacting peptides. Membrane proteins with 148-661 amino acids show stable root-mean-squared-deviations (RMSD) between 2 and 4 Å for most systems. Transmembrane helical peptides maintain helical shape and exhibit tilt angles in good agreement with experimental or other simulation data. The association of two glycophorin A (GpA) helices was simulated using replica exchange molecular dynamics simulations yielding the correct dimer structure with a crossing angle in agreement with previous studies. Finally, conformational sampling of the influenza fusion peptide also generates structures in agreement with previous studies. Overall, these findings suggest that PRIMO-M can be used to study membrane bound peptides and proteins and validates the transferable nature of the PRIMO coarse-grained force field.

Introduction

Membrane proteins play significant roles in many cellular processes, such as intercellular communications, molecular transports, and signal transductions.[1,2] About 8000 membrane proteins are encoded in the human genome (∼25% of all genes), and they are also the targets of about 60% of the currently available drugs.[3,4] Hence, understanding their structures, dynamics, and functions are of the utmost importance in the biological and pharmacological sciences.[5] However, a detailed study of membrane proteins is challenging.[6−8] Because of the complexity of lipid bilayers, it has proved difficult to map details of protein–membrane interactions using experimental techniques.[9] Molecular dynamics simulation can serve as a complementary tool to illustrate the structural and dynamic details of proteins at the atomistic level.[9−12] All-atom simulations with explicit lipid and solvent are the most accurate, but they are often limited to simulation lengths that are short compared to the slow relaxation kinetics of lipid bilayers, and hence convergence of such simulations is often a concern. To overcome such limitations, a common strategy involves the use of either coarse-graining or implicit solvent techniques.

In coarse-grained (CG) models, the idea is to reduce the degree of freedom by grouping several atoms together into single sites to reduce the degrees of freedom of the system. A wide range of CG models for biomolecules and lipids have been developed in the past decade to study protein folding, aggregation, and design;[13−18] bilayer structure and dynamics;[19−26] and protein interactions with lipid bilayers.[27−31] The most popular and widely used coarse-grained model for membrane-bound proteins and peptides is MARTINI,[21,27,32] which has been employed to a wide range of applications. The MARTINI model is based on a four-to-one mapping; i.e., on average four heavy atoms plus associated hydrogens are represented by a single CG site that interacts through an empirical interaction potential. The MARTINI model represents both lipids and the surrounding solvent in an explicit, coarse-grained fashion. This provides some degree of transferability between membrane and water environments, but the model is overall too coarse to accurately reflect amino-acid specific secondary structure propensities and requires knowledge-based secondary structure constraints. Other available CG models suffer from similar limitations.[22,30,33−38] Following the mapping scheme of MARTINI, Dal Peraro and co-workers[39−41] have developed a coarse-grained model for proteins in which they have introduced an explicit dipole defined by three consecutive Cα beads along with the treatment of nonradial dipole–dipole interactions in dynamics. This model can maintain stable secondary structure elements of unspecific proteins and sample conformational transitions. In order to improve the ability to accurately sample peptide conformations in the context of membranes, some groups have explored multiscale strategies where the peptides are represented at higher levels of resolution,[31] for example with united-atom protein models interacting with MARTINI CG lipids and solvent.[42]

Implicit membrane models focus on eliminating the solvent/lipid degrees of freedom entirely. Because of the slow relaxation of lipids in response to peptide/protein dynamics, this approach can lead to significantly accelerated sampling. In the past, implicit membrane models have been used primarily in conjunction with atomistic representations of membrane-interacting peptides or proteins.[43−46] In some cases, implicit models have also been combined with CG models, but typically they are used to model only the aqueous solvent while maintaining an explicit representation of the lipid bilayer.[47−49]

Here, we are presenting a new model where a CG representation of the peptides is integrated with a fully implicit membrane model. Such a model combines the computational speedup of CG models with the kinetic acceleration in the absence of explicit lipids. The CG model is our recently introduced PRIMO force field that performs similar to all-atom force fields in aqueous solvent but with a third to a fourth of the number of particles.[50,51] PRIMO is designed as a transferable model that does not require system-specific parametrization or constraints. It relies on a generalized Born (GB) term to describe solvation effects in aqueous solvent. This makes it possible to switch to membrane environments by simply replacing the standard GB term with an implicit membrane GB model. Here, we are describing the result of combining PRIMO with the heterogeneous dielectric generalized Born (HDGB) model[52] (PRIMO-M) to describe the energetics and dynamics of membrane-bound peptides and proteins.

In the next section, the PRIMO force field is described briefly along with the heterogeneous dielectric generalized Born implicit membrane model and the simulation methods employed in this study, followed by a presentation and discussion of results.

Model

PRIMO Force Field

Since the PRIMO force field and its parameters are described in detail in our previous paper,[50,51] we limit ourselves to a brief introduction. The CG interaction sites in PRIMO were chosen to allow an analytical reconstruction of all-atom model based on molecular bonding geometries to near-atomistic accuracies.[50,53] In order to preserve the backbone hydrogen bonding interactions, the backbone in PRIMO is represented with N, Cα, and a combined carbonyl site (CO) placed at the geometric center of the carbonyl C and O atoms. This is particularly crucial for an accurate description of the secondary structures of a protein. Nonglycine side chains (SC) are represented with one to five CG sites.

The PRIMO energy function follows an all-atom-like physically motivated force field with additional terms for a combined generalized Born/atomic solvation parameter (GB/ASP) implicit solvent, an explicit angle- and distance-based hydrogen bonding interaction potential, and spline-based bonded potentials to maintain correct bond geometries at the coarse-grained level. The PRIMO force field is expressed as in eq 1.

Bonded interactions between PRIMO sites correspond to both real covalent bonds and virtual bonds. Standard harmonic potentials are employed to model the real covalent bonds, such as 1–2 (bond) and 1–3 (angle) interactions. Otherwise, distance-dependent spline-interpolated potentials are used to reproduce nonharmonic functions and multiple minima. The bonded interactions between virtual sites and the primary CG sites are described by standard harmonic potentials. It should be noted here that the virtual atoms do not participate in nonbonded interactions, and they are reconstructed on the fly to improve local molecular geometries. In addition to the one-dimensional 1–4 terms, PRIMO also uses two-dimensional spline-interpolated CMAP potentials[54] to couple the sampling of CO–N–CA–CO and N–CA–CO–N torsions, and thereby the sampling of φ/ψ backbone torsions is controlled in PRIMO as in the atomistic CHARMM force field. Nonbonded terms in PRIMO consist of Lennard-Jones, electrostatic interactions, and an explicit hydrogen bonding term, which is the only knowledge-based term in the PRIMO force field. The force field is mainly parametrized based on direct comparison with all-atom simulations in a bottom-up fashion as described in our previous paper.[51]

HDGB Model

The HDGB model is based on the GBMV method, and as the details of the HDGB model have been described elsewhere,[52,55] we limit ourselves to a brief summary of the relevant features here. The solvation free energy of a solute in any solvent environment can be decomposed into the electrostatic contribution ΔGsolvGB and the nonpolar contribution ΔGsolvnp.In the HDGB model, the following modified expression for the polar component of the solvation free energy is used to introduce a heterogeneous dielectric environmentwhere q are charges of CG beads, n is the total number of CG particles, r are distances between particles i and j, α are the effective Born radii, and F is a dimensionless quantity set to be 8. The effective dielectric constant (ε) for each atom varies as a function of the membrane insertion depth (distance from the membrane center), along the membrane normal. The ε profile was generated initially by solving the Poisson equation for a monovalent spherical probe ion that is moved across the heterogeneous multicontinuum dielectric environments[52] and subsequently optimized by Sayadi and Feig[56] to match free energies of insertion for amino acid side chain analogues from experiment and explicit membrane simulations.

For heterogeneous environments, such as biomembranes, the nonpolar component of the solvation free energy is especially important as it allows molecules to remain in the interior of biological membranes (i.e., low dielectric regions) by compensating the polar solvation free energy. In the HDGB formalism, the nonpolar component of the solvation free energy is described by the solvent accessible surface area model along with variable surface tensions:where A is the solvent accessible surface area of the ith atom, z is the distance of atoms i from the membrane center along the membrane normal, γ is the coefficient of surface tension for different atom types, and S(z) is a switching function that is used to reflect the change of the surface tension along the membrane normal. The switching function S(z) was determined initially by matching the free energy profile of insertion of O2 into lipid bilayers obtained from explicit lipid simulations and later optimized along with the ε profile to match insertion free energies for amino acid side chain analogues.[52,56]

PRIMO-M

In the membrane version of PRIMO, PRIMO-M, the standard GB term is simply replaced by the HDGB model, and the SASA term is scaled in a z-dependent fashion according to eq 4. No further changes were made in the parameters of the original PRIMO model and, as in the original PRIMO model, the radii used in the solvation terms are the PRIMO Lennard-Jones radii.

Simulation Methods and Test Systems

PRIMO-M was applied to a number of test cases to evaluate its performance and applicability. We studied the stability of 13 membrane proteins (see Table 1) by carrying out MD simulations in the membrane environment. Furthermore, we carried out simulations of WALP peptides that have been widely studied both experimentally and theoretically and two other helical peptides (see Table 2). We paid particular attention to the tilt angle of these peptides in the membrane bilayer. We also studied the dimerization of glycophorin A (GpA), and finally, we examined the conformational sampling of the membrane-bound influenza fusion peptide (IFP) using replica exchange molecular dynamics simulations.

Table 1

List of Membrane Proteins Simulated Using PRIMO-M Methodology

no	protein name	PDB ID	code	resolution (Å)	TM sec. struct.	residues
1	bacteriorhodopsin	1QHJ	BRD7	1.9	21	228
2	V-ATPase	2BL2	VATP	2.1	4	156
3	rhomboid intramembrane protease	4H1D	GlpG	2.9	6	173
4	lactose permease	2CFQ	LacY	3.0	12	417
5	nucleobase-cation-symport 1 transporter	2JLN	Mhp1	2.9	12	463
6	Na+/H+ antiporter	1ZCD	NhaA	3.5	14	376
7	human aquaporin 5	3D9S	AQP5	2.0	8	245
8	adhesion/invasin	1K24	opcA	2.0	10	249
9	β1-adrenergic receptor GPCR	2VT4	adr	2.7	7	276
10	outer membrane protease	1I78	ompT	2.6	10	297
11	outer membrane protein X	1QJ8	ompX	1.9	8	148
12	outer membrane transporter	1KMO	FecA	2.0	22	661
13	intramembrane protease	2IC8	GlpG	2.1	6	182

Table 2

Overview of PRIMO-REMD Simulations of Helical Peptides

system	sequence	time (ns)	replicas	temp. range (K)	time for analysis (ns)
WALP23	GWW(LA)8LWWA	50	8	300–400	20–50
WALP19	GWW(LA)6LWWA	50	8	300–400	20–50
GWALP23	GGALW(LA)6LWLAGA	50	8	300–400	20–50
KWALP23	GKALW(LA)6LWLAKA	50	8	300–400	20–50
AChR M2	GSEKMSTAISVLLAQAVFLLLTSQR	50	8	300–400	20–50
NMDA M2	GSNGDALTLSAMWFSWGVLLNSGIGE	50	8	300–400	20–50
IFP	GLFGAIAGFIENGWEGMIDG	100	8	300–500	20–100
GpA	EITLIIFGVMAGVIGTILLISYGIR	400	12	270–500	200–400

In this study, periodic dielectric and nonpolar profiles were used as described previously.[57] This ensures a finite peptide concentration along the membrane normal. Furthermore, if a peptide dissociates and diffuses away from the membrane interface, because of the periodic nature of the membrane, it eventually interacts with the membrane. The thickness of the membrane was chosen to be 25 Å, measured as the distance between the center of the membrane (z = 0) and where ε = 80 is reached.

All the proteins and transmembrane peptides were capped with an acetyl group on the C-terminus and an N-methyl amide group on the N-terminus. MD simulations with PRIMO-M were carried out with the latest version of the heterogeneous dielectric generalized Born (HDGB) implicit membrane model to reflect the membrane environment[52] using GB parameters described in detail previously.[55,56] We downloaded the initial structures from the Protein Data Bank and then energy-minimized to release the side-chain strains before beginning the simulations. The minimized structures were then placed so that its principal axis coincided with the bilayer normal (z-axis) and its center of the mass was at the origin of the membrane bilayer. The Langevin thermostat with a friction coefficient of 10 ps–1 was employed to maintain a temperature of 300 K. The time step for leapfrog Verlet integrator was set to 4 fs. Nonbonded interactions were cut off at 18 Å, with smooth switching to zero starting at 16 Å. The nonbonded interaction list was maintained up to 20 Å. A similar protocol was followed for simulations with PRIMO in an implicit aqueous environment.[51] It should be noted here that it was previously found that the simulation of the bacteriorhodopsin monomer (PDB code 1QHJ) with the HDGB model and a short cutoff of only 16 Å led to serious artifacts.[55] In that case, the protein adopted a horizontal orientation during the simulation from its vertical transmembrane orientation. However, simulations of the bacteriorhodopsin monomer with larger cutoff distances of 18 and 38 Å did not lead to rotation of the monomer. Hence, we used a cutoff value of 18 Å for the nonbonded interactions in PRIMO-M simulations which did not result in any apparent artifacts (see below).

The equilibration protocol for all of the PRIMO-M/MD simulations generally consisted of initial minimization followed by stepwise heating to 300 K. The heating phase for each of the 13 membrane proteins consisted of six stages, and in each equilibration stage, the temperature of the peptide or proteins was increased by 50 K, and 40 ps MD simulations were carried out. During the heating phase, a harmonic constraint was employed to keep the Cα heavy atoms fixed. However, production simulations were carried out completely unrestrained.

To increase conformational space sampling efficiency, we have carried out replica exchange molecular dynamics (REMD) simulations of peptides listed in Table 2. For all but the GpA dimer, eight replicas were employed that were distributed over a temperature range of 300 to 400 K or 300 to 500 K, spaced exponentially. In the case of the GpA dimer, 12 replicas were used spanning a temperature range from 270 to 500 K. Langevin dynamics with a friction coefficient of 10.0 ps–1 was used to maintain the target temperatures and ensure random drifts of the peptide. In all cases, exchange moves between adjacent replicas were attempted every 5 ps, leading to an acceptance ratio for successful exchanges of 60–65%. Each replica was simulated for 50 to 400 ns for different peptides, and only the final 30 to 200 ns data were used for analysis for different peptides (see Table 2).

PRIMO-M simulations were carried out using version c36a4 of the CHARMM macromolecular modeling package[58] where the PRIMO model is implemented. The MMTSB (Multiscale Modeling Tools for Structural Biology) Tool Set[59] in combination with CHARMM was employed for all analyses.

Results and Discussion

The PRIMO-M model was implemented as described above and applied to a number of membrane-related test cases as detailed in the following.

Solvation Free Energy Profiles for Membrane Insertion of Amino Acids

The performance of the PRIMO/HDGB model was first tested by comparing amino acid insertion free energy profiles for an ideal alanine-based transmembrane (TM) helix of the type A49XA50 with the corresponding all-atom CHARMM/HDGB profiles. As discussed previously,[52,56,60] CHARMM/HDGB is in good agreement with the explicit lipid results for most residues except for the charged amino acids where membrane deformations play a dominant role. Therefore, we focus the comparison here on reproducing the profiles obtained from CHARMM/HDGB with PRIMO-M.

In the past, amino acid side chain analogues were used in explicit simulations for evaluating transfer free energies. However, contributions of the peptide bonds should be included in the hydrophobicity scale for the obvious reason that whole residues, not just side chains, partition into membranes.[61] This means that the scales must be whole-residue scales. The advantage of choosing a TM helix instead of side chain analogues is that such an approach allows us not only to account for the effect of the protein environment but also to include the contribution of peptide bonds, which have been shown to influence the transfer free energy.[61] Another reason for adopting such an approach is the parametrization strategy used in optimizing the PRIMO force field. PRIMO was parametrized against the CHARMM force field by using a variety of short peptides, rather than amino acid side chain analogue molecules where the use of a CG model becomes problematic.[51]

The HDGB dielectric and nonpolar profiles are homogeneous in the x–y plane and affect only the sampling in the z direction. Taking this into consideration, a decoy set consisting of a transmembrane helix of type A49XA50 spanning the entire membrane is constructed for each residue. We have employed 100 configurations for each residue. The free energy of insertion (ΔGinsertion) was first calculated for all the structures with the CHARMM force field with default HDGB parameters. Later, this step was repeated for PRIMO-M. The free energy of a residue (X) for a given insertion depth is calculated by subtracting the free energies ΔGinsertion (A49XA50) from ΔGinsertion (A100) of an equivalent A100 at the same insertion depth.The free energy profiles were generated by placing the center of mass of the initial structure at the center of the membrane, and then translating it along the bilayer normal until membrane surface (0, 0, 30) in steps of 0.5 Å. The resulting profiles are shown in Figure 1 for both the CHARMM/HDGB and PRIMO-M methods. The profiles obtained with PRIMO-M are generally quite similar to the results obtained with the CHARMM/HDGB model. There are only three residues where PRIMO-M deviates significantly from CHARMM/HDGB. For phenylalanine (Phe), the free energy decreases by 1 kcal/mol upon membrane insertion with PRIMO-M vs a slight increase with CHARMM/HDGB. For histidine (HSD/HSE), PRIMO-M shows ∼3 kcal/mol less of a penalty to reside at the membrane center compared to the all-atom CHARMM/HDGB model, and for glycine (GLY), PRIMO-M is less favorable in the center of the membrane by about 0.6 kcal/mol. In all other cases, the agreement is exceptionally good. In general, polar residues, such as histidine, may be more challenging because of the reduced charges in PRIMO. However, it should be noted here that the dielectric and the nonpolar profiles were highly optimized for the all-atom CHARMM force field. A further optimization of these two profiles for the PRIMO force field may improve the overall agreement with CHARMM/HDGB. Parameterization of PRIMO-M could also improve the agreement with the all-atom model, especially for histidine, but we decided to proceed without modifying the underlying PRIMO model in order to preserve its overall balance and transferability.

Figure 1

Amino acid insertion free energy profiles with PRIMO-M (solid, orange) compared to results from CHARMM/HDGB (dashed, green) simulations.

Transfer Free Energies of Amino Acids Compared to Simulation and Experiment

Transfer free energies from water to the membrane were computed as the difference between the average solvation free energy at the center of membrane (z = 0 Å) and the average energy in the aqueous environment (z = 30 Å). We used 10 neighboring structures for each averaging. Figure 2 compares the computationally derived scales (CHARMM/HDGB and PRIMO-M) with experimental apparent free energies. The PRIMO-M and CHARMM/HDGB results are correlated remarkably well (r = 0.93), again with histidine and phenylalanine being the major outliers. The accuracy of the calculated transfer free energies exemplifies the transferability of the PRIMO model. PRIMO-M values also correlate well with the biological scale[62] (r = 0.93) and with transfer free energies of side chain analogues into cyclohexane[63] (r = 0.91). The CHARMM/HDGB model also yields similar correlation with the biological and cyclohexane scales (Figure S1, Supporting Information). Interestingly, PRIMO-M histidine transfer free energies agree better with the experimental scales than with CHARMM/HDGB, suggesting that the disagreement between PRIMO-M and CHARMM/HDGB may be in part due to uncertainties in the all-atom energetics. In other studies, Ulmschneider et al.[64] also obtained similar correlations for both experimental scales with their Monte Carlo based GB membrane model using the OPLS-AA force field.

Figure 2

Insertion energy of the designed peptides: CHARMM/HDGB versus PRIMO-M (A), experimental transfer free energies of side-chain analogues from water into cyclohexane versus PRIMO-M (B), and biological hydrophobicity scale versus PRIMO (C). The corresponding correlation coefficients (r) are also provided. Charged residues were excluded from the scale since their insertion free energies are considerably overestimated by both CHARMM/HDGB and PRIMO-M methodologies.

Transfer free energies for charged residues are given in Table 3. Both CHARMM/HDGB and PRIMO-M, while agreeing with each other, greatly overestimate the experimental values. This is a result of the fixed membrane geometry in HDGB and fixed ionization states that does not consider membrane deformations or the possibility of charge neutralization upon membrane insertion. Explicit lipid simulations have shown that upon burial in the hydrophobic membrane core, a charged residue is likely to be either neutralized or accompanied by a shell of water molecules, which will lower its insertion energy significantly.[64] Tieleman et al. calculated pKa values of ionizable groups in DOPC using an all-atom model.[65] Asp and Glu were found to be neutral in the lipid core, whereas Lys and Arg prefer the charged state.

Table 3

Bulk-Solvent-to-Membrane-Center Transfer Free Energies of Charged Amino Acids Integrated into an α-Helical Conformation Compared with the CHARMM/HDGB Model and the Biological Scalea

amino acid	ΔGtransferCHARMM (kcal/mol)	ΔGtransferPRIMO (kcal/mol)	ΔGtransferPRIMO/DHDGB (kcal/mol)	biological scale (kcal/mol)	water to cyclohexane (kcal/mol)
Arg	34.9	38.8	8.9	2.6	14.9
Asp	44.8	44.7	8.7	3.5	8.7
Glu	45.7	50.7	8.4	2.7	6.8
Lys	38.6	42.9	8.1	2.7	5.6

Experimental water to cyclohexane transfer free energies of side chain analogues are also shown. The transfer free energy is computed as the difference between the average solvation free energy at the center of membrane and the average energy in the aqueous environment (ε = 80).

Water defects effectively constitute membrane deformations that can be taken into account with a deformable implicit membrane model such as in the DHDGB model that was recently developed by us.[60] We see a considerable decrease in the transfer free energy of the charged residues when the DHDGB model replaces the HDGB model in PRIMO-M to allow dynamic membrane deformations. The resulting insertion profiles for Arg, Lys, Asp, and Glu are shown in Figure 3, and the corresponding transfer free energies are reported in Table 3. The PMFs for all of the charged residues continue to increase until the center of the membrane as observed in previous explicit membrane simulations.[65] Tieleman and co-workers have estimated a free energy barrier of ∼14 kcal/mol for the side chain analogue of Arg at the center of the DOPC bilayer, while we estimate a value of ∼9 kcal/mol, which is 5 kcal/mol lower. However, PRIMO-M/DHDGB results for Asp, Glu, and Lys are comparable to explicit simulation data. We obtained a value of 8.7, 8.4, and 8.1 kcal/mol for Asp, Glu, and Lys, respectively, while the corresponding values from the explicit simulations are 7.4, 5.1, and 4.7 kcal/mol, respectively. Furthermore, PRIMO-M/DHDGB results are also comparable to CHARMM/DHDGB results for amino acid side chain analogue insertion.[60] The agreement is quite reasonable considering that DHDGB was used “as-is” without further optimization. However, since DHDGB is optimized for the atomistic CHARMM force field and the current model is only valid for single spanning helical structures, this model was not employed for subsequent studies in this paper.

Figure 3

Insertion free energy profiles for charged amino acids (Arg, Asp, Glu, and Lys) obtained with PRIMO/DHDGB.

Structural Stability and Dynamic Properties of Membrane Proteins

PRIMO-M was designed with a major goal of running stable molecular dynamics simulations of arbitrary membrane-protein systems. In contrast to other coarse-grained protein models, PRIMO does not require any bias toward known secondary structures or other structural constraints to model a given protein system. A set of 13 membrane proteins with 148 to 661 amino acids (see Table 1) and different topologies was tested. All protein simulations in the implicit membrane environment were started from the experimental structures and simulated with blocked termini for 50 ns. We analyzed the simulations with respect to thermodynamic stability and dynamic properties in comparison with experiments.

Figure 4 shows the time evolution of Cα RMSD of all the membrane proteins with respect to their initial structures. It is evident from the figure that stable conformations for most of the proteins are reached within the first 10 ns and thereafter the Cα RMSDs are kept within 4.0 Å for most proteins. However, in the case of the β1-adrenergic receptor (PDB code 2VT4), the RMSD increases steadily after 20 ns and reaches equilibrium plateau at 30 ns, fluctuating around 4.3 Å for the last 20 ns of the simulation. Table 4 lists average Cα RMSD values for all proteins during the simulation as well as Cα RMSD of the average structure over the entire trajectory. It is evident from Table 4 that the RMSD of the average structure is lower than the average instantaneous RMSD values as it corresponds more closely to the experimental scenario. Therefore, we focus our discussions on those values. Table 4 shows that the RMSD of average structures varies between 2.2 and 4.5 Å. On the other hand, Tanizaki and Feig found the RMSD of the average structures varies between 1.7 to 2.3 Å for relatively short simulations (5–10 ns) with CHARMM/HDGB of bacteriorhodopsin and BtuCD.[55] Out of the 13 proteins, the RMSD is found to be between 2 and 4 Å for 11 systems and only two systems show RMSDs above 4 Å. The smallest RMSD (2.2 Å) is obtained for the intramembrane protease, 2IC8, while the largest RMSD (4.5 Å) is obtained for the outer membrane transporter protein (FecA), 1KMO. A similar drift in RMSD is also obtained for the outer membrane protease (ompT), 1I78. Interestingly, both FecA and ompT are beta barrels. The membrane proteins with channels, such as ompT and FecA, contain water molecules inside the channel, which is neglected because the HDGB model used here assumes the presence of lipids wherever there is no solute in the membrane region. This may cause destabilization of the protein resulting in larger RMSD. Furthermore, the presence of long, more mobile loops at the extracellular part may cause somewhat higher RMSD values. To circumvent this problem, the current HDGB formalism would need to be extended to allow different dielectric environments, not just perpendicular to the membrane bilayer but also along the bilayer. One possibility would be to limit the application of a varying dielectric only to atoms facing outward while atoms facing an internal cavity would be assumed to be in an aqueous phase. While it is straightforward to implement such an approach within HDGB, this would have to be based on structural knowledge rather than first principles. Therefore, we did not pursue this option here.

Figure 4

Time evolution of Cα RMSD during PRIMO MD simulations using the HDGB methodology for selected proteins from their respective crystal structure.

Table 4

Root Mean Square Deviations from Experimental Structures in PRIMO MD Simulations in Membrane (PRIMO-M) and Aqueous (PRIMO) Environmentsa

		PRIMO-M		PRIMO
PDB	number of residues	avg. Cα RMSD (Å)	Cα RMSD of avg. struct. (Å)	avg. Cα RMSD (Å)	Cα RMSD of avg. struct. (Å)
1QHJ	228	3.9 (0.6)	3.6	6.8 (1.0)	6.5
2BL2	156	3.4 (0.5)	3.2	3.5 (0.4)	3.3
4H1D	173	3.3 (0.3)	3.2	4.0 (0.4)	3.7
2CFQ	417	4.3 (0.4)	3.9	6.5 (0.8)	6.2
2JLN	463	4.2 (0.6)	3.9	5.9 (0.8)	5.7
1ZCD	376	3.6 (0.3)	3.4	5.3 (0.7)	5.1
3D9S	245	3.7 (0.5)	3.4	7.5 (1.1)	7.2
1K24	249	4.3 (0.6)	3.9	6.7 (1.0)	6.3
2VT4	276	3.8 (0.6)	3.4	5.6 (1.1)	5.3
1I78	297	5.1 (1.1)	4.4	6.2 (1.1)	5.7
1QJ8	148	3.4 (0.6)	3.2	4.1 (0.5)	3.9
1KMO	661	4.7 (0.7)	4.5	9.7 (1.8)	9.3
2IC8	182	2.5 (0.3)	2.2	4.1 (0.7)	3.7
avg.		3.9 (0.7)	3.6 (0.6)	5.9 (1.7)	5.5 (1.7)

Standard deviations are provided in parentheses.

We also conducted PRIMO-M MD simulations of four proteins, namely 4H1D, 2BL2, 2IC8, and 1QJ8 with larger cutoff distances of 36/38 Å. The time evolutions of Cα RMSD for these four membrane proteins with respect to their initial structures are shown in the Supporting Information (see Figure S2). We found that the RMSD of the average structure for each protein remains close to the value that is observed with the simulations with the shorter cutoff values of 16/18 Å.

The final structures obtained from our PRIMO-M MD simulations are superimposed with the corresponding experimental structure and are shown in Figure 5. In general, the structural variations for all proteins in PRIMO-M are small, and only minor rearrangements of loops and helices, most notably at the flexible N- or C-termini, are observed. Overall, these results suggest that, in general, PRIMO-M can maintain experimental structures of proteins well in the membrane environment. For all proteins, stable trajectories were generated that remained close to the starting experimental structures.

Figure 5

Superposition of the final structure (green) for selected proteins obtained from PRIMO-MD simulations onto the corresponding crystal structures (red).

In addition to thermodynamic stability, we also examine dynamic properties of the simulated membrane proteins. B-factors of Cα atoms were calculated from the PRIMO-M simulations and compared with data from the experiment. Figure 6 shows the result for the monomeric bacteriorhodopsin system (PDB: 1QHJ) while other data are shown in the Supporting Information (see Figure S3). Qualitatively, the B-factors obtained from the coarse-grained PRIMO-M MD simulations agree very well with the experimental data, nicely reproducing the alternation between rigid secondary structure elements and flexible loop regions. In one loop region between residues 191 and 200 of bacteriorhodopsin, the calculated values are substantially larger, suggesting extensive motion. This is because the loop is completely exposed to the solvent and undergoes a conformational transition during the simulation.

Figure 6

B-factors of the Cα atoms were calculated from their root-mean-square fluctuations. The red line is from the PRIMO-M simulation while the green line is from the PRIMO simulation in an aqueous environment. The blue line is from the experiment.

Effect of Implicit Membrane Environment

In order to see the effect of the implicit membrane environment, we have carried out another set of simulations for the above-mentioned proteins with just an implicit aqueous environment using the GBMV model as in the original PRIMO model. The simulation protocol was the same as described in our previous paper.[51] Average Cα RMSD values during the simulation as well as Cα RMSD of the average structure over the entire trajectory are reported in Table 4. It is evident that the inclusion of the membrane is essential in maintaining stable structures in the simulations. With the aqueous environment, the RMSD values are larger for all of the cases, often much larger, corresponding effectively to denaturation. The time evolution of Cα RMSD of all proteins in aqueous solution with respect to their initial structures are shown in the Supporting Information (see Figure S4).

Figure 6 shows the B-factors of Cα atoms calculated from the implicit aqueous solvent simulation for bacteriorhodopsin. The qualitative features of B-factors in the BRD7 aqueous solvent simulations are notably different from the experimental values. Furthermore, the qualitative agreement in the aqueous solvent simulation is not as good as in the PRIMO-M simulation. In particular, the first and second helices have relatively large B-factors compared to those of other helices. The trend in the B-factors indicates increased flexibility with the aqueous solvent environment relative to the result of the implicit membrane simulation. However, the region 191–200 shows similar B-factors to that in the HDGB simulation.

Tilting of Helical Peptides in a Membrane Bilayer

In the previous section, we showed that PRIMO-M is capable of maintaining the native structure of membrane proteins in MD simulations. Here, we evaluate the applicability of the force field in predicting the insertion angle of helical transmembrane (TM) peptides into the membrane bilayer.

Significant variations of the tilt angle are observed for different TM peptides in membrane bilayers due to hydrophobic mismatch, which is defined as the difference between the lengths of the hydrophobic core of the helix and the hydrophobic width of the membrane. The hydrophobic mismatch is thought to affect the orientations of TM helices in membrane proteins and therefore their structures and functions.[66] When the lipid hydrophobic thickness is shorter than the peptide hydrophobic length (positive mismatch), the peptide tilts from the membrane normal to decrease its hydrophobic length to allow it to have better interactions with the lipids and lesser exposure of hydrophobic residues to the solvent. However, smaller tilting is observed for peptides with the hydrophobic length smaller than the hydrophobic thickness of the lipid (negative mismatch).[67] Peptides may adopt a surface orientation if the mismatch is excessively negative.

In this work, we investigated the tilting of WALP23, WALP19, GWALP23, KWALP23, AChR M2, and NMDA M2 peptides, comparing results obtained with the CG model both to experimental data and to simulations performed with various atomistic and CG force fields. We note, though, that the exact tilting of the WALP peptides is still debated, as experimental and simulation results differ. The starting structures of the peptides were fully α-helical and oriented exactly parallel to the bilayer normal. The sequences of these peptides are provided in Table 2. The tilt angle was measured as the angle between the membrane bilayer normal and the principal axis of the helix. We noticed that WALP19 sometimes left the membrane and adopted an interfacial orientation. While we are not sure about the exact origin of that behavior, we applied a weak restraint along the z-axis to keep WALP19 within the membrane since the focus here was on determining tilt angles of the membrane-inserted peptides.

Tilt angles and their distributions were calculated and are shown in Figure 7. In all cases, there is a large variation in the tilt angle. Figure 8 depicts the conformation with the most probable insertion angles of these peptides obtained from our simulations. We did not see any kink or distortion in our simulated structures. The mean tilt angles obtained from PRIMO-M are compared with experimental and PACE results in Table 5. All peptides investigated have their center of geometry close to the core of the membrane (1.0–2.8 Å) and low tilt angles (8.6°–10.5°) as reported in Table 5.

Figure 7

Distributions of insertion angles of transmembrane helical peptides at 300 K obtained from REX–PRIMO-M simulations.

Figure 8

Dominant structures of helical peptides with most probable tilt angle.

Table 5

Average Tilt Angles of WALP23, WALP19, GWALP23, KWALP23, AChR, and NAMDR Compared to Experiment (DOPC) and PACE (DOPC)a

peptide	PRIMO (deg)	PACE[42] (deg)	experimental (deg)	avg. z (Å)
WALP23	9.1 (6.0)	7.5 (5.0)	11.0[73]	1.3 (1.0)
WALP19	10.5 (6.0)	8.0 (5.7)	4.0[103]	1.4 (0.9)
GWALP23	9.6 (5.0)	6.5 (5.3)	6.0[78]	1.4 (0.9)
KWALP23	10.3 (5.2)	7.5 (5.7)	7.3[78]	1.0 (0.7)
AChR	8.6 (4.6)	NAb	11.0[79]	1.7 (1.2)
NMDAR	10.1 (5.0)	NAb	NAb	2.8 (1.8)

The mean displacement of the center of mass with respect to the membrane center (z) is also reported. The membrane center is located at z = 0, and negative z values correspond to shifts toward the extracellular side, whereas positive z values correspond to shifts toward the cytoplasm. Standard deviations are provided in parentheses.

Not available.

The WALP model peptides were designed by Killian and co-workers to investigate the hydrophobic mismatch in different membrane environments.[68−70] Figure 7A compares the distribution of tilting of different WALP peptides. The calculated average tilt angles and fluctuations are 10.5° ± 6.0° and 9.1 ± 5.0° for WALP19 and WALP23, respectively. This shows that both peptides respond to this hydrophobic mismatch by tilting its helical access relative to the membrane normal.

Our result for WALP23 agrees very well with other CG simulations and experiments also. In DMPC, a tilt angle of about 12° was measured by ATR-FTIP spectroscopy.[68] Koeppe et al. found a tilt angle of 8.1° in DLPC using the GALA method.[71] On the contrary, a recent study by fluorescence spectroscopy yielded a large tilt angle for WALP23 in DOPC (23.6°).[72] Tilt angles of WALP23 in DMPC and DLPC have been found to be 7° and 15°, respectively, by the GALA method using the dynamical model.[73]

It is interesting to compare our result with other CG and all-atom results. In a hybrid molecular dynamics study by Wan et al., a moderate tilting in DOPC (7.5°) and DLPC (17.5°) was found.[42] A tilt angle of 14° was found for WALP23 in DPPC using the MARTINI CG force field and Bond and co-workers’ protein model,[74] while the MARTINI protein and lipid model showed a tilt angle of 11.4° in DOPC.[75] Kim and Im observed a tilting of 14.9° in explicit solvent simulations of the peptide in a POPC bilayer.[76] Finally, with the CHARMM/HDGB simulations, Panahi and Feig found that the most probable insertion angles are located at 10.8° and 16.3° with the smaller angle being slightly more favorable.[60] However, with CHARMM/DHDGB simulations, the second peak was shifted to 26.4°, and the first peak was observed at 9.2°. The bimodal distribution of the tilt angle in their simulations could be due to the insufficient sampling of the peptide as they had conducted 32 independent short simulations of 20 ns using the Nosé-Hoover thermostat.

We found a larger tilt angle for WALP19 (10.5° ± 6.0°) than the experimental result (4.0°) obtained from the GALA method without considering the dynamic motion of peptides. This apparent discrepancy can be solved through a different interpretation of the result obtained with the GALA method. However, our result for WALP19 agrees very well with what was determined with the PACE force field[42] (8.0° ± 5.0° in DOPC) and the explicit membrane simulations with the CHARMM22 force field[76] (12.1° in DMPC). The GBSW implicit membrane also predicted a larger tilt angle (15.5°) for WALP19 and compared fairly with our result.[77] On the other hand, Bond et al. predicted larger tilting for shorter peptides (22° for WALP19 versus 14° for WALP23).[74]

The effect of hydrophobic mismatch is also observed for another family of peptide, GWALP. The average tilt angles of GWALP23 and its mutant, KWALP23 (with Lys residue at the two ends), were found to be 9.6° ± 5.0° and 10.3 ± 5.2°, respectively. Our results again match the experimental findings very well.[73] Tilt angles of 6.0° and 7.3° were determined in an experiment for GWALP23 and KWALP23, respectively.[78] Wan et al. obtained a tilt angle of 6.5° ± 5.3° and 7.5° ± 5.7° for GWALP23 and KWALP23, respectively, with their hybrid force field, PACE.[42] Vostrikov et al. suggested that there is no significant difference in the tilt angle between KWALP23 with a Lys residue at the both ends and GWALP23.[78] The tilt angle difference is 1.3° in the experiment. We have also obtained a similar tilt angle difference, and our result supports this finding from the experiment.

Apart from studies of artificial transmembrane peptides, we have also investigated the tilting of M2 channel-lining segments from the nicotinic acetylcholine receptor (AChR, PDB ID: 1A11) and of a glutamate receptor of the NMDA subtype (NMDAR, PDB ID: 2NR1).[79] The distribution of tilt angles is shown in Figure 6B. The average tilt angles were found to be 8.6 ± 4.6° and 10.1 ± 5.0° for AChR and NMDAR, respectively. The mean tilt angle for AChR M2 obtained from our simulation compares well to the experimental value (11.0°)[79] as well as to the previously calculated value (11.0° ± 5.0°) for a generalized Born implicit membrane.[64]

In the case of NMDAR, solid-state NMR experiments point to an inserted TM orientation, but the exact tilt angle could not be determined.[79] In our simulation, the peptide partitions into the membrane bilayer and assumes a transmembrane orientation with a tilting angle of 10.1°. Our simulations suggest that the center of mass of this peptide is located at z = 2.8 Å from the membrane center, confirming the experimental finding of TM orientation. On the contrary, Ulmschneider et al. obtained parallel surface-bound orientations with their implicit membrane model.[80]

Conformational Sampling of Influenza Fusion Peptide

Significant effort has been put into determining the membrane bound structure of the influenza fusion peptide (IFP). In micelles, IFP was found to be an inverted V-like helix–break–helix configuration, in which both the N- and C-termini insert into the membrane.[81] Computer simulations have also been employed to study the configuration adopted by the IFP in a membrane environment, and a wide variety of results are obtained. Panahi and Feig predicted[57] a predominantly helical hairpin conformation of the native influenza fusion peptide in the membrane environment using CHARMM/HDGB simulations, which was later experimentally validated by Lorieau et al.[82] Similar structural insights were revealed by a recent solid-state NMR (ssNMR) experiment.[83]

Here, we have conducted a replica exchange molecular dynamics simulation study of the monomeric IFP (PDB ID: 1IBN) using an implicit environment description (HDGB) and coarse-grained representations of the peptide. Each replica was simulated for 100 ns, and the last 80 ns of data were used for the analysis (see Table 2). The insertion angle was calculated as the angle between membrane normal (z axis) and the vector that connects the N (i + 4) and O (i) atoms of the backbone for i from 3 to 6.[57] The N-terminal insertion angle and insertion depth (the distance of center of mass from the membrane center along the membrane normal) are determined at 300 K and their distributions are shown in Figure 9. The focus on the N-terminal part of the peptide was chosen to match computational results based on CHARMM22/HDGB.[57] Figure 9 indicates that most of the conformations adopted insertion angles between 80 and 100° and insertion depths between 19 and 21 Å. This means that peptides with parallel membrane orientations and an interfacial location dominate the conformational ensemble. Other less-populated similar structures have more obliquely inserted N-termini with angles of ∼110°.

Figure 9

Distribution of insertion angle (A) and insertion depth (B) of IFP at the lowest temperature 300 K.

It is interesting to compare our results with other data. The parallel membrane orientation matches results from previous simulations.[57,84,85] The heteronuclear triple resonance NMR study suggested that the influenza fusion peptide adopts a tight helical hairpin conformation at the membrane–water interface.[82] Furthermore, previous studies suggest that the structure of the influenza fusion peptide in membranes could be a kinked helix, a straight helix, and a tight shaped helical hairpin. Interestingly, in our PRIMO-M simulations, all of these three conformations were observed. These three states were also visited in the simulations conducted by Larsson and Kasson.[84] A total of ∼19% of structures were found to be helical hairpins in our simulations (Figure 10). Nearly 60% of the total structures were found to be predominantly flat α-helical, while 20% of the total conformations were of kinked helix. Overall, our results rather agree very well with findings from recent simulations[57,84] or experiments.[82,83]

Figure 10

Representative structures of IFP at the lowest temperature (300 K).

Helix–Helix Association of Glycophorin A

Oligomerization of transmembrane (TM) helices is a key stage in the folding of membrane proteins.[86] The two-stage folding model has been proposed for the association of TM helices: each TM helix is inserted into the membrane independently followed by the assembling of helices.[87] Here, we have investigated the association of two Glycophorin A (GpA) helices, which has served as an excellent model system for studying TM membrane protein structure and stability, from both an experimental[87−89] and a computational[29,46,90−97] perspective. The GpA dimer contains a seven-residue motif (L[75]I[76]xxG[79]V[80]xxG[83]V[84]xxT[87]) that has been found to be important in the packing and dimerization of GpA helices.[98] Previous CG simulations were conducted to study such processes,[86,99] and GpA is an obvious test case for PRIMO-M.

The single TM domain of GpA consists of ∼25 residues that adopt an α-helical conformation and are sufficient for dimerization. The sequence of the monomeric TM domain is show in Table 2. Two GpA helices were inserted into the implicit lipid bilayer in a parallel fashion with an interhelix separation of ∼30 Å. We performed a REMD simulation with 16 replicas spanning a temperature range of 270–500 K. Each replica was simulated for 200 ns. The final 100 ns were used for the analysis of the structural features of the GpA dimer.

A well-defined structural feature of the GpA dimer is the crossing angle (Ω) between the two helices and is found to be right-handed with a negative Ω. The experimental crossing angles were reported to be −40° and −35° for the DPC micelle and the DMPC bilayer, respectively.[89,100] The extent of crossing angle and its distribution were calculated and are shown in Figure 11. A large variation in Ω compared to atomistic simulation[90,91] was seen in Figure 11. This reflects the soft nature of the coarse-grained free energy landscape but may also reflect the nature of the HDGB implicit membrane model. In our simulations, the probability distribution of the crossing angle shows a symmetric bimodal distribution. The most probable crossing angles were located at ±42°, and they are almost equally favorable. This result compares well with the previous study by Bu et al.[101] using the GBSW membrane model. They found that the most probable crossing angles were located at −50° and 50° with −50° being more favorable. The MARTINI protein force field[102] was also found to yield a large distribution of Ω. Previous CG simulations reported an average value of −20° to −25°, which is lower than the experimental values by 15°–20°. A positive crossing angle as observed in PRIMO was sampled with the MARTINI lipid force field and Bond and Sansom’s protein model.[99] Han et al. found that the average crossing angles vary between −32.5° and −43.5° with their hybrid coarse-grained force field.[42]

Figure 11

Distribution of the crossing angle of GpA dimer at 300 K obtained from the PRIMO-REMD simulation.

Figure 12A and B show the distribution of Cα RMSD of the simulated dimer from the solution NMR structure in a DPC micelle and interhelix separation at the lowest temperature (300 K) during the REMD simulation. The RMSD and the interhelical distance are well correlated with the interhelical crossing angle, Ω. On the basis of the distribution of crossing angles, the GpA dimer could be clustered into two distinct families of conformations: a right-handed dimer (Ω at −42°) and a left-handed dimer (Ω at 42°). The right-handed dimer has a most probable RMSD and interhelical distance value of 4.2 and 8.2 Å, respectively, whereas the most probable RMSD and the interhelical distance value of the left-handed dimer are 5.6 and 9.5 Å, respectively. The interhelical distance between two helices for the right-handed dimer conformation matches very well to an explicit solvent simulation of the GpA dimer in the DPPC bilayer, the membrane modeled here with HDGB (8.2 Å in PRIMO versus 8.1 Å in DPPC).[90] Sengupta and Marrink[102] also found in their CG-MD simulations that the optimum interhelical distance was 7.5 Å, along with a second population with an interhelical distance of 9.5 Å. The average Cα RMSD was measured to be 3.6 Å in their simulations. With the PACE force field, an interhelical distance of ∼10 Å was obtained.[42] A similar result was obtained by Psachoulia et al.[99] with the MARTINI force field and Bond–Sansom protein model. To further evaluate the structural features of the simulated GpA dimer, the contact map was calculated and shown in the Supporting Information (Figure S5). It is evident from the figure, barring a few, that all the native contacts present in the crystal structure tend to form in our simulations. The PRIMO simulations also suggest, as in previous studies, that the interactions are dominated by the key residues of the (L[75]I[76]xxG[79]V[80]xxG[83]V[84]xxT[87]) motif. Overall, our result overlaps with experiments and compares well with other computational studies.[42,99]

Figure 12

Distributions of Cα RMSD relative to the NMR structure in DPC micelle (A) and interhelix separation distance (B) at 300 K obtained from the REX–PRIMO-M simulation.

Efficiency of PRIMO-M Compared to CHARMM/HDGB

For a CG model to be attractive, computational efficiency is essential. Given its rather high resolution, it could be expected that the PRIMO-M methodology is not extremely fast compared to atomistic simulations with the HDGB methodology. In order to compare the computational efficiency, all of the above-mentioned 13 proteins were simulated for 1 ns with PRIMO and CHARMM using HDGB methodology. In the all-atom simulations, the CHARMM36 force field parameters were employed with the CMAP correction term. The bond lengths involving hydrogen atoms were fixed by using the SHAKE algorithm, so that a simulation time step of 2 fs could be used. Furthermore, both PRIMO-M and CHARMM/HDGB simulations were carried out with cutoff distances of 16/18 Å and without cutoff distances. All the PRIMO-M and CHARMM/HDGB simulations were performed in serial on an Intel E5–2680 processor (2.7 GHz). The CHARMM/HDGB simulations were carried out with a time step of 2 fs, while a time-step of 4 fs was used in the case of PRIMO-M. Table 6 lists the corresponding simulation time. In both cases, the simulation time is proportional to the system size. Compared to CHARMM/HDGB simulations, PRIMO-M can achieve about 10- to 20-fold speedups provided the simulations are carried out without any nonbonded cutoff distances. However, about 9–13 speedups are achieved with PRIMO-M with respect to the CHARMM/HDGB simulations when the simulations are conducted with the nonbonded cutoff distances of 16/18 Å. A similar speedup is achieved over all-atom explicit lipid/water simulations. A main bottleneck in PRIMO is the use of the GBMV methodology for modeling the membrane environment. Nearly 80% of the total simulation time is spent for GB calculations. One possibility for further accelerating the PRIMO-M model is to replace the GBMV-based implicit membrane model with a computationally more efficient GB implementation.

Table 6

Efficiency of Simulations of 13 Proteins with Two Different Simulation Methodologiesa

				timing (ns/day)
		sites		no cutoff		cutoff (16/18 Å)		PRIMO speedup
PDB	Res	CHARMM	PRIMO	CHARMM	PRIMO	CHARMM	PRIMO	no cutoff	cutoff
1QHJ	228	3599	1207	0.15	2.49	0.26	3.28	16.6	12.6
2BL2	156	2319	784	0.31	4.63	0.46	4.95	14.9	10.8
4H1D	173	2798	928	0.23	3.28	0.38	4.05	14.3	10.7
2CFQ	417	6632	2204	0.05	1.00	0.14	1.40	20.0	10.0
2JLN	463	7245	2420	0.05	0.85	0.12	1.03	17.0	8.6
1ZCD	376	5814	1940	0.07	1.21	0.15	1.47	17.3	9.8
3D9S	245	3712	1246	0.14	2.35	0.27	3.18	16.8	11.8
1K24	249	3926	1354	0.13	1.85	0.23	2.47	14.2	10.7
2VT4	276	4480	1493	0.18	1.74	0.21	2.34	9.7	11.1
1I78	297	4591	1617	0.10	1.42	0.20	2.12	14.2	11.4
1QJ8	148	2235	787	0.31	4.19	0.46	4.69	13.5	10.2
1KMO	661	10 094	3549	0.03	0.42	0.08	0.91	14.0	11.4
2IC8	182	2919	973	0.22	3.30	0.33	3.74	15.0	11.3

The simulations were carried out with or without the nonbonded cutoff distances. CPU times are provided in ns/day.

Discussions and Conclusion

Many biologically interesting phenomena, such as the dynamics of large proteins and self-assembly of biomolecules that occur on a time scale that is too long to be studied by fully atomistic simulations. Coarse-graining can drastically cut down the necessary simulation times. Further acceleration is obtained with implicit models of the environment, especially for membrane environments. In the present paper, we are demonstrating the performance of such a model, PRIMO-M, in the context of a number of test cases. The amino acid insertion free energy profiles obtained with the PRIMO-M model are very similar to those obtained with the CHARMM/HDGB model. The free energy of insertion for amino acids is highly correlated to the biological and water-to-cyclohexane scales. It was shown that the model could be applied successfully to obtain stable and dynamically well-behaved trajectories of membrane proteins. Simulations of membrane proteins of varying complexity remained close to the starting X-ray structure after 50 ns of simulation time, while B-factors calculated from the simulations are in good agreement with the experiment. The force field is further validated by correctly predicting the tilt angle of several transmembrane peptides, association of the GpA homo dimer, and conformational sampling of an influenza fusion peptide (IFP) using replica exchange molecular dynamics simulations. Our simulations reproduce related experimental or theoretical observations quite well, which implies that the environment in PRIMO is interchangeable between aqueous and membrane environments.

PRIMO-M was obtained by simply swapping the GBMV implicit solvent model for the HDGB implicit membrane model without any reparameterization of the underlying PRIMO model. This implies that further advances in the implicit membrane model would also directly benefit PRIMO-M. In particular, membrane proteins with internal cavities or channels remain a challenge that needs to be addressed. Another issue is a better description of the nonpolar interaction that can be addressed by separating the cost of cavity formation from solute–membrane van der Waals interactions. PRIMO (and PRIMO-M) is also especially attractive in the context of AA/CG schemes because of the close structural and energetic correspondence between PRIMO and fully atomistic models, and exploring such models to optimize the balance between accuracy and computational efficiency will be the subject of further studies.

PRIMO-M relies on PRIMO and HDGB, both of which are implemented in CHARMM c38a2 and newer versions. The PRIMO force field files are available from the authors upon request.