A QM-MD simulation approach to the analysis of FRET processes in (bio)molecular systems. A case study: complexes of E. coli purine nucleoside phosphorylase and its mutants with formycin A.

Predicting FRET pathways in proteins using computer simulation techniques is very important for reliable interpretation of experimental data. A novel and relatively simple methodology has been developed and applied to purine nucleoside phosphorylase (PNP) complexed with a fluorescent ligand - formycin A (FA). FRET occurs between an excited Tyr residue (D*) and FA (A). This study aims to interpret experimental data that, among others, suggests the absence of FRET for the PNPF159A mutant in complex with FA, based on novel theoretical methodology. MD simulations for the protein molecule containing D*, and complexed with A, are carried out. Interactions of D* with its molecular environment are accounted by including changes of the ESP charges in S1, compared to S0, and computed at the SCF-CI level. FRET probability W F depends on the inverse six-power of the D*-A distance, R da . The orientational factor 0 < k(2) < 4 between D* and A is computed and included in the analysis. Finally W F is time-averaged over the MD trajectories resulting in its mean value. The red-shift of the tyrosinate anion emission and thus lack of spectral overlap integral and thermal energy dissipation are the reasons for the FRET absence in the studied mutants at pH 7 and above. The presence of the tyrosinate anion results in a competitive energy dissipation channel and red-shifted emission, thus in consequence in the absence of FRET. These studies also indicate an important role of the phenyl ring of Phe159 for FRET in the wild-type PNP, which does not exist in the Ala159 mutant, and for the effective association of PNP with FA. In a more general context, our observations point out very interesting and biologically important properties of the tyrosine residue in its excited state, which may undergo spontaneous deprotonation in the biomolecular systems, resulting further in unexpected physical and/or biological phenomena. Until now, this observation has not been widely discussed in the literature.

Prediction of secondary ionization of the phosphate group in phosphotyrosine peptides.

A computational approach, based on a continuum molecular electrostatics model, for the calculation of the pK(a) values of secondary ionization of the phosphate group in phenyl phosphate derivatives is described. The method uses the ESP atomic charges of the mono-anionic and di-anionic forms of the ionizable phosphate group, computed with the use of the density functional method, and applies a new concept of the model group, being the reference state for the pK(a) calculations. Both conformational flexibility and tautomeric degrees of freedom are taken into account in the calculations. The method was parameterized using experimentally available pK(a) values of four derivatives of phenyl phosphates, and phosphotyrosine. Subsequently this parameterization was used to predict pK(a) of the phosphate group in a short peptide Gly-Gly-Tyr(P)-Ala, and in a longer peptide consisting of 12 residues, the latter in water, and in a complex with a protein-phospholipase. The agreement between the computed and the experimental pK(a) values is better than +/-0.3 pH units for the optimized solute dielectric constant of 11-13. This approach is promising and its extension to other phospho-amino acids is in progress.

Structural basis of oncogenic activation caused by point mutations in the kinase domain of the MET proto-oncogene: modeling studies.

Missense mutations in the tyrosine kinase domain of the MET proto-oncogene occur in selected cases of papillary renal carcinoma. In biochemical and biological assays, these mutations produced constitutive activation of the MET kinase and led to tumor formation in nude mice. Some mutations caused transformation of NIH 3T3 cells. To elucidate the mechanism of ligand-independent MET kinase activation by point mutations, we constructed several 3D models of the wild-type and mutated MET catalytic core domains. Analysis of these structures showed that some mutations (e.g., V1110I, Y1248H/D/C, M1268T) directly alter contacts between residues from the activation loop in its inhibitory conformation and those from the main body of the catalytic domain; others (e.g., M1149T, L1213V) increase flexibility at the critical points of the tertiary structure and facilitate subdomain movements. Mutation D1246N plays a role in stabilizing the active form of the enzyme. Mutation M1268T affects the S+1 and S+3 substrate-binding pockets. Models implicate that although these changes do not compromise the affinity toward the C-terminal autophosphorylation site of the MET protein, they allow for binding of the substrate for the c-Abl tyrosine kinase. We provide biochemical data supporting this observation. Mutation L1213V affects the conformation of Tyr1212 in the active form of MET. Several somatic mutations are clustered at the surface of the catalytic domain in close vicinity of the probable location of the MET C-terminal docking site for cytoplasmic effectors.

Structure-based sequence alignment for the beta-trefoil subdomain of the clostridial neurotoxin family provides residue level information about the putative ganglioside binding site.

