Revealing the mechanism of action of a first-in-class covalent inhibitor of KRASG12C (ON) and other functional properties of oncogenic KRAS by 31P NMR.

Individual oncogenic KRAS mutants confer distinct differences in biochemical properties and signaling for reasons that are not well understood. KRAS activity is closely coupled to protein dynamics and is regulated through two interconverting conformations: state 1 (inactive, effector binding deficient) and state 2 (active, effector binding enabled). Here, we use 31P NMR to delineate the differences in state 1 and state 2 populations present in WT and common KRAS oncogenic mutants (G12C, G12D, G12V, G13D, and Q61L) bound to its natural substrate GTP or a commonly used nonhydrolyzable analog GppNHp (guanosine-5'-[(ß,γ)-imido] triphosphate). Our results show that GppNHp-bound proteins exhibit significant state 1 population, whereas GTP-bound KRAS is primarily (90% or more) in state 2 conformation. This observation suggests that the predominance of state 1 shown here and in other studies is related to GppNHp and is most likely nonexistent in cells. We characterize the impact of this differential conformational equilibrium of oncogenic KRAS on RAF1 kinase effector RAS-binding domain and intrinsic hydrolysis. Through a KRAS G12C drug discovery, we have identified a novel small-molecule inhibitor, BBO-8956, which is effective against both GDP- and GTP-bound KRAS G12C. We show that binding of this inhibitor significantly perturbs state 1-state 2 equilibrium and induces an inactive state 1 conformation in GTP-bound KRAS G12C. In the presence of BBO-8956, RAF1-RAS-binding domain is unable to induce a signaling competent state 2 conformation within the ternary complex, demonstrating the mechanism of action for this novel and active-conformation inhibitor.

Machine learning-driven multiscale modeling reveals lipid-dependent dynamics of RAS signaling proteins.

RAS is a signaling protein associated with the cell membrane that is mutated in up to 30% of human cancers. RAS signaling has been proposed to be regulated by dynamic heterogeneity of the cell membrane. Investigating such a mechanism requires near-atomistic detail at macroscopic temporal and spatial scales, which is not possible with conventional computational or experimental techniques. We demonstrate here a multiscale simulation infrastructure that uses machine learning to create a scale-bridging ensemble of over 100,000 simulations of active wild-type KRAS on a complex, asymmetric membrane. Initialized and validated with experimental data (including a new structure of active wild-type KRAS), these simulations represent a substantial advance in the ability to characterize RAS-membrane biology. We report distinctive patterns of local lipid composition that correlate with interfacially promiscuous RAS multimerization. These lipid fingerprints are coupled to RAS dynamics, predicted to influence effector binding, and therefore may be a mechanism for regulating cell signaling cascades.

Clustering Protein Binding Pockets and Identifying Potential Drug Interactions: A Novel Ligand-Based Featurization Method.

Protein-ligand interactions are essential to drug discovery and drug development efforts. Desirable on-target or multitarget interactions are the first step in finding an effective therapeutic, while undesirable off-target interactions are the first step in assessing safety. In this work, we introduce a novel ligand-based featurization and mapping of human protein pockets to identify closely related protein targets and to project novel drugs into a hybrid protein-ligand feature space to identify their likely protein interactions. Using structure-based template matches from PDB, protein pockets are featured by the ligands that bind to their best co-complex template matches. The simplicity and interpretability of this approach provide a granular characterization of the human proteome at the protein-pocket level instead of the traditional protein-level characterization by family, function, or pathway. We demonstrate the power of this featurization method by clustering a subset of the human proteome and evaluating the predicted cluster associations of over 7000 compounds.

OFraMP: a fragment-based tool to facilitate the parametrization of large molecules.

An Online tool for Fragment-based Molecule Parametrization (OFraMP) is described. OFraMP is a web application for assigning atomic interaction parameters to large molecules by matching sub-fragments within the target molecule to equivalent sub-fragments within the Automated Topology Builder (ATB, atb.uq.edu.au) database. OFraMP identifies and compares alternative molecular fragments from the ATB database, which contains over 890,000 pre-parameterized molecules, using a novel hierarchical matching procedure. Atoms are considered within the context of an extended local environment (buffer region) with the degree of similarity between an atom in the target molecule and that in the proposed match controlled by varying the size of the buffer region. Adjacent matching atoms are combined into progressively larger matched sub-structures. The user then selects the most appropriate match. OFraMP also allows users to manually alter interaction parameters and automates the submission of missing substructures to the ATB in order to generate parameters for atoms in environments not represented in the existing database. The utility of OFraMP is illustrated using the anti-cancer agent paclitaxel and a dendrimer used in organic semiconductor devices. OFraMP applied to paclitaxel (ATB ID 35922).

Exploring CRD mobility during RAS/RAF engagement at the membrane.

During the activation of mitogen-activated protein kinase (MAPK) signaling, the RAS-binding domain (RBD) and cysteine-rich domain (CRD) of RAF bind to active RAS at the plasma membrane. The orientation of RAS at the membrane may be critical for formation of the RAS-RBDCRD complex and subsequent signaling. To explore how RAS membrane orientation relates to the protein dynamics within the RAS-RBDCRD complex, we perform multiscale coarse-grained and all-atom molecular dynamics (MD) simulations of KRAS4b bound to the RBD and CRD domains of RAF-1, both in solution and anchored to a model plasma membrane. Solution MD simulations describe dynamic KRAS4b-CRD conformations, suggesting that the CRD has sufficient flexibility in this environment to substantially change its binding interface with KRAS4b. In contrast, when the ternary complex is anchored to the membrane, the mobility of the CRD relative to KRAS4b is restricted, resulting in fewer distinct KRAS4b-CRD conformations. These simulations implicate membrane orientations of the ternary complex that are consistent with NMR measurements. While a crystal structure-like conformation is observed in both solution and membrane simulations, a particular intermolecular rearrangement of the ternary complex is observed only when it is anchored to the membrane. This configuration emerges when the CRD hydrophobic loops are inserted into the membrane and helices α3-5 of KRAS4b are solvent exposed. This membrane-specific configuration is stabilized by KRAS4b-CRD contacts that are not observed in the crystal structure. These results suggest modulatory interplay between the CRD and plasma membrane that correlate with RAS/RAF complex structure and dynamics, and potentially influence subsequent steps in the activation of MAPK signaling.

Improved Protein-Ligand Binding Affinity Prediction with Structure-Based Deep Fusion Inference.

Predicting accurate protein-ligand binding affinities is an important task in drug discovery but remains a challenge even with computationally expensive biophysics-based energy scoring methods and state-of-the-art deep learning approaches. Despite the recent advances in the application of deep convolutional and graph neural network-based approaches, it remains unclear what the relative advantages of each approach are and how they compare with physics-based methodologies that have found more mainstream success in virtual screening pipelines. We present fusion models that combine features and inference from complementary representations to improve binding affinity prediction. This, to our knowledge, is the first comprehensive study that uses a common series of evaluations to directly compare the performance of three-dimensional (3D)-convolutional neural networks (3D-CNNs), spatial graph neural networks (SG-CNNs), and their fusion. We use temporal and structure-based splits to assess performance on novel protein targets. To test the practical applicability of our models, we examine their performance in cases that assume that the crystal structure is not available. In these cases, binding free energies are predicted using docking pose coordinates as the inputs to each model. In addition, we compare these deep learning approaches to predictions based on docking scores and molecular mechanic/generalized Born surface area (MM/GBSA) calculations. Our results show that the fusion models make more accurate predictions than their constituent neural network models as well as docking scoring and MM/GBSA rescoring, with the benefit of greater computational efficiency than the MM/GBSA method. Finally, we provide the code to reproduce our results and the parameter files of the trained models used in this work. The software is available as open source at https://github.com/llnl/fast. Model parameter files are available at ftp://gdo-bioinformatics.ucllnl.org/fast/pdbbind2016_model_checkpoints/.

Binding Affinity Prediction by Pairwise Function Based on Neural Network.

We present a new approach to estimate the binding affinity from given three-dimensional poses of protein-ligand complexes. In this scheme, every protein-ligand atom pair makes an additive free-energy contribution. The sum of these pairwise contributions then gives the total binding free energy or the logarithm of the dissociation constant. The pairwise contribution is calculated by a function implemented via a neural network that takes the properties of the two atoms and their distance as input. The pairwise function is trained using a portion of the PDBbind 2018 data set. The model achieves good accuracy for affinity predictions when evaluated with PDBbind 2018 and with the CASF-2016 benchmark, comparing favorably to many scoring functions such as that of AutoDock Vina. The framework here may be extended to incorporate other factors to further improve its accuracy and power.

Predicting Small Molecule Transfer Free Energies by Combining Molecular Dynamics Simulations and Deep Learning.

Accurately predicting small molecule partitioning and hydrophobicity is critical in the drug discovery process. There are many heterogeneous chemical environments within a cell and entire human body. For example, drugs must be able to cross the hydrophobic cellular membrane to reach their intracellular targets, and hydrophobicity is an important driving force for drug-protein binding. Atomistic molecular dynamics (MD) simulations are routinely used to calculate free energies of small molecules binding to proteins, crossing lipid membranes, and solvation but are computationally expensive. Machine learning (ML) and empirical methods are also used throughout drug discovery but rely on experimental data, limiting the domain of applicability. We present atomistic MD simulations calculating 15,000 small molecule free energies of transfer from water to cyclohexane. This large data set is used to train ML models that predict the free energies of transfer. We show that a spatial graph neural network model achieves the highest accuracy, followed closely by a 3D-convolutional neural network, and shallow learning based on the chemical fingerprint is significantly less accurate. A mean absolute error of ∼4 kJ/mol compared to the MD calculations was achieved for our best ML model. We also show that including data from the MD simulation improves the predictions, tests the transferability of each model to a diverse set of molecules, and show multitask learning improves the predictions. This work provides insight into the hydrophobicity of small molecules and ML cheminformatics modeling, and our data set will be useful for designing and testing future ML cheminformatics methods.

ddcMD: A fully GPU-accelerated molecular dynamics program for the Martini force field.

We have implemented the Martini force field within Lawrence Livermore National Laboratory's molecular dynamics program, ddcMD. The program is extended to a heterogeneous programming model so that it can exploit graphics processing unit (GPU) accelerators. In addition to the Martini force field being ported to the GPU, the entire integration step, including thermostat, barostat, and constraint solver, is ported as well, which speeds up the simulations to 278-fold using one GPU vs one central processing unit (CPU) core. A benchmark study is performed with several test cases, comparing ddcMD and GROMACS Martini simulations. The average performance of ddcMD for a protein-lipid simulation system of 136k particles achieves 1.04 µs/day on one NVIDIA V100 GPU and aggregates 6.19 µs/day on one Summit node with six GPUs. The GPU implementation in ddcMD offloads all computations to the GPU and only requires one CPU core per simulation to manage the inputs and outputs, freeing up remaining CPU resources on the compute node for alternative tasks often required in complex simulation campaigns. The ddcMD code has been made open source and is available on GitHub at https://github.com/LLNL/ddcMD.

In Silico Insights into Protein-protein Interaction Disruptive Mutations in the PCSK9-LDLR complex.

Gain-of-function mutations in PCSK9 (proprotein convertase subtilisin/kexin type 9) lead to reduced uptake of LDL (low density lipoprotein) cholesterol and, therefore, increased plasma LDL levels. However, the mechanism by which these mutants reduce LDL reuptake is not fully understood. Here, we have used molecular dynamics simulations, MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) binding affinity calculations, and residue interaction networks, to investigate the protein-protein interaction (PPI) disruptive effects of two of PCSK9's gain-of-function mutations, Ser127Arg and Asp374Tyr on the PCSK9 and LDL receptor complex. In addition to these PPI disruptive mutants, a third, non-interface mutation (Arg496Trp) is included as a positive control. Our results indicate that Ser127Arg and Asp374Tyr confer significantly improved binding affinity, as well as different binding modes, when compared to the wild-type. These PPI disruptive mutations lie between the EGF(A) (epidermal growth factor precursor homology domain A) of the LDL receptor and the catalytic domain of PCSK9 (Asp374Tyr) and between the prodomain of PCSK9 and the ß-propeller of the LDL receptor (Ser127Arg). The interactions involved in these two interfaces result in an LDL receptor that is sterically inhibited from entering its closed conformation. This could potentially implicate the prodomain as a target for small molecule inhibitors.

