Multiscale Approach for Computing Gated Ligand Binding from Molecular Dynamics and Brownian Dynamics Simulations.

We develop an approach to characterize the effects of gating by a multiconformation protein consisting of macrostate conformations that are either accessible or inaccessible to ligand binding. We first construct a Markov state model of the apo-protein from atomistic molecular dynamics simulations from which we identify macrostates and their conformations, compute their relative macrostate populations and interchange kinetics, and structurally characterize them in terms of ligand accessibility. We insert the calculated first-order rate constants for conformational transitions into a multistate gating theory from which we derive a gating factor γ that quantifies the degree of conformational gating. Applied to HIV-1 protease, our approach yields a kinetic network of three accessible (semi-open, open, and wide-open) and two inaccessible (closed and a newly identified, "parted") macrostate conformations. The parted conformation sterically partitions the active site, suggesting a possible role in product release. We find that the binding kinetics of drugs and drug-like inhibitors to HIV-1 protease falls in the slow gating regime. However, because γ = 0.75, conformational gating only modestly slows ligand binding. Brownian dynamics simulations of the diffusional association of eight inhibitors to the protease─having a wide range of experimental association constants (∼104-1010 M-1 s-1)─yields gated rate constants in the range of ∼0.5-5.7 × 108 M-1 s-1. This indicates that, whereas the association rate of some inhibitors could be described by the model, for many inhibitors either subsequent conformational transitions or alternate binding mechanisms may be rate-limiting. For systems known to be modulated by conformational gating, the approach could be scaled computationally efficiently to screen association kinetics for a large number of ligands.

The HIV-1 Nucleocapsid Regulates Its Own Condensation by Phase-Separated Activity-Enhancing Sequestration of the Viral Protease during Maturation.

A growing number of studies indicate that mRNAs and long ncRNAs can affect protein populations by assembling dynamic ribonucleoprotein (RNP) granules. These phase-separated molecular 'sponges', stabilized by quinary (transient and weak) interactions, control proteins involved in numerous biological functions. Retroviruses such as HIV-1 form by self-assembly when their genomic RNA (gRNA) traps Gag and GagPol polyprotein precursors. Infectivity requires extracellular budding of the particle followed by maturation, an ordered processing of ∼2400 Gag and ∼120 GagPol by the viral protease (PR). This leads to a condensed gRNA-NCp7 nucleocapsid and a CAp24-self-assembled capsid surrounding the RNP. The choreography by which all of these components dynamically interact during virus maturation is one of the missing milestones to fully depict the HIV life cycle. Here, we describe how HIV-1 has evolved a dynamic RNP granule with successive weak-strong-moderate quinary NC-gRNA networks during the sequential processing of the GagNC domain. We also reveal two palindromic RNA-binding triads on NC, KxxFxxQ and QxxFxxK, that provide quinary NC-gRNA interactions. Consequently, the nucleocapsid complex appears properly aggregated for capsid reassembly and reverse transcription, mandatory processes for viral infectivity. We show that PR is sequestered within this RNP and drives its maturation/condensation within minutes, this process being most effective at the end of budding. We anticipate such findings will stimulate further investigations of quinary interactions and emergent mechanisms in crowded environments throughout the wide and growing array of RNP granules.

Conformational Analysis of Heparin-Analogue Pentasaccharides by Nuclear Magnetic Resonance Spectroscopy and Molecular Dynamics Simulations.

Elucidation and improvement of the blood coagulant properties of heparin are the focus of intense research. In this study, we performed conformational analysis using nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics (MD) simulations on the heparin pentasaccharide analogue idraparinux, its disulfonatomethyl analogue, which features a slightly improved blood coagulation property, and a trisulfonatomethyl analogue, in which the activity has been totally abolished. As the ring conformation of the G subunit has been suggested as a major determinant of the biological properties, we analyzed the sugar ring conformations and dynamics of the interglycosidic linkages. We found that the conformation of the G ring is dominated by the 2SO skewed boat next to the 1C4 chair in all three derivatives. Both the thermodynamics and the kinetics of the conformational states were found to be highly similar in the three derivatives. Molecular kinetic analysis showed that the 2SO skewed boat state of the G ring is equally favorable in the three analogues, resulting in similar 2SO populations. Also, the transition kinetics from the 1C4 chair to the 2SO skewed boat was found to be comparable in the derivatives, which indicates a similar energy barrier between the two states of the G subunit. We also identified a slower conformational transition between the dominant 4C1 chair and the boat conformations on the E subunit. Both G and E ring flips are also accompanied by changes along the interglycosidic linkages, which take place highly synchronously with the ring flips. These findings indicate that conformational plasticity of the G ring and the dominance of the 2SO skewed boat populations do not necessarily warrant the biological activity of the derivatives and hence the impact of other factors also needs to be considered.

