Implicit nonpolar solvent models

Applying Spatiotemporal Modeling of Cell Dynamics to Accelerate Drug Development

Cells act as physical computational programs that utilize input signals to orchestrate molecule-level protein-protein interactions (PPIs), generating and responding to forces, ultimately shaping all of the physiological and pathophysiological behaviors. Genome editing and molecule drugs targeting PPIs hold great promise for the treatments of diseases. Linking genes and molecular drugs with protein-performed cellular behaviors is a key yet challenging issue due to the wide range of spatial and temporal scales involved. Building predictive spatiotemporal modeling systems that can describe the dynamic behaviors of cells intervened by genome editing and molecular drugs at the intersection of biology, chemistry, physics, and computer science will greatly accelerate pharmaceutical advances. Here, we review the mechanical roles of cytoskeletal proteins in orchestrating cellular behaviors alongside significant advancements in biophysical modeling while also addressing the limitations in these models. Then, by integrating generative artificial intelligence (AI) with spatiotemporal multiscale biophysical modeling, we propose a computational pipeline for developing virtual cells, which can simulate and evaluate the therapeutic effects of drugs and genome editing technologies on various cell dynamic behaviors and could have broad biomedical applications. Such virtual cell modeling systems might revolutionize modern biomedical engineering by moving most of the painstaking wet-laboratory effort to computer simulations, substantially saving time and alleviating the financial burden for pharmaceutical industries.

Identification of novel brain penetrant GSK-3β inhibitors toward Alzheimer's disease therapy by virtual screening, molecular docking, dynamic simulation, and MMPBSA analysis

One of the significant therapeutic targets for Alzheimer's disease (AD) is Glycogen Synthase Kinase-3β (GSK-3β). Inhibition of GSK-3β can prevent hyperphosphorylation of tau, and thus prevent formation and accumulation of neurofibrillary tangles and neuropil threads that block intracellular transport, trigger unfolded protein response, and increase oxidative stress, cumulatively leading to neurodegeneration. In this study, we have performed structure-based virtual screening of two small-molecule libraries from ChemDiv CNS databases using AutoDock Vina to identify hit molecules based on their binding affinities compared to that of an established GSK-3β inhibitor, indirubin-3'-monoxime (IMO). Pharmacoinformatic screening on SwissADME and pkCSM servers enabled identification of lead molecules with favorable pharmacoinformatic properties for drug likeliness, including blood brain barrier (BBB) permeability. Further, molecular dynamic simulations identified six candidate lead molecules that show stable complex formation with GSK-3β in dynamic state under physiological conditions. Principal component analysis of the dynamic state was used to plot Free Energy Landscapes (FELs) of GSK-3β-ligand complexes. STRIDE secondary structure analysis of the lowest energy conformations identified from FEL plots, and binding free energy calculations by Molecular Mechanics Poisson-Boltzmann Surface Area ((ΔGbind)MM-PBSA) of the simulation trajectories led to the identification of two ligands as potential lead molecules having favorable free energy landscape profiles as well as significantly lower (ΔGbind)MM-PBSA in dynamic state compared to that of reference inhibitor IMO. Hence, this study identifies two novel brain penetrant GSK-3β inhibitors that are likely to have therapeutic potential against Alzheimer's disease.

Leukotriene B4 receptor 1 (BLT1) activation by leukotriene B4 (LTB4) and E resolvins (RvE1 and RvE2)

Leukotriene B4 (LTB4) is a lipid inflammatory mediator derived from arachidonic acid (AA). Leukotriene B4 receptor 1 (BLT1), a G protein-coupled receptor (GPCR), is a receptor of LTB4. Nonetheless, the resolution of inflammation is driven by specialized pro-resolving lipid mediators (SPMs) such as resolvins E1 (RvE1) and E2 (RvE2). Both resolvins are derived from omega-3 fatty acid eicosapentaenoic acid (EPA). Here, long-term molecular dynamics simulations (MD) were performed to investigate the activation of the BLT1 receptor using two pro-resolution agonists (RvE1 and RvE2) and an inflammatory agonist (LTB4). We have analyzed the receptor's activation state, electrostatic interactions, and the binding affinity the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approach. The results showed that LTB4 and RvE1 have kept the receptor in an active state by higher simulation time. MD showed that the ligand-receptor interactions occurred mainly through residues H94, R156, and R267. The MMPBSA calculations showed residues R156 and R267 were the two mainly hotspots. Our MMPBSA results were compatible with experimental results from other studies. Overall, the results from this study provide new insights into the activation mechanisms of the BLT1 receptor, reinforcing the role of critical residues and interactions in the binding of pro-resolution and pro-inflammatory agonists.

Structural insights into the allosteric effects of the antiepileptic drug topiramate on the CaV2.3 channel

The R-type voltage-gated calcium channel CaV2.3 is predominantly located in the presynapse and is implicated in distinct types of epileptic seizures. It has consequently emerged as a molecular target in seizure treatment. Here, we determined the cryo-EM structure of the CaV2.3-α2δ1-β1 complex in the topiramate-bound state at a 3.0 Å resolution. We provide a snapshot of the binding site of topiramate, a widely prescribed antiepileptic drug, on a voltage-gated ion channel. The binding site is located at an intracellular juxtamembrane hydrophilic cavity. Further structural analysis revealed that topiramate may allosterically facilitate channel inactivation. These findings provide fundamental insights into the mechanism underlying the inhibitory effect of topiramate on CaV and NaV channels, elucidating a previously unseen modulator binding site and thus pointing toward a route for the development of new drugs.

Hybrid Cluster-Continuum Method for Single-Ion Solvation Free Energy in Acetonitrile Solvent

A new hybrid discrete-continuum approach named the cluster-continuum static approximation (CCSA) has been proposed for acetonitrile solvent. The continuum part uses the conductor-like polarizable continuum model for electrostatic and a surface area-dependent term for nonelectrostatic solvation. The CCSA includes only one explicit acetonitrile solvent molecule and a damping function, which makes the CCSA method reduce to pure continuum solvation in the case of weaker potential of mean force for solute-solvent interaction. The performance of the model was tested for 22 anions and 22 cations, including challenge species that cannot be adequately described by pure continuum solvation. A comparison was done with the widely used solvent model density (SMD) model. For anions, the CCSA reduces to pure continuum solvation and the method has the same performance as the SMD model, with a standard deviation of the mean signed error (SD-MSE) of 2.7 kcal mol-1 for both models. However, the CCSA method for cations considerably outperforms the SMD model, with an SD-MSE of 3.3 kcal mol-1 for the former and 8.4 kcal mol-1 for the latter. The method can be automated, and the present study suggests that continuum solvation models could be parameterized taking into account the explicit solvation as proposed in this work.

Structural mechanism underlying variations in DNA binding by the androgen receptor

The androgen receptor (AR) is a steroid activated transcription factor which recognizes DNA motifs resembling inverted repeats of a conserved 5'-AGAACA-3'-like hexanucleotides separated by a three-nucleotide spacer from a similar, but less conserved hexanucleotide. Here, we report the structures of the human AR DNA binding domain (DBD) bound to two natural AREs (C3 and MTV) in head-to-head dimer conformations, diffracting at 2.05 Å and 2.25 Å, respectively. These structures help to explain the impact of androgen insensitivity mutations on the structure integrity, DNA binding and DBD dimerization. The binding affinity of the AR DBD to different DNA motifs were measured by the BioLayer Interferometry (BLI) and further validated by Molecular Dynamics (MD) simulations. This shows that the high binding affinity of the first DBD to the upstream 5'-AGAACA-3' motif induces the cooperative binding of the second DBD to the second hexanucleotide. Our data indicate identical interaction of the DBDs to the upstream hexanucleotides, while forming an induced closer contact of the second DBD on the non-canonical hexanucleotides. The variation in binding between the DBD monomers are the result of differences in DNA occupancy, protein-protein interactions, DNA binding affinity, and DNA binding energy profiles. We propose this has functional consequences.

