Molecular Dynamics-Markov State Model of Protein Ligand Binding and Allostery in CRIB-PDZ: Conformational Selection and Induced Fit

Capturing the Effects of Single Atom Substitutions on the Inhibition Efficiency of Glycogen Synthase Kinase-3β Inhibitors via Markov State Modeling and Experiments

Small modifications in the chemical structure of ligands are known to dramatically change their ability to inhibit the activity of a protein. Unraveling the mechanisms that govern these dramatic changes requires scrutinizing the dynamics of protein-ligand binding and unbinding at the atomic level. As an exemplary case, we have studied Glycogen Synthase Kinase-3β (GSK-3β), a multifunctional kinase that has been implicated in a host of pathological processes. As such, there is a keen interest in identifying ligands that inhibit GSK-3β activity. One family of compounds that are highly selective and potent inhibitors of GSK-3β is exemplified by a molecule termed COB-187. COB-187 consists of a five-member heterocyclic ring with a thione at C2, a pyridine substituted methyl at N3, and a hydroxyl and phenyl at C4. We have studied the inhibition of GSK-3β by COB-187-related ligands that differ in a single heavy atom from each other (either in the location of nitrogen in their pyridine ring, or with the pyridine ring replaced by a phenyl ring), or in the length of the alkyl group joining the pyridine and the N3. The inhibition experiments show a large range of half-maximal inhibitory concentration (IC50) values from 10 nM to 10 μM, implying that these ligands exhibit vastly different propensities to inhibit GSK-3β. To explain these differences, we perform Markov State Modeling (MSM) using fully atomistic simulations. Our MSM results are in excellent agreement with the experiments in that they accurately capture differences in the binding propensities of the ligands. The simulations show that the binding propensities are related to the ligands' ability to attain a compact conformation where their two aromatic rings are spatially close. We rationalize this result by sampling numerous binding and unbinding events via funnel metadynamics simulations, which show that indeed while approaching the bound state, the ligands prefer to be in their compact conformation. We find that the presence of nitrogen in the aromatic ring increases the probability of attaining the compact conformation. Protein-ligand binding is understood to be dictated by the energetics of interactions and entropic factors, like the release of bound water from the binding pockets. This work shows that changes in the conformational distribution of ligands due to atom-level modifications in the structure play an important role in protein-ligand binding.

Leveraging Bidirectional Nature of Allostery To Inhibit Protein-Protein Interactions (PPIs): A Case Study of PCSK9-LDLR Interaction

The protein PCSK9 (proprotein convertase subtilisin/Kexin type 9) negatively regulates the recycling of LDLR (low-density lipoprotein receptor), leading to an elevated plasma level of LDL. Inhibition of PCSK9-LDLR interaction has emerged as a promising therapeutic strategy to manage hypercholesterolemia. However, the large interaction surface area between PCSK9 and LDLR makes it challenging to identify a small molecule competitive inhibitor. An alternative strategy would be to identify distal cryptic sites as targets for allosteric inhibitors that can remotely modulate PCSK9-LDLR interaction. Using several microseconds long molecular dynamics (MD) simulations, we demonstrate that on binding with LDLR, there is a significant conformational change (population shift) in a distal loop (residues 211-222) region of PCSK9. Consistent with the bidirectional nature of allostery, we establish a clear correlation between the loop conformation and the binding affinity with LDLR. Using a thermodynamic argument, we establish that the loop conformations predominantly present in the apo state of PCSK9 would have lower LDLR binding affinity, and they would be potential targets for designing allosteric inhibitors. We elucidate the molecular origin of the allosteric coupling between this loop and the LDLR binding interface in terms of the population shift in a set of salt bridges and hydrogen bonds. Overall, our work provides a general strategy toward identifying allosteric hotspots: compare the conformational ensemble of the receptor between the apo and bound states of the protein and identify distal conformational changes, if any. The inhibitors should be designed to bind and stabilize the apo-specific conformations.

Reconnaissance of Allostery via the Restoration of Native p53 DNA-Binding Domain Dynamics in Y220C Mutant p53 Tumor Suppressor Protein

Allosteric regulation of protein dynamics infers a long-range deliberate propagation of information via micro- and macroscale interactions. The Y220C structural mutant is one of the most frequent cancerous p53 mutants. The mutation is distally located from the DNA-binding site of the p53 DNA-binding domain yet causes changes in DNA recognition. This system presents a unique opportunity to examine the allosteric control of mutated proteins under a drug design paradigm. We focus on the key case study of p53 Y220C mutation restoration by a series of new compounds suggested to have Y220C reactivation properties in comparison to our previous findings on the restorative potential of PK11000, a compound studied extensively for reactivation in vitro and in vivo. Previously, we implemented all-atom molecular dynamics (MD) simulations and our lab's techniques of MD-Sectors and MD-Markov state models on the wild type, the Y220C mutant, and Y220C with PK11000 to characterize the effector's restorative properties in terms of conformational dynamics and hydrogen bonding. In this study, we turn to probing the effects made by docking the battery of a new but less well-tested set of aminobenzothiazole derivative compounds reported by Baud et al., which show promise of Y220C rescue. We find that while complete and precise reconstitution of p53 WT molecular dynamics may not be observed as was the case with PK11000, dispersed local reconstitution of loop dynamics provides evidence of rescuing effects by aminobenzothiazole derivative N,2-dihydroxy-3,5-diiodo-4-(1H-pyrrol-1-yl)benzamide, Effector 22, like what we observed for PK11000. Generalizable insights into the mutation and allosteric reactivation of p53 by various effectors by reconstitution of WT dynamics observed in statistical conformational ensemble analysis and network inference are discussed, considering the development of allosteric drug design rooted in first principles.