Inter-subunit interactions drive divergent dynamics in mammalian and Plasmodium actin filaments.

Cell motility is essential for protozoan and metazoan organisms and typically relies on the dynamic turnover of actin filaments. In metazoans, monomeric actin polymerises into usually long and stable filaments, while some protozoans form only short and highly dynamic actin filaments. These different dynamics are partly due to the different sets of actin regulatory proteins and partly due to the sequence of actin itself. Here we probe the interactions of actin subunits within divergent actin filaments using a comparative dynamic molecular model and explore their functions using Plasmodium, the protozoan causing malaria, and mouse melanoma derived B16-F1 cells as model systems. Parasite actin tagged to a fluorescent protein (FP) did not incorporate into mammalian actin filaments, and rabbit actin-FP did not incorporate into parasite actin filaments. However, exchanging the most divergent region of actin subdomain 3 allowed such reciprocal incorporation. The exchange of a single amino acid residue in subdomain 2 (N41H) of Plasmodium actin markedly improved incorporation into mammalian filaments. In the parasite, modification of most subunit-subunit interaction sites was lethal, whereas changes in actin subdomains 1 and 4 reduced efficient parasite motility and hence mosquito organ penetration. The strong penetration defects could be rescued by overexpression of the actin filament regulator coronin. Through these comparative approaches we identified an essential and common contributor, subdomain 3, which drives the differential dynamic behaviour of two highly divergent eukaryotic actins in motile cells.

Equilibrium Model of Drug-Modulated GagPol-Embedded HIV-1 Reverse Transcriptase Dimerization to Enhance Premature Protease Activation.

Lack of effective strategies for killing cells latently infected with HIV-1 limits the eradication of AIDS. Unfortunately, current antiretroviral inhibitors are designed to target virus production but not latent infection. Interestingly, some non-nucleoside reverse transcriptase inhibitors (NNRTIs) have shown off-design effects, specifically, premature activation of HIV-1 protease (PR) within virus-infected cells that induces apoptosis. Here, we analyze an equilibrium model of HIV-1 reverse transcriptase (RT) binding to NNRTIs to understand the optimal binding characteristics that enhance RT dimerization within embedded GagPol dimers. This would allow NNRTIs to act as PR autoactivation enhancers (PAEs). We compute that ∼700-fold enhancement is theoretically possible by PAEs. Both a strong drug-dimer binding affinity (KD12 < 100 nM) and relatively weaker drug-monomer affinity (KD2/KD12 > 10) are required for significant enhancement (∼50-fold or more) relative to the drug-free dimer concentration within a drug concentration limit of 10 µM. Our approach rationalizes the observed effects of efavirenz on premature activation of PR and may be useful to guide the design of suitable drug candidates and their optimal dosage regimens for this therapy class.

New approaches for computing ligand-receptor binding kinetics.

The recent and growing evidence that the efficacy of a drug can be correlated to target binding kinetics has seeded the development of a multitude of novel methods aimed at computing rate constants for receptor-ligand binding processes, as well as gaining an understanding of the binding and unbinding pathways and the determinants of structure-kinetic relationships. These new approaches include various types of enhanced sampling molecular dynamics simulations and the combination of energy-based models with chemometric analysis. We assess these approaches in the light of the varying levels of complexity of protein-ligand binding processes.

The Myxobacterial Metabolite Soraphen A Inhibits HIV-1 by Reducing Virus Production and Altering Virion Composition.

Soraphen A is a myxobacterial metabolite that blocks the acetyl-coenzyme A carboxylase of the host and was previously identified as a novel HIV inhibitor. Here, we report that soraphen A acts by reducing virus production and altering the gp120 virion content, impacting entry capacity and infectivity. These effects are partially reversed by addition of palmitic acid, suggesting that inhibition of HIV envelope palmitoylation is one of the mechanisms of antiviral action.