Exploring transmembrane allostery in the MexB: DB08385 variant as a promising inhibitor-like candidate against Pseudomonas aeruginosa antibiotic resistance: a computational study

Pseudomonas aeruginosa, a formidable pathogen renowned for its antimicrobial resistance, poses a significant threat to immunocompromised individuals. In this regard, the MexAB-OprM efflux pump acts as a pivotal line of defense by extruding antimicrobials from bacterial cells. The inner membrane homotrimeric protein MexB captures antibiotics and translocates them into the outer membrane OprM channel protein connected through the MexA adaptor protein. Despite extensive efforts, competitive inhibitors targeting the tight (T) protomer of the MexB protein have not received FDA approval for medical use. Over the past few years, allosteric inhibitors have become popular as alternatives to the classical competitive inhibitor-based approach because of their higher specificity, lower dosage, and reduced toxicological effects. Hence, in this study, we unveiled the existence of a transmembrane allosteric binding pocket of MexB inspired by the recent discovery of an important allosteric inhibitor, BDM88855, for the homolog AcrB protein. While repurposing BDM88855 proved ineffective in controlling the MexB loose (L) protomer, our investigation identified a promising alternative: a chlorine-containing variant of DB08385 (2-Cl DB08385 or Variant 1). Molecular dynamics simulations, including binding free energy estimation coupled with heterogeneous dielectric implicit membrane model (implicit-membrane MM/PBSA), interaction entropy (IE) analysis and potential of mean force (PMF) calculation, demonstrated Variant 1's superior binding affinity to the transmembrane pocket, displaying the highest energy barrier in the ligand unbinding process. To elucidate the allosteric crosstalk between the transmembrane and porter domain of MexB, we employed the 'eigenvector centrality' measure in the linear mutual information obtained from the protein correlation network. Notably, this study confirmed the presence of an allosteric transmembrane site in the MexB L protomer. In addition to this, Variant 1 emerged as a potent regulator of allosteric crosstalk, inducing an 'O-L intermediate state' in the MexB L protomer. This induced state might hold the potential to diminish substrate intake into the access pocket, leading to the ineffective efflux of antibiotics.

Decoding Molecular Mechanism Underlying Human Olfactory Receptor OR8D1 Activation by Sotolone Enantiomers

Sotolone, a chiral compound, plays an important role in the food industry. Herein, (R)-/(S)-sotolone were separated to determine their odor characteristics and thresholds in air (R-form: smoky, burned, herb, and green aroma, 0.0514 μg/m3; S-form: sweet, milk, acid, and nutty aroma, 0.0048 μg/m3). OR8D1 responses to (R)-/(S)-sotolone were detected in a HEK293 cell-based luminescence assay. (S)-Sotolone was a more potent agonist than (R)-sotolone (EC50 values of 84.98 ± 1.05 and 167.20 ± 0.25 μmol/L, respectively). Molecular dynamics simulations and molecular mechanics Poisson-Boltzmann surface area analyses confirmed that the combination of (S)-sotolone and OR8D1 was more stable than that of (R)-sotolone. Odorant docking, multiple sequence alignments, site-directed mutagenesis, and functional studies with recombinant odorant receptors (ORs) in a cell-based luminescence assay identified 11 amino-acid residues that influence the enantioselectivity of OR8D1 toward sotolone significantly and that N2065.46 was indispensable to the activation of OR8D1 by (S)-sotolone.

Interferon- as a Potential Inhibitor of SARS-CoV-2 ORF6 Accessory Protein

The ORF6 protein of the SARS-CoV-2 virus plays a crucial role in blocking the innate immune response of the infected cells by inhibiting interferon pathways. Additionally, it binds to and immobilises the RAE1 protein on the cytoplasmic membranes, thereby blocking mRNA transport from the nucleus to the cytoplasm. In all these cases, the host cell proteins are tethered by the flexible C-terminus of ORF6. A possible strategy to inhibit the biological activity of ORF6 is to bind its C-terminus with suitable ligands. Our in silico experiments suggest that hIFNγ binds the ORF6 protein with high affinity, thus impairing its interactions with RAE1 and, consequently, its activity in viral invasion. The in vitro studies reported here reveal a shift of the localisation of RAE1 in ORF6 overexpressing cells upon treatment with hIFNγ from predominantly cytoplasmic to mainly nuclear, resulting in the restoration of the export of mRNA from the nucleus. We also explored the expression of GFP in transfected-with-ORF6 cells by means of fluorescence microscopy and qRT-PCR, finding that treatment with hIFNγ unblocks the mRNA trafficking and reinstates the GFP expression level. The ability of the cytokine to block ORF6 is also reflected in minimising its negative effects on DNA replication by reducing accumulated RNA-DNA hybrids. Our results, therefore, suggest hIFNγ as a promising inhibitor of the most toxic SARS-CoV-2 protein.

Malonyl-CoA is a conserved endogenous ATP-competitive mTORC1 inhibitor

Cell growth is regulated by the mammalian/mechanistic target of rapamycin complex 1 (mTORC1), which functions both as a nutrient sensor and a master controller of virtually all biosynthetic pathways. This ensures that cells are metabolically active only when conditions are optimal for growth. Notably, although mTORC1 is known to regulate fatty acid biosynthesis, how and whether the cellular lipid biosynthetic capacity signals back to fine-tune mTORC1 activity remains poorly understood. Here we show that mTORC1 senses the capacity of a cell to synthesise fatty acids by detecting the levels of malonyl-CoA, an intermediate of this biosynthetic pathway. We find that, in both yeast and mammalian cells, this regulation is direct, with malonyl-CoA binding to the mTOR catalytic pocket and acting as a specific ATP-competitive inhibitor. When fatty acid synthase (FASN) is downregulated/inhibited, elevated malonyl-CoA levels are channelled to proximal mTOR molecules that form direct protein-protein interactions with acetyl-CoA carboxylase 1 (ACC1) and FASN. Our findings represent a conserved and unique homeostatic mechanism whereby impaired fatty acid biogenesis leads to reduced mTORC1 activity to coordinately link this metabolic pathway to the overall cellular biosynthetic output. Moreover, they reveal the existence of a physiological metabolite that directly inhibits the activity of a signalling kinase in mammalian cells by competing with ATP for binding.

Revisiting MMPBSA by Adoption of MC-Based Surface Area/Volume, ANI-ML Potentials, and Two-Valued Interior Dielectric Constant

Here, we report the accuracy improvements of molecular mechanics Poisson-Boltzmann surface area (MMPBSA) calculations by adoption of ANI-ML potentials in replacement of MM terms, the use of solvent-accessible surface area (SASA) and volume (SAV) values from the Monte Carlo sampling of the probe, and introducing two different interior dielectric constants for electrostatic interactions of protein-ligand (P-L) and polar solvation term in the MMPBSA calculations. Our results show that the Pearson correlation coefficients of MMPBSA-calculated values with respect to experimental binding free energies can be drastically improved from 0.48 to 0.90 by adoption of ANI-ML potentials in replacement of MM energy terms in the equation, referred to as ANI-PBSA. Moreover, we show that the SASA/SAV-combined equation in the scaled particle theory (SPT) can be a better choice to model nonpolar solvation term, reaching nearly the same accuracy by ANI-PBSA calculations. Finally, we introduce two different values of interior dielectric constants, which could be an alternative strategy between the single and variable constant definitions.

In silico study of HASDI (high-affinity selective DNA intercalator) as a new agent capable of highly selective recognition of the DNA sequence