Molecular Mechanism for Gramicidin Dimerization and Dissociation in Bilayers of Different Thickness.

Membrane protein functions can be altered by subtle changes in the host lipid bilayer physical properties. Gramicidin channels have emerged as a powerful system for elucidating the underlying mechanisms of membrane protein function regulation through changes in bilayer properties, which are reflected in the thermodynamic equilibrium distribution between nonconducting gramicidin monomers and conducting bilayer-spanning dimers. To improve our understanding of how subtle changes in bilayer thickness alter the gramicidin monomer and dimer distributions, we performed extensive atomistic molecular dynamics simulations and fluorescence-quenching experiments on gramicidin A (gA). The free-energy calculations predicted a nonlinear coupling between the bilayer thickness and channel formation. The energetic barrier inhibiting gA channel formation was sharply increased in the thickest bilayer (1,2-dierucoyl-sn-glycero-3-phosphocholine). This prediction was corroborated by experimental results on gramicidin channel activity in bilayers of different thickness. To further explore the mechanism of channel formation, we performed extensive unbiased molecular dynamics simulations, which allowed us to observe spontaneous gA dimer formation in lipid bilayers. The simulations revealed structural rearrangements in the gA subunits and changes in lipid packing, as well as water reorganization, that occur during the dimerization process. Together, the simulations and experiments provide new, to our knowledge, insights into the process and mechanism of gramicidin channel formation, as a prototypical example of the bilayer regulation of membrane protein function.

