Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling

Revealing the key structural features promoting the helical conformation in algal polysaccharide carrageenan in solution

Carrageenans are industrially important polysaccharides with tunable viscoelastic and gelation properties. The function of polysaccharide depends on its conformation and chemical composition. However, the solution conformations of carrageenans are highly debated, and the structure-function relationship remains elusive. Here, we have studied the intrinsic conformational behavior of a series of carrageenan hexamers in solution, using extensive all-atom classical MD and enhanced sampling. Our findings comprehensively delineate that carrageenans containing the 3,6-anhydrous bridge (κ-C, ι-C, θ-C, and non-sulfated β-C) adopt compact helical structures as their predominant conformation in solution, whereas carrageenans without the bridge (μ-C, ν-C, and λ-C) remain as extended loosely packed helices; opposing the 'coil-to-helix' paradigm. Glycosidic linkages access a few allowed orientations. We hypothesize that the 3,6-anhydrous bridge, irrespective of carrageenan's sulfation pattern, is essential for stabilizing the helical conformation at the single-chain level. It provides necessary flexibility to the glycosidic linkage to sample conformations close to the experimentally derived helical structure and also prevents the sugar ring flipping. Sulfate groups mainly modify the chain stiffness due to steric and stereo-electronic effects and participate in hydrogen bonding. Such atomistic insights will be helpful for understanding the differential gelation mechanisms of carrageenans and fine-tuning polysaccharide backbone for various industrial applications.

Fractional Extended Diffusion Theory to capture anomalous relaxation from biased/accelerated molecular simulations

Biased and accelerated molecular simulations (BAMS) are widely used tools to observe relevant molecular phenomena occurring on time scales inaccessible to standard molecular dynamics, but evaluation of the physical time scales involved in the processes is not directly possible from them. For this reason, the problem of recovering dynamics from such kinds of simulations is the object of very active research due to the relevant theoretical and practical implications of dynamics on the properties of both natural and synthetic molecular systems. In a recent paper [A. Rapallo et al., J. Comput. Chem. 42, 586-599 (2021)], it has been shown how the coupling of BAMS (which destroys the dynamics but allows to calculate average properties) with Extended Diffusion Theory (EDT) (which requires input appropriate equilibrium averages calculated over the BAMS trajectories) allows to effectively use the Smoluchowski equation to calculate the orientational time correlation function of the head-tail unit vector defined over a peptide in water solution. Orientational relaxation of this vector is the result of the coupling of internal molecular motions with overall molecular rotation, and it was very well described by correlation functions expressed in terms of weighted sums of suitable time-exponentially decaying functions, in agreement with a Brownian diffusive regime. However, situations occur where exponentially decaying functions are no longer appropriate to capture the actual dynamical behavior, which exhibits persistent long time correlations, compatible with the so called subdiffusive regimes. In this paper, a generalization of EDT will be given, exploiting a fractional Smoluchowski equation (FEDT) to capture the non-exponential character observed in the relaxation of intramolecular distances and molecular radius of gyration, whose dynamics depend on internal molecular motions only. The calculation methods, proper to EDT, are adapted to implement the generalization of the theory, and the resulting algorithm confirms FEDT as a tool of practical value in recovering dynamics from BAMS, to be used in general situations, involving both regular and anomalous diffusion regimes.

Elucidating molecular mechanisms of protoxin-II state-specific binding to the human NaV1.7 channel

