Markov state models of biomolecular conformational dynamics

Potential Energy Weighted Reactive Flux and Total Rate of Change of Potential Energy: Theory and Illustrative Applications

Reactive flux can be largely nonzero in a nonequilibrium ensemble of trajectories and provide insightful information for reactive transitions from the reactant state to the product state. Based on the reactive flux, a theoretical framework is proposed here for two quantities, the potential energy weighted reactive flux and the total rate of change of potential energy, which are useful for the identification of the mechanism from a nonequilibrium ensemble. From such quantities, two multidimensional free-energy analogues can be derived in the subspace of collective variables and they are equivalent in the regions where the reactive flux is divergence-free. These free-energy analogues are assumed to be closely related to the free energy in the subspace of collective variables, and they are reduced in the one-dimensional case to be the ensemble average of the potential energy weighted with reactive flux intensity, which was proposed recently [Li, W. J. Phys. Chem. A 2022, DOI: 10.1021/acs.jpca.2c04130] and could be decomposed into energy components at the per-coordinate level. In the subspace of collective variables, the decomposition of the multidimensional free-energy analogues at the per-coordinate level is theoretically possible and is numerically difficult to be calculated. Interestingly, the total rate of change of potential energy is able to identify the location of the transition state ensemble or the stochastic separatrix, in addition to the locations of the reactant and product states. The total rate of change of potential energy can be decomposed at the per-coordinate level, and its components can quantify the contribution of a coordinate to the reactive transition in the subspace of collective variables. We then illustrated the main insights and objects that can be provided by the approach in the applications to a two-dimensional system with various diffusion anisotropies and the alanine peptide in vacuum in various nonequilibrium ensembles of short trajectories, from which the results were found to be consistent.

Markov state modelling reveals heterogeneous drug-inhibition mechanism of Calmodulin

Calmodulin (CaM) is a calcium sensor which binds and regulates a wide range of target-proteins. This implicitly enables the concentration of calcium to influence many downstream physiological responses, including muscle contraction, learning and depression. The antipsychotic drug trifluoperazine (TFP) is a known CaM inhibitor. By binding to various sites, TFP prevents CaM from associating to target-proteins. However, the molecular and state-dependent mechanisms behind CaM inhibition by drugs such as TFP are largely unknown. Here, we build a Markov state model (MSM) from adaptively sampled molecular dynamics simulations and reveal the structural and dynamical features behind the inhibitory mechanism of TFP-binding to the C-terminal domain of CaM. We specifically identify three major TFP binding-modes from the MSM macrostates, and distinguish their effect on CaM conformation by using a systematic analysis protocol based on biophysical descriptors and tools from machine learning. The results show that depending on the binding orientation, TFP effectively stabilizes features of the calcium-unbound CaM, either affecting the CaM hydrophobic binding pocket, the calcium binding sites or the secondary structure content in the bound domain. The conclusions drawn from this work may in the future serve to formulate a complete model of pharmacological modulation of CaM, which furthers our understanding of how these drugs affect signaling pathways as well as associated diseases.

Calmodulin (CaM) is a calcium-sensing protein which makes other proteins dependent on the surrounding calcium concentration by binding to these proteins. Such protein-protein interactions with CaM are vital for calcium to control many physiological pathways within the cell. The antipsychotic drug trifluoperazine (TFP) inhibits CaM’s ability to bind and regulate other proteins. Here, we use molecular dynamics simulations together with Markov state modeling and machine learning to understand the structural and dynamical features by which TFP bound to the one domain of CaM prevents association to other proteins. We find that TFP encourages CaM to adopt a conformation that is like the one stabilized in absence of calcium: depending on the binding orientation of TFP, the drug indeed either affects the CaM hydrophobic binding pocket, the calcium binding sites or the secondary structure content in the domain. Understanding TFP binding is a first step towards designing better drugs targeting CaM.

Deciphering the QR Code of the CRISPR-Cas9 System: Synergy between Gln768 (Q) and Arg976 (R)

Markov state models (MSMs) and machine learning (ML) algorithms can extrapolate the long-time-scale behavior of large biomolecules from molecular dynamics (MD) trajectories. In this study, an MD-MSM-ML scheme has been applied to probe the large endonuclease (Cas9) in the bacterial adaptive immunity CRISPR-Cas9 system. CRISPR has become a programmable and state-of-the-art powerful genome editing tool that has already revolutionized life sciences. CRISPR-Cas9 is programmed to process specific DNA sequences in the genome. However, human/biomedical applications are compromised by off-target DNA damage. Characterization of Cas9 at the structural and biophysical levels is a prerequisite for the development of efficient and high-fidelity Cas9 variants. The Cas9 wild type and two variants (R63A-R66A-R70A, R69A-R71A-R74A-R78A) are studied herein. The configurational space of Cas9 is provided with a focus on the conformations of the side chains of two residues (Gln768 and Arg976). A model for the synergy between those two residues is proposed. The results are discussed within the context of experimental literature. The results and methodology can be exploited for the study of large biomolecules in general and for the engineering of more efficient and safer Cas9 variants for applications.

Structure and Dynamics of the Isozymes II and IX of Human Carbonic Anhydrase

Human carbonic anhydrases (HCAs) are responsible for the pH control and sensing in our body and constitute key components in the central pH paradigm connected to cancer therapeutics. However, little or no molecular level studies are available on the pH-dependent stability and functional dynamics of the known isozymes of HCA. The main objective of this Article is to report the first bench-marking study on the structure and dynamics of the two most efficient isozymes, HCA II and IX, at neutral pH using classical molecular dynamics (MD) and constant pH MD (CpHMD) simulations combined with umbrella sampling, transition path sampling, and Markov state models. Starting from the known crystal structures of HCA II and the monomeric catalytic domain of HCA IX (labeled as HCA IX-c), we have generated classical MD and CpHMD trajectories (of length 1 μs each). In all cases, the overall stability, RMSD, and secondary structure segments of the two isozymes are found to be quite similar. Functionally important dynamics of these two enzymes have been probed in terms of active site hydration, coordination of the Zn(II) ion to a transient excess water, and the formation of putative proton transfer paths. The most important difference between the two isozymes is observed for the side-chain fluctuations of His-64 that is expected to shuttle an excess proton out of the active site as a part of the rate-determining intramolecular proton transfer reaction. The relative stability of the stable inward and outward conformations of the His-64 side-chain and the underlying free energy surfaces are found to depend strongly on the isozyme. In each case, a lower free energy barrier is detected between predominantly inward conformations from predominantly outward ones when simulated under constant pH conditions. The kinetic rate constants of interconversion between different free energy basins are found to span 10⁷-10⁸ s^-1 with faster conformational transitions predicted at constant pH condition. The estimated rate constants and free energies are expected to validate if the fluctuation of the His-64 side-chain in HCA IX may have a significance similar to that known in the multistep catalytic cycle of HCA II.