Cancer as an acquired genetic disease is based on changes both in the genome itself and in transcription processes. Accordingly, it is at the DNA level that it makes sense to search for and design agents capable of effective and selective anticancer action. In this study, we used an iterative approach based on a molecular dynamics simulation to design a highly selective DNA-intercalating agent called HASDI. To confirm its selective affinity to DNA, we conducted two simulation experiments: HASDI in a complex with a DNA fragment of the EBNA1 gene (it targets 16 nucleotide pairs of this gene) and HASDI in a complex with a random DNA fragment of the KCNH2 gene. The molecular dynamics simulation was carried out in the GROMACS 2019 package. The binding energy was calculated by gmx_MMPBSA 1.5.2. The further analysis was performed using the built-in utilities of GROMACS, gmx_MMPBSA and also XMGRACE and Pymol 1.8. As a result, we determined that the EBNA1-50nt/HASDI complex was stable throughout the whole simulation trajectory. HASDI, due to the presence of a linker modified depending on a specific pair of nitrogenous bases, formed an average of 32 hydrogen bonds with a sequence of 16 nucleotide pairs. Phenazine rings were stably intercalated every 2 base pairs. The root-mean-square deviation of HASDI in such a complex fluctuated around the value of 6.5 Å and had no tendency to increase. The calculated value of the binding free energy was - 235.3 ± 7.77 kcal/mol. The KCNH2-50nt/HASDI complex, as an example of the intercalation of the designed structure into a random part of the human genome, maintained the stability of its position at a level comparable to the EBNA1-50nt/HASDI complex. The phenazine rings were constantly intercalated in their original positions, and the root-mean-square deviation fluctuated around one value, although it had a tendency to chaotic changes. At the same time, this complex was characterized by 17-19 hydrogen bonds, on average, and the binding free energy was - 193.47 ± 14.09 kcal/mol. Moreover, the DNA duplex had local single-nucleotide melting in the region of the 4th linker. According to a significant decrease in the number of hydrogen bonds, a decrease in energy gain, as well as a decrease in the stability of the DNA duplex characteristic of the KCNH2-50nt/HASDI complex compared to the target EBNA1-50nt/HASDI complex, the molecule we designed can be considered a potentially selective DNA polyintercalating agent capable of relatively accurate recognition of 16 base pairs.

Conserved Conformational Dynamics Reveal a Key Dynamic Residue in the Gatekeeper Loop of Human Cyclophilins

Cyclophilins are ubiquitous human enzymes that catalyze peptidyl-prolyl cis-trans isomerization in protein substrates. Of the 17 unique isoforms, five closely related isoforms (CypA-E) are found in various environments and participate in diverse cellular processes, yet all have similar structures and the same core catalytic function. The question is what key residues are behind the conserved function of these enzymes. Here, conformational dynamics are compared across these isoforms to detect conserved dynamics essential for the catalytic activity of cyclophilins. A set of key dynamic residues, defined by the most dynamically conserved positions, are identified in the gatekeeper 2 region. The highly conserved glycine (Gly80) in this region is predicted to underlie the local flexibility, which is further tested by molecular dynamics simulations performed on mutants (G80A) of CypE and CypA. The mutation leads to decreased flexibility of CypE and CypA during substrate binding but increased flexibility during catalysis. Dynamical changes occur in the mutated region and a distal loop downstream of the mutation site in sequence. Examinations of the mutational effect on catalysis show that both mutated CypE and CypA exhibit shifted binding free energies of the substrate under distinct isomer conformations. The results suggest a loss of function in the mutated CypE and CypA. These catalytic changes by the mutation are likely independent of the substrate sequence, at least in CypA. Our work presents a method to identify function-related key residues in proteins.

EPISOL: A software package with expanded functions to perform 3D-RISM calculations for the solvation of chemical and biological molecules

Integral equation theory (IET) provides an effective solvation model for chemical and biological systems that balances computational efficiency and accuracy. We present a new software package, the expanded package for IET-based solvation (EPISOL), that performs 3D-reference interaction site model (3D-RISM) calculations to obtain the solvation structure and free energies of solute molecules in different solvents. In EPISOL, we have implemented 22 different closures, multiple free energy functionals, and new variations of 3D-RISM theory, including the recent hydrophobicity-induced density inhomogeneity (HI) theory for hydrophobic solutes and ion-dipole correction (IDC) theory for negatively charged solutes. To speed up the convergence and enhance the stability of the self-consistent iterations, we have introduced several numerical schemes in EPISOL, including a newly developed dynamic mixing approach. We show that these schemes have significantly reduced the failure rate of 3D-RISM calculations compared to AMBER-RISM software. EPISOL consists of both a user-friendly graphic interface and a kernel library that allows users to call its routines and adapt them to other programs. EPISOL is compatible with the force-field and coordinate files from both AMBER and GROMACS simulation packages. Moreover, EPISOL is equipped with an internal memory control to efficiently manage the use of physical memory, making it suitable for performing calculations on large biomolecules. We demonstrate that EPISOL can efficiently and accurately calculate solvation density distributions around various solute molecules (including a protein chaperone consisting of 120,715 atoms) and obtain solvent free energy for a wide range of organic compounds. We expect that EPISOL can be widely applied as a solvation model for chemical and biological systems. EPISOL is available at https://github.com/EPISOLrelease/EPISOL.

Development and test of highly accurate endpoint free energy methods. 1: Evaluation of ABCG2 charge model on solvation free energy prediction and optimization of atom radii suitable for more accurate solvation free energy prediction by the PBSA method

Accurate estimation of solvation free energy (SFE) lays the foundation for accurate prediction of binding free energy. The Poisson-Boltzmann (PB) or generalized Born (GB) combined with surface area (SA) continuum solvation method (PBSA and GBSA) have been widely used in SFE calculations because they can achieve good balance between accuracy and efficiency. However, the accuracy of these methods can be affected by several factors such as the charge models, polar and nonpolar SFE calculation methods and the atom radii used in the calculation. In this work, the performance of the ABCG2 (AM1-BCC-GAFF2) charge model as well as other two charge models, that is, RESP (Restrained Electrostatic Potential) and AM1-BCC (Austin Model 1-bond charge corrections), on the SFE prediction of 544 small molecules in water by PBSA/GBSA was evaluated. In order to improve the performance of the PBSA prediction based on the ABCG2 charge, we further explored the influence of atom radii on the prediction accuracy and yielded a set of atom radius parameters for more accurate SFE prediction using PBSA based on the ABCG2/GAFF2 by reproducing the thermodynamic integration (TI) calculation results. The PB radius parameters of carbon, oxygen, sulfur, phosphorus, chloride, bromide and iodine, were adjusted. New atom types, on, oi, hn1, hn2, hn3, were introduced to further improve the fitting performance. Then, we tuned the parameters in the nonpolar SFE model using the experimental SFE data and the PB calculation results. By adopting the new radius parameters and new nonpolar SFE model, the root mean square error (RMSE) of the SFE calculation for the 544 molecules decreased from 2.38 to 1.05 kcal/mol. Finally, the new radius parameters were applied in the prediction of protein-ligand binding free energies using the MM-PBSA method. For the eight systems tested, we could observe higher correlation between the experiment data and calculation results and smaller prediction errors for the absolute binding free energies, demonstrating that our new radius parameters can improve the free energy calculation using the MM-PBSA method.

Assessment of Different Parameters on the Accuracy of Computational Alanine Scanning of Protein-Protein Complexes with the Molecular Mechanics/Generalized Born Surface Area Method

Computational alanine scanning with the molecular mechanics generalized Born surface area (MM/GBSA) method constitutes a widely used approach for identifying critical residues at protein-protein interfaces. Despite its popularity, the MM/GBSA method still has certain drawbacks due to its dependence on many factors. Here, we performed a systematical study on the impact of four different parameters, namely, the internal dielectric constant, the generalized Born model, the entropic term, and the inclusion of structural waters on the accuracy of computational alanine scanning calculations with the MM/GBSA method. Our results show that the internal dielectric constant is the most critical parameter for getting accurate predictions. The introduction of entropy and interfacial water molecules decreased the quality of the predictions, while the generalized Born model had little to no effect. Considering the significance of the internal dielectric value, we proposed a methodology based on the energetic predominance of a particular set of amino acids at the protein-protein interface for selecting an appropriate value for this variable. We hope that these results serve as a guideline for future studies of protein-protein complexes using the MM/GBSA method.

