Autoencoder-based detection of the residues involved in G protein-coupled receptor signaling.

Regulator binding and mutations alter protein dynamics. The transmission of the signal of these alterations to distant sites through protein motion results in changes in protein expression and cell function. The detection of residues involved in signal transmission contributes to an elucidation of the mechanisms underlying processes as vast as cellular function and disease pathogenesis. We developed an autoencoder (AE) based method that detects residues essential for signaling by comparing the fluctuation data, particularly the time fluctuation of the side-chain distances between residues, during molecular dynamics simulations between the ligand-bound and -unbound forms or wild-type and mutant forms of proteins. Here, the AE-based method was applied to the G protein-coupled receptor (GPCR) system, particularly a class A-type GPCR, CXCR4, to detect the essential residues involved in signaling. Among the residues involved in the signaling of the homolog CXCR2, which were extracted from the literature based on the complex structures of the ligand and G protein, our method could detect more than half of the essential residues involved in G protein signaling, including those spanning the fifth and sixth transmembrane helices in the intracellular region, despite the lack of information regarding the interaction with G protein in our CXCR4 models.

Mutual information analysis of the dynamic correlation between side chains in proteins.

Protein dynamics play an essential role in function regulation. In recent years, many experimental and theoretical studies have shown that changes in protein fluctuations in the backbone and side chains fulfill a pivotal role associated with amino acid mutations, chemical modifications, and ligand binding. The dynamic correlations between protein side chains have not been sufficiently studied, and no reliable analysis method has been available so far. Therefore, we developed a method to evaluate the dynamic correlation between protein side chains using mutual information and molecular dynamics simulations. To eliminate the structural superposition errors dealing with conventional analysis methods, and to accurately extract the intrinsic fluctuation properties of the side chains, we employed distance principal component analysis (distPCA). The motion of the side chain was then projected onto the eigenvector space obtained by distPCA, and the mutual information between the projected motions was calculated. The proposed method was then applied to a small protein "eglin c" and the mutants. The results show that even a single mutation significantly changed the dynamic correlations and also suggest that the dynamic change is deeply related to the stability. Those results indicate that our developed method could be useful for analyzing the molecular mechanism of allosteric communication in proteins.

Path Ensembles for Pin1-Catalyzed Cis-Trans Isomerization of a Substrate Calculated by Weighted Ensemble Simulations.

Pin1 enzyme protein recognizes specifically phosphorylated serine/threonine (pSer/pThr) and catalyzes the slow interconversion of the peptidyl-prolyl bond between cis and trans forms. Structural dynamics between the cis and trans forms are essential to reveal the underlying molecular mechanism of the catalysis. In this study, we apply the weighted ensemble (WE) simulation method to obtain comprehensive path ensembles for the Pin1-catalyzed isomerization process. Associated rate constants for both cis-to-trans and trans-to-cis isomerization are calculated to be submicroseconds time scales, which are in good agreement with the calculated free energy landscape where the cis form is slightly less favorable. The committor-like analysis indicates the shift of the transition state toward trans form (at the isomerization angle ω ∼ 110°) compared to the intrinsic position for the isolated substrate (ω ∼ 90°). The calculated structural ensemble clarifies a role of both the dual-histidine motif, His59/His157, and the basic residues, Lys63/Arg68/Arg69, to anchor both sides of the peptidyl-prolyl bond, the aromatic ring in Pro, and the phosphate in pSer, respectively. The rotation of the torsion angle is found to be facilitated by relaying the hydrogen-bond partner of the main-chain oxygen in pSer from Cys113 in the cis form to Arg68 in the trans form, through Ser154 at the transition state, which is really the cause of the shift in the transition state. The role of Ser154 as a driving force of the isomerization is confirmed by additional WE and free energy calculations for S154A mutant where the isomerization takes place slightly slower and the free energy barrier increases through the mutation. The present study shows the usefulness of the WE simulation for substantial path samplings between the reactant and product states, unraveling the molecular mechanism of the enzyme catalysis.

Rebuilding Ring-Type Assembly of Peroxiredoxin by Chemical Modification.