Molecular mechanisms of spontaneous curvature and softening in complex lipid bilayer mixtures

Membrane reshaping is an essential biological process. The chemical composition of lipid membranes determines their mechanical properties and thus the energetics of their shape. Hundreds of distinct lipid species make up native bilayers, and this diversity complicates efforts to uncover what compositional factors drive membrane stability in cells. Simplifying assumptions, therefore, are used to generate quantitative predictions of bilayer dynamics based on lipid composition. One assumption commonly used is that "per lipid" mechanical properties are both additive and constant-that they are an intrinsic property of lipids independent of the surrounding composition. Related to this is the assumption that lipid bulkiness, or "shape," determines its curvature preference, independently of context. In this study, all-atom molecular dynamics simulations on three separate multilipid systems were used to explicitly test these assumptions, applying methodology recently developed to isolate properties of single lipids or nanometer-scale patches of lipids. The curvature preference experienced by populations of lipid conformations were inferred from their redistribution on a dynamically fluctuating bilayer. Representative populations were extracted by both structural similarity and semi-automated hidden Markov model analysis. The curvature preferences of lipid dimers were then determined and compared with an additive model that combines the monomer curvature preference of both the individual lipids. In all three systems, we identified conformational subpopulations of lipid dimers that showed non-additive curvature preference, in each case mediated by a special chemical interaction (e.g., hydrogen bonding). Our study highlights the importance of specific chemical interactions between lipids in multicomponent bilayers and the impact of interactions on bilayer stiffness. We identify two mechanisms of bilayer softening: diffusional softening, driven by the dynamic coupling between lipid distributions and membrane undulations, and conformational softening, driven by the inter-conversion between distinct dimeric conformations.

The influence of antibody humanization on shark variable domain (VNAR) binding site ensembles

Sharks and other cartilaginous fish produce new antigen receptor (IgNAR) antibodies, as key part of their humoral immune response and are the phylogenetically oldest living organisms that possess an immunoglobulin (Ig)-based adaptive immune system. IgNAR antibodies are naturally occurring heavy-chain-only antibodies, that recognize antigens with their single domain variable regions (VNARs). In this study, we structurally and biophysically elucidate the effect of antibody humanization of a previously published spiny dogfish VNAR (parent E06), which binds with high affinity to the human serum albumin (HSA). We analyze different humanization variants together with the parental E06 VNAR and the human Vκ1 light chain germline DPK9 antibody to characterize the influence of point mutations in the framework and the antigen binding site on the specificity of VNARs as reported by Kovalenko et al. We find substantially higher flexibility in the humanized variants, reflected in a broader conformational space and a higher conformational entropy, as well as population shifts of the dominant binding site ensembles in solution. A further variant, in which some mutations are reverted, largely restores the conformational stability and the dominant binding minimum of the parent E06. We also identify differences in surface hydrophobicity between the human Vκ1 light chain germline DPK9 antibody, the parent VNAR E06 and the humanized variants. Additional simulations of VNAR-HSA complexes of the parent E06 VNAR and a humanized variant reveal that the parent VNAR features a substantially stronger network of stabilizing interactions. Thus, we conclude that a structural and dynamic understanding of the VNAR binding site upon humanization is a key aspect in antibody humanization.

Transition Rates and Efficiency of Collective Variables from Time-Dependent Biased Simulations

Simulations with adaptive time-dependent bias enable an efficient exploration of the conformational space of a system. However, the dynamic information is altered by the bias. Infrequent metadynamics recovers the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we use Kramers' theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the efficiency of collective variables parameter by measuring how efficiently the bias accelerates the transitions. We present approximate analytical expressions of the survival probability, reproducing the barrier-crossing time statistics and enabling the extraction of the unbiased transition rate even for challenging cases. We explore the limits of our method and provide convergence criteria to assess its validity.

Breathing and tilting: mesoscale simulations illuminate influenza glycoprotein vulnerabilities

Influenza virus has resurfaced recently from inactivity during the early stages of the COVID-19 pandemic, raising serious concerns about the nature and magnitude of future epidemics. The main antigenic targets of influenza virus are two surface glycoproteins, hemagglutinin (HA) and neuraminidase (NA). Whereas the structural and dynamical properties of both glycoproteins have been studied previously, the understanding of their plasticity in the whole-virion context is fragmented. Here, we investigate the dynamics of influenza glycoproteins in a crowded protein environment through mesoscale all-atom molecular dynamics simulations of two evolutionary-linked glycosylated influenza A whole-virion models. Our simulations reveal and kinetically characterize three main molecular motions of influenza glycoproteins: NA head tilting, HA ectodomain tilting, and HA head breathing. The flexibility of HA and NA highlights antigenically relevant conformational states, as well as facilitates the characterization of a novel monoclonal antibody, derived from human convalescent plasma, that binds to the underside of the NA head. Our work provides previously unappreciated views on the dynamics of HA and NA, advancing the understanding of their interplay and suggesting possible strategies for the design of future vaccines and antivirals against influenza.

Binding mechanism of oseltamivir and influenza neuraminidase suggests perspectives for the design of new anti-influenza drugs

Oseltamivir is a widely used influenza virus neuraminidase (NA) inhibitor that prevents the release of new virus particles from host cells. However, oseltamivir-resistant strains have emerged, but effective drugs against them have not yet been developed. Elucidating the binding mechanisms between NA and oseltamivir may provide valuable information for the design of new drugs against NA mutants resistant to oseltamivir. Here, we conducted large-scale (353.4 μs) free-binding molecular dynamics simulations, together with a Markov State Model and an importance-sampling algorithm, to reveal the binding process of oseltamivir and NA. Ten metastable states and five major binding pathways were identified that validated and complemented previously discovered binding pathways, including the hypothesis that oseltamivir can be transferred from the secondary sialic acid binding site to the catalytic site. The discovery of multiple new metastable states, especially the stable bound state containing a water-mediated hydrogen bond between Arg118 and oseltamivir, may provide new insights into the improvement of NA inhibitors. We anticipated the findings presented here will facilitate the development of drugs capable of combating NA mutations.

Influenza virus neuraminidase (NA), a viral membrane glycoprotein, plays an important role in the interactions with host cell surface receptors. The emergence and spread of influenza mutants resistant to neuraminidase inhibitors (NAIs), such as oseltamivir, has been of great concern. Despite many improvements to NAIs, no new first-line NAIs are currently in clinical use. Although there have been previous molecular dynamics simulation studies on the binding and dissociation process of oseltamivir-NA, we discovered new binding pathways and states of oseltamivir through larger-scale simulations and more systematic analysis, which may provide new ideas for the improvement of oseltamivir and even a series of NAIs. In our study, we strongly demonstrate that a detailed understanding of the drug−receptor association process is of fundamental importance for drug design and provide methodological references for the improvement of other drugs.

