Rapid, Precise, and Reproducible Prediction of Peptide-MHC Binding Affinities from Molecular Dynamics That Correlate Well with Experiment

Optimal Molecular Design: Generative Active Learning Combining REINVENT with Precise Binding Free Energy Ranking Simulations

Active learning (AL) is a specific instance of sequential experimental design and uses machine learning to intelligently choose the next data point or batch of molecular structures to be evaluated. In this sense, it closely mimics the iterative design-make-test-analysis cycle of laboratory experiments to find optimized compounds for a given design task. Here, we describe an AL protocol which combines generative molecular AI, using REINVENT, and physics-based absolute binding free energy molecular dynamics simulation, using ESMACS, to discover new ligands for two different target proteins, 3CLpro and TNKS2. We have deployed our generative active learning (GAL) protocol on Frontier, the world's only exa-scale machine. We show that the protocol can find higher-scoring molecules compared to the baseline, a surrogate ML docking model for 3CLpro and compounds with experimentally determined binding affinities for TNKS2. The ligands found are also chemically diverse and occupy a different chemical space than the baseline. We vary the batch sizes that are put forward for free energy assessment in each GAL cycle to assess the impact on their efficiency on the GAL protocol and recommend their optimal values in different scenarios. Overall, we demonstrate a powerful capability of the combination of physics-based and AI methods which yields effective chemical space sampling at an unprecedented scale and is of immediate and direct relevance to modern, data-driven drug discovery.

PepBinding: A Workflow for Predicting Peptide Binding Structures by Combining Peptide Docking and Peptide Gaussian Accelerated Molecular Dynamics Simulations

Predicting protein-peptide interactions is crucial for understanding peptide binding processes and designing peptide drugs. However, traditional computational modeling approaches face challenges in accurately predicting peptide-protein binding structures due to the slow dynamics and high flexibility of the peptides. Here, we introduce a new workflow termed "PepBinding" for predicting peptide binding structures, which combines peptide docking, all-atom enhanced sampling simulations using the Peptide Gaussian accelerated Molecular Dynamics (Pep-GaMD) method, and structural clustering. PepBinding has been demonstrated on seven distinct model peptides. In peptide docking using HPEPDOCK, the peptide backbone root-mean-square deviations (RMSDs) of their bound conformations relative to X-ray structures ranged from 3.8 to 16.0 Å, corresponding to the medium to inaccurate quality models according to the Critical Assessment of PRediction of Interactions (CAPRI) criteria. The Pep-GaMD simulations performed for only 200 ns significantly improved the docking models, resulting in five medium and two acceptable quality models. Therefore, PepBinding is an efficient workflow for predicting peptide binding structures and is publicly available at https://github.com/MiaoLab20/PepBinding.

Introduction to Computational Biomedicine

The domain of computational biomedicine is a new and burgeoning one. Its areas of concern cover all scales of human biology, physiology, and pathology, commonly referred to as medicine, from the genomic to the whole human and beyond, including epidemiology and population health. Computational biomedicine aims to provide high-fidelity descriptions and predictions of the behavior of biomedical systems of both fundamental scientific and clinical importance. Digital twins and virtual humans aim to reproduce the extremely accurate duplicate of real-world human beings in cyberspace, which can be used to make highly accurate predictions that take complicated conditions into account. When that can be done reliably enough for the predictions to be actionable, such an approach will make an impact in the pharmaceutical industry by reducing or even replacing the extremely laboratory-intensive preclinical process of making and testing compounds in laboratories, and in clinical applications by assisting clinicians to make diagnostic and treatment decisions.

Ensemble-Based Approaches Ensure Reliability and Reproducibility

It is increasingly widely recognized that ensemble-based approaches are required to achieve reliability, accuracy, and precision in molecular dynamics calculations. The purpose of the present article is to address a frequently raised question: what is the optimal way to perform ensemble simulation to calculate quantities of interest?

Structure-Based Multilevel Descriptors for High-throughput Screening of Elastomers

To discover new materials, high-throughput screening (HTS) with machine learning (ML) requires universally available descriptors that can accurately predict the desired properties. For elastomers, experimental and simulation data in current descriptors may not be available for all candidates of interest, hindering elastomer discovery through HTS. To address this challenge, we introduce structure-based multilevel (SM) descriptors of elastomers derived solely from molecular structure that is universally available. Our SM descriptors are hierarchically organized to capture both local soft and hard segment structures as well as the global structures of elastomers. With the SM-Morgan Fingerprint (SM-MF) descriptor, one of our SM descriptors, a machine learning model accurately predicts elastomer toughness with a remarkable accuracy of 0.91. Furthermore, an HTS pipeline is established to swiftly screen elastomers with targeted toughness. We also demonstrate the generality and applicability of SM descriptors by using them to construct HTS pipelines for screening elastomers with a targeted critical strain or Young's modulus. The user-friendliness and low computational cost of SM descriptors make them a promising tool to significantly enhance HTS in the search for novel materials.

Design, synthesis, computational study and cytotoxic evaluation of some new quinazoline derivatives containing pyrimidine moiety

Quinazoline derivatives, as an important category of heterocyclic compounds, have received much attention for the design and development of new drugs due to their various pharmacological properties. Besides, there is a great deal of evidence showing pyrimidine analogs as anticancer agents. Thus, in the present study, for the design of new target compounds with cytotoxic activity, we focused on various quinazolinone and pyrimidine hybrids. A new series of quinazoline-pyrimidine hybrid derivatives (6a-6n) have been designed and synthesized as novel antiproliferative agents. All the synthesized compounds characterized based on their IR, NMR and Mass spectroscopic data. Antiproliferative activities of the new compounds were evaluated against three human cancer cell lines (MCF-7, A549, SW-480). The compounds were found to have appropriate potential with IC50 values ranging from 2.3 ± 5.91 to 176.5 ± 0.7 μM against the tested cell lines. Compound 6n exerted the highest antiproliferative activity with IC50 values of 5.9 ± 1.69 μM, 2.3 ± 5.91 μM and 5.65 ± 2.33 μM against A549, SW-480 and MCF-7 respectively. The results indicated that 6n could induce apoptosis in A549 cell line in a dose dependent manner and arrest in the S phase of cell cycle. Docking studies were also done to investigate the detailed binding pattern of the synthesized compounds against EGFR. Furthermore, molecular dynamic simulation and binding free energy calculation have been done to rescore initial docking pose of the synthesized compounds using ensemble-based MMGB/PBSA free energy method. According to the results, free energy calculation confirmed biological activity of compounds and also, Arg 817 and Lys 721 residues had the pivotal role in the high potency of 6n. Finally, the drug likeness and in silico ADME study were also predicted.

Structure and dynamics of an archetypal DNA nanoarchitecture revealed via cryo-EM and molecular dynamics simulations

