Whole-Cell Models and Simulations in Molecular Detail.

Comprehensive data about the composition and structure of cellular components have enabled the construction of quantitative whole-cell models. While kinetic network-type models have been established, it is also becoming possible to build physical, molecular-level models of cellular environments. This review outlines challenges in constructing and simulating such models and discusses near- and long-term opportunities for developing physical whole-cell models that can connect molecular structure with biological function.

Molecular dynamics in multidimensional space explains how mutations affect the association path of neomycin to a riboswitch.

Riboswitches are messenger RNA (mRNA) fragments binding specific small molecules to regulate gene expression. A synthetic N1 riboswitch, inserted into yeast mRNA controls the translation of a reporter gene in response to neomycin. However, its regulatory activity is sensitive to single-point RNA mutations, even those distant from the neomycin binding site. While the association paths of neomycin to N1 and its variants remain unknown, recent fluorescence kinetic experiments indicate a two-step process driven by conformational selection. This raises the question of which step is affected by mutations. To address this, we performed all-atom two-dimensional replica-exchange molecular dynamics simulations for N1 and U14C, U14C[Formula: see text], U15A, and A17G mutants, ensuring extensive conformational sampling of both RNA and neomycin. The obtained neomycin association and binding paths, along with multidimensional free-energy profiles, revealed a two-step binding mechanism, consisting of conformational selection and induced fit. Neomycin binds to a preformed N1 conformation upon identifying a stable upper stem and U-turn motif in the riboswitch hairpin. However, the positioning of neomycin in the binding site occurs at different RNA-neomycin distances for each mutant, which may explain their different regulatory activities. The subsequent induced fit arises from the interactions of the neomycin's N3 amino group with RNA, causing the G9 backbone to rearrange. In the A17G mutant, the critical C6-A17/G17 stacking forms at a closer RNA-neomycin distance compared to N1. These findings together with estimated binding free energies coincide with experiments and elucidate why the A17G mutation decreases and U15A enhances N1 activity in response to neomycin.

Formation of extramembrane ß-strands controls dimerization of transmembrane helices in amyloid precursor protein C99.

The 99-residue C-terminal domain of amyloid precursor protein (APP-C99), precursor to amyloid beta (Aß), is a transmembrane (TM) protein containing intrinsically disordered N- and C-terminal extramembrane domains. Using molecular dynamics (MD) simulations, we show that the structural ensemble of the C99 monomer is best described in terms of thousands of states. The C99 monomer has a propensity to form ß-strand in the C-terminal extramembrane domain, which explains the slow spin relaxation times observed in paramagnetic probe NMR experiments. Surprisingly, homodimerization of C99 not only narrows the conformational ensemble from thousands to a few states through the formation of metastable ß-strands in extramembrane domains but also stabilizes extramembrane α-helices. The extramembrane domain structure is observed to dramatically impact the homodimerization motif, resulting in the modification of TM domain conformations. Our study provides an atomic-level structural basis for communication between the extramembrane domains of the C99 protein and TM homodimer formation. This finding could serve as a general model for understanding the influence of disordered extramembrane domains on TM protein structure.

Deciphering the Multi-state Conformational Equilibrium of HDM2 in the Regulation of p53 Binding: Perspectives from Molecular Dynamics Simulation and NMR Analysis.

HDM2 negatively regulates the activity of the tumor suppressor p53. Previous NMR studies have shown that apo-HDM2 interconverts between an "open" state in which the N-terminal "lid" is disordered and a "closed" state in which the lid covers the p53-binding site in the core region. Molecular dynamics (MD) simulation studies have been performed to elucidate the conformational dynamics of HDM2, but the direct relevance of the experimental and computational analyses is unclear. In addition, how the phosphorylation of S17 in the lid contributes to the inhibition of p53 binding remains controversial. Here, we used both NMR and MD simulations to investigate the conformational dynamics of apo-HDM2. The NMR analysis revealed that apo-HDM2 exists in a fast-exchanging equilibrium within two closed states, closed 1 and closed 2, in addition to a previously demonstrated slow-exchanging "open-closed" equilibrium. MD simulations visualized two characteristic closed states, where the spatial orientation of the key residues corresponds well to the chemical shift changes of the NMR spectra. Furthermore, the phosphorylation of S17 induced an equilibrium shift toward closed 1, thereby suppressing the binding of p53 to HDM2. This study reveals a multi-state equilibrium of apo-HDM2 and provides new insights into the regulation mechanism of HDM2-p53 interactions.

