Automated Absolute Binding Free Energy Calculation Workflow for Drug Discovery.

Absolute binding free energies play a crucial role in drug development, particularly as part of the lead discovery process. In recent work, we showed how in silico predictions directly could support drug development by ranking and recommending favorable ideas over unfavorable ones. Here, we demonstrate a Python workflow that enables the calculation of ABFEs with minimal manual input effort, such as the receptor PDB and ligand SDF files, and outputs a .tsv file containing the ranked ligands and their corresponding binding free energies. The implementation uses Snakemake to structure and control the execution of tasks, allowing for dynamic control of parameters and execution patterns. We provide an example of a benchmark system that demonstrates the effectiveness of the automated workflow.

Proton Control of Transitions in an Amino Acid Transporter.

Amino acid transport into the cell is often coupled to the proton electrochemical gradient, as found in the solute carrier 36 family of proton-coupled amino acid transporters. Although no structure of a human proton-coupled amino acid transporter exists, the crystal structure of a related homolog from bacteria, GkApcT, has recently been solved in an inward-occluded state and allows an opportunity to examine how protons are coupled to amino acid transport. Our working hypothesis is that release of the amino acid substrate is facilitated by the deprotonation of a key glutamate residue (E115) located at the bottom of the binding pocket, which forms part of the intracellular gate, allowing the protein to transition from an inward-occluded to an inward-open conformation. During unbiased molecular dynamics simulations, we observed a transition from the inward-occluded state captured in the crystal structure to a much more open state, which we consider likely to be representative of the inward-open state associated with substrate release. To explore this and the role of protons in these transitions, we have used umbrella sampling to demonstrate that the transition from inward occluded to inward open is more energetically favorable when E115 is deprotonated. That E115 is likely to be protonated in the inward-occluded state and deprotonated in the inward-open state is further confirmed via the use of absolute binding free energies. Finally, we also show, via the use of absolute binding free energy calculations, that the affinity of the protein for alanine is very similar regardless of either the conformational state or the protonation of E115, presumably reflecting the fact that all the key interactions are deep within the binding cavity. Together, our results give a detailed picture of the role of protons in driving one of the major transitions in this transporter.

Ring Puckering Landscapes of Glycosaminoglycan-Related Monosaccharides from Molecular Dynamics Simulations.

The conformational flexibility of the glycosaminoglycans (GAGs) is known to be key in their binding and biological function, for example in regulating coagulation and cell growth. In this work, we employ enhanced sampling molecular dynamics simulations to probe the ring conformations of GAG-related monosaccharides, including a range of acetylated and sulfated GAG residues. We first perform unbiased MD simulations of glucose anomers and the epimers glucuronate and iduronate. These calculations indicate that in some cases, an excess of 15 µs is required for adequate sampling of ring pucker due to the high energy barriers between states. However, by applying our recently developed msesMD simulation method (multidimensional swarm-enhanced sampling molecular dynamics), we were able to quantitatively and rapidly reproduce these ring pucker landscapes. From msesMD simulations, the puckering free energy profiles were then compared for 15 further monosaccharides related to GAGs; this includes to our knowledge the first simulation study of sulfation effects on ß-GalNAc ring puckering. For the force field employed, we find that in general the calculated pucker free energy profiles for sulfated sugars were similar to the corresponding unsulfated profiles. This accords with recent experimental studies suggesting that variation in ring pucker of sulfated GAG residues is primarily dictated by interactions with surrounding residues rather than by intrinsic conformational preference. As an exception to this, however, we predict that 4-O-sulfation of ß-GalNAc leads to reduced ring rigidity, with a significant lowering in energy of the 1C4 ring conformation; this observation may have implications for understanding the structural basis of the biological function of ß-GalNAc-containing glycosaminoglycans such as dermatan sulfate.

Energetic Fingerprinting of Ligand Binding to Paralogous Proteins: The Case of the Apoptotic Pathway.