DNA can be folded into rationally designed, unique, and functional materials. To fully realise the potential of these DNA materials, a fundamental understanding of their structure and dynamics is necessary, both in simple solvents as well as more complex and diverse anisotropic environments. Here we analyse an archetypal six-duplex DNA nanoarchitecture with single-particle cryo-electron microscopy and molecular dynamics simulations in solvents of tunable ionic strength and within the anisotropic environment of biological membranes. Outside lipid bilayers, the six-duplex bundle lacks the designed symmetrical barrel-type architecture. Rather, duplexes are arranged in non-hexagonal fashion and are disorted to form a wider, less elongated structure. Insertion into lipid membranes, however, restores the anticipated barrel shape due to lateral duplex compression by the bilayer. The salt concentration has a drastic impact on the stability of the inserted barrel-shaped DNA nanopore given the tunable electrostatic repulsion between the negatively charged duplexes. By synergistically combining experiments and simulations, we increase fundamental understanding into the environment-dependent structural dynamics of a widely used nanoarchitecture. This insight will pave the way for future engineering and biosensing applications.

Computational Approaches Drive Developments in Immune-Oncology Therapies for PD-1/PD-L1 Immune Checkpoint Inhibitors

Computational approaches in immune-oncology therapies focus on using data-driven methods to identify potential immune targets and develop novel drug candidates. In particular, the search for PD-1/PD-L1 immune checkpoint inhibitors (ICIs) has enlivened the field, leveraging the use of cheminformatics and bioinformatics tools to analyze large datasets of molecules, gene expression and protein-protein interactions. Up to now, there is still an unmet clinical need for improved ICIs and reliable predictive biomarkers. In this review, we highlight the computational methodologies applied to discovering and developing PD-1/PD-L1 ICIs for improved cancer immunotherapies with a greater focus in the last five years. The use of computer-aided drug design structure- and ligand-based virtual screening processes, molecular docking, homology modeling and molecular dynamics simulations methodologies essential for successful drug discovery campaigns focusing on antibodies, peptides or small-molecule ICIs are addressed. A list of recent databases and web tools used in the context of cancer and immunotherapy has been compilated and made available, namely regarding a general scope, cancer and immunology. In summary, computational approaches have become valuable tools for discovering and developing ICIs. Despite significant progress, there is still a need for improved ICIs and biomarkers, and recent databases and web tools have been compiled to aid in this pursuit.

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.

Accuracy and precision of binding free energy prediction for a tacrine related lead inhibitor of acetylcholinesterase with an arsenal of supercomputerized molecular modelling methods: a comparative study

Nowadays, advanced computational chemistry methods offer various strategies for revealing prospective hit structures in drug development essentially through accurate binding free energy predictions. After the era of molecular docking and quantitative structure-activity relationships, much interest has been lately oriented to perturbed molecular dynamic approaches like replica exchange with solute tempering and free energy perturbation (REST/FEP) and the potential of the mean force with adaptive biasing and accelerated weight histograms (PMF/AWH). Both of these receptor-based techniques can exploit exascale CPU&GPU supercomputers to achieve high throughput performance. In this fundamental study, we have compared the predictive power of a panel of supercomputerized molecular modelling methods to distinguish the major binding modes and the corresponding binding free energies of a promising tacrine related potential antialzheimerics in human acetylcholinesterase. The binding free energies were estimated using flexible molecular docking, molecular mechanics/generalized Born surface area/Poisson-Boltzmann surface area (MM/GBSA/PBSA), transmutation REST/FEP with 12 x 5 ns/λ windows, annihilation FEP with 20 x 5 ns/λ steps, PMF with weight histogram analysis method (WHAM) and 40 x 5 ns samples, and PMF/AWH with 10 x 100 ns replicas. Confrontation of the classical approaches such as canonical molecular dynamics and molecular docking with alchemical calculations and steered molecular dynamics enabled us to show how large errors in ΔG predictions can be expected if these in silico methods are employed in the elucidation of a common case of enzyme inhibition.Communicated by Ramaswamy H. Sarma.

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.

In-silico molecular modelling, MM/GBSA binding free energy and molecular dynamics simulation study of novel pyrido fused imidazo[4,5-c]quinolines as potential anti-tumor agents

With an alarming increase in the number of cancer patients and a variety of tumors, it is high time for intensive investigation on more efficient and potent anti-tumor agents. Though numerous agents have enriched the literature, still there exist challenges, with the availability of different targets and possible cross-reactivity. Herein we have chosen the phosphoinositide 3-kinase (PI3K) as the target of interest and investigated the potential of pyrido fused imidazo[4,5-c]quinoline derivatives to bind strongly to the active site, thereby inhibiting the progression of various types of tumors. The AutoDock, Glide and the Prime-MM/GBSA analysis are used to execute the molecular docking investigation and validation for the designed compounds. The anti-tumor property evaluations were carried out by using PASS algorithm. Based on the GLIDE score, the binding affinity of the designed molecules towards the target PI3K was evaluated. The energetics associated with static interactions revealed 1j as the most potential candidate and the dynamic investigations including RMSD, RMSF, Rg, SASA and hydrogen bonding also supported the same through relative stabilization induced through ligand interactions. Subsequently, the binding free energy of the Wortmannin and 1j complex calculated using MM-PBSA analysis. Further evaluations with PASS prediction algorithm also supported the above results. The studies reveal that there is evidence for considering appropriate pyrido fused imidazo[4,5-c]quinoline compounds as potential anti-tumor agents.

Charge-based interactions through peptide position 4 drive diversity of antigen presentation by human leukocyte antigen class I molecules

Human leukocyte antigen class I (HLA-I) molecules bind and present peptides at the cell surface to facilitate the induction of appropriate CD8+ T cell-mediated immune responses to pathogen- and self-derived proteins. The HLA-I peptide-binding cleft contains dominant anchor sites in the B and F pockets that interact primarily with amino acids at peptide position 2 and the C-terminus, respectively. Nonpocket peptide-HLA interactions also contribute to peptide binding and stability, but these secondary interactions are thought to be unique to individual HLA allotypes or to specific peptide antigens. Here, we show that two positively charged residues located near the top of peptide-binding cleft facilitate interactions with negatively charged residues at position 4 of presented peptides, which occur at elevated frequencies across most HLA-I allotypes. Loss of these interactions was shown to impair HLA-I/peptide binding and complex stability, as demonstrated by both in vitro and in silico experiments. Furthermore, mutation of these Arginine-65 (R65) and/or Lysine-66 (K66) residues in HLA-A*02:01 and A*24:02 significantly reduced HLA-I cell surface expression while also reducing the diversity of the presented peptide repertoire by up to 5-fold. The impact of the R65 mutation demonstrates that nonpocket HLA-I/peptide interactions can constitute anchor motifs that exert an unexpectedly broad influence on HLA-I-mediated antigen presentation. These findings provide fundamental insights into peptide antigen binding that could broadly inform epitope discovery in the context of viral vaccine development and cancer immunotherapy.