Fewer Dimensions, More Structures for Improved Discrete Models of Dynamics of Free versus Antigen-Bound Antibody

Over the past decade, Markov State Models (MSM) have emerged as powerful methodologies to build discrete models of dynamics over structures obtained from Molecular Dynamics trajectories. The identification of macrostates for the MSM is a central decision that impacts the quality of the MSM but depends on both the selected representation of a structure and the clustering algorithm utilized over the featurized structures. Motivated by a large molecular system in its free and bound state, this paper investigates two directions of research, further reducing the representation dimensionality in a non-parametric, data-driven manner and including more structures in the computation. Rigorous evaluation of the quality of obtained MSMs via various statistical tests in a comparative setting firmly shows that fewer dimensions and more structures result in a better MSM. Many interesting findings emerge from the best MSM, advancing our understanding of the relationship between antibody dynamics and antibody–antigen recognition.

Methodology-Centered Review of Molecular Modeling, Simulation, and Prediction of SARS-CoV-2

Despite tremendous efforts in the past two years, our understanding of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), virus-host interactions, immune response, virulence, transmission, and evolution is still very limited. This limitation calls for further in-depth investigation. Computational studies have become an indispensable component in combating coronavirus disease 2019 (COVID-19) due to their low cost, their efficiency, and the fact that they are free from safety and ethical constraints. Additionally, the mechanism that governs the global evolution and transmission of SARS-CoV-2 cannot be revealed from individual experiments and was discovered by integrating genotyping of massive viral sequences, biophysical modeling of protein-protein interactions, deep mutational data, deep learning, and advanced mathematics. There exists a tsunami of literature on the molecular modeling, simulations, and predictions of SARS-CoV-2 and related developments of drugs, vaccines, antibodies, and diagnostics. To provide readers with a quick update about this literature, we present a comprehensive and systematic methodology-centered review. Aspects such as molecular biophysics, bioinformatics, cheminformatics, machine learning, and mathematics are discussed. This review will be beneficial to researchers who are looking for ways to contribute to SARS-CoV-2 studies and those who are interested in the status of the field.

Time-Lagged Flux in the Transition Path Ensemble: Flux Maximization and Relation to Transition Path Theory

The transition path ensemble is of special interest in reaction coordinate identification as it consists of reactive trajectories that start from the reactant state and end in the product one. As a theoretical framework for describing the transition path ensemble, the transition path theory has been introduced more than 10 years ago, and so far, its applications have only been illustrated in several low-dimensional systems. Given the transition path ensemble, expressions for calculating flux, current (a vector field), and principal curves are derived here in the space of collective variables from the transition path theory, and they are applicable to time series obtained from molecular dynamics simulations of high-dimensional systems, i.e., the position coordinates as a function of time in the transition path ensemble. The connection of the transition path theory is made to a density-weighted average flux, a quantity proposed in a previous work to appraise the relevance of a coordinate to the reaction coordinate [Li, W. J. Chem. Phys. 2022, 156, 054117]. Most importantly, as an extension of the existing quantities, time-lagged quantities such as flux and current are also proposed. The main insights and objects provided by these time-lagged quantities are illustrated in the application to the alanine peptide in vacuum.

A Small Molecule Stabilizes the Disordered Native State of the Alzheimer's Aβ Peptide

The stabilization of native states of proteins is a powerful drug discovery strategy. It is still unclear, however, whether this approach can be applied to intrinsically disordered proteins. Here, we report a small molecule that stabilizes the native state of the Aβ42 peptide, an intrinsically disordered protein fragment associated with Alzheimer's disease. We show that this stabilization takes place by a disordered binding mechanism, in which both the small molecule and the Aβ42 peptide remain disordered. This disordered binding mechanism involves enthalpically favorable local π-stacking interactions coupled with entropically advantageous global effects. These results indicate that small molecules can stabilize disordered proteins in their native states through transient non-specific interactions that provide enthalpic gain while simultaneously increasing the conformational entropy of the proteins.

Nanoparticle cluster formation mechanisms elucidated via Markov state modeling: Attraction range effects, aggregation pathways, and counterintuitive transition rates

Nanoparticle clusters are promising candidates for developing functional materials. However, it is still a challenging task to fabricate them in a predictable and controllable way, which requires investigation of the possible mechanisms underlying cluster formation at the nanoscale. By constructing Markov state models (MSMs) at the microstate level, we find that for highly dispersed particles to form a highly aggregated cluster, there are multiple coexisting pathways, which correspond to direct aggregation, or pathways that need to pass through partially aggregated, intermediate states. Varying the range of attraction between nanoparticles is found to significantly affect pathways. As the attraction range becomes narrower, compared to direct aggregation, some pathways that need to pass through partially aggregated intermediate states become more competitive. In addition, from MSMs constructed at the macrostate level, the aggregation rate is found to be counterintuitively lower with a lower free-energy barrier, which is also discussed.

Dynamics of peptide loading into major histocompatibility complex class I molecules chaperoned by TAPBPR

Major histocompatibility complex class I (MHC-I) molecules display antigenic peptides on the cell surface for T cell receptor scanning, thereby activating the immune response. Peptide loading into MHC-I molecules is thus a critical step during the antigen presentation process. Chaperone TAP-binding protein related (TAPBPR) plays a critical role in promoting high-affinity peptide loading into MHC-I, by discriminating against the low-affinity ones. However, the complete peptide loading dynamics into TAPBPR-bound MHC-I is still elusive. Here, we constructed kinetic network models based on hundreds of short-time MD simulations with an aggregated simulation time of ∼21.7 μs, and revealed, at atomic level, four key intermediate states of one antigenic peptide derived from melanoma-associated MART-1/Melan-A protein during its loading process into TAPBPR-bound MHC-I. We find that the TAPBPR binding at the MHC-I pocket-F can substantially reshape the distant pocket-B via allosteric regulations, which in turn promotes the following peptide N-terminal loading. Intriguingly, the partially loaded peptide could profoundly weaken the TAPBPR-MHC stability, promoting the dissociation of the TAPBPR scoop-loop (SL) region from the pocket-F to a more solvent-exposed conformation. Structural inspections further indicate that the peptide loading could remotely affect the SL binding site through both allosteric perturbations and direct contacts. In addition, another structural motif of TAPBPR, the jack hairpin region, was also found to participate in mediating the peptide editing. Our study sheds light on the detailed molecular mechanisms underlying the peptide loading process into TAPBPR-bound MHC-I and pinpoints the key structural factors responsible for dictating the peptide-loading dynamics.

Computational investigations on target-site searching and recognition mechanisms by thymine DNA glycosylase during DNA repair process

