Dihedral angle principal component analysis of molecular dynamics simulations

Repurposing of Drug Bank Compounds against Plasmodium falciparum Dihydroorotate Dehydrogenase as novel anti malarial drug candidates by Computational approaches

This study aimed to repurpose Drug Bank Compounds against P. falciparum Dihydroorotate dehydrogenase (Pf-DHODH)a potential molecular target for antimalarial drug development due to its vital role in P. falciparum survival. Initially, the MATGEN server was used to screen drugs against Pf-DHODH (PDB ID 6GJG), followed by revalidating the results through docking by Autodock Vina through PyRx. Based on the docking results, three drugs namely, Talnifumate, Sulfaphenazole, and (3S)-N-[(2S)-1-[2-(1H-indol-3-yl)ethylamino]-1-oxopropan-2-yl]-1-(4-methoxyphenyl)-5-oxopyrrolidine-3-carboxamide-were subjected to molecular dynamics simulation for 100 ns. Molecular dynamics simulation results indicate that (3S)-N-[(2S)-1-[2-(1H-indol-3-yl)ethylamino]-1-oxopropan-2-yl]-1-(4-methoxyphenyl)-5-oxopyrrolidine-3-carboxamide- and Sulfaphenazole may target Pf-DHODH by forming a stable protein-ligand complex as they showed better free binding energy -130.58 kJ/mol, and -79.84 kJ/mol, respectively as compared to the free binding energy 116.255 kJ/mol of the reference compound; 3,6-dimethyl- ~ {N}-[4-(trifluoromethyl)phenyl]-[1,2]oxazolo[5,4-d]pyrimidin-4-amine. Although the studied compounds are drugs, still we applied Lipinski's rules and ADMET analysis that reconfirmed that these drugs have favorable drug-like properties. In conclusion, the results of the study show that Talniflumate and Sulfaphenazole may be potential antimalarial drug candidates.The derivatives of these drugs could be designed and tested to develop better drugs against Plasmodium species.

Graphical Abstract

Sampling Conformational Ensembles of Highly Dynamic Proteins via Generative Deep Learning

Proteins are inherently dynamic, and their conformational ensembles are functionally important in biology. Large-scale motions may govern protein structure-function relationship, and numerous transient but stable conformations of intrinsically disordered proteins (IDPs) can play a crucial role in biological function. Investigating conformational ensembles to understand regulations and disease-related aggregations of IDPs is challenging both experimentally and computationally. In this paper we first introduced an unsupervised deep learning-based model, termed Internal Coordinate Net (ICoN), which learns the physical principles of conformational changes from molecular dynamics (MD) simulation data. Second, we selected interpolating data points in the learned latent space that rapidly identify novel synthetic conformations with sophisticated and large-scale sidechains and backbone arrangements. Third, with the highly dynamic amyloid-β 1-42 (Aβ42) monomer, our deep learning model provided a comprehensive sampling of Aβ42's conformational landscape. Analysis of these synthetic conformations revealed conformational clusters that can be used to rationalize experimental findings. Additionally, the method can identify novel conformations with important interactions in atomistic details that are not included in the training data. New synthetic conformations showed distinct sidechain rearrangements that are probed by our EPR and amino acid substitution studies. This approach is highly transferable and can be used for any available data for training. The work also demonstrated the ability for deep learning to utilize learned natural atomistic motions in protein conformation sampling.

Elucidating the structural basis of vitamin B12 derivatives as novel potent inhibitors of PTP1B: Insights from inhibitory mechanisms using Gaussian accelerated molecular dynamics (GaMD) and in vitro study

Vitamin B12 is a group of biologically active cobalamin compounds. In this study, we investigated the inhibitory effects of methylcobalamin (MeCbl) and hydroxocobalamin acetate (OHCbl Acetate) on protein tyrosine phosphatase 1B (PTP1B). MeCbl and OHCbl Acetate exhibited an IC50 of approximately 58.390 ± 2.811 μM and 8.998 ± 0.587 μM, respectively. The Ki values of MeCbl and OHCbl Acetate were 25.01 μM and 4.04 μM respectively. To elucidate the inhibition mechanism, we conducted a 500 ns Gaussian accelerated molecular dynamics (GaMD) simulation. Utilizing PCA and tICA, we constructed Markov state models (MSM) to examine secondary structure changes during motion. Our findings revealed that the α-helix at residues 37-42 remained the most stable in the PTP1B-OHCbl Acetate system. Furthermore, upon binding of OHCbl Acetate or MeCbl, the WPD loop of PTP1B moved inward to the active pocket, forming a closed conformation and potentially obstructs substrate entry. Protein-ligand interaction analysis and MM-PBSA showed that OHCbl Acetate exhibited lower binding free energy and engaged in more residue interactions with PTP1B. In summary, our study confirmed the substantial inhibitory activity of OHCbl Acetate against PTP1B, with its inhibitory potency notably surpassing that of MeCbl. We demonstrated potential molecular mechanisms of OHCbl Acetate inhibiting PTP1B.

Functional protein dynamics in a crystal

Proteins are molecular machines and to understand how they work, we need to understand how they move. New pump-probe time-resolved X-ray diffraction methods open up ways to initiate and observe protein motions with atomistic detail in crystals on biologically relevant timescales. However, practical limitations of these experiments demands parallel development of effective molecular dynamics approaches to accelerate progress and extract meaning. Here, we establish robust and accurate methods for simulating dynamics in protein crystals, a nontrivial process requiring careful attention to equilibration, environmental composition, and choice of force fields. With more than seven milliseconds of sampling of a single chain, we identify critical factors controlling agreement between simulation and experiments and show that simulated motions recapitulate ligand-induced conformational changes. This work enables a virtuous cycle between simulation and experiments for visualizing and understanding the basic functional motions of proteins.

Interpretation of autoencoder-learned collective variables using Morse-Smale complex and sublevelset persistent homology: An application on molecular trajectories

Dimensionality reduction often serves as the first step toward a minimalist understanding of physical systems as well as the accelerated simulations of them. In particular, neural network-based nonlinear dimensionality reduction methods, such as autoencoders, have shown promising outcomes in uncovering collective variables (CVs). However, the physical meaning of these CVs remains largely elusive. In this work, we constructed a framework that (1) determines the optimal number of CVs needed to capture the essential molecular motions using an ensemble of hierarchical autoencoders and (2) provides topology-based interpretations to the autoencoder-learned CVs with Morse-Smale complex and sublevelset persistent homology. This approach was exemplified using a series of n-alkanes and can be regarded as a general, explainable nonlinear dimensionality reduction method.

Characterization of Plant-Derived Natural Inhibitors of Dipeptidyl Peptidase-4 as Potential Antidiabetic Agents: A Computational Study

Diabetes, characterized by elevated blood sugar levels, poses significant health and economic risks, correlating with complications like cardiovascular disease, kidney failure, and blindness. Dipeptidyl peptidase-4 (DPP-4), also referred to as T-cell activation antigen CD26 (EC 3.4.14.5.), plays a crucial role in glucose metabolism and immune function. Inhibiting DPP-4 was anticipated as a potential new therapy for diabetes. Therefore, identification of plant-based natural inhibitors of DPP-4 would help in eradicating diabetes worldwide. Here, for the identification of the potential natural inhibitors of DPP-4, we developed a phytochemicals library consisting of over 6000 phytochemicals detected in 81 medicinal plants that exhibited anti-diabetic potency. The library has been docked against the target proteins, where isorhamnetin, Benzyl 5-Amino-5-deoxy-2,3-O-isopropyl-alpha-D-mannofuranoside (DTXSID90724586), and 5-Oxo-7-[4-(trifluoromethyl) phenyl]-4H,6H,7H-[1,2]thiazolo[4,5-b]pyridine 3-carboxylic acid (CHEMBL3446108) showed binding affinities of -8.5, -8.3, and -8.3 kcal/mol, respectively. These compounds exhibiting strong interactions with DPP-4 active sites (Glu205, Glu206, Tyr547, Trp629, Ser630, Tyr662, His740) were identified. ADME/T and bioactivity predictions affirmed their pharmacological safety. Density functional theory calculations assessed stability and reactivity, while molecular dynamics simulations demonstrated persistent stability. Analyzing parameters like RMSD, RG, RMSF, SASA, H-bonds, MM-PBSA, and FEL confirmed stable protein-ligand compound formation. Principal component analysis provided structural variation insights. Our findings suggest that those compounds might be possible candidates for developing novel inhibitors targeting DPP-4 for treating diabetes.