The performance of ensemble-based free energy protocols in computing binding affinities to ROS1 kinase

Optimization of binding affinities for compounds to their target protein is a primary objective in drug discovery. Herein we report on a collaborative study that evaluates a set of compounds binding to ROS1 kinase. We use ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling) protocols to rank the binding free energies. The predicted binding free energies from ESMACS simulations show good correlations with experimental data for subsets of the compounds. Consistent binding free energy differences are generated for TIES and ESMACS. Although an unexplained overestimation exists, we obtain excellent statistical rankings across the set of compounds from the TIES protocol, with a Pearson correlation coefficient of 0.90 between calculated and experimental activities.

Inhibition constant and stability of tripeptide inhibitors of gut trypsin-like enzyme of the soybean pest Anticarsia gemmatalis

Insects overcome the action of natural protease inhibitors (PIs) due to evolutionary adaptations through endogenous proteolysis and reprogramming proteases. Insect adaptations complicate the formulation of IP-based crop protection products. However, small peptides designed based on the active site of enzymes have shown promising results that could change this scenario. GORE1 and GORE2 are designed tripeptides that reduce the survival of Anticarsia gemmatalis when ingested orally. In this article, the stability and ability of the peptides to bind trypsin-like enzymes of A. gemmatalis were evaluated by molecular dynamics (MD) simulations. The ability of the peptides to inhibit trypsin-like enzymes in vivo was compared with the SKTI protein by feeding A. gemmatalis larvae at different concentrations, followed by an inhibition persistence assay. During the MD simulation of enzyme-ligand complexes, both peptides showed a small variation of root-mean-square deviation and root-mean-square fluctuation, suggesting that these molecules reach equilibrium when forming a complex with the trypsin-like enzyme. Furthermore, both peptides form hydrogen bonds with substrate recognition sites of A. gemmatalis trypsin-like enzyme, with GORE2 having more interactions than GORE1. Larvae of A. gemmatalis exposed to the peptides and SKTI showed a similar reduction in proteolytic activity, but the persistence of inhibition of trypsin-like enzyme was longer in peptide-fed insects. Despite their size, the peptides exhibit important active and substrate binding site interactions, stability during complex formation, and steadiness effects in vivo. The results provide fundamental information for the development of mimetic molecules and help in decision-making for the selection of delivery methods for larger-scale experiments regarding similar molecules.

Multiple-Allele MHC Class II Epitope Engineering by a Molecular Dynamics-Based Evolution Protocol

Epitopes that bind simultaneously to all human alleles of Major Histocompatibility Complex class II (MHC II) are considered one of the key factors for the development of improved vaccines and cancer immunotherapies. To engineer MHC II multiple-allele binders, we developed a protocol called PanMHC-PARCE, based on the unsupervised optimization of the epitope sequence by single-point mutations, parallel explicit-solvent molecular dynamics simulations and scoring of the MHC II-epitope complexes. The key idea is accepting mutations that not only improve the affinity but also reduce the affinity gap between the alleles. We applied this methodology to enhance a Plasmodium vivax epitope for multiple-allele binding. In vitro rate-binding assays showed that four engineered peptides were able to bind with improved affinity toward multiple human MHC II alleles. Moreover, we demonstrated that mice immunized with the peptides exhibited interferon-gamma cellular immune response. Overall, the method enables the engineering of peptides with improved binding properties that can be used for the generation of new immunotherapies.

Ensemble Simulations and Experimental Free Energy Distributions: Evaluation and Characterization of Isoxazole Amides as SMYD3 Inhibitors

Optimization of binding affinities for ligands to their target protein is a primary objective in rational drug discovery. Herein, we report on a collaborative study that evaluates various compounds designed to bind to the SET and MYND domain-containing protein 3 (SMYD3). SMYD3 is a histone methyltransferase and plays an important role in transcriptional regulation in cell proliferation, cell cycle, and human carcinogenesis. Experimental measurements using the scintillation proximity assay show that the distributions of binding free energies from a large number of independent measurements exhibit non-normal properties. We use ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and TIES (thermodynamic integration with enhanced sampling) protocols to predict the binding free energies and to provide a detailed chemical insight into the nature of ligand-protein binding. Our results show that the 1-trajectory ESMACS protocol works well for the set of ligands studied here. Although one unexplained outlier exists, we obtain excellent statistical ranking across the set of compounds from the ESMACS protocol and good agreement between calculations and experiments for the relative binding free energies from the TIES protocol. ESMACS and TIES are again found to be powerful protocols for the accurate comparison of the binding free energies.

A Bionic-Homodimerization Strategy for Optimizing Modulators of Protein-Protein Interactions: From Statistical Mechanics Theory to Potential Clinical Translation

Emerging protein-protein interaction (PPI) modulators have brought out exciting ability as therapeutics in human diseases, but its clinical translation has been greatly hampered by the limited affinity. Inspired by the homodimerize structure of antibody, the homodimerization contributes hugely to generating the optimized affinity is conjectured. Herein, a statistical-mechanics-theory-guided method is established to quantize the affinity of ligands with different topologies through analyzing the change of enthalpy and the loss of translational and rotational entropies. A peptide modulator for p53-MDM2 termed CPAP is used to homodimerize connecting, and this simple homodimerization can significantly increase the affinity. To realize the cellular internalization and tumor accumulation, Dimer CPAP and Mono CPAP are nanoengineered into gold(I)-CPAP supermolecule by the aurophilic interaction-driven self-assembly. Nano-Dimer CPAP potently suppressed tumor growth in lung cancer allograft model and a patient-derived xenograft model in more action than Nano-Mono CPAP, while keeping a favorable drug safety profile. This work not only presents a physico-mechanical method for calculating the affinity of PPI modulators, but also provides a simple yet robust homodimerization strategy to optimize the affinity of PPI modulators.

Alchemical Free Energy Methods Applied to Complexes of the First Bromodomain of BRD4

Accurate and rapid predictions of the binding affinity of a compound to a target are one of the ultimate goals of computer aided drug design. Alchemical approaches to free energy estimations follow the path from an initial state of the system to the final state through alchemical changes of the energy function during a molecular dynamics simulation. Herein, we explore the accuracy and efficiency of two such techniques: relative free energy perturbation (FEP) and multisite lambda dynamics (MSλD). These are applied to a series of inhibitors for the bromodomain-containing protein 4 (BRD4). We demonstrate a procedure for obtaining accurate relative binding free energies using MSλD when dealing with a change in the net charge of the ligand. This resulted in an impressive comparison with experiment, with an average difference of 0.4 ± 0.4 kcal mol-1. In a benchmarking study for the relative FEP calculations, we found that using 20 lambda windows with 0.5 ns of equilibration and 1 ns of data collection for each window gave the optimal compromise between accuracy and speed. Overall, relative FEP and MSλD predicted binding free energies with comparable accuracy, an average of 0.6 kcal mol-1 for each method. However, MSλD makes predictions for a larger molecular space over a much shorter time scale than relative FEP, with MSλD requiring a factor of 18 times less simulation time for the entire molecule space.