A combined ligand and target-based virtual screening strategy to repurpose drugs as putrescine uptake inhibitors with trypanocidal activity

Chagas disease, also known as American trypanosomiasis, is a neglected tropical disease caused by the protozoa Trypanosoma cruzi, affecting nearly 7 million people only in the Americas. Polyamines are essential compounds for parasite growth, survival, and differentiation. However, because trypanosomatids are auxotrophic for polyamines, they must be obtained from the host by specific transporters. In this investigation, an ensemble of QSAR classifiers able to identify polyamine analogs with trypanocidal activity was developed. Then, a multi-template homology model of the dimeric polyamine transporter of T. cruzi, TcPAT12, was created with Rosetta, and then refined by enhanced sampling molecular dynamics simulations. Using representative snapshots extracted from the trajectory, a docking model able to discriminate between active and inactive compounds was developed and validated. Both models were applied in a parallel virtual screening campaign to repurpose known drugs as anti-trypanosomal compounds inhibiting polyamine transport in T. cruzi. Montelukast, Quinestrol, Danazol, and Dutasteride were selected for in vitro testing, and all of them inhibited putrescine uptake in biochemical assays, confirming the predictive ability of the computational models. Furthermore, all the confirmed hits proved to inhibit epimastigote proliferation, and Quinestrol and Danazol were able to inhibit, in the low micromolar range, the viability of trypomastigotes and the intracellular growth of amastigotes.

Discovery of All-d-Peptide Inhibitors of SARS-CoV-2 3C-like Protease

During the replication process of SARS-CoV-2, the main protease of the virus [3-chymotrypsin-like protease (3CL^pro)] plays a pivotal role and is essential for the life cycle of the pathogen. Numerous studies have been conducted so far, which have confirmed 3CL^pro as an attractive drug target to combat COVID-19. We describe a novel and efficient next-generation sequencing (NGS) supported phage display selection strategy for the identification of a set of SARS-CoV-2 3CL^pro targeting peptide ligands that inhibit the 3CL protease, in a competitive or noncompetitive mode, in the low μM range. From the most efficient l-peptides obtained from the phage display, we designed all-d-peptides based on the retro-inverso (ri) principle. They had IC₅₀ values also in the low μM range and in combination, even in the sub-micromolar range. Additionally, the combination with Rutinprivir decreases 10-fold the IC₅₀ value of the competitive inhibitor. The inhibition modes of these d-ri peptides were the same as their respective l-peptide versions. Our results demonstrate that retro-inverso obtained all-d-peptides interact with high affinity and inhibit the SARS-CoV-2 3CL protease, thus reinforcing their potential for further development toward therapeutic agents. The here described d-ri peptides address limitations associated with current l-peptide inhibitors and are promising lead compounds. Further optimization regarding pharmacokinetic properties will allow the development of even more potent d-peptides to be used for the prevention and treatment of COVID-19.

tLyp-1: A peptide suitable to target NRP-1 receptor

Targeting vascular endothelial growth factor receptor (VEFGR) and its co-receptor neuropilin-1 (NRP-1) is an interesting vascular strategy. tLyp-1 is a tumor-homing and penetrating peptide of 7 amino acids (CGNKRTR). It is a truncated form of Lyp-1 (CGNKRTRGC), which is known to target NRP-1 receptor, with high affinity and specificity. It is mediated by endocytosis via C-end rule (CendR) internalization pathway. The aim of this study is to evaluate the importance of each amino acid in the tLyp-1 sequence through alanine-scanning (Ala-scan) technique, during which each of the amino acid in the sequence was systematically replaced by alanine to produce 7 different analogues. In silico approach through molecular docking and molecular dynamics are employed to understand the interaction between the peptide and its analogues with the NRP-1 receptor, followed by in vitro ligand binding assay study. The C-terminal Arg is crucial in the interaction of tLyp-1 with NRP-1 receptor. Substituting this residue dramatically reduces the affinity of this peptide which is clearly seen in this study. Lys-4 is also important in the interaction, which is confirmed via the in vitro study and the MM-PBSA analysis. The finding in this study supports the CendR, in which the presence of R/K-XX-R/K motif is essential in the binding of a ligand with NRP-1 receptor. This presented work will serve as a guide in the future work pertaining the development of active targeting agent towards NRP-1 receptor.

MM/PB(GB)SA benchmarks on soluble proteins and membrane proteins

Predicting protein-ligand binding free energy rapidly and accurately remains a challenging question in modern drug discovery. Molecular mechanics/Poisson-Boltzmann (Generalized Born) surface area (MM/PB(GB)SA) has emerged as an essential tool for accelerating cost-efficient binding free energy calculation. This study presents benchmarks with three membrane-bound protein systems and six soluble protein systems. Different parameters were sampled for different benchmarks to explore the highest accuracy. These include ligand charges, protein force fields, extra points, GB models, nonpolar optimization methods, internal dielectric constants and membrane dielectric constants. Comparisons of accuracy were made between MM/PB(GB)SA, docking and free energy perturbation (FEP). The results reveal a competitive performance between MM/PB(GB)SA and FEP. In summary, MM/PB(GB)SA is a powerful approach to predict ligand binding free energy rapidly and accurately. Parameters of MM/PB(GB)SA calculations, such as the GB models and membrane dielectric constants, need to be optimized for different systems. This method can be served as a powerful tool for drug design.

Electrostatic Contributions to the Binding Free Energy of Nicotine to the Acetylcholine Binding Protein

Biomolecular binding relies on specific attractive interactions between two partner molecules, including electrostatics, dispersion, hydrophobicity, and solvation. Assessing the contributions of electrostatic interactions to binding is key to the understanding of ligand binding mechanisms and the design of improved biomolecular binders. For example, nicotine is a well-known agonist of nicotinic acetylcholine receptors (nAChRs), but the molecular mechanisms for the differential action of nicotine on brain and muscle nAChRs remain elusive. In this work, we have chosen the acetylcholine binding protein (AChBP) in complex with nicotine as a model system to interrogate the electrostatic contributions to nicotine binding. Our absolute binding free energy simulations confirm that nicotine binds AChBP predominantly in its protonated (charged) form. By comparing energetic contributions from decomposed interactions for either neutral or charged nicotine, our calculations shed light on the nature of the binding of nicotine to the AChBP. The preferred binding of charged nicotine over neutral nicotine originates from its stronger electrostatic interactions with AChBP, a cation-π interaction to a tryptophan residue and a hydrogen bond between nicotine and the backbone carbonyl of the tryptophan, whereas the major force driving the binding process appears to be van der Waals interactions. The various nonelectrostatic terms can also indirectly modulate the electrostatic interactions through fine-tuning the binding pose of the ligand in the binding site, providing an explanation of why the binding specificity of nicotine to the brain versus muscle nAChRs is driven by electrostatic interaction, given that the immediate binding site residues, including the key tryptophan residue, are identical in the two receptors.

Exploring the druggability of the binding site of aurovertin, an exogenous allosteric inhibitor of FOF1-ATP synthase

