In Silico Identification of JMJD3 Demethylase Inhibitors.

In the search for new demethylase inhibitors, we have developed a multistep protocol for in silico screening. Millions of poses generated by high-throughput docking or a 3D-pharmacophore search are first minimized by a classical force field and then filtered by semiempirical quantum mechanical calculations of the interaction energy with a selected set of functional groups in the binding site. The final ranking includes solvation effects which are evaluated in the continuum dielectric approximation (finite-difference Poisson equation). Application of the multistep protocol to JMJD3 jumonji demethylase has resulted in a dozen low-micromolar inhibitors belonging to five different chemical classes. We have solved the crystal structure of JMJD3 inhibitor 8 in the complex with UTX (a demethylase in the same subfamily as JMJD3) which validates the predicted binding mode. Compound 8 is a promising candidate for future optimization as it has a favorable ligand efficiency of 0.32 kcal/mol per nonhydrogen atom.

CHARMM: the biomolecular simulation program.

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

ETNA: equilibrium transitions network and Arrhenius equation for extracting folding kinetics from REMD simulations.

It is difficult to investigate folding kinetics by conventional atomistic simulations of proteins. The replica exchange molecular dynamics (REMD) simulation technique enhances conformational sampling at the expenses of reduced kinetic information, which in REMD is directly available only for very short time scales. Here, we propose a procedure for obtaining kinetic data from REMD by making use of the equilibrium transitions network (ETN) sampled at the temperature of interest. This information is supplemented by mean folding times extracted from ETNs at higher REMD temperatures and scaled according to the Arrhenius equation. The procedure is applied to a three-stranded antiparallel beta-sheet peptide which has a very heterogeneous denatured state with a broad entropic basin and several enthalpic traps. Despite the complexity of the system and the REMD exchange time of only 0.1 ns, the procedure is able to estimate folding times (ranging from about 0.1 micros at the melting temperature of 330 K to about 8 micros at 286 K) as well as transition times from individual non-native basins to the native state.

Identification of the protein folding transition state from molecular dynamics trajectories.

The rate of protein folding is governed by the transition state so that a detailed characterization of its structure is essential for understanding the folding process. In vitro experiments have provided a coarse-grained description of the folding transition state ensemble (TSE) of small proteins. Atomistic details could be obtained by molecular dynamics (MD) simulations but it is not straightforward to extract the TSE directly from the MD trajectories, even for small peptides. Here, the structures in the TSE are isolated by the cut-based free-energy profile (cFEP) using the network whose nodes and links are configurations sampled by MD and direct transitions among them, respectively. The cFEP is a barrier-preserving projection that does not require arbitrarily chosen progress variables. First, a simple two-dimensional free-energy surface is used to illustrate the successful determination of the TSE by the cFEP approach and to explain the difficulty in defining boundary conditions of the Markov state model for an entropically stabilized free-energy minimum. The cFEP is then used to extract the TSE of a beta-sheet peptide with a complex free-energy surface containing multiple basins and an entropic region. In contrast, Markov state models with boundary conditions defined by projected variables and conventional histogram-based free-energy profiles are not able to identify the TSE of the beta-sheet peptide.

Iriomoteolides: novel chemical tools to study actin dynamics.

Despite its promising biological profile, the cellular targets of iriomoteolide-3a, a novel 15-membered macrolide isolated from Amphidinium sp., have remained unknown. A small library of non-natural iriomoteolide-3a analogues is presented here as a result of a novel, highly convergent, catalysis-based scaffold-diversification campaign, which revealed the suitable sites for chemical editing in the original core. We provide compelling experimental evidence for actin as one of iriomoteolides' primary cellular targets, establishing the ability of these secondary metabolites to inhibit cell migration, induce severe morphological changes in cells and cause a reversible cytoplasmic retraction and reduction of F-actin fibers in a time and dose dependent manner. These results are interpreted in light of the ability of iriomoteolides to stabilize F-actin filaments. Molecular dynamics simulations provide evidence for iriomoteolide-3a binding to the barbed end of G-actin. These results showcase iriomoteolides as novel and easily tunable chemical probes for the in vitro study of actin dynamics in the context of cell motility processes including cell invasion and division.

A molecular dynamics approach to the structural characterization of amyloid aggregation.

