An exa-scale high-performance molecular dynamics simulation program: MODYLAS.

A new version of the highly parallelized general-purpose molecular dynamics (MD) simulation program MODYLAS with high performance on the Fugaku computer was developed. A benchmark test using Fugaku indicated highly efficient communication, single instruction, multiple data (SIMD) processing, and on-cache arithmetic operations. The system's performance deteriorated only slightly, even under high parallelization. In particular, a newly developed minimum transferred data method, requiring a significantly lower amount of data transfer compared to conventional communications, showed significantly high performance. The coordinates and forces of 101 810 176 atoms and the multipole coefficients of the subcells could be distributed to the 32 768 nodes (1 572 864 cores) in 2.3 ms during one MD step calculation. The SIMD effective instruction rates for floating-point arithmetic operations in direct force and fast multipole method (FMM) calculations measured on Fugaku were 78.7% and 31.5%, respectively. The development of a data reuse algorithm enhanced the on-cache processing; the cache miss rate for direct force and FMM calculations was only 2.74% and 1.43%, respectively, on the L1 cache and 0.08% and 0.60%, respectively, on the L2 cache. The modified MODYLAS could complete one MD single time-step calculation within 8.5 ms for the aforementioned large system. Additionally, the program contains numerous functions for material research that enable free energy calculations, along with the generation of various ensembles and molecular constraints.

Prediction of self-diffusion coefficients of chemically diverse pure liquids by all-atom molecular dynamics simulations.

Molecular self-diffusion coefficients underlie various kinetic properties of the liquids involved in chemistry, physics, and pharmaceutics. In this study, 547 self-diffusion coefficients are calculated based on all-atom molecular dynamics (MD) simulations of 152 diverse pure liquids at various temperatures employing the OPLS4 force field. The calculated coefficients are compared with experimental data (424 extracted from the literature and 123 newly measured by pulsed-field gradient nuclear magnetic resonance). The calculations well agree with the experimental values. The determination coefficient and root mean square error between the observed and calculated logarithmic self-diffusion coefficients of the 547 entries are 0.931 and 0.213, respectively, demonstrating that the MD calculation can be an excellent industrial tool for predicting, for example, molecular transportation in liquids such as the diffusion of active ingredients in biological and pharmaceutical liquids. The self-diffusion coefficients collected in this study are compiled into a database for broad researches including artificial intelligence calculations.

Global diffusion of hydrogen molecules in the heterogeneous structure of polymer electrolytes for fuel cells: Dynamic Monte Carlo combined with molecular dynamics calculations.

Using our recently developed dynamic Monte Carlo (MC) method [Nagai et al., J. Chem. Phys. 156, 154506 (2022)], we investigated the global diffusion of hydrogen molecules over structural heterogeneities of polymer electrolyte membranes in fuel cells. The three-dimensional position-dependent free energies and the diffusion constants of the hydrogen molecules, required by the present dynamic MC calculations, were taken from our previous study [Nagai et al., J. Chem. Phys. 156, 044507 (2022)] and newly evaluated in this work, respectively. The calculations enabled evaluating the hydrogen dynamics over long-time scales, including global diffusion constants. Based on the calculated global diffusion constants and free energies, the permeability of hydrogen molecules was estimated via the solubility-diffusion model. The estimated values were in good agreement with the reported experimental data, thus validating the present methodology. The analysis of the Monte Carlo trajectories indicated that the main permeation paths are located in the polymer and interfacial phases, although the water phase may make a non-negligible contribution to mass transport.

Dynamic Monte Carlo calculation generating particle trajectories that satisfy the diffusion equation for heterogeneous systems with a position-dependent diffusion coefficient and free energy.

A series of new Monte Carlo (MC) transition probabilities was investigated that could produce molecular trajectories statistically satisfying the diffusion equation with a position-dependent diffusion coefficient and potential energy. The MC trajectories were compared with the numerical solution of the diffusion equation by calculating the time evolution of the probability distribution and the mean first passage time, which exhibited excellent agreement. The method is powerful when investigating, for example, the long-distance and long-time global transportation of a molecule in heterogeneous systems by coarse-graining them into one-particle diffusive molecular motion with a position-dependent diffusion coefficient and free energy. The method can also be applied to many-particle dynamics.

Three-dimensional free-energy landscape of hydrogen and oxygen molecules in polymer electrolyte membranes: Insight into diffusion paths.

