Interdependence across time and length scales is common in biology, where atomic interactions can impact larger-scale phenomenon. Such dependence is especially true for a well-known cancer signaling pathway, where the membrane-bound RAS protein binds an effector protein called RAF. To capture the driving forces that bring RAS and RAF (represented as two domains, RBD and CRD) together on the plasma membrane, simulations with the ability to calculate atomic detail while having long time and large length- scales are needed. The Multiscale Machine-Learned Modeling Infrastructure (MuMMI) is able to resolve RAS/RAF protein-membrane interactions that identify specific lipid-protein fingerprints that enhance protein orientations viable for effector binding. MuMMI is a fully automated, ensemble-based multiscale approach connecting three resolution scales: (1) the coarsest scale is a continuum model able to simulate milliseconds of time for a 1 μm2 membrane, (2) the middle scale is a coarse-grained (CG) Martini bead model to explore protein-lipid interactions, and (3) the finest scale is an all-atom (AA) model capturing specific interactions between lipids and proteins. MuMMI dynamically couples adjacent scales in a pairwise manner using machine learning (ML). The dynamic coupling allows for better sampling of the refined scale from the adjacent coarse scale (forward) and on-the-fly feedback to improve the fidelity of the coarser scale from the adjacent refined scale (backward). MuMMI operates efficiently at any scale, from a few compute nodes to the largest supercomputers in the world, and is generalizable to simulate different systems. As computing resources continue to increase and multiscale methods continue to advance, fully automated multiscale simulations (like MuMMI) will be commonly used to address complex science questions.

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.

A rapid response is necessary to contain emergent biological outbreaks before they can become pandemics. The novel coronavirus (SARS-CoV-2) that causes COVID-19 was first reported in December of 2019 in Wuhan, China and reached most corners of the globe in less than two months. In just over a year since the initial infections, COVID-19 infected almost 100 million people worldwide. Although similar to SARS-CoV and MERS-CoV, SARS-CoV-2 has resisted treatments that are effective against other coronaviruses. Crystal structures of two SARS-CoV-2 proteins, spike protein and main protease, have been reported and can serve as targets for studies in neutralizing this threat. We have employed molecular docking, molecular dynamics simulations, and machine learning to identify from a library of 26 million molecules possible candidate compounds that may attenuate or neutralize the effects of this virus. The viability of selected candidate compounds against SARS-CoV-2 was determined experimentally by biolayer interferometry and FRET-based activity protein assays along with virus-based assays. In the pseudovirus assay, imatinib and lapatinib had IC50 values below 10 μM, while candesartan cilexetil had an IC50 value of approximately 67 µM against Mpro in a FRET-based activity assay. Comparatively, candesartan cilexetil had the highest selectivity index of all compounds tested as its half-maximal cytotoxicity concentration 50 (CC50) value was the only one greater than the limit of the assay (>100 μM).

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/.

Structure-based Deep Fusion models were recently shown to outperform several physics- and machine learning-based protein-ligand binding affinity prediction methods. As part of a multi-institutional COVID-19 pandemic response, over 500 million small molecules were computationally screened against four protein structures from the novel coronavirus (SARS-CoV-2), which causes COVID-19. Three enhancements to Deep Fusion were made in order to evaluate more than 5 billion docked poses on SARS-CoV-2 protein targets. First, the Deep Fusion concept was refined by formulating the architecture as one, coherently backpropagated model (Coherent Fusion) to improve binding-affinity prediction accuracy. Secondly, the model was trained using a distributed, genetic hyper-parameter optimization. Finally, a scalable, high-throughput screening capability was developed to maximize the number of ligands evaluated and expedite the path to experimental evaluation. In this work, we present both the methods developed for machine learning-based high-throughput screening and results from using our computational pipeline to find SARS-CoV-2 inhibitors.

BACKGROUND:Living organisms need to allocate their limited resources in a manner that optimizes their overall fitness by simultaneously achieving several different biological objectives. Examination of these biological trade-offs can provide invaluable information regarding the biophysical and biochemical bases behind observed cellular phenotypes. A quantitative knowledge of a cell system's critical objectives is also needed for engineering of cellular metabolism, where there is interest in mitigating the fitness costs that may result from human manipulation. RESULTS:To study metabolism in photoheterotrophs, we developed and validated a genome-scale model of metabolism in Rhodopseudomonas palustris, a metabolically versatile gram-negative purple non-sulfur bacterium capable of growing phototrophically on various carbon sources, including inorganic carbon and aromatic compounds. To quantitatively assess trade-offs among a set of important biological objectives during different metabolic growth modes, we used our new model to conduct an 8-dimensional multi-objective flux analysis of metabolism in R. palustris. Our results revealed that phototrophic metabolism in R. palustris is light-limited under anaerobic conditions, regardless of the available carbon source. Under photoheterotrophic conditions, R. palustris prioritizes the optimization of carbon efficiency, followed by ATP production and biomass production rate, in a Pareto-optimal manner. To achieve maximum carbon fixation, cells appear to divert limited energy resources away from growth and toward CO2 fixation, even in the presence of excess reduced carbon. We also found that to achieve the theoretical maximum rate of biomass production, anaerobic metabolism requires import of additional compounds (such as protons) to serve as electron acceptors. Finally, we found that production of hydrogen gas, of potential interest as a candidate biofuel, lowers the cellular growth rates under all circumstances. CONCLUSIONS:Photoheterotrophic metabolism of R. palustris is primarily regulated by the amount of light it can absorb and not the availability of carbon. However, despite carbon's secondary role as a regulating factor, R. palustris' metabolism strives for maximum carbon efficiency, even when this increased efficiency leads to slightly lower growth rates.