Assessing the Performance of Screening MM/PBSA in Protein-Ligand Interactions

Accurate calculation of the binding free energies between a protein and a ligand is the primary objective of structure-based drug design, but it still remains a challenging problem. In this work, we apply the screening molecular mechanics/Poisson Boltzmann surface area (MM/PBSA) method to calculate the binding affinity of protein-ligand interactions. Our results show that the performance of the screening MM/PBSA is better than that of the standard MM/PBSA, especially in a charged-ligand system. In addition, we also investigate the effect of the solute dielectric constant on the results, and find that the optimal solute dielectric constants are different between the neutral-ligand system and the charged-ligand system. Moreover, we also evaluate the effect of the atomic-charge methods on the performance of the screening MM/PBSA. The present study demonstrates that the screening MM/PBSA should be a reliable method for calculating binding energy of biosystems.

Thermodynamic and structural insights into the repurposing of drugs that bind to SARS-CoV-2 main protease

Although researchers have been working tirelessly since the COVID-19 outbreak, so far only three drugs - remdesivir, ronapreve and molnupiravir - have been approved for use in some countries which directly target the SARS-CoV-2 virus. Given the slow pace and substantial costs of new drug discovery and development, together with the urgency of the matter, repurposing of existing drugs for the ongoing disease is an attractive proposition. In a recent study, a high-throughput X-ray crystallographic screen was performed for a selection of drugs which have been approved or are in clinical trials. Thirty-seven compounds have been identified from drug libraries all of which bind to the SARS-CoV-2 main protease (3CL^pro). In the current study, we use molecular dynamics simulation and an ensemble-based free energy approach, namely, enhanced sampling of molecular dynamics with approximation of continuum solvent (ESMACS), to investigate a subset of the aforementioned compounds. The drugs studied here are highly diverse, interacting with different binding sites and/or subsites of 3CL^pro. The predicted free energies are compared with experimental results wherever they are available and they are found to be in excellent agreement. Our study also provides detailed energetic insights into the nature of the associated drug-protein binding, in turn shedding light on the design and discovery of potential drugs.

Reliable In Silico Ranking of Engineered Therapeutic TCR Binding Affinities with MMPB/GBSA

Accurate and efficient in silico ranking of protein–protein binding affinities is useful for protein design with applications in biological therapeutics. One popular approach to rank binding affinities is to apply the molecular mechanics Poisson–Boltzmann/generalized Born surface area (MMPB/GBSA) method to molecular dynamics (MD) trajectories. Here, we identify protocols that enable the reliable evaluation of T-cell receptor (TCR) variants binding to their target, peptide-human leukocyte antigens (pHLAs). We suggest different protocols for variant sets with a few (≤4) or many mutations, with entropy corrections important for the latter. We demonstrate how potential outliers could be identified in advance and that just 5–10 replicas of short (4 ns) MD simulations may be sufficient for the reproducible and accurate ranking of TCR variants. The protocols developed here can be applied toward in silico screening during the optimization of therapeutic TCRs, potentially reducing both the cost and time taken for biologic development.

Pandemic drugs at pandemic speed: infrastructure for accelerating COVID-19 drug discovery with hybrid machine learning- and physics-based simulations on high-performance computers

The race to meet the challenges of the global pandemic has served as a reminder that the existing drug discovery process is expensive, inefficient and slow. There is a major bottleneck screening the vast number of potential small molecules to shortlist lead compounds for antiviral drug development. New opportunities to accelerate drug discovery lie at the interface between machine learning methods, in this case, developed for linear accelerators, and physics-based methods. The two in silico methods, each have their own advantages and limitations which, interestingly, complement each other. Here, we present an innovative infrastructural development that combines both approaches to accelerate drug discovery. The scale of the potential resulting workflow is such that it is dependent on supercomputing to achieve extremely high throughput. We have demonstrated the viability of this workflow for the study of inhibitors for four COVID-19 target proteins and our ability to perform the required large-scale calculations to identify lead antiviral compounds through repurposing on a variety of supercomputers.

Principles of Small-Molecule Transport through Synthetic Nanopores

Synthetic nanopores made from DNA replicate the key biological processes of transporting molecular cargo across lipid bilayers. Understanding transport across the confined lumen of the nanopores is of fundamental interest and of relevance to their rational design for biotechnological applications. Here we reveal the transport principles of organic molecules through DNA nanopores by synergistically combining experiments and computer simulations. Using a highly parallel nanostructured platform, we synchronously measure the kinetic flux across hundreds of individual pores to obtain rate constants. The single-channel transport kinetics are close to the theoretical maximum, while selectivity is determined by the interplay of cargo charge and size, the pores' sterics and electrostatics, and the composition of the surrounding lipid bilayer. The narrow distribution of transport rates implies a high structural homogeneity of DNA nanopores. The molecular passageway through the nanopore is elucidated via coarse-grained constant-velocity steered molecular dynamics simulations. The ensemble simulations pinpoint with high resolution and statistical validity the selectivity filter within the channel lumen and determine the energetic factors governing transport. Our findings on these synthetic pores' structure-function relationship will serve to guide their rational engineering to tailor transport selectivity for cell biological research, sensing, and drug delivery.

Locating ligand binding sites in G-protein coupled receptors using combined information from docking and sequence conservation

GPCRs (G-protein coupled receptors) are the largest family of drug targets and share a conserved structure. Binding sites are unknown for many important GPCR ligands due to the difficulties of GPCR recombinant expression, biochemistry, and crystallography. We describe our approach, ConDockSite, for predicting ligand binding sites in class A GPCRs using combined information from surface conservation and docking, starting from crystal structures or homology models. We demonstrate the effectiveness of ConDockSite on crystallized class A GPCRs such as the beta2 adrenergic and A2A adenosine receptors. We also demonstrate that ConDockSite successfully predicts ligand binding sites from high-quality homology models. Finally, we apply ConDockSite to predict the ligand binding sites on a structurally uncharacterized GPCR, GPER, the G-protein coupled estrogen receptor. Most of the sites predicted by ConDockSite match those found in other independent modeling studies. ConDockSite predicts that four ligands bind to a common location on GPER at a site deep in the receptor cleft. Incorporating sequence conservation information in ConDockSite overcomes errors introduced from physics-based scoring functions and homology modeling.

Understanding the avidin-biotin binding based on polarized protein-specific charge