Structural Control of Nonnative Ligand Binding in Engineered Mutants of Phosphoenolpyruvate Carboxykinase.

Protein engineering to alter recognition underlying ligand binding and activity has enormous potential. Here, ligand binding for Escherichia coli phosphoenolpyruvate carboxykinase (PEPCK), which converts oxaloacetate into CO2 and phosphoenolpyruvate as the first committed step in gluconeogenesis, was engineered to accommodate alternative ligands as an exemplary system with structural information. From our identification of bicarbonate binding in the PEPCK active site at the supposed CO2 binding site, we probed binding of nonnative ligands with three oxygen atoms arranged to resemble the bicarbonate geometry. Crystal structures of PEPCK and point mutants with bound nonnative ligands thiosulfate and methanesulfonate along with strained ATP and reoriented oxaloacetate intermediates and unexpected bicarbonate were determined and analyzed. The mutations successfully altered the bound ligand position and orientation and its specificity: mutated PEPCKs bound either thiosulfate or methanesulfonate but never both. Computational calculations predicted a methanesulfonate binding mutant and revealed that release of the active site ordered solvent exerts a strong influence on ligand binding. Besides nonnative ligand binding, one mutant altered the Mn2+ coordination sphere: instead of the canonical octahedral ligand arrangement, the mutant in question had an only five-coordinate arrangement. From this work, critical features of ligand binding, position, and metal ion cofactor geometry required for all downstream events can be engineered with small numbers of mutations to provide insights into fundamental underpinnings of protein-ligand recognition. Through structural and computational knowledge, the combination of designed and random mutations aids in the robust design of predetermined changes to ligand binding and activity to engineer protein function.

