Conservation of coevolving protein interfaces bridges prokaryote-eukaryote homologies in the twilight zone.

Protein-protein interactions are fundamental for the proper functioning of the cell. As a result, protein interaction surfaces are subject to strong evolutionary constraints. Recent developments have shown that residue coevolution provides accurate predictions of heterodimeric protein interfaces from sequence information. So far these approaches have been limited to the analysis of families of prokaryotic complexes for which large multiple sequence alignments of homologous sequences can be compiled. We explore the hypothesis that coevolution points to structurally conserved contacts at protein-protein interfaces, which can be reliably projected to homologous complexes with distantly related sequences. We introduce a domain-centered protocol to study the interplay between residue coevolution and structural conservation of protein-protein interfaces. We show that sequence-based coevolutionary analysis systematically identifies residue contacts at prokaryotic interfaces that are structurally conserved at the interface of their eukaryotic counterparts. In turn, this allows the prediction of conserved contacts at eukaryotic protein-protein interfaces with high confidence using solely mutational patterns extracted from prokaryotic genomes. Even in the context of high divergence in sequence (the twilight zone), where standard homology modeling of protein complexes is unreliable, our approach provides sequence-based accurate information about specific details of protein interactions at the residue level. Selected examples of the application of prokaryotic coevolutionary analysis to the prediction of eukaryotic interfaces further illustrate the potential of this approach.

From residue coevolution to protein conformational ensembles and functional dynamics.

The analysis of evolutionary amino acid correlations has recently attracted a surge of renewed interest, also due to their successful use in de novo protein native structure prediction. However, many aspects of protein function, such as substrate binding and product release in enzymatic activity, can be fully understood only in terms of an equilibrium ensemble of alternative structures, rather than a single static structure. In this paper we combine coevolutionary data and molecular dynamics simulations to study protein conformational heterogeneity. To that end, we adapt the Boltzmann-learning algorithm to the analysis of homologous protein sequences and develop a coarse-grained protein model specifically tailored to convert the resulting contact predictions to a protein structural ensemble. By means of exhaustive sampling simulations, we analyze the set of conformations that are consistent with the observed residue correlations for a set of representative protein domains, showing that (i) the most representative structure is consistent with the experimental fold and (ii) the various regions of the sequence display different stability, related to multiple biologically relevant conformations and to the cooperativity of the coevolving pairs. Moreover, we show that the proposed protocol is able to reproduce the essential features of a protein folding mechanism as well as to account for regions involved in conformational transitions through the correct sampling of the involved conformers.

Large-Scale Conformational Transitions and Dimerization Are Encoded in the Amino-Acid Sequences of Hsp70 Chaperones.

Hsp70s are a class of ubiquitous and highly conserved molecular chaperones playing a central role in the regulation of proteostasis in the cell. Hsp70s assist a myriad of cellular processes by binding unfolded or misfolded substrates during a complex biochemical cycle involving large-scale structural rearrangements. Here we show that an analysis of coevolution at the residue level fully captures the characteristic large-scale conformational transitions of this protein family, and predicts an evolutionary conserved-and thus functional-homo-dimeric arrangement. Furthermore, we highlight that the features encoding the Hsp70 dimer are more conserved in bacterial than in eukaryotic sequences, suggesting that the known Hsp70/Hsp110 hetero-dimer is a eukaryotic specialization built on a pre-existing template.

Insights into the conformational switching mechanism of the human vascular endothelial growth factor receptor type 2 kinase domain.

Human vascular endothelial growth factor receptor type 2 (h-VEFGR2) is a receptor tyrosine kinase involved in the angiogenesis process and regarded as an interesting target for the design of anticancer drugs. Its activation/inactivation mechanism is related to conformational changes in its cytoplasmatic kinase domain, involving first among all the αC-helix in N-lobe and the A-loop in C-lobe. Affinity of inhibitors for the active or inactive kinase form could dictate the open or closed conformation of the A-loop, thus making the different conformations of the kinase domain receptor (KDR) domain different drug targets in drug discovery. In this view, a detailed knowledge of the conformational landscape of KDR domain is of central relevance to rationalize the efficiency and selectivity of kinase inhibitors. Here, molecular dynamics simulations were used to gain insight into the conformational switching activity of the KDR domain and to identify intermediate conformations between the two limiting active and inactive conformations. Specific energy barriers have been selectively removed to induce, and hence highlight at the atomistic level, the regulation mechanism of the A-loop opening. The proposed strategy allowed to repeatedly observe the escape of the KDR domain from the DFG-out free energy basin and to identify rare intermediate conformations between the DFG-out and the DFG-in structures to be employed in a structure-based drug discovery process.

