Philippe Ferrara

Evaluation of a fast implicit solvent model for molecular dynamics simulations

A solvation term based on the solvent accessible surface area (SASA) is combined with the CHARMM polar hydrogen force field for the efficient simulation of peptides and small proteins in aqueous solution. Only two atomic solvation parameters are used: one is negative for favoring the direct solvation of polar groups and the other positive for taking into account the hydrophobic effect on apolar groups. To approximate the water screening effects on the intrasolute electrostatic interactions, a distance‐dependent dielectric function is used and ionic side chains are neutralized. The use of an analytical approximation of the SASA renders the model extremely efficient (i.e., only about 50% slower than in vacuo simulations). The limitations and range of applicability of the SASA model are assessed by simulations of proteins and structured peptides. For the latter, the present study and results reported elsewhere show that with the SASA model it is possible to sample a significant amount of folding/unfolding transitions, which permit the study of the thermodynamics and kinetics of folding at an atomic level of detail. Proteins 2002;46:24–33.

Calculation of conformational transitions and barriers in solvated systems: Application to the alanine dipeptide in water

Optimal free energy paths (OFEPs) for conformational transitions are parallel to the mean force at every nonstationary point of the free energy landscape. In contrast to adiabatic paths, which are parallel to the force, OFEPs include the effect of entropy and are relevant even for systems with diffusive degrees of freedom. In this study the OFEPs are computed for the alanine dipeptide in solution. The potential of mean force is calculated and an effective potential is derived that is used to obtain the paths with a minimization based algorithm. The comparison of the calculated paths with the adiabatic paths in vacuo shows the influence of the environment on conformational transitions. The dynamics of the alanine dipeptide in water are more complex, since there are more minima and the barriers are lower. Two simpler methods for the calculation of reaction pathways in solution are evaluated by comparing their results with the OFEPs. In the first method the mean electrostatic field of the water is approximated by an analytical continuum model. The obtained paths show qualitative agreement with the OFEPs and the height of the barriers are similar. Targeted molecular dynamics (TMD), the second approach, constrains the distance to a target conformation to accelerate the transition. In the general case, however, it is difficult to assess the physical significance of the obtained paths. Changing the initial conditions by assigning different velocities leads to different solutions for the conformational transition. Furthermore, it is shown that by performing the simulations with different reaction coordinates or in opposite directions different pathways are preferred. This result can be explained by the structure of the free energy landscape around the initial conformations. In a first approximation the physical significance of different pathways is assumed to depend mainly on the free energy at the highest saddle point. In the literature the total energy of the system has often been used to estimate the position and the height of the energy barriers in the path. By comparing the total energy with the calculated free energy it is shown that the former largely overestimates the height of the barriers. Furthermore, the positions of the maxima of the total energy do not coincide with the free energy barriers. Simple approximations to the free energy lead to good quantitative agreement.Optimal free energy paths (OFEPs) for conformational transitions are parallel to the mean force at every nonstationary point of the free energy landscape. In contrast to adiabatic paths, which are parallel to the force, OFEPs include the effect of entropy and are relevant even for systems with diffusive degrees of freedom. In this study the OFEPs are computed for the alanine dipeptide in solution. The potential of mean force is calculated and an effective potential is derived that is used to obtain the paths with a minimization based algorithm. The comparison of the calculated paths with the adiabatic paths in vacuo shows the influence of the environment on conformational transitions. The dynamics of the alanine dipeptide in water are more complex, since there are more minima and the barriers are lower. Two simpler methods for the calculation of reaction pathways in solution are evaluated by comparing their results with the OFEPs. In the first method the mean electrostatic field of the water is approximat...

Computer simulations of protein folding by targeted molecular dynamics

