Biosynthesis of eriodictyol in citrus waster by endowing P450BM3 activity of naringenin hydroxylation

The flavonoid naringenin is abundantly present in pomelo peels, and the unprocessed naringenin in wastes is not friendly for the environment once discarded directly. Fortunately, the hydroxylated product of eriodictyol from naringenin exhibits remarkable antioxidant and anticancer properties. The P450s was suggested promising for the bioconversion of the flavonoids, but less naturally existed P450s show hydroxylation activity to C3' of the naringenin. By well analyzing the catalytic mechanism and the conformations of the naringenin in P450, we proposed that the intermediate Cmpd I ((porphyrin)Fe = O) is more reasonable as key conformation for the hydrolyzation, and the distance between C3'/C5' of naringenin to the O atom of CmpdI determines the hydroxylating activity for the naringenin. Thus, the "flying kite model" that gradually drags the C-H bond of the substrate to the O atom of CmpdI was put forward for rational design. With ab initio design, we successfully endowed the self-sufficient P450-BM3 hydroxylic activity to naringenin and obtained mutant M5-5, with kcat, Km, and kcat/Km values of 230.45 min-1, 310.48 µM, and 0.742 min-1 µM-1, respectively. Furthermore, the mutant M4186 was screened with kcat/Km of 4.28-fold highly improved than the reported M13. The M4186 also exhibited 62.57% yield of eriodictyol, more suitable for the industrial application. This study provided a theoretical guide for the rational design of P450s to the nonnative compounds. KEY POINTS: •The compound I is proposed as the starting point for the rational design of the P450BM3 •"Flying kite model" is proposed based on the distance between O of Cmpd I and C3'/C5' of naringenin •Mutant M15-5 with 1.6-fold of activity than M13 was obtained by ab initio modification.

A deep learning based multi-model approach for predicting drug-like chemical compound's toxicity

Ensuring the safety and efficacy of chemical compounds is crucial in small-molecule drug development. In the later stages of drug development, toxic compounds pose a significant challenge, losing valuable resources and time. Early and accurate prediction of compound toxicity using deep learning models offers a promising solution to mitigate these risks during drug discovery. In this study, we present the development of several deep-learning models aimed at evaluating different types of compound toxicity, including acute toxicity, carcinogenicity, hERG_cardiotoxicity (the human ether-a-go-go related gene caused cardiotoxicity), hepatotoxicity, and mutagenicity. To address the inherent variations in data size, label type, and distribution across different types of toxicity, we employed diverse training strategies. Our first approach involved utilizing a graph convolutional network (GCN) regression model to predict acute toxicity, which achieved notable performance with Pearson R 0.76, 0.74, and 0.65 for intraperitoneal, intravenous, and oral administration routes, respectively. Furthermore, we trained multiple GCN binary classification models, each tailored to a specific type of toxicity. These models exhibited high area under the curve (AUC) scores, with an impressive AUC of 0.69, 0.77, 0.88, and 0.79 for predicting carcinogenicity, hERG_cardiotoxicity, mutagenicity, and hepatotoxicity, respectively. Additionally, we have used the approved drug dataset to determine the appropriate threshold value for the prediction score in model usage. We integrated these models into a virtual screening pipeline to assess their effectiveness in identifying potential low-toxicity drug candidates. Our findings indicate that this deep learning approach has the potential to significantly reduce the cost and risk associated with drug development by expediting the selection of compounds with low toxicity profiles. Therefore, the models developed in this study hold promise as critical tools for early drug candidate screening and selection.

Deep Learning-Based construction of a Drug-Like compound database and its application in virtual screening of HsDHODH inhibitors

The process of virtual screening relies heavily on the databases, but it is disadvantageous to conduct virtual screening based on commercial databases with patent-protected compounds, high compound toxicity and side effects. Therefore, this paper utilizes generative recurrent neural networks (RNN) containing long short-term memory (LSTM) cells to learn the properties of drug compounds in the DrugBank, aiming to obtain a new and virtual screening compounds database with drug-like properties. Ultimately, a compounds database consisting of 26,316 compounds is obtained by this method. To evaluate the potential of this compounds database, a series of tests are performed, including chemical space, ADME properties, compound fragmentation, and synthesizability analysis. As a result, it is proved that the database is equipped with good drug-like properties and a relatively new backbone, its potential in virtual screening is further tested. Finally, a series of seedling compounds with completely new backbones are obtained through docking and binding free energy calculations.

Substituent effects on the selectivity of ambimodal [6+4]/[4+2] cycloaddition

In this work, we report a density functional theory (DFT) study and a dynamical trajectory study of substituent effects on the ambimodal [6+4]/[4+2] cycloaddition proposed for 1,3,5,10,12-cycloheptadecapentaene, referred to as cycloheptadecapentaene. The proposed cycloaddition proceeds through an ambimodal transition state, which results in both a [6+4] adduct a [4+2] adduct directly. The [6+4] adduct can be readily converted to the [4+2] adduct via a Cope rearrangement. We study the selectivity of the reaction with regard to the position of substituents, steric effects of substituents, and electronic effects of substituents. In the dynamical trajectory study, we find that nitro-substituted reactants lead to a new product from the ambimodal transition state via the hetero Diels-Alder reaction, and this new product can then be converted to a [4+2] adduct by a hetero [3, 3]-sigmatropic rearrangement. These results may provide insights for designing more bridged heterocyclic compounds.

An Efficient Approach to the Accurate Prediction of Mutational Effects in Antigen Binding to the MHC1

The major histocompatibility complex (MHC) can recognize and bind to external peptides to generate effective immune responses by presenting the peptides to T cells. Therefore, understanding the binding modes of peptide-MHC complexes (pMHC) and predicting the binding affinity of pMHCs play a crucial role in the rational design of peptide vaccines. In this study, we employed molecular dynamics (MD) simulations and free energy calculations with an Alanine Scanning with Generalized Born and Interaction Entropy (ASGBIE) method to investigate the protein-peptide interaction between HLA-A*02:01 and the G9209 peptide derived from the melanoma antigen gp100. The energy contribution of individual residue was calculated using alanine scanning, and hotspots on both the MHC and the peptides were identified. Our study shows that the pMHC binding is dominated by the van der Waals interactions. Furthermore, we optimized the ASGBIE method, achieving a Pearson correlation coefficient of 0.91 between predicted and experimental binding affinity for mutated antigens. This represents a significant improvement over the conventional MM/GBSA method, which yields a Pearson correlation coefficient of 0.22. The computational protocol developed in this study can be applied to the computational screening of antigens for the MHC1 as well as other protein-peptide binding systems.

A precise swaying map for how promiscuous cellobiose-2-epimerase operate bi-reaction

Promiscuous enzymes play a crucial role in organism survival and new reaction mining. However, comprehensive mapping of the catalytic and regulatory mechanisms hasn't been well studied due to the characteristic complexity. The cellobiose 2-epimerase from Caldicellulosiruptor saccharolyticus (CsCE) with complex epimerization and isomerization was chosen to comprehensively investigate the promiscuous mechanisms. Here, the catalytic frame of ring-opening, cis-enediol mediated catalysis and ring-closing was firstly determined. To map the full view of promiscuous CE, the structure of CsCE complex with the isomerized product glucopyranosyl-β1,4-fructose was determined. Combined with computational calculation, the promiscuity was proved a precise cooperation of the double subsites, loop rearrangement, and intermediate swaying. The flexible loop was like a gear, whose structural reshaping regulates the sway of the intermediates between the two subsites of H377-H188 and H377-H247, and thus regulates the catalytic directions. The different protonated states of cis-enediol intermediate catalyzed by H188 were the key point for the catalysis. The promiscuous enzyme tends to utilize all elements at hand to carry out the promiscuous functions.

Screening Power of End-Point Free-Energy Calculations in Cucurbituril Host-Guest Systems