SPANA: Spatial decomposition analysis for cellular-scale molecular dynamics simulations.

The rapid increase in computational power with the latest supercomputers has enabled atomistic molecular dynamics (MDs) simulations of biomolecules in biological membrane, cytoplasm, and other cellular environments. These environments often contain a million or more atoms to be simulated simultaneously. Therefore, their trajectory analyses involve heavy computations that can become a bottleneck in the computational studies. Spatial decomposition analysis (SPANA) is a set of analysis tools in the Generalized-Ensemble Simulation System (GENESIS) software package that can carry out MD trajectory analyses of large-scale biological simulations using multiple CPU cores in parallel. SPANA applies the spatial decomposition of a large biological system to distribute structural and dynamical analyses into individual CPU cores, which reduces the computational time and the memory size, significantly. SPANA opens new possibilities for detailed atomistic analyses of biomacromolecules as well as solvent water molecules, ions, and metabolites in MD simulation trajectories of very large biological systems containing more than millions of atoms in cellular environments.

Micelle-like clusters in phase-separated Nanog condensates: A molecular simulation study.

The phase separation model for transcription suggests that transcription factors (TFs), coactivators, and RNA polymerases form biomolecular condensates around active gene loci and regulate transcription. However, the structural details of condensates remain elusive. In this study, for Nanog, a master TF in mammalian embryonic stem cells known to form protein condensates in vitro, we examined protein structures in the condensates using residue-level coarse-grained molecular simulations. Human Nanog formed micelle-like clusters in the condensate. In the micelle-like cluster, the C-terminal disordered domains, including the tryptophan repeat (WR) regions, interacted with each other near the cluster center primarily via hydrophobic interaction. In contrast, hydrophilic disordered N-terminal and DNA-binding domains were exposed on the surface of the clusters. Electrostatic attractions of these surface residues were responsible for bridging multiple micelle-like structures in the condensate. The micelle-like structure and condensate were dynamic and liquid-like. Mutation of tryptophan residues in the WR region which was implicated to be important for a Nanog function resulted in dissolution of the Nanog condensate. Finally, to examine the impact of Nanog cluster to DNA, we added DNA fragments to the Nanog condensate. Nanog DNA-binding domains exposed to the surface of the micelle-like cluster could recruit more than one DNA fragments, making DNA-DNA distance shorter.

Similarity scores of vibrational spectra reveal the atomistic structure of pentapeptides in multiple basins.

Vibrational spectroscopy combined with theoretical calculations is a powerful tool for analyzing the interaction and conformation of peptides at the atomistic level. Nonetheless, identifying the structure becomes increasingly difficult as the peptide size grows large. One example is acetyl-SIVSF-N-methylamide, a capped pentapeptide, whose atomistic structure has remained unknown since its first observation [T. Sekiguchi, M. Tamura, H. Oba, P. Çarçarbal, R. R. Lozada-Garcia, A. Zehnacker-Rentien, G. Grégoire, S. Ishiuchi and M. Fujii, Angew. Chem., Int. Ed., 2018, 57, 5626-5629]. Here, we propose a novel conformational search method, which exploits the structure-spectrum correlation using a similarity score that measures the agreement of theoretical and experimental spectra. Surprisingly, the two conformers have distinctly different energy and geometry. The second conformer is 25 kJ mol-1 higher in energy than the other, lowest-energy conformer. The result implies that there are multiple pathways in the early stage of the folding process: one to the global minimum and the other to a different basin. Once such a structure is established, the second conformer is unlikely to overcome the barrier to produce the most stable structure due to a vastly different hydrogen bond network of the backbone. Our proposed method can characterize the lowest-energy conformer and kinetically trapped, high-energy conformers of complex biomolecules.

