Observations on AMBER Force Field Performance under the

Subscriber access provided by Gothenburg University Library

Article

Observations on AMBER Force Field Performance under the Conditions of Low pH and High Salt Concentrations Hanzhong Liu, Qingzhe Tan, Li Han, and Shuanghong Huo J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b07528 • Publication Date (Web): 29 Sep 2017 Downloaded from http://pubs.acs.org on October 1, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Observations on AMBER Force Field Performance under the Conditions of Low pH and High Salt Concentrations Hanzhong Liu,1 Qingzhe Tan,1 Li Han,2* Shuanghong Huo1* 1

Gustaf H. Carlson School of Chemistry and Biochemistry, Clark University, 950 Main

Street, Worcester, Massachusetts 01610, USA 2

Department of Mathematics and Computer Science, Clark University, 950 Main Street,

Worcester, Massachusetts 01610, USA *Corresponding authors: [email protected], [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Molecular dynamics simulations have become an important tool for the study of structures, dynamics, and functions of biomolecules. Timescales and force field accuracy are key factors for the reliability of these simulations. With the advancement of computational platforms and simulation technologies, all-atom simulations of proteins in explicitly represented aqueous solutions can reach as long as milliseconds. However, there are indications of minor force field flaws in literature. In this work we present our observations on force field accuracy under uncommon conditions. We performed molecular dynamics simulations of BBL refolding in aqueous solutions of acidic pH and high salt concentrations (up to 6 M) with the AMBER99SBILDN force field for a microsecond timescale. The reliability of the simulations relies on the accuracy of the physical models of protein, water, and ions. Our simulations show the same trend as experiments: higher salt concentration facilities refolding. However, we have observed the presence of β−sheet in the native helical region as well as α−helix and β−sheet in the native loop region. Some of the nonnative secondary structures are even more stable than native helices. Aside from the secondary structure issue under the uncommon conditions, we have found that the rigidity of glycine dihedral angles in the loop region leads to low root-mean-square deviations but large topological differences from the native structure. Whether this is due to a force field deficiency or not needs further investigations. Recently developed ion parameters exhibit evident liquid features even in the 6 M LiCl solution. However, cation-anion interactions in TIP3P water still seem too strong, leading to high fractions of contact ion pairs. We do not find any specific ion-binding motif, thus we conclude that the effect of salt is a nonspecific electrostatic screening in our simulations. Our observations on the AMBER force field performance under acidic conditions and high salt concentrations show that simulations under extreme conditions can provide informative tests of physical models.


Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Molecular dynamics (MD) simulation is a valuable tool in studying structures, thermodynamics, and kinetics of biomolecules. New computational platforms and algorithms make it feasible to perform long simulations or collect thousands of shorter trajectories. For example, with the development of a special-purpose supercomputer, the timescale of simulating a protein (< 100 amino acids) with an all-atom model and explicit solvent has reached milliseconds;1 cloud-based simulations of G-protein-coupled receptor activation pathways have also been in the millisecond timescale;2 with inexpensive hardware and implicit solvation, folding simulations of proteins with diverse sizes, secondary structures, and topologies can be finished in days.3 Folding@home distributed computing platform4 is a viable alternative, which enables one to collect thousands of shorter trajectories for a total duration of milliseconds.5 Aside from the simulation timescale, another important factor that affects the reliability of simulations and computed thermodynamic and kinetic properties for MD-based studies is the accuracy of force fields.6-7 Typical force field validations include validating secondary structural propensities,8-12 examining force fields’ capabilities to reproduce NMR measurements,8,

10, 13-14

and folding

simulations of various α and β proteins.3, 8, 15-16 In general, the structures and dynamics of native states as well as overall folding rates are better predicted than detailed kinetics, unfolded structures, and folding enthalpies.7,

17

Newly developed residue-specific force field based on

AMBER99SB showed more significant improvements in predicting folding enthalpies and entropies than in reproducing the 3JHNHα coupling constants.18,19 Folding simulations of helix bundle proteins suggest minor force field deficiencies,16 which may be more easily pinpointed in some systems than others. For example, minor differences in force fields may result in significant changes in the obtained thermodynamic and kinetic quantities as well as structural features of proteins with relatively low Tm (thermal denaturation midpoint). Thus, folding simulations of



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

such proteins could be an informative test of force field accuracy, which will in turn benefit simulations of all proteins. In this work we chose BBL as such a model system. The Tm of BBL ranges from 295 K to 335 K, depending on the protein constructs, experimental conditions, and biophysical techniques used,20-22 with a marginal folding barrier estimated as 1–2 kBT.23 Recent studies showed that salts (LiCl, NaCl, and CsCl) at concentrations above 1 M can induce complete refolding of acid-denatured BBL.24 Simulations of protein refolding under acidic conditions and high salt concentrations pose even greater challenges to the current force fields. And to our knowledge no tests have been done under such conditions. We simulated the refolding of BBL under acidic conditions in the 2 M, 4 M, and 6 M solutions of LiCl with explicit water in microsecond timescales. We chose AMBER99SB-ILDN25 force field because it is one of the recently improved force fields and its performance is comparable to CHARMM22*.8, 14 In the AMBER99SB-ILDN force field corrections have been made on side chain potentials of ILE, LEU, ASP, and ASN of AMBER ff99SB26 which is in turn built on AMBER ff99.27 The newly developed alkali and halide monovalent ion parameters were employed in this work.28 These ion parameters have been tuned to balance crystal and solution properties with different water models, including TIP3P29 and TIP4PEWD.30 The ion parameter set has given excellent agreements with reference values with both TIP3P and TIP4PEWD models.28 In this paper we present our observations on the force field deficiencies and analysis approaches that can be used in future force field assessments.


Page 4 of 29

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Methods System preparation To mimic the low pH environment, both GLU and ASP of BBL (PDB ID: 1W4H) (Fig. 1A) were protonated, and HIS was doubly protonated. The N-terminus was set to −NH3+ and the Cterminus was patched with the N-methyl group. The pH under this condition, therefore, is lower than the experimental value (pH = 3). The unfolded structure (Fig. 1B) was generated by an 1 nsunfolding simulation at 400 K starting from the protonated structure with the AMBER 11 package31 and AMBER99SB-ILDN force field. Because the purpose was to obtain an unfolded structure, the generalized Born (igb= 5) solvation model was used to facilitate the unfolding simulation.

Figure 1. Native structure (A) of BBL and an unfolded structure (B). The two helices in the native states are in the regions from residue 8 to 17 (helix 1) and from residue 35 to 44 (helix 2). A 310 helix is from residue 20-23. All the protein conformations in this work were depicted using UCSF Chimera.32 The color scheme is rainbow with the N-terminus in blue and the C-terminus in red.

The unfolded protein was then solvated in 2 M and 4 M solutions of LiCl using Maestro33 (version 9.3). The AMBER99SB-ILDN force field for proteins and the TIP3P explicit water model were employed. The recently developed ion (Li+ and Cl- ) parameters28 that were tuned to better balance crystal and solution properties in Ewald simulations with TIP3P and some other explicit water models were employed. The distance between the protein and the closest edge of the water box was set to 16 Å. Ten Cl¯ ions were added as counter ions to neutralize the system. 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Additionally, 431 and 922 LiCl were added to reach the desired salt concentration of 2 M and 4 M, respectively. The final systems contain 35,357 atoms/ions and 34,243 atoms/ions for the salt concentrations of 2 M and 4 M solutions, respectively. However, Maestro failed to solvate the protein with a 6 M solution of LiCl, probably due to the extremely high concentration. GROMACS34 (version 4.5.5), instead, was used to set up the system for a total of 27,059 atoms/ions including 1054 LiCl with the same force field. To minimize the proteins in the 2 M and 4 M solutions of LiCl, a hybrid of the steepest decent (10 steps) and the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) (2000 steps at most) implemented in Desmond (3.0.11)35 was used. There were two stages of minimization: the positions of protein atoms were restrained by harmonic forces with the force constant equal to 50 kcal/(mol Å ), followed by the minimization without any restraints. The protein in the 6 M solution of LiCl was minimized with the steepest descent algorithm until the maximum iteration step reached 50,000 or reached the convergence criterion (1000 kJ/mol/nm), whichever came first by GROMACS. To equilibrate protein in the 2 M and 4 M solutions of LiCl, Desmond standard NPT relaxation protocol was used. The system was first equilibrated for 1.2 ns in the NVT ensemble at 10 K using the Berendsen thermostat36 with harmonic restraints (force constant = 50.0 kcal/(mol Å2)) on the protein heavy atoms, followed by another round equilibration for 1.2 ns in the NPT ensemble at 10 K using the Berendsen thermostat and barostat with the same restraints. After that, the temperature of system was increased to 300 K linearly over the time period of 1.2 ns with restraints. The protein in the 6 M solution of LiCl was equilibrated for 1 ns using GROMACS with the same restraints using Berendsen thermostat in the NVT ensemble at 300 K followed by another equilibration in the NPT ensemble for 1 ns. Then the system was equilibrated without any restraints in the NPT ensemble for 10 ns at 300 K. Finally, all the systems were equilibrated for


