Markov State Models Research Articles

Uncertainty in Markov state models for small macromolecules

Markov state models (MSMs) are a tool to describe and analyse protein dynamics. They are useful in particular to determine characteristic timescales of protein motion. For larger biomolecules such as proteins it is challenging to obtain sufficient sampling for the timescales to converge, which is thus a frequent concern. There are, however, also several other sources of uncertainty, such as choice of lag time and the number of dimensions in the dimension reduction preprocessing step (e.g. time-lagged independent component analysis), the number of Markov states, or choice of the MSM lag time.

Characterizing blebbistatin pocket conformational dynamics with Markov state models

Allosteric modulators of protein function hold great promise as therapeutics. As they do not target protein active sites, these compounds are much more likely to achieve isoform specificity, a key challenge in drug design especially when targeting closely related families of proteins like G-protein-coupled receptors, kinases, or myosins. There are many known allosteric modulators that selectively bind a subset of isoforms within a family of proteins. Currently, compound specificity is attributed to sequence differences within binding sites discovered through crystallization experiments.

Exploring the myosin active/auto-inhibited state equilibrium by Markov state modeling

Myosins are a ubiquitous class of molecular motors found in eukaryotes. Myosins play specialized roles in virtually every human cell. In fact, there are over 38 different myosin genes in the human genome. Due to their ubiquity, myosins have extensive clinical relevance as mutations in these proteins are often deleterious. For example, mutations in human β-cardiac myosin (MYH7) result in a spectrum of heart diseases called cardiomyopathies. Hypertrophic cardiomyopathy (HCM) is one of the leading causes of sudden cardiac death in the United States.

Insights of hSERT transport mechanism revealed by Markov state models

The human serotonin transporter (hSERT) belongs to the solute carrier 6 (SLC6) family and is expressed in neurons. hSERT is essential for regulating serotonin levels in the synaptic cleft by transporting the neurotransmitter into the pre-synaptic neuron. Many psychopathological diseases are linked to serotonergic regulation, and pharmacological treatments frequently target hSERT. However, available medications suffer from limited efficacy and severe side effects. Essential for overcoming these limitations is a comprehensive understanding of the transport mechanism at the molecular level.

The substrate import mechanism of the human serotonin transporter

The serotonin transporter (SERT) initiates the reuptake of extracellular serotonin in the synapse to terminate neurotransmission. The cryogenic electron microscopy structures of SERT bound to ibogaine and the physiological substrate serotonin resolved in different states have provided a glimpse of the functional conformations at atomistic resolution. However, the conformational dynamics and structural transitions to intermediate states are not fully understood. Furthermore, the molecular basis of how serotonin is recognized and transported remains unclear. In this study, we performed unbiased microsecond-long simulations of the human SERT to investigate the structural dynamics to various intermediate states and elucidated the complete substrate import pathway. Using Markov state models, we characterized a sequential order of conformational-driven ion-coupled substrate binding and transport events and calculated the free energy barriers of conformation transitions associated with the import mechanism. We find that the transition from the occluded to inward-facing state is the rate-limiting step for substrate import and that the substrate decreases the free energy barriers to achieve the inward-facing state. Our study provides insights on the molecular basis of dynamics-driven ion-substrate recognition and transport of SERT that can serve as a model for other closely related neurotransmitter transporters.

Whole virion simulations unveiling the hemagglutinin-neuraminidase interplay

Despite decades of research, the mechanism of how influenza virions release from a cell is still murky. We know the membrane protein neuraminidase must cleave a bond, releasing another membrane protein, hemagglutinin, from the cell surface. We examine this interplay between hemagglutinin and neuraminidase through simulations. We created an influenza virion and ran all-atom molecular dynamics. We analyzed the ensuing protein motions through Markov state models and protein contacts to determine how these two membrane proteins interact to free influenza virions from host cells.

Membrane-like interfaces can facilitate the conformational closure of cyclic peptides for passive membrane permeability