ORAC: a molecular dynamics simulation program to explore free energy surfaces in biomolecular systems at the atomistic level.

We present the new release of the ORAC engine (Procacci et al., Comput Chem 1997, 18, 1834), a FORTRAN suite to simulate complex biosystems at the atomistic level. The previous release of the ORAC code included multiple time steps integration, smooth particle mesh Ewald method, constant pressure and constant temperature simulations. The present release has been supplemented with the most advanced techniques for enhanced sampling in atomistic systems including replica exchange with solute tempering, metadynamics and steered molecular dynamics. All these computational technologies have been implemented for parallel architectures using the standard MPI communication protocol. ORAC is an open-source program distributed free of charge under the GNU general public license (GPL) at http://www.chim.unifi.it/orac.

Fragment 101-108 of myelin oligodendrocyte glycoprotein: a possible lead compound for multiple sclerosis.

Multiple Sclerosis (MS) is a highly invalidating autoimmune disease of the central nervous system, leading to progressive paralysis and, sometimes, to premature death. One of the potential targets of the autoimmune reaction is the myelin protein MOG (Myelin Oligodendrocyte Glycoprotein). Since the 101-108 fragment of MOG plays a key role in the interaction with the MS-autoantibody 8-18C5, we performed an analysis of the equilibrium conformations of this peptide using the Replica Exchange Molecular Dynamics technique in conjunction with the Generalized Born continuum solvent model. Four variants of the peptide, stabilized by a disulfide bond, were also studied. We found that a significant fraction of the equilibrium population retains the original beta-hairpin conformation, and the amount of crystal-like conformations increases in the disulfide-closed analogues. When the equilibrium structures were used in docking simulations with the 8-18C5 autoantibody, we discovered the existence of a docking funnel whose bottom is populated by stable complexes where the peptide occupies the same region of space that was occupied in the crystal. It follows that the MOG 101-108 fragment represents a promising starting point for the design of a drug capable of blocking the 8-18C5 antibody. The molecule may also be used for the development of a diagnostic assay for multiple sclerosis.

Calculation of the potential of mean force from nonequilibrium measurements via maximum likelihood estimators.

We present an approach to the estimate of the potential of mean force along a generic reaction coordinate based on maximum likelihood methods and path-ensemble averages in systems driven far from equilibrium. Following similar arguments, various free energy estimators can be recovered, all providing comparable computational accuracy. The method, applied to the unfolding process of the alpha -helix form of an alanine decapeptide, gives results in good agreement with thermodynamic integration.

Generalization of the Jarzynski and Crooks nonequilibrium work theorems in molecular dynamics simulations.

The Jarzynski identity [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)] and the Crooks equation [G. E. Crooks, J. Stat. Phys. 90, 1481 (1998)] relate thermodynamic free energy differences to the work done on a system during a collection of nonequilibrium transformations. In the present Rapid Communication we provide generalized versions of these nonequilibrium work theorems, which hold for dissipative transformations where the system may undergo simultaneously mechanical work and pressure-temperature or volume-temperature changes. The proof is valid in the context of dynamic systems that evolve with NPT -based equations of motion according to the Martyna-Tobias-Klein algorithm [Martyna J. Chem. Phys. 101, 4177 (1994)]. An extension of the proof to dynamic systems that evolve through NVT -based equations of motion is also provided. The theorems may be effectively used in non-Hamiltonian molecular dynamics simulations for evaluating Helmholtz or Gibbs free energy differences, or the ratio of partition functions at different temperatures to be eventually used in thermodynamic cycles.

Self-healing umbrella sampling: a non-equilibrium approach for quantitative free energy calculations.

We propose a new approach for the umbrella sampling method in molecular dynamics simulations of complex systems. An accelerated sampling of the slow degrees of freedom is achieved by generating a single self-adaptive trajectory that tends to span uniformly the reaction coordinate using a time dependent bias potential derived from the preceding history of the system. To show the convergent behavior and the efficiency of the method, we present the free energy surface of alanine dipeptide in water as a function of the backbone dihedral angles.