Page 6 of 29

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50 ns to 150 ns using Desmond to ensure that the potential energy fluctuation is less than 0.5 % within the last 50 ns. The production runs were all carried out on Anton with the same force field, water model, and ion parameters. In total, four trajectories of 28.8 µs were obtained: one 7.2 µs−trajectory was collected for the protein in the 2 M and 4 M solutions of LiCl, respectively, and two 7.2 µs−trajectories were generated for the protein in the 6 M solution of LiCl. The Berendsen thermostat and barostat were employed to regulate temperature (= 300 K) and pressure (= 1 bar) in the NPT ensemble. A multiple time step approach37 (bonded, near, and far time steps = 2 fs, 2 fs, and 6 fs) was used for the integration of the equations of motion. The long-range electrostatic interactions were computed using the k-space Gaussian split Ewald method38 with the grid size of 32 × 32 × 32. The cutoff of van der Waals and short-range electrostatic interactions was optimized by the MD package for each system. Depending on the box size and grid size, the corresponding cutoffs of BBL in the 2 M, 4 M, and 6 M solutions of LiCl were 15.7 Å, 15.9 Å and 14.9 Å, respectively. A snapshot was saved every 240 ps, resulted in 29743, 29867, 30019, and 30022 proteins conformations in the 2 M, 4 M, and 6 M (two trajectories) solutions of LiCl, respectively. Analysis To characterize the structural features of the collected conformations, we calculated the secondary structures, RMSD (root-mean-square deviation), total surface area, and surface area of hydrophobic residues using Visual Molecular Dynamics (VMD).39 For the RMSD analysis the structure of BBL resolved at neutral pH was used as the native structure (PDB ID: 1W4H). Because of the flexibility of the N-terminus, we selected Cα atoms from residue 8 to 45 to calculate Cα−RMSD with respect to the native structure. The secondary structure assignment of a given snapshot was carried out with STRIDE40 implemented in VMD. The helicity of a particular residue is defined as the percentage of the number of snapshots of which the residue adopts a 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