DNA glycosylase, as one member of DNA repair machineries, plays an essential role in correcting mismatched/damaged DNA nucleotides by cleaving the N-glycosidic bond between the sugar and target nucleobase through the base excision repair (BER) pathways. Efficient corrections of these DNA lesions are critical for maintaining genome integrity and preventing premature aging and cancers. The target-site searching/recognition mechanisms and the subsequent conformational dynamics of DNA glycosylase, however, remain challenging to be characterized using experimental techniques. In this review, we summarize our recent studies of sequential structural changes of thymine DNA glycosylase (TDG) during the DNA repair process, achieved mostly by molecular dynamics (MD) simulations. Computational simulations allow us to reveal atomic-level structural dynamics of TDG as it approaches the target-site, and pinpoint the key structural elements responsible for regulating the translocation of TDG along DNA. Subsequently, upon locating the lesions, TDG adopts a base-flipping mechanism to extrude the mispaired nucleobase into the enzyme active-site. The constructed kinetic network model elucidates six metastable states during the base-extrusion process and suggests an active role of TDG in flipping the intrahelical nucleobase. Finally, the molecular mechanism of product release dynamics after catalysis is also summarized. Taken together, we highlight to what extent the computational simulations advance our knowledge and understanding of the molecular mechanism underlying the conformational dynamics of TDG, as well as the limitations of current theoretical work.

Nanobody Paratope Ensembles in Solution Characterized by MD Simulations and NMR

Variable domains of camelid antibodies (so-called nanobodies or VHH) are the smallest antibody fragments that retain complete functionality and therapeutic potential. Understanding of the nanobody-binding interface has become a pre-requisite for rational antibody design and engineering. The nanobody-binding interface consists of up to three hypervariable loops, known as the CDR loops. Here, we structurally and dynamically characterize the conformational diversity of an anti-GFP-binding nanobody by using molecular dynamics simulations in combination with experimentally derived data from nuclear magnetic resonance (NMR) spectroscopy. The NMR data contain both structural and dynamic information resolved at various timescales, which allows an assessment of the quality of protein MD simulations. Thus, in this study, we compared the ensembles for the anti-GFP-binding nanobody obtained from MD simulations with results from NMR. We find excellent agreement of the NOE-derived distance maps obtained from NMR and MD simulations and observe similar conformational spaces for the simulations with and without NOE time-averaged restraints. We also compare the measured and calculated order parameters and find generally good agreement for the motions observed in the ps-ns timescale, in particular for the CDR3 loop. Understanding of the CDR3 loop dynamics is especially critical for nanobodies, as this loop is typically critical for antigen recognition.

Structural Rearrangement of the Serotonin Transporter Intracellular Gate Induced by Thr276 Phosphorylation

The reuptake of the neurotransmitter serotonin from the synaptic cleft by the serotonin transporter, SERT, is essential for proper neurological signaling. Biochemical studies have shown that Thr276 of transmembrane helix 5 is a site of PKG-mediated SERT phosphorylation, which has been proposed to shift the SERT conformational equilibria to promote inward-facing states, thus enhancing 5-HT transport. Recent structural and simulation studies have provided insights into the conformation transitions during substrate transport but have not shed light on SERT regulation via post-translational modifications. Using molecular dynamics simulations and Markov state models, we investigate how Thr276 phosphorylation impacts the SERT mechanism and its role in enhancing transporter stability and function. Our simulations show that Thr276 phosphorylation alters the hydrogen-bonding network involving residues on transmembrane helix 5. This in turn decreases the free energy barriers for SERT to transition to the inward-facing state, thus facilitating 5-HT import. The results provide atomistic insights into in vivo SERT regulation and can be extended to other pharmacologically important transporters in the solute carrier family.

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.

From Data to Knowledge: Systematic Review of Tools for Automatic Analysis of Molecular Dynamics Output

An increasing number of crystal structures available on one side, and the boost of computational power available for computer-aided drug design tasks on the other, have caused that the structure-based drug design tools are intensively used in the drug development pipelines. Docking and molecular dynamics simulations, key representatives of the structure-based approaches, provide detailed information about the potential interaction of a ligand with a target receptor. However, at the same time, they require a three-dimensional structure of a protein and a relatively high amount of computational resources. Nowadays, as both docking and molecular dynamics are much more extensively used, the amount of data output from these procedures is also growing. Therefore, there are also more and more approaches that facilitate the analysis and interpretation of the results of structure-based tools. In this review, we will comprehensively summarize approaches for handling molecular dynamics simulations output. It will cover both statistical and machine-learning-based tools, as well as various forms of depiction of molecular dynamics output.

A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls

The string method with swarms of trajectories (SMwST) is an algorithm that identifies a physically meaningful transition pathway─a one-dimensional curve, embedded within a high-dimensional space of selected collective variables. The SMwST algorithm leans on a series of short, unbiased molecular dynamics simulations spawned at different locations of the discretized path, from whence an average dynamic drift is determined to evolve the string toward an optimal pathway. However conceptually simple in both its theoretical formulation and practical implementation, the SMwST algorithm is computationally intensive and requires a careful choice of parameters for optimal cost-effectiveness in applications to challenging problems in chemistry and biology. In this contribution, the SMwST algorithm is presented in a self-contained manner, discussing with a critical eye its theoretical underpinnings, applicability, inherent limitations, and use in the context of path-following free-energy calculations and their possible extension to kinetics modeling. Through multiple simulations of a prototypical polypeptide, combining the search of the transition pathway and the computation of the potential of mean force along it, several practical aspects of the methodology are examined with the objective of optimizing the computational effort, yet without sacrificing accuracy. In light of the results reported here, we propose some general guidelines aimed at improving the efficiency and reliability of the computed pathways and free-energy profiles underlying the conformational transitions at hand.

Intrinsically disordered proteins: Ensembles at the limits of Anfinsen's dogma

Intrinsically disordered proteins (IDPs) are proteins that lack rigid 3D structure. Hence, they are often misconceived to present a challenge to Anfinsen's dogma. However, IDPs exist as ensembles that sample a quasi-continuum of rapidly interconverting conformations and, as such, may represent proteins at the extreme limit of the Anfinsen postulate. IDPs play important biological roles and are key components of the cellular protein interaction network (PIN). Many IDPs can interconvert between disordered and ordered states as they bind to appropriate partners. Conformational dynamics of IDPs contribute to conformational noise in the cell. Thus, the dysregulation of IDPs contributes to increased noise and "promiscuous" interactions. This leads to PIN rewiring to output an appropriate response underscoring the critical role of IDPs in cellular decision making. Nonetheless, IDPs are not easily tractable experimentally. Furthermore, in the absence of a reference conformation, discerning the energy landscape representation of the weakly funneled IDPs in terms of reaction coordinates is challenging. To understand conformational dynamics in real time and decipher how IDPs recognize multiple binding partners with high specificity, several sophisticated knowledge-based and physics-based in silico sampling techniques have been developed. Here, using specific examples, we highlight recent advances in energy landscape visualization and molecular dynamics simulations to discern conformational dynamics and discuss how the conformational preferences of IDPs modulate their function, especially in phenotypic switching. Finally, we discuss recent progress in identifying small molecules targeting IDPs underscoring the potential therapeutic value of IDPs. Understanding structure and function of IDPs can not only provide new insight on cellular decision making but may also help to refine and extend Anfinsen's structure/function paradigm.