Patient-specific models of the ventricular myocardium, combined with the computational power to run rapid simulations, are approaching the level where they could be used for personalized cardiovascular medicine. A major remaining challenge is determining model parameters from available patient data, especially for models of the Purkinje-myocardial junctions (PMJs): the sites of initial ventricular electrical activation. There are no non-invasive methods for localizing PMJs in patients, and the relationship between the standard clinical ECG and PMJ model parameters is underexplored. Thus, this study aimed to determine the sensitivity of the QRS complex of the ECG to the anatomical location and regional number of PMJs. The QRS complex was simulated using an image-based human torso and biventricular model, and cardiac electrophysiology was simulated using Cardioid. The PMJs were modeled as discrete current injection stimuli, and the location and number of stimuli were varied within initial activation regions based on published experiments. Results indicate that the QRS complex features were most sensitive to the presence or absence of four "seed" stimuli, and adjusting locations of nearby "regional" stimuli provided finer tuning. Decreasing number of regional stimuli by an order of magnitude resulted in virtually no change in the QRS complex. Thus, a minimal 12-stimuli configuration was identified that resulted in physiological excitation, defined by QRS complex feature metrics and ventricular excitation pattern. Overall, the sensitivity results suggest that parameterizing PMJ location, rather than number, be given significantly higher priority in future studies creating personalized ventricular models from patient-derived ECGs.

Late-stage or post-market identification of adverse drug reactions (ADRs) is a significant public health issue and a source of major economic liability for drug development. Thus, reliable in silico screening of drug candidates for possible ADRs would be advantageous. In this work, we introduce a computational approach that predicts ADRs by combining the results of molecular docking and leverages known ADR information from DrugBank and SIDER. We employed a recently parallelized version of AutoDock Vina (VinaLC) to dock 906 small molecule drugs to a virtual panel of 409 DrugBank protein targets. L1-regularized logistic regression models were trained on the resulting docking scores of a 560 compound subset from the initial 906 compounds to predict 85 side effects, grouped into 10 ADR phenotype groups. Only 21% (87 out of 409) of the drug-protein binding features involve known targets of the drug subset, providing a significant probe of off-target effects. As a control, associations of this drug subset with the 555 annotated targets of these compounds, as reported in DrugBank, were used as features to train a separate group of models. The Vina off-target models and the DrugBank on-target models yielded comparable median area-under-the-receiver-operating-characteristic-curves (AUCs) during 10-fold cross-validation (0.60-0.69 and 0.61-0.74, respectively). Evidence was found in the PubMed literature to support several putative ADR-protein associations identified by our analysis. Among them, several associations between neoplasm-related ADRs and known tumor suppressor and tumor invasiveness marker proteins were found. A dual role for interstitial collagenase in both neoplasms and aneurysm formation was also identified. These associations all involve off-target proteins and could not have been found using available drug/on-target interaction data. This study illustrates a path forward to comprehensive ADR virtual screening that can potentially scale with increasing number of CPUs to tens of thousands of protein targets and millions of potential drug candidates.

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.

An approach to catalyst design is presented in which local potential energy surface models are first built to elucidate design principles and then used to identify larger scaffold motifs that match the target geometries. Carbon sequestrationviahydration is used as the model reaction, and three- and four-coordinatesp2orsp3nitrogen-ligand motifs are considered for ZnIImetals. The comparison of binding, activation and product release energies over a large range of interaction distances and angles suggests that four-coordinate short ZnII—Nsp3bond distances favor a rapid turnover for CO2hydration. This design strategy is then confirmed by computationally characterizing the reactivity of a known mimic over a range of metal–nitrogen bond lengths. A search of existing catalysts in a chemical database reveals structures that match the target geometry from model calculations, and subsequent calculations have identified these structures as potentially effective for CO2hydration and sequestration.

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.

Product regioselectivity as influenced by molecular recognition is a key aspect of enzyme catalysis. We applied large-scale two-dimensional (2D) umbrella sampling (USP) simulations to characterize acetaminophen (APAP) binding in the active sites of the family of Cytochrome P450 (CYP) enzymes as a case study to show the different regioselectivity exhibited by a single substrate in comparative enzymes. Our results successfully explain the experimentally observed product regioselectivity for all five human CYPs included in this study, demonstrating that binding events play an important role in determining regioselectivity. In CYP2C9 and CYP3A4, weak interactions in an overall large active site cavity result in a fairly small binding free energy difference between APAP reactive binding states, consistent with experimental results that show little preference for resulting metabolites. In contrast, in CYP1A2 and CYP2E1, APAP is strongly restrained by a compact binding pocket, leading to a preferred binding conformation. The calculated binding equilibrium of APAP within the compact active site of CYP2A6 is able to predict the experimentally documented product ratios and is also applied to explain APAP regioselectivity in CYP1A2 and CYP2C9. APAP regioselectivity seems to be related to the selectivity for one binding conformation over another binding conformation as dictated by the size and shape of the active site. Additionally, unlike docking and molecular dynamics (MD), our free energy calculations successfully reproduced a unique APAP pose in CYP3A4 that had been reported experimentally, suggesting this approach is well suited to find the realistic binding pose and the lowest-energy starting structure for studying the chemical reaction step in the future.

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.