To establish a procedure for screening compounds that inhibit ligand–receptor binding, we employed a multidimensional virtual-system coupled molecular dynamics (mD-VcMD), a generalized ensemble method recently developed by our group. In this approach, each compound is initially placed far from the receptor. Both receptor and compound were fully flexible in explicit solvent during sampling. The mD-VcMD generated a free-energy landscape of the compound–receptor interactions, in which a probability of existence was assigned to each sampled conformation. We examined four compounds that bind to the papain-like protease (PLpro) of SARS-CoV-2. The resultant free-energy landscapes were funnel-like for all compounds. The probabilities assigned to the free-energy basins correlated well with the measured dissociation constants. Furthermore, structural clustering revealed two types of binding modes within the free-energy basin. The probabilities assigned to the binding modes correlated well with the measured enzyme inhibitory activity. These results suggest that the proposed procedure is effective for selecting candidate inhibitors among the examined compounds.

Ligand-receptor docking simulation is difficult when the biomolecules have high intrinsic flexibility. If some knowledge on the ligand-receptor complex structure or inter-molecular contact sites are presented in advance, the difficulty of docking problem considerably decreases. This paper proposes a generalized-ensemble method "cartesian-space division mD-VcMD" (or CSD-mD-VcMD), which calculates stable complex structures without assist of experimental knowledge on the complex structure. This method is an extension of our previous method that requires the knowledge on the ligand-receptor complex structure in advance. Both the present and previous methods enhance the conformational sampling, and finally produce a binding free-energy landscape starting from a completely dissociated conformation, and provide a free-energy landscape. We applied the present method to same system studied by the previous method: A ligand (ribocil A or ribocil B) binding to an RNA (the aptamer domain of the FMN riboswitch). The two methods produced similar results, which explained experimental data. For instance, ribocil B bound to the aptamer's deep binding pocket more strongly than ribocil A did. However, this does not mean that two methods have a similar performance. Note that the present method did not use the experimental knowledge of binding sites although the previous method was supported by the knowledge. The RNA-ligand binding site could be a cryptic site because RNA and ligand are highly flexible in general. The current study showed that CSD-mD-VcMD is actually useful to obtain a binding free-energy landscape of a flexible system, i.e., the RNA-ligand interacting system.

To investigate the effects of phosphorylation on the function of the human positive cofactor 4 (PC4), an enhanced molecular dynamics (MD) simulation was performed. The simulation system consists of the N-terminal intrinsic disordered region (IDR) of PC4 and a complex that comprises the C-terminal acidic activation domain of a herpes simplex virion protein 16 (VP16ad) and a homodimer of the C-terminal structured core domain of PC4 (PC4ctd). An earlier report of an experimental study reported that the PC4-VP16ad interaction is modulated by incremental phosphorylation of the IDR. That report also proposed a dynamic model where phosphorylated serine residues of a segment (SEAC) in the IDR contact positively charged residues (lysin and arginine) of another segment (K-rich) in the IDR. This contact formation induced by the phosphorylation results in variation of PC4-VP16ad interaction. However, this contact formation has not yet been measured directly because it is transiently formed and because the SEAC and K-rich segments are unstructured with high flexibility. We performed two simulations to mimic the incremental phosphorylation: the IDR was not phosphorylated in one simulation and only partially phosphorylated in the other. Our simulation results indicate that the phosphorylation weakens the IDR-VP16ad contact considerably with the induction of a compact structure in the IDR. This structure was stabilized by electrostatic interactions between the phosphorylated serine residues of a segment and lysine or arginine residues of another segment in the IDR, but the conformational fluctuation of this compact structure was considerably large. Consequently, the present study supports the experimentally proposed dynamic model. Results of this study can be important for computational elucidation of the functional modulation of PC4.