Computational Lipidomics of the Neuronal Plasma Membrane.

Membrane lipid composition varies greatly within submembrane compartments, different organelle membranes, and also between cells of different cell stage, cell and tissue types, and organisms. Environmental factors (such as diet) also influence membrane composition. The membrane lipid composition is tightly regulated by the cell, maintaining a homeostasis that, if disrupted, can impair cell function and lead to disease. This is especially pronounced in the brain, where defects in lipid regulation are linked to various neurological diseases. The tightly regulated diversity raises questions on how complex changes in composition affect overall bilayer properties, dynamics, and lipid organization of cellular membranes. Here, we utilize recent advances in computational power and molecular dynamics force fields to develop and test a realistically complex human brain plasma membrane (PM) lipid model and extend previous work on an idealized, "average" mammalian PM. The PMs showed both striking similarities, despite significantly different lipid composition, and interesting differences. The main differences in composition (higher cholesterol concentration and increased tail unsaturation in brain PM) appear to have opposite, yet complementary, influences on many bilayer properties. Both mixtures exhibit a range of dynamic lipid lateral inhomogeneities ("domains"). The domains can be small and transient or larger and more persistent and can correlate between the leaflets depending on lipid mixture, Brain or Average, as well as on the extent of bilayer undulations.

An Electrostatic Funnel in the GABA-Binding Pathway.

The γ-aminobutyric acid type A receptor (GABAA-R) is a major inhibitory neuroreceptor that is activated by the binding of GABA. The structure of the GABAA-R is well characterized, and many of the binding site residues have been identified. However, most of these residues are obscured behind the C-loop that acts as a cover to the binding site. Thus, the mechanism by which the GABA molecule recognizes the binding site, and the pathway it takes to enter the binding site are both unclear. Through the completion and detailed analysis of 100 short, unbiased, independent molecular dynamics simulations, we have investigated this phenomenon of GABA entering the binding site. In each system, GABA was placed quasi-randomly near the binding site of a GABAA-R homology model, and atomistic simulations were carried out to observe the behavior of the GABA molecules. GABA fully entered the binding site in 19 of the 100 simulations. The pathway taken by these molecules was consistent and non-random; the GABA molecules approach the binding site from below, before passing up behind the C-loop and into the binding site. This binding pathway is driven by long-range electrostatic interactions, whereby the electrostatic field acts as a 'funnel' that sweeps the GABA molecules towards the binding site, at which point more specific atomic interactions take over. These findings define a nuanced mechanism whereby the GABAA-R uses the general zwitterionic features of the GABA molecule to identify a potential ligand some 2 nm away from the binding site.