Revolutionizing Peptide-Based Drug Discovery: Advances in the Post-AlphaFold Era

Peptide-based drugs offer high specificity, potency, and selectivity. However, their inherent flexibility and differences in conformational preferences between their free and bound states create unique challenges that have hindered progress in effective drug discovery pipelines. The emergence of AlphaFold (AF) and Artificial Intelligence (AI) presents new opportunities for enhancing peptide-based drug discovery. We explore recent advancements that facilitate a successful peptide drug discovery pipeline, considering peptides' attractive therapeutic properties and strategies to enhance their stability and bioavailability. AF enables efficient and accurate prediction of peptide-protein structures, addressing a critical requirement in computational drug discovery pipelines. In the post-AF era, we are witnessing rapid progress with the potential to revolutionize peptide-based drug discovery such as the ability to rank peptide binders or classify them as binders/non-binders and the ability to design novel peptide sequences. However, AI-based methods are struggling due to the lack of well-curated datasets, for example to accommodate modified amino acids or unconventional cyclization. Thus, physics-based methods, such as docking or molecular dynamics simulations, continue to hold a complementary role in peptide drug discovery pipelines. Moreover, MD-based tools offer valuable insights into binding mechanisms, as well as the thermodynamic and kinetic properties of complexes. As we navigate this evolving landscape, a synergistic integration of AI and physics-based methods holds the promise of reshaping the landscape of peptide-based drug discovery.

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

Molecular dynamics of mismatch detection-How MutS uses indirect readout to find errors in DNA

The mismatch repair protein MutS safeguards genomic integrity by finding and initiating repair of basepairing errors in DNA. Single-molecule studies show MutS diffusing on DNA, presumably scanning for mispaired/unpaired bases, and crystal structures show a characteristic "mismatch-recognition" complex with DNA enclosed within MutS and kinked at the site of error. But how MutS goes from scanning thousands of Watson-Crick basepairs to recognizing rare mismatches remains unanswered, largely because atomic-resolution data on the search process are lacking. Here, 10 μs all-atom molecular dynamics simulations of Thermus aquaticus MutS bound to homoduplex DNA and T-bulge DNA illuminate the structural dynamics underlying the search mechanism. MutS-DNA interactions constitute a multistep mechanism to check DNA over two helical turns for its 1) shape, through contacts with the sugar-phosphate backbone, 2) conformational flexibility, through bending/unbending engineered by large-scale motions of the clamp domain, and 3) local deformability, through basepair destabilizing contacts. Thus, MutS can localize a potential target by indirect readout due to lower energetic costs of bending mismatched DNA and identify a site that distorts easily due to weaker base stacking and pairing as a mismatch. The MutS signature Phe-X-Glu motif can then lock in the mismatch-recognition complex to initiate repair.

Deconstructing allostery by computational assessment of the binding determinants of allosteric PTP1B modulators

Fragment-based drug discovery is an established methodology for finding hit molecules that can be elaborated into lead compounds. However it is currently challenging to predict whether fragment hits that do not bind to an orthosteric site could be elaborated into allosteric modulators, as in these cases binding does not necessarily translate into a functional effect. We propose a workflow using Markov State Models (MSMs) with steered molecular dynamics (sMD) to assess the allosteric potential of known binders. sMD simulations are employed to sample protein conformational space inaccessible to routine equilibrium MD timescales. Protein conformations sampled by sMD provide starting points for seeded MD simulations, which are combined into MSMs. The methodology is demonstrated on a dataset of protein tyrosine phosphatase 1B ligands. Experimentally confirmed allosteric inhibitors are correctly classified as inhibitors, whereas the deconstructed analogues show reduced inhibitory activity. Analysis of the MSMs provide insights into preferred protein-ligand arrangements that correlate with functional outcomes. The present methodology may find applications for progressing fragments towards lead molecules in FBDD campaigns.

The quaternary question: Determining allostery in spastin through dynamics classification learning and bioinformatics