A novel computational approach to the structural analysis of ordered beta-aggregation is presented and validated on three known amyloidogenic polypeptides. The strategy is based on the decomposition of the sequence into overlapping stretches and equilibrium implicit solvent molecular dynamics (MD) simulations of an oligomeric system for each stretch. The structural stability of the in-register parallel aggregates sampled in the implicit solvent runs is further evaluated using explicit water simulations for a subset of the stretches. The beta-aggregation propensity along the sequence of the Alzheimer's amyloid-beta peptide (Abeta(42)) is found to be highly heterogeneous with a maximum in the segment V(12)HHQKLVFFAE(22) and minima at S(8)G(9), G(25)S(26), G(29)A(30), and G(38)V(39), which are turn-like segments. The simulation results suggest that these sites may play a crucial role in determining the aggregation tendency and the fibrillar structure of Abeta(42). Similar findings are obtained for the human amylin, a 37-residue peptide that displays a maximal beta-aggregation propensity at Q(10)RLANFLVHSSNN(22) and two turn-like sites at G(24)A(25) and G(33)S(34). In the third application, the MD approach is used to identify beta-aggregation "hot-spots" within the N-terminal domain of the yeast prion Ure2p (Ure2p(1-94)) and to design a double-point mutant (Ure2p-N4748S(1-94)) with lower beta-aggregation propensity. The change in the aggregation propensity of Ure2p-N4748S(1-94) is verified in vitro using the thioflavin T binding assay.

Structural details of urea binding to barnase: a molecular dynamics analysis.

BACKGROUND: The molecular mechanism of urea-induced protein unfolding has not been established. It is generally thought that denaturation results from the stabilizing interactions of urea with portions of the protein that are buried in the native state and become exposed upon unfolding of the protein. RESULTS: We have performed molecular dynamics simulations of barnase (a 110 amino acid RNase from Bacillus amyloliquefaciens) with explicit water and urea molecules at 300 K and 360 K. The native conformation was unaffected in the 300 K simulations at neutral and low pH. Two of the three runs at 360 K and low pH showed some denaturation, with partial unfolding of the hydrophobic core 2. The first solvation shell has a much higher density of urea molecules (water/urea ratio ranging from 2.07 to 2.73) than the bulk (water/urea ratio of 4.56). About one half of the first-shell urea molecules are involved in hydrogen bonds with polar or charged groups on the barnase surface, and between 15% and 18% of the first-shell urea molecules participate in multiple hydrogen bonds with barnase. The more stably bound urea molecules tend to be in crevices or pockets on the barnase surface. CONCLUSIONS: The simulation results indicate that an aqueous urea solution solvates the surface of a polypeptide chain more favorably than pure water. Urea molecules interact more favorably with nonpolar groups of the protein than water does, and the presence of urea improves the interactions of water molecules with the hydrophilic groups of the protein. The results suggest that urea denaturation involves effects on both nonpolar and polar groups of proteins.

Acid and thermal denaturation of barnase investigated by molecular dynamics simulations.

The transition in barnase from the native state to a partially unfolded conformation has been studied by molecular dynamics simulations with explicit water molecules at 360 K and low pH(450 ps), and at 600 K and neutral pH (three simulations of 120, 250 and 200 ps each). The use of several simulations provides evidence that the results are not sensitive to initial conditions. To mimic low pH conditions, the acidic sidechains in barnase were neutralized and the two histidine residues were doubly protonated. Runs at 300 K showed that the solvated structures at low pH (300 ps) and neutral pH (310 ps) are very similar. The main structural differences involved the acidic residues, histidine residues, and the beta-turn connecting strands 4 and 5. When the temperature is raised to 360 K at low pH and to 600 K at neutral pH the barnase molecule begins to unfold. The molecule rapidly expands (Rg changes from 13.9 A to 15.3 A in 450 ps at 360 K and from 13.7 A to between 15.1 and 15.5 A in 120 ps at 600 K). However, the expansion is not uniform. In all the simulations, the chain termini, loops and the N-terminal parts of the main alpha-helix (helix 1) show a continuous and progressive unfolding. An essential step in the denaturation process is that the major alpha-helix (helix 1) separates from the beta-sheet; this is coupled to the exposure of the principal hydrophobic core, many of whose non-polar side chains become solvated by hydrogen-bonded water molecules. The barnase-water interaction energy improves during unfolding at the expense of the barnase self-energy. The deterioration of the intramolecular van der Waals energy suggests that the rupture of the tight packing during the initial unfolding phase contributes to the energy barrier of the denaturation process. The mutationally well-analyzed Asp8-Arg110-Asp12 double salt-bridge on the barnase surface is found to be marginally stable in the folded form in the simulations. A Poisson-Boltzmann calculation indicates that the salt-bridge is unstable; this is probably due to an overestimate of the solvation energy. A detailed analysis of the main hydrophobic core reveals that increase in solvent-accessible surface area and penetration of water molecules are simultaneous in the high-temperature simulation; at lower temperatures there is significant cavity formation and the entrance of the water molecules is somewhat delayed. The cavities occur in the neighborhood of the hydrophobic sidechains; the region formed by the sidechains of Val10, Leu14, Leu20, Tyr24, Ala74, Ile76 and Tyr90 is involved. The loosening of the core packing is coupled to an increase in the number of dihedral transitions.(ABSTRACT TRUNCATED AT 400 WORDS)

