Alanine-scanning mutagenesis has proven to be a key experimental tool to identify the residues that are most important for mediating interactions between a protein and a ligand (which could be a small molecule, a nucleic acid, a different protein, or a different region within the same protein, etc.). However, mutating individual / combinations of amino acids to alanine experimentally is both time-consuming and costly. Consequently, several computational alanine scanning (CAS) programs have been developed that enable users to identify hotspot residues that make the greatest contribution to the interaction of interest in silico.

In two recent papers ^{[1, 2]}, the Sessions and Woolfson groups at the University of Bristol, in collaboration with the Wilson group at the University of Leeds, introduce BUDE Alanine Scan (BAlaS), a command-line tool to run CAS within the BUDE force field ^{[3, 4]}. BAlaS enables high-throughput CAS of both individual X-ray crystal structures plus ensembles of structures obtained from NMR / MD analysis. In Ibarra et al. (2019), the authors introduce the BAlaS software, and benchmark its performance relative to a range of pre-existing CAS software, firstly on the SKEMPI database (a database of 3047 binding free energy changes upon mutation, collated from the literature, for PPIs with known structure), and secondly on newly acquired alanine-scanning mutagenesis data collected for three diverse protein-protein interactions. They find that BAlaS predicts the experimentally identified hotspot residues with comparable accuracy to the best-performing CAS software, whilst running considerably faster. However, the authors also find that averaging the results from all CAS software tested achieves better results than any of the individual programs, and hence recommend CAS users adopt such a meta-scoring approach.

Figure 1: Comparison of experimentally measured vs. computationally predicted (via 6 different CAS software packages) DDG values for individual alanine mutations of interface residues in three diverse protein-protein interactions: a) NOXA-B/MCL-1 (PDB model 1); b) SIMS/SUMO (PDB model 1); and c) GKAP/SHANK-PDZ (chains A and C).