The nanomachine from the ATPases associated with various cellular activities superfamily, called spastin, severs microtubules during cellular processes. To characterize the functionally important allostery in spastin, we employed methods from evolutionary information, to graph-based networks, to machine learning applied to atomistic molecular dynamics simulations of spastin in its monomeric and the functional hexameric forms, in the presence or absence of ligands. Feature selection, using machine learning approaches, for transitions between spastin states recognizes all the regions that have been proposed as allosteric or functional in the literature. The analysis of the composition of the Markov State Model macrostates in the spastin monomer, and the analysis of the direction of change in the top machine learning features for the transitions, indicate that the monomer favors the binding of ATP, which primes the regions involved in the formation of the inter-protomer interfaces for binding to other protomer(s). Allosteric path analysis of graph networks, built based on the cross-correlations between residues in simulations, shows that perturbations to a hub specific for the pre-hydrolysis hexamer propagate throughout the structure by passing through two obligatory regions: the ATP binding pocket, and pore loop 3, which connects the substrate binding site to the ATP binding site. Our findings support a model where the changes in the terminal protomers due to the binding of ligands play an active role in the force generation in spastin. The secondary structures in spastin, which are found to be highly degenerative within the network paths, are also critical for feature transitions of the classification models, which can guide the design of allosteric effectors to enhance or block allosteric signaling.

Allosteric Signaling in PDZ Energetic Networks: Embedding Error Analysis

Allosteric signaling in proteins has been known for some half a century, yet how the signal traverses the protein remains an active area of research. Recently, the importance of electrostatics to achieve long-range signaling has become increasingly appreciated. Our laboratory has been working on developing network approaches to capture such interactions. In this study, we turn our attention to the well-studied allosteric model protein, PDZ. We study the allosteric dynamics on a per-residue basis in key constructs involving the PDZ domain, its allosteric effector, and its peptide ligand. We utilize molecular dynamics trajectories to create the networks for the constructs to explore the allosteric effect by plotting the heat kernel results onto axes defined by principal components. We introduce a new metric to quantitate the volume sampled by a residue in the latent space. We relate our findings to PDZ and the greater field of allostery.

Network Analysis of Molecular Dynamics Sectors in the p53 Protein

Design of allosteric regulators is an emergent field in the area of drug discovery holding promise for currently untreated diseases. Allosteric regulators bind to a protein in one location and affect a distant site. The ubiquitous presence of allosteric effectors in biology and the success of serendipitously identified allosteric compounds point to the potential they hold. Although the mechanism of transmission of an allosteric signal is not unequivocally determined, one hypothesis suggests that groups of evolutionarily covarying residues within a protein, termed sectors, are conduits. A long-term goal of our lab is to allosterically modulate the activity of proteins by binding small molecules at points of allosteric control. However, methods to consistently identify such points remain unclear. Sector residues on the surfaces of proteins are a promising source of allosteric targets. Recently, we introduced molecular dynamics (MD)-based sectors; MD sectors capitalize on covariance of motion, in place of evolutionary covariance. By focusing on motional covariance, MD sectors tap into the framework of statistical mechanics afforded by the Boltzmann ensemble of structural conformations comprising the underlying data set. We hypothesized that the method of MD sectors can be used to identify a cohesive network of motionally covarying residues capable of transmitting an allosteric signal in a protein. While our initial qualitative results showed promise for the method to predict sectors, that a network of cohesively covarying residues had been produced remained an untested assumption. In this work, we apply network theory to rigorously analyze MD sectors, allowing us to quantitatively assess the biologically relevant property of network cohesiveness of sectors in the context of the tumor suppressor protein, p53. We revised the methodology for assessing and improving MD sectors. Specifically, we introduce a metric to calculate the cohesive properties of the network. Our new approach separates residues into two categories: sector residues and non-sector residues. The relatedness within each respective group is computed with a distance metric. Cohesive sector networks are identified as those that have high relatedness among the sector residues which exceeds the relatedness of the residues to the non-sector residues in terms of the correlation of motions. Our major finding was that the revised means of obtaining sectors was more efficacious than previous iterations, as evidenced by the greater cohesion of the networks. These results are discussed in the context of the development of allosteric regulators of p53 in particular and the expected applicability of the method to the drug design field in general.

Investigating the allosteric response of the PICK1 PDZ domain to different ligands with all-atom simulations

The PDZ family is comprised of small modular domains that play critical roles in the allosteric modulation of many cellular signaling processes by binding to the C-terminal tail of different proteins. As dominant modular proteins that interact with a diverse set of peptides, it is of particular interest to explore how different binding partners induce different allosteric effects on the same PDZ domain. Because the PICK1 PDZ domain can bind different types of ligands, it is an ideal test case to answer this question and explore the network of interactions that give rise to dynamic allostery. Here, we use all-atom molecular dynamics simulations to explore dynamic allostery in the PICK1 PDZ domain by modeling two PICK1 PDZ systems: PICK1 PDZ-DAT and PICK1 PDZ-GluR2. Our results suggest that ligand binding to the PICK1 PDZ domain induces dynamic allostery at the αA helix that is similar to what has been observed in other PDZ domains. We found that the PICK1 PDZ-ligand distance is directly correlated with both dynamic changes of the αA helix and the distance between the αA helix and βB strand. Furthermore, our work identifies a hydrophobic core between DAT/GluR2 and I35 as a key interaction in inducing such dynamic allostery. Finally, the unique interaction patterns between different binding partners and the PICK1 PDZ domain can induce unique dynamic changes to the PICK1 PDZ domain. We suspect that unique allosteric coupling patterns with different ligands may play a critical role in how PICK1 performs its biological functions in various signaling networks.