Networks of biological molecules are key to cellular function, governing processes ranging from signal cascade propagation to metabolic pathway regulation. Genetic duplication processes give rise to sets of regulatory proteins that have evolved from a common ancestor, leading to interactomes whose dysregulation is often associated with disease. A better understanding of the determinants of specificity at interfaces shared by functionally related proteins is crucial to the rational design of novel pharmacotherapeutic agents. To this end, a comprehensive data set of drug and drug-like binders was assembled for the Bcl-xL and Bcl-2 antiapoptotic proteins-archetypal examples of regulatory systems governed by evolutionarily conserved protein-protein interactions. These were first used to derive a two-dimensional quantitative structure-activity relationship (2D QSAR) model, predicting ligand specificity for these homologous proteins. The strengths and weaknesses of high-throughput 2D QSAR were then compared and contrasted to those of theoretically rigorous thermodynamic integration calculations performed on 14 complexes of Bcl-xL-specific, Bcl-2-specific, and potent dual binders bound to the Bcl-xL and Bcl-2 proteins. We demonstrate that free energy calculations provide an added layer of essential information, which traditional QSAR cannot capture. Moreover, we show that protein energetic responses to different ligands, expressed as per-residue energy values, can be used to fingerprint the protein-ligand interaction, extending the framework of four-dimensional molecular dynamics/quantitative structure-activity relationships (4D-MD/QSAR) toward the facilitation of future drug design strategies.

Free Energy Calculations using a Swarm-Enhanced Sampling Molecular Dynamics Approach.

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm-enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT-TI) method. Free energy changes for transitions computed by using IT-TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm-enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.

Kartograf: A Geometrically Accurate Atom Mapper for Hybrid-Topology Relative Free Energy Calculations.

