Self-consistent mean-field model based on molecular dynamics: application to lipid-cholesterol bilayers

We have developed a dynamic self-consistent mean-field model, based on molecular-dynamics simulations, to study lipid-cholesterol bilayers. In this model the lipid bilayer is represented as a two-dimensional lattice field in the lipid chain order parameters, while cholesterol molecules are represented by hard rods. The motion of rods in the system is continuous and is not confined to lattice cells. The statistical mechanics of chain ordering is described by a mean field derived from an extension of a model due to Marcelja. The time evolution of the system is governed by stochastic equations. The ensemble of chain configurations required in partition sums, and the energies of interaction, are taken from atomistic level molecular-dynamics simulations of lipid bilayers. The model allows us to simulate systems 500 nm in lateral size for 20 micros time scales, or greater. We have applied the model to dipalmitoyl-phosphatidylcholine-cholesterol (Chol) bilayers at 50 degrees C for Chol concentrations between 2% and 33%. At low concentrations of Chol (2%-4%), the model predicts the formation of isolated clusters of Chol surrounded by relatively ordered lipid chains, randomly dispersed in the disordered bilayer. With increasing Chol composition, regions of Chol-induced order begin to overlap. Starting from about 11% Chol this ordering effect becomes system wide and regions unaffected by Chol are no longer detectable. From the analysis of properties of the model we conclude that the change in lipid chain order with increasing Chol concentration is continuous over the 20-mus scale of the simulations. We also conclude that at 50 degrees C no large-scale Chol-rich and Chol-depleted coexisting phase-separated regions form at any concentration. At no point in any of the simulations do we observe a higher degree of lateral organization, such as Chol-based superlattice structures.     

3D-QSAR and molecular modeling of HIV-1 integrase inhibitors

Three-dimensional quantitative structure-activity relationship (3D QSAR) methods were applied on a series of inhibitors of HIV-1 integrase with respect to their inhibition of 3-processing and 3-end joining steps in vitro.The training set consisted of 27 compounds belonging to the class of thiazolothiazepines. The predictive ability of each model was evaluated using test set I consisting of four thiazolothiazepines and test set II comprised of seven compounds belonging to an entirely different structural class of coumarins. Maximum Common Substructure (MCS) based method was used to align the molecules and this was compared with other known methods of alignment. Two methods of 3D QSAR: comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) were analyzed in terms of their predictive abilities. CoMSIA produced significantly better results for all correlations. The results indicate a strong correlation between the inhibitory activity of these compounds and the steric and electrostatic fields around them. CoMSIA models with considerable internal as well as external predictive ability were obtained. A poor correlation obtained with hydrophobic field indicates that the binding of thiazolothiazepines to HIV-1 integrase is mainly enthalpic in nature. Further the most active compound of the series was docked into the active site using the crystal structure of integrase. The binding site was formed by the amino acid residues 64-67, 116, 148, 151-152, 155-156, and 159. The comparison of coefficient contour maps with the steric and electrostatic properties of the receptor shows high level of compatibility.     

Eigen value analysis of HIV-1 integrase inhibitors.

A three-dimensional quantitative structure activity relationship using the eigen value analysis (EVA) paradigm applied to 41 HIV-1 integrase inhibitors that inhibit integrase mediated cleavage (3'-processing step) and integration (3'-strand transfer step) in vitro was performed. The training set consisted of 35 molecules from five structurally diverse classes: salicylhydrazines, lichen acids, coumarins, quinones, and thiazolothiazepines. Models derived using semiempirical (MOPAC AM1 and PM3) calculated normal-mode frequencies were compared. The predictive ability of each resultant model was evaluated using a test set comprised of six molecules belonging to a different structural class: hydrazides. Models derived using AM1 method showed considerable internal as well as external predictivity (r(2)(cv) = 0.806, r(2)(pred) = 0.761 for 3'-processing and r(2)(cv) = 0.677, r(2)(pred) = 0.591 for 3'-strand transfer).     

Theory for trivial trajectory parallelization of multicanonical molecular dynamics and application to a polypeptide in water