Structure-guided engineering and molecular simulations to design a potent monoclonal antibody to target aP2 antigen for adaptive immune response instigation against type 2 diabetes

Introduction

Diabetes mellitus (DM) is recognized as one of the oldest chronic diseases and has become a significant public health issue, necessitating innovative therapeutic strategies to enhance patient outcomes. Traditional treatments have provided limited success, highlighting the need for novel approaches in managing this complex disease.

Methods

In our study, we employed graph signature-based methodologies in conjunction with molecular simulation and free energy calculations. The objective was to engineer the CA33 monoclonal antibody for effective targeting of the aP2 antigen, aiming to elicit a potent immune response. This approach involved screening a mutational landscape comprising 57 mutants to identify modifications that yield significant enhancements in binding efficacy and stability.

Results

Analysis of the mutational landscape revealed that only five substitutions resulted in noteworthy improvements. Among these, mutations T94M, A96E, A96Q, and T94W were identified through molecular docking experiments to exhibit higher docking scores compared to the wild-type. Further validation was provided by calculating the dissociation constant (KD), which showed a similar trend in favor of these mutations. Molecular simulation analyses highlighted T94M as the most stable complex, with reduced internal fluctuations upon binding. Principal components analysis (PCA) indicated that both the wild-type and T94M mutant displayed similar patterns of constrained and restricted motion across principal components. The free energy landscape analysis underscored a single metastable state for all complexes, indicating limited structural variability and potential for high therapeutic efficacy against aP2. Total binding free energy (TBE) calculations further supported the superior performance of the T94M mutation, with TBE values demonstrating the enhanced binding affinity of selected mutants over the wild-type.

Discussion

Our findings suggest that the T94M substitution, along with other identified mutations, significantly enhances the therapeutic potential of the CA33 antibody against DM by improving its binding affinity and stability. These results not only contribute to a deeper understanding of antibody-antigen interactions in the context of DM but also provide a valuable framework for the rational design of antibodies aimed at targeting this disease more effectively.

Dramatic Differences between the Structural Susceptibility of the S1 Pre- and S2 Postfusion States of the SARS-CoV-2 Spike Protein to External Electric Fields Revealed by Molecular Dynamics Simulations

In its prefusion state, the SARS-CoV-2 spike protein (similarly to other class I viral fusion proteins) is metastable, which is considered to be an important feature for optimizing or regulating its functions. After the binding process of its S1 subunit (S1) with ACE2, the spike protein (S) undergoes a dramatic conformational change where S1 splits from the S2 subunit, which then penetrates the membrane of the host cell, promoting the fusion of the viral and cell membranes. This results in the infection of the host cell. In a previous work, we showed-using large-scale molecular dynamics simulations-that the application of external electric fields (EFs) induces drastic changes and damage in the receptor-binding domain (RBD) of the wild-type spike protein, as well of the Alpha, Beta, and Gamma variants, leaving a structure which cannot be recognized anymore by ACE2. In this work, we first extend the study to the Delta and Omicron variants and confirm the high sensitivity and extreme vulnerability of the RBD of the prefusion state of S to moderate EF (as weak as 104 V/m), but, more importantly, we also show that, in contrast, the S2 subunit of the postfusion state of the spike protein does not suffer structural damage even if electric field intensities four orders of magnitude higher are applied. These results provide a solid scientific basis to confirm the connection between the prefusion-state metastability of the SARS-CoV-2 spike protein and its susceptibility to be damaged by EF. After the virus docks to the ACE2 receptor, the stable and robust postfusion conformation develops, which exhibits a similar resistance to EF (damage threshold higher than 108 V/m) like most globular proteins.

Probing the mechanisms of hydrazide-based HDAC inhibitors binding to HDAC3 using Gaussian accelerated molecular dynamics (GaMD) simulations

Histone deacetylases (HDACs) have emerged as promising targets for anticancer drug development. They regulate gene expression by removing acetyl groups from lysine residues on histone tails, leading to chromatin condensation. A hydrazide-based HDAC inhibitor, N-(4-(2-Propylhydrazine-1-carbonyl)benzyl)-1H-indole-2-carboxamide (11h), has been reported to exhibit significant in vivo antitumor activity. In comparison to the lead compound N-(4-(2-Propylhydrazine-1-carbonyl)benzyl)cinnamamide (17), compound 11h demonstrates 2- to 5-fold higher HDAC inhibition and cell-based antitumor activity. However, the inhibitory mechanism of 11h remains insufficiently explored. In this study, we conducted 500 ns Gaussian Accelerated Molecular Dynamics (GaMD) simulations on Histone deacetylase 3 (HDAC3) and two complex systems (HDAC3-17 and HDAC3-11h). Our findings revealed that upon inhibitor binding, the active pocket volume of HDAC3 undergone alterations, and the movement of the L6-loop toward the active site impeded substrate entry. Moreover, we observed a destabilization of the α-helix in the aa75-89 region of HDAC3 compared to the two complex systems, indicating partial unwinding. Notably, 11h exhibited a closer proximity of its carbonyl oxygen to the active pocket's Zn2+ metal compared to 17, increasing the likelihood of coordination with the Zn2+ metal. The analysis of protein-ligand interactions highlighted a greater number of hydrogen bonds and other interactions between 11h and the receptor protein when compared to 17, underscoring the stronger binding of 11h to HDAC3. In conclusion, our study provided theoretical insights into the inhibitory mechanism of hydrazide-based HDAC inhibitors on HDAC3, thereby contributing to the development of improved drug targets for cancer therapy.Communicated by Ramaswamy H. Sarma.

Toward a Benchmark for Markov State Models: The Folding of HP35

Adopting a 300 μs long MD trajectory of the folding of villin headpiece (HP35) by D. E. Shaw Research, we recently constructed a Markov state model (MSM) based on inter-residue contacts. The model reproduces the folding time and predicts that the native basin and unfolded region consist of metastable substates that are structurally well-characterized. Recognizing the need to establish well-defined benchmark problems, we study to what extent and in what sense this MSM can be employed as a reference model. Hence, we test the robustness of the MSM by comparing it to models that use alternative combinations of features, dimensionality reduction methods, and clustering schemes. The study suggests some main characteristics of the folding of HP35 that should be reproduced by other competitive models. Moreover, the discussion reveals which parts of the MSM workflow matter most for the considered problem and illustrates the promises and pitfalls of state-based models for the interpretation of biomolecular simulations.

SARS-CoV-2 Spike Protein Interaction Space

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a +sense single-strand RNA virus. The virus has four major surface proteins: spike (S), envelope (E), membrane (M), and nucleocapsid (N), respectively. The constitutive proteins present a high grade of symmetry. Identifying a binding site is difficult. The virion is approximately 50-200 nm in diameter. Angiotensin-converting enzyme 2 (ACE2) acts as the cell receptor for the virus. SARS-CoV-2 has an increased affinity to human ACE2 compared with the original SAR strain. Topological space, and its symmetry, is a critical component in molecular interactions. By exploring this space, a suitable ligand space can be characterized accordingly. A spike protein (S) computational model in a complex with ACE 2 was generated using silica methods. Topological spaces were probed using high computational throughput screening techniques to identify and characterize the topological space of both SARS and SARS-CoV-2 spike protein and its ligand space. In order to identify the symmetry clusters, computational analysis techniques, together with statistical analysis, were utilized. The computations are based on crystallographic protein data bank PDB-based models of constitutive proteins. Cartesian coordinates of component atoms and some cluster maps were generated and analyzed. Dihedral angles were used in order to compute a topological receptor space. This computational study uses a multimodal representation of spike protein interactions with some fragment proteins. The chemical space of the receptors (a dimensional volume) suggests the relevance of the receptor as a drug target. The spike protein S of SARS and SARS-CoV-2 is analyzed and compared. The results suggest a mirror symmetry of SARS and SARS-CoV-2 spike proteins. The results show thatSARS-CoV-2 space is variable and has a distinct topology. In conclusion, surface proteins grant virion variability and symmetry in interactions with a potential complementary target (protein, antibody, ligand). The mirror symmetry of dihedral angle clusters determines a high specificity of the receptor space.

MDiGest: A Python package for describing allostery from molecular dynamics simulations