In this study, charge updating schemes based on the local polarized protein-specific charge (LPPC) were introduced to vary the atomic charges of the biotin molecule and the residues in close contact during the simulation of the avidin-biotin complexes. The need of the charge variation of the ligand in response to changes in its surroundings was thoroughly studied. The results show that the calculated binding energy difference between biotin (BTN1) and 2'-iminobiotin (BTN2) and avidin is in excellent agreement with the experimental value, thus verifying the feasibility of updating the atomic charges of ligands during the simulation.

The effect of protein mutations on drug binding suggests ensuing personalised drug selection

The advent of personalised medicine promises a deeper understanding of mechanisms and therefore therapies. However, the connection between genomic sequences and clinical treatments is often unclear. We studied 50 breast cancer patients belonging to a population-cohort in the state of Qatar. From Sanger sequencing, we identified several new deleterious mutations in the estrogen receptor 1 gene (ESR1). The effect of these mutations on drug treatment in the protein target encoded by ESR1, namely the estrogen receptor, was achieved via rapid and accurate protein-ligand binding affinity interaction studies which were performed for the selected drugs and the natural ligand estrogen. Four nonsynonymous mutations in the ligand-binding domain were subjected to molecular dynamics simulation using absolute and relative binding free energy methods, leading to the ranking of the efficacy of six selected drugs for patients with the mutations. Our study shows that a personalised clinical decision system can be created by integrating an individual patient's genomic data at the molecular level within a computational pipeline which ranks the efficacy of binding of particular drugs to variant proteins.

Artificial intelligence in drug design: algorithms, applications, challenges and ethics

The discovery paradigm of drugs is rapidly growing due to advances in machine learning (ML) and artificial intelligence (AI). This review covers myriad faces of AI and ML in drug design. There is a plethora of AI algorithms, the most common of which are summarized in this review. In addition, AI is fraught with challenges that are highlighted along with plausible solutions to them. Examples are provided to illustrate the use of AI and ML in drug discovery and in predicting drug properties such as binding affinities and interactions, solubility, toxicology, blood–brain barrier permeability and chemical properties. The review also includes examples depicting the implementation of AI and ML in tackling intractable diseases such as COVID-19, cancer and Alzheimer’s disease. Ethical considerations and future perspectives of AI are also covered in this review.

Uncertainty quantification in classical molecular dynamics

Molecular dynamics simulation is now a widespread approach for understanding complex systems on the atomistic scale. It finds applications from physics and chemistry to engineering, life and medical science. In the last decade, the approach has begun to advance from being a computer-based means of rationalizing experimental observations to producing apparently credible predictions for a number of real-world applications within industrial sectors such as advanced materials and drug discovery. However, key aspects concerning the reproducibility of the method have not kept pace with the speed of its uptake in the scientific community. Here, we present a discussion of uncertainty quantification for molecular dynamics simulation designed to endow the method with better error estimates that will enable it to be used to report actionable results. The approach adopted is a standard one in the field of uncertainty quantification, namely using ensemble methods, in which a sufficiently large number of replicas are run concurrently, from which reliable statistics can be extracted. Indeed, because molecular dynamics is intrinsically chaotic, the need to use ensemble methods is fundamental and holds regardless of the duration of the simulations performed. We discuss the approach and illustrate it in a range of applications from materials science to ligand-protein binding free energy estimation. This article is part of the theme issue 'Reliability and reproducibility in computational science: implementing verification, validation and uncertainty quantification in silico'.

MHCBI: a pipeline for calculating peptide-MHC binding energy using semi-empirical quantum mechanical methods with explicit/implicit solvent models

Experimentally estimating peptide-major histocompatibility complex (pMHC) binding affinity has been quite challenging due to the many receptors and the many potential ligands implicated in it. We have thus proposed a straightforward computational methodology considering the different mechanisms involved in pMHC binding to facilitate studying such receptor-ligand interactions. We have developed a pipeline using semi-empirical quantum mechanical methods for calculating pMHC class I and II molecules' binding energy (BE). This pipeline can systematize the methodology for calculating pMHC system BE, enabling the rational design of T-cell epitopes to be used as pharmaceuticals and vaccines.

Peptide Gaussian accelerated molecular dynamics (Pep-GaMD): Enhanced sampling and free energy and kinetics calculations of peptide binding

Peptides mediate up to 40% of known protein-protein interactions in higher eukaryotes and play an important role in cellular signaling. However, it is challenging to simulate both binding and unbinding of peptides and calculate peptide binding free energies through conventional molecular dynamics, due to long biological timescales and extremely high flexibility of the peptides. Based on the Gaussian accelerated molecular dynamics (GaMD) enhanced sampling technique, we have developed a new computational method "Pep-GaMD," which selectively boosts essential potential energy of the peptide in order to effectively model its high flexibility. In addition, another boost potential is applied to the remaining potential energy of the entire system in a dual-boost algorithm. Pep-GaMD has been demonstrated on binding of three model peptides to the SH3 domains. Independent 1 µs dual-boost Pep-GaMD simulations have captured repetitive peptide dissociation and binding events, which enable us to calculate peptide binding thermodynamics and kinetics. The calculated binding free energies and kinetic rate constants agreed very well with available experimental data. Furthermore, the all-atom Pep-GaMD simulations have provided important insights into the mechanism of peptide binding to proteins that involves long-range electrostatic interactions and mainly conformational selection. In summary, Pep-GaMD provides a highly efficient, easy-to-use approach for unconstrained enhanced sampling and calculations of peptide binding free energies and kinetics.

Unraveling the Allosteric Communication Mechanisms in T-Cell Receptor-Peptide-Loaded Major Histocompatibility Complex Dynamics Using Molecular Dynamics Simulations: An Approach Based on Dynamic Cross Correlation Maps and Residue Interaction Energy Calculations

Antigen presentation by major histocompatibility complex (MHC) proteins to T-cell receptors (TCRs) plays a crucial role in triggering the adaptive immune response. Most of our knowledge on TCR-peptide-loaded major histocompatibility complex (pMHC) interaction stemmed from experiments yielding static structures, yet the dynamic aspects of this molecular interaction are equally important to understand the underlying molecular mechanisms and to develop treatment strategies against diseases such as cancer and autoimmune diseases. To this end, computational biophysics studies including all-atom molecular dynamics simulations have provided useful insights; however, we still lack a basic understanding of an overall allosteric mechanism that results in conformational changes in the TCR and subsequent T-cell activation. Previous hydrogen-deuterium exchange and nuclear magnetic resonance studies provided clues regarding these molecular mechanisms, including global rigidification and allosteric effects on the constant domain of TCRs away from the pMHC interaction site. Here, we show that molecular dynamics simulations can be used to identify how this overall rigidification may be related to the allosteric communication within TCRs upon pMHC interaction via essential dynamics and nonbonded residue-residue interaction energy analyses. The residues taking part in the rigidification effect are highlighted with an intricate analysis on residue interaction changes, which lead to a detailed outline of the complex formation event. Our results indicate that residues of the Cβ domain of TCRs show significant differences in their nonbonded interactions upon complex formation. Moreover, the dynamic cross correlations between these residues are also increased, in line with their nonbonded interaction energy changes. Altogether, our approach may be valuable for elucidating intramolecular allosteric changes in the TCR structure upon pMHC interaction in molecular dynamics simulations.