In addition to playing a central role in the mitochondria as the main producer of ATP, F_OF₁-ATP synthase performs diverse key regulatory functions in the cell membrane. Its malfunction has been linked to a growing number of human diseases, including hypertension, atherosclerosis, cancer, and some neurodegenerative, autoimmune, and aging diseases. Furthermore, inhibition of this enzyme jeopardizes the survival of several bacterial pathogens of public health concern. Therefore, F_OF₁-ATP synthase has emerged as a novel drug target both to treat human diseases and to combat antibiotic resistance. In this work, we carried out a computational characterization of the binding sites of the fungal antibiotic aurovertin in the bovine F₁ subcomplex, which shares a large identity with the human enzyme. Molecular dynamics simulations showed that although the binding sites can be described as preformed, the inhibitor hinders inter-subunit communications and exerts long-range effects on the dynamics of the catalytic site residues. End-point binding free energy calculations revealed hot spot residues for aurovertin recognition. These residues were also relevant to stabilize solvent sites determined from mixed-solvent molecular dynamics, which mimic the interaction between aurovertin and the enzyme, and could be used as pharmacophore constraints in virtual screening campaigns. To explore the possibility of finding species-specific inhibitors targeting the aurovertin binding site, we performed free energy calculations for two bacterial enzymes with experimentally solved 3D structures. Finally, an analysis of bacterial sequences was carried out to determine conservation of the aurovertin binding site. Taken together, our results constitute a first step in paving the way for structure-based development of new allosteric drugs targeting F_OF₁-ATP synthase sites of exogenous inhibitors.

Identification of Novel Ribonucleotide Reductase Inhibitors for Therapeutic Application in Bile Tract Cancer: An Advanced Pharmacoinformatics Study

Biliary tract cancer (BTC) is constituted by a heterogeneous group of malignant tumors that may develop in the biliary tract, and it is the second most common liver cancer. Human ribonucleotide reductase M1 (hRRM1) has already been proven to be a potential BTC target. In the current study, a de novo design approach was used to generate novel and effective chemical therapeutics for BTC. A set of comprehensive pharmacoinformatics approaches was implemented and, finally, seventeen potential molecules were found to be effective for the modulation of hRRM1 activity. Molecular docking, negative image-based ShaEP scoring, absolute binding free energy, in silico pharmacokinetics, and toxicity assessments corroborated the potentiality of the selected molecules. Almost all molecules showed higher affinity in comparison to gemcitabine and naphthyl salicylic acyl hydrazone (NSAH). On binding interaction analysis, a number of critical amino acids was found to hold the molecules at the active site cavity. The molecular dynamics (MD) simulation study also indicated the stability between protein and ligands. High negative MM-GBSA (molecular mechanics generalized Born and surface area) binding free energy indicated the potentiality of the molecules. Therefore, the proposed molecules might have the potential to be effective therapeutics for the management of BTC.

Extracellular pH sensing by plant cell-surface peptide-receptor complexes

The extracellular pH is a vital regulator of various biological processes in plants. However, how plants perceive extracellular pH remains obscure. Here, we report that plant cell-surface peptide-receptor complexes can function as extracellular pH sensors. We found that pattern-triggered immunity (PTI) dramatically alkalinizes the acidic extracellular pH in root apical meristem (RAM) region, which is essential for root meristem growth factor 1 (RGF1)-mediated RAM growth. The extracellular alkalinization progressively inhibits the acidic-dependent interaction between RGF1 and its receptors (RGFRs) through the pH sensor sulfotyrosine. Conversely, extracellular alkalinization promotes the alkaline-dependent binding of plant elicitor peptides (Peps) to its receptors (PEPRs) through the pH sensor Glu/Asp, thereby promoting immunity. A domain swap between RGFR and PEPR switches the pH dependency of RAM growth. Thus, our results reveal a mechanism of extracellular pH sensing by plant peptide-receptor complexes and provide insights into the extracellular pH-mediated regulation of growth and immunity in the RAM.

Solvation Thermodynamics of Solutes in Water and Ionic Liquids Using the Multiscale Solvation-Layer Interface Condition Continuum Model

Molecular assembly processes are generally driven by thermodynamic properties in solutions. Atomistic modeling can be very helpful in designing and understanding complex systems, except that bulk solvent is very inefficient to treat explicitly as discrete molecules. In this work, we develop and assess two multiscale solvation models for computing solvation thermodynamic properties. The new SLIC/CDC model combines continuum solvent electrostatics based on the solvent layer interface condition (SLIC) with new statistical thermodynamic models for hydrogen bonding and nonpolar modes: cavity formation, dispersion interactions, combinatorial mixing (CDC). Given the structures of 500 solutes, the SLIC/CDC model predicts Gibbs energies of solvation in water with an average accuracy better than 1 kcal/mol, when compared to experimental measurements, and better than 0.8 kcal/mol, when compared to explicit-solvent molecular dynamics simulations. The individual SLIC/CDC energy mode values agree quantitatively with those computed from explicit-solvent molecular dynamics. The previously published SLIC/SASA multiscale model combines the SLIC continuum electrostatic model with the solvent-accessible surface area (SASA) nonpolar energy mode. With our new, improved parametrization method, the SLIC/SASA model now predicts Gibbs energies of solvation with better than 1.4 kcal/mol average accuracy in aqueous systems, compared to experimental and explicit-solvent molecular dynamics, and better than 1.6 kcal/mol average accuracy in ionic liquids, compared to explicit-solvent molecular dynamics. Both models predict solvation entropies, and are the first implicit-solvation models capable of predicting solvation heat capacities.

Identification of Potential Cytochrome P450 3A5 Inhibitors: An Extensive Virtual Screening through Molecular Docking, Negative Image-Based Screening, Machine Learning and Molecular Dynamics Simulation Studies

Cytochrome P450 3A5 (CYP3A5) is one of the crucial CYP family members and has already proven to be an important drug target for cardiovascular diseases. In the current study, the PubChem database was screened through molecular docking and high-affinity molecules were adopted for further assessment. A negative image-based (NIB) model was used for a similarity search by considering the complementary shape and electrostatics of the target and small molecules. Further, the molecules were segregated into active and inactive groups through six machine learning (ML) matrices. The active molecules found in each ML model were used for in silico pharmacokinetics and toxicity assessments. A total of five molecules followed the acceptable pharmacokinetics and toxicity profiles. Several potential binding interactions between the proposed molecules and CYP3A5 were observed. The dynamic behavior of the selected molecules in the CYP3A5 was explored through a molecular dynamics (MD) simulation study. Several parameters obtained from the MD simulation trajectory explained the stability of the protein–ligand complexes in dynamic states. The high binding affinity of each molecule was revealed by the binding free energy calculation through the MM-GBSA methods. Therefore, it can be concluded that the proposed molecules might be potential CYP3A5 molecules for therapeutic application in cardiovascular diseases subjected to in vitro/in vivo validations.

A strategy to optimize the peptide-based inhibitors against different mutants of the spike protein of SARS-CoV-2

The SARS-CoV-2 virus has caused high-priority health concerns at a global level. Vaccines have stalled the proliferation of viruses to some extent. Yet, the emergence of newer, potentially more infectious, and dangerous mutants such as Delta and Omicron are among the major challenges in finding a more permanent solution for this pandemic. The effectiveness of antivirals Molnupiravir and Paxlovid, authorized for emergency use by the FDA, are yet to be assessed on a larger population. Patients with a high risk of disease progression have received treatment with antibody-cocktail. Most of the mutations leading to the new lineage of SARS-CoV-2 are found in the spike protein of this virus that plays a key role in facilitating host entry. The current study has investigated how to modify a promising peptide-based inhibitor of spike protein, LCB3, against common mutations, N501Y and K417N in the target protein so that it retains its efficacy against the spike protein. LCB3 being a prototype for protein-based inhibitors is an ideal testing system to learn about protein-based inhibitors. This study proposes the substitutions of amino acid residues of LCB3 inhibitor using a structure-based approach that considers free energy decomposition of residues, the distance between atoms, and the interaction among amino acids. The binding free energy calculations suggest a possible improvement in the binding affinity of existing inhibitor LCB3 to the mutant forms of the S-protein using simple substitutions at specific positions of the inhibitor. This approach, being general, can be used in different inhibitors and other mutations and help in fighting against SARS-CoV-2.Communicated by Ramaswamy H. Sarma.