Many biological processes are regulated by allosteric mechanisms that communicate with distant sites in the protein responsible for functionality. The binding of a small molecule at an allosteric site typically induces conformational changes that propagate through the protein along allosteric pathways regulating enzymatic activity. Elucidating those communication pathways from allosteric sites to orthosteric sites is, therefore, essential to gain insights into biochemical processes. Targeting the allosteric pathways by mutagenesis can allow the engineering of proteins with desired functions. Furthermore, binding small molecule modulators along the allosteric pathways is a viable approach to target reactions using allosteric inhibitors/activators with temporal and spatial selectivity. Methods based on network theory can elucidate protein communication networks through the analysis of pairwise correlations observed in molecular dynamics (MD) simulations using molecular descriptors that serve as proxies for allosteric information. Typically, single atomic descriptors such as α-carbon displacements are used as proxies for allosteric information. Therefore, allosteric networks are based on correlations revealed by that descriptor. Here, we introduce a Python software package that provides a comprehensive toolkit for studying allostery from MD simulations of biochemical systems. MDiGest offers the ability to describe protein dynamics by combining different approaches, such as correlations of atomic displacements or dihedral angles, as well as a novel approach based on the correlation of Kabsch-Sander electrostatic couplings. MDiGest allows for comparisons of networks and community structures that capture physical information relevant to allostery. Multiple complementary tools for studying essential dynamics include principal component analysis, root mean square fluctuation, as well as secondary structure-based analyses.

A deep encoder-decoder framework for identifying distinct ligand binding pathways

The pathway(s) that a ligand would adopt en route to its trajectory to the native pocket of the receptor protein act as a key determinant of its biological activity. While Molecular Dynamics (MD) simulations have emerged as the method of choice for modeling protein-ligand binding events, the high dimensional nature of the MD-derived trajectories often remains a barrier in the statistical elucidation of distinct ligand binding pathways due to the stochasticity inherent in the ligand's fluctuation in the solution and around the receptor. Here, we demonstrate that an autoencoder based deep neural network, trained using an objective input feature of a large matrix of residue-ligand distances, can efficiently produce an optimal low-dimensional latent space that stores necessary information on the ligand-binding event. In particular, for a system of L99A mutant of T4 lysozyme interacting with its native ligand, benzene, this deep encoder-decoder framework automatically identifies multiple distinct recognition pathways, without requiring user intervention. The intermediates involve the spatially discrete location of the ligand in different helices of the protein before its eventual recognition of native pose. The compressed subspace derived from the autoencoder provides a quantitatively accurate measure of the free energy and kinetics of ligand binding to the native pocket. The investigation also recommends that while a linear dimensional reduction technique, such as time-structured independent component analysis, can do a decent job of state-space decomposition in cases where the intermediates are long-lived, autoencoder is the method of choice in systems where transient, low-populated intermediates can lead to multiple ligand-binding pathways.

Relationships between Molecular Structural Order Parameters and Equilibrium Water Dynamics in Aqueous Mixtures

Water's unique thermophysical properties and how it mediates aqueous interactions between solutes have long been interpreted in terms of its collective molecular structure. The seminal work of Errington and Debenedetti [Nature 2001, 409, 318-321] revealed a striking hierarchy of relationships among the thermodynamic, dynamic, and structural properties of water, motivating many efforts to understand (1) what measures of water structure are connected to different experimentally accessible macroscopic responses and (2) how many such structural metrics are adequate to describe the collective structural behavior of water. Diffusivity constitutes a particularly interesting experimentally accessible equilibrium property to investigate such relationships because advanced NMR techniques allow the measurement of bulk and local water dynamics in nanometer proximity to molecules and interfaces, suggesting the enticing possibility of measuring local diffusivities that report on water structure. Here, we apply statistical learning methods to discover persistent structure-dynamic correlations across a variety of simulated aqueous mixtures, from alcohol-water to polypeptoid-water systems. We investigate a variety of molecular water structure metrics and find that an unsupervised statistical learning algorithm (namely, sequential feature selection) identifies only two or three independent structural metrics that are sufficient to predict water self-diffusivity accurately. Surprisingly, the translational diffusivity of water across all mixed systems studied here is strongly correlated with a measure of tetrahedral order given by water's triplet angle distribution. We also identify a separate small number of structural metrics that well predict an important thermodynamic property, the excess chemical potential of an idealized methane-sized hydrophobe in water. Ultimately, we offer a Bayesian method of inferring water structure by using only structure-dynamics linear regression models with experimental Overhauser dynamic nuclear polarization (ODNP) measurements of water self-diffusivity. This study thus quantifies the relationships among several distinct structural order parameters in water and, through statistical learning, reveals the potential to leverage molecular structure to predict fundamental thermophysical properties. In turn, these findings suggest a framework for solving the inverse problem of inferring water's molecular structure using experimental measurements such as ODNP studies that probe local water properties.

Selecting Features for Markov Modeling: A Case Study on HP35

Markov state models represent a popular means to interpret molecular dynamics trajectories in terms of memoryless transitions between metastable conformational states. To provide a mechanistic understanding of the considered biomolecular process, these states should reflect structurally distinct conformations and ensure a time scale separation between fast intrastate and slow interstate dynamics. Adopting the folding of villin headpiece (HP35) as a well-established model problem, here we discuss the selection of suitable input coordinates or "features", such as backbone dihedral angles and interresidue distances. We show that dihedral angles account accurately for the structure of the native energy basin of HP35, while the unfolded region of the free energy landscape and the folding process are best described by tertiary contacts of the protein. To construct a contact-based model, we consider various ways to define and select contact distances and introduce a low-pass filtering of the feature trajectory as well as a correlation-based characterization of states. Relying on input data that faithfully account for the mechanistic origin of the studied process, the states of the resulting Markov model are clearly discriminated by the features, describe consistently the hierarchical structure of the free energy landscape, and─as a consequence─correctly reproduce the slow time scales of the process.

A Retrospective on the Development of Methods for the Analysis of Protein Conformational Ensembles

Analysing protein conformational ensembles whether from molecular dynamics (MD) simulation or other sources for functionally relevant conformational changes can be very challenging. In the nineteen nineties dimensional reduction methods were developed primarily for analysing MD trajectories to determine dominant motions with the aim of understanding their relationship to function. Coarse-graining methods were also developed so the conformational change between two structures could be described in terms of the relative motion of a small number of quasi-rigid regions rather than in terms of a large number of atoms. When these methods are combined, they can characterize the large-scale motions inherent in a conformational ensemble providing insight into possible functional mechanism. The dimensional reduction methods first applied to protein conformational ensembles were referred to as Quasi-Harmonic Analysis, Principal Component Analysis and Essential Dynamics Analysis. A retrospective on the origin of these methods is presented, the relationships between them explained, and more recent developments reviewed.

Unveiling the Inhibitory Potentials of Peptidomimetic Azanitriles and Pyridyl Esters towards SARS-CoV-2 Main Protease: A Molecular Modelling Investigation

The severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) is responsible for COVID-19, which was declared a global pandemic in March 2020 by the World Health Organization (WHO). Since SARS-CoV-2 main protease plays an essential role in the virus’s life cycle, the design of small drug molecules with lower molecular weight has been a promising development targeting its inhibition. Herein, we evaluated the novel peptidomimetic azatripeptide and azatetrapeptide nitriles against SARS-CoV-2 main protease. We employed molecular dynamics (MD) simulations to elucidate the selected compounds’ binding free energy profiles against SARS-CoV-2 and further unveil the residues responsible for the drug-binding properties. Compound 8 exhibited the highest binding free energy of −49.37 ± 0.15 kcal/mol, followed by compound 7 (−39.83 ± 0.19 kcal/mol), while compound 17 showed the lowest binding free energy (−23.54 ± 0.19 kcal/mol). In addition, the absorption, distribution, metabolism, and excretion (ADME) assessment was performed and revealed that only compound 17 met the drug-likeness parameters and exhibited high pharmacokinetics to inhibit CYP1A2, CYP2C19, and CYP2C9 with better absorption potential and blood-brain barrier permeability (BBB) index. The additional intermolecular evaluations suggested compound 8 as a promising drug candidate for inhibiting SARS-CoV-2 Mpro. The substitution of isopropane in compound 7 with an aromatic benzene ring in compound 8 significantly enhanced the drug’s ability to bind better at the active site of the SARS-CoV-2 Mpro.