Deeptime: a Python library for machine learning dynamical models from time series data

Generation and analysis of time-series data is relevant to many quantitative fields ranging from economics to fluid mechanics. In the physical sciences, structures such as metastable and coherent sets, slow relaxation processes, collective variables, dominant transition pathways or manifolds and channels of probability flow can be of great importance for understanding and characterizing the kinetic, thermodynamic and mechanistic properties of the system. Deeptime is a general purpose Python library offering various tools to estimate dynamical models based on time-series data including conventional linear learning methods, such as Markov state models (MSMs), Hidden Markov Models and Koopman models, as well as kernel and deep learning approaches such as VAMPnets and deep MSMs. The library is largely compatible with scikit-learn, having a range of Estimator classes for these different models, but in contrast to scikit-learn also provides deep Model classes, e.g. in the case of an MSM, which provide a multitude of analysis methods to compute interesting thermodynamic, kinetic and dynamical quantities, such as free energies, relaxation times and transition paths. The library is designed for ease of use but also easily maintainable and extensible code. In this paper we introduce the main features and structure of the deeptime software. Deeptime can be found under https://deeptime-ml.github.io/.

Deciphering the Effect of Lysine Acetylation on the Misfolding and Aggregation of Human Tau Fragment 171IPAKTPPAPK180 Using Molecular Dynamic Simulation and the Markov State Model

The formation of neurofibrillary tangles (NFT) with β-sheet-rich structure caused by abnormal aggregation of misfolded microtubule-associated protein Tau is a hallmark of tauopathies, including Alzheimer’s Disease. It has been reported that acetylation, especially K174 located in the proline-rich region, can largely promote Tau aggregation. So far, the mechanism of the abnormal acetylation of Tau that affects its misfolding and aggregation is still unclear. Therefore, revealing the effect of acetylation on Tau aggregation could help elucidate the pathogenic mechanism of tauopathies. In this study, molecular dynamics simulation combined with multiple computational analytical methods were performed to reveal the effect of K174 acetylation on the spontaneous aggregation of Tau peptide ¹⁷¹IPAKTPPAPK¹⁸⁰, and the dimerization mechanism as an early stage of the spontaneous aggregation was further specifically analyzed by Markov state model (MSM) analysis. The results showed that both the actual acetylation and the mutation mimicking the acetylated state at K174 induced the aggregation of the studied Tau fragment; however, the effect of actual acetylation on the aggregation was more pronounced. In addition, acetylated K174 plays a major contributing role in forming and stabilizing the antiparallel β-sheet dimer by forming several hydrogen bonds and side chain van der Waals interactions with residues I171, P172, A173 and T175 of the corresponding chain. In brief, this study uncovered the underlying mechanism of Tau peptide aggregation in response to the lysine K174 acetylation, which can deepen our understanding on the pathogenesis of tauopathies.

Molecular Dynamics Simulations of Protein Aggregation: Protocols for Simulation Setup and Analysis with Markov State Models and Transition Networks

Protein disorder and aggregation play significant roles in the pathogenesis of numerous neurodegenerative diseases, such as Alzheimer's and Parkinson's diseases. The end products of the aggregation process in these diseases are highly structured amyloid fibrils. Though in most cases, small, soluble oligomers formed during amyloid aggregation are the toxic species. A full understanding of the physicochemical forces that drive protein aggregation is thus required if one aims for the rational design of drugs targeting the formation of amyloid oligomers. Among a multitude of biophysical and biochemical techniques that are employed for studying protein aggregation, molecular dynamics (MD) simulations at the atomic level provide the highest temporal and spatial resolution of this process, capturing key steps during the formation of amyloid oligomers. Here we provide a step-by-step guide for setting up, running, and analyzing MD simulations of aggregating peptides using GROMACS. For the analysis, we provide the scripts that were developed in our lab, which allow to determine the oligomer size and inter-peptide contacts that drive the aggregation process. Moreover, we explain and provide the tools to derive Markov state models and transition networks from MD data of peptide aggregation.

Bispecific antibodies-effects of point mutations on CH3-CH3 interface stability

A new format of therapeutic proteins is bispecific antibodies, in which two different heavy chains heterodimerize to obtain two different binding sites. Therefore, it is crucial to understand and optimize the third constant domain (CH3-CH3) interface to favor heterodimerization over homodimerization, and to preserve the physicochemical properties, as thermal stability. Here, we use molecular dynamics simulations to investigate the dissociation process of 19 CH3-CH3 crystal structures that differ from each other in few point mutations. We describe the dissociation of the dimeric interface as a two-steps mechanism. As confirmed by a Markov state model, apart from the bound and the dissociated state, we observe an additional intermediate state, which corresponds to an encounter complex. The analysis of the interdomain contacts reveals key residues that stabilize the interface. We expect that our results will improve the understanding of the CH3-CH3 interface interactions and thus advance the developability and design of new antibodies formats.

Ligand Binding Path Sampling Based on Parallel Cascade Selection Molecular Dynamics: LB-PaCS-MD

Parallel cascade selection molecular dynamics (PaCS-MD) is a rare-event sampling method that generates transition pathways between a reactant and product. To sample the transition pathways, PaCS-MD repeats short-time MD simulations from important configurations as conformational resampling cycles. In this study, PaCS-MD was extended to sample ligand binding pathways toward a target protein, which is referred to as LB-PaCS-MD. In a ligand-concentrated environment, where multiple ligand copies are randomly arranged around the target protein, LB-PaCS-MD allows for the frequent sampling of ligand binding pathways. To select the important configurations, we specified the center of mass (COM) distance between each ligand and the relevant binding site of the target protein, where snapshots generated by the short-time MD simulations were ranked by their COM distance values. From each cycle, snapshots with smaller COM distance values were selected as the important configurations to be resampled using the short-time MD simulations. By repeating conformational resampling cycles, the COM distance values gradually decreased and converged to constants, meaning that a set of ligand binding pathways toward the target protein was sampled by LB-PaCS-MD. To demonstrate relative efficiency, LB-PaCS-MD was applied to several proteins, and their ligand binding pathways were sampled more frequently than conventional MD simulations.

gmxapi: A GROMACS-native Python interface for molecular dynamics with ensemble and plugin support