Kinetics for Drug Discovery: an industry-driven effort to target drug residence time.

A considerable number of approved drugs show non-equilibrium binding characteristics, emphasizing the potential role of drug residence times for in vivo efficacy. Therefore, a detailed understanding of the kinetics of association and dissociation of a target-ligand complex might provide crucial insight into the molecular mechanism-of-action of a compound. This deeper understanding will help to improve decision making in drug discovery, thus leading to a better selection of interesting compounds to be profiled further. In this review, we highlight the contributions of the Kinetics for Drug Discovery (K4DD) Consortium, which targets major open questions related to binding kinetics in an industry-driven public-private partnership.

Reaction-diffusion basis of retroviral infectivity.

Retrovirus particle (virion) infectivity requires diffusion and clustering of multiple transmembrane envelope proteins (Env3) on the virion exterior, yet is triggered by protease-dependent degradation of a partially occluding, membrane-bound Gag polyprotein lattice on the virion interior. The physical mechanism underlying such coupling is unclear and only indirectly accessible via experiment. Modelling stands to provide insight but the required spatio-temporal range far exceeds current accessibility by all-atom or even coarse-grained molecular dynamics simulations. Nor do such approaches account for chemical reactions, while conversely, reaction kinetics approaches handle neither diffusion nor clustering. Here, a recently developed multiscale approach is considered that applies an ultra-coarse-graining scheme to treat entire proteins at near-single particle resolution, but which also couples chemical reactions with diffusion and interactions. A model is developed of Env3 molecules embedded in a truncated Gag lattice composed of membrane-bound matrix proteins linked to capsid subunits, with freely diffusing protease molecules. Simulations suggest that in the presence of Gag but in the absence of lateral lattice-forming interactions, Env3 diffuses comparably to Gag-absent Env3 Initial immobility of Env3 is conferred through lateral caging by matrix trimers vertically coupled to the underlying hexameric capsid layer. Gag cleavage by protease vertically decouples the matrix and capsid layers, induces both matrix and Env3 diffusion, and permits Env3 clustering. Spreading across the entire membrane surface reduces crowding, in turn, enhancing the effect and promoting infectivity.This article is part of the themed issue 'Multiscale modelling at the physics-chemistry-biology interface'.

Computing the role of near attack conformations in an enzyme-catalyzed nucleophilic bimolecular reaction.

Near attack conformations (NACs) are conformations extending from the ground state (GS) that lie on the transition path of a chemical reaction. Here, we develop a method for computing the thermodynamic contribution to catalysis due to NAC formation in bimolecular reactions, within the limit of a classical molecular dynamics force field. We make use of the Bürgi­Dunitz theory applied to large-scale unbiased all-atom ensemble molecular dynamics simulations. We apply this to HIV-1 protease peptide hydrolysis, known to achieve a rate enhancement of ∼1011 (ΔGcat⧧ ∼ 15 kcal/mol) over the uncatalyzed bimolecular reaction (ΔGnon⧧ ∼ 30 kcal/mol). The ground state consists of a nucleophilic water molecule bound to an octapeptide substrate in the active site. We first observe multiple and reversible binding of a nucleophilic water molecule into the active site giving a free energy of binding of ΔG = −1 kcal/mol to form the GS. The free energy barriers for catalyzed and uncatalyzed NAC formation are both equivalent: ΔGNAC⧧ = 4.6 kcal/mol, constituting ∼30% and ∼15% of the overall barriers, respectively. Therefore, not only does adoption of NACs only account for minor progress along the transition path in both catalyzed and uncatalyzed reactions, but there is no preferential formation of them in the catalyzed reaction. Analysis of the catalytic hydrogen bond network reveals interactions that stabilize the GS; however, subsequent NAC formation does not preferentially favor any of the possible hydrogen bond configurations. This supports the view that the catalytic power of HIV-1 protease is not due to NAC formation.

Gag-Pol processing during HIV-1 virion maturation: a systems biology approach.