Variable Regions of p53 Isoforms Allosterically Hard Code DNA Interaction

Allosteric regulation of protein activity pervades biology as the "second secret of life." We have been examining the allosteric regulation and mutant reactivation of the tumor suppressor protein p53. We have found that generalizing the definition of allosteric effector to include entire proteins and expanding the meaning of binding site to include the interface of a transcription factor with its DNA to be useful in understanding the modulation of protein activity. Here, we cast the variable regions of p53 isoforms as allosteric regulators of p53 interactions with its consensus DNA. We implemented molecular dynamics simulations and our lab's new techniques of molecular dynamics (MD) sectors and MD-Markov state models to investigate the effects of nine naturally occurring splice variant isoforms of p53. We find that all of the isoforms differ from wild type in their dynamic properties and how they interact with the DNA. We consider the implications of these findings on allostery and cancer treatment.

Mutagenic Activation of Glutathione Peroxidase-4: Approaches toward Rational Design of Allosteric Drugs

Glutathione peroxidase 4 (GPX4) reduces lipid hydroperoxides in lipid membranes, effectively inhibiting iron-dependent cell death or ferroptosis. The upregulation of the enzyme by the mutations at residues D21 and D23 has been suggested to be associated with higher protein activity, which confers more protection against neurodegenerative diseases such as Alzheimer's, Parkinson's, and Huntington's diseases. Therefore, it has become an attractive target for treating and preventing neurodegenerative diseases. However, identifying means of mimicking the beneficial effects of these mutations distant from the active site constitutes a formidable challenge in moving toward therapeutics. In this study, we explore using molecular dynamics simulations to computationally map the conformational and energetic landscape of the wild-type GPX4 protein and three mutant variants to identify the allosteric networks of the enzyme. We present the conformational dynamic profile providing the desired signature behavior of the enzyme. We also discuss the implications of these findings for drug design efforts.

A litmus test for classifying recognition mechanisms of transiently binding proteins

Partner recognition in protein binding is critical for all biological functions, and yet, delineating its mechanism is challenging, especially when recognition happens within microseconds. We present a theoretical and experimental framework based on straight-forward nuclear magnetic resonance relaxation dispersion measurements to investigate protein binding mechanisms on sub-millisecond timescales, which are beyond the reach of standard rapid-mixing experiments. This framework predicts that conformational selection prevails on ubiquitin’s paradigmatic interaction with an SH3 (Src-homology 3) domain. By contrast, the SH3 domain recognizes ubiquitin in a two-state binding process. Subsequent molecular dynamics simulations and Markov state modeling reveal that the ubiquitin conformation selected for binding exhibits a characteristically extended C-terminus. Our framework is robust and expandable for implementation in other binding scenarios with the potential to show that conformational selection might be the design principle of the hubs in protein interaction networks.

The authors provide a litmus test for the recognition mechanism of transiently binding proteins based on nuclear magnetic resonance and find a conformational selection binding mechanism through concentration-dependent kinetics of ubiquitin and SH3.

The CAR-mRNA Interaction Surface Is a Zipper Extension of the Ribosome A Site

The ribosome CAR interaction surface behaves as an extension of the decoding center A site and has H-bond interactions with the +1 codon, which is next in line to enter the A site. Through molecular dynamic simulations, we investigated the codon sequence specificity of this CAR–mRNA interaction and discovered a strong preference for GCN codons, suggesting that there may be a sequence-dependent layer of translational regulation dependent on the CAR interaction surface. Dissection of the CAR–mRNA interaction through nucleotide substitution experiments showed that the first nucleotide of the +1 codon dominates over the second nucleotide position, consistent with an energetically favorable zipper-like activity that emanates from the A site through the CAR–mRNA interface. Moreover, the CAR/+1 codon interaction is affected by the identity of nucleotide 3 of +1 GCN codons, which influences the stacking of G and C. Clustering analysis suggests that the A-site decoding center adopts different neighborhood substates that depend on the identity of the +1 codon.

Insights into Rational Design of a New Class of Allosteric Effectors with Molecular Dynamics Markov State Models and Network Theory

The development of drugs to restore protein function has been a major advance facilitated by molecular medicine. Allosteric regulation, a phenomenon widely observed in nature, in which a molecule binds to control a distance active site, holds great promise for regulating proteins, yet how to rationally design such a molecule remains a mystery. Over the past few years, we and others have developed several techniques based on molecular dynamics (MD) simulations: MD-Markov state models to capture global conformational substates, and network theory approach utilizing the interaction energy within the protein to confer local allosteric control. We focus on the key case study of the p53 Y220C mutation restoration by PK11000, a compound experimentally shown to reactivate p53 native function in Y220C mutant present tumors. We gain insights into the mutation and allosteric reactivation of the protein, which we anticipate will be applicable to de novo design to engineer new compounds not only for this mutation, but in other macromolecular systems as well.