Learning QM/MM potential using equivariant multiscale model.

The machine learning (ML) method emerges as an efficient and precise surrogate model for high-level electronic structure theory. Its application has been limited to closed chemical systems without considering external potentials from the surrounding environment. To address this limitation and incorporate the influence of external potentials, polarization effects, and long-range interactions between a chemical system and its environment, the first two terms of the Taylor expansion of an electrostatic operator have been used as extra input to the existing ML model to represent the electrostatic environments. However, high-order electrostatic interaction is often essential to account for external potentials from the environment. The existing models based only on invariant features cannot capture significant distribution patterns of the external potentials. Here, we propose a novel ML model that includes high-order terms of the Taylor expansion of an electrostatic operator and uses an equivariant model, which can generate a high-order tensor covariant with rotations as a base model. Therefore, we can use the multipole-expansion equation to derive a useful representation by accounting for polarization and intermolecular interaction. Moreover, to deal with long-range interactions, we follow the same strategy adopted to derive long-range interactions between a target system and its environment media. Our model achieves higher prediction accuracy and transferability among various environment media with these modifications.

Structural and energetic analysis of metastable intermediate states in the E1P-E2P transition of Ca2+-ATPase.

Sarcoplasmic reticulum (SR) Ca2+-ATPase transports two Ca2+ ions from the cytoplasm to the SR lumen against a large concentration gradient. X-ray crystallography has revealed the atomic structures of the protein before and after the dissociation of Ca2+, while biochemical studies have suggested the existence of intermediate states in the transition between E1P⋅ADP⋅2Ca2+ and E2P. Here, we explore the pathway and free energy profile of the transition using atomistic molecular dynamics simulations with the mean-force string method and umbrella sampling. The simulations suggest that a series of structural changes accompany the ordered dissociation of ADP, the A-domain rotation, and the rearrangement of the transmembrane (TM) helices. The luminal gate then opens to release Ca2+ ions toward the SR lumen. Intermediate structures on the pathway are stabilized by transient sidechain interactions between the A- and P-domains. Lipid molecules between TM helices play a key role in the stabilization. Free energy profiles of the transition assuming different protonation states suggest rapid exchanges between Ca2+ ions and protons when the Ca2+ ions are released toward the SR lumen.

Structural effects of spike protein D614G mutation in SARS-CoV-2.

A single mutation from aspartate to glycine at position 614 has dominated all circulating variants of the severe acute respiratory syndrome coronavirus 2. D614G mutation induces structural changes in the spike (S) protein that strengthen the virus infectivity. Here, we use molecular dynamics simulations to dissect the effects of mutation and 630-loop rigidification on S-protein structure. The introduction of the mutation orders the 630-loop structure and thereby induces global structural changes toward the cryoelectron microscopy structure of the D614G S-protein. The ordered 630-loop weakens local interactions between the 614th residue and others in contrast to disordered structures in the wild-type protein. The mutation allosterically alters global interactions between receptor-binding domains, forming an asymmetric and mobile down conformation and facilitating transitions toward up conformation. The loss of salt bridge between D614 and K854 upon the mutation generally stabilizes S-protein protomer, including the fusion peptide proximal region that mediates membrane fusion. Understanding the molecular basis of D614G mutation is crucial as it dominates in all variants of concern, including Delta and Omicron.

Acceleration of generalized replica exchange with solute tempering simulations of large biological systems on massively parallel supercomputer.