Polymer electrolyte membranes, for example, the Nafion™ membranes, used in the fuel cells are responsible for separating reactive gas molecules as well as for the efficient exchange of protons. Although control of the permeation of the gases is important to enhance the fuel cell performance, the mechanism by which hydrogen and oxygen molecules permeate through the membranes remains unclear. To clarify the mechanism, we investigated the three-dimensional free-energy landscape of hydrogen and oxygen molecules in Nafion membranes with various water contents focusing on relevant diffusion paths. Low-free-energy paths are found mainly in the polymer phase and the interfacial region between the polymer and water phases. Thus, the path of the transportation may be attributed to the polymer phase and interfacial phases. However, the free-energy value in the aqueous phase is only slightly higher (∼1-2 kBT) than that in the other two phases, which indicates that a secondary contribution from the aqueous phase is expected. The free-energy landscape in the polymer and interfacial phases was found rugged, while it is comparatively flat in the water phase. We also found that an increase in water content brings about a smoother free-energy landscape in the polymer and interfacial phases. The decreased ruggedness may facilitate the gas diffusivity. These observations help understand the molecular mechanism of the gas diffusion in the membranes.

Momentum removal to obtain the position-dependent diffusion constant in constrained molecular dynamics simulation.

The position-dependent diffusion coefficient along with free energy profile are important parameters needed to study mass transport in heterogeneous systems such as biological and polymer membranes, and molecular dynamics (MD) calculation is a popular tool to obtain them. Among many methodologies, the Marrink-Berendsen (MB) method is often employed to calculate the position-dependent diffusion coefficient, in which the autocorrelation function of the force on a fixed molecule is related to the friction on the molecule. However, the diffusion coefficient is shown to be affected by the period of the removal of the center-of-mass velocity, τv0 , which is necessary when performing MD calculations using the Ewald method for Coulombic interaction. We have clarified theoretically in this study how this operation affects the diffusion coefficient calculated by the MB method, and the theoretical predictions are proven by MD calculations. Therefore, we succeeded in providing guidance on how to select an appropriate τv0 value in estimating the position-dependent diffusion coefficient by the MB method. This guideline is applicable also to the Woolf-Roux method.

Cryo-Cooling Effect on DHFR Crystal Studied by Replica-Exchange Molecular Dynamics Simulations.

Cryo-cooling is routinely performed before x-ray diffraction image collection to reduce the damage to crystals due to ionizing radiation. It has been suggested that although backbone structures are usually very similar between room temperature and cryo-temperature, cryo-cooling may hamper biologically relevant dynamics. In this study, the crystal of Escherichia coli dihydrofolate reductase is studied with replica-exchange molecular dynamics simulation, and the results are compared with the crystal structure determined at cryo-temperature and room temperature with the time-averaged ensemble method. Although temperature dependence of unit cell compaction and root mean-square fluctuation of Cα is found in accord with experiment, it is found that the protein structure at low temperature can be more heterogeneous than the ensemble of structures reported by using the time-averaged ensemble method, encouraging further development of the time-averaged ensemble method and indicating that data should be examined carefully to avoid overinterpretation of one average structure.

Exploring the impact of proteins on the line tension of a phase-separating ternary lipid mixture.

The separation of lipid mixtures into thermodynamically stable phase-separated domains is dependent on lipid composition, temperature, and system size. Using molecular dynamics simulations, the line tension between thermodynamically stable lipid domains formed from ternary mixtures of di-C16:0 PC:di-C18:2 PC:cholesterol at 40:40:20 mol. % ratio was investigated via two theoretical approaches. The line tension was found to be 3.1 ± 0.2 pN by capillary wave theory and 4.7 ± 3.7 pN by pressure tensor anisotropy approaches for coarse-grained models based on the Martini force field. Using an all-atom model of the lipid membrane based on the CHARMM36 force field, the line tension was found to be 3.6 ± 0.9 pN using capillary wave theory and 1.8 ± 2.2 pN using pressure anisotropy approaches. The discrepancy between estimates of the line tension based on capillary wave theory and pressure tensor anisotropy methods is discussed. Inclusion of protein in Martini membrane lipid mixtures was found to reduce the line tension by 25%-35% as calculated by the capillary wave theory approach. To further understand and predict the behavior of proteins in phase-separated membranes, we have formulated an analytical Flory-Huggins model and parameterized it against the simulation results. Taken together these results suggest a general role for proteins in reducing the thermodynamic cost associated with domain formation in lipid mixtures and quantifies the thermodynamic driving force promoting the association of proteins to domain interfaces.

Gaussian mixture model for coarse-grained modeling from XFEL.