DFT and TD-DFT study of hydrogen bonded complexes of aspartic acid and n water (n = 1 and 2)

CONTEXT

Hydrogen bonds (HB) influence the conformational preferences of biomolecules and their optical and electronic properties. The directional interaction of molecules of water can be a prototype to understand the effects of HBs on biomolecules. Among the neurotransmitters (NT), L-aspartic acid (ASP) stands out due to its importance in health and as a precursor of several biomolecules. As it presents different functional groups and readily forms inter- and intramolecular HBs, ASP can be considered a prototype for understanding the behavior of NTs when interacting by HB with other substances. Although several theoretical studies have been performed in the past on isolated ASP and its formed complexes with water, both in gas and liquid phases, using DFT and TD-DFT formalisms, these works did not perform large basis set calculations or study electronic transitions of ASP-water complexes. We investigated the HB interactions in complexes of ASP and water molecules. The results show that the interactions between the carboxylic groups of ASP with water molecules, forming cyclic structures with two HBs, lead to more stable and less polar complexes than other conformers formed between water and the NH₂ group. It was observed that there is a relationship between the deviation in the UV-Vis absorption band of the ASP and the interactions of water with the HOMO and LUMO orbitals with the stabilization/destabilization of the S₁ state to the S₀ of the complexes. However, in some cases, such as 1:1 complex ASP-W2, this analysis may be inaccurate due to small changes in ΔE.

METHODS

We studied the landscapes of the ground state surface of different conformers of isolated L-ASP and the L-ASP-(H₂O)_n complexes (n = 1 and 2) using the DFT formalism, with the B3LYP functional, and six different basis sets: 6-31 +  + G(d,p), 6-311 +  + G(d,p), D95 +  + (d,p), D95V +  + (d,p), cc-pVDZ, and, cc-pVTZ basis sets. The cc-pVTZ basis set provides the minimum energy of all conformers, and therefore, we performed the analysis with this basis set. We evaluated the stabilization of the ASP and complexes using the minimum ground state energy, corrected by the zero point energy and the interaction energy between the ASP and the water molecules. We also calculated the vertical electronic transitions S₁ ← S₀, and their properties using the TD-DFT formalism at B3LYP/cc-pVTZ level with the optimized geometries for S₀ state with the same basis set. For the analysis of the vertical transitions of isolated ASP and the ASP-(H₂O)_n complexes, we calculated the electrostatic energy in the S₀ and S₁ states. We performed the calculations with the Gaussian 09 software package. We used the VMD software package to visualize the geometries and shapes of the molecule and complexes.

Revealing Free Energy Landscape From MD Data via Conditional Angle Partition Tree

Deciphering the free energy landscape of biomolecular structure space is crucial for understanding many complex molecular processes, such as protein-protein interaction, RNA folding, and protein folding. A major source of current dynamic structure data is Molecular Dynamics (MD) simulations. Several methods have been proposed to investigate the free energy landscape from MD data, but all of them rely on the assumption that kinetic similarity is associated with global geometric similarity, which may lead to unsatisfactory results. In this paper, we proposed a new method called Conditional Angle Partition Tree to reveal the hierarchical free energy landscape by correlating local geometric similarity with kinetic similarity. Its application on the benchmark alanine dipeptide MD data showed a much better performance than existing methods in exploring and understanding the free energy landscape. We also applied it to the MD data of Villin HP35. Our results are more reasonable on various aspects than those from other methods and very informative on the hierarchical structure of its energy landscape.

Conformational ensemble of the NSP1 CTD in SARS-CoV-2: Perspectives from the free energy landscape

The nonstructural protein-1 (NSP1) of the severe acute respiratory syndrome-associated coronavirus 2 plays a crucial role in the translational shutdown and immune evasion inside host cells. Despite its known intrinsic disorder, the C-terminal domain (CTD) of NSP1 has been reported to form a double α-helical structure and block the 40S-ribosomal channel for mRNA translation. Experimental studies indicate that NSP1 CTD functions independently from the globular N-terminal region separated with a long linker domain, underscoring the necessity of exploring the standalone conformational ensemble. In this contribution, we utilize exascale computing resources to yield unbiased molecular dynamics simulation of NSP1 CTD in all-atom resolution starting from multiple initial seed structures. A data-driven approach elicits collective variables (CVs) that are significantly superior to conventional descriptors in capturing the conformational heterogeneity. The free energy landscape as a function of the CV space is estimated using the modified expectation maximized molecular dynamics. Originally developed by us for small peptides, here, we establish the efficacy of expectation maximized molecular dynamics in conjunction with data-driven CV space for a more complex and relevant biomolecular system. The results reveal the existence of two disordered metastable populations in the free energy landscape that are separated from the conformation resembling ribosomal subunit bound state by high kinetic barriers. Chemical shift correlation and secondary structure analysis capture significant differences among key structures of the ensemble. Altogether, these insights can underpin drug development studies and mutational experiments that help induce population shifts to alter translational blocking and understand its molecular basis in further detail.

A new approach for extracting information from protein dynamics

Increased ability to predict protein structures is moving research focus towards understanding protein dynamics. A promising approach is to represent protein dynamics through networks and take advantage of well-developed methods from network science. Most studies build protein dynamics networks from correlation measures, an approach that only works under very specific conditions, instead of the more robust inverse approach. Thus, we apply the inverse approach to the dynamics of protein dihedral angles, a system of internal coordinates, to avoid structural alignment. Using the well-characterized adhesion protein, FimH, we show that our method identifies networks that are physically interpretable, robust, and relevant to the allosteric pathway sites. We further use our approach to detect dynamical differences, despite structural similarity, for Siglec-8 in the immune system, and the SARS-CoV-2 spike protein. Our study demonstrates that using the inverse approach to extract a network from protein dynamics yields important biophysical insights.

Quantum modeling of dimethoxyl-indaceno dithiophene based acceptors for the development of semiconducting acceptors with outstanding photovoltaic potential

In the current DFT study, seven dimethoxyl-indaceno dithiophene based semiconducting acceptor molecules (ID1-ID7) are designed computationally by modifying the parent molecule (IDR). Here, based on a DFT exploration at a carefully selected level of theory, we have compiled a list of the optoelectronic properties of ID1-ID7 and IDR. In light of these results, all newly designed molecules, except ID5 have shown a bathochromic shift in their highest absorbance (λ _max). ID1-ID4, ID6 and ID7 molecules have smaller band gap (E _gap) and excitation energy (E _x). IP of ID5 is the smallest and EA of ID1 is the largest among all others. Compared to the parent molecule, ID1-ID3 have increased electron mobility, with ID1 being the most improved in hole mobility. ID4 had the best light harvesting efficiency in this investigation, due to its strongest oscillator. The acceptor molecules' open-circuit voltages (V _OC) were computed after being linked to the PTB7-Th donor molecule. Fill factor (FF) and normalized V _OC of ID1-ID7 were calculated and compared to the parent molecule. Based on the outcomes of this study, the modified acceptors may be further scrutinised for empirical usage in the production of organic solar cells with enhanced photovoltaic capabilities.

Biochemical Properties of Naturally Occurring Human Bloom Helicase Variants

Bloom syndrome helicase (BLM) is a RecQ-family helicase implicated in a variety of cellular processes, including DNA replication, DNA repair, and telomere maintenance. Mutations in human BLM cause Bloom syndrome (BS), an autosomal recessive disorder that leads to myriad negative health impacts including a predisposition to cancer. BS-causing mutations in BLM often negatively impact BLM ATPase and helicase activity. While BLM mutations that cause BS have been well characterized both in vitro and in vivo , there are other less studied BLM mutations that exist in the human population that do not lead to BS. Two of these non-BS mutations, encoding BLM P868L and BLM G1120R, when homozygous, increase sister chromatid exchanges in human cells. To characterize these naturally occurring BLM mutant proteins in vitro , we purified the BLM catalytic core (BLM core , residues 636-1298) with either the P868L or G1120R substitution. We also purified a BLM core K869A K870A mutant protein, which alters a lysine-rich loop proximal to the P868 residue. We found that BLM core P868L and G1120R proteins were both able to hydrolyze ATP, bind diverse DNA substrates, and unwind G-quadruplex and duplex DNA structures. Molecular dynamics simulations suggest that the P868L substitution weakens the DNA interaction with the winged-helix domain of BLM and alters the orientation of one lobe of the ATPase domain. Because BLM core P868L and G1120R retain helicase function in vitro , it is likely that the increased genome instability is caused by specific impacts of the mutant proteins in vivo . Interestingly, we found that BLM core K869A K870A has diminished ATPase activity, weakened binding to duplex DNA structures, and less robust helicase activity compared to wild-type BLM core . Thus, the lysine-rich loop may have an important role in ATPase activity and specific binding and DNA unwinding functions in BLM.

