J. Phys. Chem. B 2008, 112, 11137–11146
11137
Systematic Comparison of Empirical Forcefields for Molecular Dynamic Simulation of Insulin Nevena Todorova,† F. Sue Legge,† Herbert Treutlein,‡ and Irene Yarovsky*,† Applied Physics, School of Applied Sciences, RMIT UniVersity, GPO Box 2476V, Melbourne, Victoria, 3001, Australia, and Cytopia Research Pty. Ltd., PO Box 6492, St. Kilda Road Central, Melbourne, Victoria, 8008, Australia ReceiVed: August 26, 2007; ReVised Manuscript ReceiVed: June 11, 2008
The use of atomistic simulation methodologies based on empirical forcefields has enhanced our understanding of many physical processes governing protein structure and dynamics. However, the forcefields used in classical modeling studies are often designed for a particular class of proteins and rely on continuous improvement and validation by comparison of simulations with experimental data. We present a comprehensive comparison of five popular forcefields for simulating insulin. The effect of each forcefield on the conformational evolution and structural properties of the peptide is analyzed in detail and compared with available experimental results. In this study we observed that different forcefields favor different structural trends. However, the all-atom forcefield CHARMM27 and the united-atom forcefield GROMOS 43A1 delivered the best representation of the experimentally observed dynamic behavior of chain B of insulin. Introduction In the last few decades molecular dynamics (MD) simulations have emerged as a powerful tool for the characterization of biomolecular structure and dynamics at the atomic level. This notable technique has helped us understand complex molecular processes associated with protein conformational changes, ranging from studies of enzyme-reaction mechanisms and ligand binding to problems of protein refolding and denaturation (see references within ref 1). Fundamental to such simulations are the forces that govern the atomic motions, derived from atomatom interactions usually given in the form of empirical potential energy function or a so-called forcefield. One functional form of such forcefield can be represented as follows:
ki ki Vn (li - li,0)2 + Σ (θi - θi,0)2 + Σ (1 + bonds 2 angles 2 dihedrals 2 qiqj σij 6 σij 12 N N + (1) cos(nω - δ)) + Σ Σ 4εij i)1j)i+1 rij rij 4πε0rij
V(rN) ) Σ
( [( ) ( ) ]
)
where V(rN) denotes the potential energy, which is a function of the positions (r) of N particles. The first three terms in Equation 1, model the bonded or intramolecular interactions, where the bonds and angles are represented by a harmonic potential and the dihedrals by torsional potential. The forth contribution is the nonbonded term, which in a simple forcefield is usually modeled using a Coulomb potential term for the electrostatic interactions and a Lennard-Jones potential for the van der Waals interactions. The choice of mathematical function describing a forcefield is important, since it will ultimately determine the quality of the results. The in-depth theory and history of development of different forcefields will not be covered in great detail in this article. The interested reader is referred to the papers by Ponder and Case1 and MacKerell2 for more information. * Corresponding author. Telephone: +61 3 9925 2571. Fax: +61 3 9925 5290. E-mail address:
[email protected]. † Applied Physics, School of Applied Sciences, RMIT University. ‡ Cytopia Research Pty. Ltd.
The biomolecular forcefields selected for our study include the recent versions of the all-atom CHARMM27,3 AMBER03,4 OPLS,5 and the united-atom GROMOS 43A16 and GROMOS 53A67 forcefields. The parameters for all five forcefields were extensively optimized with particular emphasis on the treatment of proteins. There are several variations available of the AMBER forcefield. Duan et al. 4 extensively modified the previously available AMBER948 and AMBER99,9 now called AMBER03, in which the derivation of the partial atomic charges was fundamentally different. Instead of relying on the Hartree-Fock/ 6-31G* approach to provide aqueous-phase charges, a lowdielectric continuum model corresponding to an organic solvent environment was directly included in the quantum mechanical (QM) calculation of the dihedral parameters and electrostatic potential, from which the charges are obtained. AMBER03 was selected for this study, because it showed the best agreement with experimental data in an extensive forcefield investigation between the readily available AMBER forcefields,10 not considering some custom modified AMBER versions. In OPLS and CHARMM27 forcefield the partial charges are based on HF/ 6-31G* supramolecular data.3,5 In contrast to the parametrization of other biomolecular forcefields, the parametrization of the united-atom GROMOS 53A6 forcefield is based primarily on reducing the free enthalpies of hydration and apolar solvation for a range of compounds. An essential aspect of the use of a forcefield for biomolecular simulations is the accurate treatment of the condensed aqueous environment. When selecting a water model to use for a particular study, the most important consideration is its compatibility with the biomolecular forcefield used. The reason for this is that most forcefields have been developed in conjunction with a specific water model. For example, AMBER, CHARMM and OPLS have been developed with the TIP3P water model, OPLS also with TIP4P while GROMOS with SPC and F3C models. Hess and van der Vegt11 performed an extensive study of hydration thermodynamic properties on several side chains, using AMBER99, GROMOS 53A6, and OPLS, with five water models, SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew (see
10.1021/jp076825d CCC: $40.75 2008 American Chemical Society Published on Web 08/12/2008
11138 J. Phys. Chem. B, Vol. 112, No. 35, 2008 references within ref 11). They found differences in hydration free energies between the water models to be small, however with AMBER99 and OPLS, TIP3P performs slightly better than the other water models. The GROMOS 53A6 forcefield performs stronger than the other forcefields for the free energy, which is not surprising, as it was parametrized on this property. Recent simulations12 of three proteins, performed using the wellknown AMBER94, CHARMM, and OPLS forcefields, concluded that these forcefields behave comparably or at least that it was not possible to distinguish between the forcefields in the 2 ns time scale investigated. In a recent article by Villa et al.,13 the sensitivity of molecular dynamics simulation to different forcefields was investigated. The three parameter sets, 43A1, 53A5, and 53A6 of the GROMOS forcefield were used to simulate 36 structures of 31 proteins. In their investigation, no major differences were detected in a wide range of structural properties such as the root-mean square deviation from the experimental structure, radii of gyration, solvent accessible surface, secondary structure, or hydrogen bond propensities on a 5-10 ns time scale, despite the differences in the forcefield parameters. On the contrary, Mu et al.14 showed that the recent version of the AMBER, CHARMM, GROMOS, and OPLS forcefields used in their study of the structure of trialanine differ considerably in the description of the dynamical properties of the system, where it was found that the lifetime of some conformational states differ by more then an order of magnitude, depending on which model was used. Several empirical forcefields and quantum mechanics/molecular mechanics (QM/MM) forcefields were used by Hu et al.15 to investigate the sampled conformations of unfolded polypeptide chains in aqueous solution. This group also found insufficient accuracy in their results and the conclusion followed from the widespread of results with different forcefields, rather than from a direct comparison with experimental data. The findings by Sorin and Pande16 on helix-coil systems in conjunction with Zeman et al. work on alanine dimers and trimers17 suggest that a dependence of forcefield accuracy on peptide length exists even for simple models such as polyalanine based systems. Discrepancies in simulation results due to the choice of forcefield have been well documented for a significant number of systems.18 The protein we selected for the present forcefield investigation is insulin. Insulin is a widely studied protein and there is an extensive experimental and computational literature available for comparison. The protein insulin is composed of two chains, A and B, containing 21 and 30 amino acids, respectively, and linked together by disulfide bonds. There are many structures currently available of insulin and the general binding mode of the insulin-receptor complex is known.19 Chain B of insulin is believed to retain much of its structure independently of chain A20-22 and is the subject of our investigation. Structure activity studies of insulin indicate the C-terminus of chain B is integral to receptor interaction,23-25 and also suggest that the conformation of the C-terminus is influenced by the structure of the N-terminus.26-28 Various studies of insulin,29,30 including a preliminary crystallographic structure of the native insulin monomer at a low pH,31 and also solution structures of isolated chain B determined by NMR spectroscopy,20,22 confirm both termini’s mobility. The secondary structure characteristics of chain B of insulin are generally defined as the N-terminal (residues 1 to 8), a central R-helix (residues 9 to 19), a characteristic chain B fold (β-turn, residues 20 to 23) and the extended C-terminus (residues 24 to 30). There are several known conformational states for the N-terminus of chain B.26,27,32-34 The principal conformations
Todorova et al.
Figure 1. Overlay of X-ray structures, showing experimentally observed states of insulin Chain B.
have been designated as the R-state and T-state,35 although additional conformations have also been identified.36 A schematic showing several conformational states is presented in Figure 1. The T-state is associated with insulin’s activity as it is believed to be the monomeric solution state.37 In the T-state, chain B consists of the R-helix (9 to 19) and the extended Nand C-termini regions.38,39 In the R-state, all N-terminal residues form a helix, joining with the central helix.40,41 Another known variation is the “frayed” or Rf-state, where the R-helix is only present for some of the N-terminal residues (4 to 8) in addition to the central helix. Molecular dynamics enables sampling of conformational states of a protein under controlled conditions and therefore delivers additional understanding of the proteins dynamics. Previous MD simulations using the GROMOS 37c forcefield, show conformational flexibility of insulin chain B, reporting a high degree of movement in aqueous solution of both the monomer and dimer.42 Simulations of 5-10 ns length using the CHARMM22 forcefield were performed by Zoete et al.43 of monomeric insulin in the T-state conformation, described the flexibility of the N- and C-termini regions of chain B. This inherent flexibility was observed in our previous computational study of possible thermal and chemical effects on the dynamics of chain B44 using the NAMD simulation package and the CHARMM27 parameter set. The effects of static and oscillating electric field on the flexibility of the chain was also investigated.45 The most recent of these studies performed by Legge et al.,46 yielded information clarifying uncertainties about the structure and dynamics of insulin with respect to its biological behavior by performing multiple MD simulations with CHARMM27 forcefield in explicit solvent. In this article we attempt to identify and compare possible systematic effects of multiple forcefields on the structure and dynamics of chain B of insulin in a course of identical molecular dynamics simulations. We also aim to investigate which forcefield gives the best representation of the experimentally observed conformational features and behavior of chain B within our selected time frame. The assessment of the forcefields has been performed by comparison of results with experimental data which illustrates the 3-dimensional morphology of the protein, such as X-ray crystallographic structures,35,36,40 nuclear Overhauser enhancement (NOE) distance restraints from the solution NMR structure of isolated chain B of insulin,20 and previous molecular dynamics study on chain B which reproduced known experimental structural features.46 To the best of our knowledge this is the only information available in the literature for direct comparison of conformational features observed in the simulations with different forcefields. Dynamic and structural properties such as the secondary structure evolution, solvent accessible
Comparison of Empirical Forcefields
J. Phys. Chem. B, Vol. 112, No. 35, 2008 11139
Figure 2. Secondary structure evolution of insulin chain B. Simulations for each forcefield are labeled as described in the Computational Methods section. Secondary structure color codes: magenta, R helix; red, π helix; cyan, turn; white, coil; yellow, extended conformation; green, bridge.
surface area (SASA), radius of gyration, interproton distance violations, have been obtained from each forcefield simulation for comparison. Computational Methods This study utilizes a classical molecular dynamics algorithm47 as implemented in the GROMACS 3.348 simulation package. To investigate the effect of different potential energy functions
on the atomic interactions the widely used all-atom CHARMM27, OPLS, AMBER03, and the united-atom GROMOS 43A1 and GROMOS 53A6 forcefields were selected. In this article they will be referred to as, CHARMM, OPLS, AMBER, GROMOSA1, and GROMOSA6, respectively. For each forcefield, two separate simulations were performed, for 50 ns each, starting from the same initial structure. All the simulations are labeled individually, based on the forcefield used and the number of
11140 J. Phys. Chem. B, Vol. 112, No. 35, 2008 the run, such as, for the calculations using AMBER, the simulations are labeled a1 and a2, with CHARMM, b1 and b2, with OPLS, c1 and c2, with GROMOSA1, d1 and d2 and with GROMOSA6, e1 and e2. Each system was immersed in explicit solvent (water), for which the model was selected according to that chosen by the developers at the time of parametrization. Specifically, for the CHARMM, AMBER and OPLS forcefields, the TIP3P water model was used, while for the united atom GROMOS forcefields we employed the SPC model. The starting structure was the Rf conformational state obtained from the crystal structure of porcine insulin (PDB code 1ZNI40). The His residues were protonated, and both termini regions and titrable sidechains were charged, resulting in a neutral to acidic pH environment. The protein was enclosed in a periodic box of 60 Å × 60 Å × 60 Å size, then solvated with ∼7075 TIP3P or SPC water molecules (depending on the forcefield chosen), corresponding to water density of ∼1.0 g/cm3. The whole system was energy minimized to remove steric clashes using the steepest descent algorithm. Two equilibration stages were then performed, each lasting 600 ps. First, equilibration was performed in the constant volume (NVT) ensemble, and then followed by the constant pressure (NPT) ensemble simulation at 1 bar. Constant pressure was achieved by coupling the system to a Parrinello-Rahman barostat49 with a relaxation time of 4 ps. The data collection stage was performed at the constant pressure (NPT) for 50 ns for each forcefield simulation. Constant temperature of 300 K was maintained using the Berendsen thermostat.50 The Particle Mesh Ewald51 (PME) summation was applied to treat long-range electrostatic interactions. Atom based cutoff of 12 Å with switching at 10 Å was used for nonbonded van der Waals interactions. All bond lengths were constrained to their equilibrium value using the SHAKE52 algorithm. A recent study by Schulten et al.53 highlighted today’s difficulties in obtaining sufficient sampling for studying folding using classical MD, and the possible effects forcefield might have on obtaining the correct folded structure. True thermodynamic equilibrium is not feasible from simulations such as these, however adequate equillibration can be demonstrated by monitoring the statistical deviations of the energy of the system. We monitored this property by calculating the deviations in potential energy from the mean value for each simulation. No systematic drifts or fluctuations higher than 3% were observed, signifying that the system was properly equilibrated for the statistical data collection.
Todorova et al.
Figure 3. The average rmsd structure of the two most populated clusters of each forcefield simulation, labeled as described in the Computational Methods section. The coloring of each protein is by its residue type: blue, basic; red, acidic; green, polar; silver, nonpolar. Hydrogen bonds of structure d2 are colored purple.
Results and Discussions Secondary Structure Evolution. The secondary structure dynamics reveal details of the conformational changes experienced by the protein over the simulation period. The evolution of the secondary structure of insulin chain B as a function of time for each forcefield system is presented in Figure 2. The secondary structure was classified by VMD using the algorithm STRIDE which utilizes hydrogen bond energy and mainchain dihedral angles in addition to hydrogen bond distances.54 The average rmsd structure of the two most populated clusters of each simulation were superimposed over the R-helix region, residues 9 to 19, and color coded by their residue type in Figure 3. There is a variation of secondary structures featured in different forcefield simulations, however, some similar trends can be observed. All systems showed an early change from the starting structure at residue Cys 19, most evident in the AMBER (Figure 2a) and GROMOSA1 (Figure 2d) simulations. This change can be expected due to the loss of the disulfide bonds
with chain A and the associated destabilization around this residue. Some helical content was maintained in each simulation although with great variation, as the close inspection of Figure 2 demonstrates. Specifically, the AMBER forcefield simulations preserved the most helical content, the helix ranging from Asn 3 to Cys 19, maintaining the majority of structural elements of the initial conformation (Figure 2a). Slight helical extension was seen after 11 ns in simulation a1, between Glu 13 and Val 18, however within the next 13 ns the helix was fully restored and stayed stable for the rest of the simulation, demonstrating the AMBER forcefield tendency to “favor” the helical conformation.10 The turn region between residues 19 and 22 was largely maintained throughout the simulation. Interestingly, this type of secondary structure was observed in a modeling study using CHARMM27 forcefield45 where chain B was simulated in a static electric field of 108 V/m, which stabilized the helix by aligning the peptide along the field direction. AMBER’s tendency to overstabilize
Comparison of Empirical Forcefields helical conformers has been observed in other studies.10 The two most populated clusters in our simulations (Figure 3a) confirm this preference, showing very little variance in the helix with a slight indication of mobility in the C-terminus. It should be noted that the flexibility of the carboxy terminus of chain B was proposed to be an important feature in receptor binding33,55,56 of insulin and should be featured in a reliable simulation. In comparison to AMBER, CHARMM produced more structural changes within our simulation time (Figure 2b). A break in the helix at the Gly 8 residue has been seen within the first 5 ns of the two independent calculations, followed by an attempt to retain the helical conformation between residues 4 and 7. These results are consistent with experimental studies showing that Gly residues are often associated with perturbations in structure and are known as “helix-breakers”.57 The break in the helix contributed to an increased flexibility of the Nterminus, which enabled the N- and C-termini to move close together, producing hydrogen-bonded bridges between residues 2 and 26. This feature can also be seen in the cluster structures of system b1 (Figure 3b), where the turn region between residues 19 and 22 transformed into a π-helix after 12 ns and maintained a stable π-helical conformation for the rest of the simulation. It has previously been reported that the formation of π-helix may be an artifact of a specific forcefield applied to a MD simulation with implicit solvent.58 While an implicit solvent simulation cannot be directly compared to simulations in explicit solvent such as those reported here, we must stress that π-helix has been suggested to be 10 times more prevalent than previously believed.59 Transitions between R- and π- helix have been previously reported in MD simulations and are believed to be genuine physical transitions and represent an important stage in the helix-to-coil transition.60 Interestingly, the formation of π-helices can also be observed in the GROMOSA6 forcefield simulations (Figure 2e). In contrast to system b1, there is very little π-helix formation seen in system b2, where most of the region between residues 19 and 22, retains its turn conformation. Noticeable structural changes were seen in both simulations using OPLS (Figure 2c), starting with a break at Ser 9, which continued into dissipation of the helix. Most helical content was lost by the end of our simulation time, with the majority of the peptide structure adopting a β-turn conformation. The loss in helical content can be observed in the cluster structures of both simulations using OPLS shown in Figure 3c, particularly in the second system c2, where very small amount of helical elements are seen. An interesting structural feature observed in the GROMOSA1 systems (Figure 2d) is the break in the helix at Gly 8, which leads to an attempt of helix reformation between residues 4 and 7. Evidence of helical conformation in residues 4 to 7 was observed previously in the NMR solution structure of an insulin chain B mutant22 and in our previous computational study using the CHARMM27 forcefield.46 In the present study, the helical region between residues 9 and 18 was well conserved throughout both simulations using GROMOSA1, with a slight fraying at residue 18, near the end of simulation d2 (Figure 2d). This was likely caused by the hydrogen bonding between the two termini, reflected by bridge in the secondary structure evolution plot of this system. This feature can be seen in the cluster structures of system d2 in Figure 3, presented as small parallel β-sheet with the hydrogen bonds colored in purple. The GROMOSA6 forcefield produced many changes in the secondary structure of the protein from its initial state (Figure 2e). A break in the helix at residue 10 occurred during the equilibration stage, which destabilized the protein and totally
J. Phys. Chem. B, Vol. 112, No. 35, 2008 11141 dissolved the helical conformation. A transition from R-helix to π-helix to turn conformation can also be seen at the beginning of the simulations. The formation of extended conformation between residues 9, 10 and 26, 27, represented in yellow in system e2 of Figure 3 is associated with the loss of the conserved helical region observed in most of the previous simulations, however no hydrogen bonds were detected. The total loss of helical content can be clearly seen in the most populated clusters of both simulations with GROMOSA6 shown in Figure 3e. NOE Energies and Violations. Simulated molecular dynamics trajectories of proteins and nucleic acids are often compared with nuclear magnetic resonance (NMR) data when assessing the quality of forcefield used, or in the interpretation of ambiguous experimental data.61 The atom-atom distances derived from the nuclear Overhauser enhancement (NOE) intensities are particularly useful, as they are integral in NMR-based protein structure determination61 and the finding of specific conformational features. The only experimental data currently available for the description of the 3D conformational features being investigated in this study are the NOE interproton distances. The data derived from the NMR structure of bovine insulin chain B in solution20 provide a useful benchmark for evaluating the conformations observed during the MD calculations. Hawkins et al.20 derived NOE constraints from a 250 ms NOESY spectrum of the oxidized chain B at 500 MHz, 300 K and pH of 2.2 to 2.5. As coordinates for the structure were not available, NOE distance restraints were used for comparison in this study. Using the square-well potential function62 implemented in the program XPLOR,63,64 the interproton distance restraints were used to calculate the NOE energy (ENOE) values for the MD conformations from the trajectories of each forcefield system. The square-well potential function is defined as:
{
(rij - rijupper)2,
ENOE ) 0, rlower ij
if rij > rijupper
if rijlower e rij e rijupper
(rij - rijlower)2,
if
(2)
rij < rijlower
rupper ij
where and are the lower and upper limit of the target distance restraint, respectively. The distance between selected set of atoms is averaged according to
(〈rij-6〉)-1/6
(3)
where rij runs through all possible combinations of distance restraints between atoms “i” and atoms “j” of the selected molecule. The NOE energy is widely used by experimentalists as a distance restraining function in addition to the standard empirical energy, eq 1, for simulating annealing, restrained minimization, and molecular dynamics simulations. This energy is a crucial inclusion in various methodologies used for NMR structure determination and refinement (see references within ref 65). In this study the NOE energy was used for evaluation of the conformity of the structures sampled by our simulations with the NMR distance restraints. To evaluate our results, we compared these pseudo energy values with the ENOE of the original PDB structure of each forcefield system, presented in Figure 4. The total numbers of interproton distance violations, along with the number of “severe” violations (1 Å above the upper bound value of the NOE restraint) are presented in Table 1. Additionally, several interesting structures were selected from the forcefield simulations and are shown in Figure 5, highlighting their satisfied and violated distances, along with the region of their occurrence.
11142 J. Phys. Chem. B, Vol. 112, No. 35, 2008
Todorova et al.
Figure 5. Important protein conformations and their NOE distance violations. The satisfied distance restraints are shown in yellow and the violated distances in black. The coloring of each protein is by its secondary structure as in Figure 2.
Figure 4. NOE energy derived from the distance restraints of the sampled structures. The region below the line at ∼41.5 × 103 au represents the MD structures with NOE energy less than that of the starting structure. The plot shows a moving average of 250 ps.
The NOE energies calculated from the simulations using AMBER were mostly higher than the energy of the starting structure (Figure 4a). A decrease in energy is seen after 10 ns for system a1 indicating favorable conformational changes. This region of the secondary structure plot corresponds to the turn conformation forming near the β-turn of the peptide. The conformations sampled at ∼12.8 ns and their NOE distance violations were closely analyzed for both simulations using AMBER (Figure 5 a1 and a2) because of their significant NOE energy variance. The structural differences in the turn region of both structures are evident, with the extension of turn conformation observed in simulation a1 and the fully retained helical region in simulation a2. These features produced a large difference in the NOE distances between residues 22 and 26 of both structures, in particular, that of structure a2 which was more than double in length, ∼12 Å, of the upper bound value of 6 Å of that restraint, compared to ∼8 Å of structure a1. The following increase in energy is evoked by the rebuilding of the helix. The turn conformation extension at the C-terminus produced a decrease in energy near the end of the simulation. Interestingly, very similar NOE violations were observed in a folding study of chain B of insulin using the AMBER03 forcefield (to be published elsewhere), where ∼1 ms of data was accumulated using the Bias-Exchange Metadynamics technique.66 This finding further highlights that adequate sampling has been acquired for this study. Overall, the protein
Comparison of Empirical Forcefields
J. Phys. Chem. B, Vol. 112, No. 35, 2008 11143
TABLE 1: Average NOE Interproton Distance Violations Calculated for Each Forcefield Simulationa forcefield AMBER03 CHARMM27 OPLS-AA GROMOS 43A1 GROMOS 53A6
system
total violations (309)
strong 2.7 Å (24)
medium 3.5 Å (141)
weak 5 Å (88)
very weak 6 Å (56)
a1 a2 b1 b2 c1 c2 d1 d2 e1 e2
13 9 11 11 17 16 7 8 14 16
2 2 1 0 4 0 1 0 1 1
1 1 1 1 1 1 1 1 1 1
3 2 3 3 3 3 3 3 8 8
1 1 1 1 4 3 1 1 2 2
a The total average number of violations, as well as the severe violations larger than 1 Å from the upper bound value of the classified cross-peak intensities by Hawkins et al.20 are indicated.
behavior produced by the AMBER forcefield, such as the restrictions of atomic movement and the tendency to keep helices well conserved is not in a good agreement with the NOE restraints of chain B of insulin. In the CHARMM simulations there is a variance in the NOE energy between the two separate systems (Figure 4b), reflected by the structural differences observed. Initially, the NOE energy is below that of the minimized original structure in both simulations, however an increase is seen in system b1 after ∼10 ns, where the formation of a π-helix and a bridge between residues 24 and 28 is observed. This indicates a less favorable secondary structure, which may be due to the bridge formation near the C-terminus. A decrease in energy is seen near the end of the simulation as the bridge disappears and the C-terminus extends and the π-helix reforms. Overall, system b2 produced the lowest NOE energies, indicating that the turn conformation of the chain at residues 17-19 can be more energetically favorable than a π-helix. This finding is in agreement with the distance violations experienced by the structures sampled at ∼30.2 ns of the two simulations using CHARMM in Figure 5. It is evident that the π-helix formation in structure b1 produced more severe interproton distance violations than the extended turn conformation of structure b2. Turn-like structures were also apparent over residues 17 to 21 of chain B at the carboxy terminal in the solution state study performed by Hawkins et al.20 The simulations using OPLS forcefield produced the least favorable results, with most of the energies above the NOE energy value of the initial structure (Figure 4c). These results correlate well with the calculated interproton violations (Table 1), which show that the reduction of helical content is the main contributor to the largest number of violations. The stretching of the protein structure and the violations experienced can be visualized in the conformation sampled at ∼20 ns of simulation c1 in Figure 5. The helix in chain B of insulin is a well conserved feature observed experimentally while in the OPLS calculations most helical content is lost. The simulations using the united-atom GROMOSA1 forcefield agreed with the experimental data better than all the other simulations (Figure 4d). In both simulations the NOE energies are below the energy of the initial structure, also having the lowest number of interproton distance violations (Table 1). The conformational changes such as the Cys 8 helix break are in agreement with the NOE restraints, showing a decrease in energy. This event was followed by the formation of a T-like state (associated with insulin’s activity37), represented by structure sampled at ∼10.24 ns of simulation d1 in Figure 5. The small magnitude of NOE distance violations indicates favorable structure. The increase in energy of system d2, near
the end of the simulation can be correlated with the close interaction of the termini, however this is probably not an unfavorable event because the NOE energies are still well below that of the starting structure. In contrast to GROMOSA1, GROMOSA6 produced less favorable results (Figure 4e). The increase in conformational changes produced large fluctuation in energy with most of the values being above that of the minimized structure. A high NOE energy conformation sampled at ∼ 40 ns of simulation e1 is presented in Figure 5, where the extent of the structural deformation and large number of violated distance restraints are clearly seen. As observed with OPLS, the dissipation of the helix is not a favored change in the chain B structure. The calculated NOE energies enabled us to compare the change in conformation of chain B relative to the starting structure with respect to the NOE intensities during each simulation. The widely used averaging method described by Zagrovic et al.61 was implemented to give the total violations within each forcefield system. The highest number of violations are found in the OPLS (17 and 16), and GROMOSA6 (14 and 16) simulations, and the least in the GROMOSA1 (7 and 8). Most violations are found for the weak intensities (Table 1) where the same NOEs found in the C-termini’s region of the protein are violated in each simulation. We have observed that all simulations produced a range of structures, however the number of severe violations compared to the total number of restraints, 309, is quite small. This seems to indicate that the NOE restraints allow some degree of flexibility in the molecule, which is in agreement with the range of experimental structures observed.20,37 These findings highlight the usefulness of performing multiple simulations in cases like this to explore the range of bioavailable conformations. Root Mean Square Deviation (Rmsd). The central helix region plays a key role in the insulin’s activity. To investigate the effect of each forcefield on the stability of the helix, rmsd calculations have been performed on this region. The central helix residues (Ser 9 and Cys 19) in each simulation were aligned and compared to the minimized PDB structure. The rmsd of the central helix for the ten systems is shown in Figure 6. The plots show a moving average of 250 ps to enable effective comparison of the overall trends between the systems. In the AMBER and GROMOSA1 systems, most of the helical region has retained its structure throughout the simulation, as illustrated by the most populated cluster structures in Figure 3. A small increase in rmsd is seen at the end of the simulation d2 of GROMOSA1, generated by the formation of turn conformation at the end of the helix. Transformation of part of the R-helical region into a π-helix (system b1) or turn conformation (system b2) contributed to the sharp increase in
11144 J. Phys. Chem. B, Vol. 112, No. 35, 2008
Todorova et al. TABLE 2: Average Radius of Gyration, Rg, with Standard Deviations Shown in Parentheses forcefield AMBER03 CHARMM27 OPLS-AA GROMOS 43A1 GROMOS 53A6
Figure 6. Rmsd of the helical region of chain B, from residues 9 to 19, compared to the starting structure as a function of time for all systems. The plot shows a moving average of 250 ps.
system
Rg (nm)
a1 a2 b1 b2 c1 c2 d1 d2 e1 e2
1.026 (0.027) 1.007 (0.042) 0.957 (0.059) 0.957 (0.049) 0.976 (0.047) 1.004 (0.045) 0.897 (0.049) 0.891 (0.051) 1.008 (0.055) 1.006 (0.061)
rmsd for the CHARMM simulations. However, this was stabilized for the rest of the simulations by the retained helix. The slow dissipation of the helix observed in the OPLS and GROMOSA6 simulations is obvious from the continuous increase in rmsd. This is almost double the value of rmsd obtained in our previous simulation of thermally stressed chain B at 400 K46 and points out the limitations of the forcefields in reproducing the experimentally observed behavior of the protein. Radius of Gyration. The radius of gyration (Rg) was calculated for each forcefield system to give an indication of the distribution of atoms relative to the proteins center of mass. It is used to measure the change in shape of the protein by comparing it to its starting structure. The average Rg for each forcefield system, together with its standard deviation, is presented in Table 2. The lowest average radius of gyration was produced by the GROMOSA1 forcefield, giving a value of 0.897 and 0.891 nm for the two trajectories, respectively. These results indicate that the structures formed are closely packed, as can be seen from the most populated cluster structures displayed in Figure 3. In addition, the secondary structure evolution plots show continuous interaction between the termini, producing a closed loop of the protein which results in a compact structure. Similar results were produced by the CHARMM forcefield, though with a slightly higher Rg average, likely caused by some loss of the helical content in the protein. The highest value of the radius of gyration was produced by the AMBER forcefield, closely followed by GROMOSA6 and OPLS. For AMBER this is not surprising because of the stability of the inherently elongated helical secondary structure observed in our simulations. For the other two forcefields the relatively high Rg indicates the disordering of the structure, including the loss of helical content and the lack of specific termini binding of the C-terminus with the helix, contributed to the higher average radius. Solvent Accessible Surface Area (SASA). The SASA was calculated for each forcefield system using a probe radius of 1.4 Å and Lennard-Jones hard-shell radii for each atom to define the surface. The average SASA of the hydrophobic and hydrophilic residues of chain B, along with the total area, including their standard deviations is presented in Table 3. The lowest SASA was produced by the GROMOSA1 forcefield, ∼4 nm2 smaller area than other systems. This result correlated well with the smallest radius of gyration also observed by this forcefield, confirming the compactness of the structures explored. The close interaction between the termini which reduced the accessible area to the solvent is the main reason for the low value. The other forcefield systems produced very similar values of the total SASA. The only significant difference was observed for the hydrophobic/hydrophilic SASA of the CHARMM simulations. Specifically, the SASA for the hydro-
Comparison of Empirical Forcefields
J. Phys. Chem. B, Vol. 112, No. 35, 2008 11145
TABLE 3: Average Solvent Accessible Surface Area, SASA, with Standard Deviations Shown in Parentheses forcefield AMBER03 CHARMM27 OPLS-AA GROMOS 43A1 GROMOS 53A6
system
hydrophobic (nm2)
hydrophilic (nm2)
total (nm2)
a1 a2 b1 b2 c1 c2 d1 d2 e1 e2
17.47 (0.75) 17.33 (0.87) 10.88 (0.77) 11.28 (0.76) 16.14 (0.70) 16.21 (0.87) 15.80 (1.01) 16.27 (0.98) 17.69 (0.80) 17.24 (0.84)
11.05 (0.54) 10.81 (0.66) 17.75 (0.86) 18.29 (0.88) 12.10 (0.59) 12.41 (0.65) 8.96 (0.56) 8.70 (0.66) 11.04 (0.63) 11.15 (0.72)
28.52 (1.04) 28.13 (1.35) 28.63 (1.38) 29.57 (1.37) 28.24 (0.99) 28.61 (1.30) 24.76 (1.26) 24.98 (1.28) 28.73 (1.06) 28.38 (1.07)
phobic residues is ∼ 6 nm2 lower compared to the corresponding structures produced by the other forcefields, while the exposed area of the hydrophilic residues is ∼6 nm2 higher. This result is produced mainly due to the compactness of the structure and the orientation of the side chains. This outcome is in agreement with accepted theory that proteins will fold to minimize the hydrophobic surface that is exposed to water.67 As can be seen from the cluster analysis of the CHARMM simulations (Figure 3b), the most populated structures are “closed” structures, where a strong interaction between the C- and N-terminii has led to their close proximity, thus reducing the amount of solvent accessing the hydrophobic residues. The outward orientation of some of the hydrophilic residues, such as Asn 3, His 5, Ser 9 and Lys 29, contribute to the slightly larger hydrophilic SASA. Cross-Forcefield Simulations. To further illustrate the possible systematic biases of the forcefields toward particular elements of secondary structure cross-forcefield simulations were performed on two predominant structures identified by different forcefields. First, an unfolded structure was taken from a GROMOSA6 simulation (Figure 3, e1) and was simulated using the AMBER03 forcefield for 15 ns. During this time, the secondary structure of the protein remained unfolded without significant structural change, and no helical formation was observed. This finding, together with the two previously described simulations performed on the native structure of chain B, suggests that the AMBER03 forcefield may have very steep and rugged free energy surface, where a protein can get easily trapped in a local minimum and even a long simulation may fail to sample the remaining conformational space. Second, one of the most frequently sampled structure (Figure 3, b2) by the CHARMM forcefield was taken and simulated using the OPLS for 15 ns. Significant change in the conformation of the protein was seen after 7.5 ns of MD, where reduction in helical content became most evident. This result confirms possible favoring of the OPLS forcefield of the turn conformation as it is in agreement with the behavior observed in the two simulations of natively folded chain B, described above. Conclusion Consistent systematic comparison of multiple simulations of insulin chain B using five different forcefields was performed to gain an improved understanding of the forcefield influences on the representation of the conformational behavior of the protein. The effect of these widely used forcefields on the secondary structure of insulin and its dynamics were investigated in detail by calculating the conformational evolution, solvent accessible surface area, radius of gyration and interproton distance violations for each forcefield simulation. Comparison of our results with X-ray crystallographic structures and NOE
distance restraints is, to the best of our knowledge, the only experimental data available for the description of conformational features we are interested to identify. We have observed that different forcefields favor different conformational trends, which is important to be aware of for the analysis of any classical simulation results. The all-atom forcefield CHARMM27 and the united-atom forcefield GROMOS 43A1 gave the best representation of the experimentally observed dynamic behavior of chain B of insulin. Similar structural trends were observed with both forcefields, including breaking of the helix at Gly 8, and the formation of some biologically important conformational states, such as the T and Rf states. The majority of the simulated structures of these two forcefields satisfied the NMR distance restraints, although GROMOS 43A1 had lower NOE energy and fewer violations. Alternatively, CHARMM27 was most effective in reducing the exposure of the hydrophobic surface of the peptide to the solvent. The AMBER03 forcefield produced well conserved helical regions, and did not lead to many interproton distance violations. However, the resultant loss of inherent flexibility of the chain limited the exploration of the conformational space, preventing the sampling of the biologically important structures. The OPLS-AA and GROMOS 53A6 forcefields also produced unfavorable results, where the simulated structures had very high NOE energies, indicating a large number of distance restraint violations, mostly contributed by the loss of the helical structure. It must be stressed, however, that even though some forcefields considered here may not have reproduced all features of the normal behavior of chain B within our simulated time, it does not mean that this behavior will not be observed if longer simulations are conducted. Overall, we can conclude that the CHARMM27 and GROMOS 43A1 forcefields appear to be the most suitable for our further studies of insulin behavior. Acknowledgment. The authors thank the Australian Research Council and Cytopia Pty. Ltd. for providing funding for this project, the Australian Partnership for Advanced Computing and Victorian Partnership for Advanced Computing for providing grants of computer time. We thank, Dr Stefano Piana, Dr Andrew Hung and Dr Akin Budi for the helpful discussions and Dr Mark Abraham from Australian National University for the help with the implementation of the CHARMM27 forcefield in the Gromacs simulation package. References and Notes (1) Ponder, J. W.; Case, D. A. AdV. Protein Chem. 2003, 66, 27. (2) MacKerell Jr, A. D. J. Comput. Chem. 2004, 25, 1584. (3) MacKerell Jr, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wio´rkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (4) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999. (5) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474. (6) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hunenberger, P. H.; Kruger, P.; Mark, A. E.; Scott., W. R. P.; Tironi, I. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; ETH Zurich and BIOMOS b.V; Zurich: Gronigen, 1996. (7) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2004, 25, 1656. (8) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (9) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049.
11146 J. Phys. Chem. B, Vol. 112, No. 35, 2008 (10) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins: Structure, Function, and Bioinformatics 2006, 65, 712. (11) Hess, B.; van der Vegt, N. F. A. J. Phys. Chem. B 2006, 110, 17616. (12) Price, D. J.; Brooks, C. L., III J. Comput. Chem. 2002, 23, 1045. (13) Villa, A.; Fan, H.; Wassenaar, T.; Mark, A. E. J. Phys. Chem. B 2007. (14) Mu, Y. G.; Kosov, D. S.; Stock, G. J. Phys. Chem. B 2003, 107, 5064. (15) Hu, H.; Elstner, M.; Hermans, J. Proteins: Struct., Funct. Genetics 2003, 50, 451. (16) Sorin, E. J.; Pande, V. S. J. Comput. Chem. 2005, 26, 682. (17) Zaman, M. H.; Shen, M.-Y.; Berry, R. S.; Freed, K. F.; Sosnick, T. R. J. Mol. Biol. 2003, 331, 693. (18) Yoda, T.; Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 2004, 386, 460. (19) Luo, R. Z.; Beniac, D. R.; Fernandes, A.; Yip, C. C.; Ottensmeyer, F. P. Science 1999, 285, 1077. (20) Hawkins, B.; Cross, K.; Craik, D. Int. J. Peptide Protein Res. 1995, 46, 424. (21) Qiao, Z. S.; Min, C. Y.; Hua, Q. X.; Weiss, M. A.; Feng, Y. M. J. Biol. Chem. 2003, 278, 17800. (22) Dupradeau, F. Y.; Richard, T.; Le Flem, G.; Oulyadi, H.; Prigent, Y.; Monti, J. P. J. Peptide Res. 2002, 60, 56. (23) Pullen, R. A.; Lindsay, D. G.; Wood, S. P.; Tickle, I. J.; Blundell, T. L.; Wollmer, A.; Krail, G.; Brandenburg, D.; Zahn, H.; Gliemann, J.; Gammeltoft, S. Nature 1976, 259, 369. (24) Mirmira, R. G.; Nakagawa, S. H.; Tager, H. S. J. Biol. Chem. 1991, 266, 1428. (25) Derewenda, U.; Derewenda, Z.; Dodson, E. J.; Dodson, G. G.; Bing, X.; Markussen, J. J. Mol. Biol. 1991, 220, 425. (26) Smith, G. D.; Dodson, G. G. Biopolymers 1992, 32, 441. (27) Ciszak, E.; Smith, G. D. Biochemistry 1994, 33, 1512. (28) Pittman, I. t.; Tager, H. S. Biochemistry 1995, 34, 10578. (29) Derewenda, U.; Derewenda, Z.; Dodson, G. G.; Hubbard, R. E.; Korber, F. Br. Med. Bull. 1989, 45, 4. (30) Ludvigsen, S.; Olsen, H. B.; Kaarsholm, N. C. J. Mol. Biol. 1998, 279, 1. (31) Zhang, Y.; Whittingham, J. L.; Turkenburg, J. P.; Dodson, E. J.; Brange, J.; Dodson, G. G. Acta Crystallogr. Sect. D: Biol. Crystallogr. 2002, 58, 186. (32) Derewenda, U.; Derewenda, Z.; Dodson, E. J.; Dodson, G. G.; Reynolds, C. D.; Smith, G. D.; Sparks, C.; Swenson, D. Nature 1989, 338, 594. (33) Baker, E. N.; Blundell, T. L.; Cutfield, J. F.; Cutfield, S. M.; Dodson, E. J.; Dodson, G. G.; Hodgkin, D. M. C.; Hubbard, R. E.; Isaacs, N. W.; Reynolds, C. D.; Sakabe, K.; Sakabe, N.; Vijayan, N. M. Philos. Trans. R. Soc. London, Ser. B: Biol. Sci. 1988, 319, 369. (34) Smith, G. D.; Swenson, D. C.; Dodson, E. J.; Dodson, G. G.; Reynolds, C. D. Proc. Natl. Acad. Sci. U.S.A.: Biol. Sci. 1984, 81, 7093. (35) Kaarsholm, N. C.; Ko, H. C.; Dunn, M. F. Biochemistry 1989, 28, 4427. (36) Yao, Z. P.; Zeng, Z. H.; Li, H. M.; Zhang, Y.; Feng, Y. M.; Wang, D. C. Acta Crystallogr., Sect. D: Biol. Crystallogr. 1999, 55, 1524.
Todorova et al. (37) Hua, Q. X.; Weiss, M. A. Biochemistry 1991, 30, 5505. (38) Adams, M. J.; Dodson, E.; Dodson, G. G.; Vijayan, M.; Baker, E. N. Nature 1969, 224, 491. (39) Peking, I. S. G. Peking ReV. 1971, 40, 11. (40) Bentley, G.; Dodson, E.; Dodson, G.; Hodgkin, D.; Mercola, D. Nature 1976, 261, 166. (41) Bentley, G. A.; Brange, J.; Derewenda, Z.; Dodson, E. J.; Dodson, G. G.; Markussen, J.; Wilkinson, A. J.; Wollmer, A.; Xiao, B. J. Mol. Biol. 1992, 228, 1163. (42) Falconi, M.; Cambria, M. T.; Cambria, A.; Desideri, A. J. Biomol. Struct. Dyn. 2001, 18, 761. (43) Zoete, V.; Meuwly, M.; Karplus, M. J. Mol. Biol. 2004, 342, 913. (44) Budi, A.; Legge, F. S.; Treutlein, H.; Yarovsky, I. Eur. Biophys. J. Biophys. Lett. 2004, 33, 121. (45) Budi, A.; Legge, F. S.; Treutlein, H.; Yarovsky, I. J. Phys. Chem. B 2005, 109, 22641. (46) Legge, F. S.; Budi, A.; Treutlein, H.; Yarovsky, I. Biophys. Chem. 2006, 119, 146. (47) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Prentice Hall: Harlow, England, 2001. (48) van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (49) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (50) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (51) Cheatham, T. E.; Miller, J. L.; Fox, T.; Darden, T. A.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 4193. (52) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (53) Freddolino, P. L.; Liu, F.; Gruebele, M.; Schulten, K. Biophys. J. 2008, 94, L75. (54) Frishman, D.; Argos, P. Proteins 1995, 23, 566. (55) Nakagawa, S. H.; Tager, H. S. J. Biol. Chem. 1986, 261, 7332. (56) Nakagawa, S. H.; Tager, H. S. J. Biol. Chem. 1987, 262, 12054. (57) Okamoto, Y.; Hansmann, U. H. E.; Nakazawa, T. Chem. Lett. 1995, 391. (58) Feig, M., Jr.; Brooks, C. L. J. Phys. Chem. B 2003, 107, 2831. (59) Fodje, M. N.; Al-Karadaghi, S. Protein Eng. 2002, 15, 353. (60) Armen, R.; Alonso, D. O.; Daggett, V. Protein Sci. 2003, 12, 1145. (61) Zagrovic, B.; van Gunsteren, W. F. Proteins: Struct., Funct., Bioinform. 2006, 63, 210. (62) Clore, G. M.; Nilges, M.; Sukumaran, D. K.; Brunger, A. T.; Karplus, M.; Gronenborn, A. M. Embo J. 1986, 5, 2729. (63) Brunger, Q. ReV, Biophys, 1993, 26, 49. (64) Brunger. X-PLOR, a system for crystallography and NMR.; Yale University: New Haven, CT, 1992. (65) Schwieters, C. D.; Kuszewski, J. J.; Clore, G. M. Prog. Nucl. Magn. Reson. Spectrosc. 2006, 48, 47. (66) Piana, S.; Laio, A. J, Phys, Chem, B 2007, 111, 4553. (67) Dyson, H. J.; Wright, P. E.; Scheraga, H. A. PNAS 2006, 103, 13057.
JP076825D