Clostridial neurotoxins embrace a family of extremely potent toxins comprised of tetanus toxin (TeNT) and seven different serotypes of botulinum toxin (BoNT/A-G). The beta-trefoil subdomain of the C-terminal part of the heavy chain (H(C)), responsible for ganglioside binding, is the most divergent region in clostridial neurotoxins with sequence identity as low as 15%. We re-examined the alignment between family sequences within this subdomain, since in this region all alignments published to date show obvious inconsistencies with the beta-trefoil fold. The final alignment was obtained by considering the general constraints imposed by this fold, and homology modeling studies based on the TeNT structure. Recently solved structures of BoNT/A confirm the validity of this structure-based approach. Taking into account biochemical data and crystal structures of TeNT and BoNT/A, we also re-examined the location of the putative ganglioside binding site and, using the new alignment, characterized this site in other BoNT serotypes.

Quantum-dynamical picture of a multistep enzymatic process: reaction catalyzed by phospholipase A(2).

A quantum-classical molecular dynamics model (QCMD), applying explicit integration of the time-dependent Schrödinger equation (QD) and Newtonian equations of motion (MD), is presented. The model is capable of describing quantum dynamical processes in complex biomolecular systems. It has been applied in simulations of a multistep catalytic process carried out by phospholipase A(2) in its active site. The process includes quantum-dynamical proton transfer from a water molecule to histidine localized in the active site, followed by a nucleophilic attack of the resulting OH(-) group on a carbonyl carbon atom of a phospholipid substrate, leading to cleavage of an adjacent ester bond. The process has been simulated using a parallel version of the QCMD code. The potential energy function for the active site is computed using an approximate valence bond (AVB) method. The dynamics of the key proton is described either by QD or classical MD. The coupling between the quantum proton and the classical atoms is accomplished via Hellmann-Feynman forces, as well as the time dependence of the potential energy function in the Schrödinger equation (QCMD/AVB model). Analysis of the simulation results with an Advanced Visualization System revealed a correlated rather than a stepwise picture of the enzymatic process. It is shown that an sp(2)--> sp(3) configurational change at the substrate carbonyl carbon is mostly responsible for triggering the activation process.

A simple model for predicting the free energy of binding between anthracycline antibiotics and DNA.

A theoretical model for predicting the free energy of binding between anthracycline antibiotics and DNA was developed using the electron density functional (DFT) and molecular mechanics (MM) methods. Partial DFT-ESP charges were used in calculating the MM binding energies for complexes formed between anthracycline antibiotics and oligodeoxynucleotides. These energies were then compared with experimental binding free energies. The good correlation between the experimental and theoretical energies allowed us to propose a model for predicting the binding free energy for derivatives of anthracycline antibiotics and for quickly screening new anthracycline derivatives.

A mezoscopic model of nucleic acids. Part 1. Lagrangian and quaternion molecular dynamics.

This study presents a model for mezoscopic molecular dynamics simulations with objects of different scale and properties e.e. atoms, pseudoatoms, rigid and pseudo-elastic bodies, described by the external coordinates and internal degrees of freedom. The Lagrangian approach is used to derive equations of motion and a quaternion representation is used for the description of the dynamics of rigid and pseudo-elastic molecular elements. Stability of the LQMD algorithm was tested for a 10-base pair deoxynucleotide. The total energy, momentum and angular momentum are conserved for time-steps up to 20 fs.

A mezoscopic model of nucleic acids. Part 2. An effective potential energy function for DNA.

In this study we present an effective Potential of Mean Force (PMF) designed for Lagrangian and Quaternion Molecular Dynamics (LQMD) of DNA. The DNA model is built from pseudoatoms as well as rigid and pseudo-elastic bodies described by a limited number of selected Cartesian and internal degrees of freedom. Phosphate groups, deoxyribose rings and nucleic acid bases are represented by pseudoparticles, some of them with internal degrees of freedom. PMF is defined as the sum of effective bonded and long-range potentials. The potentials were fitted to numerical free energy surfaces. Over 50 free energy surfaces, each depending on a conformational variable (pseudobond length, angle or dihedral angle) and the pseudorotation phase of the nearest neighbour deoxribose ring, were computed. The numerical free energy surfaces were obtained from probability distributions derived from a 1.5 ns conventional, microscopic MD simulation of the d(GpC)9 double helical DNA molecule. An umbrella sampling method was used to simulate transitions between the A and B DNA forms, and PMF reproduces these transitions.

Poisson-Boltzmann model studies of molecular electrostatic properties of the cAMP-dependent protein kinase.