Naturally Occurring 8ß,13ß-kaur-15-en-17-al and Anti-Malarial Activity from Podocarpus polystachyus Leaves

Despite much interest and studies toward the genus Podocarpus, the anti-malarial evaluation of Podocarpus polystachyus’s phytoconstituents remains lacking. Herein, the phytoconstituents of P. polystachyus leaves and their anti-malarial effect against Plasmodium falciparum were investigated for the first time. One new natural product, 8ß,13ß-kaur-15-en-17-al (1), along with three known compounds, 8ß,13ß-kaur-15-en-17-ol (2) and 13ß-kaur-16-ene (3), and α-tocopherol hydroquinone (4) were isolated via HR-ESI-MS and NMR analyses. Compounds 1 and 2 inhibited P. falciparum growth at 12 and 52 µM of IC₅₀, respectively. Their anti-malarial activity was associated with the in silico P. falciparum lactate dehydrogenase (PfLDH) inhibition. Molecular docking of ligands 1 and 2 with the putative target PfLDH revealed ~−2 kcal/mol of binding energies more negative than the control. Molecular dynamic simulations (100 ns) showed equal or smaller deviation values (RMSD, RMSF, Rg) and stronger interactions of PfLDH-1 and PfLDH-2 complexes via at least one consistent H-bond than the control. Additionally, a slightly increased PfLDH H-bond profile in their interactions improved the PfLDH dynamic and structural stabilities. Overall, this study supports the relevance of 1 and 2 as plasmodial growth inhibitors with their putative anti-PfLDH activity, which could be a potential scaffold for developing anti-malarial drugs.

Calcium Ion Binding to the Mutants of Calmodulin: A Structure-Based Computational Predictive Model of Binding Affinity Using a Charge Scaling Approach in Molecular Dynamics Simulation

The binding of calcium ions (Ca²⁺) to the calcium-binding proteins (CBPs) controls a plethora of regulatory processes. Among the roles played by CBPs in several diseases, the onset and progress of some cardiovascular diseases are caused by mutations in calmodulin (CaM), an important member of CBPs. Rationalization and prediction of the binding affinity of Ca²⁺ ions to the CaM can play important roles in understanding the origin of cardiovascular diseases. However, there is no robust structure-based computational method for predicting the binding affinity of Ca²⁺ ions to the different forms of CBPs in general and CaM in particular. In the current work, we have devised a fast yet accurate computational technique to accurately calculate the binding affinity of Ca²⁺ to the different forms of CaM. This method combines the well-known molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method and a charge scaling approach developed by previous authors that takes care of the polarization of CaM and Ca²⁺ ions. Our detailed analysis of the different components of binding free energy shows that subtle changes in electrostatics and van der Waals contribute to the difference in the binding affinity of mutants from that of the wild type (WT), and the charge scaling approach is superior in calculating these subtle changes in electrostatics as compared to the nonpolarizable force field used in this work. A statistically significant regression model made from our binding free energy calculations gives a correlation coefficient close to 0.8 to the experimental results. This structure-based predictive model can open up a new strategy to understand and predict the binding of Ca²⁺ to the mutants of CBPs, in general.

Computational design of stapled peptide inhibitor against SARS‐CoV ‐2 receptor binding domain

Since its first detection in 2019, the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS‐CoV‐2) has been the cause of millions of deaths worldwide. Despite the development and administration of different vaccines, the situation is still worrisome as the virus is constantly mutating to produce newer variants some of which are highly infectious. This raises an urgent requirement to understand the infection mechanism and thereby design therapeutic‐based treatment for COVID‐19. The gateway of the virus to the host cell is mediated by the binding of the receptor binding domain (RBD) of the virus spike protein to the angiotensin‐converting enzyme 2 (ACE2) of the human cell. Therefore, the RBD of SARS‐CoV‐2 can be used as a target to design therapeutics. The α1 helix of ACE2, which forms direct contact with the RBD surface, has been used as a template in the current study to design stapled peptide therapeutics. Using computer simulation, the mechanism and thermodynamics of the binding of six stapled peptides with RBD have been estimated. Among these, the one with two lactam stapling agents has shown binding affinity, sufficient to overcome RBD‐ACE2 binding. Analyses of the mechanistic detail reveal that a reorganization of amino acids at the RBD‐ACE2 interface produces favorable enthalpy of binding whereas conformational restriction of the free peptide reduces the loss in entropy to result higher binding affinity. The understanding of the relation of the nature of the stapling agent with their binding affinity opens up the avenue to explore stapled peptides as therapeutic against SARS‐CoV‐2.

QM/MM Energy Decomposition Using the Interacting Quantum Atoms Approach

The interacting quantum atoms (IQA) method decomposes the quantum mechanical (QM) energy of a molecular system in terms of one- and two-center (atomic) contributions within the context of the quantum theory of atoms in molecules. Here, we demonstrate that IQA, enhanced with molecular mechanics (MM) and Poisson-Boltzmann surface-area (PBSA) solvation methods, is naturally extended to the realm of hybrid QM/MM methodologies, yielding intra- and inter-residue energy terms that characterize all kinds of covalent and noncovalent bonding interactions. To test the robustness of this approach, both metal-water interactions and QM/MM boundary artifacts are characterized in terms of the IQA descriptors derived from QM regions of varying size in Zn(II)- and Mg(II)-water clusters. In addition, we analyze a homologous series of inhibitors in complex with a matrix metalloproteinase (MMP-12) by carrying out QM/MM-PBSA calculations on their crystallographic structures followed by IQA energy decomposition. Overall, these applications not only show the advantages of the IQA QM/MM approach but also address some of the challenges lying ahead for expanding the QM/MM methodology.

Computational investigation of the impact of core sequence on immobile DNA four-way junction structure and dynamics

Immobile four-way junctions (4WJs) are core structural motifs employed in the design of programmed DNA assemblies. Understanding the impact of sequence on their equilibrium structure and flexibility is important to informing the design of complex DNA architectures. While core junction sequence is known to impact the preferences for the two possible isomeric states that junctions reside in, previous investigations have not quantified these preferences based on molecular-level interactions. Here, we use all-atom molecular dynamics simulations to investigate base-pair level structure and dynamics of four-way junctions, using the canonical Seeman J1 junction as a reference. Comparison of J1 with equivalent single-crossover topologies and isolated nicked duplexes reveal conformational impact of the double-crossover motif. We additionally contrast J1 with a second junction core sequence termed J24, with equal thermodynamic preference for each isomeric configuration. Analyses of the base-pair degrees of freedom for each system, free energy calculations, and reduced-coordinate sampling of the 4WJ isomers reveal the significant impact base sequence has on local structure, isomer bias, and global junction dynamics.

gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS

Molecular mechanics/Poisson-Boltzmann (Generalized-Born) surface area is one of the most popular methods to estimate binding free energies. This method has been proven to balance accuracy and computational efficiency, especially when dealing with large systems. As a result of its popularity, several programs have been developed for performing MM/PB(GB)SA calculations within the GROMACS community. These programs, however, present several limitations. Here we present gmx_MMPBSA, a new tool to perform end-state free energy calculations from GROMACS molecular dynamics trajectories. gmx_MMPBSA provides the user with several options, including binding free energy calculations with different solvation models (PB, GB, or 3D-RISM), stability calculations, computational alanine scanning, entropy corrections, and binding free energy decomposition. Noteworthy, several promising methodologies to calculate relative binding free energies such as alanine scanning with variable dielectric constant and interaction entropy have also been implemented in gmx_MMPBSA. Two additional tools-gmx_MMPBSA_test and gmx_MMPBSA_ana-have been integrated within gmx_MMPBSA to improve its usability. Multiple illustrating examples can be accessed through gmx_MMPBSA_test, while gmx_MMPBSA_ana provides fast, easy, and efficient access to different graphics plotted from gmx_MMPBSA output files. The latest version (v1.4.3, 26/05/2021) is available free of charge (documentation, test files, and tutorials included) at https://github.com/Valdes-Tresanco-MS/gmx_MMPBSA.