Direct control of the protein quaternary structure (QS) is challenging owing to the complexity of the protein structure. As a protein with a characteristic QS, peroxiredoxin from Aeropyrum pernix K1 (ApPrx) forms a decamer, wherein five dimers associate to form a ring. Here, we disrupted and reconstituted ApPrx QS via amino acid mutations and chemical modifications targeting hot spots for protein assembly. The decameric QS of an ApPrx* mutant, wherein all cysteine residues in wild-type ApPrx were mutated to serine, was destructed to dimers via an F80C mutation. The dimeric ApPrx*F80C mutant was then modified with a small molecule and successfully assembled as a decamer. Structural analysis confirmed that an artificially installed chemical moiety potentially facilitates suitable protein-protein interactions to rebuild a native structure. Rebuilding of dodecamer was also achieved through an additional amino acid mutation. This study describes a facile method to regulate the protein assembly state.

Atomistic detailed free-energy landscape of intrinsically disordered protein studied by multi-scale divide-and-conquer molecular dynamics simulation.

Calcineurin (CaN) is a eukaryotic serine/threonine protein phosphatase activated by both Ca2+ and calmodulin (CaM), including intrinsically disordered region (IDR). The region undergoes folding into an α-helix form in the presence Ca2+ -loaded CaM. To sample the ordered structure of the IDR by conventional all atom model (AAM) molecular dynamics (MD) simulation, the IDR and Ca2+ -loaded CaM must be simultaneously treated. However, it is time-consuming task because the coupled folding and binding should include repeated binding and dissociation. Then, in this study, we propose novel multi-scale divide-and-conquer MD (MSDC-MD), which combines AAM-MD and coarse-grained model MD (CGM-MD). To speed up the conformation sampling, MSDC-MD simulation first treats the IDR by CGM to sample conformations from wide conformation space; then, multiple AAM-MD in a limited area is initiated using the resultant CGM conformation, which is reconstructed by homology modeling method. To investigate performance, we sampled the ordered conformation of the IDR using MSDC-MD; the root-mean-square distance (RMSD) with respect to the experimental structure was 2.23 Å.

Epigenetic Protection of Vertebrate Lymphoid Progenitor Cells by Dnmt1.

DNA methylation is a universal epigenetic mechanism involved in regulation of gene expression and genome stability. The DNA maintenance methylase DNMT1 ensures that DNA methylation patterns are faithfully transmitted to daughter cells during cell division. Because loss of DNMT1 is lethal, a pan-organismic analysis of DNMT1 function is lacking. We identified new recessive dnmt1 alleles in medaka and zebrafish and, guided by the structures of mutant proteins, generated a recessive variant of mouse Dnmt1. Each of the three missense mutations studied here distorts the catalytic pocket and reduces enzymatic activity. Because all three DNMT1 mutant animals are viable, it was possible to examine their phenotypes throughout life. The consequences of genome-wide hypomethylation of DNA of somatic tissues in the Dnmt1 mutants are surprisingly mild but consistently affect the development of the lymphoid lineage. Our findings indicate that developing lymphocytes in vertebrates are sensitive to perturbations of DNA maintenance methylation.

Loosening of Side-Chain Packing Associated with Perturbations in Peripheral Dynamics Induced by the D76N Mutation of ß2-Microglobulin Revealed by Pressure-NMR and Molecular Dynamic Simulations.

ß2-Microglobulin (ß2m) is the causative protein of dialysis-related amyloidosis, and its D76N variant is less stable and more prone to aggregation. Since their crystal structures are indistinguishable from each other, enhanced amyloidogenicity induced by the mutation may be attributed to changes in the structural dynamics of the molecule. We examined pressure and mutation effects on the ß2m molecule by NMR and MD simulations, and found that the mutation induced the loosening of the inter-sheet packing of ß2m, which is relevant to destabilization and subsequent amyloidogenicity. On the other hand, this loosening was coupled with perturbed dynamics at some peripheral regions. The key result for this conclusion was that both the mutation and pressure induced similar reductions in the mobility of these residues, suggesting that there is a common mechanism underlying the suppression of inherent fluctuations in the ß2m molecule. Analyses of data obtained under high pressure conditions suggested that the network of dynamically correlated residues included not only the mutation site, but also distal residues, such as those of the C- and D-strands. Reductions in these local dynamics correlated with the loosening of inter-sheet packing.

Autoencoder-Based Detection of Dynamic Allostery Triggered by Ligand Binding Based on Molecular Dynamics.