End-point free-energy methods as an indispensable component in virtual screening are commonly recognized as a tool with a certain level of screening power in pharmaceutical research. While a huge number of records could be found for end-point applications in protein-ligand, protein-protein, and protein-DNA complexes from academic and industrial reports, up to now, there is no large-scale benchmark in host-guest complexes supporting the screening power of end-point free-energy techniques. A good benchmark requires a data set of sufficient coverage of pharmaceutically relevant chemical space, a long-time sampling length supporting the trajectory approximation of the ensemble average, and a sufficient sample size of receptor-acceptor pairs to stabilize the performance statistics. In this work, selecting a popular family of macrocyclic hosts named cucurbiturils, we construct a large data set containing 154 host-guest pairs, perform extensive end-point sampling of several hundred nanosecond lengths for each system, and extract the free-energy estimates with a variety of end-point free-energy techniques, including the advanced three-trajectory dielectric-constant-variable regime proposed in our recent work. The best-performing end-point protocol employs GAFF2 for solute descriptions, the three-trajectory end-point sampling regime, and the MM/GBSA Hamiltonian in free-energy extraction, achieving a high ranking metrics of Kendall τ > 0.6, a Pearlman predictive index of ∼0.8, and a high scoring power of Pearson r > 0.8. The current project as the first large-scale systematic benchmark of end-point methods in host-guest complexes in academic publications provides solid evidence of the applicability of end-point techniques and direct guidance of computational setups in practical host-guest systems.

Entropy driven cooperativity effect in multi-site drug optimization targeting SARS-CoV-2 papain-like protease

Papain-like protease (PLpro), a non-structural protein encoded by SARS-CoV-2, is an important therapeutic target. Regions 1 and 5 of an existing drug, GRL0617, can be optimized to produce cooperativity with PLpro binding, resulting in stronger binding affinity. This work investigated the origin of the cooperativity using molecular dynamics simulations combined with the interaction entropy (IE) method. The regions' improvement exhibits cooperativity by calculating the binding free energies between the complex of PLpro-inhibitor. The thermodynamic integration method further verified the cooperativity generated in the drug improvement. To further determine the specific source of cooperativity, enthalpy and entropy in the complexes were calculated using molecular mechanics/generalized Born surface area and IE. The results show that the entropic change is an important contributor to the cooperativity. Our study also identified residues P248, Q269, and T301 that play a significant role in cooperativity. The optimization of the inhibitor stabilizes these residues and minimizes the entropic loss, and the cooperativity observed in the binding free energy can be attributed to the change in the entropic contribution of these residues. Based on our research, the application of cooperativity can facilitate drug optimization, and provide theoretical ideas for drug development that leverage cooperativity by reducing the contribution of entropy through multi-locus binding.

Assessment of an Electrostatic Energy-Based Charge Model for Modeling the Electrostatic Interactions in Water Solvent

The protein force field based on the restrained electrostatic potential (RESP) charges has limitations in accurately describing hydrogen bonding interactions in proteins. To address this issue, we propose an alternative approach called the electrostatic energy-based charges (EEC) model, which shows improved performance in describing electrostatic interactions (EIs) of hydrogen bonds in proteins. In this study, we further investigate the performance of the EEC model in modeling EIs in water solvent. Our findings demonstrate that the fixed EEC model can effectively reproduce the quantum mechanics/molecular mechanics (QM/MM)-calculated EIs between a water molecule and various water solvent environments. However, to achieve the same level of computational accuracy, the electrostatic potential (ESP) charge model needs to fluctuate according to the electrostatic environment. Our analysis indicates that the requirement for charge adjustments depends on the specific mathematical and physical representation of EIs as a function of the environment for deriving charges. By comparing with widely used empirical water models calibrated to reproduce experimental properties, we confirm that the performance of the EEC model in reproducing QM/MM EIs is similar to that of general purpose TIP4P-like water models such as TIP4P-Ew and TIP4P/2005. When comparing the computed 10,000 distinct EI values within the range of -40 to 0 kcal/mol with the QM/MM results calculated at the MP2/aug-cc-pVQZ/TIP3P level, we noticed that the mean unsigned error (MUE) for the EEC model is merely 0.487 kcal/mol, which is remarkably similar to the MUE values of the TIP4P-Ew (0.63 kcal/mol) and TIP4P/2005 (0.579 kcal/mol) models. However, both the RESP method and the TIP3P model exhibit a tendency to overestimate the EIs, as evidenced by their higher MUE values of 1.761 and 1.293 kcal/mol, respectively. EEC-based molecular dynamics simulations have demonstrated that, when combined with appropriate van der Waals parameters, the EEC model can closely reproduce oxygen-oxygen radial distribution function and density of water, showing a remarkable similarity to the well-established TIP4P-like empirical water models. Our results demonstrate that the EEC model has the potential to build force fields with comparable accuracy to more sophisticated empirical TIP4P-like water models.

Combined Antibodies Evusheld against the SARS-CoV-2 Omicron Variants BA.1.1 and BA.5: Immune Escape Mechanism from Molecular Simulation

The Omicron lineage of SARS-CoV-2, which was first reported in November 2021, has spread globally and become dominant, splitting into several sublineages. Experiments have shown that Omicron lineage has escaped or reduced the activity of existing monoclonal antibodies, but the origin of escape mechanism caused by mutation is still unknown. This work uses molecular dynamics and umbrella sampling methods to reveal the escape mechanism of BA.1.1 to monoclonal antibody (mAb) Tixagevimab (AZD1061) and BA.5 to mAb Cilgavimab (AZD8895), both mAbs were combined to form antibody cocktail, Evusheld (AZD7442). The binding free energy of BA.1.1-AZD1061 and BA.5-AZD8895 has been severely reduced due to multiple-site mutated Omicron variants. Our results show that the two Omicron variants, which introduce a substantial number of positively charged residues, can weaken the electrostatic attraction between the receptor binding domain (RBD) and AZD7442, thus leading to a decrease in affinity. Additionally, using umbrella sampling along dissociation pathway, we found that the two Omicron variants severely impaired the interaction between the RBD of SARS-CoV-2's spike glycoprotein (S protein) and complementary determining regions (CDRs) of mAbs, especially in CDR3H. Although mAbs AZD8895 and AZD1061 are knocked out by BA.5 and BA.1.1, respectively, our results confirm that the antibody cocktail AZD7442 retains activity against BA.1.1 and BA.5 because another antibody is still on guard. The study provides theoretical insights for mAbs interacting with BA.1.1 and BA.5 from both energetic and dynamic perspectives, and we hope this will help in developing new monoclonals and combinations to protect those unable to mount adequate vaccine responses.

DeepBindGCN: Integrating Molecular Vector Representation with Graph Convolutional Neural Networks for Protein-Ligand Interaction Prediction

The core of large-scale drug virtual screening is to select the binders accurately and efficiently with high affinity from large libraries of small molecules in which non-binders are usually dominant. The binding affinity is significantly influenced by the protein pocket, ligand spatial information, and residue types/atom types. Here, we used the pocket residues or ligand atoms as the nodes and constructed edges with the neighboring information to comprehensively represent the protein pocket or ligand information. Moreover, the model with pre-trained molecular vectors performed better than the one-hot representation. The main advantage of DeepBindGCN is that it is independent of docking conformation, and concisely keeps the spatial information and physical-chemical features. Using TIPE3 and PD-L1 dimer as proof-of-concept examples, we proposed a screening pipeline integrating DeepBindGCN and other methods to identify strong-binding-affinity compounds. It is the first time a non-complex-dependent model has achieved a root mean square error (RMSE) value of 1.4190 and Pearson r value of 0.7584 in the PDBbind v.2016 core set, respectively, thereby showing a comparable prediction power with the state-of-the-art affinity prediction models that rely upon the 3D complex. DeepBindGCN provides a powerful tool to predict the protein-ligand interaction and can be used in many important large-scale virtual screening application scenarios.

Revealing the Sequence Characteristics and Molecular Mechanisms of ACE Inhibitory Peptides by Comprehensive Characterization of 160,000 Tetrapeptides