Human voltage-gated sodium (hNaV) channels are responsible for initiating and propagating action potentials in excitable cells, and mutations have been associated with numerous cardiac and neurological disorders. hNaV1.7 channels are expressed in peripheral neurons and are promising targets for pain therapy. The tarantula venom peptide protoxin-II (PTx2) has high selectivity for hNaV1.7 and is a valuable scaffold for designing novel therapeutics to treat pain. Here, we used computational modeling to study the molecular mechanisms of the state-dependent binding of PTx2 to hNaV1.7 voltage-sensing domains (VSDs). Using Rosetta structural modeling methods, we constructed atomistic models of the hNaV1.7 VSD II and IV in the activated and deactivated states with docked PTx2. We then performed microsecond-long all-atom molecular dynamics (MD) simulations of the systems in hydrated lipid bilayers. Our simulations revealed that PTx2 binds most favorably to the deactivated VSD II and activated VSD IV. These state-specific interactions are mediated primarily by PTx2's residues R22, K26, K27, K28, and W30 with VSD and the surrounding membrane lipids. Our work revealed important protein-protein and protein-lipid contacts that contribute to high-affinity state-dependent toxin interaction with the channel. The workflow presented will prove useful for designing novel peptides with improved selectivity and potency for more effective and safe treatment of pain.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

Elucidating Molecular Mechanisms of Protoxin-2 State-specific Binding to the Human NaV1.7 Channel

Human voltage-gated sodium (hNaV) channels are responsible for initiating and propagating action potentials in excitable cells and mutations have been associated with numerous cardiac and neurological disorders. hNaV1.7 channels are expressed in peripheral neurons and are promising targets for pain therapy. The tarantula venom peptide protoxin-2 (PTx2) has high selectivity for hNaV1.7 and serves as a valuable scaffold to design novel therapeutics to treat pain. Here, we used computational modeling to study the molecular mechanisms of the state-dependent binding of PTx2 to hNaV1.7 voltage-sensing domains (VSDs). Using Rosetta structural modeling methods, we constructed atomistic models of the hNaV1.7 VSD II and IV in the activated and deactivated states with docked PTx2. We then performed microsecond-long all-atom molecular dynamics (MD) simulations of the systems in hydrated lipid bilayers. Our simulations revealed that PTx2 binds most favorably to the deactivated VSD II and activated VSD IV. These state-specific interactions are mediated primarily by PTx2's residues R22, K26, K27, K28, and W30 with VSD as well as the surrounding membrane lipids. Our work revealed important protein-protein and protein-lipid contacts that contribute to high-affinity state-dependent toxin interaction with the channel. The workflow presented will prove useful for designing novel peptides with improved selectivity and potency for more effective and safe treatment of pain.

Selectivity and Ranking of Tight-Binding JAK-STAT Inhibitors Using Markovian Milestoning with Voronoi Tessellations

Janus kinases (JAK), a group of proteins in the nonreceptor tyrosine kinase (NRTKs) family, play a crucial role in growth, survival, and angiogenesis. They are activated by cytokines through the Janus kinase-signal transducer and activator of a transcription (JAK-STAT) signaling pathway. JAK-STAT signaling pathways have significant roles in the regulation of cell division, apoptosis, and immunity. Identification of the V617F mutation in the Janus homology 2 (JH2) domain of JAK2 leading to myeloproliferative disorders has stimulated great interest in the drug discovery community to develop JAK2-specific inhibitors. However, such inhibitors should be selective toward JAK2 over other JAKs and display an extended residence time. Recently, novel JAK2/STAT5 axis inhibitors (N-(1H-pyrazol-3-yl)pyrimidin-2-amino derivatives) have displayed extended residence times (hours or longer) on target and adequate selectivity excluding JAK3. To facilitate a deeper understanding of the kinase-inhibitor interactions and advance the development of such inhibitors, we utilize a multiscale Markovian milestoning with Voronoi tessellations (MMVT) approach within the Simulation-Enabled Estimation of Kinetic Rates v.2 (SEEKR2) program to rank order these inhibitors based on their kinetic properties and further explain the selectivity of JAK2 inhibitors over JAK3. Our approach investigates the kinetic and thermodynamic properties of JAK-inhibitor complexes in a user-friendly, fast, efficient, and accurate manner compared to other brute force and hybrid-enhanced sampling approaches.

DeepWEST: Deep Learning of Kinetic Models with the Weighted Ensemble Simulation Toolkit for Enhanced Sampling