Dynamic allostery on proteins, triggered by regulator binding or chemical modifications, transmits information from the binding site to distant regions, dramatically altering protein function. It is accompanied by subtle changes in side-chain conformations of the protein, indicating that the changes in dynamics, and not rigid or large conformational changes, are essential to understand regulation of protein function. Although a lot of experimental and theoretical studies have been dedicated to investigate this issue, the regulation mechanism of protein function is still being debated. Here, we propose an autoencoder-based method that can detect dynamic allostery. The method is based on the comparison of time fluctuations of protein structures, in the form of distance matrices, obtained from molecular dynamics simulations in ligand-bound and -unbound forms. Our method detected that the changes in dynamics by ligand binding in the PDZ2 domain led to the reorganization of correlative fluctuation motions among residue pairs, which revealed a different view of the correlated motions from the PCA and DCCM. In addition, other correlative motions were also found as a result of the dynamic perturbation from the ligand binding, which may lead to dynamic allostery. This autoencoder-based method would be usefully applied to the signal transduction and mutagenesis systems involved in protein functions and severe diseases.

Mutational effects of Cys113 on structural dynamics of Pin1.

Pin1 is a peptidyl-prolyl isomerase (PPIase) which catalyzes cis/trans isomerization of pS/pT-P bond. Its activity is related to various cellular functions including suppression of Alzheimer's disease. A cysteine residue C113 is known to be important for its PPIase activity; a mutation C113A reduced the activity by 130-fold. According to various nuclear magnetic resonance experiments for mutants of C113 and molecular dynamics (MD) simulation of wild-type Pin1, the protonation sate of Sγ of C113 regulates the hydrogen-bonding network of the dual-histidine motif (H59, H157) whose dynamics may affect substrate binding ability. However, it was still unclear why such local dynamic changes altered the PPIase activity of Pin1. In this study, we performed 500 ns of MD simulations of full-length wild-type Pin1 and C113A mutant in order to elucidate why the mutation C113A drastically reduced the PPIase activity of Pin1. The principal component analysis for both MD trajectories clearly elucidated that the mutation C113A suppressed the dynamics of Pin1 because it stabilized a hydrogen-bond between Nɛ of H59 and Oγ of S115. In the dynamics of wild-type protein, the phosphate binding loop (K63-S71) as well as the interdomain hinge showed the closed-open dynamics which correlated with the change of the hydrogen-bonding network of the dual-histidine motif. In contrast, in the dynamics of C113A mutant, the phosphate binding loop took only the closed conformation together with the interdomain hinge. Such closed-open dynamics must be essential for the PPIase activity of Pin1.

Allosteric activation of cytochrome P450 3A4 by efavirenz facilitates midazolam binding.

1. The purpose of this study is to investigate the heteroactivation mechanism of CYP3A4 by efavirenz, which enhances metabolism of midazolam in vivo, in terms of its binding to CYP3A4 with in vitro spectroscopic methods. 2. Efavirenz exhibited a type II spectral change with binding to CYP3A4 indicating a possible inhibitor. Although dissociation constant (K d) was approximated as 520 µM, efavirenz enhanced binding affinity of midazolam as a co-existing drug with an estimated iK d value of 5.6 µM which is comparable to a clinical concentration. 3. Efavirenz stimulated the formation of 1'-hydroxymidazolam, and the product formation rate (V max) concentration-dependently increased without changing the K m. Besides, an efavirenz analogue, [6-chloro-1,4-dihydro-4-(1-pentynyl)-4-(trifluoromethyl)-2H-3,1-benzoxazin-2-one] (efavirenz impurity) slightly facilitated the binding affinity of midazolam in a concentration-dependent manner. These results propose that efavirenz affects midazolam-binding via binding to the peripheral site which is apart from the active site of CYP3A4. 4. A molecular dynamics simulation also suggested the bound-efavirenz was repositioned to effector-binding site. As a consequence, our spectroscopic studies clarified the heteroactivation of CYP3A4 caused by efavirenz with a proper affinity to the peripheral site, and we concluded the method can be a useful tool for characterising the potential for drug-drug interactions.

On-the-fly analysis of molecular dynamics simulation trajectories of proteins using the Bayesian inference method.

Robust and reliable analyses of long trajectories from molecular dynamics simulations are important for investigations of functions and mechanisms of proteins. Structural fitting is necessary for various analyses of protein dynamics, thus removing time-dependent translational and rotational movements. However, the fitting is often difficult for highly flexible molecules. Thus, to address the issues, we proposed a fitting algorithm that uses the Bayesian inference method in combination with rotational fitting-weight improvements, and the well-studied globular protein systems trpcage and lysozyme were used for investigations. The present method clearly identified rigid core regions that fluctuate less than other regions and also separated core regions from highly fluctuating regions with greater accuracy than conventional methods. Our method also provided simultaneous variance-covariance matrix elements composed of atomic coordinates, allowing us to perform principle component analysis and prepare domain cross-correlation map during molecular dynamics simulations in an on-the-fly manner.