Generalizing Continuum Solvation in Crystal to Nonaqueous Solvents: Implementation, Parametrization, and Application to Molecules and Surfaces

We present an extension of a generalized finite-difference Poisson-Boltzmann (FDPB) continuum solvation model based on a self-consistent reaction field treatment to nonaqueous solvents. Implementation and reparametrization of the cavitation, dispersion, and structural (CDS) effects nonelectrostatic model are presented in CRYSTAL, with applications to both finite and infinite periodic systems. For neutral finite systems, computed errors with respect to available experimental data on free energies of solvation of 2523 solutes in 91 solvents, as well as 144 transfer energies from water to 14 organic solvents are on par with the reference SM12 solvation model for which the CDS parameters have been developed. Calculations performed on a TiO₂ anatase surface and compared to VASPsol data revealed an overall very good agreement of computed solvation energies, surface energies, as well as band structure changes upon solvation in three different solvents, validating the general applicability of the reparametrized FDPB approach to neutral nonperiodic and periodic solutes in aqueous and nonaqueous solvents. For ionic species, while the reparametrized CDS model led to large errors on free energies of solvation of anions, addition of a corrective term based on Abraham's acidity of the solvent significantly improved the accuracy of the proposed continuum solvation model, leading to errors on aqueous pK_a of a test set of 83 solutes divided by a factor of 4 compared to the reference solvation model based on density (SMD). Overall, therefore, these encouraging results demonstrate that the generalized FDPB continuum solvation model can be applied to a broad range of solutes in various solvents, ranging from finite neutral or charged solutes to extended periodic surfaces.

Design and Characterization of a Cell-Penetrating Peptide Derived from the SOX2 Transcription Factor

SOX2 is an oncogenic transcription factor overexpressed in nearly half of the basal-like triple-negative breast cancers associated with very poor outcomes. Targeting and inhibiting SOX2 is clinically relevant as high SOX2 mRNA levels are positively correlated with decreased overall survival and progression-free survival in patients affected with breast cancer. Given its key role as a master regulator of cell proliferation, SOX2 represents an important scaffold for the engineering of dominant-negative synthetic DNA-binding domains (DBDs) that act by blocking or interfering with the oncogenic activity of the endogenous transcription factor in cancer cells. We have synthesized an interference peptide (iPep) encompassing a truncated 24 amino acid long C-terminus of SOX2 containing a potential SOX-specific nuclear localization sequence, and the determinants of the binding of SOX2 to the DNA and to its transcription factor binding partners. We found that the resulting peptide (SOX2-iPep) possessed intrinsic cell penetration and promising nuclear localization into breast cancer cells, and decreased cellular proliferation of SOX2 overexpressing cell lines. The novel SOX2-iPep was found to exhibit a random coil conformation predominantly in solution. Molecular dynamics simulations were used to characterize the interactions of both the SOX2 transcription factor and the SOX2-iPep with FGF4-enhancer DNA in the presence of the POU domain of the partner transcription factor OCT4. Predictions of the free energy of binding revealed that the iPep largely retained the binding affinity for DNA of parental SOX2. This work will enable the future engineering of novel dominant interference peptides to transport different therapeutic cargo molecules such as anti-cancer drugs into cells.

Recent Developments in Free Energy Calculations for Drug Discovery

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Optimized Halogen Atomic Radii for PBSA Calculations Using Off-Center Point Charges

In force-field methods, the usage of off-center point charges, also called extra points (EPs), is a common strategy to tackle the anisotropy of the electrostatic potential of covalently bonded halogens (X), thus allowing the description of halogen bonds (XBs) at the molecular mechanics/molecular dynamics (MM/MD) level. Diverse EP implementations exist in the literature differing on the charge sets and/or the X-EP distances. Poisson-Boltzmann and surface area (PBSA) calculations can be used to obtain solvation free energies (ΔG_solv) of small molecules, often to compute binding free energies (ΔG_bind) at the MM-PBSA level. This method depends, among other parameters, on the empirical assignment of atomic radii (PB radii). Given the multiplicity of off-center point-charge models and the lack of specific PB radii for halogens compatible with such implementations, in this work, we assessed the performance of PBSA calculations for the estimation of ΔG_solv values in water (ΔG_hyd), also conducting an optimization of the halogen PB radii (Cl, Br, and I) for each EP model. We not only expand the usage of EP models in the scope of the general AMBER force field (GAFF) but also provide the first optimized halogen PB radii in the context of the CHARMM general force field (CGenFF), thus contributing to improving the description of halogenated compounds in PBSA calculations.

Computational design of ultrashort peptide inhibitors of the receptor-binding domain of the SARS-CoV-2 S protein

Targeting the interaction between severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2)-receptor-binding domain (RBD) and angiotensin-converting enzyme 2 (ACE2) is believed to be an effective strategy for drug design to inhibit the infection of SARS-CoV-2. Herein, several ultrashort peptidase inhibitors against the RBD-ACE2 interaction were obtained by a computer-aided approach based on the RBD-binding residues on the protease domain (PD) of ACE2. The designed peptides were tested on a model coronavirus GX_P2V, which has 92.2 and 86% amino acid identity to the SARS-CoV-2 spike protein and RBD, respectively. Molecular dynamics simulations and binding free energy analysis predicted a potential binding pocket on the RBD of the spike protein, and this was confirmed by the specifically designed peptides SI5α and SI5α-b. They have only seven residues, showing potent antiviral activity and low cytotoxicity. Enzyme-linked immunosorbent assay result also confirmed their inhibitory ability against the RBD-ACE2 interaction. The ultrashort peptides are promising precursor molecules for the drug development of Corona Virus Disease 2019, and the novel binding pocket on the RBD may be helpful for the design of RBD inhibitors or antibodies against SARS-CoV-2.

Cationic Dendrimers for siRNA Delivery: Computational Approaches for Characterization

Nowadays, computer simulations have been established as a fundamental tool in the design and development of new dendrimer-based nanocarriers for drug and gene delivery. Moreover, the level of detail contained in the information that can be gathered by performing atomistic-scale simulations cannot be obtained with any other available experimental technique. In this chapter we describe the main computational toolbox that can be exploited in the different stages of novel dendritic nanocarrier production-from the initial conception to the stage of biological intermolecular interactions.

Discovery of Selective Inhibitor Leads by Targeting an Allosteric Site in Insulin-Regulated Aminopeptidase

Insulin-Regulated aminopeptidase (IRAP) is a zinc-dependent aminopeptidase with several important biological functions and is an emerging pharmaceutical target for cognitive enhancement and immune system regulation. Aiming to discover lead-like IRAP inhibitors with enhanced selectivity versus homologous enzymes, we targeted an allosteric site at the C-terminal domain pocket of IRAP. We compiled a library of 2.5 million commercially available compounds from the ZINC database, and performed molecular docking at the target pocket of IRAP and the corresponding pocket of the homologous endoplasmic reticulum aminopeptidase 1 (ERAP1). Of the top compounds that showed high selectivity, 305 were further analyzed by molecular dynamics simulations and free energy calculations, leading to the selection of 33 compounds for in vitro evaluation. Two orthogonal functional assays were employed: one using a small fluorogenic substrate and one following the degradation of oxytocin, a natural peptidic substrate of IRAP. In vitro evaluation suggested that several of the compounds tested can inhibit IRAP, but the inhibition profile was dependent on substrate size, consistent with the allosteric nature of the targeted site. Overall, our results describe several novel leads as IRAP inhibitors and suggest that the C-terminal domain pocket of IRAP is a promising target for developing highly selective IRAP inhibitors.