Generalized replica exchange with solute tempering (gREST) is one of the enhanced sampling algorithms for proteins or other systems with rugged energy landscapes. Unlike the replica-exchange molecular dynamics (REMD) method, solvent temperatures are the same in all replicas, while solute temperatures are different and are exchanged frequently between replicas for exploring various solute structures. Here, we apply the gREST scheme to large biological systems containing over one million atoms using a large number of processors in a supercomputer. First, communication time on a multi-dimensional torus network is reduced by matching each replica to MPI processors optimally. This is applicable not only to gREST but also to other multi-copy algorithms. Second, energy evaluations, which are necessary for the multistate bennet acceptance ratio (MBAR) method for free energy estimations, are performed on-the-fly during the gREST simulations. Using these two advanced schemes, we observed 57.72 ns/day performance in 128-replica gREST calculations with 1.5 million atoms system using 16,384 nodes in Fugaku. These schemes implemented in the latest version of GENESIS software could open new possibilities to answer unresolved questions on large biomolecular complex systems with slow conformational dynamics.

Development of hidden Markov modeling method for molecular orientations and structure estimation from high-speed atomic force microscopy time-series images.

High-speed atomic force microscopy (HS-AFM) is a powerful technique for capturing the time-resolved behavior of biomolecules. However, structural information in HS-AFM images is limited to the surface geometry of a sample molecule. Inferring latent three-dimensional structures from the surface geometry is thus important for getting more insights into conformational dynamics of a target biomolecule. Existing methods for estimating the structures are based on the rigid-body fitting of candidate structures to each frame of HS-AFM images. Here, we extend the existing frame-by-frame rigid-body fitting analysis to multiple frames to exploit orientational correlations of a sample molecule between adjacent frames in HS-AFM data due to the interaction with the stage. In the method, we treat HS-AFM data as time-series data, and they are analyzed with the hidden Markov modeling. Using simulated HS-AFM images of the taste receptor type 1 as a test case, the proposed method shows a more robust estimation of molecular orientations than the frame-by-frame analysis. The method is applicable in integrative modeling of conformational dynamics using HS-AFM data.

Implementation of residue-level coarse-grained models in GENESIS for large-scale molecular dynamics simulations.

Residue-level coarse-grained (CG) models have become one of the most popular tools in biomolecular simulations in the trade-off between modeling accuracy and computational efficiency. To investigate large-scale biological phenomena in molecular dynamics (MD) simulations with CG models, unified treatments of proteins and nucleic acids, as well as efficient parallel computations, are indispensable. In the GENESIS MD software, we implement several residue-level CG models, covering structure-based and context-based potentials for both well-folded biomolecules and intrinsically disordered regions. An amino acid residue in protein is represented as a single CG particle centered at the Cα atom position, while a nucleotide in RNA or DNA is modeled with three beads. Then, a single CG particle represents around ten heavy atoms in both proteins and nucleic acids. The input data in CG MD simulations are treated as GROMACS-style input files generated from a newly developed toolbox, GENESIS-CG-tool. To optimize the performance in CG MD simulations, we utilize multiple neighbor lists, each of which is attached to a different nonbonded interaction potential in the cell-linked list method. We found that random number generations for Gaussian distributions in the Langevin thermostat are one of the bottlenecks in CG MD simulations. Therefore, we parallelize the computations with message-passing-interface (MPI) to improve the performance on PC clusters or supercomputers. We simulate Herpes simplex virus (HSV) type 2 B-capsid and chromatin models containing more than 1,000 nucleosomes in GENESIS as examples of large-scale biomolecular simulations with residue-level CG models. This framework extends accessible spatial and temporal scales by multi-scale simulations to study biologically relevant phenomena, such as genome-scale chromatin folding or phase-separated membrane-less condensations.

Towards de novo design of transmembrane α-helical assemblies using structural modelling and molecular dynamics simulation.