Trivial trajectory parallelization of multicanonical molecular dynamics (TTP-McMD) explores the conformational space of a biological system with multiple short runs of McMD starting from various initial structures. This method simply connects (i.e., trivially parallelizes) the short trajectories and generates a long trajectory. First, we theoretically prove that the simple trajectory connection satisfies a detailed balance automatically. Thus, the resultant long trajectory is regarded as a single multicanonical trajectory. Second, we applied TTP-McMD to an alanine decapeptide with an all-atom model in explicit water to compute a free-energy landscape. The theory imposes two requirements on the multiple trajectories. We have demonstrated that TTP-McMD naturally satisfies the requirements. The TTP-McMD produces the free-energy landscape considerably faster than a single-run McMD does. We quantitatively showed that the accuracy of the computed landscape increases with increasing the number of multiple runs. Generally, the free-energy landscape of a large biological system is unknown a priori. The current method is suitable for conformational sampling of such a large system to reduce the waiting time to obtain a canonical ensemble statistically reliable.     

A simple clustering technique to improve QSAR model selection and predictivity: application to a receptor independent 4D-QSAR analysis of cyclic urea derived inhibitors of HIV-1 protease

A training set of 50 tetrahydropyrimidine-2-one based inhibitors of HIV-1 protease, for which the -log K(i) values were measured, was used to construct receptor independent 4D-QSAR models. A novel clustering technique was employed to facilitate and improve model selection as well as test set predictions. Following the manifold model theory, five unique models were chosen by the clustering algorithm (q(2) = 0.81-0.84). The models were used to map the atom type morphology of the inhibitor binding site of HIV-1 protease as well as to predict the potencies (-log K(i)) of 10 test set compounds. The rank-difference correlation coefficient was used to evaluate the quality of the test set predictions, which was improved from 0.39 to 0.68 when the clustering technique was applied. The set of five models, collectively, identify the important binding characteristics of the HIV protease receptor site. This study demonstrates that the selected simple clustering technique provides a discrete algorithm for model selection, as well as improving the quality of test set, or unknown, compound prediction as determined by the rank-difference correlation coefficient.     

Thermodynamic properties in the molecular dynamics ensemble applied to the gaussian core model fluid

The thermodynamic properties of pressure, energy, isothermal pressure coefficient, thermal expansion coefficient, isothermal and adiabatic compressibilities, isobaric and isochoric heat capacities, Joule-Thomson coefficient, and speed of sound are considered in a classical molecular dynamics ensemble. These properties were obtained using the treatment of Lustig [J. Chem. Phys. 100, 3048 (1994)] and Meier and Kabelac [J. Chem. Phys. 124, 064104 (2006)], whereby thermodynamic state variables are expressible in terms of phase-space functions determined directly from molecular dynamics simulations. The complete thermodynamic information about an equilibrium system can be obtained from this general formalism. We apply this method to the gaussian core model fluid because the complex phase behavior of this simple model provides a severe test for this treatment. Waterlike and other anomalies are observed for some of the thermodynamic properties of the gaussian core model fluid.     

Novel virtual screening protocol based on the combined use of molecular modeling and electron-ion interaction potential techniques to design HIV-1 integrase inhibitors

HIV-1 integrase (IN) is an essential enzyme for viral replication and represents an intriguing target for the development of new drugs. Although a large number of compounds have been reported to inhibit IN in biochemical assays, no drug active against this enzyme has been approved by the FDA so far. In this study, we report, for the first time, the use of the electron-ion interaction potential (EIIP) technique in combination with molecular modeling approaches for the identification of new IN inhibitors. An innovative virtual screening approach, based on the determination of both short- and long-range interactions between interacting molecules, was employed with the aim of identifying molecules able to inhibit the binding of IN to viral DNA. Moreover, results from a database screening on the commercial Asinex Gold Collection led to the selection of several compounds. One of them showed a significant inhibitory potency toward IN in the overall integration assay. Biological investigations also showed, in agreement with modeling studies, that these compounds prevent recognition of DNA by IN in a fluorescence fluctuation assay, probably by interacting with the DNA binding domain of IN.     

A structural study of model peptides derived from HIV-1 integrase central domain