For the design and development of innovative carbon nanotube (CNT)-based tools and applications, an understanding of the molecular interactions between CNTs and biological systems is essential. In this study, a three-dimensional protein-structure-based in silico screen identified the paired immune receptors, sialic acid immunoglobulin-like binding lectin-5 (Siglec-5) and Siglec-14, as CNT-recognizing receptors. Molecular dynamics simulations showed the spatiotemporally stable association of aromatic residues on the extracellular loop of Siglec-5 with CNTs. Siglec-14 mediated spleen tyrosine kinase (Syk)-dependent phagocytosis of multiwalled CNTs and the subsequent secretion of interleukin-1β from human monocytes. Ectopic in vivo expression of human Siglec-14 on mouse alveolar macrophages resulted in enhanced recognition of multiwalled CNTs and exacerbated pulmonary inflammation. Furthermore, fostamatinib, a Syk inhibitor, blocked Siglec-14-mediated proinflammatory responses. These results indicate that Siglec-14 is a human activating receptor recognizing CNTs and that blockade of Siglec-14 and the Syk pathway may overcome CNT-induced inflammation.

The pivotal importance of intrinsically disordered proteins (IDPs) in cellular biochemistry is widely accepted due to extensive studies, especially on liquid-liquid phase separation (LLPS). Elucidating their molecular mechanisms is essential for the comprehensive analysis of the molecular and structural biology of the cell. However, the high flexibility and dynamic features of IDPs make it challenging to analyze their structures. A molecular dynamics (MD) simulation is a promising approach to investigate highly complicated interactions among IDPs and their complexes. This chapter presents an overview of recent advances in MD simulations for higher-order complex formation of IDPs and LLPS using a coarse-grained model.

Nonstructural protein 1 (nsp1) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a 180-residue protein that blocks translation of host mRNAs in SARS-CoV-2-infected cells. Although it is known that SARS-CoV-2’s own RNA evades nsp1’s host translation shutoff, the molecular mechanism underlying the evasion was poorly understood. We performed an extended ensemble molecular dynamics simulation to investigate the mechanism of the viral RNA evasion. Simulation results suggested that the stem loop structure of the SARS-CoV-2 RNA 5’-untranslated region (SL1) binds to both nsp1’s N-terminal globular region and intrinsically disordered region. The consistency of the results was assessed by modeling nsp1-40S ribosome structure based on reported nsp1 experiments, including the X-ray crystallographic structure analysis, the cryo-EM electron density map, and cross-linking experiments. The SL1 binding region predicted from the simulation was open to the solvent, yet the ribosome could interact with SL1. Cluster analysis of the binding mode and detailed analysis of the binding poses suggest residues Arg124, Lys47, Arg43, and Asn126 may be involved in the SL1 recognition mechanism, consistent with the existing mutational analysis.

Several basic leucine zipper (bZIP) transcription factors have accessory motifs in their DNA-binding domains, such as the CNC motif of CNC family or the EHR motif of small Maf (sMaf) proteins. CNC family proteins heterodimerize with sMaf proteins to recognize CNC-sMaf binding DNA elements (CsMBEs) in competition with sMaf homodimers, but the functional role of the CNC motif remains elusive. In this study, we report the crystal structures of Nrf2/NFE2L2, a CNC family protein regulating anti-stress transcriptional responses, in a complex with MafG and CsMBE. The CNC motif restricts the conformations of crucial Arg residues in the basic region, which form extensive contact with the DNA backbone phosphates. Accordingly, the Nrf2-MafG heterodimer has approximately a 200-fold stronger affinity for CsMBE than canonical bZIP proteins, such as AP-1 proteins. The high DNA affinity of the CNC-sMaf heterodimer may allow it to compete with the sMaf homodimer on target genes without being perturbed by other low-affinity bZIP proteins with similar sequence specificity.