Specificity of bispecific T cell receptors and antibodies targeting peptide-HLA

Tumor-associated peptide-human leukocyte antigen complexes (pHLAs) represent the largest pool of cell surface-expressed cancer-specific epitopes, making them attractive targets for cancer therapies. Soluble bispecific molecules that incorporate an anti-CD3 effector function are being developed to redirect T cells against these targets using 2 different approaches. The first achieves pHLA recognition via affinity-enhanced versions of natural TCRs (e.g., immune-mobilizing monoclonal T cell receptors against cancer [ImmTAC] molecules), whereas the second harnesses an antibody-based format (TCR-mimic antibodies). For both classes of reagent, target specificity is vital, considering the vast universe of potential pHLA molecules that can be presented on healthy cells. Here, we made use of structural, biochemical, and computational approaches to investigate the molecular rules underpinning the reactivity patterns of pHLA-targeting bispecifics. We demonstrate that affinity-enhanced TCRs engage pHLA using a comparatively broad and balanced energetic footprint, with interactions distributed over several HLA and peptide side chains. As ImmTAC molecules, these TCRs also retained a greater degree of pHLA selectivity, with less off-target activity in cellular assays. Conversely, TCR-mimic antibodies tended to exhibit binding modes focused more toward hot spots on the HLA surface and exhibited a greater degree of crossreactivity. Our findings extend our understanding of the basic principles that underpin pHLA selectivity and exemplify a number of molecular approaches that can be used to probe the specificity of pHLA-targeting molecules, aiding the development of future reagents.

Effects of side chain length of 10-methyl-aplog-1, a simplified analog of debromoaplysiatoxin, on PKC binding, anti-proliferative, and pro-inflammatory activities

10-Methyl-aplog-1 (1), a simplified analog of debromoaplysiatoxin, exhibits a high binding affinity for protein kinase C (PKC) isozymes and potent antiproliferative activity against several cancer cells with few adverse effects. A recent study has suggested that its phenol group in the side chain is involved in hydrogen bonding and CH/π interactions with the binding cleft-forming loops in the PKCδ-C1B domain. To clarify the effects of the side chain length on these interactions, four analogs of 1 with various lengths of side chains (2-5) were prepared. The maximal PKC binding affinity and antiproliferative activity were observed in 1. Remarkably, the introduction of a bromine atom into the phenol group of 2 increased not only these activities but also proinflammatory activity. These results indicated that 1 has the optimal side chain length as an anticancer seed. This conclusion was supported by docking simulations of 1-5 to the PKCδ-C1B domain.

Rapid, accurate, precise and reproducible ligand-protein binding free energy prediction

A central quantity of interest in molecular biology and medicine is the free energy of binding of a molecule to a target biomacromolecule. Until recently, the accurate prediction of binding affinity had been widely regarded as out of reach of theoretical methods owing to the lack of reproducibility of the available methods, not to mention their complexity, computational cost and time-consuming procedures. The lack of reproducibility stems primarily from the chaotic nature of classical molecular dynamics (MD) and the associated extreme sensitivity of trajectories to their initial conditions. Here, we review computational approaches for both relative and absolute binding free energy calculations, and illustrate their application to a diverse set of ligands bound to a range of proteins with immediate relevance in a number of medical domains. We focus on ensemble-based methods which are essential in order to compute statistically robust results, including two we have recently developed, namely thermodynamic integration with enhanced sampling and enhanced sampling of MD with an approximation of continuum solvent. Together, these form a set of rapid, accurate, precise and reproducible free energy methods. They can be used in real-world problems such as hit-to-lead and lead optimization stages in drug discovery, and in personalized medicine. These applications show that individual binding affinities equipped with uncertainty quantification may be computed in a few hours on a massive scale given access to suitable high-end computing resources and workflow automation. A high level of accuracy can be achieved using these approaches.

Hit-to-lead and lead optimization binding free energy calculations for G protein-coupled receptors

We apply the hit-to-lead ESMACS (enhanced sampling of molecular dynamics with approximation of continuum solvent) and lead-optimization TIES (thermodynamic integration with enhanced sampling) methods to compute the binding free energies of a series of ligands at the A1 and A2A adenosine receptors, members of a subclass of the GPCR (G protein-coupled receptor) superfamily. Our predicted binding free energies, calculated using ESMACS, show a good correlation with previously reported experimental values of the ligands studied. Relative binding free energies, calculated using TIES, accurately predict experimentally determined values within a mean absolute error of approximately 1 kcal mol-1. Our methodology may be applied widely within the GPCR superfamily and to other small molecule-receptor protein systems.

Large-scale binding affinity calculations on commodity compute clouds

In recent years, it has become possible to calculate binding affinities of compounds bound to proteins via rapid, accurate, precise and reproducible free energy calculations. This is imperative in drug discovery as well as personalized medicine. This approach is based on molecular dynamics (MD) simulations and draws on sequence and structural information of the protein and compound concerned. Free energies are determined by ensemble averages of many MD replicas, each of which requires hundreds of cores and/or GPU accelerators, which are now available on commodity cloud computing platforms; there are also requirements for initial model building and subsequent data analysis stages. To automate the process, we have developed a workflow known as the binding affinity calculator. In this paper, we focus on the software infrastructure and interfaces that we have developed to automate the overall workflow and execute it on commodity cloud platforms, in order to reliably predict their binding affinities on time scales relevant to the domains of application, and illustrate its application to two free energy methods.

Markov state modeling reveals alternative unbinding pathways for peptide-MHC complexes

Peptide binding to major histocompatibility complexes (MHCs) is a central component of the immune system, and understanding the mechanism behind stable peptide-MHC binding will aid the development of immunotherapies. While MHC binding is mostly influenced by the identity of the so-called anchor positions of the peptide, secondary interactions from nonanchor positions are known to play a role in complex stability. However, current MHC-binding prediction methods lack an analysis of the major conformational states and might underestimate the impact of secondary interactions. In this work, we present an atomically detailed analysis of peptide-MHC binding that can reveal the contributions of any interaction toward stability. We propose a simulation framework that uses both umbrella sampling and adaptive sampling to generate a Markov state model (MSM) for a coronavirus-derived peptide (QFKDNVILL), bound to one of the most prevalent MHC receptors in humans (HLA-A24:02). While our model reaffirms the importance of the anchor positions of the peptide in establishing stable interactions, our model also reveals the underestimated importance of position 4 (p4), a nonanchor position. We confirmed our results by simulating the impact of specific peptide mutations and validated these predictions through competitive binding assays. By comparing the MSM of the wild-type system with those of the D4A and D4P mutations, our modeling reveals stark differences in unbinding pathways. The analysis presented here can be applied to any peptide-MHC complex of interest with a structural model as input, representing an important step toward comprehensive modeling of the MHC class I pathway.