Protonation equilibria of residues important in the catalytic mechanism of a protein kinase were analyzed on the basis of the Poisson-Boltzmann electrostatic model along with a cluster-based treatment of the multiple titration state problem. Calculations were based upon crystallographic structures of the mammalian cAMP-dependent protein kinase, one representing the so called closed form of the enzyme and the other representing an open conformation. It was predicted that at pH 7 the preferred form of the phosphate group at the catalytically essential threonine 197 (P-Thr197) in the closed form is dianionic, whereas in the open form a monoanionic ionization state is preferred. This dianionic state of P-Thr197, in the closed form, is stabilized by interactions with ionizable residues His87, Arg165, and Lys189. Our calculations predict that the hydroxyl of the Ser residue in the peptide substrate is very difficult to ionize, both in the closed and open structures of the complex. Also, the supposed catalytic base, Asp166, does not seem to have a pK(a) appropriate to remove the hydroxyl group proton of the peptide substrate. However, when Ser of the peptide substrate is forced to remain ionized, the predicted pK(a) of Asp166 increases strongly, which suggests that the Asp residue is a likely candidate to attract the proton if the Ser residue becomes deprotonated, possibly during some structural change preceding formation of the transition state. Finally, in accord with suggestions made on the basis of the pH-dependence of kinase kinetics, our calculations predict that Glu230 and His87 are the residues responsible for the molecular pK(a) values of 6.2 and 8.5, observed in the experiment.

Modelling of insulin receptor tyrosine kinase in its active form: a case study for validation of theoretical methods.

An active form of an insulin receptor tyrosine kinase (IRK) catalytic core was modelled based on its experimentally known inactive form and the active form of a serine/threonine kinase, protein kinase A (PKA). This theoretical model was compared with the crystallographic structure of the active form of IRK reported later. The structures are very similar, which shows that all the most important features and interactions have been taken into account in the modelling procedure. The elaborated procedure can be applied to other tyrosine kinases. This would allow designing of a wide class of tyrosine kinase inhibitors, very important potential anti-cancer and/or anti-viral drugs.

Binding of SPXK- and APXK-peptide motifs to AT-rich DNA. Experimental and theoretical studies.

The binding properties of the SPXK- and APXK-type peptides to the AT-rich DNA fragments of different length were studied by measuring the competition of peptides with Hoechst 33258 dye for DNA binding and by the gel shift assay analysis. In parallel to the experimental studies, molecular modeling techniques were used to analyze possible binding modes of the SPXZ and APXK motifs to the AT-rich DNA. The results of the competition measurements and gel shift assays suggest that serine at the i-1 position (i is proline) can be replaced by alanine without affecting the binding properties of the motif. Thus, the presence of the conserved serine in this motif in many DNA-binding proteins is probably not dictated by structural requirements. Based on the results of molecular modeling studies we propose that the binding mode of the SPXK-type motifs to the AT-rich DNA resembles closely that between the N-terminal arm of the homeodomain and DNA. This model confirms that serine in the SPXK motifs is not essential for the DNA binding. The model also indicates that if X in the motif is glutamic acid, this residue is probably protonated in the complex with DNA.

Modelling of active forms of protein kinases: p38--a case study.

An active form of p38 protein kinase, belonging to the mitogen-activated protein kinases subfamily, has been designed based on crystallographically known structures of two other kinases, an active form of protein kinase A (PKA) and an inactive form of extracellular signal-regulated kinase 2 (ERK2). The modelling procedure is described. Its general scheme can also be applied to other kinases. The structure of the active forms of p38 and PKA is very similar in the region which binds the substrate. The ATP-binding mode is very similar in the active forms of all the three studied kinases. Models of the active forms allow for further studies on transphosphorylation processes at the molecular level, and modelling of inhibitors competitive with ATP and/or substrates.

Lagrangian molecular dynamics using selected conformational degrees of freedom, with application to the pseudorotation dynamics of furanose rings.

Using internal conformational degrees of freedom for biopolymers as natural variables, and introducing a Lagrangian dynamics approach, one can simulate time-dependent processes over a much longer time scale than in classical Newtonian molecular dynamics (MD) techniques. Two factors contribute to this: a substantial reduction in the number of degrees of freedom and a very large increase in the size of the time step. We present the Lagrangian equations of motion for repuckering transitions in model furanose (F), ribose (R), and 2'-deoxyribose (dR) ring systems using the pseudorotation phase angle as the single dynamic variable. As in most Lagrangian analyses, the effective masses for the R and dR models are dependent on conformation, and we test the behavior of this variable mass (VM) model. Since the variation in effective mass is small, the VM model is compared with a simplified constant mass (CM) model, which is shown to be an excellent approximation. The equations of motion for the CM and VM models are integrated with the leapfrog and the iterative leapfrog algorithms, respectively. The Lagrangian dynamics approach reduces the number of degrees of freedom from about 40 to 1, and allows the use of time steps on the order of 20 fs, about an order of magnitude greater than is used in conventional MD simulations.

Molecular modeling methods. Basic techniques and challenging problems.