Plant sphingolipids mostly possess 2-hydroxy fatty acids (HFA), the synthesis of which is catalyzed by FA 2-hydroxylases (FAHs). In Arabidopsis (Arabidopsis thaliana), two FAHs (FAH1 and FAH2) have been identified. However, the functions of FAHs and sphingolipids with HFAs (2-hydroxy sphingolipids) are still unknown because of the lack of Arabidopsis lines with the complete deletion of FAH1. In this study, we generated a FAH1 mutant (fah1c) using CRISPR/Cas9-based genome editing. Sphingolipid analysis of fah1c, fah2, and fah1cfah2 mutants revealed that FAH1 hydroxylates very long-chain FAs (VLCFAs), whereas the substrates of FAH2 are VLCFAs and palmitic acid. However, 2-hydroxy sphingolipids are not completely lost in the fah1cfah2 double mutant, suggesting the existence of other enzymes catalyzing the hydroxylation of sphingolipid FAs. Plasma membrane (PM) analysis and molecular dynamics simulations revealed that hydroxyl groups of sphingolipid acyl chains play a crucial role in the organization of nanodomains, which are nanoscale liquid-ordered domains mainly formed by sphingolipids and sterols in the PM, through hydrogen bonds. In the PM of the fah1cfah2 mutant, the expression levels of 26.7% of the proteins, including defense-related proteins such as the pattern recognition receptors (PRRs) brassinosteroid insensitive 1-associated receptor kinase 1 and chitin elicitor receptor kinase 1, NADPH oxidase respiratory burst oxidase homolog D (RBOHD), and heterotrimeric G proteins, were lower than that in the wild-type. In addition, reactive oxygen species (ROS) burst was suppressed in the fah1cfah2 mutant after treatment with the pathogen-associated molecular patterns flg22 and chitin. These results indicated that 2-hydroxy sphingolipids are necessary for the organization of PM nanodomains and ROS burst through RBOHD and PRRs during pattern-triggered immunity.

Elucidating the principles of sequence–structure relationships of proteins is a long-standing issue in biology. The nature of a short segment of a protein is determined by both the subsequence of the segment itself and its environment. For example, a type of subsequence, the so-called chameleon sequences, can form different secondary structures depending on its environments. Chameleon sequences are considered to have a weak tendency to form a specific structure. Although many chameleon sequences have been identified, they are only a small part of all possible subsequences in the proteome. The strength of the tendency to take a specific structure for each subsequence has not been fully quantified. In this study, we comprehensively analyzed subsequences consisting of four to nine amino acid residues, or N-gram (4≤N≤9), observed in non-redundant sequences in the Protein Data Bank (PDB). Tendencies to form a specific structure in terms of the secondary structure and accessible surface area are quantified as information quantities for each N-gram. Although the majority of observed subsequences have low information quantity due to lack of samples in the current PDB, thousands of N-grams with strong tendencies, including known structural motifs, were found. In addition, machine learning partially predicted the tendency of unknown N-grams, and thus, this technique helps to extract knowledge from the limited number of samples in the PDB.

Prediction of ligand-receptor complex structure is important in both the basic science and the industry such as drug discovery. We report various computation molecular docking methods: fundamental in silico (virtual) screening, ensemble docking, enhanced sampling (generalized ensemble) methods, and other methods to improve the accuracy of the complex structure. We explain not only the merits of these methods but also their limits of application and discuss some interaction terms which are not considered in the in silico methods. In silico screening and ensemble docking are useful when one focuses on obtaining the native complex structure (the most thermodynamically stable complex). Generalized ensemble method provides a free-energy landscape, which shows the distribution of the most stable complex structure and semi-stable ones in a conformational space. Also, barriers separating those stable structures are identified. A researcher should select one of the methods according to the research aim and depending on complexity of the molecular system to be studied.