Role of native topology investigated by multiple unfolding simulations of four SH3 domains.

The relative importance of amino acid sequence and native topology in the unfolding process of two SH3 domains and two circular permutants was investigated by 120 molecular dynamics runs at 375 K for a total simulation time of 0.72 micros. The alpha-spectrin (aSH3) and src SH3 (sSH3) domains, which have the same topology and a sequence identity of only 34%, show similar unfolding pathways. The disappearance of the three-stranded antiparallel beta-sheet is the last unfolding event, in agreement with a large repertoire of kinetic data derived from point mutations as well as glycine insertions and disulfide crosslinks. Two alternative routes of beta-sheet unfolding have emerged from the analysis of the trajectories. One is statistically preferred in aSH3 (n-src loop breaks before distal hair-pin) and the inverse in sSH3. An elongation of the beta2-beta3 hairpin was observed during the unfolding of sSH3 at 375 K and in 300 K simulations started from the putative transition state of sSH3 in accord with unusual kinetic data for point mutations at the n-src loop. The change of connectivity in the permutants influenced the sequence of unfolding events mainly at the permutation site. Regions where the connectivity remained unaffected showed the same chronology of contact disappearance. Taken together with previous folding simulations of two designed three-stranded antiparallel beta-sheet peptides, these results indicate that, at least for small beta-sheet proteins, the folding mechanism is primarily defined by the native state topology, whilst specific interactions determine the statistically predominant folding route.

Native topology or specific interactions: what is more important for protein folding?

Fifty-five molecular dynamics runs of two three-stranded antiparallel beta-sheet peptides were performed to investigate the relative importance of amino acid sequence and native topology. The two peptides consist of 20 residues each and have a sequence identity of 15 %. One peptide has Gly-Ser (GS) at both turns, while the other has d-Pro-Gly ((D)PG). The simulations successfully reproduce the NMR solution conformations, irrespective of the starting structure. The large number of folding events sampled along the trajectories at 360 K (total simulation time of about 5 micros) yield a projection of the free-energy landscape onto two significant progress variables. The two peptides have compact denatured states, similar free-energy surfaces, and folding pathways that involve the formation of a beta-hairpin followed by consolidation of the unstructured strand. For the GS peptide, there are 33 folding events that start by the formation of the 2-3 beta-hairpin and 17 with first the 1-2 beta-hairpin. For the (D)PG peptide, the statistical predominance is opposite, 16 and 47 folding events start from the 2-3 beta-hairpin and the 1-2 beta-hairpin, respectively. These simulation results indicate that the overall shape of the free-energy surface is defined primarily by the native-state topology, in agreement with an ever-increasing amount of experimental and theoretical evidence, while the amino acid sequence determines the statistically predominant order of the events.

Forces and energetics of hapten-antibody dissociation: a biased molecular dynamics simulation study.

The unbinding of fluorescein from the single-chain Fv fragment of the 4D5Flu antibody is investigated by biased molecular dynamics with an implicit solvation model. To obtain statistically meaningful results, a large number of unbinding trajectories are calculated; they involve a total simulation time of more than 200 ns. Simulations are carried out with a time-dependent perturbation and in the presence of a constant force. The two techniques, which provide complementary information, induce unbinding by favoring an increase in the distance between the ligand and the antibody. This distance is an appropriate progress variable for the dissociation reaction and permits direct comparison of the unbinding forces in the simulations with data from atomic force microscopy (AFM). The time-dependent perturbation generates unfolding pathways that are close to equilibrium and can be used to reconstruct the mean force; i.e. the derivative of the potential of mean force, along the reaction coordinate. This is supported by an analysis of the overall unbinding profile and the magnitude of the mean force, which are similar to those of the unbinding force (i.e. the external force due to the time-dependent perturbation) averaged over several unbinding events. The multiple simulations show that unbinding proceeds along a rather well-defined pathway for a broad range of effective pulling speeds. Initially, there is a distortion of the protein localized in the C-terminal region followed by the fluorescein exit from the binding site. This occurs in steps that involve breaking of specific electrostatic and van der Waals interactions. It appears that the simulations do not explore the same barriers as those measured in the AFM experiments because of the much higher unfolding speed in the former. The dependence of the force on the logarithm of the loading rate is linear and the slope is higher than in the AFM, in agreement with experiment in other systems, where different slopes were observed for different regimes. Based on the unbinding events, mutations in the 4D5Flu antigen binding site are predicted to result in significant changes in the unbinding force.