Ligands-induced open-close conformational change during DapE catalysis: Insights from molecular dynamics simulations

The microbial enzyme DapE plays a critical role in the lysine biosynthetic pathway and is considered as a potentially safe antibiotic target. In this study, atomistic simulations are employed to identify the modes of essential dynamics that define the conformational response of the enzyme to ligand binding and unbinding. The binding modes and the binding affinities of the products to the DapE enzyme are estimated from the MM-PBSA method, and the residues contributing to the ligand binding are identified. Various structural analyses and the principal component analysis of the molecular dynamics trajectories reveal that the removal of products from the active site causes a significant change in the overall enzyme structure. Both Cartesian and dihedral principal component analyses are used to characterize the structural changes in terms of domain unfolding and domain twisting motions. In the most dominant mode, that is, the domain unfolding motion, the two catalytic domains move away from the two dimerization domains of the dimeric enzyme, representing a closed-to-open conformational change. The conformational changes are initiated by the coordinated movement of three loops (Asp75-Pro82, Gly240-Asn244, and Thr347-Glu353) that trigger a domain-level movement. From multiple short trajectories, the time constant associated with the domain opening motion is estimated as 43.6 ns. Physiologically, this close-to-open conformational change is essential for the regeneration of the initial state of the enzyme for the subsequent cycle of catalytic action and provides the apo enzyme enough flexibility for efficient substrate binding.

Chasing collective variables using temporal data-driven strategies

The convergence of free-energy calculations based on importance sampling depends heavily on the choice of collective variables (CVs), which in principle, should include the slow degrees of freedom of the biological processes to be investigated. Autoencoders (AEs), as emerging data-driven dimension reduction tools, have been utilised for discovering CVs. AEs, however, are often treated as black boxes, and what AEs actually encode during training, and whether the latent variables from encoders are suitable as CVs for further free-energy calculations remains unknown. In this contribution, we review AEs and their time-series-based variants, including time-lagged AEs (TAEs) and modified TAEs, as well as the closely related model variational approach for Markov processes networks (VAMPnets). We then show through numerical examples that AEs learn the high-variance modes instead of the slow modes. In stark contrast, time series-based models are able to capture the slow modes. Moreover, both modified TAEs with extensions from slow feature analysis and the state-free reversible VAMPnets (SRVs) can yield orthogonal multidimensional CVs. As an illustration, we employ SRVs to discover the CVs of the isomerizations of N-acetyl-N'-methylalanylamide and trialanine by iterative learning with trajectories from biased simulations. Last, through numerical experiments with anisotropic diffusion, we investigate the potential relationship of time-series-based models and committor probabilities.

Structure-activity correlations for peptaibols obtained from clade Longibrachiatum of Trichoderma: A combined experimental and computational approach

Integrated disease management and plant protection have been discussed with much fervor in the past decade due to the rising environmental concerns of using industrially produced pesticides. Members of the genus Trichoderma are a subject of considerable research today due to their several properties as biocontrol agents. In our study, the peptaibol production of Trichoderma longibrachiatum SZMC 1775, T. longibrachiatum f. bissettii SZMC 12546, T. reesei SZMC 22616, T. reesei SZMC 22614, T. saturnisporum SZMC 22606 and T. effusum SZMC 22611 were investigated to elucidate structure-activity relationships (SARs) between the properties of peptaibols and their 3D structures. The effects of peptaibol mixtures obtained from every Trichoderma strain were examined against nine commonly known bacteria. The lowest minimum inhibitory concentrations (MIC, mg ml^-1) were exerted by T. longibrachiatum f. bissettii SZMC 12546 against Gram-positive bacteria, which was also able to inhibit the plant pathogenic Gram-negative Rhizobium radiobacter. Accelerated molecular dynamics (aMD) simulations were performed in aqueous solvent to explore the folding dynamics of 12 selected peptaibol sequences. The most characteristic difference between the peptaibols from group A and B relies in the 'Gly-Leu-Aib-Pro' and 'Gly-Aib-Aib-Pro' motifs ('Aib' stands for α-aminoisobutyric acid), which imparted a significant effect on the folding dynamics in water and might be correlated with their expressed bioactivity. In our aMD simulation experiments, Group A peptaibols showed more restricted folding dynamics with well-folded helical conformations as the most stable representative structures. This structural stability and dynamics may contribute to their bioactivity against the selected bacterial species.

ALS mutations in the TIA-1 prion-like domain trigger highly condensed pathogenic structures

T cell intracellular antigen-1 (TIA-1) plays a central role in stress granule (SG) formation by self-assembly via the prion-like domain (PLD). In the TIA-1 PLD, amino acid mutations associated with neurodegenerative diseases, such as amyotrophic lateral sclerosis (ALS) or Welander distal myopathy (WDM), have been identified. However, how these mutations affect PLD self-assembly properties has remained elusive. In this study, we uncovered the implicit pathogenic structures caused by the mutations. NMR analysis indicated that the dynamic structures of the PLD are synergistically determined by the physicochemical properties of amino acids in units of five residues. Molecular dynamics simulations and three-dimensional electron crystallography, together with biochemical assays, revealed that the WDM mutation E384K attenuated the sticky properties, whereas the ALS mutations P362L and A381T enhanced the self-assembly by inducing β-sheet interactions and highly condensed assembly, respectively. These results suggest that the P362L and A381T mutations increase the likelihood of irreversible amyloid fibrillization after phase-separated droplet formation, and this process may lead to pathogenicity.

Fundamental Redesign of the TIGER2hs Kernel to Address Severe Parameter Sensitivity

Replica exchange molecular dynamics simulations are one of the most popular approaches to enhance conformational sampling of molecular systems. Applications range from protein folding to protein-protein or other host-guest interactions, as well as binding free energy calculations. While these methods are computationally expensive, highly accurate results can be obtained. We recently developed TIGER2hs, an improved version of the temperature intervals with global exchange of replicas (TIGER2) algorithm. This method combines the replica-based enhanced sampling in an explicit solvent with a hybrid solvent energy evaluation. During the exchange attempts, bulk water is replaced by an implicit solvent model, allowing sampling with significantly less replicas than parallel tempering (REMD). This enables accurate enhanced sampling calculations with only a fraction of computational resources compared to REMD. Our latest results highlight several issues with sampling imbalance and parameter sensitivity within the original TIGER2 exchange algorithms that affect the overall state populations. A high sensitivity on replica number and maximum temperature is eliminated by changing to a pairwise exchange kernel (PE) without additional sorting. Simulations are controlled by adjusting the average temperature change per exchange ⟨ΔT/χ⟩ to below 30 K to mimic a controlled temperature mixing of replicas similar to REMD. Thus, this parameter provides an applicable property for selecting combinations of replica number and maximum temperature to adjust simulations for best accuracy, with flexible resource investment. This increases the robustness of the method and ensures results in excellent agreement with REMD, as demonstrated for three different peptides.

Identification of a HTT-specific binding motif in DNAJB1 essential for suppression and disaggregation of HTT

Huntington’s disease is a neurodegenerative disease caused by an expanded polyQ stretch within Huntingtin (HTT) that renders the protein aggregation-prone, ultimately resulting in the formation of amyloid fibrils. A trimeric chaperone complex composed of Hsc70, DNAJB1 and Apg2 can suppress and reverse the aggregation of HTTExon1Q₄₈. DNAJB1 is the rate-limiting chaperone and we have here identified and characterized the binding interface between DNAJB1 and HTTExon1Q₄₈. DNAJB1 exhibits a HTT binding motif (HBM) in the hinge region between C-terminal domains (CTD) I and II and binds to the polyQ-adjacent proline rich domain (PRD) of soluble as well as aggregated HTT. The PRD of HTT represents an additional binding site for chaperones. Mutation of the highly conserved H244 of the HBM of DNAJB1 completely abrogates the suppression and disaggregation of HTT fibrils by the trimeric chaperone complex. Notably, this mutation does not affect the binding and remodeling of any other protein substrate, suggesting that the HBM of DNAJB1 is a specific interaction site for HTT. Overexpression of wt DNAJB1, but not of DNAJB1^H244A can prevent the accumulation of HTTExon1Q₉₇ aggregates in HEK293 cells, thus validating the biological significance of the HBM within DNAJB1.