Gmxapi provides an integrated, native Python API for both standard and advanced molecular dynamics simulations in GROMACS. The Python interface permits multiple levels of integration with the core GROMACS libraries, and legacy support is provided via an interface that mimics the command-line syntax, so that all GROMACS commands are fully available. Gmxapi has been officially supported since the GROMACS 2019 release and is enabled by default in current versions of the software. Here we describe gmxapi 0.3 and later. Beyond simply wrapping GROMACS library operations, the API permits several advanced operations that are not feasible using the prior command-line interface. First, the API allows custom user plugin code within the molecular dynamics force calculations, so users can execute custom algorithms without modifying the GROMACS source. Second, the Python interface allows tasks to be dynamically defined, so high-level algorithms for molecular dynamics simulation and analysis can be coordinated with loop and conditional operations. Gmxapi makes GROMACS more accessible to custom Python scripting while also providing support for high-level data-flow simulation algorithms that were previously feasible only in external packages.

Searching for Low Probability Opening Events in a DNA Sliding Clamp

The β subunit of E. coli DNA polymererase III is a DNA sliding clamp associated with increasing the processivity of DNA synthesis. In its free form, it is a circular homodimer structure that can accomodate double-stranded DNA in a nonspecific manner. An open state of the clamp must be accessible before loading the DNA. The opening mechanism is still a matter of debate, as is the effect of bound DNA on opening/closing kinetics. We use a combination of atomistic, coarse-grained, and enhanced sampling strategies in both explicit and implicit solvents to identify opening events in the sliding clamp. Such simulations of large nucleic acid and their complexes are becoming available and are being driven by improvements in force fields and the creation of faster computers. Different models support alternative opening mechanisms, either through an in-plane or out-of-plane opening event. We further note some of the current limitations, despite advances, in modeling these highly charged systems with implicit solvent.

WESTPA 2.0: High-Performance Upgrades for Weighted Ensemble Simulations and Analysis of Longer-Timescale Applications

The weighted ensemble (WE) family of methods is one of several statistical mechanics-based path sampling strategies that can provide estimates of key observables (rate constants and pathways) using a fraction of the time required by direct simulation methods such as molecular dynamics or discrete-state stochastic algorithms. WE methods oversee numerous parallel trajectories using intermittent overhead operations at fixed time intervals, enabling facile interoperability with any dynamics engine. Here, we report on the major upgrades to the WESTPA software package, an open-source, high-performance framework that implements both basic and recently developed WE methods. These upgrades offer substantial improvements over traditional WE methods. The key features of the new WESTPA 2.0 software enhance the efficiency and ease of use: an adaptive binning scheme for more efficient surmounting of large free energy barriers, streamlined handling of large simulation data sets, exponentially improved analysis of kinetics, and developer-friendly tools for creating new WE methods, including a Python API and resampler module for implementing both binned and "binless" WE strategies.

Distinct Binding Mechanisms for Allosteric Sodium Ion in Cannabinoid Receptors

The therapeutic potential of cannabinoid receptors is not fully explored due to psychoactive side effects and lack of selectivity associated with orthosteric ligands. Allosteric modulators have the potential to become selective therapeutics for cannabinoid receptors. Biochemical experiments have shown the effects of the allosteric Na⁺ binding on cannabinoid receptor activity. However, the Na⁺ coordination site and binding pathway are still unknown. Here, we perform molecular dynamic simulations to explore Na⁺ binding in the cannabinoid receptors, CB₁ and CB₂. Simulations reveal that Na⁺ binds to the primary binding site from different extracellular sites for CB₁ and CB₂. A distinct secondary Na⁺ coordination site is identified in CB₁ that is not present in CB₂. Furthermore, simulations also show that intracellular Na⁺ could bind to the Na⁺ binding site in CB₁. Constructed Markov state models show that the standard free energy of Na⁺ binding is similar to the previously calculated free energy for other class A GPCRs.

Broadly neutralizing antibodies target a haemagglutinin anchor epitope

Broadly neutralizing antibodies that target epitopes of haemagglutinin on the influenza virus have the potential to provide near universal protection against influenza virus infection1. However, viral mutants that escape broadly neutralizing antibodies have been reported2,3. The identification of broadly neutralizing antibody classes that can neutralize viral escape mutants is critical for universal influenza virus vaccine design. Here we report a distinct class of broadly neutralizing antibodies that target a discrete membrane-proximal anchor epitope of the haemagglutinin stalk domain. Anchor epitope-targeting antibodies are broadly neutralizing across H1 viruses and can cross-react with H2 and H5 viruses that are a pandemic threat. Antibodies that target this anchor epitope utilize a highly restricted repertoire, which encodes two public binding motifs that make extensive contacts with conserved residues in the fusion peptide. Moreover, anchor epitope-targeting B cells are common in the human memory B cell repertoire and were recalled in humans by an oil-in-water adjuvanted chimeric haemagglutinin vaccine4,5, which is a potential universal influenza virus vaccine. To maximize protection against seasonal and pandemic influenza viruses, vaccines should aim to boost this previously untapped source of broadly neutralizing antibodies that are widespread in the human memory B cell pool.

Decisive role of water and protein dynamics in residence time of p38α MAP kinase inhibitors

Target residence time plays a crucial role in the pharmacological activity of small molecule inhibitors. Little is known, however, about the underlying causes of inhibitor residence time at the molecular level, which complicates drug optimization processes. Here, we employ all-atom molecular dynamics simulations (~400 μs in total) to gain insight into the binding modes of two structurally similar p38α MAPK inhibitors (type I and type I½) with short and long residence times that otherwise show nearly identical inhibitory activities in the low nanomolar IC₅₀ range. Our results highlight the importance of protein conformational stability and solvent exposure, buried surface area of the ligand and binding site resolvation energy for residence time. These findings are further confirmed by simulations with a structurally diverse short residence time inhibitor SB203580. In summary, our data provide guidance in compound design when aiming for inhibitors with improved target residence time.

The molecular determinants of the residence time of a small molecule inhibitor at its target protein are not well understood. Here, Pantsar et al. show that the target protein’s conformational stability and solvent exposure are key factors governing the target residence time of kinase inhibitors.

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.

Delineating the activation mechanism and conformational landscape of a class B G protein-coupled receptor glucagon receptor