Relative binding free energy (RBFE) calculations have emerged as a powerful tool that supports ligand optimization in drug discovery. Despite many successes, the use of RBFEs can often be limited by automation problems, in particular, the setup of such calculations. Atom mapping algorithms are an essential component in setting up automatic large-scale hybrid-topology RBFE calculation campaigns. Traditional algorithms typically employ a 2D subgraph isomorphism solver (SIS) in order to estimate the maximum common substructure. SIS-based approaches can be limited by time-intensive operations and issues with capturing geometry-linked chemical properties, potentially leading to suboptimal solutions. To overcome these limitations, we have developed Kartograf, a geometric-graph-based algorithm that uses primarily the 3D coordinates of atoms to find a mapping between two ligands. In free energy approaches, the ligand conformations are usually derived from docking or other previous modeling approaches, giving the coordinates a certain importance. By considering the spatial relationships between atoms related to the molecule coordinates, our algorithm bypasses the computationally complex subgraph matching of SIS-based approaches and reduces the problem to a much simpler bipartite graph matching problem. Moreover, Kartograf effectively circumvents typical mapping issues induced by molecule symmetry and stereoisomerism, making it a more robust approach for atom mapping from a geometric perspective. To validate our method, we calculated mappings with our novel approach using a diverse set of small molecules and used the mappings in relative hydration and binding free energy calculations. The comparison with two SIS-based algorithms showed that Kartograf offers a fast alternative approach. The code for Kartograf is freely available on GitHub (https://github.com/OpenFreeEnergy/kartograf). While developed for the OpenFE ecosystem, Kartograf can also be utilized as a standalone Python package.

Evaluating the use of absolute binding free energy in the fragment optimisation process.

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman's r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely Nwat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

Atomistic mechanisms of human TRPA1 activation by electrophile irritants through molecular dynamics simulation and mutual information analysis.

The ion channel TRPA1 is a promiscuous chemosensor, with reported response to a wide spectrum of noxious electrophilic irritants, as well as cold, heat, and mechanosensation. It is also implicated in the inception of itch and pain and has hence been investigated as a drug target for novel analgesics. The mechanism of electrophilic activation for TRPA1 is therefore of broad interest. TRPA1 structures with the pore in both open and closed states have recently been published as well as covalent binding modes for electrophile agonists. However, the detailed mechanism of coupling between electrophile binding sites and the pore remains speculative. In addition, while two different cysteine residues (C621 and C665) have been identified as critical for electrophile bonding and activation, the bound geometry has only been resolved at C621. Here, we use molecular dynamics simulations of TRPA1 in both pore-open and pore-closed states to explore the allosteric link between the electrophile binding sites and pore stability. Our simulations reveal that an open pore is structurally stable in the presence of open 'pockets' in the C621/C665 region, but rapidly collapses and closes when these pockets are shut. Binding of electrophiles at either C621 or C665 provides stabilisation of the pore-open state, but molecules bound at C665 are shown to be able to rotate in and out of the pocket, allowing for immediate stabilisation of transient open states. Finally, mutual information analysis of trajectories reveals an informational path linking the electrophile binding site pocket to the pore via the voltage-sensing-like domain, giving a detailed insight into the how the pore is stabilized in the open state.

The structure of nontypeable Haemophilus influenzae SapA in a closed conformation reveals a constricted ligand-binding cavity and a novel RNA binding motif.

Nontypeable Haemophilus influenzae (NTHi) is a significant pathogen in respiratory disease and otitis media. Important for NTHi survival, colonization and persistence in vivo is the Sap (sensitivity to antimicrobial peptides) ABC transporter system. Current models propose a direct role for Sap in heme and antimicrobial peptide (AMP) transport. Here, the crystal structure of SapA, the periplasmic component of Sap, in a closed, ligand bound conformation, is presented. Phylogenetic and cavity volume analysis predicts that the small, hydrophobic SapA central ligand binding cavity is most likely occupied by a hydrophobic di- or tri- peptide. The cavity is of insufficient volume to accommodate heme or folded AMPs. Crystal structures of SapA have identified surface interactions with heme and dsRNA. Heme binds SapA weakly (Kd 282 µM) through a surface exposed histidine, while the dsRNA is coordinated via residues which constitute part of a conserved motif (estimated Kd 4.4 µM). The RNA affinity falls within the range observed for characterized RNA/protein complexes. Overall, we describe in molecular-detail the interactions of SapA with heme and dsRNA and propose a role for SapA in the transport of di- or tri-peptides.

Identification of Rare Lewis Oligosaccharide Conformers in Aqueous Solution Using Enhanced Sampling Molecular Dynamics.

Determining the conformations accessible to carbohydrate ligands in aqueous solution is important for understanding their biological action. In this work, we evaluate the conformational free-energy surfaces of Lewis oligosaccharides in explicit aqueous solvent using a multidimensional variant of the swarm-enhanced sampling molecular dynamics (msesMD) method; we compare with multi-microsecond unbiased MD simulations, umbrella sampling, and accelerated MD approaches. For the sialyl Lewis A tetrasaccharide, msesMD simulations in aqueous solution predict conformer landscapes in general agreement with the other biased methods and with triplicate unbiased 10 µs trajectories; these simulations find a predominance of closed conformer and a range of low-occupancy open forms. The msesMD simulations also suggest closed-to-open transitions in the tetrasaccharide are facilitated by changes in ring puckering of its GlcNAc residue away from the 4C1 form, in line with previous work. For sialyl Lewis X tetrasaccharide, msesMD simulations predict a minor population of an open form in solution corresponding to a rare lectin-bound pose observed crystallographically. Overall, from comparison with biased MD calculations, we find that triplicate 10 µs unbiased MD simulations may not be enough to fully sample glycan conformations in aqueous solution. However, the computational efficiency and intuitive approach of the msesMD method suggest potential for its application in glycomics as a tool for analysis of oligosaccharide conformation.

'Dual' peptidyl-oligonucleotide conjugates: Role of conformational flexibility in catalytic cleavage of RNA.

Traditional therapeutic interventions against abnormal gene expression in disease states at the level of expressed proteins are becoming increasingly difficult due to poor selectivity, off-target effects and associated toxicity. Upstream catalytic targeting of specific RNA sequences offers an alternative platform for drug discovery to achieve more potent and selective treatment through antisense interference with disease-relevant RNAs. We report a novel class of catalytic biomaterials, comprising amphipathic RNA-cleaving peptides placed between two RNA recognition motifs, here demonstrated to target the TΨC loop and 3'- acceptor stem of tRNAPhe. These unique peptidyl-oligonucleotide 'dual' conjugates (DCs) were created by phosphoramidate or thiol-maleimide conjugation chemistry of a TΨC-targeting oligonucleotide to the N-terminus of the amphipathic peptide sequence, followed by amide coupling of a 3'-acceptor stem-targeting oligonucleotide to the free C-terminal carboxylic acid functionality of the same peptide. Hybridization of the DCs bearing two spatially-separated recognition motifs with the target tRNAPhe placed the peptide adjacent to a single-stranded RNA region and promoted cleavage within the 'action radius' of the catalytic peptide. Up to 100% cleavage of the target tRNAPhe was achieved by the best candidate (i.e. DC6) within 4 h, when conformational flexibility was introduced into the linker regions between the peptide and oligonucleotide components. This study provides the strong position for future development of highly selective RNA-targeting agents that can potentially be used for disease-selective treatment at the level of messenger, micro, and genomic viral RNA.