A GA-guided multidimensional virtual-system coupled molecular dynamics (GA-mD-VcMD) simulation was conducted to elucidate binding mechanisms of a middle-sized flexible molecule, bosentan, to a GPCR protein, human endothelin receptor type B (hETB). GA-mD-VcMD is a generalized ensemble method that produces a free-energy landscape of the ligand-receptor binding by searching large-scale motions accompanied with stable maintenance of the fragile cell-membrane structure. All molecular components (bosentan, hETB, membrane, and solvent) were represented with an all-atom model. Then sampling was conducted from conformations where bosentan was distant from the binding site in the hETB binding pocket. The deepest basin in the resultant free-energy landscape was assigned to native-like complex conformation. The following binding mechanism was inferred. First, bosentan fluctuating randomly in solution is captured using a tip region of the flexible N-terminal tail of hETB via nonspecific attractive interactions (fly casting). Bosentan then slides occasionally from the tip to the root of the N-terminal tail (ligand–sliding). During this sliding, bosentan passes the gate of the binding pocket from outside to inside of the pocket with an accompanying rapid reduction of the molecular orientational variety of bosentan (orientational selection). Last, in the pocket, ligand–receptor attractive native contacts are formed. Eventually, the native-like complex is completed. The bosentan-captured conformations by the tip-region and root-region of the N-terminal tail correspond to two basins in the free-energy landscape. The ligand-sliding corresponds to overcoming of a free-energy barrier between the basins.

Nonstructural protein 1 (nsp1) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a 180-residue protein that blocks translation of host mRNAs in SARS-CoV-2-infected cells. Although it is known that SARS-CoV-2’s own RNA evades nsp1’s host translation shutoff, the molecular mechanism underlying the evasion was poorly understood. We performed an extended ensemble molecular dynamics simulation to investigate the mechanism of the viral RNA evasion. Simulation results showed that the stem loop structure of the SARS-CoV-2 RNA 5’-untranslated region (SL1) is recognized by both nsp1’s globular region and intrinsically disordered region. The recognition presumably enables selective translation of viral RNAs. Cluster analysis of the binding mode and detailed analysis of the binding poses revealed several residues involved in the SL1 recognition mechanism. The simulation results imply that the nsp1 C-terminal helices are lifted from the 40S ribosome upon the binding of SL1 to nsp1, unblocking translation of the viral RNA.

A preceding experiment suggested that a compound, which inhibits binding of the REST/NRSF segment to the cleft of a receptor protein mSin3B, can be a potential drug candidate to ameliorate many neuropathies. We have recently developed an enhanced conformational sampling method, genetic-algorithm-guided multi-dimensional virtual-system-coupled canonical molecular dynamics, and in the present study, applied it to three systems consisting of mSin3B and one of three compounds, sertraline, YN3, and acitretin. Other preceding experiments showed that only sertraline inhibits the binding of REST/NRSF to mSin3B. The current simulation study produced the spatial distribution of the compounds around mSin3B, and showed that sertraline and YN3 bound to the cleft of mSin3B with a high propensity, although acitretin did not. Further analyses of the simulation data indicated that only the sertraline-mSin3B complex produced a hydrophobic core similar to that observed in the molecular interface of the REST/NRSF-mSin3B complex: An aromatic ring of sertraline sunk deeply in the mSin3B's cleft forming a hydrophobic core contacting to hydrophobic amino-acid residues located at the bottom of the cleft. The present study proposes a step to design a compound that inhibits competitively the binding of a ligand to its receptor.

Protein-protein interactions between transmembrane helices are essential elements for membrane protein structures and functions. To understand the effects of peptide sequences and lipid compositions on these interactions, single-molecule experiments using model systems comprising artificial peptides and membranes have been extensively performed. However, their dynamic behavior at the atomic level remains largely unclear. In this study, we applied the all-atom molecular dynamics (MD) method to simulate the interactions of single-transmembrane helical peptide dimers in membrane environments, which has previously been analyzed by single-molecule experiments. The simulations were performed with two peptides (Ala- and Leu-based artificially designed peptides, termed "host peptide", and the host peptide added with the GXXXG motif, termed "GXXXG peptide"), two membranes (pure-POPC and POPC mixed with 30% cholesterols), and two dimer directions (parallel and antiparallel), consistent with those in the previous experiment. As a result, the MD simulations with parallel dimers reproduced the experimentally observed tendency that introducing cholesterols weakened the interactions in the GXXXG dimer and facilitated those in the host dimer. Our simulation suggested that the host dimer formed hydrogen bonds but the GXXXG dimer did not. However, some discrepancies were also observed between the experiments and simulations. Limitations in the space and time scales of simulations restrict the large-scale undulation and peristaltic motions of the membranes, resulting in differences in lateral pressure profiles. This effect could cause a discrepancy in the rotation angles of helices against the membrane normal.