The HIV-1 integrase (IN) catalyzes the integration of viral DNA in the human genome. In vitro the enzyme displays an equilibrium of monomers, dimers, tetramers and larger oligomers. However, its functional oligomeric form in vivo is not known. We report a study of the auto-associative properties of three peptides denoted K156, E156 and E159. These derive from the alpha4 helix of the IN catalytic core. The alpha4 helix is an amphipatic helix exposed at the surface of the protein and could be involved in the oligomerization process through its hydrophobic face. The peptides were obtained from the replacement of several amino acid residues by more helicogenic ones in the alpha4 helix peptide. K156 carries the basic residues Lys156 and Lys159, which have been shown important for the binding of IN to viral DNA. In E156 and E159 they are replaced with the acidic residue Glu. A fourth peptide K(E)156 obtained from the replacement of hydrophobic residues with Glu in K156 in order to abolish the auto-associative properties is used as a negative control. The capacity shown by peptides for alpha-helical formation is demonstrated by circular dichroism (CD) analysis performed in aqueous solution and in aqueous trifluoroethanol (TFE) mixtures. Both electrospray ionization mass spectrometry (ESI-MS) and glutaraldehyde chemical cross-linking show that peptides adopt different solvent-dependent equilibriums of monomers, dimers, trimers and tetramers. Oligomerization of peptides in aqueous solution is related to their ability to form helical structures. Addition of a small amount of TFE (<10%) stimulates helix stabilization and the interhelical hydrophobic contacts. Higher amounts of TFE alter the hydrophobic contacts and disrupt the oligomeric species. In addition to hydrophobic interactions, the patterns indicate that the biologically important Lys156 and Lys159 residues also participate in helix association. K(E)156 despite its ability to adopt a helical structure is unable to associate into oligomers, demonstrating the importance of hydrophobic contacts for oligomerization. Thus, the designed peptides provide us information on the functional properties of the alpha4 IN that seems to hold a dual role in DNA recognition and protein oligomerization.     

De novo design based identification of potential HIV-1 integrase inhibitors: A pharmacoinformatics study

In the present study, pharmacoinformatics paradigms include receptor-based de novo design, virtual screening through molecular docking and molecular dynamics (MD) simulation are implemented to identify novel and promising HIV-1 integrase inhibitors. The de novodrug/ligand/molecule design is a powerful and effective approach to design a large number of novel and structurally diverse compounds with the required pharmacological profiles. A crystal structure of HIV-1 integrase bound with standard inhibitor BI-224436 is used and a set of 80,000 compounds through the de novo approach in LigBuilder is designed. Initially, a number of criteria including molecular docking, in-silico toxicity and pharmacokinetics profile assessments are implied to reduce the chemical space. Finally, four de novo designed molecules are proposed as potential HIV-1 integrase inhibitors based on comparative analyses. Notably, strong binding interactions have been identified between a few newly identified catalytic amino acid residues and proposed HIV-1 integrase inhibitors. For evaluation of the dynamic stability of the protein-ligand complexes, a number of parameters are explored from the 100 ns MD simulation study. The MD simulation study suggested that proposed molecules efficiently retained their molecular interaction and structural integrity inside the HIV-1 integrase. The binding free energy is calculated through the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) approach for all complexes and it also explains their thermodynamic stability. Hence, proposed molecules through de novo design might be critical to inhibiting the HIV-1 integrase.     

A quantum mechanics/molecular mechanics study of the protein-ligand interaction for inhibitors of HIV-1 integrase

Human immunodeficiency virus type-1 integrase (HIV-1 IN) is an essential enzyme for effective viral replication. Diketo acids such as L-731,988 and S-1360 are potent and selective inhibitors of HIV-1 IN. In this study, we used molecular dynamics simulations, within the hybrid quantum mechanics/molecular mechanics (QM/MM) approach, to determine the protein-ligand interaction energy between HIV-1 IN and L-731,988 and 10 of its derivatives and analogues. This hybrid methodology has the advantage that it includes quantum effects such as ligand polarisation upon binding, which can be very important when highly polarisable groups are embedded in anisotropic environments, as for example in metal-containing active sites. Furthermore, an energy decomposition analysis was performed to determine the contributions of individual residues to the enzyme-inhibitor interactions on averaged structures obtained from rather extensive conformational sampling. Analysis of the results reveals first that there is a correlation between protein-ligand interaction energy and experimental strand transfer into human chromosomes and secondly that the Asn-155, Lys-156 and Lys-159 residues and the Mg(2+) ion are crucial to anti-HIV IN activity. These results may explain the available experimental data.     

Independent trajectory implementation of the semiclassical Liouville method: application to multidimensional reaction dynamics