Class B G protein-coupled receptors (GPCRs) are important targets in the treatment of metabolic syndrome and diabetes. Although multiple structures of class B GPCRs-G protein complexes have been elucidated, the detailed activation mechanism of the receptors remains unclear. Here, we combine Gaussian accelerated molecular dynamics simulations and Markov state models (MSM) to investigate the activation mechanism of a canonical class B GPCR, human glucagon receptor-GCGR, including the negative allosteric modulator-bound inactive state, the agonist glucagon-bound active state, and both glucagon- and Gs-bound fully active state. The free-energy landscapes of GCGR show the conformational ensemble consisting of three activation-associated states: inactive, active, and fully active. The structural analysis indicates the high dynamics of GCGR upon glucagon binding with both active and inactive conformations in the ensemble. Significantly, the H8 and TM6 exhibits distinct features from the inactive to the active states. The additional simulations demonstrate the role of H8 in the recruitment of Gs. Gs binding presents a crucial function of stabilizing the glucagon binding site and MSM highlights the absolute requirement of Gs to help the GCGR reach the fully active state. Together, our results reveal the detailed activation mechanism of GCGR from the view of conformational dynamics.

Enhanced sampling without borders: on global biasing functions and how to reweight them

Molecular dynamics (MD) simulations are a powerful tool to follow the time evolution of biomolecular motions in atomistic resolution. However, the high computational demand of these simulations limits the timescales of motions that can be observed. To resolve this issue, so called enhanced sampling techniques are developed, which extend conventional MD algorithms to speed up the simulation process. Here, we focus on techniques that apply global biasing functions. We provide a broad overview of established enhanced sampling methods and promising new advances. As the ultimate goal is to retrieve unbiased information from biased ensembles, we also discuss benefits and limitations of common reweighting schemes. In addition to concisely summarizing critical assumptions and implications, we highlight the general application opportunities as well as uncertainties of global enhanced sampling.

RPnet: a reverse-projection-based neural network for coarse-graining metastable conformational states for protein dynamics

The Markov State Model (MSM) is a powerful tool for modeling long timescale dynamics based on numerous short molecular dynamics (MD) simulation trajectories, which makes it a useful tool for elucidating the conformational changes of biological macromolecules. By partitioning the phase space into discretized states and estimating the probabilities of inter-state transitions based on short MD trajectories, one can construct a kinetic network model that could be used to extrapolate long-timescale kinetics if the Markovian condition is met. However, meeting the Markovian condition often requires hundreds or even thousands of states (microstates), which greatly hinders the comprehension of the conformational dynamics of complex biomolecules. Kinetic lumping algorithms can coarse grain numerous microstates into a handful of metastable states (macrostates), which would greatly facilitate the elucidation of biological mechanisms. In this work, we have developed a reverse-projection-based neural network (RPnet) to lump microstates into macrostates, by making use of a physics-based loss function that is based on the projection operator framework of conformational dynamics. By recognizing that microstate and macrostate transition modes can be related through a projection process, we have developed a reverse-projection scheme to directly compare the microstate and macrostate dynamics. Based on this reverse-projection scheme, we designed a loss function that allows the effective assessment of the quality of a given kinetic lumping. We then make use of a neural network to efficiently minimize this loss function to obtain an optimized set of macrostates. We have demonstrated the power of our RPnet in analyzing the dynamics of a numerical 2D potential, alanine dipeptide, and the clamp opening of an RNA polymerase. In all these systems, we have illustrated that our method could yield comparable or better results than competing methods in terms of state partitioning and reproduction of slow dynamics. We expect that our RPnet holds promise in analyzing the conformational dynamics of biological macromolecules.

Unprotected peptide macrocyclization and stapling via a fluorine-thiol displacement reaction

We report the discovery of a facile peptide macrocyclization and stapling strategy based on a fluorine thiol displacement reaction (FTDR), which renders a class of peptide analogues with enhanced stability, affinity, cellular uptake, and inhibition of cancer cells. This approach enabled selective modification of the orthogonal fluoroacetamide side chains in unprotected peptides in the presence of intrinsic cysteines. The identified benzenedimethanethiol linker greatly promoted the alpha helicity of a variety of peptide substrates, as corroborated by molecular dynamics simulations. The cellular uptake of benzenedimethanethiol stapled peptides appeared to be universally enhanced compared to the classic ring-closing metathesis (RCM) stapled peptides. Pilot mechanism studies suggested that the uptake of FTDR-stapled peptides may involve multiple endocytosis pathways in a distinct pattern in comparison to peptides stapled by RCM. Consistent with the improved cell permeability, the FTDR-stapled lead Axin and p53 peptide analogues demonstrated enhanced inhibition of cancer cells over the RCM-stapled analogues and the unstapled peptides.

Strategies capable of stapling unprotected peptides in a straightforward, chemoselective, and clean manner, as well as promoting cellular uptake are of great interest. Here the authors report a peptide macrocyclization and stapling strategy which satisfies those criteria, based on a fluorine thiol displacement reaction.

Transient exposure of a buried phosphorylation site in an autoinhibited protein

Autoinhibition is a mechanism used to regulate protein function, often by making functional sites inaccessible through the interaction with a cis-acting inhibitory domain. Such autoinhibitory domains often display a substantial degree of structural disorder when unbound, and only become structured in the inhibited state. These conformational dynamics make it difficult to study the structural origin of regulation, including effects of regulatory post-translational modifications. Here, we study the autoinhibition of the Dbl Homology domain in the protein Vav1 by the so-called acidic inhibitory domain. We use molecular simulations to study the process by which a mostly unstructured inhibitory domain folds upon binding and how transient exposure of a key buried tyrosine residue makes it accessible for phosphorylation. We show that the inhibitory domain, which forms a helix in the bound and inhibited stated, samples helical structures already before binding and that binding occurs via a molten-globule-like intermediate state. Together, our results shed light on key interactions that enable the inhibitory domain to sample a finely tuned equilibrium between an inhibited and a kinase-accessible state.

Long-Timescale Simulations Revealed Critical Non-Conserved Residues of Phosphodiesterases Affecting Selectivity of BAY60-7550

A major obstacle of the selective inhibitor design for specific human phosphodiesterase (PDE) is that highly conserved catalytic pockets are difficult to be distinguished by inhibitor molecules. To overcome this, a feasible path is to understand the molecular determinants underlying the selectivity of current inhibitors. BAY60-7550 (BAY for short; IC₅₀ = 4.7 nM) is a highly selective inhibitor targeting PDE2A which is a dual-specificity PDE and an attractive target for therapeutic intervention of the central nervous system (CNS) disorders. Recent studies suggest that molecular determinants may be in binding processes of BAY. However, a detailed understanding of these processes are still lacking. To explore these processes, High-Throughput Molecular Dynamics (HTMD) simulations were performed to reproduce the spontaneous association of BAY with catalytic pockets of 4 PDE isoforms; Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD) simulations were performed to reproduce the unbinding-rebinding processes of FKG and 10.13039/100016266MC2, two pyrazolopyrimidinone PDE2A selective inhibitors, in the PDE2A system. The produced molecular trajectories were analyzed by the Markov state model (MSM) and the molecular mechanics/generalized Born surface area (MM/GBSA). The results showed that the non-covalent interactions between the non-conserved residues and BAY, especially the hydrogen bonds, determined the unique binding pathways of BAY on the surface of PDE2A. These pathways were different from those of BAY on the surface of the other three PDE isoforms and the binding pathways of the other two PDE2A inhibitors in PDE2A systems. These differences were ultimately reflected in the high selectivity of this inhibitor for PDE2A. As a result, this study demonstrates the critical role of the binding processes in the selectivity of BAY, and also identifies the key non-conserved residues affecting the binding processes of BAY. Thus, this study provides a new perspective and data support for the further development of BAY-derived inhibitors targeting PDE2A.

