Search

The COVID-19 pandemic continues to pose a substantial threat to human lives and is likely to do so for years to come. Despite the availability of vaccines, searching for efficient small-molecule drugs that are widely available, including in low- and middle-income countries, is an ongoing challenge. In this work, we report the results of an open science community effort, the "Billion molecules against COVID-19 challenge", to identify small-molecule inhibitors against SARS-CoV-2 or relevant human receptors. Participating teams used a wide variety of computational methods to screen a minimum of 1 billion virtual molecules against 6 protein targets. Overall, 31 teams participated, and they suggested a total of 639,024 molecules, which were subsequently ranked to find 'consensus compounds'. The organizing team coordinated with various contract research organizations (CROs) and collaborating institutions to synthesize and test 878 compounds for biological activity against proteases (Nsp5, Nsp3, TMPRSS2), nucleocapsid N, RdRP (only the Nsp12 domain), and (alpha) spike protein S. Overall, 27 compounds with weak inhibition/binding were experimentally identified by binding-, cleavage-, and/or viral suppression assays and are presented here. Open science approaches such as the one presented here contribute to the knowledge base of future drug discovery efforts in finding better SARS-CoV-2 treatments., (© 2023 The Authors. Molecular Informatics published by Wiley-VCH GmbH.)

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 (∼10 4 -10 10 M -1 s -1 )─yields gated rate constants in the range of ∼0.5-5.7 × 10 8 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.

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.

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 2 S O skewed boat next to the 1 C 4 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 2 S O skewed boat state of the G ring is equally favorable in the three analogues, resulting in similar 2 S O populations. Also, the transition kinetics from the 1 C 4 chair to the 2 S O 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 4 C 1 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 2 S O 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.

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 (K D 12  < 100 nM) and relatively weaker drug-monomer affinity (K D 2 /K D 12 > 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.

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., Competing Interests: The authors have declared that no competing interests exist.

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., (Copyright © 2017 Elsevier Ltd. All rights reserved.)

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., (Copyright © 2017 American Society for Microbiology.)

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., (Copyright © 2017. Published by Elsevier Ltd.)

Retrovirus particle (virion) infectivity requires diffusion and clustering of multiple transmembrane envelope proteins (Env 3 ) 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 Env 3 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, Env 3 diffuses comparably to Gag-absent Env 3 Initial immobility of Env 3 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 Env 3 diffusion, and permits Env 3 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'., (© 2016 The Author(s).)

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.

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.

Functioning of G protein-coupled receptors (GPCRs) is tightly linked to the membrane environment, but a molecular level understanding of the modulation of GPCR by membrane lipids is not available. However, specific receptor-lipid interactions as well as unspecific effects mediated by the bulk properties of the membrane (thickness, curvature, etc.) have been proposed to be key regulators of GPCR modulation. In this review, we examine computational efforts made towards modeling and simulation of (i) the complex behavior of membrane lipids, (ii) membrane lipid-GPCR interactions as well as membrane lipid-mediated effects on GPCRs and (iii) GPCR oligomerization in a native-like membrane environment. We propose that, from the perspective of computational modeling, all three of these components need to be addressed in order to achieve a deeper understanding of GPCR functioning. Presently, we are able to simulate numerous lipid properties applying advanced computational techniques, although some barriers, such as the time-length of these simulations, need to be overcome. Implementing three-dimensional structures of GPCRs in such validated membrane systems can give novel insights in membrane-dependent receptor modulation and formation of higher order receptor complexes. Finally, more realistic GPCR-membrane models would provide a very useful tool in studying receptor behavior and its modulation by small drug-like ligands, a relevant issue for drug discovery.

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.

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.

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.

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.

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.

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., (Proteins 2010. © 2010 Wiley-Liss, Inc.)

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.

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.

Patient-specific medical simulation holds the promise of determining tailored medical treatment based on the characteristics of an individual patient (for example, using a genotypic assay of a sequence of DNA). Decision-support systems based on patient-specific simulation can potentially revolutionize the way that clinicians plan courses of treatment for various conditions, ranging from viral infections to arterial abnormalities. Basing medical decisions on the results of simulations that use models derived from data specific to the patient in question means that the effectiveness of a range of potential treatments can be assessed before they are actually administered, preventing the patient from experiencing unnecessary or ineffective treatments. We illustrate the potential for patient-specific simulation by first discussing the scale of the evolving international grid infrastructure that is now available to underpin such applications. We then consider two case studies, one concerned with the treatment of patients with HIV/AIDS and the other addressing neuropathologies associated with the intracranial vasculature. Such patient-specific medical simulations require access to both appropriate patient data and the computational resources on which to perform potentially very large simulations. Computational infrastructure providers need to furnish access to a wide range of different types of resource, typically made available through heterogeneous computational grids, and to institute policies that facilitate the performance of patient-specific simulations on those resources. To support these kinds of simulations, where life and death decisions are being made, computational resource providers must give urgent priority to such jobs, for example by allowing them to pre-empt the queue on a machine and run straight away. We describe systems that enable such priority computing.

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.

To explain drug resistance by computer simulations at the molecular level, we first have to assess the accuracy of theoretical predictions. Herein we report an application of the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) technique to the ranking of binding affinities of the inhibitor saquinavir with the wild type (WT) and three resistant mutants of HIV-1 protease: L90M, G48V, and G48V/L90M. For each ligand-protein complex we report 10 ns of fully unrestrained molecular dynamics (MD) simulations with explicit solvent. We investigate convergence, internal consistency, and model dependency of MM/PBSA ligand binding energies. Converged enthalpy and entropy estimates produce ligand binding affinities within 1.5 kcal/mol of experimental values, with a remarkable level of correlation to the experimentally observed ranking of resistance levels. A detailed analysis of the enthalpic/entropic balance of drug-protease interactions explains resistance in L90M in terms of a higher vibrational entropy than in the WT complex, while G48V disrupts critical hydrogen bonds at the inhibitor's binding site and produces an altered, more unfavorable balance of Coulomb and polar desolvation energies.

We provide insight into the first stages of a kinetic mechanism of lateral drug expulsion from the active site of HIV-1 protease, by conducting all atom molecular dynamics simulations with explicit solvent over a time scale of 24 ns for saquinavir bound to the wildtype, G48V, L90M and G48V/L90M mutant proteases. We find a consistent escape mechanism associated with the G48V mutation. First, increased hydrophilic and hydrophobic flap coupling and water mediated disruption of catalytic dyad hydrogen bonding induce drug motion away from the dyad and promote protease flap transition to the semi-open form. Conversely, flap-inhibitor motion is decoupled in the wildtype. Second, the decrease of total interactions causes unidirectional lateral inhibitor translation by up to 4 A toward the P3 subsite exit of the active site, increased P3 subsite exposure to solvent and a complete loss of hydrophobic interactions with the opposite end of the active site. The P1 subsite moves beyond the hydrophobic active site side pocket, the only remaining steric barrier to complete expulsion being the "breathable" residue, P81. Significant inhibitor deviation is reported over 24 ns, and subsequent complete expulsion, implemented using steered molecular dynamics simulations, is shown to occur most easily for the G48V-containing mutants. Our simulations thus provide compelling support for lateral drug escape from a protease in a semi-open flap conformation. It is likely that some mutations take advantage of this escape mechanism to increase the rate of inhibitor dissociation from the protease. Finally, unidirectional translation may be countered by designing inhibitors with terminal subsites that provide sufficient anchoring to the flaps, thus increasing the steric barrier for translation in either direction.