We explore the advantage of Gaussian mixture model (GMM) for interpretation of single particle diffraction patterns from X-ray free electron laser (XFEL) experiments. GMM approximates a biomolecular shape by the superposition of Gaussian distributions. As the Fourier transformation of GMM can be quickly performed, we can efficiently simulate XFEL diffraction patterns from approximated structure models. We report that the resolution that GMM can accurately reproduce is proportional to the cubic root of the number of Gaussians used in the modeling. This behavior can be attributed to the correspondence between the number of adjustable parameters in GMM and the amount of sampling points in diffraction space. Furthermore, GMMs can successfully be used to perform angular assignment and to detect conformational variation. These results demonstrate that GMMs serve as useful coarse-grained models for hybrid approach in XFEL single particle experiments.

Role of Computational Methods in Going beyond X-ray Crystallography to Explore Protein Structure and Dynamics.

Protein structural biology came a long way since the determination of the first three-dimensional structure of myoglobin about six decades ago. Across this period, X-ray crystallography was the most important experimental method for gaining atomic-resolution insight into protein structures. However, as the role of dynamics gained importance in the function of proteins, the limitations of X-ray crystallography in not being able to capture dynamics came to the forefront. Computational methods proved to be immensely successful in understanding protein dynamics in solution, and they continue to improve in terms of both the scale and the types of systems that can be studied. In this review, we briefly discuss the limitations of X-ray crystallography in studying protein dynamics, and then provide an overview of different computational methods that are instrumental in understanding the dynamics of proteins and biomacromolecular complexes.

Critical size dependence of domain formation observed in coarse-grained simulations of bilayers composed of ternary lipid mixtures.

Model cellular membranes are known to form micro- and macroscale lipid domains dependent on molecular composition. The formation of macroscopic lipid domains by lipid mixtures has been the subject of many simulation investigations. We present a critical study of system size impact on lipid domain phase separation into liquid-ordered and liquid-disordered macroscale domains in ternary lipid mixtures. In the popular di-C16:0 PC:di-C18:2 PC:cholesterol at 35:35:30 ratio mixture, we find systems with a minimum of 1480 lipids to be necessary for the formation of macroscopic phase separated domains and systems of 10 000 lipids to achieve structurally converged conformations similar to the thermodynamic limit. To understand these results and predict the behavior of any mixture forming two phases, we develop and investigate an analytical Flory-Huggins model which is recursively validated using simulation and experimental data. We find that micro- and macroscale domains can coexist in ternary mixtures. Additionally, we analyze the distributions of specific lipid-lipid interactions in each phase, characterizing domain structures proposed based on past experimental studies. These findings offer guidance in selecting appropriate system sizes for the study of phase separations and provide new insights into the nature of domain structure for a popular ternary lipid mixture.

On the use of mass scaling for stable and efficient simulated tempering with molecular dynamics.

Simulated tempering (ST) is a generalized-ensemble algorithm that employs trajectories exploring a range of temperatures to effectively sample rugged energy landscapes. When implemented using the molecular dynamics method, ST can require the use of short time steps for ensuring the stability of trajectories at high temperatures. To address this shortcoming, a mass-scaling ST (MSST) method is presented in which the particle mass is scaled in proportion to the temperature. Mass scaling in the MSST method leads to velocity distributions that are independent of temperature and eliminates the need for velocity scaling after the accepted temperature updates that are required in conventional ST simulations. The homogeneity in time scales with changing temperature improves the stability of simulations and allows for the use of longer time steps at high temperatures. As a result, the MSST is found to be more efficient than the standard ST method, particularly for cases in which a large temperature range is employed. © 2016 Wiley Periodicals, Inc.

Mass-scaling replica-exchange molecular dynamics optimizes computational resources with simpler algorithm.

We develop a novel method of replica-exchange molecular dynamics (REMD) simulation, mass-scaling REMD (MSREMD) method, which improves numerical stability of simulations. In addition, the MSREMD method can also simplify a replica-exchange routine by eliminating velocity scaling. As a pilot system, a Lennard-Jones fluid is simulated with the new method. The results suggest that the MSREMD method improves the numerical stability at high temperatures compared with the conventional REMD method. For the Nosé-Hoover thermostats, we analytically demonstrate that the MSREMD simulations can reproduce completely the same trajectories of the conventional REMD ones with shorter time steps at high temperatures. Accordingly, we can easily compare the computational costs of the REMD and MSREMD simulations. We conclude that the MSREMD method decreases the instability and optimizes the computational resources with simpler algorithm.

A hybrid approach to study large conformational transitions of biomolecules from single particle XFEL diffraction data.