Analysis of Molecular Dynamics Simulations of Protein Folding

Unbiased molecular dynamics simulations of proteins can now capture spontaneous folding events. This provides a wealth of data reflecting information on folding mechanism, but raises the challenge of interpreting it in a meaningful way. Here, I describe how such simulations can be used to identify reactive states and reaction coordinates for describing folding, and how folding dynamics can be captured by projection onto those coordinates. Methods are described for quantifying the interactions important for defining the folding mechanism, and for comparison of simulations with experimental mechanistic probes, such as ϕ-values.

CDR loop interactions can determine heavy and light chain pairing preferences in bispecific antibodies

As the current biotherapeutic market is dominated by antibodies, the design of different antibody formats, like bispecific antibodies, is critical to the advancement of the field. In contrast to monovalent antibodies, which consist of two identical antigen-binding sites, bispecific antibodies can target two different epitopes by containing two different antigen-binding sites. Thus, the rise of new formats as successful therapeutics has reignited the interest in advancing and facilitating the efficient production of bispecific antibodies. Here, we investigate the influence of point mutations in the antigen-binding site, the paratope, on heavy and light chain pairing preferences by using molecular dynamics simulations. In agreement with experiments, we find that specific residues in the antibody variable domain (Fv), i.e., the complementarity-determining region (CDR) L3 and H3 loops, determine heavy and light chain pairing preferences. Excitingly, we observe substantial population shifts in CDR-H3 and CDR-L3 loop conformations in solution accompanied by a decrease in bispecific IgG yield. These conformational changes in the CDR3 loops induced by point mutations also influence all other CDR loop conformations and consequentially result in different CDR loop states in solution. However, besides their effect on the obtained CDR loop ensembles, point mutations also lead to distinct interaction patterns in the VH-VL interface. By comparing the interaction patterns among all investigated variants, we observe specific contacts in the interface that drive heavy and light chain pairing. Thus, these findings have broad implications in the field of antibody engineering and design because they provide a mechanistic understanding of antibody interfaces, by identifying critical factors driving the pairing preferences, and thus can help to advance the design of bispecific antibodies.

Chasing Collective Variables Using Autoencoders and Biased Trajectories

Free energy biasing methods have proven to be powerful tools to accelerate the simulation of important conformational changes of molecules by modifying the sampling measure. However, most of these methods rely on the prior knowledge of low-dimensional slow degrees of freedom, i.e., collective variables (CVs). Alternatively, such CVs can be identified using machine learning (ML) and dimensionality reduction algorithms. In this context, approaches where the CVs are learned in an iterative way using adaptive biasing have been proposed: at each iteration, the learned CV is used to perform free energy adaptive biasing to generate new data and learn a new CV. In this paper, we introduce a new iterative method involving CV learning with autoencoders: Free Energy Biasing and Iterative Learning with AutoEncoders (FEBILAE). Our method includes a reweighting scheme to ensure that the learning model optimizes the same loss at each iteration and achieves CV convergence. Using the alanine dipeptide system and the solvated chignolin mini-protein system as examples, we present results of our algorithm using the extended adaptive biasing force as the free energy adaptive biasing method.

Harnessing Reversed Allosteric Communication: A Novel Strategy for Allosteric Drug Discovery

Allostery is a fundamental and extensive mechanism of intramolecular signal transmission. Allosteric drugs possess several unique pharmacological advantages over traditional orthosteric drugs, including greater selectivity, better physicochemical properties, and lower off-target toxicity. However, owing to the complexity of allosteric regulation, experimental approaches for the development of allosteric modulators are traditionally serendipitous. Recently, the reversed allosteric communication theory has been proposed, providing a feasible tool for the unbiased detection of allosteric sites. Herein, we review the latest research on the reversed allosteric communication effect using the examples of sirtuin 6, epidermal growth factor receptor, 3-phosphoinositide-dependent protein kinase 1, and Related to A and C kinases (RAC) serine/threonine protein kinase B and recapitulate the methodologies of reversed allosteric communication strategy. The novel reversed allosteric communication strategy greatly expands the horizon of allosteric site identification and allosteric mechanism exploration and is expected to accelerate an end-to-end framework for drug discovery.

Modeling of biomolecular machines in non-equilibrium steady states

Numerical computations have become a pillar of all modern quantitative sciences. Any computation involves modeling-even if often this step is not made explicit-and any model has to neglect details while still being physically accurate. Equilibrium statistical mechanics guides both the development of models and numerical methods for dynamics obeying detailed balance. For systems driven away from thermal equilibrium, such a universal theoretical framework is missing. For a restricted class of driven systems governed by Markov dynamics and local detailed balance, stochastic thermodynamics has evolved to fill this gap and to provide fundamental constraints and guiding principles. The next step is to advance stochastic thermodynamics from simple model systems to complex systems with tens of thousands or even millions of degrees of freedom. Biomolecules operating in the presence of chemical gradients and mechanical forces are a prime example for this challenge. In this Perspective, we give an introduction to isothermal stochastic thermodynamics geared toward the systematic multiscale modeling of the conformational dynamics of biomolecular and synthetic machines, and we outline some of the open challenges.

Hidden allosteric sites and De-Novo drug design

INTRODUCTION

Hidden allosteric sites are not visible in apo-crystal structures, but they may be visible in holo-structures when a certain ligand binds and maintains the ligand intended conformation. Several computational and experimental techniques have been used to investigate these hidden sites but identifying them remains a challenge.

AREAS COVERED

This review provides a summary of the many theoretical approaches for predicting hidden allosteric sites in disease-related proteins. Furthermore, promising cases have been thoroughly examined to reveal the hidden allosteric site and its modulator.

EXPERT OPINION

In the recent past, with the development in scientific techniques and bioinformatics tools, the number of drug targets for complex human diseases has significantly increased but unfortunately most of these targets are undruggable due to several reasons. Alternative strategies such as finding cryptic (hidden) allosteric sites are an attractive approach for exploitation of the discovery of new targets. These hidden sites are difficult to recognize compared to allosteric sites, mainly due to a lack of visibility in the crystal structure. In our opinion, after many years of development, MD simulations are finally becoming successful for obtaining a detailed molecular description of drug-target interaction.