Computational de novo protein design involves iterative processes consisting of amino acid sequence design, structural modelling and scoring, and design validation by synthesis and experimental characterisation. Recent advances in protein structure prediction and modelling methods have enabled the highly efficient and accurate design of water-soluble proteins. However, the design of membrane proteins remains a major challenge. To advance membrane protein design, considering the higher complexity of membrane protein folding, stability, and dynamic interactions between water, ions, lipids, and proteins is an important task. For introducing explicit solvents and membranes to these design methods, all-atom molecular dynamics (MD) simulations of designed proteins provide useful information that cannot be obtained experimentally. In this review, we first describe two major approaches to designing transmembrane α-helical assemblies, consensus and de novo design. We further illustrate recent MD studies of membrane protein folding related to protein design, as well as advanced treatments in molecular models and conformational sampling techniques in the simulations. Finally, we discuss the possibility to introduce MD simulations after the existing static modelling and screening of design decoys as an additional step for refinement of the design, which considers membrane protein folding dynamics and interactions with explicit membranes.

Allosteric regulation of ß-reaction stage I in tryptophan synthase upon the α-ligand binding.

Tryptophan synthase (TRPS) is a bifunctional enzyme consisting of α- and ß-subunits that catalyzes the last two steps of L-tryptophan (L-Trp) biosynthesis. The first stage of the reaction at the ß-subunit is called ß-reaction stage I, which converts the ß-ligand from an internal aldimine [E(Ain)] to an α-aminoacrylate [E(A-A)] intermediate. The activity is known to increase 3-10-fold upon the binding of 3-indole-D-glycerol-3'-phosphate (IGP) at the α-subunit. The effect of α-ligand binding on ß-reaction stage I at the distal ß-active site is not well understood despite the abundant structural information available for TRPS. Here, we investigate the ß-reaction stage I by carrying out minimum-energy pathway searches based on a hybrid quantum mechanics/molecular mechanics (QM/MM) model. The free-energy differences along the pathway are also examined using QM/MM umbrella sampling simulations with QM calculations at the B3LYP-D3/aug-cc-pVDZ level of theory. Our simulations suggest that the sidechain orientation of ßD305 near the ß-ligand likely plays an essential role in the allosteric regulation: a hydrogen bond is formed between ßD305 and the ß-ligand in the absence of the α-ligand, prohibiting a smooth rotation of the hydroxyl group in the quinonoid intermediate, whereas the dihedral angle rotates smoothly after the hydrogen bond is switched from ßD305-ß-ligand to ßD305-ßR141. This switch could occur upon the IGP-binding at the α-subunit, as evidenced by the existing TRPS crystal structures.

Extension of the iSoLF implicit-solvent coarse-grained model for multicomponent lipid bilayers.

iSoLF is a coarse-grained (CG) model for lipid molecules with the implicit-solvent approximation used in molecular dynamics (MD) simulations of biological membranes. Using the original iSoLF (iSoLFv1), MD simulations of lipid bilayers consisting of either POPC or DPPC and these bilayers, including membrane proteins, can be performed. Here, we improve the original model, explicitly treating the electrostatic interactions between different lipid molecules and adding CG particle types. As a result, the available lipid types increase to 30. To parameterize the potential functions of the new model, we performed all-atom MD simulations of each lipid at three different temperatures using the CHARMM36 force field and the modified TIP3P model. Then, we parameterized both the bonded and non-bonded interactions to fit the area per lipid and the membrane thickness of each lipid bilayer by using the multistate Boltzmann Inversion method. The final model reproduces the area per lipid and the membrane thickness of each lipid bilayer at the three temperatures. We also examined the applicability of the new model, iSoLFv2, to simulate the phase behaviors of mixtures of DOPC and DPPC at different concentrations. The simulation results with iSoLFv2 are consistent with those using Dry Martini and Martini 3, although iSoLFv2 requires much fewer computations. iSoLFv2 has been implemented in the GENESIS MD software and is publicly available.

Short disordered protein segment regulates cross-species transmission of a yeast prion.