Estradiol induces allosteric coupling and partitioning of sex-hormone-binding globulin monomers among conformational states

Sex-hormone-binding globulin (SHBG) regulates the transport and bioavailability of estradiol. The dynamics of estradiol's binding to SHBG are incompletely understood, although it is believed that estradiol binds to each monomer of SHBG dimer with identical affinity (Kd ∼2 nM). Contrary to the prevalent view, we show that estradiol's binding to SHBG is nonlinear, and the "apparent" Kd changes with varying estradiol and SHBG concentrations. Estradiol's binding to each SHBG monomer influences residues in the ligand-binding pocket of both monomers and differentially alters the conformational and energy landscapes of both monomers. Monomers are not energetically or conformationally equivalent even in fully bound state. Estradiol's binding to SHBG involves bidirectional, inter-monomeric allostery that changes the distribution of both monomers among various energy and conformational states. Inter-monomeric allostery offers a mechanism to extend the binding range of SHBG and regulate hormone bioavailability as estradiol concentrations vary widely during life.

The Speed of Allosteric Signaling Within a Single-Domain Protein

While much is known about different allosteric regulation mechanisms, the nature of the allosteric signal and the time scale on which it propagates remains elusive. The PDZ3 domain from postsynaptic density-95 protein is a small protein domain with a terminal third α-helix, i.e., the α3-helix, which is known to be allosterically active. By cross-linking the allosteric helix with an azobenzene moiety, we obtained a photocontrollable PDZ3 variant. Photoswitching triggers its allosteric transition, resulting in a change in binding affinity of a peptide to the remote binding pocket. Using time-resolved infrared and UV/vis spectroscopy, we follow the allosteric signal transduction and reconstruct the timeline in which the allosteric signal propagates through the protein within 200 ns.

Identifying Conformational Isomers of Organic Molecules in Solution via Unsupervised Clustering

We present a systematic approach for the identification of statistically relevant conformational macrostates of organic molecules from molecular dynamics trajectories. The approach applies to molecules characterized by an arbitrary number of torsional degrees of freedom and enables the transferability of the macrostates definition across different environments. We formulate a dissimilarity measure between molecular configurations that incorporates information on the characteristic energetic cost associated with transitions along all relevant torsional degrees of freedom. Such metric is employed to perform unsupervised clustering of molecular configurations based on the Fast Search and Find of Density Peaks algorithm. We apply this method to investigate the equilibrium conformational ensemble of Sildenafil, a conformationally complex pharmaceutical compound, in different environments including the crystal bulk, the gas phase, and three different solvents (acetonitrile, 1-butanol, and toluene). We demonstrate that while Sildenafil can adopt more than 100 metastable conformational configurations, only 12 are significantly populated across all of the environments investigated. Despite the complexity of the conformational space, we find that the most abundant conformers in solution are the closest to the conformers found in the most common Sildenafil crystal phase.

Molecular dynamics of the immune checkpoint programmed cell death protein I, PD-1: conformational changes of the BC-loop upon binding of the ligand PD-L1 and the monoclonal antibody nivolumab

BACKGROUND

The immune checkpoint receptor programmed cell death protein I (PD-1) has been identified as a key target in immunotherapy. PD-1 reduces the risk of autoimmunity by inducing apoptosis in antigen-specific T cells upon interaction with programmed cell death protein ligand I (PD-L1). Various cancer types overexpress PD-L1 to evade the immune system by inducing apoptosis in tumor-specific CD8+ T cells. The clinically used blocking antibody nivolumab binds to PD-1 and inhibits the immunosuppressive interaction with PD-L1. Even though PD-1 is already used as a drug target, the exact mechanism of the receptor is still a matter of debate. For instance, it is hypothesized that the signal transduction is based on an active conformation of PD-1.

RESULTS

Here we present the results of the first molecular dynamics simulations of PD-1 with a complete extracellular domain with a focus on the role of the BC-loop of PD-1 upon binding PD-L1 or nivolumab. We could demonstrate that the BC-loop can form three conformations. Nivolumab binds to the BC-loop according to the conformational selection model whereas PD-L1 induces allosterically a conformational change of the BC-loop.

CONCLUSION

Due to the structural differences of the BC-loop, a signal transduction based on active conformation cannot be ruled out. These findings will have an impact on drug design and will help to refine immunotherapy blocking antibodies.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

Ranking of Ligand Binding Kinetics Using a Weighted Ensemble Approach and Comparison with a Multiscale Milestoning Approach