X-ray free-electron laser (XFEL) is the latest generation of the X-ray source that could become an invaluable technique in structural biology. XFEL has ultrashort pulse duration, extreme peak brilliance, and high spatial coherence, which could enable the observation of the biological molecules in near nature state at room temperature without crystallization. However, for biological systems, due to their low diffraction power and complexity of sample delivery, experiments and data analysis are not straightforward, making it extremely challenging to reconstruct three-dimensional (3D) structures from single particle XFEL data. Given the current limitations to the amount and resolution of the data from such XFEL experiments, we propose a new hybrid approach for characterizing biomolecular conformational transitions by using a single 2D low-resolution XFEL diffraction pattern in combination with another known conformation. In our method, we represent the molecular structure with a coarse-grained model, the Gaussian mixture model, to describe large conformational transitions from low-resolution XFEL data. We obtain plausible 3D structural models that are consistent with the XFEL diffraction pattern by deforming an initial structural model to maximize the similarity between the target pattern and the simulated diffraction patterns from the candidate models. We tested the proposed algorithm on two biomolecules of different sizes with different complexities of conformational transitions, adenylate kinase, and elongation factor 2, using synthetic XFEL data. The results show that, with the proposed algorithm, we can successfully describe the conformational transitions by flexibly fitting the coarse-grained model of one conformation to become consistent with an XFEL diffraction pattern simulated from another conformation. In addition, we showed that the incident beam orientation has some effect on the accuracy of the 3D structure modeling and discussed the reasons for the inaccuracies for certain orientations. The proposed method could serve as an alternative approach for retrieving information on 3D conformational transitions from the XFEL diffraction patterns to interpret experimental data. Since the molecules are represented by Gaussian kernels and no atomic structure is needed in principle, such a method could also be used as a tool to seek initial models for 3D reconstruction algorithms.

Synchronization of circadian oscillation of phosphorylation level of KaiC in vitro.

In recent experimental reports, robust circadian oscillation of the phosphorylation level of KaiC has been reconstituted by incubating three cyanobacterial proteins, KaiA, KaiB, and KaiC, with ATP in vitro. This reconstitution indicates that protein-protein interactions and the associated ATP hydrolysis suffice to generate the oscillation, and suggests that the rhythm arising from this protein-based system is the circadian clock pacemaker in cyanobacteria. The mechanism of this reconstituted oscillation, however, remains elusive. In this study, we extend our previous model of oscillation by explicitly taking two phosphorylation sites of KaiC into account and we apply the extended model to the problem of synchrony of two oscillatory samples mixed at different phases. The agreement between the simulated and observed data suggests that the combined mechanism of the allosteric transition of KaiC hexamers and the monomer shuffling between them plays a key role in synchronization among KaiC hexamers and hence underlies the population-level oscillation of the ensemble of Kai proteins. The predicted synchronization patterns in mixtures of unequal amounts of two samples provide further opportunities to experimentally check the validity of the proposed mechanism. This mechanism of synchronization should be important in vivo for the persistent oscillation when Kai proteins are synthesized at random timing in cyanobacterial cells.

Position-Dependent Diffusion Constant of Molecules in Heterogeneous Systems as Evaluated by the Local Mean Squared Displacement.

The authors propose a novel method to evaluate the position-dependent diffusion constant by analyzing unperturbed segments of a trajectory determined by the additional flat-bottom potential. The accuracy of this novel method is first established by studying homogeneous systems, where the reference value can be obtained by the Einstein relation. The applicability of this new method to heterogeneous systems is then demonstrated by studying a hydrophobic solute near a hydrophobic wall. The proposed method is also comprehensively compared with popular conventional methods, whereby the significance of the present method is illustrated. The novel method is powerful and useful for studying kinetics in heterogeneous systems based on molecular dynamics calculations.

Simulated tempering and magnetizing: application of two-dimensional simulated tempering to the two-dimensional Ising model and its crossover.

We have performed two-dimensional simulated tempering (ST) simulations of the two-dimensional Ising model with different lattice sizes in order to investigate the two-dimensional ST's applicability to dealing with phase transitions and study the crossover of critical scaling behavior. The external field, as well as the temperature, was treated as a dynamical variable updated during the simulations. Thus this simulation can be referred to as simulated tempering and magnetizing (STM). We also performed simulated magnetizing (SM) simulations, in which the external field was considered as a dynamical variable and temperature was not. As discussed in previous studies, the ST method is not always compatible with first-order phase transitions. This is also true in the magnetizing process. Flipping of the entire magnetization did not occur in the SM simulations under the critical temperature T{c} in large-lattice-size simulations; however, the phase changed through the high-temperature region in the STM simulations. Thus the dimensional extension let us eliminate the difficulty of the first-order phase transitions and study a wide area of the phase space. We discuss how frequently parameter-updating attempts should be made for optimal convergence. The results favor frequent attempts. We finally study the crossover behavior of the phase transitions with respect to the temperature and external field. The crossover behavior is clearly observed in the simulations, in agreement with the theoretical implications.