Epigenomic Co-localization and Co-evolution Reveal a Key Role for 5hmC as a Communication Hub in the Chromatin Network of ESCs.

Epigenetic communication through histone and cytosine modifications is essential for gene regulation and cell identity. Here, we propose a framework that is based on a chromatin communication model to get insight on the function of epigenetic modifications in ESCs. The epigenetic communication network was inferred from genome-wide location data plus extensive manual annotation. Notably, we found that 5-hydroxymethylcytosine (5hmC) is the most-influential hub of this network, connecting DNA demethylation to nucleosome remodeling complexes and to key transcription factors of pluripotency. Moreover, an evolutionary analysis revealed a central role of 5hmC in the co-evolution of chromatin-related proteins. Further analysis of regions where 5hmC co-localizes with specific interactors shows that each interaction points to chromatin remodeling, stemness, differentiation, or metabolism. Our results highlight the importance of cytosine modifications in the epigenetic communication of ESCs.

Energy-Driven Undocking (EDU-HREM) in Solute Tempering Replica Exchange Simulations.

We present a new computational strategy for calculating the absolute binding free energy for protein ligand association in the context of atomistic simulation in explicit solvent. The method is based on an appropriate definition of a solute tempering scheme enforced via Hamilton replica exchange method (HREM). The definition of "solute" includes both the ligand and the active site, with the remainder of the systems defined as "solvent". The hydrophilicity of the solute and the solute torsional plus nonbonded intrasolute interactions are increased and decreased, respectively, along the replica progression, thus favoring the extrusion of the drug form the active site in the scaled states of the generalized ensemble. The proposed technique, named "Energy Driven Undocking" (EDU-HREM), completely bypasses the need for defining and/or identifying the relevant reaction coordinates in a ligand receptor interactions and allows the calculation of the absolute binding free energy in one single generalized simulation of the drug-receptor system. The methodology is applied, with encouraging results, to the calculation of the absolute binding free energy of some FK506-related ligands of the peptidyl prolyl cis-trans isomerase protein (FKBP12) with known dissociation constants. Aspects of the binding/inhibition mechanism in FKBP12 are also analyzed and discussed.

Intraligand hydrophobic interactions rationalize drug affinities for peptidyl-prolyl cis-trans isomerase protein.

The conformational landscape of three FK506-related drugs with disparate inhibition constants is determined in bulk solution using a replica exchange simulation method with solute torsional tempering. Energetic fitness of important drug conformations with respect to the FKBP12 protein is evaluated by molecular mechanics. Results show that the experimental affinity toward peptidyl-prolyl cis-trans isomerase protein (FKBP12) of the analyzed ligands appears to be positively correlated to the observed population of specific chair structures of the drug piperidinic ring in bulk solution. This observation is rationalized on the basis that such structures, stabilized by stereospecific intramolecular hydrophobic interactions, allows the formation of a pair of protein-ligand hydrogen bonds upon binding.

Conformational structure of the MOG-derived peptide 101-108 in solution.

One of the most important targets in the autoimmune attack in experimental autoimmune encephalomyielitis is the myelin oligodendrocyte glycoprotein (MOG). The complex with demyelinating 8-18C5 antibody was recently resolved by X-ray crystallography, showing a remarkable adhesion of the 101-108 MOG subsequence to the heavy chain of the autoantibody. In this study, we have determined, using replica exchange molecular dynamics methods, the structure of the MOG-derived peptide 101-108 in solution at ambient conditions. According to the simulation, the peptide exhibits, with significant probability, a distorted beta-turn structure highly similar to that of the corresponding subsequence in the crystal in complex with 8-18C5 antibody. Such results are found to be fully consistent with circular dichroism spectra of the peptide in solution, suggesting the use of the MOG-derived 101-108 peptide as a potential lead compound for designing decoy targets for the autoimmune attack in multiple sclerosis.

Free energy reconstruction in bidirectional force spectroscopy experiments: the effect of the device stiffness.