Multiple copy simultaneous search and construction of ligands in binding sites: application to inhibitors of HIV-1 aspartic proteinase.

Rational ligand design is a complex problem that can be divided into three parts: the search for optimal positions and orientations of functional groups in the binding site, the connection of such positions to form candidate ligands, and the estimation of their binding constants. Approaches for addressing the first two parts of the problem are described in the present work. They are applied to the construction of peptide ligands in the binding site of the human immunodeficiency virus 1 (HIV-1) proteinase. The primary objective is to test the method by comparison of the results with the MVT-101 complex structure for which coordinates are available; the results obtained with the liganded and unliganded proteinase structure are used to examine the utility of the latter for binding studies. A secondary objective is to show how to find new inhibitor candidates. The multiple copy simultaneous search (MCSS) method is utilized to search for optimal positions and orientations of a set of functional groups. For peptide ligands, functional groups corresponding to the protein main chain (N-methylacetamide) and to protein side chains (e.g., methanol, ethyl guanidinium) are used. The resulting N-methylacetamide minima are connected to form hexapeptide main chains with a simple pseudoenergy function that permits a complete search of all possible ways of connecting the minima. Side chains are added to the main-chain candidates by application of the same pseudoenergy function to the appropriate functional group minima. A set of 15 hexapeptides with the sequence of MVT-101 is then minimized by a Monte Carlo scheme, which allows for escape from local minima. Comparison of the MCSS results with the structure of MVT-101 in the HIV-1 binding site showed that all of its functional group positions correspond (within 2.4 A) to some (usually more than one) MCSS minima. There were also many other low-energy MCSS minima which do not appear in any known inhibitors, e.g., methyl ammonium minima in the neighborhood of the catalytic aspartates. Among the 15 lowest minima are seven hexapeptides with the same main-chain orientation as the one found by X-ray crystallography for the inhibitor MVT-101 in the binding site and eight with the main chain oriented in the opposite direction; the latter tend to be more stable. [Addendum: These results are in agreement with recent high-resolution crystallographic data provided after the study was completed.(ABSTRACT TRUNCATED AT 400 WORDS)

Computational ligand design.

A variety of computational tools that are used to assist drug design are reviewed. Particular emphasis is given to the limitations and merits of different methodologies. Recently, a number of general methods have been proposed for clustering compounds in classes of drug-like and non-drug-like molecules. The usefulness of this classification for drug design is discussed. The estimation of (relative) binding affinities is from a theoretical point of view the most challenging part of ligand design. We review three methods for the estimation of binding energies. Firstly, quantitative structure-activity relationships (QSAR) are presented. These have gained significantly from recent developments of experimental techniques for combinatorial synthesis and high-throughput screening as well as the use of powerful computational procedures like genetic algorithms and neural networks for the derivation of models. Secondly, empirical energy functions are shown to lead to more general models than standard QSAR, since they are fitted to a variety of complexes. They have been used recently with considerable success. Thirdly, we briefly outline free energy calculations based on molecular dynamics simulations, the method with the most sound theoretical foundation. Recent developments are reestablishing the interest in this approach. In the last part of this review structure-based ligand design programs are described. These are closely related to docking, with the difference that in design, unlike in most docking procedures, ligands are built on a fragment-by-fragment basis. Finally, a short description of our approach to computational combinatorial ligand design is given.

An evolutionary approach for structure-based design of natural and non-natural peptidic ligands.

A new computational approach (PEP) is presented for the structure-based design of linear polymeric ligands consisting of any type of amino acid. Ligands are grown from a seed by iteratively adding amino acids to the growing construct. The search in chemical space is performed by a build-up approach which employs all of the residues of a user-defined library. At every growing step, a genetic algorithm is used for conformational optimization of the last added monomer inside the binding site of a rigid target protein. The binding energy with electrostatic solvation is evaluated to select sequences for further growing. PEP is tested on the peptide substrate binding site of the insulin receptor tyrosine kinase and farnesyltransferase. In both test cases, the peptides designed by PEP correspond to the sequence motifs of known substrates. For tyrosine kinase, tyrosine residues are suggested at position P and P+2. While the tyrosine at P is in agreement with the experimental data, the one at P+2 is a prediction which awaits experimental validation. For farnesyltransferase, it is shown that electrostatic solvation is necessary for the correct design of peptidic inhibitors.

Solution conformation of phakellistatin 8 investigated by molecular dynamics simulations.

Phakellistatin 8 is a cyclic decapeptide that inhibits cancer cell growth and has sequence and structure similar to antamanide. In molecular dynamics simulations of phakellistatin 8 in water, the decapeptide ring undergoes a conformational change from the saddle-like crystal structure to a more elongated conformation by a transition of the Tyr9 main chain from the alpha L to an extended structure. This is coupled to the loss of the NH9-O6 beta-turn hydrogen bond and the transient dissociation of the Pro7-Tyr9 side-chain packing. Furthermore, the water molecule acting as a transannular bridge forms an additional hydrogen bond with phakellistatin 8, namely with the NH group of Val5 besides those already present in the crystal structure, i.e., with the NH of Ile10 and the CO of Leu6. The alpha-turn hydrogen bond between the Phe4 amide hydrogen and the Ile10 carbonyl oxygen is always present. The solution conformations of the two cyclic decapeptides are similar, in particular in the region involving the NH4-O10 alpha turn of phakellistatin 8 and the NH5-O1 alpha turn of antamanide. The simulation results suggest that in aqueous solution the conformation of phakellistatin 8 is more extended than in the crystalline state, and on a nanosecond time scale phakellistatin 8 is more flexible than antamanide.

Hydrophobicity and functionality maps of farnesyltransferase.

Farnesyltransferase (FTase) catalyzes the attachment of a 15-carbon isoprenoid moiety, farnesyl, through a thioether linkage to a cysteine near the C-terminus of oncogenic Ras proteins. These transform animal cells to a malignant phenotype when farnesylated. Hence, FTase is an interesting target for the development of antitumor agents. In this work we first investigate the active site of FTase by mapping its hydrophobic patches. Then the program SEED is used to dock functional groups into the active site by an exhaustive search and efficient evaluation of the binding energy with solvation. The electrostatic energy is SEED is based on the continuum dielectric approximation and consists of screened intermolecular energy and protein and fragment desolvation terms. The results are found to be consistent with the sequence variability of the tetrapeptide substrate. The distribution of functional groups (functionality maps) on the substrate binding site allows for identification of modifications of the tetrapeptide sequence that are consistent with potent peptidic inhibitors. Furthermore, the best minima of benzene match corresponding moieties of an inhibitor in clinical trials. The functionality maps are also used to design a library of disubstituted indoles that might prevent the binding of the protein substrates.

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.

Mechanism and Kinetics of Acetyl-Lysine Binding to Bromodomains.

Bromodomains are four-helix bundle proteins that specifically recognize acetylation of lysine side chains on histones. The available X-ray structures of bromodomain/histone tail complexes show that the conserved Asn residue in the loop between helices B and C is involved in a hydrogen bond with the acetyl-lysine side chain. Here we analyze the spontaneous binding of acetyl-lysine to the bromodomain TAF1(2) by the first molecular dynamics simulations of histone mark binding to an epigenetic reader protein. Multiple events of reversible association sampled along the unbiased simulations allow us to determine the pathway and kinetics of binding. The simulations show that acetyl-lysine has two major binding modes in TAF1(2) one of which corresponds to the available crystal structures and is stabilized by a hydrogen bond to the conserved Asn side chain. The other major binding mode is more buried than in the crystal structures and is stabilized by two hydrogen bonds with conserved residues of the loop between helices Z and A. In the more buried binding conformation, three of the six structured water molecules at the bottom of the binding pocket are displaced by the acetyl-lysine side chain. The kinetic analysis shows that the two binding modes interconvert on a faster time scale with respect to the association/dissociation process. The atomic-level description of the binding pathway and binding modes is useful for the design of small molecule modulators of histone binding to bromodomains.

Slow folding of cross-linked alpha-helical peptides due to steric hindrance.

The folding process of a 16-residue alpha-helical peptide with an azobenzene cross-linker (covalently bound to residues Cys3 and Cys14) is investigated by 50 molecular dynamics simulations of 4 micros each. The folding kinetics at 281 K show a stretched exponential behavior but become simpler and much faster when a distance restraint is used to emulate a nonbulky cross-linker. The free-energy basin of the helical state is divided into two subbasins by a barrier that separates helical conformations with opposite orientations of the Arg10 side chain with respect to the azobenzene cross-linker. In contrast, such barrier is not present in the helical basin of the peptide with the nonbulky cross-linker, which folds with speed similar to the unrestrained peptide. These results indicate that the cross-linker slows down folding because of steric hindrance rather than its restraining effect on the two ends of the helical segment.