To improve lead optimization efforts in finding the right ligand, pharmaceutical industries need to know the ligand's binding kinetics, such as binding and unbinding rate constants, which often correlate with the ligand's efficacy in vivo. To predict binding kinetics efficiently, enhanced sampling methods, such as milestoning and the weighted ensemble (WE) method, have been used in molecular dynamics (MD) simulations of these systems. However, a comparison of these enhanced sampling methods in ranking ligands has not been done. Hence, a WE approach called the concurrent adaptive sampling (CAS) algorithm that uses MD simulations was used to rank seven ligands for β-cyclodextrin, a system in which a multiscale milestoning approach called simulation enabled estimation of kinetic rates (SEEKR) was also used, which uses both MD and Brownian dynamics simulations. Overall, the CAS algorithm can successfully rank ligands using the unbinding rate constant k_off values and binding free energy ΔG values, as SEEKR did, with reduced computational cost that is about the same as SEEKR. We compare the CAS algorithm simulations with different parameters and discuss the impact of parameters in ranking ligands and obtaining rate constant and binding free energy estimates. We also discuss similarities and differences and advantages and disadvantages of SEEKR and the CAS algorithm for future use.

A disordered encounter complex is central to the yeast Abp1p SH3 domain binding pathway

Protein-protein interactions are involved in a wide range of cellular processes. These interactions often involve intrinsically disordered proteins (IDPs) and protein binding domains. However, the details of IDP binding pathways are hard to characterize using experimental approaches, which can rarely capture intermediate states present at low populations. SH3 domains are common protein interaction domains that typically bind proline-rich disordered segments and are involved in cell signaling, regulation, and assembly. We hypothesized, given the flexibility of SH3 binding peptides, that their binding pathways include multiple steps important for function. Molecular dynamics simulations were used to characterize the steps of binding between the yeast Abp1p SH3 domain (AbpSH3) and a proline-rich IDP, ArkA. Before binding, the N-terminal segment 1 of ArkA is pre-structured and adopts a polyproline II helix, while segment 2 of ArkA (C-terminal) adopts a 310 helix, but is far less structured than segment 1. As segment 2 interacts with AbpSH3, it becomes more structured, but retains flexibility even in the fully engaged state. Binding simulations reveal that ArkA enters a flexible encounter complex before forming the fully engaged bound complex. In the encounter complex, transient nonspecific hydrophobic and long-range electrostatic contacts form between ArkA and the binding surface of SH3. The encounter complex ensemble includes conformations with segment 1 in both the forward and reverse orientation, suggesting that segment 2 may play a role in stabilizing the correct binding orientation. While the encounter complex forms quickly, the slow step of binding is the transition from the disordered encounter ensemble to the fully engaged state. In this transition, ArkA makes specific contacts with AbpSH3 and buries more hydrophobic surface. Simulating the binding between ApbSH3 and ArkA provides insight into the role of encounter complex intermediates and nonnative hydrophobic interactions for other SH3 domains and IDPs in general.

Correlation functions, mean first passage times, and the Kemeny constant

Markov processes are widely used models for investigating kinetic networks. Here, we collate and present a variety of results pertaining to kinetic network models in a unified framework. The aim is to lay out explicit links between several important quantities commonly studied in the field, including mean first passage times (MFPTs), correlation functions, and the Kemeny constant. We provide new insights into (i) a simple physical interpretation of the Kemeny constant, (ii) a relationship to infer equilibrium distributions and rate matrices from measurements of MFPTs, and (iii) a protocol to reduce the dimensionality of kinetic networks based on specific requirements that the MFPTs in the coarse-grained system should satisfy. We prove that this protocol coincides with the one proposed by Hummer and Szabo [J. Phys. Chem. B 119, 9029 (2014)], and it leads to a variational principle for the Kemeny constant. Finally, we introduce a modification of this protocol, which preserves the Kemeny constant. Our work underpinning the theoretical aspects of kinetic networks will be useful in applications including milestoning and path sampling algorithms in molecular simulations.

Spectral analysis of molecular dynamics simulations on PDZ: MD sectors

The idea of protein "sectors" posits that sparse subsets of amino acid residues form cooperative networks that are key elements of protein stability, ligand binding, and allosterism. To date, protein sectors have been calculated by the statistical coupling analysis (SCA) method of Ranganathan and co-workers via the spectral analysis of conservation-weighted evolutionary covariance matrices obtained from a multiple sequence alignments of homologous families of proteins. SCA sectors, a knowledge-based protocol, have been indentified with functional properties and allosterism for a number of systems. In this study, we investigate the utility of the sector idea for the analysis of physics-based molecular dynamics (MD) trajectories of proteins. Our test case for this procedure is PSD95- PDZ3, one of the smallest proteins for which allosterism has been observed. It has served previously as a model system for a number of prediction algorithms, and is well characterized by X-ray crystallography, NMR spectroscopy and site specific mutagenisis. All-atom MD simulations were performed for a total of 500 nanoseconds using AMBER, and MD-calculated covariance matrices for the fluctuations of residue displacements and non-bonded interaction energies were subjected to spectral analysis in a manner analogous to that of SCA. The composition of MD sectors was compared with results from SCA, site specific mutagenesis, and allosterism. The concordance indicates that MD sectors are a viable protocol for analyzing MD trajectories and provide insight into the physical origin of the phenomenon.Communicated by Ramaswamy H. Sarma.