Proteolytic processing of Gag and Gag-Pol polyproteins by the viral protease (PR) is crucial for the production of infectious HIV-1, and inhibitors of the viral PR are an integral part of current antiretroviral therapy. The process has several layers of complexity (multiple cleavage sites and substrates; multiple enzyme forms; PR auto-processing), which calls for a systems level approach to identify key vulnerabilities and optimal treatment strategies. Here we present the first full reaction kinetics model of proteolytic processing by HIV-1 PR, taking into account all canonical cleavage sites within Gag and Gag-Pol, intermediate products and enzyme forms, enzyme dimerization, the initial auto-cleavage of full-length Gag-Pol as well as self-cleavage of PR. The model allows us to identify the rate limiting step of virion maturation and the parameters with the strongest effect on maturation kinetics. Using the modelling framework, we predict interactions and compensatory potential between individual cleavage rates and drugs, characterize the time course of the process, explain the steep dose response curves associated with PR inhibitors and gain new insights into drug action. While the results of the model are subject to limitations arising from the simplifying assumptions used and from the uncertainties in the parameter estimates, the developed framework provides an extendable open-access platform to incorporate new data and hypotheses in the future.

Rapid Conformational Fluctuations of Disordered HIV-1 Fusion Peptide in Solution.

The conformationally flexible fusion peptide (FP) of HIV-1 is indispensible for viral infection of host cells, due to its ability to insert into and tightly couple with phospholipid membranes. There are conflicting reports on the membrane-associated structure of FP, and solution structure information is limited, yet such a structure is the target for a novel class of antiretroviral inhibitors. An ensemble of explicit solvent molecular dynamics simulations, initiated from a disordered HIV-1 FP (aggregate time of ∼30 µs), revealed that while the vast majority of conformations predominantly lack secondary structure, both spontaneous formation and rapid interconversion of local secondary structure elements occur, highlighting the structural plasticity of the peptide. Therefore, even at this rapid time scale, FP constitutes a diverse and flexible conformational ensemble in solution. Secondary structure clustering reveals that the most prominent ordered elements are α- and 3-10-helical subsets of membrane-bound conformations, while trace populations within 2 Å RMSD of all complete membrane-bound conformations are found to pre-exist in the solution ensemble. Since inhibitor bound conformations of FP are only rarely found, FP inhibitors could function by modulating the conformational ensemble and binding to nonfusogenic FP structures. A thermodynamic characterization of the most prominent ordered nonfusogenic structures could facilitate the future design of improved FP inhibitors.

Kinetic characterization of the critical step in HIV-1 protease maturation.

HIV maturation requires multiple cleavage of long polyprotein chains into functional proteins that include the viral protease itself. Initial cleavage by the protease dimer occurs from within these precursors, and yet only a single protease monomer is embedded in each polyprotein chain. Self-activation has been proposed to start from a partially dimerized protease formed from monomers of different chains binding its own N termini by self-association to the active site, but a complete structural understanding of this critical step in HIV maturation is missing. Here, we captured the critical self-association of immature HIV-1 protease to its extended amino-terminal recognition motif using large-scale molecular dynamics simulations, thus confirming the postulated intramolecular mechanism in atomic detail. We show that self-association to a catalytically viable state requires structural cooperativity of the flexible ß-hairpin "flap" regions of the enzyme and that the major transition pathway is first via self-association in the semiopen/open enzyme states, followed by enzyme conformational transition into a catalytically viable closed state. Furthermore, partial N-terminal threading can play a role in self-association, whereas wide opening of the flaps in concert with self-association is not observed. We estimate the association rate constant (k(on)) to be on the order of ∼1 × 10(4) s(-1), suggesting that N-terminal self-association is not the rate-limiting step in the process. The shown mechanism also provides an interesting example of molecular conformational transitions along the association pathway.

Thumbs down for HIV: domain level rearrangements do occur in the NNRTI-bound HIV-1 reverse transcriptase.