Recent advances in computational power and algorithms have enabled molecular dynamics (MD) simulations to reach greater time scales. However, for observing conformational transitions associated with biomolecular processes, MD simulations still have limitations. Several enhanced sampling techniques seek to address this challenge, including the weighted ensemble (WE) method, which samples transitions between metastable states using many weighted trajectories to estimate kinetic rate constants. However, initial sampling of the potential energy surface has a significant impact on the performance of WE, i.e., convergence and efficiency. We therefore introduce deep-learned kinetic modeling approaches that extract statistically relevant information from short MD trajectories to provide a well-sampled initial state distribution for WE simulations. This hybrid approach overcomes any statistical bias to the system, as it runs short unbiased MD trajectories and identifies meaningful metastable states of the system. It is shown to provide a more refined free energy landscape closer to the steady state that could efficiently sample kinetic properties such as rate constants.

Weighted ensemble: Recent mathematical developments

Weighted ensemble (WE) is an enhanced sampling method based on periodically replicating and pruning trajectories generated in parallel. WE has grown increasingly popular for computational biochemistry problems due, in part, to improved hardware and accessible software implementations. Algorithmic and analytical improvements have played an important role, and progress has accelerated in recent years. Here, we discuss and elaborate on the WE method from a mathematical perspective, highlighting recent results that enhance the computational efficiency. The mathematical theory reveals a new strategy for optimizing trajectory management that approaches the best possible variance while generalizing to systems of arbitrary dimension.

Engineering and Design of Programmable Genome Editors

Programmable genome editors are enzymes that can be targeted to a specific location in the genome for making site-specific alterations or deletions. The engineering, design, and development of sequence-specific editors has resulted in a dramatic increase in the precision of editing for nucleotide sequences. These editors can target specific locations in a genome, in vivo. The genome editors are being deployed for the development of genetically modified organisms for agriculture and industry, and for gene therapy of inherited human genetic disorders, cancer, and immunotherapy. Experimental and computational studies of structure, binding, activity, dynamics, and folding, reviewed here, have provided valuable insights that have the potential for increasing the functional efficiency of these gene/genome editors. Biochemical and biophysical studies of the specificities of natural and engineered genome editors reveal that increased binding affinity can be detrimental because of the increase of off-target effects and that the engineering and design of genome editors with higher specificity may require modulation and control of the conformational dynamics.

Protein-Protein Interaction-Gaussian Accelerated Molecular Dynamics (PPI-GaMD): Characterization of Protein Binding Thermodynamics and Kinetics

Protein-protein interactions (PPIs) play key roles in many fundamental biological processes such as cellular signaling and immune responses. However, it has proven challenging to simulate repetitive protein association and dissociation in order to calculate binding free energies and kinetics of PPIs due to long biological timescales and complex protein dynamics. To address this challenge, we have developed a new computational approach to all-atom simulations of PPIs based on a robust Gaussian accelerated molecular dynamics (GaMD) technique. The method, termed "PPI-GaMD", selectively boosts interaction potential energy between protein partners to facilitate their slow dissociation. Meanwhile, another boost potential is applied to the remaining potential energy of the entire system to effectively model the protein's flexibility and rebinding. PPI-GaMD has been demonstrated on a model system of the ribonuclease barnase interactions with its inhibitor barstar. Six independent 2 μs PPI-GaMD simulations have captured repetitive barstar dissociation and rebinding events, which enable calculations of the protein binding thermodynamics and kinetics simultaneously. The calculated binding free energies and kinetic rate constants agree well with the experimental data. Furthermore, PPI-GaMD simulations have provided mechanistic insights into barstar binding to barnase, which involves long-range electrostatic interactions and multiple binding pathways, being consistent with previous experimental and computational findings of this model system. In summary, PPI-GaMD provides a highly efficient and easy-to-use approach for binding free energy and kinetics calculations of PPIs.