Molecular Mechanism for Attractant Signaling to DHMA by E. coli Tsr

The attractant chemotaxis response of Escherichia coli to norepinephrine requires that it be converted to 3,4-dihydroxymandelic acid (DHMA) by the monoamine oxidase TynA and the aromatic aldehyde dehydrogenase FeaB. DHMA is sensed by the serine chemoreceptor Tsr, and the attractant response requires that at least one subunit of the periplasmic domain of the Tsr homodimer (pTsr) has an intact serine-binding site. DHMA that is generated in vivo by E. coli is expected to be a racemic mixture of the (R) and (S) enantiomers, so it has been unclear whether one or both chiral forms are active. Here, we used a combination of state-of-the-art tools in molecular docking and simulations, including an in-house simulation-based docking protocol, to investigate the binding properties of (R)-DHMA and (S)-DHMA to E. coli pTsr. Our studies computationally predicted that (R)-DHMA should promote a stronger attractant response than (S)-DHMA because of a consistently greater-magnitude piston-like pushdown of the pTsr α-helix 4 toward the membrane upon binding of (R)-DHMA than upon binding of (S)-DHMA. This displacement is caused primarily by interaction of DHMA with Tsr residue Thr156, which has been shown by genetic studies to be critical for the attractant response to L-serine and DHMA. These findings led us to separate the two chiral species and test their effectiveness as chemoattractants. Both the tethered cell and motility migration coefficient assays validated the prediction that (R)-DHMA is a stronger attractant than (S)-DHMA. Our study demonstrates that refined computational docking and simulation studies combined with experiments can be used to investigate situations in which subtle differences between ligands may lead to diverse chemotactic responses.

Identification of a distal allosteric ligand binding pocket in HtrA3

Human HtrA3 (High temperature requirement protease A3) is a trimeric PDZ bearing propapoptotic serine protease, which is involved in various diseases including cancer and pre-eclampsia. Proposed to be a tumor suppressor, its role as a potential therapeutic target is strongly advocated. Therefore, it becomes imperative to gain insights into its mechanism of action and regulation. Allostery is a well-known mechanism of catalytic activation for many HtrA3 homologs, which opens up avenues for manipulating enzyme functions for therapeutic intervention. In our study, through in silico and biochemical approaches, we have reported for the first time that HtrA3 shows allosteric behaviour. We identified a novel selective binding pocket, which triggers conformational reorientations through signal propagation to the distantly situated active-site pocket via the functionally important loop regions. Using molecular docking, simulation studies and biochemical studies we have identified the regulatory movements at and around the active site pocket. Our study is the first one to report a non-classical binding site for HtrA3, which is instrumental for formation of a catalytically efficient orthosteric pocket upon substrate binding.

Functional Relevance of Interleukin-1 Receptor Inter-domain Flexibility for Cytokine Binding and Signaling

The interleukin 1 (IL-1) receptor family, whose members contain three immunoglobulin-like domains (D1-D3) in the extracellular region, is responsible for transmitting pleiotropic signals of IL-1 cytokines. The inter-domain flexibility of IL-1 receptors and its functional roles have not been fully elucidated. In this study, we used small-angle X-ray scattering to show that ligand-binding primary receptors and co-receptors in the family all have inherent inter-domain flexibility due to the D2/D3 linker. Variants of the IL-1RAcP and IL-18Rβ co-receptors with mutated D2/D3 linkers cannot form a cytokine-receptor complex and mediate signaling. Our analysis further revealed that these mutated co-receptors exhibited a changed conformational ensemble, suggesting that loss of function is due to the alteration of receptor dynamics. Taken together, our results demonstrate that the D2/D3 linker is a critical functional determinant of IL-1 receptor and underscore the important roles of the inter-domain flexibility in cytokine/receptor binding and signaling.

Markov models for the elucidation of allosteric regulation

Allosteric regulation refers to the process where the effect of binding of a ligand at one site of a protein is transmitted to another, often distant, functional site. In recent years, it has been demonstrated that allosteric mechanisms can be understood by the conformational ensembles of a protein. Molecular dynamics (MD) simulations are often used for the study of protein allostery as they provide an atomistic view of the dynamics of a protein. However, given the wealth of detailed information hidden in MD data, one has to apply a method that allows extraction of the conformational ensembles underlying allosteric regulation from these data. Markov state models are one of the most promising methods for this purpose. We provide a short introduction to the theory of Markov state models and review their application to various examples of protein allostery studied by MD simulations. We also include a discussion of studies where Markov modelling has been employed to analyse experimental data on allosteric regulation. We conclude our review by advertising the wider application of Markov state models to elucidate allosteric mechanisms, especially since in recent years it has become straightforward to construct such models thanks to software programs like PyEMMA and MSMBuilder.This article is part of a discussion meeting issue 'Allostery and molecular machines'.