Ayala Mariscal et al have identified and characterized the interface of pathogenic Huntingtin and the molecular chaperone DNAJB1. Histidine-244 of the C-terminal domain of DNAJB1 is a key residues for binding to the poly-proline region of HTT. This binding site is specific for the interaction with Huntingtin.

Investigating the Bioactive Conformation of Angiotensin II Using Markov State Modeling Revisited with Web-Scale Clustering

Molecular dynamics simulation is a powerful technique for studying the structure and dynamics of biomolecules in atomic-level detail by sampling their various conformations in real time. Because of the long timescales that need to be sampled to study biomolecular processes and the big and complex nature of the corresponding data, relevant analyses of important biophysical phenomena are challenging. Clustering and Markov state models (MSMs) are efficient computational techniques that can be used to extract dominant conformational states and to connect those with kinetic information. In this work, we perform Molecular Dynamics simulations to investigate the free energy landscape of Angiotensin II (AngII) in order to unravel its bioactive conformations using different clustering techniques and Markov state modeling. AngII is an octapeptide hormone, which binds to the AT1 transmembrane receptor, and plays a vital role in the regulation of blood pressure, conservation of total blood volume, and salt homeostasis. To mimic the water-membrane interface as AngII approaches the AT1 receptor and to compare our findings with available experimental results, the simulations were performed in water as well as in water-ethanol mixtures. Our results show that in the water-ethanol environment, AngII adopts more compact U-shaped (folded) conformations than in water, which resembles its structure when bound to the AT1 receptor. For clustering of the conformations, we validate the efficiency of an inverted-quantized k-means algorithm, as a fast approximate clustering technique for web-scale data (millions of points into thousands or millions of clusters) compared to k-means, on data from trajectories of molecular dynamics simulations with reasonable trade-offs between time and accuracy. Finally, we extract MSMs using various clustering techniques for the generation of microstates and macrostates and for the selection of the macrostate representatives.

Inhibitor induced conformational changes in SARS-COV-2 papain-like protease

SARS-CoV-2's papain-like protease (PL^pro) interaction with ligands has recently been explored with a myriad of crystal structures. We used molecular dynamics (MD) simulations to study different PL^pro-ligand complexes, their ligand-induced conformational changes, and interactions. We focused on inhibitors reported with known IC₅₀ against PL^pro, namely GRL-0617, XR8-89, PLP_Snyder530, and Sander's recently published compound 7 (CPD7), and compared these trajectories against the apostructure (Apo), with a total of around 60 µs worth simulation data. We aimed to study the conformational changes using molecular dynamics simulations for the inhibitors in the PL^pro. PCA analyses and the MSM models revealed distinct conformations of PL^pro in the absence/presence of ligands and proposed that BL2-loop contributes to the accessibility of these inhibitors. Further, bulkier substituents closer to Tyr268 and Gln269 could improve inhibition of SARS-CoV-2 PL^pro by occupying the region between BL2-groove and BL2-loop, but we also expand on the relevance of exploring multiple PL^pro sub-pockets to improve inhibition.

Structural Investigations and Binding Mechanisms of Oseltamivir Drug Resistance Conferred by the E119V Mutation in Influenza H7N9 Virus

The use of vaccinations and antiviral medications have gained popularity in the therapeutic management of avian influenza H7N9 virus lately. Antiviral medicines are more popular due to being readily available. The presence of the neuraminidase protein in the avian influenza H7N9 virus and its critical role in the cleavage of sialic acid have made it a target drug in the development of influenza virus drugs. Generally, the neuraminidase proteins have common conserved amino acid residues and any mutation that occurs around or within these conserved residues affects the susceptibility and replicability of the influenza H7N9 virus. Herein, we investigated the interatomic and intermolecular dynamic impacts of the experimentally reported E119V mutation on the oseltamivir resistance of the influenza H7N9 virus. We extensively employed molecular dynamic (MD) simulations and subsequent post-MD analyses to investigate the binding mechanisms of oseltamivir-neuraminidase wildtype and E119V mutant complexes. The results revealed that the oseltamivir-wildtype complex was more thermodynamically stable than the oseltamivir-E119V mutant complex. Oseltamivir exhibited a greater binding affinity for wildtype (−15.46 ± 0.23 kcal/mol) relative to the E119V mutant (−11.72 ± 0.21 kcal/mol). The decrease in binding affinity (−3.74 kcal/mol) was consistent with RMSD, RMSF, SASA, PCA, and hydrogen bonding profiles, confirming that the E119V mutation conferred lower conformational stability and weaker protein–ligand interactions. The findings of this oseltamivir-E119V mutation may further assist in the design of compounds to overcome E119V mutation in the treatment of influenza H7N9 virus patients.

Additive eigenvectors as optimal reaction coordinates, conditioned trajectories, and time-reversible description of stochastic processes

A fundamental way to analyze complex multidimensional stochastic dynamics is to describe it as diffusion on a free energy landscape—free energy as a function of reaction coordinates (RCs). For such a description to be quantitatively accurate, the RC should be chosen in an optimal way. The committor function is a primary example of an optimal RC for the description of equilibrium reaction dynamics between two states. Here, additive eigenvectors (addevs) are considered as optimal RCs to address the limitations of the committor. An addev master equation for a Markov chain is derived. A stationary solution of the equation describes a sub-ensemble of trajectories conditioned on having the same optimal RC for the forward and time-reversed dynamics in the sub-ensemble. A collection of such sub-ensembles of trajectories, called stochastic eigenmodes, can be used to describe/approximate the stochastic dynamics. A non-stationary solution describes the evolution of the probability distribution. However, in contrast to the standard master equation, it provides a time-reversible description of stochastic dynamics. It can be integrated forward and backward in time. The developed framework is illustrated on two model systems—unidirectional random walk and diffusion.

Correlation-Based Feature Selection to Identify Functional Dynamics in Proteins

To interpret molecular dynamics simulations of biomolecular systems, systematic dimensionality reduction methods are commonly employed. Among others, this includes principal component analysis (PCA) and time-lagged independent component analysis (TICA), which aim to maximize the variance and the time scale of the first components, respectively. A crucial first step of such an analysis is the identification of suitable and relevant input coordinates (the so-called features), such as backbone dihedral angles and interresidue distances. As typically only a small subset of those coordinates is involved in a specific biomolecular process, it is important to discard the remaining uncorrelated motions or weakly correlated noise coordinates. This is because they may exhibit large amplitudes or long time scales and therefore will be erroneously considered important by PCA and TICA, respectively. To discriminate collective motions underlying functional dynamics from uncorrelated motions, the correlation matrix of the input coordinates is block-diagonalized by a clustering method. This strategy avoids possible bias due to presumed functional observables and conformational states or variation principles that maximize variance or time scales. Considering several linear and nonlinear correlation measures and various clustering algorithms, it is shown that the combination of linear correlation and the Leiden community detection algorithm yields excellent results for all considered model systems. These include the functional motion of T4 lysozyme to demonstrate the successful identification of collective motion, as well as the folding of the villin headpiece to highlight the physical interpretation of the correlated motions in terms of a functional mechanism.

Principal Component Analysis and Related Methods for Investigating the Dynamics of Biological Macromolecules

Principal component analysis (PCA) is used to reduce the dimensionalities of high-dimensional datasets in a variety of research areas. For example, biological macromolecules, such as proteins, exhibit many degrees of freedom, allowing them to adopt intricate structures and exhibit complex functions by undergoing large conformational changes. Therefore, molecular simulations of and experiments on proteins generate a large number of structure variations in high-dimensional space. PCA and many PCA-related methods have been developed to extract key features from such structural data, and these approaches have been widely applied for over 30 years to elucidate macromolecular dynamics. This review mainly focuses on the methodological aspects of PCA and related methods and their applications for investigating protein dynamics.

Folding molecular dynamics simulation of T-peptide, a HIV viral entry inhibitor: Structure, dynamics, and comparison with the experimental data