helical conformation out of the total number of conformations. The radius of probe used to calculate the solvent accessible surface area was set to 1.4 Å. When each of the structural descriptors was plotted against time, the average values in every 24 ns were used to obtain a smooth curve. SHIFTX241 was used to back calculate the native BBL 1H and 15N chemical shifts as well as to predict the ensemble average chemical shifts of the low-RMSD (< 5 Å to the native) conformations in the 4 M and 6 M LiCl solutions. The distances between side-chain protons of the hydrophobic core (formed by Leu14, His17, Leu19, Ala21, Ile24, Leu33, Val38, His41, and Leu42) were calculated to assess the formation of tertiary structures. The Cl--protein side chain radial distribution function is defined as −

Cl 1 δN prot ( r ) Cl g prot ( r ) = Cl − , δV ρ −

where ρ

Cl−

(1)

is the density of chloride ion in the bulk solution. When the radial distribution

functions are integrated to their corresponding minima following the respective first peaks, the average number of ions accumulated around protein side chains in the first shell is obtained −

( N Cl prot ). The difference between the integrations of the second peak and the first peak gives the 42 average number of ions assembled in the second shell around the protein ( ∆N Cl prot ) side chains. −

We also used Scaled Gauss Metric (SGM)43-44 to measure the similarity between protein topologies. This metric has been used for automatic classification of protein structures44 and description of formation of native topology during protein folding.16 In mathematics, topology concerns with the properties that are preserved under continuous deformations through twisting, bending, and stretching in contrast to discontinuous changes such as tearing or merging. One key property used in the study of topological properties of curves is crossing, which can be used to characterize protein conformations. We briefly describe the concept here and refer readers to the Ref.43-45 for more details. Each protein conformation is considered as a polygonal space curve µ 8 ACS Paragon Plus Environment

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that connects all the Cα atoms, with one segment for every two adjacent Cα atoms. If two line segments cross each other when they are viewed from a certain direction, a sign is assigned to this crossing according to the right-hand rule (Fig. 2). Given two segments and of a polygonal protein curve, a topological characterization of these two segments, denoted by ( , ), is the signed probability of crossing averaged over all directions in space. For a protein with N Cα atoms, the sum of all ( , )s along the protein chain is defined as the first-order Gauss integral writhe of the curve (µ),