Cyclic peptides have the potential to vastly extend the scope of druggable proteins and to lead to new therapeutics for currently untreatable diseases. However, the design of cyclic peptides with oral bioavailability is challenging and often only achieved via a tedious trial-and-error process. We use molecular dynamics (MD) simulations and Markov state modeling of cyclic peptides to rationalize the relationship between their structure and passive membrane permeability. Previous work has shown that an important determinant for permeability is a peptide's conformational behavior in different environments.

Optimization of hydration sites in plant hormone receptors for agrochemical design

Phytohormones are endogenous small molecules derived from common metabolic pathways that regulate key processes in plants, such as defense, development, and growth. Plant hormones are specifically recognized by the binding pocket in their receptors. Phytohormone receptors' agonists and antagonists have the potential to become useful agrochemicals. For this reason, it is necessary to characterize the factors determining the binding affinity and selectivity of these ligands to their receptors. In this work, we employed extensive molecular dynamics simulations coupled with Markov state models to study the binding of eight major classes of plant hormones (auxin, jasmonate, gibberellin, strigolactone, brassinosteroid, cytokinin, salicylic acid, and abscisic acid) to their receptors.

Iterative steady state restarting to accelerate weighted ensemble convergence

Although the weighted ensemble (WE) algorithm provides an efficient, unbiased framework for rare-event sampling in molecular dynamics (MD) simulations, WE convergence timescales may still be limiting for complex systems. In particular, unbiased estimates of key observables require relaxation to steady state, which can be slow in standard WE simulations. Previously, history-augmented Markov state models (haMSMs) were demonstrated to provide estimates of steady state properties from transient, relaxation-phase WE data.

Trajectory data for the publication "Investigating the bioactive conformation of Angiotensin II using Markov State Modeling revisited with web-scale clustering"

Trajectory data for the publication "Investigating the bioactive conformation of Angiotensin II using Markov State Modeling revisited with web-scale clustering"

Prediction of Protein-Protein Interactions Between Alsin DH/PH and Rac1 and Resulting Protein Dynamics.

Alsin is a protein of 1,657 amino acids known for its crucial role in vesicular trafficking in neurons thanks to its ability to interact with two guanosine triphosphatases, Rac1 and Rab5. Evidence suggests that Rac1 can bind Alsin central region, composed by a Dbl Homology (DH) domain followed by a Pleckstrin Homology (PH) domain, leading to Alsin relocalization. However, Alsin three-dimensional structure and its relationship with known biological functions of this protein are still unknown. In this work, a homology model of the Alsin DH/PH domain was developed and studied through molecular dynamics both in the presence and in the absence of its binding partner, Rac1. Due to different conformations of DH domain, the presence of Rac1 seems to stabilize an open state of the protein, while the absence of its binding partner results in closed conformations. Furthermore, Rac1 interaction was able to reduce the fluctuations in the second conserved region of DH motif, which may be involved in the formation of a homodimer. Moreover, the dynamics of DH/PH was described through a Markov State Model to study the pathways linking the open and closed states. In conclusion, this work provided an all-atom model for the DH/PH domain of Alsin protein; moreover, molecular dynamics investigations suggested underlying molecular mechanisms in the signal transduction between Rac1 and Alsin, providing the basis for a deeper understanding of the whole structure–function relationship for Alsin protein.

Divide and Conquer. Pocket-Opening Mixed-Solvent Simulations in the Perspective of Docking Virtual Screening Applications for Drug Discovery.