Peptide T is a synthetic octapeptide fragment, which corresponds to the region 185-192 of the gp120 HIV coat protein and functions as a viral entry inhibitor. In this work, a folding molecular dynamics simulation of peptide T in a membrane-mimicking (DMSO) solution was performed with the aim of characterizing the peptide's structural and dynamical properties. We show that peptide T is highly flexible and dynamic. The main structural characteristics observed were rapidly interconverting short helical stretches and turns, with a notable preference for the formation of β-turns. The simulation also indicated that the C-terminal part appears to be more stable than the rest of the peptide, with the most preferred conformation for residues 5-8 being a β-turn. In order to validate the accuracy of the simulations, we compared our results with the experimental NMR data obtained for the T-peptide in the same solvent. In agreement with the simulation, the NMR data indicated the presence of a preferred structure in solution that was consistent with a β-turn comprising the four C-terminal residues. An additional comparison between the experimental and simulation-derived chemical shifts also showed a reasonable agreement between experiment and simulation, further validating the simulation-derived structural characterization of the T-peptide. We conclude that peptide folding simulations produce physically relevant results even when performed with organic solvents that were not part of the force field parameterization procedure.

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.

Structure of the stress-related LHCSR1 complex determined by an integrated computational strategy

Light-harvesting complexes (LHCs) are pigment-protein complexes whose main function is to capture sunlight and transfer the energy to reaction centers of photosystems. In response to varying light conditions, LH complexes also play photoregulation and photoprotection roles. In algae and mosses, a sub-family of LHCs, light-harvesting complex stress-related (LHCSR), is responsible for photoprotective quenching. Despite their functional and evolutionary importance, no direct structural information on LHCSRs is available that can explain their unique properties. In this work, we propose a structural model of LHCSR1 from the moss P. patens, obtained through an integrated computational strategy that combines homology modeling, molecular dynamics, and multiscale quantum chemical calculations. The model is validated by reproducing the spectral properties of LHCSR1. Our model reveals the structural specificity of LHCSR1, as compared with the CP29 LH complex, and poses the basis for understanding photoprotective quenching in mosses.

The structure of the moss P. patens light-harvesting complex stress-related 1 (LHCSR1) is determined using a multi-scale computational approach for investigations of its photoprotective function.

Human Learning for Molecular Simulations: The Collective Variables Dashboard in VMD

The Collective Variables Dashboard is a software tool for real-time, seamless exploration of molecular structures and trajectories in a customizable space of collective variables. The Dashboard arises from the integration of the Collective Variables Module (also known as Colvars) with the visualization software VMD, augmented with a fully discoverable graphical interface offering interactive workflows for the design and analysis of collective variables. Typical use cases include a priori design of collective variables for enhanced sampling and free energy simulations as well as analysis of any type of simulation or collection of structures in a collective variable space. A combination of those cases commonly occurs when preliminary simulations, biased or unbiased, reveal that an optimized set of collective variables is necessary to improve sampling in further simulations. Then the Dashboard provides an efficient way to intuitively explore the space of likely collective variables, validate them on existing data, and use the resulting collective variable definitions directly in further biased simulations using the Collective Variables Module. Visualization of biasing energies and forces is proposed to help analyze or plan biased simulations. We illustrate the use of the Dashboard on two applications: discovering coordinates to describe ligand unbinding from a protein binding site and designing volume-based variables to bias the hydration of a transmembrane pore.

Probing Protein Aggregation Using the Coarse-Grained UNRES Force Field

Protein aggregation is the cause of many, often lethal, diseases, including the Alzheimer's, Parkinson's, and Huntington's diseases, and familial amyloidosis. Theoretical investigation of the mechanism of this process, including the structures of the oligomeric intermediates which are the most toxic, is difficult because of long time scale of aggregation. Coarse-grained models, which enable us to extend the simulation time scale by three or more orders of magnitude, are, therefore, of great advantage in such studies. In this chapter, we describe the application of the physics-based UNited RESidue (UNRES) force field developed in our laboratory to study protein aggregation, in both free simulations and simulations of aggregation propagation from an existing template (seed), and illustrate it with the examples of Aβ-peptide aggregation and Aβ-peptide-assisted aggregation of the peptides derived from the repeat domains of tau (TauRD).

Toward a Computational NMR Procedure for Modeling Dipeptide Side-Chain Conformation

Theoretical relationships between the vicinal spin–spin coupling constants (SSCCs) and the χ₁ torsion angles have been studied to predict the conformations of protein side chains. An efficient computational procedure is developed to obtain the conformation of dipeptides through theoretical and experimental SSCCs, Karplus equations, and quantum chemistry methods, and it is applied to three aliphatic hydrophobic residues (Val, Leu, and Ile). Three models are proposed: unimodal-static, trimodal-static-stepped, and trimodal-static-trigonal, where the most important factors are incorporated (coupled nuclei, nature and orientation of the substituents, and local geometric properties). Our results are validated by comparison with NMR and X-ray empirical data described in the literature, obtaining successful results on the 29 residues considered. Using out trimodal residue treatment, it is possible to detect and resolve residues with a simple conformation and those with two or three staggered conformers. In four residues, a deeper analysis explains that they do not have a unique conformation and that the population of each conformation plays an important role.

MLCV: Bridging Machine-Learning-Based Dimensionality Reduction and Free-Energy Calculation

Importance-sampling algorithms leaning on the definition of a model reaction coordinate (RC) are widely employed to probe processes relevant to chemistry and biology alike, spanning time scales not amenable to common, brute-force molecular dynamics (MD) simulations. In practice, the model RC often consists of a handful of collective variables (CVs) chosen on the basis of chemical intuition. However, constructing manually a low-dimensional RC model to describe an intricate geometrical transformation for the purpose of free-energy calculations and analyses remains a daunting challenge due to the inherent complexity of the conformational transitions at play. To solve this issue, remarkable progress has been made in employing machine-learning techniques, such as autoencoders, to extract the low-dimensional RC model from a large set of CVs. Implementation of the differentiable, nonlinear machine-learned CVs in common MD engines to perform free-energy calculations is, however, particularly cumbersome. To address this issue, we present here a user-friendly tool (called MLCV) that facilitates the use of machine-learned CVs in importance-sampling simulations through the popular Colvars module. Our approach is critically probed with three case examples consisting of small peptides, showcasing that through hard-coded neural network in Colvars, deep-learning and enhanced-sampling can be effectively bridged with MD simulations. The MLCV code is versatile, applicable to all the CVs available in Colvars, and can be connected to any kind of dense neural networks. We believe that MLCV provides an effective, powerful, and user-friendly platform accessible to experts and nonexperts alike for machine-learning (ML)-guided CV discovery and enhanced-sampling simulations to unveil the molecular mechanisms underlying complex biochemical processes.

A different perspective for nonphotochemical quenching in plant antenna complexes

Light-harvesting complexes of plants exert a dual function of light-harvesting (LH) and photoprotection through processes collectively called nonphotochemical quenching (NPQ). While LH processes are relatively well characterized, those involved in NPQ are less understood. Here, we characterize the quenching mechanisms of CP29, a minor LHC of plants, through the integration of two complementary enhanced-sampling techniques, dimensionality reduction schemes, electronic calculations and the analysis of cryo-EM data in the light of the predicted conformational ensemble. Our study reveals that the switch between LH and quenching state is more complex than previously thought. Several conformations of the lumenal side of the protein occur and differently affect the pigments' relative geometries and interactions. Moreover, we show that a quenching mechanism localized on a single chlorophyll-carotenoid pair is not sufficient but many chlorophylls are simultaneously involved. In such a diffuse mechanism, short-range interactions between each carotenoid and different chlorophylls combined with a protein-mediated tuning of the carotenoid excitation energies have to be considered in addition to the commonly suggested Coulomb interactions.

Fast and Accurate Multidimensional Free Energy Integration