Soluble prion proteins contingently encounter foreign prion aggregates, leading to cross-species prion transmission. However, how its efficiency is regulated by structural fluctuation of the host soluble prion protein remains unsolved. In the present study, through the use of two distantly related yeast prion Sup35 proteins, we found that a specific conformation of a short disordered segment governs interspecies prion transmissibility. Using a multidisciplinary approach including high-resolution NMR and molecular dynamics simulation, we identified critical residues within this segment that allow interspecies prion transmission in vitro and in vivo, by locally altering dynamics and conformation of soluble prion proteins. Remarkably, subtle conformational differences caused by a methylene group between asparagine and glutamine sufficed to change the short segment structure and substantially modulate the cross-seeding activity. Thus, our findings uncover how conformational dynamics of the short segment in the host prion protein impacts cross-species prion transmission. More broadly, our study provides mechanistic insights into cross-seeding between heterologous proteins.

Modified Hamiltonian in FEP Calculations for Reducing the Computational Cost of Electrostatic Interactions.

The free-energy perturbation (FEP) method predicts relative and absolute free-energy changes of biomolecules in solvation and binding with other molecules. FEP is, therefore, one of the most essential tools in in silico drug design. In conventional FEP, to smoothly connect two thermodynamic states, the potential energy is modified as a linear combination of the end-state potential energies by introducing scaling factors. When the particle mesh Ewald is used for electrostatic calculations, conventional FEP requires two reciprocal-space calculations per time step, which largely decreases the computational performance. To overcome this problem, we propose a new FEP scheme by introducing a modified Hamiltonian instead of interpolation of the end-state potential energies. The scheme introduces nonuniform scaling into the electrostatic potential as used in Replica Exchange with Solute Tempering 2 (REST2) and does not require additional reciprocal-space calculations. We tested this modified Hamiltonian in FEP calculations in several biomolecular systems. In all cases, the calculated free-energy changes with the current scheme are in good agreement with those from conventional FEP. The modified Hamiltonian in FEP greatly improves the computational performance, which is particularly marked for large biomolecular systems whose reciprocal-space calculations are the major bottleneck of total computational time.

Protein folding intermediates on the dimensionality reduced landscape with UMAP and native contact likelihood.

To understand protein folding mechanisms from molecular dynamics (MD) simulations, it is important to explore not only folded/unfolded states but also representative intermediate structures on the conformational landscape. Here, we propose a novel approach to construct the landscape using the uniform manifold approximation and projection (UMAP) method, which reduces the dimensionality without losing data-point proximity. In the approach, native contact likelihood is used as feature variables rather than the conventional Cartesian coordinates or dihedral angles of protein structures. We tested the performance of UMAP for coarse-grained MD simulation trajectories of B1 domain in protein G and observed on-pathway transient structures and other metastable states on the UMAP conformational landscape. In contrast, these structures were not clearly distinguished on the dimensionality reduced landscape using principal component analysis or time-lagged independent component analysis. This approach is also useful to obtain dynamical information through Markov state modeling and would be applicable to large-scale conformational changes in many other biomacromolecules.

Clustering and dynamics of crowded proteins near membranes and their influence on membrane bending.

Atomistic molecular dynamics simulations of concentrated protein solutions in the presence of a phospholipid bilayer are presented to gain insights into the dynamics and interactions at the cytosol-membrane interface. The main finding is that proteins that are not known to specifically interact with membranes are preferentially excluded from the membrane, leaving a depletion zone near the membrane surface. As a consequence, effective protein concentrations increase, leading to increased protein contacts and clustering, whereas protein diffusion becomes faster near the membrane for proteins that do occasionally enter the depletion zone. Since protein-membrane contacts are infrequent and short-lived in this study, the structure of the lipid bilayer remains largely unaffected by the crowded protein solution, but when proteins do contact lipid head groups, small but statistically significant local membrane curvature is induced, on average.