One of the principal targets in human immunodeficiency virus type 1 (HIV-1) therapy is the reverse transcriptase (RT) enzyme. Non-nucleoside RT inhibitors (NNRTIs) are a class of highly specific drugs which bind to a pocket approximately 10 Šfrom the polymerase active site, inhibiting the enzyme allosterically. It is widely believed that NNRTIs function as "molecular wedges", disrupting the region between thumb and palm subdomains of the p66 subunit and locking the thumb in a wide-open conformation. Crystal structure data suggest that the binding of NNRTIs forces RT into a wide-open conformation in which the separation between the thumb and fingers subdomains is much higher than in the apo structure. Using ensemble molecular dynamics simulations (aggregate sampling ∼600 ns), we have captured RT bound to the NNRTI efavirenz in a closed conformation similar to that of the apo enzyme, suggesting the constraint of thumb motion is not as complete as previously believed. Rather, our investigation confirms that a conformational distribution across open and closed states must exist in the drug-bound enzyme and that allosteric modulation is effected via the alteration of the kinetic landscape of conformational transitions upon drug-binding. A more detailed understanding of the mechanism of NNRTI inhibition and the effect of binding upon domain motion could aid the design of more effective inhibitors and help identify novel allosteric sites.

Reaction kinetics of catalyzed competitive heteropolymer cleavage.

A theoretical formulation for complete heteropolymer degradation is developed in terms of Michaelis-Menten reaction kinetics under the quasi-steady-state approximation. This allows the concentration of the entire intermediate decomposition cascade to be accounted for as well as each species of emerging final product. The formulation is implemented computationally and results in stable reaction kinetics across a range of orders of magnitude for K(M) and k(cat). The model is compared with experiment, specifically in vitro HIV-1 protease-catalyzed retroviral Gag-polyprotein processing. Using an experimentally determined cleavage-polypeptide parameter set, good qualitative agreement is reached with Gag degradation kinetics, given the difference in experimental conditions. A parameter search within 1 order of magnitude of variation of the experimental set results in the determination of an optimal parameter set in complete agreement with experiment which allows the time evolution of each individual as well as intermediate species in Gag to be accurately followed. Future investigations that determine the required enzymatic parameters to populate such a scheme will allow for the model to be refined in order to track the time for viral maturation and infectivity.

Optimized Potential of Mean Force Calculations for Standard Binding Free Energies.

The prediction of protein-ligand binding free energies is an important goal of computational biochemistry, yet accuracy, reproducibility, and cost remain a problem. Nevertheless, these are essential requirements for computational methods to become standard binding prediction tools in discovery pipelines. Here, we present the results of an extensive search for an optimal method based on an ensemble of umbrella sampling all-atom molecular simulations tested on the phosphorylated tetrapeptide, pYEEI, binding to the SH2 domain, resulting in an accurate and converged binding free energy of -9.0 ± 0.5 kcal/mol (compared to an experimental value of -8.0 ± 0.1 kcal/mol). We find that a minimum of 300 ns of sampling is required for every prediction, a target easily achievable using new generation accelerated MD codes. Convergence is obtained by using an ensemble of simulations per window, each starting from different initial conformations, and by optimizing window-width, orthogonal restraints, reaction coordinate harmonic potentials, and window-sample time. The use of uncorrelated initial conformations in neighboring windows is important for correctly sampling conformational transitions from the unbound to bound states that affect significantly the precision of the calculations. This methodology thus provides a general recipe for reproducible and practical computations of binding free energies for a class of semirigid protein-ligand systems, within the limit of the accuracy of the force field used.

Explicit solvent dynamics and energetics of HIV-1 protease flap opening and closing.

An accurate description of the conformational dynamics of the ß-hairpin flaps of HIV-1 protease is of central importance in elucidating the functional recognition of the enzyme by ligands. Using all-atom molecular dynamics simulations in explicit solvent, with a total of 461 trajectories of ∼50 ns each, we report the closed, semiopen, open, and wide-open flap conformation of the free wild-type protease. The free energy of flap opening and closing from the semiopen state is 0.9 ± 0.2 and 2.4 ± 0.4 kcal/mol, respectively. The mean relaxation time of opening is ∼8 ns, in good agreement with NMR data. The explicit solvent simulations quantitatively confirm the hypothesis that the semiopen state is the dominant population in the free protease whilst fast flap tip fluctuations lead frequently to an open state. More pronounced flap rearrangements lead to a rare wide-open state with the catalytic site completely exposed to the solvent. The structures of the different flap conformations provided herein are of general interest for improved drug design of HIV-1 protease, in particular, the wide-open conformation could be favored by the large Gag and GagPol polyprotein chains. Strategies that take into account multiple flap-gating mechanisms may lead to more effective inhibitors.