Chronic diseases, such as hypertension, cause great harm to human health. Conventional drugs have promising therapeutic effects, but also cause significant side effects. Food-sourced angiotensin-converting enzyme (ACE) inhibitory peptides are an excellent therapeutic alternative to pharmaceuticals, as they have fewer side effects. However, there is no systematic and effective screening method for ACE inhibitory peptides, and the lack of understanding of the sequence characteristics and molecular mechanism of these inhibitory peptides poses a major obstacle to the development of ACE inhibitory peptides. Through systematically calculating the binding effects of 160,000 tetrapeptides with ACE by molecular docking, we found that peptides with Tyr, Phe, His, Arg, and especially Trp were the characteristic amino acids of ACE inhibitory peptides. The tetrapeptides of WWNW, WRQF, WFRV, YYWK, WWDW, and WWTY rank in the top 10 peptides exhibiting significantly high ACE inhibiting behaviors, with IC₅₀ values between 19.98 ± 8.19 μM and 36.76 ± 1.32 μM. Salt bridges, π-π stacking, π-cations, and hydrogen bonds contributed to the high binding characteristics of the inhibitors and ACE. Introducing eight Trp into rabbit skeletal muscle protein (no Trp in wide sequence) endowed the protein with a more than 90% ACE inhibition rate, further suggesting that meat with a high content of Trp could have potential utility in hypertension regulation. This study provides a clear direction for the development and screening of ACE inhibitory peptides.

Binding of berberine derivates to G-quadruplex: insight from a computational study

Human telomerase exhibits significant activity in cancer cells relative to normal cells, which contributes to the immortal proliferation of cancer cells. To counter this, the stabilization of G-quadruplexes formed in the guanine-rich sequence of the cancer cell chromosome has emerged as a promising avenue for anti-cancer therapy. Berberine (BER), an alkaloid that is derived from traditional Chinese medicines, has shown potential for stabilizing G-quadruplexes. To investigate the atomic interactions between G-quadruplexes and BER and its derivatives, molecular dynamics simulations were conducted. Modeling the interactions between G-quadruplexes and ligands accurately is challenging due to the strong negative charge of nucleic acids. Thus, various force fields and charge models for the G-quadruplex and ligands were tested to obtain precise simulation results. The binding energies were calculated by a combination of molecular mechanics/generalized Born surface area and interaction entropy methods, and the calculated results correlated well with experimental results. B-factor and hydrogen bond analyses demonstrated that the G-quadruplex was more stable in the presence of ligands than in the absence of ligands. Calculation of the binding free energy showed that the BER derivatives bind to a G-quadruplex with higher affinity than that of BER. The breakdown of the binding free energy to per-nucleotide energies suggested that the first G-tetrad played a primary role in binding. Additionally, energy and geometric properties analyses indicated that van der Waals interactions were the most favorable interactions between the derivatives and the G-quadruplexes. Overall, these findings provide crucial atomic-level insights into the binding of G-quadruplexes and their inhibitors.

MolTaut: A Tool for the Rapid Generation of Favorable Tautomer in Aqueous Solution

Fast and proper treatment of the tautomeric states for drug-like molecules is critical in computer-aided drug discovery since the major tautomer of a molecule determines its pharmacophore features and physical properties. We present MolTaut, a tool for the rapid generation of favorable states of drug-like molecules in water. MolTaut works by enumerating possible tautomeric states with tautomeric transformation rules, ranking tautomers with their relative internal energies and solvation energies calculated by AI-based models, and generating preferred ionization states according to predicted microscopic pK_a. Our test shows that the ranking ability of the AI-based tautomer scoring approach is comparable to the DFT method (wB97X/6-31G*//M062X/6-31G*/SMD) from which the AI models try to learn. We find that the substitution effect on tautomeric equilibrium is well predicted by MolTaut, which is helpful in computer-aided ligand design. The source code of MolTaut is freely available to researchers and can be accessed at https://github.com/xundrug/moltaut. To facilitate the usage of MolTaut by medicinal chemists, we made a free web server, which is available at http://moltaut.xundrug.cn. MolTaut is a handy tool for investigating the tautomerization issue in drug discovery.

Comprehensive Evaluation of End-Point Free Energy Techniques in Carboxylated-Pillar[6]arene Host–Guest Binding: III. Force-Field Comparison, Three-Trajectory Realization and Further Dielectric Augmentation

Host–guest binding, despite the relatively simple structural and chemical features of individual components, still poses a challenge in computational modelling. The extreme underperformance of standard end-point methods in host–guest binding makes them practically useless. In the current work, we explore a potentially promising modification of the three-trajectory realization. The alteration couples the binding-induced structural reorganization into free energy estimation and suffers from dramatic fluctuations in internal energies in protein–ligand situations. Fortunately, the relatively small size of host–guest systems minimizes the magnitude of internal fluctuations and makes the three-trajectory realization practically suitable. Due to the incorporation of intra-molecular interactions in free energy estimation, a strong dependence on the force field parameters could be incurred. Thus, a term-specific investigation of transferable GAFF derivatives is presented, and noticeable differences in many aspects are identified between commonly applied GAFF and GAFF2. These force-field differences lead to different dynamic behaviors of the macrocyclic host, which ultimately would influence the end-point sampling and binding thermodynamics. Therefore, the three-trajectory end-point free energy calculations are performed with both GAFF versions. Additionally, due to the noticeable differences between host dynamics under GAFF and GAFF2, we add additional benchmarks of the single-trajectory end-point calculations. When only the ranks of binding affinities are pursued, the three-trajectory realization performs very well, comparable to and even better than the regressed PBSA_E scoring function and the dielectric constant-variable regime. With the GAFF parameter set, the TIP3P water in explicit solvent sampling and either PB or GB implicit solvent model in free energy estimation, the predictive power of the three-trajectory realization in ranking calculations surpasses all existing end-point methods on this dataset. We further combine the three-trajectory realization with another promising modified end-point regime of varying the interior dielectric constant. The combined regime does not incur sizable improvements for ranks and deviations from experiment exhibit non-monotonic variations.

Deep Learning-Based Bioactive Therapeutic Peptide Generation and Screening

Many bioactive peptides demonstrated therapeutic effects over complicated diseases, such as antiviral, antibacterial, anticancer, etc. It is possible to generate a large number of potentially bioactive peptides using deep learning in a manner analogous to the generation of de novo chemical compounds using the acquired bioactive peptides as a training set. Such generative techniques would be significant for drug development since peptides are much easier and cheaper to synthesize than compounds. Despite the limited availability of deep learning-based peptide-generating models, we have built an LSTM model (called LSTM_Pep) to generate de novo peptides and fine-tuned the model to generate de novo peptides with specific prospective therapeutic benefits. Remarkably, the Antimicrobial Peptide Database has been effectively utilized to generate various kinds of potential active de novo peptides. We proposed a pipeline for screening those generated peptides for a given target and used the main protease of SARS-COV-2 as a proof-of-concept. Moreover, we have developed a deep learning-based protein-peptide prediction model (DeepPep) for rapid screening of the generated peptides for the given targets. Together with the generating model, we have demonstrated that iteratively fine-tuning training, generating, and screening peptides for higher-predicted binding affinity peptides can be achieved. Our work sheds light on developing deep learning-based methods and pipelines to effectively generate and obtain bioactive peptides with a specific therapeutic effect and showcases how artificial intelligence can help discover de novo bioactive peptides that can bind to a particular target.

Rational Design of a Highly Specific Prolyl Endopeptidase To Activate the Antihypertensive Effect of Peptides

Enzymatic hydrolysis of food-derived proteins to produce bioactive peptides could activate food functions such as antihypertension. However, the diversity of enzymatic hydrolysis products can reduce bioactive peptides' efficacy. Highly specific proteases can homogenize the hydrolysis products to reduce the production of impotent peptides. In this study, we successfully obtained M. xanthus prolyl endopeptidase mutant Y451M by constraint/free molecular dynamics simulations and binding energy calculations. The specificity of Y451M for proline was increased by 286 % compared to WT, while its activity was almost unchanged. Milk-derived substrates processed with Y451M showed an antihypertensive effect that was 567 % higher than without enzymes. The ability to activate food antihypertension increased 152 % and the use of enzyme by 192 % compared with WT. Specific proteases are thus valuable tools in the processing of complex substrates to obtain bioactive peptides.

Discovery of novel inhibitors of SARS-CoV-2 main protease

Corona Virus Disease 2019 (COVID-19), referred to as 'New Coronary Pneumonia', is a type of acute infectious disease caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) infection. M^pro is one of the main targets for treating COVID-19. The current research on M^pro mainly focuses on the repurposing of old drugs, and there are only a few novel ligands that inhibit M^pro. In this research, we used computational free energy calculation to screen a compound library against M^pro, and discovered four novel compounds with the two best compounds (AG-690/13507628 and AG-690/13507724) having experimental measured IC₅₀ of just under 3 μM and low cell toxicity. Detailed decomposition of the interactions between the inhibitors and M^pro reveals key interacting residues and interactions that determine the activity. The results from this study should provide a basis for further development of anti-SARS-CoV-2 drugs.Communicated by Ramaswamy H. Sarma.