An efficient algorithm to find the binding position and mode of small ligands bound at an active site of protein is proposed based on the spatial distribution function (SDF) obtained from the three-dimensional reference interaction site model (3D-RISM) theory with the Kovalenko-Hirata (KH) closure relation. The ligand examined includes hydrophobic, acidic, and basic molecules and zwitterions. Eighteen different types of proteins, which serve as targets for those ligands, are selected to examine the robustness of the algorithm. An imaginary atom, referred to as an "anchor site", is defined at the center of geometry of a ligand molecule that serves as a center for searching the binding position and mode of the ligand molecule in the translational and rotational spaces. The probable binding sites (PBSs) are identified based on the SDFs of the ligand molecules around the protein, and the PBS is ranked according to the peak height of SDF. The deviations from the mean height of the peak values of SDFs for 50 PBSs are analyzed based on the z-score, which is a measure of prominence of the site. The PBS found at the closest distance from the anchor site of the crystal structure is referred to as the "nearest site". The orientation of the ligand molecule at each PBS is explored by changing the Euler angles, and the most probable binding mode is determined based on the superposition approximation. The binding position of ligand molecules is successfully predicted as one of the distinct peaks in SDF of the anchor site, with a few exceptions. The binding mode of the ligand molecule predicted based on the superposition approximation is consistent with the X-ray crystal structure in nine systems, a half of the systems investigated. The significance of the results is discussed in detail. An application of the new protocol to fragment-based drug discovery is suggested.

The molecular dynamics (MD) method is a promising approach for investigating the molecular mechanisms of microscopic phenomena. In particular, generalized ensemble MD methods can efficiently explore the conformational space with a rugged free-energy surface. However, the implementation and acquisition of technical knowledge for each generalized ensemble MD method are not straightforward for end-users. Here, we present a new version of the myPresto/omegagene software, which is an MD simulation engine tailored for a series of generalized ensemble methods, which are virtual-system coupled multicanonical MD (V-McMD), virtual-system coupled adaptive umbrella sampling (V-AUS), and virtual-system coupled canonical MD (VcMD). This program has been applied in several studies analyzing free-energy landscapes of a variety of molecular systems with all-atom simulations. The updated version provides new functionality for coarse-grained simulations powered by the hydrophobicity scale method. The software package includes a step-by-step tutorial document for enhanced conformational sampling of the poly-glutamine (poly-Q) oligomer expressed as a one-bead per residue model. The myPresto/omegagene software is freely available at the following URL: https://github.com/kotakasahara/omegagene under the Apache2 license.

The molecular dynamics (MD) method is a promising technique to dissect the atomistic details of water dynamics around a solute. However, the quantitative predictions of experimentally measured kinetic properties, e.g. translational diffusion coefficient (D) and rotational relaxation time (?), are not straightforward. Current water models have failed to reproduce these properties quantitatively; therefore, the fine-tuning of water models is required. In this study, we examined the effects of ion?water Lennard-Jones (LJ) potentials on the water dynamics around a monovalent atomic ion. For the TIP5P water model, we introduced new LJ potentials for the ion?hydrogen and ion?pseudoatom interactions, which were zero in the original TIP5P model. The hydration properties, i.e. D, ?, and the radius of the first hydration shell (r(MO)), were examined for various parameter settings. As a result, the new parameters certainly improved the reproducibility of the hydration properties in correspondence with experimental values. However, it is still difficult to reproduce faster rotational relaxation of hydration water than that of bulk water. In addition, we found that the three hydration properties (D, ?, and r(MO)) were artificially correlated in the MD simulations.