(,) (µ) =

( , ) (2),

which is interpreted as the signed average number of crossings when averaged over all directions in space with lying downstream of . The unsigned number of crossings is called average crossing number and defined as

|(,)| () = ∑ |( , )| (3). Both writhe and the average crossing number were used to distinguish protein conformations previously.46-47 In addition to the first-order measures, there are also higher-order measures43 and the following equation defines one second-order measure

(,)(,) =

( , )( , ) (4),

where the condition “0 < < < < < ” essentially enumerates all four tuples of Cα segments with ( , ) as one pair and ( , ) as the other pair. The order of these segments is from the N-terminal to the C-terminal.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Illustrations of line segment crossings of the Cα chain of BBL. The arrows point from the upstream residue to the downstream residue of the corresponding segments. The sign of a crossing is determined by the right-hand rule, screwing from the top segment to the one underneath. For the shown crossing, the view is in front of the paper and a negative crossing is seen. “N” and “C” denote the Nterminus and C-terminus, respectively.

Røgen defined a thirty-dimensional vector descriptor of topological properties: the first element of the vector is the number of Cα atoms and the rest are the crossing-based measures, including the first-order and higher order Gauss integral writhe and the unsigned number of crossings.43 The SGM between two conformations of a protein was computed as the Euclidean distances between their corresponding vectors.44 In contrast to RMSD, this process does not need pair-wise structure alignment, which leads to substantially reduced computational time. Topological properties can provide additional information about structure similarity that might not be captured by geometric measures, such as RMSD. Thus, we used both SGM and RMSD with respect to the native structure in our analysis of MD trajectories.


Page 10 of 29

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results and Discussion Ion parameters In the development of ion parameters, Joung et al. focused on fitting thermodynamic properties.28 In this work, we have investigated ion local structures that provide important information about ion interactions in an aqueous solution but not explicitly considered in the parameterization. We calculated the Li+ ̶ Cl- radial distribution functions for the aqueous solutions containing BBL and the 2 M, 4 M, and 6 M LiCl. As seen in Fig. 3A, the radial distribution functions show the characteristics of a liquid state: the distinct first peak (at 2.35 Å) and the second peak (at 4.65 Å) at short distance and a constant value at long distance. The peak positions are well in line with other simulation results using different ion parameters and water model.48-49 Recent success in the application of electron-transfer-mediated decay has provided important information on local chemical environment of ions in aqueous solutions.49 Features of radial distribution function can be attributed to different ion pair structures based on the electrontransfer-mediated decay spectra combined with ab initio calculation and MD simulations of LiCl aqueous solutions.49 The first peak of the Li+ ̶ Cl- radial distribution function is attributed to a contact ion pair arrangement that is formed when one of the water molecules nearest to Li+ is replaced by Cl−. A solvent-shared ion pair structure, corresponding to the second peak, contains as few as only two bridging water molecules separating the two counter ions. In a solventseparated ion pair configuration, Cl− is in the presence of the third hydration shell of Li+, much further away than in a solvent-shared ion pair. Fig. 3B shows a histogram of the distance between any given Li+ and the closest Cl−. The integration of this histogram up to a certain distance gives the percentage of the ion pairs within the range: contact ion pairs are 42.0 %, 62.0 %, and 73.1 % in the 2 M, 4 M, and 6 M LiCl solutions, respectively (Fig. 3C). This result is different from that of a previous simulation that showed the fraction of contact ion pairs was 40 % even in 6 M LiCl solution using SPC/E (extended simple point charge model) water50 and different ion



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