Insights into small molecule inhibitor bindings to PD-L1 with residue-specific binding free energy calculation

Targeting the immunological checkpoint PD-1/PD-L1 with antibodies has shown opportunities to improve cancer treatment in recent years. However, antibody therapy is a double-edged sword with high cost, low patient tolerance, lack of oral bioavailability, and a reaction to most solid tumors that prevents the adoption of antibodies. Advancement of small-molecule PD-1/PD-L1 inhibitors that could overwhelm these drawbacks is sluggish because of the poor pharmacodynamic properties and shallow pocket of the PD-1/PD-L1 binding interface. Recently, a number of compounds have been discovered to bind the PD-L1/PD-L1 dimer interface, providing an excellent alternative to inhibit the interaction between PD-1/PD-L1 and small molecules. Quantitative characterization of PD-L1 interactions with these inhibitors will advance the design of novel and efficient inhibitors in the future. Here, the binding free energies of 35 PD-L1 dimer inhibitors have been calculated using the alanine-scanning-interaction-entropy (AS-IE) method. Hotspot residues on PD-L1 and potential modification groups on the inhibitors were identified. The experimental results for the AS-IE method were better correlated than the classical MM/GBSA method. These results may set the stage for the design the more powerful PD-L1 inhibitors.Communicated by Ramaswamy H. Sarma.

An efficient method to predict protein thermostability in alanine mutation

The relationship between protein sequence and its thermodynamic stability is a critical aspect of computational protein design. In this work, we present a new theoretical method to calculate the free energy change (ΔΔG) resulting from a single-point amino acid mutation to alanine in a protein sequence. The method is derived based on physical interactions and is very efficient in estimating the free energy changes caused by a series of alanine mutations from just a single molecular dynamics (MD) trajectory. Numerical calculations are carried out on a total of 547 alanine mutations in 19 diverse proteins whose experimental results are available. The comparison between the experimental ΔΔG_exp and the calculated values shows a generally good correlation with a correlation coefficient of 0.67. Both the advantages and limitations of this method are discussed. This method provides an efficient and valuable tool for protein design and engineering.

Comprehensive evaluation of end-point free energy techniques in carboxylated-pillar[6]arene host-guest binding: II. regression and dielectric constant

End-point free energy calculations as a powerful tool have been widely applied in protein-ligand and protein-protein interactions. It is often recognized that these end-point techniques serve as an option of intermediate accuracy and computational cost compared with more rigorous statistical mechanic models (e.g., alchemical transformation) and coarser molecular docking. However, it is observed that this intermediate level of accuracy does not hold in relatively simple and prototypical host-guest systems. Specifically, in our previous work investigating a set of carboxylated-pillar[6]arene host-guest complexes, end-point methods provide free energy estimates deviating significantly from the experimental reference, and the rank of binding affinities is also incorrectly computed. These observations suggest the unsuitability and inapplicability of standard end-point free energy techniques in host-guest systems, and alteration and development are required to make them practically usable. In this work, we consider two ways to improve the performance of end-point techniques. The first one is the PBSA_E regression that varies the weights of different free energy terms in the end-point calculation procedure, while the second one is considering the interior dielectric constant as an additional variable in the end-point equation. By detailed investigation of the calculation procedure and the simulation outcome, we prove that these two treatments (i.e., regression and dielectric constant) are manipulating the end-point equation in a somehow similar way, i.e., weakening the electrostatic contribution and strengthening the non-polar terms, although there are still many detailed differences between these two methods. With the trained end-point scheme, the RMSE of the computed affinities is improved from the standard ~ 12 kcal/mol to ~ 2.4 kcal/mol, which is comparable to another altered end-point method (ELIE) trained with system-specific data. By tuning PBSA_E weighting factors with the host-specific data, it is possible to further decrease the prediction error to ~ 2.1 kcal/mol. These observations along with the extremely efficient optimized-structure computation procedure suggest the regression (i.e., PBSA_E as well as its GBSA_E extension) as a practically applicable solution that brings end-point methods back into the library of usable tools for host-guest binding. However, the dielectric-constant-variable scheme cannot effectively minimize the experiment-calculation discrepancy for absolute binding affinities, but is able to improve the calculation of affinity ranks. This phenomenon is somehow different from the protein-ligand case and suggests the difference between host-guest and biomacromolecular (protein-ligand and protein-protein) systems. Therefore, the spectrum of tools usable for protein-ligand complexes could be unsuitable for host-guest binding, and numerical validations are necessary to screen out really workable solutions in these 'prototypical' situations.

Molecular Mechanism of the Non-Covalent Orally Targeted SARS-CoV-2 Mpro Inhibitor S-217622 and Computational Assessment of Its Effectiveness against Mainstream Variants

Convenient and efficient therapeutic agents are urgently needed to block the continued spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Here, the mechanism for the novel orally targeted SARS-CoV-2 main protease (M^pro) inhibitor S-217622 is revealed through a molecular dynamics simulation. The difference in the movement modes of the S-217622-M^pro complex and apo-M^pro suggested S-217622 could inhibit the motility intensity of M^pro, thus maintaining their stable binding. Subsequent energy calculations showed that the P2 pharmacophore possessed the highest energy contribution among the three pharmacophores of S-217622. Additionally, hot-spot residues H41, M165, C145, E166, and H163 have strong interactions with S-217622. To further investigate the resistance of S-217622 to six mainstream variants, the binding modes of S-217622 with these variants were elucidated. The subtle differences in energy compared to that of the wild type implied that the binding patterns of these systems were similar, and S-217622 still inhibited these variants. We hope this work will provide theoretical insights for optimizing novel targeted M^pro drugs.

Comprehensive evaluation of end-point free energy techniques in carboxylated-pillar[6]arene host-guest binding: I. Standard procedure

Despite the massive application of end-point free energy methods in protein-ligand and protein-protein interactions, computational understandings about their performance in relatively simple and prototypical host-guest systems are limited. In this work, we present a comprehensive benchmark calculation with standard end-point free energy techniques in a recent host-guest dataset containing 13 host-guest pairs involving the carboxylated-pillar[6]arene host. We first assess the charge schemes for solutes by comparing the charge-produced electrostatics with many ab initio references, in order to obtain a preliminary albeit detailed view of the charge quality. Then, we focus on four modelling details of end-point free energy calculations, including the docking procedure for the generation of initial condition, the charge scheme for host and guest molecules, the water model used in explicit-solvent sampling, and the end-point methods for free energy estimation. The binding thermodynamics obtained with different modelling schemes are compared with experimental references, and some practical guidelines on maximizing the performance of end-point methods in practical host-guest systems are summarized. Further, we compare our simulation outcome with predictions in the grand challenge and discuss further developments to improve the prediction quality of end-point free energy methods. Overall, unlike the widely acknowledged applicability in protein-ligand binding, the standard end-point calculations cannot produce useful outcomes in host-guest binding and thus are not recommended unless alterations are performed.

An Efficient Modern Strategy to Screen Drug Candidates Targeting RdRp of SARS-CoV-2 With Potentially High Selectivity and Specificity