Accurate ensemble molecular dynamics binding free energy ranking of multidrug-resistant HIV-1 proteases.

Accurate calculation of important thermodynamic properties, such as macromolecular binding free energies, is one of the principal goals of molecular dynamics simulations. However, single long simulation frequently produces incorrectly converged quantitative results due to inadequate sampling of conformational space in a feasible wall-clock time. Multiple short (ensemble) simulations have been shown to explore conformational space more effectively than single long simulations, but the two methods have not yet been thermodynamically compared. Here we show that, for end-state binding free energy determination methods, ensemble simulations exhibit significantly enhanced thermodynamic sampling over single long simulations and result in accurate and converged relative binding free energies that are reproducible to within 0.5 kcal/mol. Completely correct ranking is obtained for six HIV-1 protease variants bound to lopinavir with a correlation coefficient of 0.89 and a mean relative deviation from experiment of 0.9 kcal/mol. Multidrug resistance to lopinavir is enthalpically driven and increases through a decrease in the protein-ligand van der Waals interaction, principally due to the V82A/I84V mutation, and an increase in net electrostatic repulsion due to water-mediated disruption of protein-ligand interactions in the catalytic region. Furthermore, we correctly rank, to within 1 kcal/mol of experiment, the substantially increased chemical potency of lopinavir binding to the wild-type protease compared to saquinavir and show that lopinavir takes advantage of a decreased net electrostatic repulsion to confer enhanced binding. Our approach is dependent on the combined use of petascale computing resources and on an automated simulation workflow to attain the required level of sampling and turn around time to obtain the results, which can be as little as three days. This level of performance promotes integration of such methodology with clinical decision support systems for the optimization of patient-specific therapy.

Real science at the petascale.

We describe computational science research that uses petascale resources to achieve scientific results at unprecedented scales and resolution. The applications span a wide range of domains, from investigation of fundamental problems in turbulence through computational materials science research to biomedical applications at the forefront of HIV/AIDS research and cerebrovascular haemodynamics. This work was mainly performed on the US TeraGrid 'petascale' resource, Ranger, at Texas Advanced Computing Center, in the first half of 2008 when it was the largest computing system in the world available for open scientific research. We have sought to use this petascale supercomputer optimally across application domains and scales, exploiting the excellent parallel scaling performance found on up to at least 32 768 cores for certain of our codes in the so-called 'capability computing' category as well as high-throughput intermediate-scale jobs for ensemble simulations in the 32-512 core range. Furthermore, this activity provides evidence that conventional parallel programming with MPI should be successful at the petascale in the short to medium term. We also report on the parallel performance of some of our codes on up to 65 636 cores on the IBM Blue Gene/P system at the Argonne Leadership Computing Facility, which has recently been named the fastest supercomputer in the world for open science.

Automated molecular simulation based binding affinity calculator for ligand-bound HIV-1 proteases.

The successful application of high throughput molecular simulations to determine biochemical properties would be of great importance to the biomedical community if such simulations could be turned around in a clinically relevant timescale. An important example is the determination of antiretroviral inhibitor efficacy against varying strains of HIV through calculation of drug-protein binding affinities. We describe the Binding Affinity Calculator (BAC), a tool for the automated calculation of HIV-1 protease-ligand binding affinities. The tool employs fully atomistic molecular simulations alongside the well established molecular mechanics Poisson-Boltzmann solvent accessible surface area (MMPBSA) free energy methodology to enable the calculation of the binding free energy of several ligand-protease complexes, including all nine FDA approved inhibitors of HIV-1 protease and seven of the natural substrates cleaved by the protease. This enables the efficacy of these inhibitors to be ranked across several mutant strains of the protease relative to the wildtype. BAC is a tool that utilizes the power provided by a computational grid to automate all of the stages required to compute free energies of binding: model preparation, equilibration, simulation, postprocessing, and data-marshaling around the generally widely distributed compute resources utilized. Such automation enables the molecular dynamics methodology to be used in a high throughput manner not achievable by manual methods. This paper describes the architecture and workflow management of BAC and the function of each of its components. Given adequate compute resources, BAC can yield quantitative information regarding drug resistance at the molecular level within 96 h. Such a timescale is of direct clinical relevance and can assist in decision support for the assessment of patient-specific optimal drug treatment and the subsequent response to therapy for any given genotype.