parameters.49 Practically, Pluhařová et al. rescaled the ionic charges to obtain a better fitting of neutron scattering data.51

Figure 3. Radial distribution function and ion pair percentage. A) Li+−Cl− radial distribution functions for the aqueous solutions containing BBL and the 2 M, 4 M, and 6 M LiCl; B) Histogram of the distance between any given Li+ and the closest Cl−; C) Integration of this histogram in B) up to 8 Å gives the percentage of ion pairs within the range. Color scheme: BBL in 2 M solution of LiCl (black curve); in 4 M solution of LiCl (red curve); in 6 M solution of LiCl (blue: traj. 1 and green: traj. 2).

Effect of salt concentrations on BBL refolding To study the salt effect on the refolding of BBL, we analyzed the change of structural properties along the trajectories. The plot of the Cα RMSD with respect to the native structure against simulation time (Fig. 4A) shows that the protein refolds to ca. 5 Å RMSD away from the native state in the 4 M and 6 M solutions of LiCl whereas it stays at least 8 Å RMSD to the native in the 2 M solution. This trend is consistent with that observed in the experiments: higher LiCl concentration facilitates refolding. Note that the actual pH of our simulation is lower than 3 due to the protonation of all Asp, Glu, and His as well as the N-terminus. Therefore, the fact that BBL does not refold to low RMSD value in the 2 M LiCl solution in our simulation does not contradict the experimental observation that in the presence of 2 M LiCl, BBL refolds at pH 3 to the native structure at neutral pH.24 It is also not surprising that the refolding does not reach a substantially low RMSD to the native state even in the 6 M LiCl solution. As shown in the previous folding simulation of BBL at neutral pH for hundreds of microseconds, the most populated state is 4.8 Å away from its native structure with partially formed helix 2, suggesting that there are minor flaws in the force field.16


Page 12 of 29

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It was found experimentally that higher salt concentration results in higher helicity.24 Our results show that the helicity in the helix 1 region increases when the LiCl concentration is raised (Fig. 4B). The helicity grows from less than 20 % in the helix 1 region in the 2 M solution to an average of 75 % in the 4 M solution and reaches the full length in the 6 M solution. Both experimental and simulation data show that helix 2 is less stable than helix 1 and forms in the late stage of folding.16, 21 The calculation of solvent accessible surface area shows that most of the conformations are as compact as the native state conformation or even exceeds the native compactness (Fig. 4C). Non-polar solvation free energy is proportional to the solvent accessible surface area with the proportional constant described as the microscopic surface tension that increases with ionic strength.52-54 To minimize the non-polar solvation free energy, the protein surface area will decrease in solutions of high salt concentrations to compensate the increase in the microscopic surface tension. Therefore, we have observed that the surface areas of most conformations in the 2 M-6 M solutions of LiCl are less than that of the native state. However, the smaller-than-native surface areas do not necessarily lead to smaller-than-native hydrophobic surface areas. The conformations collected under the conditions of the 2 M and 4 M solutions and the conformations near the end of the trajectory 2 in the 6 M solution exhibit larger hydrophobic surface areas than the native state (Fig. 4D).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Structural properties as a function of simulation time. (A) Cα RMSD of residues 8 to 45 with respect to the native structure of BBL (PDB ID: 1W4H) at neutral pH; (B) Helicity; (C) Surface area; (D) Surface area of hydrophobic residues. In (B) The native helical regions and the β-sheet regions predicted by PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred/) are shown on the top of the figure. In (C) and (D) the dashed lines denote the corresponding values of the native state. Color scheme: black and red are for the simulations in the 2 M and 4 M solutions of LiCl, respectively; blue and green are for trajectory 1 and trajectory 2 of BBL in the 6 M solution, respectively.