Free Energy Reconstruction from Logarithmic Mean-Force Dynamics Using Multiple Nonequilibrium Trajectories.

Efficient and reliable estimation of the mean force (MF), the derivatives of the free energy with respect to a set of collective variables (CVs), has been a challenging problem because free energy differences are often computed by integrating the MF. Among various methods for computing free energy differences, logarithmic mean-force dynamics (LogMFD) [ Morishita et al., Phys. Rev. E 2012 , 85 , 066702 ] invokes the conservation law in classical mechanics to integrate the MF, which allows us to estimate the free energy profile along the CVs on-the-fly. Here, we present a method called parallel dynamics, which improves the estimation of the MF by employing multiple replicas of the system and is straightforwardly incorporated in LogMFD or a related method. In the parallel dynamics, the MF is evaluated by a nonequilibrium path-ensemble using the multiple replicas based on the Crooks-Jarzynski nonequilibrium work relation. Thanks to the Crooks relation, realizing full-equilibrium states is no longer mandatory for estimating the MF. Additionally, sampling in the hidden subspace orthogonal to the CV space is highly improved with appropriate weights for each metastable state (if any), which is hardly achievable by typical free energy computational methods. We illustrate how to implement parallel dynamics by combining it with LogMFD, which we call logarithmic parallel dynamics (LogPD). Biosystems of alanine dipeptide and adenylate kinase in explicit water are employed as benchmark systems to which LogPD is applied to demonstrate the effect of multiple replicas on the accuracy and efficiency in estimating the free energy profiles using parallel dynamics.

Substrate recognition of N,N'-diacetylchitobiose deacetylase from Pyrococcus horikoshii.

Enzymes of carbohydrate esterase (CE) family 14 catalyze hydrolysis of N-acetyl groups at the non-reducing end of the N-acetylglucosamine (GlcNAc) residue of chitooligosaccharides or related compounds. N,N'-diacetylchitobiose deacetylase (Dac) belongs to the CE-14 family and plays a role in the chitinolytic pathway in archaea by deacetylating N,N'-diacetylchitobiose (GlcNAc2), which is the end product of chitinase. In this study, we revealed the structural basis of reaction specificity in CE-14 deacetylases by solving a crystal structure of Dac from Pyrococcus horikoshii (Ph-Dac) in complex with a novel reaction intermediate analog. We developed 2-deoxy-2-methylphosphoramido-d-glucose (MPG) as the analog of the tetrahedral oxyanion intermediate of the monosaccharide substrate GlcNAc. The crystal structure of Ph-Dac in complex with MPG demonstrated that Arg92, Asp115, and His152 side chains interact with hydroxyl groups of the glucose moiety of the non-reducing-end GlcNAc residue. The amino acid residues responsible for recognition of the MPG glucose moiety are spatially conserved in other CE-14 deacetylases. Molecular dynamics simulation of the structure of the Ph-Dac-GlcNAc2 complex indicated that the reducing GlcNAc residue is placed in a large intermolecular cleft and is not involved with specific interactions with the enzyme. This observation was consistent with results indicating that Ph-Dac displayed similar kinetic parameters for both GlcNAc and GlcNAc2. This study provides the structural basis of reaction-site specificity of Dac and related CE-14 enzymes.

Development of massive multilevel molecular dynamics simulation program, Platypus (PLATform for dYnamic Protein Unified Simulation), for the elucidation of protein functions.

A massively parallel program for quantum mechanical-molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc-pVDZ and B3LYP/cc-pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6-31G** calculations. We also performed excited QM/MM-MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH-insensitive and photo-stable ultramarine fluorescent protein. Platypus accelerated on-the-fly excited-state QM/MM-MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50-ps (200,000-step) on-the-fly excited-state QM/MM-MD simulations for the SIRIUS in water.

Direct Interaction between the Voltage Sensors Produces Cooperative Sustained Deactivation in Voltage-gated H+ Channel Dimers.

The voltage-gated H(+) channel (Hv) is a voltage sensor domain-like protein consisting of four transmembrane segments (S1-S4). The native Hv structure is a homodimer, with the two channel subunits functioning cooperatively. Here we show that the two voltage sensor S4 helices within the dimer directly cooperate via a π-stacking interaction between Trp residues at the middle of each segment. Scanning mutagenesis showed that Trp situated around the original position provides the slow gating kinetics characteristic of the dimer's cooperativity. Analyses of the Trp mutation on the dimeric and monomeric channel backgrounds and analyses with tandem channel constructs suggested that the two Trp residues within the dimer are functionally coupled during Hv deactivation but are less so during activation. Molecular dynamics simulation also showed direct π-stacking of the two Trp residues. These results provide new insight into the cooperative function of voltage-gated channels, where adjacent voltage sensor helices make direct physical contact and work as a single unit according to the gating process.

