membrane 
 proteins play significant roles in many cellular processes , 
 such as intercellular communications , molecular transports , and signal 
 transductions . 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 . hence ,
understanding their structures , dynamics , and functions are 
 of the utmost importance in the biological and pharmacological sciences .
however , a detailed study of membrane proteins 
 is challenging . because of the complexity of lipid bilayers ,
it has 
 proved difficult to map details of protein  membrane interactions 
 using experimental techniques .
molecular 
 dynamics simulation can serve as a complementary tool to illustrate 
 the structural and dynamic details of proteins at the atomistic level .
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 ; bilayer structure and dynamics ; and protein interactions with lipid bilayers .
the most popular and widely used coarse - grained model 
 for membrane - bound proteins and peptides is martini , 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 .
following the mapping scheme of martini , dal peraro and co - workers 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 
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 , for example 
 with united - atom protein models interacting with martini cg lipids 
 and solvent .
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 . 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 . 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 .
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 ( 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 .
since the primo 
 force field and its parameters are described in detail in our previous 
 paper , 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 . 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.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 12 ( bond ) 
 and 13 ( 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 
 14 terms , primo also uses two - dimensional spline - interpolated 
 cmap potentials to couple the sampling 
 of 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 .
the hdgb model is based 
 on the gbmv method , and as the details of the hdgb model have been 
 described elsewhere , 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 gsolv and the nonpolar contribution gsolv.2 in the hdgb model , the following modified 
 expression for the polar component of the solvation free energy is 
 used to introduce a heterogeneous dielectric environment3where qi are charges 
 of cg beads , n is the total number 
 of cg particles , rij are 
 distances between particles i and j , i 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 and subsequently optimized 
 by sayadi and feig 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:4where ai is the solvent accessible surface area 
 of the ith atom , zi is the distance 
 of atoms i from the membrane center along the membrane 
 normal , i is the coefficient of 
 surface tension for different atom types , and s(zi ) is a switching function 
 that is used to reflect the change of the surface tension along the 
 membrane normal .
the switching function s(zi ) 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 . 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 .
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 . in this study , periodic 
 dielectric and nonpolar profiles were used 
 as described previously .
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 using gb parameters described in detail previously .
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 was employed to maintain a temperature 
 of 300 k. the time step for leapfrog verlet integrator
nonbonded interactions were cut off at 18  , with smooth 
 switching to zero starting at 16  .
a similar protocol was followed 
 for simulations with primo in an implicit aqueous environment .
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 . 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 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 6065% .
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 where the primo model is implemented .
the mmtsb ( multiscale modeling 
 tools for structural biology ) tool set in combination with charmm was employed for all analyses .
the primo - m model was implemented 
 as described above and applied 
 to a number of membrane - related test cases as detailed in the following .
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 , 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 .
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 .
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 .
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 .
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.5the 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 .
amino acid 
 insertion free energy profiles with primo - m ( solid , 
 orange ) compared to results from charmm / hdgb ( dashed , green ) simulations .
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  ) .
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 ( r = 0.93 ) and with transfer 
 free energies of side chain analogues into cyclohexane ( 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 . also obtained similar 
 correlations for both experimental scales with their monte carlo based 
 gb membrane model using the opls - aa force field .
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 ) .
charged residues were excluded 
 from the scale since their insertion free energies are considerably 
 overestimated by both charmm / hdgb and primo - m methodologies .
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 .
calculated pka values of ionizable groups in dopc using an all - atom 
 model .
asp and glu were found to be neutral 
 in the lipid core , whereas lys and arg prefer the charged state .
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 .
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 .
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 .
primo - m / dhdgb results are also comparable to charmm / dhdgb results 
 for amino acid side chain analogue insertion .
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 . insertion free energy 
 profiles for charged amino acids ( arg , asp , 
 glu , and lys ) obtained with primo / dhdgb .
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 .
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 ( 510 ns ) with charmm / hdgb 
 of bacteriorhodopsin and btucd . 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 .
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 .
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 .
time evolution of c rmsd 
 during primo md simulations 
 using the hdgb methodology for selected proteins from their respective 
 crystal structure .
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 ,
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 ,
this is because the loop is completely 
 exposed to the solvent and undergoes a conformational transition during 
 the simulation .
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 . 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 .
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 191200 shows 
 similar b - factors to that in the hdgb simulation . 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 .
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 ) .
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 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 .
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.02.8  ) 
 and low tilt angles ( 8.610.5 ) as reported in table 5 .