In Wood et al. (2020), the authors introduce an interactive web application to run BAlaS (https://pragmaticproteindesign.bio.ed.ac.uk/balas/). This interface is simple and intuitive, allowing non-expert users to easily run CAS on their own protein-ligand interactions of interest. In addition to allowing mutation of individual residues to alanine, the interface also allows the user to simultaneously mutate constellations of residues (selected either manually or via all permutations of a specified constellation size within a selected subset of residues) to alanine in order to identify “hotspot patches” of residues within the ligand. The capability to perform mutations to residues other than alanine, available in the command-line tool, will be introduced in a future version of the web application.

Figure 2: The BAlaS web interface. Here DDG values are calculated for residues in the transactivation domain of p53 (the “ligand”) bound to MDM2 (the “receptor”) (PDB structure 1YCR).

References

1. A. A. Ibarra et al., Predicting and experimentally validating hot-spot residues at protein–protein interfaces. ACS Chem Biol 14, 2252-2263 (2019).
2. C. W. Wood et al., BAlaS: fast, interactive and accessible computational alanine-scanning using BudeAlaScan. Bioinformatics (2020).
3. S. McIntosh-Smith, T. Wilson, A. Á. Ibarra, J. Crisp, R. B. Sessions, Benchmarking energy efficiency, power costs and carbon emissions on heterogeneous systems. The Computer Journal 55, 192-205 (2012).
4. S. McIntosh-Smith, J. Price, R. B. Sessions, A. A. Ibarra, High performance in silico virtual drug screening on many-core processors. The International Journal of High Performance Computing Applications 29, 119-134 (2015).

Intricate molecular choreography separates the bound and unbound states of proteins, resulting in slow binding kinetics that make simulating these events difficult. Interactive molecular dynamics (iMD) allows simulations to be driven to a desired state using a combination of human chemical and spatial intuition and current advances present virtual reality (VR) as a novel strategy for simulation biasing. Previous work shows VR interfaces to iMD (iMD-VR) are sufficiently performant for basic molecular manipulation tasks, but this begs the question of whether these interfaces can be used for more complex tasks, namely correctly establishing bound poses between proteins and ligands.

Our iMD-VR framework for interactively chaperoning drug unbinding and rebinding events was tested on a cohort of non-expert users (n=10), defined as being unfamiliar with the iMD-VR interface. Three example protein-ligand systems were chosen, and participants were asked to undock and redock a ligand from these proteins twice, first with the aid of a trace representation of the ligand in the correct position to guide them, and once again with no trace guide present. Below is a representation of two people in VR undocking a ligand from a protein, and also the three protein-ligand systems used in the study, with key interactions highlighted.

We show that it is possible to unbind and rebind ligands from three protein targets in a simulation timescale totalling less than 50ps, simulated at a rate of 4.5 ps/min of real time. Where non-expert users had trace atoms showing them the correct pose, all users were able to establish a docking pose within 2Å of the starting structure (1Å for two out of the three tasks). Where no trace atoms were present, binding poses understandably had higher variation, however, participants were still able to get within the same range of RMSD for all three systems. These results were achieved within a single hour-long training session with each participant and are presented below.

Our results show that, even with very limited experience, iMD-VR is performant enough to allow (a) unbinding of a ligand from a protein binding pocket and (b) re-establishment of the original binding pose, demonstrating iMD-VR as a useful tool for sampling states where a ligand is both bound and unbound from a protein.

Helen M. Deeks, Rebecca K. Walters, Stephanie R. Hare, Michael B. O’Connor, Adrian J. Mulholland and David R. Glowacki

PLOSONE, 2020

doi: to be added very soon!

In this recent publication included in the ‘Women in Computational Chemistry’ special issue of the Journal of Chemical Information and Modeling, we investigate the dependence of predicted reaction energetics on the choice of density functional and QM region size used in QM/MM calculations. DFT methods potentially offer a good combination of accuracy and computational cost but suffer from some well-known limitations and are not systematically improvable. The Claisen rearrangement of chorismate to prephenate, a simple unimolecular reaction catalysed by the enzyme chorismate mutase, was used as a model system to investigate how these choices for QM/MM calculations impact the predicted energetics of reaction.

There are many different density functionals available and the best density functional to choose for any application is not obvious. Commonly used density functionals (e.g. B3LYP, BH&HLYP, MO6-2X) predict barrier heights for the reaction in the enzyme that differ by as much 13 kcal/mol. When different density functionals give different results for the same system, which result should be preferred? Ab initio methods are potentially highly accurate and are systematically improvable but often restricted by their computational cost. Projector-based embedding allows the incorporation of ab initio methods into QM/MM calculations at reasonable computational cost as a small number of atoms can be selected for treatment at the highest levels. Embedding SCS-MP2 in DFT in these QM/MM calculations reduced the spread in predicted barrier height from 13 kcal/mol to 0.3 kcal/mol, essentially eliminating the dependence on the density functional, as shown below (top: standard DFT, bottom: SCS-MP2-in-DFT).

The optimum size of the QM region in QM/MM calculations is a matter of much debate in the literature. The effect of the size of the QM region was tested in solution (by adding additional water molecules to the QM region) and in the enzyme (by adding additional arginine side chains) showing that the difference in barrier heights predicted by common density functionals for any size of the QM region is significantly larger than the change in barrier caused by increasing the size of the QM region.

Ranaghan, K. E.,  Shchepanovska, D., Bennie, S. J.,  Lawan, N., Macrae, S. J., Zurek, J., Manby, F. R. and Mulholland, A. J.
J. Chem. Inf. Model., 2019
DOI: 10.1021/acs.jcim.8b00940

Stefano Artin Serapian and Marc W. van der Kamp

In this recent ACS Catalysis publication, Stefano Serapian and Marc van der Kamp used a variety of computational tools to shed light on an enticing problem in biocatalysis that has proven very difficult to solve experimentally.

Our enzyme of interest—actinorhodin ketoreductase (actKR)—is found in the soil bacterium Streptomyces coelicolor, where it is normally implicated in the biosynthesis of the antibiotic actinorhodin.  In addition to actKR’s natural scope, as is the case with other ketoreductases, some of its re-engineered variants are highly attractive to synthetic chemists by virtue of their high stereoselectivity in reducing small-molecule achiral ketones to chiral alcohols.

Experimental work on actKR featuring the small model substrate trans-1-decalone (a bicyclic aliphatic ketone) and similarly-sized chiral alcohols suggests that whereas the wild-type enzyme is mildly S-selective, some variants (e.g. the P94L mutant) were entirely S- selective, whereas others (e.g. the V151L mutant) only exhibited R- selectivity.

Featuring classical MD, MM/PBSA and hybrid QM/MM MD with umbrella sampling, our study successfully unravels the causes of such remarkable behaviour in wild-type, P94L, and V151L actKR towards trans-1-decalone.  Explicitly examining both enantiomers of this naturally racemic substrate (something difficult to achieve in vitro), we conclude that changes in stereocontrol across actKR variants can be dictated by a subtle interplay of different causes (including reaction barrier height and accessibility of reactive poses). Interestingly, however, which factor is dominant in conferring stereoselectivity differs per variant.

The protocols we have used were chosen such that they (1) require input of the WT structure only; (2) use relatively limited computational resources (short simulations and semiempirical QM treatment); and (3) can be automated. Our study is thereby a good example of how computational biochemistry can become a practical, useful and efficient tool in biocatalyst engineering, offering perspectives that might otherwise be difficult to explore in vitro.