A method for predicting protein conformational pathways by using molecular dynamics simulations guided by difference distance matrices.

Here, an efficient method that predicts natural transition pathways between two endpoint states of an allosteric protein has been proposed. This method helps create structures that bridge these endpoints through multiple iterative and unbiased molecular dynamics simulations with explicit water. Difference distance matrices provide an approach for identifying states involving concerted slow motion. A series of structures are readily generated along the transition pathways of adenylate kinase. Predicted structures may be useful for an initial pathway to evaluate free energy landscapes via umbrella sampling and chain-of-states methods. © 2016 Wiley Periodicals, Inc.

Molecular dynamics study of the phosphorylation effect on the conformational states of the C-terminal domain of RNA polymerase II.

The carboxyl-terminal domain (CTD) of RNA polymerase II in eukaryotes regulates mRNA processing processes by recruiting various regulation factors. A main function of the CTD relies on the heptad consensus sequence (YSPTSPS). The CTD dynamically changes its conformational state to recognize and bind different regulation factors. The dynamical conformation changes are caused by modifications, mainly phosphorylation and dephosphorylation, to the serine residues. In this study, we investigate the conformational states of the unit consensus CTD peptide with various phosphorylation patterns of the serine residues by extended ensemble simulations. The results show that the CTD without phosphorylation has a flexible disordered structure distributed between twisted and extended states, but phosphorylation tends to reduce the conformational space. It was found that phosphorylation induces a ß-turn around the phosphorylated serine residue and the cis conformation of the proline residue significantly inhibits the ß-turn formation. The ß-turn should contribute to specific CTD binding of the different regulation factors by changing the conformation propensity combined with induced fit.

Simple and accurate scheme to compute electrostatic interaction: zero-dipole summation technique for molecular system and application to bulk water.

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 Å for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.

A long-range electrostatic potential based on the Wolf method charge-neutral condition.

Molecular simulations rely heavily on a long range electrostatic Coulomb interaction. The Coulomb potential decays inversely with distance, indicating infinite effective range. In practice, molecular simulations do not directly take into account such an infinite interaction. Therefore, the Ewald, fast multipole, and cutoff methods are frequently used. Although cutoff methods are implemented easily and the calculations are fast, it has been pointed out that they produce serious artifacts. Wolf and coworkers recently discovered one source of the artifacts. They found that when the total charge in a cutoff sphere disappeared, the cutoff error is dramatically suppressed. The Wolf method uses the charge-neutral principle combined with a potential damping that is realized using a complementary error function. To date, many molecular simulation studies have demonstrated the accuracy and reliability of the Wolf method. We propose a novel long-range potential that is constructed only from the charge-neutral condition of the Wolf method without potential damping. We also show that three simulation systems, in which involve liquid sodium-chloride, TIP3P water, and a charged protein in explicit waters with neutralized ions using the new potential, provide accurate statistical and dielectric properties when compared with the particle mesh Ewald method.

Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle.

We propose a novel idea, zero-dipole summation, for evaluating the electrostatic energy of a classical particle system, and have composed an algorithm for effectively utilizing the idea for molecular dynamics. It conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. The resulting energy formula is nevertheless represented by a simple pairwise function sum, which enables facile application to high-performance computation. By following a heuristic approach to derive the current electrostatic energy formula, we developed an axiomatic approach to construct the method consistently. Explorations of the theoretical details of our method revealed the structure of the generated error, and we analyzed it by comparisons with other methods. A numerical simulation using liquid sodium chloride confirmed that the current method with a small damping factor yielded sufficient accuracy with a practical cutoff distance region. The current energy function also conducts stable numerical integration in a liquid MD simulation. Our method is an extension of the charge neutralized summation developed by Wolf et al. [J. Chem. Phys. 110, 8254 (1999)]. Furthermore, we found that the current method becomes a generalization of the preaveraged potential method proposed by Yakub and Ronchi [J. Chem. Phys. 119, 11556 (2003)], which is based on a viewpoint different from the neutrality. The current study presents these relationships and suggests possibilities for their further applications.