Predicting Cross-Reactivity and Antigen Specificity of T Cell Receptors

Adaptive immune recognition is mediated by specific interactions between heterodimeric T cell receptors (TCRs) and their cognate peptide-MHC (pMHC) ligands, and the methods to accurately predict TCR:pMHC interaction would have profound clinical, therapeutic and pharmaceutical applications. Herein, we review recent developments in predicting cross-reactivity and antigen specificity of TCR recognition. We discuss current experimental and computational approaches to investigate cross-reactivity and antigen-specificity of TCRs and highlight how integrating kinetic, biophysical and structural features may offer valuable insights in modeling immunogenicity. We further underscore the close inter-relationship of these two interconnected notions and the need to investigate each in the light of the other for a better understanding of T cell responsiveness for the effective clinical applications.

Combining Three-Dimensional Modeling with Artificial Intelligence to Increase Specificity and Precision in Peptide-MHC Binding Predictions

The reliable prediction of the affinity of candidate peptides for the MHC is important for predicting their potential antigenicity and thus influences medical applications, such as decisions on their inclusion in T cell-based vaccines. In this study, we present a rapid, predictive computational approach that combines a popular, sequence-based artificial neural network method, NetMHCpan 4.0, with three-dimensional structural modeling. We find that the ensembles of bound peptide conformations generated by the programs MODELLER and Rosetta FlexPepDock are less variable in geometry for strong binders than for low-affinity peptides. In tests on 1271 peptide sequences for which the experimental dissociation constants of binding to the well-characterized murine MHC allele H-2Db are known, by applying thresholds for geometric fluctuations the structure-based approach in a standalone manner drastically improves the statistical specificity, reducing the number of false positives. Furthermore, filtering candidates generated with NetMHCpan 4.0 with the structure-based predictor led to an increase in the positive predictive value (PPV) of the peptides correctly predicted to bind very strongly (i.e., K d < 100 nM) from 40 to 52% (p = 0.027). The combined method also significantly improved the PPV when tested on five human alleles, including some with limited data for training. Overall, an average increase of 10% in the PPV was found over the standalone sequence-based method. The combined method should be useful in the rapid design of effective T cell-based vaccines.

Molecular Rules Underpinning Enhanced Affinity Binding of Human T Cell Receptors Engineered for Immunotherapy

Immuno-oncology approaches that utilize T cell receptors (TCRs) are becoming highly attractive because of their potential to target virtually all cellular proteins, including cancer-specific epitopes, via the recognition of peptide-human leukocyte antigen (pHLA) complexes presented at the cell surface. However, because natural TCRs generally recognize cancer-derived pHLAs with very weak affinities, efforts have been made to enhance their binding strength, in some cases by several million-fold. In this study, we investigated the mechanisms underpinning human TCR affinity enhancement by comparing the crystal structures of engineered enhanced affinity TCRs with those of their wild-type progenitors. Additionally, we performed molecular dynamics simulations to better understand the energetic mechanisms driving the affinity enhancements. These data demonstrate that supra-physiological binding affinities can be achieved without altering native TCR-pHLA binding modes via relatively subtle modifications to the interface contacts, often driven through the addition of buried hydrophobic residues. Individual energetic components of the TCR-pHLA interaction governing affinity enhancements were distinct and highly variable for each TCR, often resulting from additive, or knock-on, effects beyond the mutated residues. This comprehensive analysis of affinity-enhanced TCRs has important implications for the future rational design of engineered TCRs as efficacious and safe drugs for cancer treatment.

Large-Scale Structure-Based Prediction of Stable Peptide Binding to Class I HLAs Using Random Forests

Prediction of stable peptide binding to Class I HLAs is an important component for designing immunotherapies. While the best performing predictors are based on machine learning algorithms trained on peptide-HLA (pHLA) sequences, the use of structure for training predictors deserves further exploration. Given enough pHLA structures, a predictor based on the residue-residue interactions found in these structures has the potential to generalize for alleles with little or no experimental data. We have previously developed APE-Gen, a modeling approach able to produce pHLA structures in a scalable manner. In this work we use APE-Gen to model over 150,000 pHLA structures, the largest dataset of its kind, which were used to train a structure-based pan-allele model. We extract simple, homogenous features based on residue-residue distances between peptide and HLA, and build a random forest model for predicting stable pHLA binding. Our model achieves competitive AUROC values on leave-one-allele-out validation tests using significantly less data when compared to popular sequence-based methods. Additionally, our model offers an interpretation analysis that can reveal how the model composes the features to arrive at any given prediction. This interpretation analysis can be used to check if the model is in line with chemical intuition, and we showcase particular examples. Our work is a significant step toward using structure to achieve generalizable and more interpretable prediction for stable pHLA binding.

HLA-Arena: A Customizable Environment for the Structural Modeling and Analysis of Peptide-HLA Complexes for Cancer Immunotherapy

PURPOSE

HLA protein receptors play a key role in cellular immunity. They bind intracellular peptides and display them for recognition by T-cell lymphocytes. Because T-cell activation is partially driven by structural features of these peptide-HLA complexes, their structural modeling and analysis are becoming central components of cancer immunotherapy projects. Unfortunately, this kind of analysis is limited by the small number of experimentally determined structures of peptide-HLA complexes. Overcoming this limitation requires developing novel computational methods to model and analyze peptide-HLA structures.

METHODS

Here we describe a new platform for the structural modeling and analysis of peptide-HLA complexes, called HLA-Arena, which we have implemented using Jupyter Notebook and Docker. It is a customizable environment that facilitates the use of computational tools, such as APE-Gen and DINC, which we have previously applied to peptide-HLA complexes. By integrating other commonly used tools, such as MODELLER and MHCflurry, this environment includes support for diverse tasks in structural modeling, analysis, and visualization.

RESULTS

To illustrate the capabilities of HLA-Arena, we describe 3 example workflows applied to peptide-HLA complexes. Leveraging the strengths of our tools, DINC and APE-Gen, the first 2 workflows show how to perform geometry prediction for peptide-HLA complexes and structure-based binding prediction, respectively. The third workflow presents an example of large-scale virtual screening of peptides for multiple HLA alleles.

CONCLUSION

These workflows illustrate the potential benefits of HLA-Arena for the structural modeling and analysis of peptide-HLA complexes. Because HLA-Arena can easily be integrated within larger computational pipelines, we expect its potential impact to vastly increase. For instance, it could be used to conduct structural analyses for personalized cancer immunotherapy, neoantigen discovery, or vaccine development.

A comprehensive review and performance evaluation of bioinformatics tools for HLA class I peptide-binding prediction