ABSTRACTEnhanced conformational sampling, a genetic-algorithm-guided multi-dimensional virtual-system coupled molecular dynamics, can provide equilibrated conformational distributions of a receptor protein and a flexible ligand at room temperature. The distributions provide not only the most stable but also semi-stable complex structures, and propose a ligand–receptor binding process. This method was applied to a system consisting of a receptor protein, 14-3-3ε, and a flexible peptide, phosphorylated Myeloid leukemia factor 1 (pMLF1). The results present comprehensive binding pathways of pMLF1 to 14-3-3ε. We identified four thermodynamically stable clusters of MLF1 on the 14-3-3ε surface, and free-energy barriers among some clusters. The most stable cluster includes two high-density spots connected by a narrow corridor. When pMLF1 passes the corridor, a salt-bridge relay (switching) related to the phosphorylated residue of pMLF1 occurs. Conformations in one high-density spots are similar to the experimentally determined complex structure. Three-dimensional distributions of residues in the intermolecular interface rationally explain the binding-constant changes resultant from alanine–mutation experiment for the residues. We performed a simulation of non-phosphorylated peptide and 14-3-3ε, which demonstrated that the complex structure was unstable, suggesting that phosphorylation of the peptide is crucially important for binding to 14-3-3ε.

We introduced a conformational sampling method in an earlier report: The multi-dimensional virtual-system coupled molecular dynamics (mD-VcMD) enhances conformational sampling of a biomolecular system by computer simulations. Herein, new sampling method, a subzone-based mD-VcMD, is presented as an extension of mD-VcMD. Then, the subzone-based method is extended further using a genetic algorithm (GA) named the GA-guided mD-VcMD. In these methods, iterative simulation runs are performed to increase the sampled region gradually. The new methods have the following benefits: (1) They are free from a production run: i.e., all snapshots from all iterations are useful for analyses. (2) They are free from fine tuning of a weight function (probability distribution function or potential of mean force). (3) A canonical ensemble (i.e., a thermally equilibrated ensemble) is generated from a simple procedure. A thermodynamic weight is assigned to each snapshot. (4) Selective sampling can be performed for particularly addressing a poorly sampled region without breaking the proportion of the canonical ensemble if the poorly sampled conformational region emerges in sampling. By applying the methods to a simple system that involves an energy barrier between potential-energy minima, we demonstrated that the new methods have considerably higher sampling efficiency than the original mD-VcMD does.

Intrinsically disordered regions (IDRs) of a protein employ a flexible binding manner when recognizing a partner molecule. Moreover, it is recognized that binding of IDRs to a partner molecule is accompanied by folding, with a variety of bound conformations often being allowed in formation of the complex. In this study, we investigated a fragment of the disordered p53 C-terminal domain (CTDf) that interacts with one of its partner molecules, S100B, as a representative IDR. Although the 3D structure of CTDf in complex with S100B has been previously reported, the specific interactions remained controversial. To clarify these interactions, we performed generalized ensemble molecular dynamics (MD) simulations (virtual-system coupled multicanonical MD, termed V-McMD), which enable effective conformational sampling beyond that provided by conventional MD. These simulations generated a multimodal structural distribution for our system including CTDf and S100B, indicating that CTDf forms a variety of complex structures upon binding to S100B. We confirmed that our results are consistent with chemical shift perturbations and nuclear Overhauser effects that were observed in previous studies. Furthermore, we calculated the conformational entropy of CTDf in bound and isolated (free) states. Comparison of these CTDf entropies indicated that the disordered CTDf shows further increase in conformational diversity upon binding to S100B. Such entropy gain by binding may comprise an important feature of complex formation for IDRs.