In force spectroscopy single-molecule experiments, an individual molecule, usually a polymer, is mechanically stretched by means of an externally controlled driving potential. Typically, the stiffness of this potential is much smaller than the stiffness of the potential of mean force along the molecular extension coordinate. Here we discuss how such a disparity alters the free energy and the reversibility of the driven system, with respect to the pristine molecular system under examination. In particular, by simulating unfolding/refolding experiments of a small protein, we examine the traits of the potential of mean force that are responsible for the dramatic amount of work dissipated in experiments using a soft device. Finally, we show that in bidirectional experiments the free energy of the free molecular system can be easily recovered by appropriate reweighting methods.

Thermodynamics of stacking interactions in proteins.

Using a database of 6166 experimental structures taken from the Protein Data Bank, we have studied pair interactions between planar residues (Phe, Tyr, His, Arg, Glu and Asp) in proteins, known as pi-pi interactions. On the basis of appropriate coordinates defining the mutual arrangement of two residues, we have calculated 2-D potentials of mean force aimed at determining the stability of the most probable structures for aromatic-aromatic, aromatic-cation and aromatic-anion bound pairs. Our analysis reveals the thermodynamic relevance and the ubiquity of stacked complexes in proteins.

Numerical verification of the generalized Crooks nonequilibrium work theorem for non-Hamiltonian molecular dynamics simulations.

The generalized Crooks theorem (GCT) for deterministic non-Hamiltonian molecular dynamics simulations [Phys. Rev. E 75, 050101 (2007)] connects the probabilities of nonequilibrium realizations switching the system between two thermodynamic states, to the partition functions of these states. In comparison to the "classical" Crooks nonequilibrium work theorem [J. Stat. Phys. 90, 1481 (1998)], which deals with realizations involving only mechanical work, the GCT also accounts for additional work resulting from changes of the intensive and extensive thermodynamic variables of the system. In this article we present a numerical verification of the GCT using a Lennard-Jones fluid model where two particles are subject to a time-dependent external potential. Moreover, in order to switch the system between different thermodynamic states, the temperature and the pressure (or volume), which are controlled through the Martyna-Tobias-Klein equations of motion [J. Chem. Phys. 101, 4177 (1994)], are also varied externally. The free energy difference between states characterized by different distances of the target particles is evaluated using both a standard methodology (pair radial distribution functions) and the GCT. In order to exploit the various options provided by the GCT approach, i.e., the possibility of temperature/pressure/volume changes during the realizations, the free energy difference is recovered via arbitrary thermodynamic cycles. In all tests, the GCT is quantitatively verified.

Recovering the Crooks equation for dynamical systems in the isothermal-isobaric ensemble: a strategy based on the equations of motion.

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)] relates the work done on a system during a nonequilibrium transformation to the free energy difference between the final and the initial state of the transformation. Recently, the authors have derived the Crooks equation for systems in the canonical ensemble thermostatted by the Nose-Hoover or Nose-Hoover chain method [P. Procacci et al., J. Chem. Phys. 125, 164101 (2006)]. That proof is essentially based on the fluctuation theorem by Evans and Searles [Adv. Phys. 51, 1529 (2002)] and on the equations of motion. Following an analogous approach, the authors derive here the Crooks equation in the context of molecular dynamics simulations of systems in the isothermal-isobaric (NPT) ensemble, whose dynamics is regulated by the Martyna-Tobias-Klein algorithm [J. Chem. Phys. 101, 4177 (1994)]. Their present derivation of the Crooks equation correlates to the demonstration of the Jarzynski identity for NPT systems recently proposed by Cuendet [J. Chem. Phys. 125, 144109 (2006)].

Crooks equation for steered molecular dynamics using a Nosé-Hoover thermostat.

The Crooks equation [Eq. (10) in J. Stat. Phys. 90, 1481 (1998)], originally derived for microscopically reversible Markovian systems, relates the work done on a system during an irreversible transformation to the free energy difference between the final and the initial state of the transformation. In the present work we provide a theoretical proof of the Crooks equation in the context of constant volume, constant temperature steered molecular dynamics simulations of systems thermostated by means of the Nosé-Hoover method (and its variant using a chain of thermostats). As a numerical test we use the folding and unfolding processes of decaalanine in vacuo at finite temperature. We show that the distribution of the irreversible work for the folding process is markedly non-Gaussian thereby implying, according to Crooks equation, that also the work distribution of the unfolding process must be inherently non-Gaussian. The clearly asymmetric behavior of the forward and backward irreversible work distributions is a signature of a non-Markovian regime for the folding/unfolding of decaalanine.