We have performed 128 folding and 45 unfolding molecular dynamics runs of chymotrypsin inhibitor 2 (CI2) with an implicit solvation model for a total simulation time of 0.4 microseconds. Folding requires that the three‐dimensional structure of the native state is known. It was simulated at 300 K by supplementing the force field with a harmonic restraint which acts on the root‐mean‐square deviation and allows to decrease the distance to the target conformation. High temperature and/or the harmonic restraint were used to induce unfolding. Of the 62 folding simulations started from random conformations, 31 reached the native structure, while the success rate was 83% for the 66 trajectories which began from conformations unfolded by high‐temperature dynamics. A funnel‐like energy landscape is observed for unfolding at 475 K, while the unfolding runs at 300 K and 375 K as well as most of the folding trajectories have an almost flat energy landscape for conformations with less than about 50% of native contacts formed. The sequence of events, i.e., secondary and tertiary structure formation, is similar in all folding and unfolding simulations, despite the diversity of the pathways. Previous unfolding simulations of CI2 performed with different force fields showed a similar sequence of events. These results suggest that the topology of the native state plays an important role in the folding process. Proteins 2000;39:252–260.

Weak temperature dependence of the free energy surface and folding pathways of structured peptides

The thermodynamics and energetics of a 20‐residue synthetic peptide with a stable three‐stranded antiparallel β‐sheet fold are investigated by implicit solvent molecular dynamics (MD) at 330 K (slightly above the melting temperature in the model) and compared with previous simulation results at 360 K. At both temperature values, the peptide folds reversibly to the NMR solution conformation, irrespective of the starting conformation. The sampling of the conformational space (2.3 μs and 25 folding events at 330 K, and 3 μs and 50 folding events at 360 K) is sufficient to obtain a thermodynamic description of minima and transition states on the free energy surface, which is determined near equilibrium by counting populations. The free energy surface, plotted as a function of two‐order parameters that monitor formation of either of the β‐hairpins, is similar at both temperature values. The statistically predominant folding pathway and its frequency (about two‐thirds of the folding events) are the same at 330 K and 360 K. Furthermore, the main unfolding route is the reverse of the predominant folding pathway. The effective energy and its electrostatic and van der Waals contributions show a downhill profile at both temperatures, implying that the free energy barrier is of entropic origin and corresponds to the freezing of about two‐thirds of the chain into a β‐hairpin conformation. The average folding rate is nearly the same at 330 K and 360 K, while the unfolding rate is about four times slower at 330 K than at 360 K. Taken together with previous MD analysis of α‐helices and β‐hairpins, the present simulation results indicate that the free energy surface and folding mechanism of structured peptides have a weak temperature dependence. Proteins 2002;47:305–314.

Flexibility of the murine prion protein and its Asp178Asn mutant investigated by molecular dynamics simulations

Inherited forms of transmissible spongiform encephalopathy, e.g. familial Creutzfeldt-Jakob disease, Gerstmann-Sträussler-Scheinker syndrome and fatal familial insomnia, segregate with specific point mutations of the prion protein. It has been proposed that the pathologically relevant Asp178Asn (D178N) mutation might destabilize the structure of the prion protein because of the loss of the Arg164-Asp178 salt bridge. Molecular dynamics simulations of the structured C-terminal domain of the murine prion protein and the D178N mutant were performed to investigate this hypothesis. The D178N mutant did not deviate from the NMR conformation more than the wild type on the nanosecond time scale of the simulations. In agreement with CD spectroscopy experiments, no major structural rearrangement could be observed for the D178N mutant, apart from the N-terminal elongation of helix 2. The region of structure around the disulfide bridge deviated the least from the NMR conformation and showed the smallest fluctuations in all simulations in agreement with hydrogen exchange data of the wild type prion protein. Large deviations and flexibility were observed in the segments which are ill-defined in the NMR conformation. Moreover, helix 1 showed an increased degree of mobility, especially at its N-terminal region. The dynamic behavior of the D178N mutant and its minor deviation from the folded conformation suggest that the salt bridge between Arg164 and Asp178 might not be crucial for the stability of the prion protein.

Calculation of protein ionization equilibria with conformational sampling: pKa of a model leucine zipper, GCN4 and barnase †

The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pKa calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pKa estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve‐averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pKa calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pKa predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pKa predictions for the model leucine zipper was reduced to as low as approximately 0.25 pKa units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pKa calculations, are discussed. Proteins 2002;46:41–60.