Desired drug candidates should have both a high potential binding chance and high specificity. Recently, many drug screening strategies have been developed to screen compounds with high possible binding chances or high binding affinity. However, there is still no good solution to detect whether those selected compounds possess high specificity. Here, we developed a reverse DFCNN (Dense Fully Connected Neural Network) and a reverse docking protocol to check a given compound’s ability to bind diversified targets and estimate its specificity with homemade formulas. We used the RNA-dependent RNA polymerase (RdRp) target as a proof-of-concept example to identify drug candidates with high selectivity and high specificity. We first used a previously developed hybrid screening method to find drug candidates from an 8888-size compound database. The hybrid screening method takes advantage of the deep learning-based method, traditional molecular docking, molecular dynamics simulation, and binding free energy calculated by metadynamics, which should be powerful in selecting high binding affinity candidates. Also, we integrated the reverse DFCNN and reversed docking against a diversified 102 proteins to the pipeline for assessing the specificity of those selected candidates, and finally got compounds that have both predicted selectivity and specificity. Among the eight selected candidates, Platycodin D and Tubeimoside III were confirmed to effectively inhibit SARS-CoV-2 replication in vitro with EC₅₀ values of 619.5 and 265.5 nM, respectively. Our study discovered that Tubeimoside III could inhibit SARS-CoV-2 replication potently for the first time. Furthermore, the underlying mechanisms of Platycodin D and Tubeimoside III inhibiting SARS-CoV-2 are highly possible by blocking the RdRp cavity according to our screening procedure. In addition, the careful analysis predicted common critical residues involved in the binding with active inhibitors Platycodin D and Tubeimoside III, Azithromycin, and Pralatrexate, which hopefully promote the development of non-covalent binding inhibitors against RdRp.

Immune Escape Mechanisms of SARS-CoV-2 Delta and Omicron Variants against Two Monoclonal Antibodies That Received Emergency Use Authorization

Multiple-site mutated SARS-CoV-2 Delta and Omicron variants may trigger immune escape against existing monoclonal antibodies. Here, molecular dynamics simulations combined with the interaction entropy method reveal the escape mechanism of Delta/Omicron variants to Bamlanivimab/Etesevimab. The result shows the significantly reduced binding affinity of the Omicron variant for both antibodies, due to the introduction of positively charged residues that greatly weaken their electrostatic interactions. Meanwhile, significant structural deflection induces fewer atomic contacts and an unstable binding mode. As for the Delta variant, the reduced binding affinity for Bamlanivimab is owing to the alienation of the receptor-binding domain to the main part of this antibody, and the binding mode of the Delta variant to Etesevimab is similar to that of the wild type, suggesting that Etesevimab could still be effective against the Delta variant. We hope this work will provide timely theoretical insights into developing antibodies to prevalent and possible future variants of SARS-CoV-2.

Generating and screening de novo compounds against given targets using ultrafast deep learning models as core components

Deep learning is an artificial intelligence technique in which models express geometric transformations over multiple levels. This method has shown great promise in various fields, including drug development. The availability of public structure databases prompted the researchers to use generative artificial intelligence models to narrow down their search of the chemical space, a novel approach to chemogenomics and de novo drug development. In this study, we developed a strategy that combined an accelerated LSTM_Chem (long short-term memory for de novo compounds generation), dense fully convolutional neural network (DFCNN), and docking to generate a large number of de novo small molecular chemical compounds for given targets. To demonstrate its efficacy and applicability, six important targets that account for various human disorders were used as test examples. Moreover, using the M protease as a proof-of-concept example, we find that iteratively training with previously selected candidates can significantly increase the chance of obtaining novel compounds with higher and higher predicted binding affinities. In addition, we also check the potential benefit of obtaining reliable final de novo compounds with the help of MD simulation and metadynamics simulation. The generation of de novo compounds and the discovery of binders against various targets proposed here would be a practical and effective approach. Assessing the efficacy of these top de novo compounds with biochemical studies is promising to promote related drug development.

A fast-slow method to treat solute dynamics in explicit solvent

Aiming to reduce the computational cost in the current explicit solvent molecular dynamics (MD) simulation, this paper proposes a fast-slow method for the fast MD simulation of biomolecules in explicit solvent. This fast-slow method divides the entire system into two parts: a core layer (typically solute or biomolecule) and a peripheral layer (typically solvent molecules). The core layer is treated using standard MD method but the peripheral layer is treated by a slower dynamics method to reduce the computational cost. We compared four different simulation models in testing calculations for several small proteins. These include gas-phase, implicit solvent, fast-slow explicit solvent and standard explicit solvent MD simulations. Our study shows that gas-phase and implicit solvent models do not provide a realistic solvent environment and fail to correctly produce reliable dynamic structures of proteins. On the other hand, the fast-slow method can essentially reproduce the same solvent effect as the standard explicit solvent model while gaining an order of magnitude in efficiency. This fast-slow method thus provides an efficient approach for accelerating the MD simulation of biomolecules in explicit solvent.

Rational Design of P450 aMOx for Improving Anti-Markovnikov Selectivity Based on the “Butterfly” Model

Aromatic aldehydes are important industrial raw materials mainly synthesized by anti-Markovnikov (AM) oxidation of corresponding aromatic olefins. The AM product selectivity remains a big challenge. P450 aMOx is the first reported enzyme that could catalyze AM oxidation of aromatic olefins. Here, we reported a rational design strategy based on the “butterfly” model of the active site of P450 aMOx. Constrained molecular dynamic simulations and a binding energy analysis of key residuals combined with an experimental alanine scan were applied. As a result, the mutant A275G showed high AM selectivity of >99%. The results also proved that the “butterfly” model is an effective design strategy for enzymes.

Ab initio neural network MD simulation of thermal decomposition of a high energy material CL-20/TNT

CL-20 (2,4,6,8,10,12-hexanitro-2,4,6,8,10,12-hexaazaisowurtzitane, also known as HNIW) is one of the most powerful energetic materials. However, its high sensitivity to environmental stimuli greatly reduces its safety and severely limits its application. In this work, ab initio based neural network potential (NNP) energy surfaces for both β-CL-20 and CL-20/TNT co-crystals were constructed. To accurately simulate the thermal decomposition processes of these two crystal systems, reactive molecular dynamics simulations based on the NNPs were performed. Many important intermediate species and their associated reaction paths during the decomposition had been identified in the simulations and the direct results on detonation temperatures of both systems were provided. The simulations also showed clearly that 2,4,6-trinitrotoluene (TNT) molecules in the co-crystal act as a buffer to slow down the chain reactions triggered by nitrogen dioxide and this effect is more significant at lower temperatures. Specifically, the addition of TNT molecules in the CL-20/TNT co-crystal introduces intermolecular hydrogen bonds between CL-20 and TNT molecules in the system, thereby increasing the thermal stability of the co-crystal. The current reactive molecular dynamics simulation is performed based on the NNP which helps in accelerating the speed of ab initio molecular dynamics (AIMD) simulation by more than 3 orders of magnitude while preserving the accuracy of density functional theory (DFT) calculations. This enabled us to perform longer-time simulations at more realistic temperatures that traditional AIMD methods cannot achieve. With the advantage of the NNP in its powerful fitting ability and transferability, the NNP-based MD simulation can be widely applied to energetic material systems.

AA-Score: a New Scoring Function Based on Amino Acid-Specific Interaction for Molecular Docking

The protein-ligand scoring function plays an important role in computer-aided drug discovery and is heavily used in virtual screening and lead optimization. In this study, we developed a new empirical protein-ligand scoring function with amino acid-specific interaction components for hydrogen bond, van der Waals, and electrostatic interactions. In addition, hydrophobic, π-stacking, π-cation, and metal-ligand interactions are also included in the new scoring function. To better evaluate the performance of the AA-Score, we generated several new test sets for evaluation of scoring, ranking, and docking performances, respectively. Extensive tests show that AA-Score performs well on scoring, docking, and ranking as compared to other widely used traditional scoring functions. The performance improvement of AA-Score benefits from the decomposition of individual interaction into amino acid-specific types. To facilitate applications, we developed an easy-to-use tool to analyze protein-ligand interaction fingerprint and predict binding affinity using the AA-Score. The source code and associated running examples can be found at https://github.com/xundrug/AA-Score-Tool.

Mutational Effect of Some Major COVID-19 Variants on Binding of the S Protein to ACE2

COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which has many variants that accelerated the spread of the virus. In this study, we investigated the quantitative effect of some major mutants of the spike protein of SARS-CoV-2 binding to the human angiotensin-converting enzyme 2 (ACE2). These mutations are directly related to the Variant of Concern (VOC) including Alpha, Beta, Gamma, Delta and Omicron. Our calculations show that five major mutations (N501Y, E484K, L452R, T478K and K417N), first reported in Alpha, Beta, Gamma and Delta variants, all increase the binding of the S protein to ACE2 (except K417N), consistent with the experimental findings. We also studied an additional eight mutations of the Omicron variant that are located on the interface of the receptor binding domain (RDB) and have not been reported in other VOCs. Our study showed that most of these mutations (except Y505H and G446S) enhance the binding of the S protein to ACE2. The computational predictions helped explain why the Omicron variant quickly became dominant worldwide. Finally, comparison of several different computational methods for binding free energy calculation of these mutants was made. The alanine scanning method used in the current calculation helped to elucidate the residue-specific interactions responsible for the enhanced binding affinities of the mutants. The results show that the ASGB (alanine scanning with generalized Born) method is an efficient and reliable method for these binding free energy calculations due to mutations.

Discovery of Novel Inhibitors of CDK2 Using Docking and Physics-based Binding Free Energy Calculation

Cyclin-dependent kinase (CDK) is a serine/threonine protein kinase family that cooperates with cyclin and plays an important role in the regulation of cell cycle. Cyclin-dependent kinase 2 is an important member of the CDK family and holds great promise as an anti-cancer drug target. In this study, we used molecular docking and physics-based binding free energy calculation method AS-IE that explicitly calculated protein-ligand binding entropy to discover novel inhibitors of CDK2. A total of 17 inhibitors were discovered with the best IC₅₀ reaching ~2 μM. Decomposition of the binding free energy using AS-IE reveals key protein-ligand interactions that determines the activity. These results provided a good example of drug design using physics-based free energy calculation method such as AS-IE and the novel compounds offered a good start point for further development of CDK2 inhibitors.

Computational Alanine Scanning Reveals Common Features of TCR/pMHC Recognition in HLA-DQ8-Associated Celiac Disease

In HLA-DQ8-associated celiac disease, Gliadin-γ1 or Gliadin-α1 peptide is presented to the cell surface and recognized by several types of T-cell receptor (TCR), but it is still unclear how the TCR, peptide, and the major histocompatibility complex (MHC) act together to trigger celiac disease. For now, most of the analysis is based on static crystal structures. And the detailed information about these structures based on energetic interaction is still lacking. Here, we took four types of celiac disease-related MHC-peptide-TCR structures from three patients to perform computational alanine scanning calculations using the molecular mechanics generalized born surface area (MM/GBSA) approach combined with a recently developed interaction entropy (IE) method to identify the key residues on TCR, peptide, and MHC. Our study aims to shed some light on the interaction mechanism of this complex protein interaction system. Based on detailed computational analysis and mutational calculations, important binding interactions in these triple-interaction complexes are analyzed, and critical residues responsible for TCR/pMHC recognition pattern in HLA-DQ8-associated celiac disease are presented. These detailed analysis and computational result should help shed light on our understanding of the celiac disease and the development of the medical treatment.

Targeting mechanism for SARS-CoV-2 in silico: interaction and key groups of TMPRSS2 toward four potential drugs

The global dissemination of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has seriously endangered human health. The number of confirmed cases is still increasing; however, treatment options are limited. Transmembrane protease serine 2 (TMPRSS2), as a key protease that primes the binding of SARS-CoV-2 spike protein and angiotensin-converting enzyme 2 (ACE2), has become an attractive target and received widespread attention. Thus, four potential drugs (bromhexine, camostat, gabexate, and nafamostat) were used to explore the mechanism of binding with TMPRSS2 in this work. A 65 ns molecular dynamics simulation was performed three times for each drug-TMPRSS2 system for reliable energy calculation and conformational analysis, of which the simulations of nafamostat-TMPRSS2 systems were further extended to 150 ns three times due to the discovery of two binding modes. Through the results of calculating binding free energy by nine methods, the binding affinity of camostat, gabexate, and nafamostat to TMPRSS2 showed great advantages compared with bromhexine, where the nafamostat was surprisingly found to present two reasonable binding conformations (forward and reverse directions). Two negatively charged amino acids (Asp435 and Glu299) can clamp the two positively charged groups (guanidinium group and amidinium group) in either forward or reverse fashion, and the forward one is more stable than the reverse. In addition, compared with gabexate, the dimethylamino group in camostat forms more van der Waals interactions with surrounding hot-spots His296 and Val280, resulting in a stronger affinity to TMPRSS2. For bromhexine, multiple binding sites are displayed in the binding pocket due to its small molecular structure, and van der Waals interactions play the dominant role in the binding process. In particular, six typical hot-spots were identified in the last three serine protease inhibitor systems, i.e., Asp435, Ser436, Gln438, Trp461, Ser463, and Gly464. The guanidinium groups of the drugs have powerful interactions with adjacent residues due to the formation of more hydrogen bonds, suggesting that this may be the critical site for drug design against TMPRSS2. This work provides valuable molecular insight into these four drug-TMPRSS2 binding mechanisms and is helpful for designing and screening drugs targeting TMPRSS2.

Automated Construction of Neural Network Potential Energy Surface: The Enhanced Self-Organizing Incremental Neural Network Deep Potential Method

In recent years, the use of deep learning (neural network) potential energy surface (NNPES) in molecular dynamics simulation has experienced explosive growth as it can be as accurate as quantum chemistry methods while being as efficient as classical mechanic methods. However, the development of NNPES is highly nontrivial. In particular, it has been troubling to construct a dataset that is as small as possible yet can cover the target chemical space. In this work, an ESOINN-DP method is developed, which has the enhanced self-organizing incremental neural network (ESOINN) and a newly proposed error indicator at its core. With ESOINN-DP, one can construct the NNPES with little human intervention, and this method ensures that the constructed reference dataset covers the target chemical space with minimum redundancy. The performance of the ESOINN-DP method has been well validated by developing neural network potential energy surfaces for water clusters, tripeptides, and by de-redundancy of a sub-dataset of the ANI-1 database. We believe that the ESOINN-DP method provides a novel idea for the construction of NNPES and, especially, the reference datasets, and it can be used for molecular dynamics (MD) simulations of various gas-phase and condensed-phase chemical systems.

Computational Insights into the Binding Mechanism of OxyS sRNA with Chaperone Protein Hfq

Under the oxidative stress condition, the small RNA (sRNA) OxyS that acts as essential post-transcriptional regulators of gene expression is produced and plays a regulatory function with the assistance of the RNA chaperone Hfq protein. Interestingly, experimental studies found that the N48A mutation of Hfq protein could enhance the binding affinity with OxyS while resulting in the defection of gene regulation. However, how the Hfq protein interacts with sRNA OxyS and the origin of the stronger affinity of N48A mutation are both unclear. In this paper, molecular dynamics (MD) simulations were performed on the complex structure of Hfq and OxyS to explore their binding mechanism. The molecular mechanics generalized born surface area (MM/GBSA) and interaction entropy (IE) method were combined to calculate the binding free energy between Hfq and OxyS sRNA, and the computational result was correlated with the experimental result. Per-residue decomposition of the binding free energy revealed that the enhanced binding ability of the N48A mutation mainly came from the increased van der Waals interactions (vdW). This research explored the binding mechanism between Oxys and chaperone protein Hfq and revealed the origin of the strong binding affinity of N48A mutation. The results provided important insights into the mechanism of gene expression regulation affected by protein mutations.

A fixed multi-site interaction charge model for an accurate prediction of the QM/MM interactions