GABAA receptor target of tetramethylenedisulfotetramine.

Use of the highly toxic and easily prepared rodenticide tetramethylenedisulfotetramine (TETS) was banned after thousands of accidental or intentional human poisonings, but it is of continued concern as a chemical threat agent. TETS is a noncompetitive blocker of the GABA type A receptor (GABAAR), but its molecular interaction has not been directly established for lack of a suitable radioligand to localize the binding site. We synthesized [(14)C]TETS (14 mCi/mmol, radiochemical purity >99%) by reacting sulfamide with H(14)CHO and s-trioxane then completion of the sequential cyclization with excess HCHO. The outstanding radiocarbon sensitivity of accelerator mass spectrometry (AMS) allowed the use of [(14)C]TETS in neuroreceptor binding studies with rat brain membranes in comparison with the standard GABAAR radioligand 4'-ethynyl-4-n-[(3)H]propylbicycloorthobenzoate ([(3)H]EBOB) (46 Ci/mmol), illustrating the use of AMS for characterizing the binding sites of high-affinity (14)C radioligands. Fourteen noncompetitive antagonists of widely diverse chemotypes assayed at 1 or 10 µM inhibited [(14)C]TETS and [(3)H]EBOB binding to a similar extent (r(2) = 0.71). Molecular dynamics simulations of these 14 toxicants in the pore region of the α1ß2γ2 GABAAR predict unique and significant polar interactions for TETS with α1T1' and γ2S2', which are not observed for EBOB or the GABAergic insecticides. Several GABAAR modulators similarly inhibited [(14)C]TETS and [(3)H]EBOB binding, including midazolam, flurazepam, avermectin Ba1, baclofen, isoguvacine, and propofol, at 1 or 10 µM, providing an in vitro system for recognizing candidate antidotes.

Catalytic site identification--a web server to identify catalytic site structural matches throughout PDB.

The catalytic site identification web server provides the innovative capability to find structural matches to a user-specified catalytic site among all Protein Data Bank proteins rapidly (in less than a minute). The server also can examine a user-specified protein structure or model to identify structural matches to a library of catalytic sites. Finally, the server provides a database of pre-calculated matches between all Protein Data Bank proteins and the library of catalytic sites. The database has been used to derive a set of hypothesized novel enzymatic function annotations. In all cases, matches and putative binding sites (protein structure and surfaces) can be visualized interactively online. The website can be accessed at http://catsid.llnl.gov.

A method to predict blood-brain barrier permeability of drug-like compounds using molecular dynamics simulations.

The blood-brain barrier (BBB) is formed by specialized tight junctions between endothelial cells that line brain capillaries to create a highly selective barrier between the brain and the rest of the body. A major problem to overcome in drug design is the ability of the compound in question to cross the BBB. Neuroactive drugs are required to cross the BBB to function. Conversely, drugs that target other parts of the body ideally should not cross the BBB to avoid possible psychotropic side effects. Thus, the task of predicting the BBB permeability of new compounds is of great importance. Two gold-standard experimental measures of BBB permeability are logBB (the concentration of drug in the brain divided by concentration in the blood) and logPS (permeability surface-area product). Both methods are time-consuming and expensive, and although logPS is considered the more informative measure, it is lower throughput and more resource intensive. With continual increases in computer power and improvements in molecular simulations, in silico methods may provide viable alternatives. Computational predictions of these two parameters for a sample of 12 small molecule compounds were performed. The potential of mean force for each compound through a 1,2-dioleoyl-sn-glycero-3-phosphocholine bilayer is determined by molecular dynamics simulations. This system setup is often used as a simple BBB mimetic. Additionally, one-dimensional position-dependent diffusion coefficients are calculated from the molecular dynamics trajectories. The diffusion coefficient is combined with the free energy landscape to calculate the effective permeability (Peff) for each sample compound. The relative values of these permeabilities are compared to experimentally determined logBB and logPS values. Our computational predictions correlate remarkably well with both logBB (R(2) = 0.94) and logPS (R(2) = 0.90). Thus, we have demonstrated that this approach may have the potential to provide reliable, quantitatively predictive BBB permeability, using a relatively quick, inexpensive method.