distributions of insertion angles of transmembrane helical 
 peptides 
 at 300 k obtained from rex  primo - m simulations .
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 .
the walp 
 model peptides were designed by killian and co - workers 
 to investigate the hydrophobic mismatch in different membrane environments .
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 ,
koeppe 
 et al . found a tilt angle of 8.1 in dlpc using the gala method . on the contrary , a recent study by fluorescence 
 spectroscopy yielded a large tilt angle for walp23 in dopc ( 23.6 ) .
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 .
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 . a tilt angle of 14
was found for walp23 
 in dppc using the martini cg force field and bond and co - workers 
 protein model , while the martini protein 
 and lipid model showed a tilt angle of 11.4 in dopc .
kim and i m observed a tilting of 14.9 
 in explicit solvent simulations of the peptide in a popc bilayer .
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 .
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 ( 8.0  5.0 in dopc ) and the 
 explicit membrane simulations with the charmm22 force field ( 12.1 in dmpc ) .
the gbsw implicit membrane 
 also predicted a larger tilt angle ( 15.5 ) for walp19 and compared 
 fairly with our result . on the other 
 hand , bond et al . predicted larger tilting for shorter peptides ( 22 
 for walp19 versus 14 for walp23 ) .
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 . 

tilt angles of 6.0 and 7.3 were 
 determined in an experiment for gwalp23 and kwalp23 , respectively .
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 .
suggested that there is no significant 
 difference in the tilt angle between kwalp23 with a lys residue at 
 the both ends and gwalp23 .
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 i d : 1a11 ) and of a glutamate 
 receptor of the nmda subtype ( nmdar , pdb i d : 2nr1 ) .
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 ) as well as to the previously calculated value ( 11.0  
 5.0 ) for a generalized born implicit membrane . in the case of nmdar ,
solid - state nmr experiments 
 point to an inserted 
 tm orientation , but the exact tilt angle could not be determined . 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 .
significant effort has been put into determining 
 the membrane bound structure of the influenza fusion peptide ( ifp ) . 
 in micelles ,
helix 
 configuration , in which both the n- and c - termini insert into the 
 membrane .
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 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 .
similar structural insights were revealed by a recent solid - state 
 nmr ( ssnmr ) experiment . here , we 
 have conducted a replica exchange molecular dynamics simulation study 
 of the monomeric ifp ( pdb i d : 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 .
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 .
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. 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 heteronuclear triple resonance 
 nmr study suggested that the influenza fusion peptide adopts a tight 
 helical hairpin conformation at the membrane  water interface .
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 . 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 or experiments .
oligomerization of transmembrane ( tm ) helices is a key 
 stage in the folding of membrane proteins .
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 . 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 and a computational perspective .
the gpa dimer contains a seven - residue motif ( lixxgvxxgvxxt ) that has been found to be important in the 
 packing and dimerization of gpa helices .
previous cg simulations were conducted to study such processes , 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 .
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 270500 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 .
the extent of crossing 
 angle and its distribution were calculated and are shown in figure 11 .
a large variation in  compared to atomistic 
 simulation 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 
 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 . 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 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 
 1520. a positive crossing angle as observed in 
 primo was sampled with the martini lipid force field and bond and 
 sansom s protein model .
found that the average crossing angles vary between 32.5 
 and 43.5 with their hybrid coarse - grained force field .
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 ) .
sengupta and marrink 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 .
a similar result was obtained by psachoulia 
 et al . 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 ( lixxgvxxgvxxt ) motif .
overall , our result overlaps with experiments 
 and compares well with other computational studies .
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 . 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 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 e52680 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 .
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 
 913 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 .
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 .
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 . 

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 .
 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 
 148661 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
Model
Results and Discussion
Discussions and Conclusion