A fixed multi-site interaction charge (FMIC) model was proposed for the accurate prediction of intermolecular electrostatic interactions based on the quantum mechanical linear response of a molecule to an external electric field. In such a model, some additional off-center interaction sites were added for capturing multipole interactions for a given molecule. By multivariate least-square fitting analysis of the calculated QM/MM interactions of a given molecule with the electrostatic environment and the electrostatic potentials of the environment at the pre-defined distributed interaction sites, the FMIC of the molecule was obtained. The model system of CO in myoglobin (Mb) was utilized to demonstrate the derivation of the FMIC. The accuracy of FMIC in predicting the electrostatic interactions between CO and the Mb environment was investigated using 10 000 different Mb-CO configurations generated from the 400 ps QM/MM MD simulation. In comparison to the QM/MM calculations at the B3LYP/aug-cc-pVTZ/ff99SB level, the mean unsigned error (MUE) of the results based on the FMIC model was merely 0.10 kcal mol^-1, and the root mean square error (RMSE) was only 0.13 kcal mol^-1, which are significantly lower than the results predicted by the ESP charge model (MUE = 1.45 kcal mol^-1, and RMSE = 1.7 kcal mol^-1, respectively). The transferability of FMIC was tested by applying the obtained FMIC in the wild type Mb-CO system to the mutants of V68F and H64L Mb-CO systems. The MUEs of the obtained results for 10 000 different configurations are both smaller than 0.2 kcal mol^-1 for the V68F and H64L Mb-CO systems in comparison to the B3LYP/aug-cc-pVTZ/ff99SB calculations, and the RMSEs are also lower than 0.2 kcal mol^-1 for both mutants. The applications of FMIC were extended to model the electrostatic interactions between a hydrogen fluoride molecule and 492 waters in a truncated octahedron box; our study showed that the FMIC could give satisfactory results with a MUE of 0.12 kcal mol^-1 and a RMSE of 0.16 kcal mol^-1 in comparison to the B3LYP/aug-cc-pVDZ/TIP3P calculations for 10 000 different configurations generated using the 10 ns classical MD simulation. Therefore, the FMIC method provides an accurate and efficient tool for predicting intermolecular electrostatic interactions, which can be utilized in the future development of molecular force fields.

Residue-specific binding mechanisms of PD-L1 to its monoclonal antibodies by computational alanine scanning

Programmed cell death 1 receptor (PD-1) on the surface of T cells and its ligand 1 (PD-L1) are immune checkpoint proteins. Treating cancer patients with inhibitors blocking this checkpoint has significantly prolonged the survival rate of patients. In this study, we examined several monoclonal antibodies (mAbs) of PD-L1 and studied their detailed binding mechanism to PD-L1. An efficient computational alanine scanning method was used to perform quantitative analysis of hotspot residues that are important for PD-1/PD-L1 binding. A total of five PD-L1/mAb complexes were investigated and hotspots on both PD-L1 and mAbs were predicted. Our result shows that PD-L1M115 and PD-L1Y123 are two relatively important hotspots in all the five PD-L1/mAb binding complexes. It is also found that the important residues of mAbs binding to PD-L1M115 and PD-L1Y123 are similar to each other. The computational alanine scanning result is compared to the experimental measurements that are available for two of the mAbs (KN035 and atezolizumab). The calculated alanine scanning result is in good agreement with the experimental data with a correlation coefficient of 0.87 for PD-L1/KN035 and 0.6 for PD-L1/atezolizumab. Our computation found more hotspots on PD-L1 in the PD-L1/KN035 complex than those in the PD-L1/atezolizumab system, indicating stronger binding affinity in the former than the latter, which is in good agreement with the experimental finding. The present work provides important insights for the design of new mAbs targeting PD-L1.

Molecular mechanism related to the binding of fluorophores to Mango-II revealed by multiple-replica molecular dynamics simulations

Recently, RNA aptamers activating small-molecule fluorophores have been successfully applied to tag and track RNAs in vivo. It is of significance to investigate the molecular mechanism of the fluorophore-RNA aptamer bindings at the atomic level to seek a possible pathway to enhance the fluorescence efficiency of fluorophores. In this work, multiple replica molecular dynamics (MRMD) simulations, essential dynamics (ED) analysis, and hierarchical clustering analysis were coupled to probe the effect of A22U mutation on the binding of two fluorophores, TO1-Biotin (TO1) and TO3-Biotin (TO3), to the Mango-II RNA aptamer (Mango-II). ED analysis reveals that A22U induces alterations in the binding pocket and sites of TO1 and TO3 to the Mango-II, which in turn tunes the fluorophore-RNA interface and changes the interactions of TO1 and TO3 with separate nucleotides of Mango-II. Dynamics analyses also uncover that A22U exerts the opposite impact on the molecular surface areas of the Mango-II and sugar puckers of nucleotides 22 and 23 in Mango-II complexed with TO1 and TO3. Moreover, the calculations of binding free energies suggest that A22U strengthens the binding ability of TO1 to the mutated Mango-II but weakens TO3 to the mutated Mango-II when compared with WT. These findings imply that point mutation in nucleotides possibly tune the fluorescence of fluorophores binding to RNA aptamers, providing a possible scheme to enhance the fluorescence of fluorophores.

Anchor-Locker Binding Mechanism of the Coronavirus Spike Protein to Human ACE2: Insights from Computational Analysis

COVID-19 has emerged as the most serious international pandemic in early 2020 and the lack of comprehensive knowledge in the recognition and transmission mechanisms of this virus hinders the development of suitable therapeutic strategies. The specific recognition during the binding of the spike glycoprotein (S protein) of coronavirus to the angiotensin-converting enzyme 2 (ACE2) in the host cell is widely considered the first step of infection. However, detailed insights on the underlying mechanism of dynamic recognition and binding of these two proteins remain unknown. In this work, molecular dynamics simulation and binding free energy calculation were carried out to systematically compare and analyze the receptor-binding domain (RBD) of six coronavirus' S proteins. We found that affinity and stability of the RBD from SARS-CoV-2 under the binding state with ACE2 are stronger than those of other coronaviruses. The solvent-accessible surface area (SASA) and binding free energy of different RBD subunits indicate an "anchor-locker" recognition mechanism involved in the binding of the S protein to ACE2. Loop 2 (Y473-F490) acts as an anchor for ACE2 recognition, and Loop 3 (G496-V503) locks ACE2 at the other nonanchoring end. Then, the charged or long-chain residues in the β-sheet 1 (N450-F456) region reinforce this binding. The proposed binding mechanism was supported by umbrella sampling simulation of the dissociation process. The current computational study provides important theoretical insights for the development of new vaccines against SARS-CoV-2.

MolGpka: A Web Server for Small Molecule pKa Prediction Using a Graph-Convolutional Neural Network

pK_a is an important property in the lead optimization process since the charge state of a molecule in physiologic pH plays a critical role in its biological activity, solubility, membrane permeability, metabolism, and toxicity. Accurate and fast estimation of small molecule pK_a is vital during the drug discovery process. We present MolGpKa, a web server for pK_a prediction using a graph-convolutional neural network model. The model works by learning pK_a related chemical patterns automatically and building reliable predictors with learned features. ACD/pK_a data for 1.6 million compounds from the ChEMBL database was used for model training. We found that the performance of the model is better than machine learning models built with human-engineered fingerprints. Detailed analysis shows that the substitution effect on pK_a is well learned by the model. MolGpKa is a handy tool for the rapid estimation of pK_a during the ligand design process. The MolGpKa server is freely available to researchers and can be accessed at https://xundrug.cn/molgpka.

Quantitative analysis of ACE2 binding to coronavirus spike proteins: SARS-CoV-2 vs. SARS-CoV and RaTG13

The global outbreak of the COVID-19 pandemic is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Bat virus RaTG13 and SARS-CoV are also members of the coronavirus family and SARS-CoV caused a world-wide pandemic in 2003. SARS-CoV-2, SARS-CoV and RaTG13 bind to angiotensin-converting enzyme 2 (ACE2) through their receptor-binding domain (RBD) of the spike protein. SARS-CoV-2 binds ACE2 with a higher binding affinity than SARS-CoV and RaTG13. Here we performed molecular dynamics simulation of these binding complexes and calculated their binding free energies using a computational alanine scanning method. Our MD simulation and hotspot residue analysis showed that the lower binding affinity of SARS-CoV to ACE2 vs. SARS-CoV-2 to ACE2 can be explained by different hotspot interactions in these two systems. We also found that the lower binding affinity of RaTG13 to ACE2 is mainly due to a mutated residue (D501) which resulted in a less favorable complex formation for binding. We also calculated an important mutation of N501Y in SARS-CoV-2 using both alanine scanning calculation and a thermodynamic integration (TI) method. Both calculations confirmed a significant increase of the binding affinity of the N501Y mutant to ACE2 and explained its molecular mechanism. The present work provides an important theoretical basis for understanding the molecular mechanism in coronavirus spike protein binding to human ACE2.