The existence of a druggable binding pocket is a prerequisite for computational drug-target interaction studies including virtual screening. Retrospective studies have shown that extended sampling methods like Markov State Modeling and mixed-solvent simulations can identify cryptic pockets relevant for drug discovery. Here, we apply a combination of mixed-solvent molecular dynamics (MD) and time-structure independent component analysis (TICA) to four retrospective case studies: NPC2, the CECR2 bromodomain, TEM-1, and MCL-1. We compare previous experimental and computational findings to our results. It is shown that the successful identification of cryptic pockets depends on the system and the cosolvent probes. We used alternative TICA internal features such as the unbiased backbone coordinates or backbone dihedrals versus biased interatomic distances. We found that in the case of NPC2, TEM-1, and MCL-1, the use of unbiased features is able to identify cryptic pockets, although in the case of the CECR2 bromodomain, more specific features are required to properly capture a pocket opening. In the perspective of virtual screening applications, it is shown how docking studies with the parent ligands depend critically on the conformational state of the targets.

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.

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, CB1 and CB2. Simulations reveal that Na+ binds to the primary binding site from different extracellular sites for CB1 and CB2. A distinct secondary Na+ coordination site is identified in CB1 that is not present in CB2. Furthermore, simulations also show that intracellular Na+ could bind to the Na+ binding site in CB1. 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.

Polar/apolar interfaces modulate the conformational behavior of cyclic peptides with impact on their passive membrane permeability.

Cyclic peptides have the potential to vastly extend the scope of druggable proteins and lead to new therapeutics for currently untreatable diseases. However, cyclic peptides often suffer from poor bioavailability. To uncover design principles for permeable cyclic peptides, a promising strategy is to analyze the conformational dynamics of the peptides using molecular dynamics (MD) and Markov state models (MSMs). Previous MD studies have focused on the conformational dynamics in pure aqueous or apolar environments to rationalize membrane permeability. However, during the key steps of the permeation through the membrane, cyclic peptides are exposed to interfaces between polar and apolar regions. Recent studies revealed that these interfaces constitute the free energy minima of the permeation process. Thus, a deeper understanding of the behavior of cyclic peptides at polar/apolar interfaces is desired. Here, we investigate the conformational and kinetic behavior of cyclic decapeptides at a water/chloroform interface using unbiased MD simulations and MSMs. The distinct environments at the interface alter the conformational equilibrium as well as the interconversion kinetics of cyclic peptide conformations. For peptides with low population of the permeable conformation in aqueous solution, the polar/apolar interface facilitates the interconversion to the closed conformation, which is required for membrane permeation. Comparison to unbiased MD simulations with a POPC bilayer reveals that not only the conformations but also the orientations are relevant in a membrane system. These findings allow us to propose a permeability model that includes both ‘prefolding’ and ‘non-prefolding’ cyclic peptides – an extension that can lead to new design considerations for permeable cyclic peptides.

Understanding and design of non-conservative optical matter systems using Markov state models

Non-conservative and permutationally-invariant Markov state models inform understanding and control of self-assembling optical matter systems.

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.

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; IC50 = 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 MC2, 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.

Studying noncovalent or covalent bond problem between smoothened and cholesterol by molecular dynamics simulation and Markov state model.

Smoothened (SMO) is an attractive therapeutic target for the treatment and prevention of several malignant tumors of the nervous system. The crystal structure of SMO shows cholesterol interacts with residue Asp95 via the noncovalent bond. However, some studies indicate that cholesterol covalently binds to residue Asp95 of SMO. To study these contradictory results, we performed molecular dynamics (MD) simulations and Markov state model (MSM) on SMO in complex with noncovalent-bound and covalent-bound cholesterol. The MD simulated results showed that the noncovalent-bound cholesterol was extremely unstable around the position of residue Asp95 of SMO, while the covalent-bound cholesterol could keep the stable connection with residue Asp95 of SMO. The free energy landscape showed that noncovalent-bound cholesterol had more deep energy wells than covalent-bound cholesterol when it dynamically interacted with the extracellular domain of SMO crystal structure. The MSM results showed the noncovalent-bound cholesterol had more dynamic configuration transformation pathways than the covalent-bound cholesterol. These results theoretically revealed cholesterol should have a covalent bond with residue Asp95 if cholesterol could be stable in the near position of residue Asp95 of SMO. Our studies not only elucidate the covalent binding contradictory issue between cholesterol and residue Asp95 of SMO, but also supply helpful information for antagonists design of SMO.

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.