Human leukocyte antigen class I (HLA-I) molecules are encoded by major histocompatibility complex (MHC) class I loci in humans. The binding and interaction between HLA-I molecules and intracellular peptides derived from a variety of proteolytic mechanisms play a crucial role in subsequent T-cell recognition of target cells and the specificity of the immune response. In this context, tools that predict the likelihood for a peptide to bind to specific HLA class I allotypes are important for selecting the most promising antigenic targets for immunotherapy. In this article, we comprehensively review a variety of currently available tools for predicting the binding of peptides to a selection of HLA-I allomorphs. Specifically, we compare their calculation methods for the prediction score, employed algorithms, evaluation strategies and software functionalities. In addition, we have evaluated the prediction performance of the reviewed tools based on an independent validation data set, containing 21 101 experimentally verified ligands across 19 HLA-I allotypes. The benchmarking results show that MixMHCpred 2.0.1 achieves the best performance for predicting peptides binding to most of the HLA-I allomorphs studied, while NetMHCpan 4.0 and NetMHCcons 1.1 outperform the other machine learning-based and consensus-based tools, respectively. Importantly, it should be noted that a peptide predicted with a higher binding score for a specific HLA allotype does not necessarily imply it will be immunogenic. That said, peptide-binding predictors are still very useful in that they can help to significantly reduce the large number of epitope candidates that need to be experimentally verified. Several other factors, including susceptibility to proteasome cleavage, peptide transport into the endoplasmic reticulum and T-cell receptor repertoire, also contribute to the immunogenicity of peptide antigens, and some of them can be considered by some predictors. Therefore, integrating features derived from these additional factors together with HLA-binding properties by using machine-learning algorithms may increase the prediction accuracy of immunogenic peptides. As such, we anticipate that this review and benchmarking survey will assist researchers in selecting appropriate prediction tools that best suit their purposes and provide useful guidelines for the development of improved antigen predictors in the future.

On Restraints in End-Point Protein-Ligand Binding Free Energy Calculations

The impact of harmonic restraints on protein heavy atoms and ligand atoms on end-point free energy calculations is systematically characterized for 54 protein-ligand complexes. We observe that stronger restraints reduce the equilibration time and statistical inefficiency, suppress conformational sampling, influence correlation with experiment, and monotonically decrease the estimated loss of entropy upon binding, leading to stronger estimated binding free energies in most systems. A statistical estimator that reweights for the biasing potential and includes data prior to the estimated equilibration time has the highest correlation with experiment. A spring constant of 20 cal mol-1 Å-2 maintains a near-native energy landscape and suppresses artifactual energy minima while minimally limiting thermal fluctuations about the crystal structure. © 2019 Wiley Periodicals, Inc.

Molecular Basis for Metabolic Regioselectivity and Mechanism of Cytochrome P450s toward Carcinogenic 4-(Methylnitrosamino)-(3-pyridyl)-1-butanone

As an abundantly present tobacco component, carcinogen 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanone (NNK) has also been detected in atmospheric particulate matter, suggesting the ineluctable exposure risk of this contaminant. NNK metabolic activation by cytochrome P450 enzymes (CYPs) is a prerequisite to exerting its genotoxicity, but the metabolic regioselectivity and mechanism are still unknown. Here the binding feature and regioselectivity of CYPs 1A1, 1A2, 2A6, 2A13, 2B6, and 3A4 toward NNK are unraveled through molecular docking and molecular dynamics (MD) simulations. Binding mode analyses reveal that 1A2 and 2B6 have definite preferences for NNK α-methyl hydroxylation, while the other four CYPs preferentially catalyze α-methylene hydroxylation. The binding affinities between NNK and CYPs evaluated by the binding free energies follow the order 2A13 > 2B6 > 1A2 > 2A6 > 1A1 > 3A4. Density functional theory (DFT) calculations are further performed to characterize the mechanism of NNK biotransformation. Results show that the α-hydroxyNNK generated from α-hydroxylation may undergo nonenzymatic decomposition to form genotoxic diazohydroxide and aldehyde, and further oxidation by P450 to yield nitrosamide, which mainly contributes to NNK toxification capacity. Meanwhile the pyridine N-oxidation and denitrosation of C_α-radical intermediate play an important role in detoxifying NNK. Overall, the present study provides the molecular basis for CYP-catalyzed regioselectivity and mechanism of NNK biotransformation, which can enable the identification of metabolites for assessing the health risk of individual NNK exposure.

Exploring the different states of wild-type T-cell receptor and mutant conformational changes towards understanding the antigen recognition

Recognition of proteolytic peptide fragments presented by major histocompatibility complex (MHC) on target cells by T-cell receptor (TCR) is among the most important interactions in the adaptive immune system. Several computational studies have been performed to investigate conformational and dynamical properties of TCRs for enhanced immunogenicity. Here, we present the large-scale molecular dynamics (MD) simulation studies of the two comprehensive systems consisting of the wild-type and mutant IG4 TCR in complex with the tumor epitope NY-ESO peptide (SLLMWITQC) and analyzed for mapping conformational changes of TCR in the states prior to antigen binding, upon antigen binding and after the antigen was released. All of the simulations were performed with different states of TCRs for each 1000 ns of simulation time, providing six simulations for time duration of 6000 ns (6µs). We show that rather than undergoing most critical conformational changes upon antigen binding, the high proportion of complementarity-determining region (CDR) loops change by comparatively small amount. The hypervariable CDRα3 and CDRβ3 loops showed significant structural changes. Interestingly, the TCR β chain loops showed the least changes, which is reliable with recent implications that β domain of TCR may propel antigen interaction. The mutant shows higher rigidity than wild-type even in released state; expose an induced fit mechanism occurring from the re-structuring of CDRα3 loop and can allow enhanced binding affinity of the peptide antigen. Additionally, we show that CDRα3 loop and peptide contacts are an adaptive feature of affinity enhanced mutant TCR.Communicated by Ramaswamy H. Sarma.

Accuracy and Precision of Alchemical Relative Free Energy Predictions with and without Replica-Exchange

A systematic and statistically robust protocol is applied for the evaluation of free energy calculations with and without replica-exchange. The protocol is based on ensemble averaging to generate accurate assessments of the uncertainties in the predictions. Comparison is made between FEP+ and TIES-free energy perturbation and thermodynamic integration with enhanced sampling-the latter with and without the so-called "enhanced sampling" based on replica-exchange protocols. Standard TIES performs best for a reference set of targets and compounds; no benefits accrue from replica-exchange methods. Evaluation of FEP+ and TIES with REST-replica-exchange with solute tempering-reveals a systematic and significant underestimation of free energy differences in FEP+, which becomes increasingly large for long duration simulations, is confirmed by extensive analysis of previous publications, and raises a number of questions pertaining to the accuracy of the predictions with the REST technique not hitherto discussed.