Toward fully automated high performance computing drug discovery: a massively parallel virtual screening pipeline for docking and molecular mechanics/generalized Born surface area rescoring to improve enrichment.

In this work we announce and evaluate a high throughput virtual screening pipeline for in-silico screening of virtual compound databases using high performance computing (HPC). Notable features of this pipeline are an automated receptor preparation scheme with unsupervised binding site identification. The pipeline includes receptor/target preparation, ligand preparation, VinaLC docking calculation, and molecular mechanics/generalized Born surface area (MM/GBSA) rescoring using the GB model by Onufriev and co-workers [J. Chem. Theory Comput. 2007, 3, 156-169]. Furthermore, we leverage HPC resources to perform an unprecedented, comprehensive evaluation of MM/GBSA rescoring when applied to the DUD-E data set (Directory of Useful Decoys: Enhanced), in which we selected 38 protein targets and a total of ∼0.7 million actives and decoys. The computer wall time for virtual screening has been reduced drastically on HPC machines, which increases the feasibility of extremely large ligand database screening with more accurate methods. HPC resources allowed us to rescore 20 poses per compound and evaluate the optimal number of poses to rescore. We find that keeping 5-10 poses is a good compromise between accuracy and computational expense. Overall the results demonstrate that MM/GBSA rescoring has higher average receiver operating characteristic (ROC) area under curve (AUC) values and consistently better early recovery of actives than Vina docking alone. Specifically, the enrichment performance is target-dependent. MM/GBSA rescoring significantly out performs Vina docking for the folate enzymes, kinases, and several other enzymes. The more accurate energy function and solvation terms of the MM/GBSA method allow MM/GBSA to achieve better enrichment, but the rescoring is still limited by the docking method to generate the poses with the correct binding modes.

Membrane lipids drive formation of KRAS4b-RAF1 RBDCRD nanoclusters on the membrane.

The oncogene RAS, extensively studied for decades, presents persistent gaps in understanding, hindering the development of effective therapeutic strategies due to a lack of precise details on how RAS initiates MAPK signaling with RAF effector proteins at the plasma membrane. Recent advances in X-ray crystallography, cryo-EM, and super-resolution fluorescence microscopy offer structural and spatial insights, yet the molecular mechanisms involving protein-protein and protein-lipid interactions in RAS-mediated signaling require further characterization. This study utilizes single-molecule experimental techniques, nuclear magnetic resonance spectroscopy, and the computational Machine-Learned Modeling Infrastructure (MuMMI) to examine KRAS4b and RAF1 on a biologically relevant lipid bilayer. MuMMI captures long-timescale events while preserving detailed atomic descriptions, providing testable models for experimental validation. Both in vitro and computational studies reveal that RBDCRD binding alters KRAS lateral diffusion on the lipid bilayer, increasing cluster size and decreasing diffusion. RAS and membrane binding cause hydrophobic residues in the CRD region to penetrate the bilayer, stabilizing complexes through ß-strand elongation. These cooperative interactions among lipids, KRAS4b, and RAF1 are proposed as essential for forming nanoclusters, potentially a critical step in MAP kinase signal activation.

Message passing interface and multithreading hybrid for parallel molecular docking of large databases on petascale high performance computing machines.

A mixed parallel scheme that combines message passing interface (MPI) and multithreading was implemented in the AutoDock Vina molecular docking program. The resulting program, named VinaLC, was tested on the petascale high performance computing (HPC) machines at Lawrence Livermore National Laboratory. To exploit the typical cluster-type supercomputers, thousands of docking calculations were dispatched by the master process to run simultaneously on thousands of slave processes, where each docking calculation takes one slave process on one node, and within the node each docking calculation runs via multithreading on multiple CPU cores and shared memory. Input and output of the program and the data handling within the program were carefully designed to deal with large databases and ultimately achieve HPC on a large number of CPU cores. Parallel performance analysis of the VinaLC program shows that the code scales up to more than 15K CPUs with a very low overhead cost of 3.94%. One million flexible compound docking calculations took only 1.4 h to finish on about 15K CPUs. The docking accuracy of VinaLC has been validated against the DUD data set by the re-docking of X-ray ligands and an enrichment study, 64.4% of the top scoring poses have RMSD values under 2.0 Å. The program has been demonstrated to have good enrichment performance on 70% of the targets in the DUD data set. An analysis of the enrichment factors calculated at various percentages of the screening database indicates VinaLC has very good early recovery of actives.