The hydrophobic core of BBL is formed by Leu14, Leu19, Ala21, Ile24, Leu33, Val38, and Leu42 as well as the two histidines (His17 and His41) in the native structure. The NOESY spectrum showed the same long range NOEs at pH 3 in the LiCl solution as expected from the native structure at neutral pH,24 indicating that salt has induced native tertiary contacts. For the conformations with low RMSDs (< 5 Å to the native) in the 4 M and 6 M LiCl solutions, we calculated the average proton-proton distances for all the residue pairs that exhibit the long range NOEs in Ref.24 Like the trend seen in the secondary structure formation in different solutions, higher salt concentration promotes native-like tertiary contacts. Particularly, among the helix 1, 310 helix, and loop regions the protein conformations in the 6 M solution exhibit more native-like


Page 14 of 29

Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tertiary contacts than those in the 4 M solution (Fig. 5). Nevertheless, the native contacts between helix 1 and helix 2 are not formed in either solution.

Figure 5. Average proton-proton distances in the hydrophobic core. The side chain protons used are Hδ of Leu, Hδ of His, Hβ of Ala, Hδ of ILE, and Hγ of Val. The residue names and numbers are color-coded for helix 1 (cyan) and helix 2 (orange), and those in the 310 helix (residues 20-23) and its downstream loop are in black. The black curve denotes the distances in the native structure. The red and blue curves represent the distances in the conformations with low RMSDs with respect to the native structure (< 5 Å ) in the 4 M and 6 M LiCl solutions, respectively. Some of the error bars are too small to show.

We analyzed the secondary structure of each snapshot along the trajectories (Fig. 6). As shown in Fig. 6A, no helical segment is maintained for more than 1 µs in the native helix 1 region in the simulation of BBL in the 2 M solution of LiCl, whereas two stable β strands have been present for five µs in the native helix 2 region. Short β strands are also seen in the native 310 helical region and its upstream and downstream. As seen in Fig. 7, the content of a two-strand βsheet in the native helix 2 region remains for more than 60% of the trajectory simulated in the 2 M solution of LiCl. For the simulation under the 4 M salt concentration, β-strand content is at most 20% occupancy in two regions (Fig. 6B and Fig. 7): from residue 22 to 25 and from residue 32 to 36. These β-strand regions overlap with the β-sheet found in a prior simulation work55 employing AMBER parm9656 force field. And PSIPRED57 predicts that this region has certain βsheet propensity. A β-sheet formed by residues 22-26 and residues 29-34 was found to be on the pathways of folding of BBL,58 but the β conformations are not verified by experiments. We



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

collected two independent trajectories of BBL in the 6 M solution of LiCl. While in both trajectories helix 1 is well preserved (Fig. 6C and D), differences are seen in its downstream regions. In trajectory 1, two segments of helices are present in the native loop region which includes the PSIPRED predicted β-sheet region, whereas in trajectory 2 a two-strand β-sheet with the first strand in the predicted β-sheet region and the second strand in the native helix 2 region has stabilized for 5 µs (Fig. 6D).

Figure 6. Secondary structures as a function of time in 2M (A), 4M (B), and 6M solutions of LiCl (trajectory 1 (C) and trajectory 2(D)). The pink color denotes helix and blue is for ! sheet. The native helix 1, helix 2, and 310 helix regions are denoted as H1, H2, and 310, respectively.

Figure 7. β-sheets content along different trajectories. Color Scheme: black and red are for the simulations in the 2 M and 4 M solutions of LiCl, respectively; blue and green are for trajectory 1 and trajectory 2 of BBL in the 6 M solution of LiCl, respectively.

The conformations that have been randomly selected from the time interval when RMSD, secondary structure, and surface area are relatively stable are depicted in Fig. 8. In addition to the 16 ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