Single-Molecule FRET Kinetics of the Mn2+ Riboswitch: Evidence for Allosteric Mg2+ Control of "Induced-Fit" vs "Conformational Selection" Folding Pathways

Gene expression in bacteria is often regulated dynamically by conformational changes in a riboswitch upon ligand binding, a detailed understanding of which is very much in its infancy. For example, the manganese riboswitch is a widespread RNA motif that conformationally responds in regulating bacterial gene expression to micromolar levels of its eponymous ligand, Mn²⁺, but the mechanistic pathways are poorly understood. In this work, we quantitatively explore the dynamic folding behavior of the manganese riboswitch by single-molecule fluorescence resonance energy transfer spectroscopy as a function of cation/ligand conditions. From the detailed analysis of the kinetics, the Mn²⁺ is shown to fold the riboswitch by a "bind-then-fold" (i.e., "induced-fit", IF) mechanism, whereby the ligand binds first and then promotes folding. On the other hand, the data also clearly reveal the presence of a folded yet ligand-free structure predominating due to the addition of physiological Mg²⁺ to a nonselective metal ion binding site. Of particular kinetic interest, such a Mg²⁺ "prefolded" conformation of the riboswitch is shown to exhibit a significantly increased affinity for Mn²⁺ and further stabilization by subsequent binding of the ligand, thereby promoting efficient riboswitch folding by a "fold-then-bind" (i.e., "conformational selection", CS) mechanism. Our results not only demonstrate Mg²⁺-controlled switching between IF and CS riboswitch folding pathways but also suggest a novel heterotropic allosteric control in the manganese riboswitch activity co-regulated by Mg²⁺ binding.

Rationalization of the Membrane Permeability Differences in a Series of Analogue Cyclic Decapeptides

Cyclization and selected backbone N-methylations are found to be often necessary but not sufficient conditions for peptidic drugs to have a good bioavailability. Thus, the design of cyclic peptides with good passive membrane permeability and good solubility remains a challenge. The backbone scaffold of a recently published series of cyclic decapeptides with six selected backbone N-methylations was designed to favor the adoption of a closed conformation with β-turns and four transannular hydrogen bonds. Although this conformation was indeed adopted by the peptides as determined by NMR measurements, substantial differences in the membrane permeability were observed. In this work, we aim to rationalize the impact of discrete side chain modifications on membrane permeability for six of these cyclic decapeptides. The thermodynamic and kinetic properties were investigated using molecular dynamics simulations and Markov state modeling in water and chloroform. The study highlights the influence that side-chain modifications can have on the backbone conformation. Peptides with a d-proline in the β-turns were more likely to adopt, even in water, the closed conformation with transannular hydrogen bonds, which facilitates transition through the membrane. The population of the closed conformation in water was found to correlate positively with PAMPA log P_e.

Exploring Configuration Space and Path Space of Biomolecules Using Enhanced Sampling Techniques-Searching for Mechanism and Kinetics of Biomolecular Functions

To understand functions of biomolecules such as proteins, not only structures but their conformational change and kinetics need to be characterized, but its atomistic details are hard to obtain both experimentally and computationally. Here, we review our recent computational studies using novel enhanced sampling techniques for conformational sampling of biomolecules and calculations of their kinetics. For efficiently characterizing the free energy landscape of a biomolecule, we introduce the multiscale enhanced sampling method, which uses a combined system of atomistic and coarse-grained models. Based on the idea of Hamiltonian replica exchange, we can recover the statistical properties of the atomistic model without any biases. We next introduce the string method as a path search method to calculate the minimum free energy pathways along a multidimensional curve in high dimensional space. Finally we introduce novel methods to calculate kinetics of biomolecules based on the ideas of path sampling: one is the Onsager⁻Machlup action method, and the other is the weighted ensemble method. Some applications of the above methods to biomolecular systems are also discussed and illustrated.

Markov-state model for CO2 binding with carbonic anhydrase under confinement

Enzyme immobilization with a nanostructure material can enhance its stability and facilitate reusability. However, the apparent activity is often compromised due to additional diffusion barriers and complex interactions with the substrates and solvent molecules. The present study elucidates the effects of the surface hydrophobicity of nano-confinement on CO₂ diffusion to the active site of human carbonic anhydrase II (CA), an enzyme that is able to catalyze CO₂ hydration at extremely high turnover rates. Using the Markov-state model in combination with coarse-grained molecular dynamics simulations, we demonstrate that a hydrophobic cage increases CO₂ local density but hinders its diffusion towards the active site of CA under confinement. By contrast, a hydrophilic cage hinders CO₂ adsorption but promotes its binding with CA. An optimal surface hydrophobicity can be identified to maximize both the CO₂ occupation probability and the diffusion rate. The simulation results offer insight into understanding enzyme performance under nano-confinement and help us to advance broader applications of CA for CO₂ absorption and recovery.