We describe an independent trajectory implementation of semiclassical Liouville method for simulating quantum processes using classical trajectories. In this approach, a single ensemble of trajectories describes all semiclassical density matrix elements of a coupled electronic state problem, with the ensemble evolving classically under a single reference Hamiltonian chosen on the basis of physical grounds. In this paper, we introduce an additional uncoupled trajectory approximation, allowing the members of the ensemble to evolve independently of one another and eliminating the major computational costs of our previous coupled trajectory implementation. The accuracy of the method is demonstrated for model one-dimensional problems. In addition, the approach is applied to the chemical reaction dynamics of a collinear triatomic system, yielding excellent agreement with exact calculations. This method allows molecular dynamics involving coupled electronic surfaces to be modeled with essentially the same effort as classical molecular dynamics and ensemble averaging.     

An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.     

Coarse master equation from Bayesian analysis of replica molecular dynamics simulations

We use Bayesian inference to derive the rate coefficients of a coarse master equation from molecular dynamics simulations. Results from multiple short simulation trajectories are used to estimate propagators. A likelihood function constructed as a product of the propagators provides a posterior distribution of the free coefficients in the rate matrix determining the Markovian master equation. Extensions to non-Markovian dynamics are discussed, using the trajectory "paths" as observations. The Markovian approach is illustrated for the filling and emptying transitions of short carbon nanotubes dissolved in water. We show that accurate thermodynamic and kinetic properties, such as free energy surfaces and kinetic rate coefficients, can be computed from coarse master equations obtained through Bayesian inference.     

Wavelet transform analysis of ab initio molecular dynamics simulation: application to core-excitation dynamics of BF3

We propose a novel analysis method of ab initio molecular dynamics (AIMD) simulation using a continuous wavelet transform (c-WT) technique. The c-WT technique, one of the time-frequency signal analysis methods, provides a clear view of the dynamical information in time developments. Combined with the auto-correlation function of velocity by AIMD simulation, c-WT analysis enables us to well understand dynamical distribution, such as the vibrational properties following a change of electronic structure in a molecular system. As a practical application, AIMD simulation of core-excited BF(3) (B1s --> 2a(2) (')) is illustrated. AIMD simulation leads to the change of vibrational motion as well as structural deformation by core-excitation. The c-WT analysis clarifies the relationship between structural deformation and the related significant vibrational modes in core-excitation within 50 fs.     

Elastic bag model for molecular dynamics simulations of solvated systems: application to liquid argon

A new approach is developed to study the dynamics of the localized process in solutions and other condensed phase systems. The approach employs a fluctuating elastic boundary (FEB) model which encloses the simulated system in an elastic bag that mimics the effects of the bulk solvent. This alleviates the need for periodic boundary conditions and allows for a reduction in the number of solvent molecules that need to be included in the simulation. The boundary bag is modeled as a mesh of quasi-particles connected by elastic bonds. The FEB model allows for volume and density fluctuations characteristic of the bulk system, and the shape of the boundary fluctuates during the course of the simulation to adapt to the configuration fluctuations of the explicit solute-solvent system inside. The method is applied to the simulation of a Lennard-Jones model of liquid argon. Various structural and dynamical quantities are computed and compared with those obtained from conventional periodic boundary simulations. The agreement between the two is excellent in most cases, thus validating the viability of the FEB method.     

Insights into the dynamics of HIV-1 protease: a kinetic network model constructed from atomistic simulations

The conformational dynamics in the flaps of HIV-1 protease plays a crucial role in the mechanism of substrate binding. We develop a kinetic network model, constructed from detailed atomistic simulations, to determine the kinetic mechanisms of the conformational transitions in HIV-1 PR. To overcome the time scale limitation of conventional molecular dynamics (MD) simulations, our method combines replica exchange MD with transition path theory (TPT) to study the diversity and temperature dependence of the pathways connecting functionally important states of the protease. At low temperatures the large-scale flap opening is dominated by a small number of paths; at elevated temperatures the transition occurs through many structurally heterogeneous routes. The expanded conformation in the crystal structure 1TW7 is found to closely mimic a key intermediate in the flap-opening pathways at low temperature. We investigated the different transition mechanisms between the semi-open and closed forms. The calculated relaxation times reveal fast semi-open ? closed transitions, and infrequently the flaps fully open. The ligand binding rate predicted from this kinetic model increases by 38-fold from 285 to 309 K, which is in general agreement with experiments. To our knowledge, this is the first application of a network model constructed from atomistic simulations together with TPT to analyze conformational changes between different functional states of a natively folded protein.