Biotherapeutic optimization, whether to improve general properties or to engineer specific attributes, is a time-consuming process with uncertain outcomes. Conversely, Consensus Protein Design has been shown to be a viable approach to enhance protein stability while retaining function. In adapting this method for a more limited number of protein sequences, we studied 21 consensus single-point variants from eight publicly available CD3 binding sequences with high similarity but diverse biophysical and pharmacological properties. All single-point consensus variants retained CD3 binding and performed similarly in cell-based functional assays. Using Ridge regression analysis, we identified the variants and sequence positions with overall beneficial effects on developability attributes of the CD3 binders. A second round of sequence generation that combined these substitutions into a single molecule yielded a unique CD3 binder with globally optimized developability attributes. In this first application to therapeutic antibodies, adapted Consensus Protein Design was found to be highly beneficial within lead optimization, conserving resources and minimizing iterations. Future implementations of this general strategy may help accelerate drug discovery and improve success rates in bringing novel biotherapeutics to market.

Supplementary Information: "Evaluating the use of absolute binding free energy in the fragment optimization process" Included are the ABFE raw free energy samples for multiple replicates (labelled by `run` number) of 19 ligands to bound MCL-1. These ligands are originally detailed by Friberg et al. (https://doi.org/10.1021/jm301448p). All samples are provided as a set of `.xvg` files as generated by GROMACS 2021 (https://doi.org/10.5281/zenodo.5849961). The `.xvg` files are labelled as dhdl.N.xvg where N represents the λ state the free energy values were sampled from. The `.xvg` files contain both ΔH/Δλ and ΔH values, please see the header of each files for more information. Samples detailing the partial decoupling of the ligand from the protein-ligand complex are contained within the `complex` folder. These consist of an orientational restraint addition step (found within the `restraints-xvg` folders), charge annihilation step (found within the `coul-xvg` folders), and Van der Waals decoupling step (found within the `vdw-xvg` folders). Samples detailing the partial decoupling of the ligand from solvent are contained within the `ligand` folder and consist of a charge annihilation step (found within the individual `coul-xvg` folders) and a Van der Waals decoupling step (found within the individual `vdw-xvg` folders).

Key to the fragment optimization process 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 optimization 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 optimization decisions can be supported by the ABFE calculations. Cases that were not supported were often limited by large uncertainties in the free energy estimates, however generally the right direction in ΔΔG is still observed. Comparing against cheaper endpoint methods, namely Nwat-MM/GBSA, we find that ABFEs offer better outcomes in ranking binders, improving correlation metrics, although a similar confidence in retrospective synthetic decisions is achieved. Our results indicate that ABFE calculations are currently at the level of accuracy that can be usefully employed to gauge which fragment elaborations are likely to offer the best gains in affinity.

Supplementary Information: "Evaluating the use of absolute binding free energy in the fragment optimization process" Provided here are the various scripts, input files, and results necessary to reproduce the outcomes of the above mentioned publication. Please see the provided README.md files for further information on the contents of this dataset.

Computational methods assisting drug discovery and development are routine in the pharmaceutical industry. Digital recording of ADMET assays has provided a rich source of data for development of predictive models. Despite the accumulation of data and the public availability of advanced modeling algorithms, the utility of prediction in ADMET research is not clear. Here, we present a critical evaluation of the relationships between data volume, modeling algorithm, chemical representation and grouping, and temporal aspect (time sequence of assays) using an in-house ADMET database. We find no large difference in prediction algorithms nor any systemic and substantial gain from increasingly large datasets. Temporal-based data enlargement led to performance improvement in only in a limited number of assays, and with fractional improvement at best. Assays that are well-, intermediately-, or poorly-suited for ADMET predictions and reasons for such behavior are systematically identified, generating realistic expectations for areas in which computational models can be used to guide decision making in molecular design and development.

In the past decade, the pharmaceutical industry has paid closer attention to covalent drugs. Differently from standard noncovalent drugs, these compounds can exhibit peculiar properties, such as higher potency or longer duration of target inhibition with a potentially lower dosage. These properties are mainly driven by the reactive functional group present in the compound, the so-called warhead that forms a covalent bond with a specific nucleophilic amino-acid on the target. In this work, we report the possibility to combine ab initio activation energies with machine-learning to estimate covalent compound intrinsic reactivity. The idea behind this approach is to have a precise estimation of the transition state barriers, and thus of the compound reactivity, but with the speed of a machine-learning algorithm. We call this method "BIreactive". Here, we demonstrate this approach on acrylamides and 2-chloroacetamides, two warhead classes that possess different reaction mechanisms. In combination with our recently implemented truncation algorithm, we also demonstrate the possibility to use BIreactive not only for fragments but also for lead-like molecules. The generic nature of this approach allows also the extension to several other warheads. The combination of these factors makes BIreactive a valuable tool for the covalent drug discovery process in a pharmaceutical context.

Fragment-based drug discovery (FBDD) permits efficient sampling of the vast chemical space for hit identification. Libraries are screened biophysically and fragment:protein co-structures are determined by X-ray crystallography. In parallel, computational methods can derive pharmacophore models or screen virtual libraries. We screened 15 very small fragments (VSFs) (HA ≤ 11) computationally, using site identification by ligand competitive saturation (SILCS), and experimentally, by X-ray crystallography, to map potential interaction sites on the FKBP51 FK1 domain. We identified three hot spots and obtained 6 X-ray co-structures, giving a hit rate of 40%. SILCS FragMaps overlapped with X-ray structures. The compounds had millimolar affinities as determined by 15N HSQC NMR. VSFs identified the same interactions as known FK1 binder and provide new chemical starting points. We propose a hybrid screening strategy starting with SILCS, followed by a pharmacophore-derived X-ray screen and 15N HSQC NMR based KD determination to rapidly identify hits and their binding poses.

Predicting the costructure of small-molecule ligands and their respective target proteins has been a long-standing problem in drug discovery. For weak binding compounds typically identified in fragment-based screening (FBS) campaigns, determination of the correct binding site and correct binding mode is usually done experimentally via X-ray crystallography. For many targets of pharmaceutical interest, however, establishing an X-ray system which allows for sufficient throughput to support a drug discovery project is not possible. In this case, exploration of fragment hits becomes a very laborious and consequently slow process with the generation of protein/ligand cocrystal structures as the bottleneck of the entire process. In this work, we introduce a computational method which is able to reliably predict binding sites and binding modes of fragment-like small molecules using solely the structure of the apoprotein and the ligand's chemical structure as input information. The method is based on molecular dynamics simulations and Markov-state models and can be run as a fully automated protocol requiring minimal human intervention. We describe the application of the method to a representative subset of different target classes and fragments from historical FBS efforts at Boehringer Ingelheim and discuss its potential integration into the overall fragment-based drug discovery workflow.

Ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design. However, their broad uptake and impact is held back by the notoriously complex setup of the calculations. Only a few tools other than the free energy perturbation approach by Schrödinger Inc. (referred to as FEP+) currently enable end-to-end application. Here, we present for the first time an approach based on the open-source software pmx that allows to easily set up and run alchemical calculations for diverse sets of small molecules using the GROMACS MD engine. The method relies on theoretically rigorous non-equilibrium thermodynamic integration (TI) foundations, and its flexibility allows calculations with multiple force fields. In this study, results from the Amber and Charmm force fields were combined to yield a consensus outcome performing on par with the commercial FEP+ approach. A large dataset of 482 perturbations from 13 different protein-ligand datasets led to an average unsigned error (AUE) of 3.64 ± 0.14 kJ mol-1, equivalent to Schrödinger's FEP+ AUE of 3.66 ± 0.14 kJ mol-1. For the first time, a setup is presented for overall high precision and high accuracy relative protein-ligand alchemical free energy calculations based on open-source software.

The target residence time (RT) for a given ligand is one of the important parameters that have to be optimized during drug design. It is well established that shielding the receptor-ligand hydrogen bond (H-bond) interactions from water has been one of the factors in increasing ligand RT. Building on this foundation, here we report that shielding an intra-protein H-bond, which confers rigidity to the binding pocket and which is not directly involved in drug-receptor interactions, can strongly influence RT for CCR2 antagonists. Based on our recently solved CCR2 structure with MK-0812 and molecular dynamics (MD) simulations, we show that the RT for this and structurally related ligands is directly dependent on the shielding of the Tyr120-Glu291 H-bond from the water. If solvated this H-bond is often broken, making the binding pocket flexible and leading to shorter RT.

Background and Purpose The bronchodilator tiotropium binds not only to its main binding site on the M-3 muscarinic receptor but also to an allosteric site. Here, we have investigated the functional relevance of this allosteric binding and the potential contribution of this behaviour to interactions with long-acting beta-adrenoceptor agonists, as combination therapy with anticholinergic agents and beta-adrenoceptor agonists improves lung function in chronic obstructive pulmonary disease. Experimental Approach ACh, tiotropium, and atropine binding to M-3 receptors were modelled using molecular dynamics simulations. Contractions of bovine and human tracheal smooth muscle strips were studied. Key Results Molecular dynamics simulation revealed extracellular vestibule binding of tiotropium, and not atropine, to M-3 receptors as a secondary low affinity binding site, preventing ACh entry into the orthosteric binding pocket. This resulted in a low (allosteric binding) and high (orthosteric binding) functional affinity of tiotropium in protecting against methacholine-induced contractions of airway smooth muscle, which was not observed for atropine and glycopyrrolate. Moreover, antagonism by tiotropium was insurmountable in nature. This behaviour facilitated functional interactions of tiotropium with the beta-agonist olodaterol, which synergistically enhanced bronchoprotective effects of tiotropium. This was not seen for glycopyrrolate and olodaterol or indacaterol but was mimicked by the interaction of tiotropium and forskolin, indicating no direct beta-adrenoceptor-M-3 receptor crosstalk in this effect. Conclusions and Implications We propose that tiotropium has two binding sites at the M-3 receptor that prevent ACh action, which, together with slow dissociation kinetics, may contribute to insurmountable antagonism and enhanced functional interactions with beta-adrenoceptor agonists.

Protein aggregation remains a major area of focus in the production of monoclonal antibodies. Improving the intrinsic properties of antibodies can improve manufacturability, attrition rates, safety, formulation, titers, immunogenicity, and solubility. Here, we explore the potential of predicting and reducing the aggregation propensity of monoclonal antibodies, based on the identification of aggregation-prone regions and their contribution to the thermodynamic stability of the protein. Although aggregation-prone regions are thought to occur in the antigen binding region to drive hydrophobic binding with antigen, we were able to rationally design variants that display a marked decrease in aggregation propensity while retaining antigen binding through the introduction of artificial aggregation gatekeeper residues. The reduction in aggregation propensity was accompanied by an increase in expression titer, showing that reducing protein aggregation is beneficial throughout the development process. The data presented show that this approach can significantly reduce liabilities in novel therapeutic antibodies and proteins, leading to a more efficient path to clinical studies.

Understanding protein function requires detailed knowledge about protein dynamicsProtein dynamics , i.e. the different conformational states the system can adopt. Despite substantial experimental progress, simulation techniques such as molecular dynamicsMolecular dynamics (MD) currently provide the only routine means to obtain dynamical information at an atomic level on timescales of nano- to microseconds. Even with the current development of computational power, sampling techniques beyond MD are necessary to enhance conformational samplingConformational sampling of large proteins and assemblies thereof. The use of collective coordinatesCollective coordinates has proven to be a promising means in this respect, either as a tool for analysis or as part of new sampling algorithms. Starting from MD simulations, several enhanced samplingEnhanced sampling algorithms for biomolecular simulations are reviewed in this chapter. Examples are given throughout illustrating how consideration of the dynamic properties of a protein sheds light on its function.

The bronchodilator drug tiotropium has been suggested to bind an allosteric site on the M3 muscarinic receptor, but its functional relevance is unclear. We hypothesized that binding to this allosteric site would have functional implications for its bronchoprotective effects and its interactions with LABAs. Methods: We used molecular dynamics simulation and contraction studies in bovine and human smooth muscle strips to study the interactions between acetylcholine and tiotropium. Results: Using molecular dynamics, we show that tiotropium binds to the allosteric site and thereby prevents acetylcholine entry into the orthosteric binding pocket. Thus, tiotropium has two receptor binding sites, both of which antagonize acetylcholine action. This resulted in a low (allosteric binding) and high (orthosteric binding) functional affinity of tiotropium in studies of methacholine induced contractions of airway smooth muscle, that were characterized by rapid versus slow kinetics of interaction between methacholine and tiotropium. Antagonism by tiotropium was insurmountable in nature, with methacholine failing to reach maximal response. This behaviour facilitated functional interactions of tiotropium with the β-agonist olodaterol, which substantially potentiated the bronchoprotective effect of tiotropium to a much greater extent than additive. This biphasic behaviour was not observed for glycopyrrolate or atropine. Conclusions: Tiotropium has two binding sites at the M3 receptor that both prevent acetylcholine action, which together with the slow dissociation kinetics contributes to insurmountable antagonism and enhanced functional interactions with β-agonists.

Monoclonal antibody (mAb)-based therapeutics often require high-concentration formulations. Unfortunately, highly concentrated antibody solutions often have biophysical properties that are disadvantageous for therapeutic development, such as high viscosity, solubility limitations, precipitation issues, or liquid liquid phase separation. In this work, we present a computational rational design principle for improving the thermodynamic stability of mAb solutions through targeted point mutations. Two publicly available IgG1 monoclonal antibodies that exhibit high viscosity at high concentrations were used as model systems. Guided by a computationally efficient approach that combines molecular dynamics simulations with three-dimensional reference interaction site model theory, point mutations of charged residues were introduced in the variable Fv regions in such a manner that the hydration free energy was optimized. Two selected point mutants were then produced by transient expression and characterized experimentally. Both engineered mAbs have reduced viscosity at high concentration, less negative second virial coefficient, and improved solubility compared to the respective wild-types. The results obtained with the suggested straightforward design principle underline the relevance of solvation effects for understanding, and ultimately optimizing, the properties of highly concentrated mAb solutions, with possible implications also for other biomolecular systems.

AbstractThe prediction of mutation‐induced free‐energy changes in protein thermostability or protein–protein binding is of particular interest in the fields of protein design, biotechnology, and bioengineering. Herein, we achieve remarkable accuracy in a scan of 762 mutations estimating changes in protein thermostability based on the first principles of statistical mechanics. The remaining error in the free‐energy estimates appears to be due to three sources in approximately equal parts, namely sampling, force‐field inaccuracies, and experimental uncertainty. We propose a consensus force‐field approach, which, together with an increased sampling time, leads to a free‐energy prediction accuracy that matches those reached in experiments. This versatile approach enables accurate free‐energy estimates for diverse proteins, including the prediction of changes in the melting temperature of the membrane protein neurotensin receptor 1.

G protein-coupled receptors (GPCRs) are the largest superfamily of membrane proteins in the human genome, mediating the propagation of extracellular ligand binding information into intracellular signal transduction cascades. Crystal structures have revealed a water-filled hydrophilic internal pocket within their transmembrane domain, extending from the orthosteric ligand-binding site to regions near the G protein binding site. Recent high-resolution structures have identified a sodium ion near the base of this pocket, coordinated by highly conserved residues (1,2). Owing to this conservation level, sodium binding may play a role for most rhodopsin-like GPCRs (3). Recently, functional or conformational effects have been shown to be elicited by physiological transmembrane voltage in several GPCRs, and voltage-clamp recordings have determined a gating charge of ∼0.85 e in the M1 and M2 muscarinic receptors (4). However, the nature of the voltage sensor has remained enigmatic. Here we show by MD and computational electrophysiology simulations (5) that the internal sodium ion in the delta-opioid and muscarinic receptors is highly mobile under the influence of transmembrane voltage changes. By calculating the energetics of ion movement along the pocket axis, we find that the free energy barriers to migration can be overcome by the physiological voltage range. Furthermore, we demonstrate that the motion along the axis creates a gating charge in excellent agreement with the experimentally measured gating currents. Thus, voltage-induced movements of Na+ in GPCRs may represent an attractive mechanism for voltage-regulation of GPCRs. 1. Liu, W. et al., Science 337: 232-236, 2012. 2. Fenalti, G. et al., Nature 506: 191-196, 2014. 3.Katritch, V. et al., TiBS 39: 233-244, 2014. 4. Ben-Chaim, Y. et al., Nature 444: 106-9, 2006. 5. Kutzner, C. et al., Biophysical Journal 101: 809-817, 2011.

Proteins under evolutionary pressure have acquired a subtle balance of thermodynamic properties that optimizes functioning at a physiological range of the environmental conditions. The fields of biotechnology and bioengineering, however, aim to recalibrate, and sometimes design, biological molecules by tailoring certain physico-chemical properties. Enhancing protein thermal stability is one of the frequently desired features both by the industry and scientific community. Increased thermostability may be achieved by introducing mutations into a protein, and identification of such mutations is thus of key interest. In the current study we utilize a recently developed automated approach [1] to perform a large scale protein mutation scan assessing changes in thermostability from free energy molecular dynamics simulations. A remarkable agreement with experimental data was achieved and further enhanced by adopting a consensus force field approach. Subsequently, we demonstrate the strength of the first principles based methods in capturing mutation induced free energy changes in a membrane protein (neurotensin receptor 1). Finally, gaining access to accurate protein thermostabilities enables us to estimate free energy changes in protein-protein binding upon an amino acid mutation. We report on a number of cases examplifying the use of first principles approaches for mutations at a protein-protein binding interface. [1] V. Gapsys, S. Michielssens, D. Seeliger, and B. L. de Groot. pmx: Automated protein structure and topology generation for alchemical perturbations. J. Comput. Chem., 36(5):348-354, 2015.

G-protein-coupled receptors (GPCRs) form the largest superfamily of membrane proteins and one-third of all drug targets in humans. A number of recent studies have reported evidence for substantial voltage regulation of GPCRs. However, the structural basis of GPCR voltage sensing has remained enigmatic. Here, we present atomistic simulations on the δ-opioid and M2 muscarinic receptors, which suggest a structural and mechanistic explanation for the observed voltage-induced functional effects. The simulations reveal that the position of an internal Na(+) ion, recently detected to bind to a highly conserved aqueous pocket in receptor crystal structures, strongly responds to voltage changes. The movements give rise to gating charges in excellent agreement with previous experimental recordings. Furthermore, free energy calculations show that these rearrangements of Na(+) can be induced by physiological membrane voltages. Due to its role in receptor function and signal bias, the repositioning of Na(+) has important general implications for signal transduction in GPCRs.