low helical content and long lasting β-sheet, the overall shape of the protein in the 2 M solution of LiCl is totally different from that of the native state, leading to roughly 8 Å of RMSD away from the native structure (Fig. 8A). The illustrated structure of BBL in the 4 M solution of LiCl (Fig. 8B) resembles the overall shape of the native state even though its helical content is low. In trajectory 1 of BBL in the 6 M solution of LiCl, the protein folds to the same RMSD range as that in the 4 M solution (Fig. 4A) and the helix 1 is fully formed. However, the loop winds in the opposite direction to that of the native structure. In trajectory 2, the protein is trapped in an α/β conformation.

Figure 8. Selected conformations during the simulations in the 2M (A), 4M (B), and 6M solutions of LiCl (trajectory 1 (C), trajectory 2 (D)). The snapshots were selected from the time intervals of which RMSD, secondary structure, and surface area are relatively stable (in Fig. 4 and Fig. 6). Snapshot (A) is at 7.14 µs, RMSD = 8.1 Å; Snapshot (B) is at 3.35 µs, RMSD = 4.6 Å; Snapshot (C) is at 6.26 µs, RMSD = 4.9 Å; Snapshot (D) is at 5.42 µs, RMSD = 5.8 Å. All the RMSDs refer to the Cα RMSD of residues 8-45 with respective to the native structure (PDB ID: 1W4H, model 1). The color scheme is rainbow with the Nterminus in blue and the C-terminus in red.

Topological descriptor as a similarity measure for protein conformations Topological and geometrical measures can provide different and complementary information. We compare SGM with RMSD along each trajectory in Fig. 9. Overall, both measures give a similar trend, but differences are also significant in some time intervals. On one hand, if both values are small (or large) for two conformations, it indicates that these two



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

conformations are similar (or different). For example, in the simulation in the 4 M solution of LiCl the RMSD values are about 5 Å away from the native structure from 2.7 µs to 4.4 µs and the corresponding SGM values drop to their low range (Fig. 9B). Thus, the topological measure and geometric measure give a congruous description of the structure similarity in this time interval. On the other hand, if one value is small and the other is large, it indicates that neither measure is sufficient on its own and further investigation is needed to identify the cause of the difference and then draw conclusions on structure similarity accordingly.

Figure 9. Comparison between scaled Gauss metric (SGM) (black curve) and RMSD (red curve) with respect to the native structure along simulations. A-D: simulations in 2 M, 4 M, and 6 M (trajectory 1 and trajectory 2) solutions of LiCl, respectively.

In this work, we are more interested in the conformations that are near native in terms of RMSD (< 5 Å with respect to the native state in our case) while topologically they are not because these conformations and the cause of the difference in the similarity measures may provide useful information about the partially folded/misfolded conformations. As shown in Fig. 9C, when the RMSDs remain less than 5 # from 4.8 µs to 7.2 µs during the simulation in the 6 M 18 ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


solution of LiCl, the corresponding SGM values are around 110, exceeding that of the initial unfolded structure shown in Fig. 1B (SGM = 102). This suggests that the topological property of the structures are substantially different from that of the native.

Figure 10. Ramachandran plot of G28 of the low-RMSD conformations. (A) For the snapshots collected from 2.7 to 4.4 μs of the simulation in the 4 M solution of LiCl, for a total of 7200 conformations; (B) For the snapshots collected from 4.8 to 7.2 μs of traj. 1 of the simulation in the 6 M solution of LiCl, for a total of 10019 conformations. The snapshots in these two time periods correspond to the structures in the low RMSD range in Fig. 9B and Fig. 9C, respectively. (C) Supposition of the native structure and the structures with the smallest SGM in the low RMSD range of the simulations in the 4M and 6M solutions of LiCl, respectively. Color scheme: Native (brown); snapshot at 3.35 µs in the 4 M solution of LiCl (cyan); snapshot at 6.26 µs in the 6 M solution of LiCl (traj. 1, lavender). G28 of these structures are highlighted in red. “N” and “C” denote the N and C terminus, respectively.

To find out the origin of the topological differences between the low-RMSD ensembles (

Observations on AMBER Force Field Performance under the

Recommend Documents