An overview is presented of computer modeling and simulation methods that play an increasing role in drug design: quantum chemical methods, molecular mechanics, molecular dynamics and Brownian dynamics. The application of molecular dynamics for the prediction of thermodynamic properties like free energy differences and binding constants is discussed. The Brownian dynamics method is presented in connection with the calculation of effective electrostatic forces using the Poisson-Boltzmann equation, which allows one to sample ligand-binding geometries and to predict the kinetics of diffusion-limited enzyme reactions. New techniques that have recently been extensively developed, such as the global energy minimization and quantum-classical dynamics methods, are also introduced. The molecular modeling methods are illustrated with selected examples.

Thermodynamic cycle-perturbation study of the binding of trifluoroacetyl dipeptide anilide inhibitors with porcine pancreatic elastase.

The variety of results of crystallographic studies of the serine proteases complexed with isocoumarin inhibitors presents a challenging problem to modeling methods and molecular energetics. Therefore, the thermodynamic cycle-perturbation technique has been used to study a model system of elastase and two peptidic inhibitors. Using the program AMBER, the technique correctly predicts changes of the binding constants for the trifluoroacetyl dipeptide inhibitors in comparison with available experimental (kinetic and crystallographic) data. However, the absolute values obtained are shown to be sensitive to the specific electrostatic interaction potential parameters used in the simulations. The reader and user are cautioned that thermodynamic cycle-perturbation results may be too optimistic by underestimating the accuracy of free energy values. This is especially a matter of concern for those cases where a direct comparison with experimental values is not possible, viz., (1) the stimulation of binding of novel compounds, (2) structurally uncertain binding sites, or (3) structurally different binding modes. With our best 4-31G* ESP (electrostatic potential) charges we were able to reproduce experimentally determined free energy differences (delta delta A) with an accuracy of about 1.5 kcal/mol. Dynamically induced structural changes in the binding site of elastase, and particularly changes in hydrogen-bond patterns of the binding site, are also reported.

Dynamic properties of the first enzymatic reaction steps of porcine pancreatic elastase. How rigid is the active site of the native enzyme? Molecular dynamics simulation.

Two molecular dynamics simulations (100 and 50 ps) of native porcine pancreatic elastase i.e., without bound substrate and with the active site hydrated by a dome of water (630 molecules) have been performed. Dynamical properties of the catalytic tetrad have been examined. While relative conformations of the Asp 102, His 57, and Ser 214 are rather stable in time, the side chain of Ser 195 undergoes several conformational changes. No preferences are observed for the formation of a hydrogen bond between the O gamma-H group (Ser 195) and nitrogen N, (His 57). A cluster of ordered water molecules effectively competes with the H-O gamma group (Ser 195) and thereby prevents the formation of this H bond, which is generally agreed to be crucial for catalysis.

A model for the hydrogen-bond-length probability distributions in the crystal structures of small-molecule components of the nucleic acids.

The probability distributions of the N-H...O = C and O-H...O = C hydrogen-bond lengths observed in the crystal structures of the purines, pyrimidines, nucleosides and nucleotides have been fitted to a one-dimensional hydrogen-bond potential-energy function. In order to obtain a quantitative correspondence between the experimental and theoretical distributions, it is necessary to include with the usual hydrogen-bond-type potential-energy function, an effective crystal-packing force and two thermodynamical parameters of the crystal lattice, the Debye temperature and the Gruneisen constant.

Energy minimization and molecular dynamics studies of Asn-102 elastase.

Four isomeric forms of the Asn-102 PPE (D102N mutant according to the emerging protocol, [Knowles, Science, 236 (1987) 1252-1258]) have been investigated using energy minimization (EM) and molecular dynamics (MD) techniques. MD simulation data for 175 ps are reported for each form (in total 700 ps for about 2500 atoms). The His-57 N epsilon-protonated forms are calculated to be more stable than the N delta-protonated ones. The active site region of the most stable form is very similar to that found in the D102N rat trypsin enzyme [Craik et al., Science, 237 (1987) 909-913]. Conformations of the active sites and their hydrogen bond patterns are presented for each of these forms and are compared with the structure of the native enzyme active site. The pH dependent activity of the D102N derivative is discussed.

Theoretical investigations on the conformation of 1,5,N(4),N(4)-tetramethylcytosine.

Theoretical investigations (Perturbative Configuration Interaction over Localized Orbitals (PCILO) and Intermediate Neglect of Differential Overlap (INDO) methods) of the conformation of tetramethylcystosine, an overcrowded molecule with planar structure, have been carried out. The physical features of rotational bending potentials of the dimethylamino group are discussed. Particularly, the controversial problem concerning the planarity of the molecule is investigated. The obtained results show that the planarity of tetramethylcytosine is an intrinsic property of the molecule. Nevertheless, because of the repulsion between the methyl groups, the planar structure of tetramethylcytosine is slightly destabilized. Further, the functional dependence of the dipole moment of tetramethylcytosine on rotation and bending of the dimethylamino group has been analyzed.