Implicit Solvents for the Polarizable Atomic Multipole AMOEBA Force Field

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the nonlinear and linearized versions of the Poisson-Boltzmann equation (PBE), the domain-decomposition conductor-like screening model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatics models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least-squares style target function based on a library of 103 small-molecule solvation free energy differences. Mean signed errors for the adaptive Poisson-Boltzmann solver (APBS), ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.58 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields opens the door to their use for folding and design applications.

Estimating the Roles of Protonation and Electronic Polarization in Absolute Binding Affinity Simulations

Accurate prediction of binding free energies is critical to streamlining the drug development and protein design process. With the advent of GPU acceleration, absolute alchemical methods, which simulate the removal of ligand electrostatics and van der Waals interactions with the protein, have become routinely accessible and provide a physically rigorous approach that enables full consideration of flexibility and solvent interaction. However, standard explicit solvent simulations are unable to model protonation or electronic polarization changes upon ligand transfer from water to the protein interior, leading to inaccurate prediction of binding affinities for charged molecules. Here, we perform extensive simulation totaling ∼540 μs to benchmark the impact of modeling conditions on predictive accuracy for absolute alchemical simulations. Binding to urokinase plasminogen activator (UPA), a protein frequently overexpressed in metastatic tumors, is evaluated for a set of 10 inhibitors with extended flexibility, highly charged character, and titratable properties. We demonstrate that the alchemical simulations can be adapted to utilize the MBAR/PBSA method to improve the accuracy upon incorporating electronic polarization, highlighting the importance of polarization in alchemical simulations of binding affinities. Comparison of binding energy prediction at various protonation states indicates that proper electrostatic setup is also crucial in binding affinity prediction of charged systems, prompting us to propose an alternative binding mode with protonated ligand phenol and Hid-46 at the binding site, a testable hypothesis for future experimental validation.

Evaluating the Multifunctionality of a New Modulator of Zinc-Induced Aβ Aggregation Using a Novel Computational Approach

The high concentration of zinc metal ions in Aβ aggregations is one of the most cited hallmarks of Alzheimer's disease (AD), and several substantial pieces of evidence emphasize the key role of zinc metal ions in the pathogenesis of AD. In this study, while designing a multifunctional peptide for simultaneous targeting Aβ aggregation and chelating the zinc metal ion, a novel and comprehensive approach is introduced for evaluating the multifunctionality of a multifunctional drugs based on computational methods. The multifunctional peptide consists of inhibitor and chelator domains, which are included in the C-terminal hydrophobic region of Aβ, and the first four amino acids of human albumin. The ability of the multifunctional peptide in zinc ion chelation has been investigated using molecular dynamics (MD) simulations of the peptide-zinc interaction for 300 ns, and Bennett's acceptance ratio (BAR) method has been used to accurately calculate the chelation free energy. Data analysis demonstrates that the peptide chelating domain can be stably linked to the zinc ion. Besides, the introduced method used for evaluating chelation and calculating the free energy of peptide binding to zinc ions was successfully validated by comparison with previous experimental and theoretical published data. The results indicate that the multifunctional peptide, coordinating with the zinc metal ion, can be effective in Aβ inhibition by preserving the native helical structure of the Aβ₄₂ monomer as well as disrupting the β-sheet structure of Aβ₄₂ aggregates. Detailed assessments of the Aβ₄₂-peptide interactions elucidate that the inhibition of Aβ is achieved by considerable hydrophobic interactions and hydrogen bonding between the multifunctional peptide and the hydrophobic Aβ regions, along with interfering in stable bridges formed inside the Aβ aggregate.

How Parkinson's disease-related mutations disrupt the dimerization of WD40 domain in LRRK2: a comparative molecular dynamics simulation study

The multidomain kinase enzyme leucine-rich-repeat kinase 2 (LRRK2), activated through a homodimerization manner, has been identified as an important pathogenic factor in Parkinson's disease (PD), the second most common neurodegenerative disease wordwide. The Trp-Asp-40 (WD40) domain, located in the C-terminal LRRK2, harbours one of the most frequent PD-related variants, G2385R. However, the detailed dynamics of WD40 during LRRK2 dimerization and the underlying mechanism through which the pathogenic mutations disrupt the formation of the WD40 dimer have remained elusive. Here, microsecond-scale molecular dynamics simulations were employed to provide a mechanistic view underlying the WD40 dimerization and unveil the structural basis by which the interface-based mutations G2385R, H2391D and R2394E compromise the corresponding process. The simulation results identified important residues, D2351, R2394, E2395, R2413, and R2443, involved in establishing the complex binding network along the dimerization interface, which was significantly weakened in the presence of interfacial mutations. A "sag-bulge" model was proposed to explain the unfavorable dimer formation in the mutant systems. In addition, mutations altered the community configuration in the wild-type system in which inter-monomeric interplay is prominent, leading to the destabilization of the WD40 dimer under mutation.

Protein–ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange–correlation functional

Quantum mechanical binding free energies based on thousands of full-protein DFT calculations are tractable, reproducible and converge well.

Cell Type-Dependent Escape of Capsid Inhibitors by Simian Immunodeficiency Virus SIVcpz

Pandemic human immunodeficiency virus type 1 (HIV-1) is the result of the zoonotic transmission of simian immunodeficiency virus (SIV) from the chimpanzee subspecies Pan troglodytes troglodytes (SIVcpzPtt). The related subspecies Pan troglodytes schweinfurthii is the host of a similar virus, SIVcpzPts, which did not spread to humans. We tested these viruses with small-molecule capsid inhibitors (PF57, PF74, and GS-CA1) that interact with a binding groove in the capsid that is also used by CPSF6. While HIV-1 was sensitive to capsid inhibitors in cell lines, human macrophages, and peripheral blood mononuclear cells (PBMCs), SIVcpzPtt was resistant in rhesus FRhL-2 cells and human PBMCs but was sensitive to PF74 in human HOS and HeLa cells. SIVcpzPts was insensitive to PF74 in FRhL-2 cells, HeLa cells, PBMCs, and macrophages but was inhibited by PF74 in HOS cells. A truncated version of CPSF6 (CPSF6-358) inhibited SIVcpzPtt and HIV-1, while in contrast, SIVcpzPts was resistant to CPSF6-358. Homology modeling of HIV-1, SIVcpzPtt, and SIVcpzPts capsids and binding energy estimates suggest that these three viruses bind similarly to the host proteins cyclophilin A (CYPA) and CPSF6 as well as the capsid inhibitor PF74. Cyclosporine treatment, mutation of the CYPA-binding loop in the capsid, or CYPA knockout eliminated the resistance of SIVcpzPts to PF74 in HeLa cells. These experiments revealed that the antiviral capacity of PF74 is controlled by CYPA in a virus- and cell type-specific manner. Our data indicate that SIVcpz viruses can use infection pathways that escape the antiviral activity of PF74. We further suggest that the antiviral activity of PF74 capsid inhibitors depends on cellular cofactors.IMPORTANCE HIV-1 originated from SIVcpzPtt but not from the related virus SIVcpzPts, and thus, it is important to describe molecular infection by SIVcpzPts in human cells to understand the zoonosis of SIVs. Pharmacological HIV-1 capsid inhibitors (e.g., PF74) bind a capsid groove that is also a binding site for the cellular protein CPSF6. SIVcpzPts was resistant to PF74 in HeLa cells but sensitive in HOS cells, thus indicating cell line-specific resistance. Both SIVcpz viruses showed resistance to PF74 in human PBMCs. Modulating the presence of cyclophilin A or its binding to capsid in HeLa cells overcame SIVcpzPts resistance to PF74. These results indicate that early cytoplasmic infection events of SIVcpzPts may differ between cell types and affect, in an unknown manner, the antiviral activity of capsid inhibitors. Thus, capsid inhibitors depend on the activity or interaction of currently uncharacterized cellular factors.