Automatically Constructed Neural Network Potentials for Molecular Dynamics Simulation of Zinc Proteins

The development of accurate and efficient potential energy functions for the molecular dynamics simulation of metalloproteins has long been a great challenge for the theoretical chemistry community. An artificial neural network provides the possibility to develop potential energy functions with both the efficiency of the classical force fields and the accuracy of the quantum chemical methods. In this work, neural network potentials were automatically constructed by using the ESOINN-DP method for typical zinc proteins. For the four most common zinc coordination modes in proteins, the potential energy, atomic forces, and atomic charges predicted by neural network models show great agreement with quantum mechanics calculations and the neural network potential can maintain the coordination geometry correctly. In addition, MD simulation and energy optimization with the neural network potential can be readily used for structural refinement. The neural network potential is not limited by the function form and complex parameterization process, and important quantum effects such as polarization and charge transfer can be accurately considered. The algorithm proposed in this work can also be directly applied to proteins containing other metal ions.

DeepBSP-a Machine Learning Method for Accurate Prediction of Protein-Ligand Docking Structures

In recent years, machine-learning-based scoring functions have significantly improved the scoring power. However, many of these methods do not perform well in distinguishing the native structure from docked decoy poses due to the lack of decoy structural information in their training data. Here, we developed a machine-learning model, named DeepBSP, that can directly predict the root mean square deviation (rmsd) of a ligand docking pose with reference to its native binding pose. Unlike the binding affinity, the rmsd between the docking poses with reference to their native structures can be straightforwardly determined. By training on a generated data set with 11,925 native complexes and more than 165,000 docked poses, our model shows excellent docking power on our test set and also on the CASF-2016 docking decoy set compared to other major scoring functions. Thus, by combining molecular dockings that generate many poses with the application of DeepBSP, one can more accurately predict the best binding pose that is closest to the native complex structure. This DeepBSP model shall be very useful in picking out poses close to their natives from many poses generated from a dock application.

Inhibition mechanism and hot-spot prediction of nine potential drugs for SARS-CoV-2 Mpro by large-scale molecular dynamic simulations combined with accurate binding free energy calculations

Coronavirus disease 2019 (COVID-19), which is caused by a new coronavirus known as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is spreading around the world. However, a universally effective treatment regimen has not been developed to date. The main protease (M^pro), a key enzyme of SARS-CoV-2, plays a crucial role in the replication and transcription of this virus in cells and has become the ideal target for rational antiviral drug design. In this study, we performed molecular dynamics simulations three times for these complexes of M^pro (monomeric and dimeric) and nine potential drugs that have a certain effect on the treatment of COVID-19 to explore their binding mechanism. In addition, a total of 12 methods for calculating binding free energy were employed to determine the optimal drug. Ritonavir, Arbidol, and Chloroquine consistently showed an outstanding binding ability to monomeric M^pro under various methods. Ritonavir, Arbidol, and Saquinavir presented the best performance when binding to a dimer, which was independent of the protonated state of Hie41 (protonated at N_ε) and Hid41 (protonated at N_δ), and these findings suggest that Chloroquine may not effectively inhibit the activity of dimeric M^proin vivo. Furthermore, three common hot-spot residues of Met165, Hie41, and Gln189 of monomeric M^pro systems dominated the binding of Ritonavir, Arbidol, and Chloroquine. In dimeric M^pro, Gln189, Met165, and Met49 contributed significantly to binding with Ritonavir, Arbidol, and Saquinavir; therefore, Gln189 and Met165 might serve as the focus in the discovery and development of anti-COVID-19 drugs. In addition, the van der Waals interaction played a significant role in the binding process, and the benzene ring of the drugs showed an apparent inhibitory effect on the normal function of M^pro. The binding cavity had great flexibility to accommodate these different drugs. The results would be notably helpful for enabling a detailed understanding of the binding mechanisms for these important drug-M^pro interactions and provide valuable guidance for the design of potent inhibitors.

Analysis of the binding modes of the first- and second-generation antiandrogens with respect to F876L mutation

Androgen receptor (AR) is an important target for the treatment of prostate cancer, and mutations in the AR have an important impact on the resistance of existing drugs. In this work, we performed molecular dynamics simulations of the existing marketed antiandrogens flutamide, nilutamide, bicalutamide, enzalutamide, apalutamide, darolutamide, and its main metabolite ORM15341 in complex with the wild-type and F876L mutant AR. We calculated the residue-specific binding free energy contribution of the wild-type and mutant ARs with the AS-IE method and analyzed the hotspot residues and the binding free energy contributions of specific residues before and after the mutation. In addition, we analyzed the total binding obtained by adding residue binding energy contributions and compared the results with experimental values. The obtained residue-specific binding information should be very helpful in understanding the mechanism of drug resistance with respect to specific mutations and in the design of new generation drugs against possible new mutations.

An electrostatic energy-based charge model for molecular dynamics simulation

The interactions of the polar chemical bonds such as C=O and N-H with an external electric field were investigated, and a linear relationship between the QM/MM interaction energies and the electric field along the chemical bond is established in the range of weak to intermediate electrical fields. The linear relationship indicates that the electrostatic interactions of a polar group with its surroundings can be described by a simple model of a dipole with constant moment under the action of an electric field. This relationship is employed to develop a general approach to generating an electrostatic energy-based charge (EEC) model for molecules containing single or multiple polar chemical bonds. Benchmark test studies of this model were carried out for (CH₃)₂-CO and N-methyl acetamide in explicit water, and the result shows that the EEC model gives more accurate electrostatic energies than those given by the widely used charge model based on fitting to the electrostatic potential (ESP) in direct comparison to the energies computed by the QM/MM method. The MD simulations of the electric field at the active site of ketosteroid isomerase based on EEC demonstrated that EEC gave a better representation of the electrostatic interaction in the hydrogen-bonding environment than the Amber14SB force field by comparison with experiment. The current study suggests that EEC should be better suited for molecular dynamics study of molecular systems with polar chemical bonds such as biomolecules than the widely used ESP or RESP (restrained ESP) charge models.

Engineering the biomimetic cofactors of NMNH for cytochrome P450 BM3 based on binding conformation refinement

Cytochrome P450 BM3 (BM3) is an important oxidoreductase that is widely used in drug synthesis, chemical synthesis, and other industries. However, as BM3 unquestionably increases costs by consuming a natural cofactor that unstably provides electrons, an alternative biomimetic cofactor with simpler structures represented by nicotinamide mononucleotide (NMNH) has been utilized. Currently, few reports exist on artificially modified BM3 enzymes using NMNH, especially regarding theoretical simulation and calculation. With the cognition of the mechanism in mind, we propose a strategy that optimizes and refines catalytic conformation. Based on constrained molecular dynamics simulation, the distance between N-5 of FAD flavin and C-4 of NMNH is used as a cue for the determination of improved conformation, and the potential positive mutants are subsequently screened virtually in accordance with binding free energy requirements. As a result, the K_cat/K_M values of the favorable mutant S848R increased to 205.38% compared to the wild-type BM3 with NMNH. These data indicate that our strategy can be applied for the specific utilization of biomimetic cofactors by oxidoreductases represented by BM3.

A rational design strategy was proposed to improve the efficient utilization of alternative biomimetic cofactor by P450 BM3 enzyme.

Computational Analysis of Residue-Specific Binding Free Energies of Androgen Receptor to Ligands

Androgen receptor (AR) is an important therapeutic target for the treatment of diseases such as prostate cancer, hypogonadism, muscle wasting, etc. In this study, the complex structures of the AR ligand-binding domain (LBD) with fifteen ligands were analyzed by molecular dynamics simulations combined with the alanine-scanning-interaction-entropy method (ASIE). The quantitative free energy contributions of the pocket residues were obtained and hotspot residues are quantitatively identified. Our calculation shows that that these hotspot residues are predominantly hydrophobic and their interactions with binding ligands are mainly van der Waals interactions. The total binding free energies obtained by summing over binding contributions by individual residues are in good correlation with the experimental binding data. The current quantitative analysis of binding mechanism of AR to ligands provides important insight on the design of future inhibitors.