Enhanced sampling and free energy calculation algorithms of the thermodynamic integration family (such as the adaptive biasing force (ABF) method) are not based on the direct computation of a free energy surface but rather of its gradient. Integrating the free energy surface is nontrivial in dimensions higher than one. Here, the author introduces a flexible, portable implementation of a Poisson equation formalism to integrate free energy surfaces from estimated gradients in dimensions 2 and 3 using any combination of periodic and nonperiodic (Neumann) boundary conditions. The algorithm is implemented in portable C++ and provided as a standalone tool that can be used to integrate multidimensional gradient fields estimated on a grid using any algorithm, such as umbrella integration as a post-treatment of umbrella sampling simulations. It is also included in the implementation of ABF (and its extended-system variant eABF) in the Collective Variables Module, enabling the seamless computation of multidimensional free energy surfaces within ABF and eABF simulations. A Python-based analysis toolchain is provided to easily plot and analyze multidimensional ABF simulation results, including metrics to assess their convergence. The Poisson integration algorithm can also be used to perform Helmholtz decomposition of noisy gradient estimates on the fly, resulting in an efficient implementation of the projected ABF (pABF) method proposed by Leliévre and co-workers. In numerical tests, pABF is found to lead to faster convergence with respect to ABF in simple cases of low intrinsic dimension but seems detrimental to convergence in a more realistic case involving degenerate coordinates and hidden barriers due to slower exploration. This suggests that variance reduction schemes do not always yield convergence improvements when applied to enhanced sampling methods.

Model-Free or Not?

Relaxation in nuclear magnetic resonance is a powerful method for obtaining spatially resolved, timescale-specific dynamics information about molecular systems. However, dynamics in biomolecular systems are generally too complex to be fully characterized based on NMR data alone. This is a familiar problem, addressed by the Lipari-Szabo model-free analysis, a method that captures the full information content of NMR relaxation data in case all internal motion of a molecule in solution is sufficiently fast. We investigate model-free analysis, as well as several other approaches, and find that model-free, spectral density mapping, LeMaster's approach, and our detector analysis form a class of analysis methods, for which behavior of the fitted parameters has a well-defined relationship to the distribution of correlation times of motion, independent of the specific form of that distribution. In a sense, they are all "model-free." Of these methods, only detectors are generally applicable to solid-state NMR relaxation data. We further discuss how detectors may be used for comparison of experimental data to data extracted from molecular dynamics simulation, and how simulation may be used to extract details of the dynamics that are not accessible via NMR, where detector analysis can be used to connect those details to experiments. We expect that combined methodology can eventually provide enough insight into complex dynamics to provide highly accurate models of motion, thus lending deeper insight into the nature of biomolecular dynamics.

CATBOSS: Cluster Analysis of Trajectories Based on Segment Splitting

Molecular dynamics (MD) simulations are an exceedingly and increasingly potent tool for molecular behavior prediction and analysis. However, the enormous wealth of data generated by these simulations can be difficult to process and render in a human-readable fashion. Cluster analysis is a commonly used way to partition data into structurally distinct states. We present a method that improves on the state of the art by taking advantage of the temporal information of MD trajectories to enable more accurate clustering at a lower memory cost. To date, cluster analysis of MD simulations has generally treated simulation snapshots as a mere collection of independent data points and attempted to separate them into different clusters based on structural similarity. This new method, cluster analysis of trajectories based on segment splitting (CATBOSS), applies density-peak-based clustering to classify trajectory segments learned by change detection. Applying the method to a synthetic toy model as well as four real-life data sets–trajectories of MD simulations of alanine dipeptide and valine dipeptide as well as two fast-folding proteins–we find CATBOSS to be robust and highly performant, yielding natural-looking cluster boundaries and greatly improving clustering resolution. As the classification of points into segments emphasizes density gaps in the data by grouping them close to the state means, CATBOSS applied to the valine dipeptide system is even able to account for a degree of freedom deliberately omitted from the input data set. We also demonstrate the potential utility of CATBOSS in distinguishing metastable states from transition segments as well as promising application to cases where there is little or no advance knowledge of intrinsic coordinates, making for a highly versatile analysis tool.

Molecular Mechanism of Food-Derived Polyphenols on PD-L1 Dimerization: A Molecular Dynamics Simulation Study

In cancer immunotherapy, an emerging approach is to block the interactions of programmed cell death-1 (PD-1) and programmed cell death-ligand 1 (PD-L1) using small-molecule inhibitors. The food-derived polyphenols curcumin (CC), resveratrol (RSV) and epigallocatechin gallate (EGCG) have anticancer immunologic functions, which, recently, have been proposed to act via the downregulation of PD-L1 expression. However, it remains unclear whether they can directly target PD-L1 dimerization and, thus, interrupt the PD-1/PD-L1 pathway. To elucidate the molecular mechanism of such compounds on PD-L1 dimerization, molecular docking and nanosecond molecular dynamics simulations were performed. Binding free energy calculations show that the affinities of CC, RSV and EGCG to the PD-L1 dimer follow a trend of CC > RSV > EGCG. Hence, CC is the most effective inhibitor of the PD-1/PD-L1 pathway. Analysis on contact numbers, nonbonded interactions and residue energy decomposition indicate that such compounds mainly interact with the C-, F- and G-sheet fragments of the PD-L1 dimer, which are involved in interactions with PD-1. More importantly, nonpolar interactions between these compounds and the key residues Ile54, Tyr56, Met115, Ala121 and Tyr123 play a dominant role in binding. Free energy landscape and secondary structure analyses further demonstrate that such compounds can stably interact with the binding domain of the PD-L1 dimer. The results provide evidence that CC, RSV and EGCG can inhibit PD-1/PD-L1 interactions by directly targeting PD-L1 dimerization. This provides a novel approach to discovering food-derived small-molecule inhibitors of the PD-1/PD-L1 pathway with potential applications in cancer immunotherapy.

Time-dependent principal component analysis: A unified approach to high-dimensional data reduction using adiabatic dynamics

Systematic reduction of the dimensionality is highly demanded in making a comprehensive interpretation of experimental and simulation data. Principal component analysis (PCA) is a widely used technique for reducing the dimensionality of molecular dynamics (MD) trajectories, which assists our understanding of MD simulation data. Here, we propose an approach that incorporates time dependence in the PCA algorithm. In the standard PCA, the eigenvectors obtained by diagonalizing the covariance matrix are time independent. In contrast, they are functions of time in our new approach, and their time evolution is implemented in the framework of Car-Parrinello or Born-Oppenheimer type adiabatic dynamics. Thanks to the time dependence, each of the step-by-step structural changes or intermittent collective fluctuations is clearly identified, which are often keys to provoking a drastic structural transformation but are easily masked in the standard PCA. The time dependence also allows for reoptimization of the principal components (PCs) according to the structural development, which can be exploited for enhanced sampling in MD simulations. The present approach is applied to phase transitions of a water model and conformational changes of a coarse-grained protein model. In the former, collective dynamics associated with the dihedral-motion in the tetrahedral network structure is found to play a key role in crystallization. In the latter, various conformations of the protein model were successfully sampled by enhancing structural fluctuation along the periodically optimized PC. Both applications clearly demonstrate the virtue of the new approach, which we refer to as time-dependent PCA.

A deep autoencoder framework for discovery of metastable ensembles in biomacromolecules

Biomacromolecules manifest dynamic conformational fluctuation and involve mutual interconversion among metastable states. A robust mapping of their conformational landscape often requires the low-dimensional projection of the conformational ensemble along optimized collective variables (CVs). However, the traditional choice for the CV is often limited by user-intuition and prior knowledge about the system, and this lacks a rigorous assessment of their optimality over other candidate CVs. To address this issue, we propose an approach in which we first choose the possible combinations of inter-residue Cα-distances within a given macromolecule as a set of input CVs. Subsequently, we derive a non-linear combination of latent space embedded CVs via auto-encoding the unbiased molecular dynamics simulation trajectories within the framework of the feed-forward neural network. We demonstrate the ability of the derived latent space variables in elucidating the conformational landscape in four hierarchically complex systems. The latent space CVs identify key metastable states of a bead-in-a-spring polymer. The combination of the adopted dimensional reduction technique with a Markov state model, built on the derived latent space, reveals multiple spatially and kinetically well-resolved metastable conformations for GB1 β-hairpin. A quantitative comparison based on the variational approach-based scoring of the auto-encoder-derived latent space CVs with the ones obtained via independent component analysis (principal component analysis or time-structured independent component analysis) confirms the optimality of the former. As a practical application, the auto-encoder-derived CVs were found to predict the reinforced folding of a Trp-cage mini-protein in aqueous osmolyte solution. Finally, the protocol was able to decipher the conformational heterogeneities involved in a complex metalloenzyme, namely, cytochrome P450.