Optimal estimates of self-diffusion coefficients from molecular dynamics simulations

Computational and Machine Learning Methods for CO2 Capture Using Metal-Organic Frameworks

Machine learning (ML) using data sets of atomic and molecular force fields (FFs) has made significant progress and provided benefits in the fields of chemistry and material science. This work examines the interactions between chemistry and materials computational science at the atomic and molecular scales for metal-organic framework (MOF) adsorbent development toward carbon dioxide (CO2) capture. Herein, a connection will be drawn between atomic forces predicted by ML algorithms and the structures of MOFs for CO2 adsorption. Our study also takes into account the successes of atomic computational screening in the field of materials science, especially quantum ML, and its relationship to ML algorithms that clarify advancements in the area of CO2 adsorption by MOFs. Additionally, we reviewed the processes for supplying data to ML algorithms for algorithm training, including text mining from scientific articles, and MOF's formula processing linked to the chemical properties of MOFs. To create ML algorithms for future research, we recommend that the digitization of scientific records can help efficiently synthesize advanced MOFs. Finally, a future vision for developing pioneer MOF synthesis routes for CO2 capture is presented in this review article.

From Latent Dynamics to Meaningful Representations

While representation learning has been central to the rise of machine learning and artificial intelligence, a key problem remains in making the learned representations meaningful. For this, the typical approach is to regularize the learned representation through prior probability distributions. However, such priors are usually unavailable or are ad hoc. To deal with this, recent efforts have shifted toward leveraging the insights from physical principles to guide the learning process. In this spirit, we propose a purely dynamics-constrained representation learning framework. Instead of relying on predefined probabilities, we restrict the latent representation to follow overdamped Langevin dynamics with a learnable transition density─a prior driven by statistical mechanics. We show that this is a more natural constraint for representation learning in stochastic dynamical systems, with the crucial ability to uniquely identify the ground truth representation. We validate our framework for different systems including a real-world fluorescent DNA movie data set. We show that our algorithm can uniquely identify orthogonal, isometric, and meaningful latent representations.

Vitrification-enabled enhancement of proton conductivity in hydrogen-bonded organic frameworks

Hydrogen-bonded organic frameworks (HOFs) are versatile materials with potential applications in proton conduction. Traditional approaches involve incorporating humidity control to address grain boundary challenges for proton conduction. This study finds vitrification as an alternative strategy to eliminate grain boundary effect in HOFs by rapidly melt quenching the kinetically stable HOF-SXU-8 to glassy state HOF-g. Notably, a remarkable enhancement in proton conductivity without humidity was achieved after vitrification, from 1.31 × 10-7 S cm-1 to 5.62× 10-2 S cm-1 at 100 °C. Long term stability test showed negligible performance degradation, and even at 30 °C, the proton conductivity remained at high level of 1.2 × 10-2 S cm-1. Molecule dynamics (MD) simulations and X-ray total scattering experiments reveal the HOF-g system is consisted of three kinds of clusters, i.e., 1,5-Naphthalenedisulfonic acid (1,5-NSA) anion clusters, N,N-dimethylformamide (DMF) molecule clusters, and H+-H2O clusters. In which, the H+ plays an important role to bridge these clusters and the high conductivity is mainly related to the H+ on H3O+. These findings provide valuable insights for optimizing HOFs, enabling efficient proton conduction, and advancing energy conversion and storage devices.

Experimentally Probing the Effect of Confinement Geometry on Lipid Diffusion

The lateral mobility of molecules within the cell membrane is ultimately governed by the local environment of the membrane. Confined regions induced by membrane structures, such as protein aggregates or the actin meshwork, occur over a wide range of length scales and can impede or steer the diffusion of membrane components. However, a detailed picture of the origins and nature of these confinement effects remains elusive. Here, we prepare model lipid systems on substrates patterned with confined domains of varying geometries constructed with different materials to explore the influences of physical boundary conditions and specific molecular interactions on diffusion. We demonstrate a platform that is capable of significantly altering and steering the long-range diffusion of lipids by using simple oxide deposition approaches, enabling us to systematically explore how confinement size and shape impact diffusion over multiple length scales. While we find that a "boundary condition" description of the system captures underlying trends in some cases, we are also able to directly compare our systems to analytical models, revealing the unexpected breakdown of several approximate solutions. Our results highlight the importance of considering the length scale dependence when discussing properties such as diffusion.

On the Statistical Mechanics of Mass Accommodation at Liquid-Vapor Interfaces

We propose a framework for describing the dynamics associated with the adsorption of small molecules to liquid-vapor interfaces using an intermediate resolution between traditional continuum theories that are bereft of molecular detail and molecular dynamics simulations that are replete with them. In particular, we develop an effective single particle equation of motion capable of describing the physical processes that determine thermal and mass accommodation probabilities. The effective equation is parametrized with quantities that vary through space away from the liquid-vapor interface. Of particular importance in describing the early time dynamics is the spatially dependent friction, for which we propose a numerical scheme to evaluate from molecular simulation. Taken together with potentials of mean force computable with importance sampling methods, we illustrate how to compute the mass accommodation coefficient and residence time distribution. Throughout, we highlight the case of ozone adsorption in aqueous solutions and its dependence on electrolyte composition.

Determination of self-diffusion coefficients in mixtures with benchtop 13C NMR spectroscopy via polarization transfer

Nuclear magnetic resonance (NMR) is an established method to determine self-diffusion coefficients in liquids with high precision. The development of benchtop NMR spectrometers makes the method accessible to a wider community. In most cases, 1H NMR spectroscopy is used to determine self-diffusion coefficients due to its high sensitivity. However, especially when using benchtop NMR spectrometers for the investigation of complex mixtures, the signals in 1H NMR spectra can overlap, hindering the precise determination of self-diffusion coefficients. In 13C NMR spectroscopy, the signals of different compounds are generally well resolved. However, the sensitivity of 13C NMR is significantly lower than that of 1H NMR spectroscopy leading to very long measurement times, which makes diffusion coefficient measurements based on 13C NMR practically infeasible with benchtop NMR spectrometers. To circumvent this problem, we have combined two known pulse sequences, one for polarization transfer from 1H to the 13C nuclei (PENDANT) and one for the measurement of diffusion coefficients (PFG). The new method (PENPFG) was used to measure the self-diffusion coefficients of three pure solvents (acetonitrile, ethanol and 1-propanol) as well as in all their binary mixtures and the ternary mixture at various compositions. For comparison, also measurements of the same systems were carried out with a standard PFG-NMR routine on a high-field NMR instrument. The results are in good agreement and show that PENPFG is a useful tool for the measurement of the absolute value of the self-diffusion coefficients in complex liquid mixtures with benchtop NMR spectrometers.

On-chip light sheet illumination for nanoparticle tracking in microfluidic channels

A simple and inexpensive method is presented to efficiently integrate light sheet illumination in a microfluidic chip for dark-field microscopic tracking and sizing of nanoparticles. The basic idea is to insert an optical fiber inside a polydimethylsiloxane (PDMS) elastomer microfluidic chip and use it as a cylindrical lens. The optical fiber is in this case no longer seen as only an optical waveguide but as a ready-made micro-optical component that is inexpensive and easy to source. Upon insertion, the optical fiber stretches the PDMS microchannel walls, which has two effects. The first effect is to tone down the intrinsic ripples in the PDMS that would otherwise create inhomogeneities in the light sheet illumination. The second effect is to remove any obliqueness of the channel wall and constrain it to be strictly perpendicular to the propagation of the illumination, avoiding the formation of a prismatic diopter. Through calculations, numerical simulations and measurements, we show that the optimal configuration consists in creating a slowly converging light sheet so that its axial thickness is almost uniform along the tracked area. The corresponding thickness was estimated at 12 μm, or 10 times the depth of field of the optical system. This leads to an at least six-fold increase in the signal-to-noise ratio compared to the case without the cylindrical lens. This original light-sheet configuration is used to track and size spherical gold nanoparticles with diameters of 80 nm and 50 nm.

Temperature-Induced Restructuring of Mycolic Acid Bilayers Modeling the Mycobacterium tuberculosis Outer Membrane: A Molecular Dynamics Study

The emergence of new drug-resistant strains of the tuberculosis pathogen Mycobacterium tuberculosis (Mtb) is a new challenge for modern medicine. Its resistance capacity is closely related to the properties of the outer membrane of the Mtb cell wall, which is a bilayer membrane formed by mycolic acids (MAs) and their derivatives. To date, the molecular mechanisms of the response of the Mtb outer membrane to external factors and, in particular, elevated temperatures have not been sufficiently studied. In this work, we consider the temperature-induced changes in the structure, ordering, and molecular mobility of bilayer MA membranes of various chemical and conformational compositions. Using all-atom long-term molecular dynamics simulations of various MA membranes, we report the kinetic parameters of temperature-dependent changes in the MA self-diffusion coefficients and conformational compositions, including the apparent activation energies of these processes, as well as the characteristic times of ordering changes and the features of phase transitions occurring over a wide range of elevated temperatures. Understanding these effects could be useful for the prevention of drug resistance and the development of membrane-targeting pharmaceuticals, as well as in the design of membrane-based materials.

Efficient single-run implementation of generalized Einstein relation to compute transport coefficients: A binary-based time sampling

A robust and simple implementation of the generalized Einstein formulation using single equilibrium molecular dynamics simulation is introduced to compute diffusion and shear viscosity. The unique features underlying this framework are as follows: (1) The use of a simple binary-based method to sample time-dependent transport coefficients results in a uniform distribution of data on a logarithmic time scale. Although we sample "on-the-fly," the algorithm is readily applicable for post-processing analysis. Overlapping same-length segments are not sampled as they indicate strong correlations. (2) Transport coefficients are estimated using a power law fitting function, a generalization of the standard linear relation, that accurately describes the long-time plateau. (3) The use of a generalized least squares (GLS) fitting estimator to explicitly consider correlations between fitted data points results in a reliable estimate of the statistical uncertainties in a single run. (4) The covariance matrix for the GLS method is estimated analytically using the Wiener process statistics and computed variances. (5) We provide a Python script to perform the fits and automate the procedure to determine the optimal fitting domain. The framework is applied to two fluids, binary hard sphere and a Lennard-Jones near the triple point, and the validity of the single-run estimates is verified against multiple independent runs. The approach should be applicable to other transport coefficients since the diffusive limit is universal to all of them. Given its rigor and simplicity, this methodology can be readily incorporated into standard molecular dynamics packages using on-the-fly or post-processing analysis.

Trajectory Statistical Learning of the Potential Mean of Force and Diffusion Coefficient from Molecular Dynamics Simulations

Central to studying the conformational changes of a complex protein is understanding the dynamics and energetics involved. Phenomenologically, structural dynamics can be formulated using an overdamped Langevin model along an observable, e.g., the distance between two residues in the protein. The Langevin model is specified by the deterministic force (the potential of mean force, PMF) and stochastic force (characterized by the diffusion coefficient, D). It is therefore of great interest to be able to extract both PMF and D from an observable time series but under the same computational framework. Here, we approach this challenge in molecular dynamics (MD) simulations by treating it as a missing-data Bayesian estimation problem. An important distinction in our methodology is that the entire MD trajectory, as opposed to the individual data elements, is used as the statistical variable in Bayesian imputation. This idea is implemented through an eigen-decomposition procedure for a time-symmetrized Fokker-Planck equation, followed by maximizing the likelihood for parameter estimation. The mathematical expressions for the functional derivatives used in learning PMF and D also provide new physical insights for the manner by which the information on both the deterministic and stochastic forces is encoded in the dynamics data. An all-atom MD simulation of a nontrivial biomolecule case is used to illustrate the application of this approach. We show that, interestingly, the results of trajectory statistical learning can motivate new order parameters for an improved description of the kinetic bottlenecks in conformational changes. Complementing purely data-driven or black-box methods, this work underscores the advantages of physics-based machine learning in gaining chemical insights from quantitative parameter estimation.

Understanding the Mechanical Properties of Ultradeformable Liposomes Using Molecular Dynamics Simulations

Improving drug delivery efficiency to solid tumor sites is a central challenge in anticancer therapeutic research. Our previous experimental study (Guo et al., Nat. Commun. 2018, 9, 130) showed that soft, elastic liposomes had increased uptake and accumulation in cancer cells and tumors in vitro and in vivo respectively, relative to rigid particles. As a first step toward understanding how liposomes' molecular structure and composition modulates their elasticity, we performed all-atom and coarse-grained classical molecular dynamics (MD) simulations of lipid bilayers formed by mixing a long-tailed unsaturated phospholipid with a short-tailed saturated lipid with the same headgroup. The former types of phospholipids considered were 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) and 1,2-dipalmitoleoyl-sn-glycero-3-phosphocholine (termed here DPMPC). The shorter saturated lipids examined were 1,2-diheptanoyl-sn-glycero-3-phosphocholine (DHPC), 1,2-didecanoyl-sn-glycero-3-phosphocholine (DDPC), 1,2-dilauroyl-sn-glycero-3-phosphocholine (DLPC), and 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC). Several lipid concentrations and surface tensions were considered. Our results show that DOPC or DPMPC systems having 25-35 mol % of the shortest lipids DHPC or DDPC are the least rigid, having area compressibility moduli KA that are ∼10% smaller than the values observed in pure DOPC or DPMPC bilayers. These results agree with experimental measurements of the stretching modulus and lysis tension in liposomes with the same compositions. These mixed systems also have lower areas per lipid and form more uneven x-y interfaces with water, the tails of both primary and secondary lipids are more disordered, and the terminal methyl groups in the tails of the long lipid DOPC or DPMPC wriggle more in the vertical direction, compared to pure DOPC or DPMPC bilayers or their mixtures with the longer saturated lipid DLPC or DMPC. These observations confirm our hypothesis that adding increasing concentrations of the short unsaturated lipid DHPC or DDPC to DOPC or DPMPC bilayers alters lipid packing and thus makes the resulting liposomes more elastic and less rigid. No formation of lipid nanodomains was noted in our simulations, and no clear trends were observed in the lateral diffusivities of the lipids as the concentration, type of secondary lipid, and surface tension were varied.

Elucidating anticancer drugs release from UiO-66 as a carrier through the computational approaches

The computational analysis of drug release from metal-organic frameworks (MOFs), specifically UiO-66, is the primary focus of this research. MOFs are recognized as nanocarriers due to their crystalline structure, porosity, and potential for added functionalities. The research examines the release patterns of three drugs: temozolomide, alendronate, and 5-fluorouracil, assessing various factors such as the drugs' distance from the UiO-66 centers, the interaction of drug functional groups with Zr metal ions, and the drug density throughout the nanocarrier. Findings reveal that 5-fluorouracil is located furthest from the UiO-66 center and exhibits the highest positive energy compared to the other drugs. Alendronate's density is observed to shift to the carrier surface, while 5-fluorouracil's density significantly decreases within the system. The drug density diminishes as the distance from the UiO-66 center of mass increases, suggesting a stronger positive interaction between the drugs and the nanocarrier. Moreover, Monte Carlo calculations were employed to load drugs onto the UiO-66 surface, leading to a substantial release of 5-fluorouracil from UiO-66. Quantum and Monte Carlo adsorption localization calculations were also conducted to gather data on the compounds' energy and geometry. This research underscores the potential of MOFs as nanocarriers for drug delivery and highlights the crucial role of temperature in regulating drug release from UiO-66. It provides insights into the complex dynamics of drug release and the factors influencing it, thereby emphasizing the promise of UiO-66 as a viable candidate for drug delivery. This work contributes to our understanding of UiO-66's role and sets the stage for improved performance optimization in the cancer treatment.

Understanding the Anomalous Diffusion of Water in Aqueous Electrolytes Using Machine Learned Potentials

The diffusivity of water in aqueous cesium iodide solutions is larger than that in neat liquid water and vice versa for sodium chloride solutions. Such peculiar ion-specific behavior, called anomalous diffusion, is not reproduced in typical force field based molecular dynamics (MD) simulations due to inadequate treatment of ion-water interactions. Herein, this hurdle is tackled by using machine learned atomic potentials (MLPs) trained on data from density functional theory calculations. MLP based atomistic MD simulations of aqueous salt solutions reproduce experimentally determined thermodynamic, structural, dynamical, and transport properties, including their varied trends in water diffusivities across salt concentration. This enables an examination of their intermolecular structure to unravel the microscopic underpinnings of the differences in their transport properties. While both ions in CsI solutions contribute to the faster diffusion of water molecules, the competition between the heavy retardation by Na ions and the slight acceleration by Cl ions in NaCl solutions reduces their water diffusivity.

Unwrapping NPT Simulations to Calculate Diffusion Coefficients

In molecular dynamics simulations in the NPT ensemble at constant pressure, the size and shape of the periodic simulation box fluctuate with time. For particle images far from the origin, the rescaling of the box by the barostat results in unbounded position displacements. Special care is thus required when a particle trajectory is unwrapped from a projection into the central box under periodic boundary conditions to a trajectory in full three-dimensional space, e.g., for the calculation of translational diffusion coefficients. Here, we review and compare different schemes in use for trajectory unwrapping. We also specify the corresponding rewrapping schemes to put an unwrapped trajectory back into the central box. On this basis, we then identify a scheme for the calculation of diffusion coefficients from NPT simulations, which is a primary application of trajectory unwrapping. In this scheme, the wrapped and unwrapped trajectory are mutually consistent and their statistical properties are preserved. We conclude with advice on best practice for the consistent unwrapping of constant-pressure simulation trajectories and the calculation of accurate translational diffusion coefficients.

H-Bonds in Crambin: Coherence in an α-Helix

We applied coherence analysis-used by engineers to identify linear interactions in stochastic systems-to molecular dynamics simulations of crambin, a thionin storage protein found in Abyssinian cabbage. A key advantage of coherence over other analyses is that it is robust, independent of the properties, or even the existence of probability distributions often relied on in statistical mechanics. For frequencies between 0.391 and 5.08 GHz (corresponding reciprocally to times of 2.56 and 0.197 ns), the displacements of oxygen and nitrogen atoms across α-helix H-bonds are strongly correlated, with a coherence greater than 0.9; the secondary structure causes the H-bonds to effectively act as a spring. Similar coherence behavior is observed for covalent bonds and other noncovalent interactions including H-bonds in β-sheets and salt bridges. In contrast, arbitrary pairs of atoms that are physically distant have uncorrelated motions and negligible coherence. These results suggest that coherence may be used to objectively identify atomic interactions without subjective thresholds such as H-bond lengths angles and angles. Strong coherence is also observed between the average position of adjacent leaves (groups of atoms) in an α-helix, suggesting that the harmonic analysis of classical molecular dynamics can successfully describe the propagation of allosteric interactions through the structure.

Real Space and Time Imaging of Collective Headgroup Dipole Motions in Zwitterionic Lipid Bilayers

Lipid bilayers are supramolecular structures responsible for a range of processes, such as transmembrane transport of ions and solutes, and sorting and replication of genetic materials, to name just a few. Some of these processes are transient and currently, cannot be visualized in real space and time. Here, we developed an approach using 1D, 2D, and 3D Van Hove correlation functions to image collective headgroup dipole motions in zwitterionic phospholipid bilayers. We show that both 2D and 3D spatiotemporal images of headgroup dipoles are consistent with commonly understood dynamic features of fluids. However, analysis of the 1D Van Hove function reveals lateral transient and re-emergent collective dynamics of the headgroup dipoles-occurring at picosecond time scales-that transmit and dissipate heat at longer times, due to relaxation processes. At the same time, the headgroup dipoles also generate membrane surface undulations due a collective tilting of the headgroup dipoles. A continuous intensity band of headgroup dipole spatiotemporal correlations-at nanometer length and nanosecond time scales-indicates that dipoles undergo stretching and squeezing elastic deformations. Importantly, the above mentioned intrinsic headgroup dipole motions can be externally stimulated at GHz-frequency scale, enhancing their flexoelectric and piezoelectric capabilities (i.e., increased conversion efficiency of mechanical energy into electric energy). In conclusion, we discuss how lipid membranes can provide molecular-level insights about biological learning and memory, and as platforms for the development of the next generation of neuromorphic computers.

Sequence-to-Sequence Change-Point Detection in Single-Particle Trajectories via Recurrent Neural Network for Measuring Self-Diffusion

A recurrent neural network is developed for segmenting between anomalous and normal diffusion in single-particle trajectories. Accurate segmentation infers a distinct change point that is used to approximate an Einstein linear regime in the mean-squared displacement curve via the transition density function, a unique physical descriptor for short-lived and delayed transiency. Through several artificial and simulated scenarios, we demonstrate the compelling accuracy of our model for dissecting linear and nonlinear behaviour. The inherent practicality of our model lies in its ability to substantiate the self-diffusion coefficient through offline trajectory segmentation, which is opposed to the common ‘best-guess’ linear fitting standard. Additionally, we show that the transition density function has fundamental implications and correspondence to underlying mechanisms that influence transition. In particular, we show that the known proportionality between salt concentration and diffusion of water also influences delayed anomalous behaviour.

Optimal Bond Constraint Topology for Molecular Dynamics Simulations of Cholesterol

We recently observed artificial temperature gradients in molecular dynamics (MD) simulations of phase-separating ternary lipid mixtures using the Martini 2 force field. We traced this artifact to insufficiently converged bond length constraints with typical time steps and default settings for the linear constraint solver (LINCS). Here, we systematically optimize the constraint scaffold of cholesterol. With massive virtual sites in an equimomental arrangement, we accelerate bond constraint convergence while preserving the original cholesterol force field and dynamics. The optimized model does not induce nonphysical temperature gradients even at relaxed LINCS settings and is at least as fast as the original model at the strict LINCS settings required for proper thermal sampling. We provide a python script to diagnose possible problems with constraint convergence for other molecules and force fields. Equimomental constraint topology optimization can also be used to boost constraint convergence in atomistic MD simulations of molecular systems.

Diffusion of small molecule drugs is affected by surface interactions and crowder proteins

Crowded environments are known to affect the diffusion of macromolecules, but their effects on the diffusion of small molecules are largely uncharacterized. We investigate how three protein crowders, bovine serum albumin (BSA), hen egg-white lysozyme, and myoglobin, influence the diffusion rates and interactions of four small molecules: fluorescein, and three drugs, doxorubicin, glycogen synthase kinase-3 inhibitor SB216763, and quinacrine. Using Line-FRAP measurements, Brownian dynamics simulations, and molecular docking, we find that the diffusion rates of the small molecules are highly affected by self-aggregation, interactions with the proteins, and surface adsorption. The diffusion of fluorescein is decreased because of its interactions with the protein crowders and their surface adsorption. Protein crowders increase the diffusion rates of doxorubicin and SB216763 by reducing surface interactions and self-aggregation, respectively. Quinacrine diffusion was not affected by protein crowders. The mechanistic insights gained here may assist in optimization of compounds for higher mobility in complex macromolecular environments.

A Multi-Scale Approach to Model K+ Permeation Through the KcsA Channel

K+ channels allow a very efficient passage of K+ ions through the membrane while excluding Na+ ions, and these properties are essential for life. The 3D structure of the KcsA K+ channel, solved more than 20 years ago, allows to address many relevant aspects of K+ permeation and selectivity mechanisms at the molecular level. Recent crystallographic data and molecular dynamics (MD) studies suggest that no water is normally present inside the selectivity filter (SF), which can instead accommodate four adjacent K+ ions. Using a multi-scale approach, whereby information taken from a low-level simulation approach is used to feed a high-level model, we studied the mechanism of K+ permeation through KcsA channels. More specifically, we used MD to find stable ion configurations under physiological conditions. They were characterized by two adjacent K+ ions occupying the more central positions of the SF (sites S2 and S3), while the other two K+ ions could be found at the external and internal entrances to the SF. Sites S1 and S4 were instead not occupied by K+. A continuum Bikerman–Poisson–Boltzmann model that takes into account the volume of the ions and their dehydration when entering the SF fully confirmed the MD results, showing peaks of K+ occupancy at S2, S3, and the external and internal entrances, with S1 and S4 sites being virtually never occupied by K+. Inspired by the newly found ion configuration in the SF at equilibrium, we developed a simple kinetic permeation model which, fed with kinetic rate constants assessed from molecular meta-dynamics, reproduced the main permeation properties of the KcsA channel found experimentally, including sublinear current-voltage and saturating conductance-concentration relationships. This good agreement with the experimental data also implies that the ion configuration in the SF we identified at equilibrium would also be a key configuration during permeation.

Small ionic radii limit time step in Martini 3 molecular dynamics simulations

Among other improvements, the Martini~3 coarse-grained force field provides a more accurate description of the solvation of protein pockets and channels through the consistent use of various bead types and sizes. Here, we show that the representation of Na+ and Cl- ions as ``tiny' (TQ5) beads limits the accessible time step to 25~fs. By contrast, with Martini~2, time steps of 30-40~fs were possible for lipid bilayer systems without proteins. This limitation is relevant for systems that require long equilibration times. We derive a quantitative kinetic model of time-integration instabilities in molecular dynamics as a function of time step, ion concentration and mass, system size, and simulation time. We demonstrate that ion-water interactions are the main source of instability at physiological conditions, followed closely by ion-ion interactions. We show that increasing the ionic masses makes it possible to use time steps up to 40~fs with minimal impact on static equilibrium properties and on dynamical quantities such as lipid and ion diffusion coefficients. Increasing the size of the bead representing the ions (and thus changing their hydration) also permits longer time steps. For a soluble protein, we find that increasing the mass of tiny beads also on the protein permits simulations with 30-fs time step. The use of larger time steps in Martini~3 results in a more efficient exploration of configuration space. The kinetic model of MD simulation crashes can be used to determine the maximum allowed time step upfront for an efficient use of resources and whenever sampling efficiency is critical.

Systematic Evaluation of Ion Diffusion and Water Exchange

As a fundamental property of all fluids, diffusion plays myriad roles in both science and our daily lives. Diffusive properties of many liquids including water have been extensively studied both experimentally and theoretically, while for transition metal ions, there exist significant experimental data that have not been extensively studied theoretically. Hence, high-confidence predictions for challenging systems like radioactive ions that are biohazardous cannot be reliably generated. In this work, a workflow named ISAIAH (Ion Simulation using AMBER for dIffusion Action when Hydrated) was designed to accurately simulate the diffusion coefficients of 15 monoatomic ions with charges varying from -1 to +3 in four water models. As the results indicate, good agreement with experimental values was achieved, leading us to select ²³⁹Pu⁴⁺ (for which no experimental data are available) as a candidate ion to make a theoretical prediction of its diffusion coefficient in water. Among all the force field parameter sets, the ones parametrized using an augmented 12-6-4 Lennard-Jones (LJ) potential showed lower average unsigned errors (AUE) for ions of various radii and electron configurations relative to some 12-6 LJ parameters. This observation agrees well with the fact that diffusion is affected by both the hydration free energy (HFE) and the ion-oxygen distance (IOD) between solute and solvent molecules, both of which are handled well by the 12-6-4 model.

Bayesian estimation of the hydroxyl radical diffusion coefficient at low temperature and high pressure from atomistic molecular dynamics

The hydroxyl radical is the primary reactive oxygen species produced by the radiolysis of water and is a significant source of radiation damage to living organisms. Mobility of the hydroxyl radical at low temperatures and/or high pressures is hence a potentially important factor in determining the challenges facing psychrophilic and/or barophilic organisms in high-radiation environments (e.g., ice-interface or undersea environments in which radiative heating is a potential heat and energy source). Here, we estimate the diffusion coefficient for the hydroxyl radical in aqueous solution using a hierarchical Bayesian model based on atomistic molecular dynamics trajectories in TIP4P/2005 water over a range of temperatures and pressures.

Multiscale molecular kinetics by coupling Markov state models and reaction-diffusion dynamics

A novel approach to simulate simple protein-ligand systems at large time and length scales is to couple Markov state models (MSMs) of molecular kinetics with particle-based reaction-diffusion (RD) simulations, MSM/RD. Currently, MSM/RD lacks a mathematical framework to derive coupling schemes, is limited to isotropic ligands in a single conformational state, and lacks multiparticle extensions. In this work, we address these needs by developing a general MSM/RD framework by coarse-graining molecular dynamics into hybrid switching diffusion processes. Given enough data to parameterize the model, it is capable of modeling protein-protein interactions over large time and length scales, and it can be extended to handle multiple molecules. We derive the MSM/RD framework, and we implement and verify it for two protein-protein benchmark systems and one multiparticle implementation to model the formation of pentameric ring molecules. To enable reproducibility, we have published our code in the MSM/RD software package.

Multiscale Investigation of the Diffusion Mechanism within the Solid-Electrolyte Interface Layer: Coupling Quantum Mechanics, Molecular Dynamics, and Macroscale Mathematical Modeling

The solid-electrolyte interface (SEI) layer has a critical role in Li-ion batteries' (LIBs) life span. The SEI layer, even in modern commercial LIBs, is responsible for more than 50% of capacity loss. Due to the inherent complexity in studying the SEI layer, many aspects of its performance and characteristics, including diffusion mechanisms in this layer, are unknown. As a result, most mathematical models use a constant value of the diffusion coefficient, instead of a variable formulation, to predict LIBs' properties and performance such as capacity fading and the SEI growth rate. In this work, by employing a multiscale investigation using a combination of quantum mechanics, molecular dynamics, and macroscale mathematical modeling, some equations are presented to evaluate the energy barrier against diffusion and the diffusion coefficient in different crystal structures in the inner section of the SEI layer. The equations are evaluated as a function of temperature and concentration and can be used to study the diffusion mechanism in the SEI layer. They can also be integrated with other mathematical models of LIBs to increase the accuracy of the latter.

Maximum likelihood estimates of diffusion coefficients from single-particle tracking experiments

Single-molecule localization microscopy allows practitioners to locate and track labeled molecules in biological systems. When extracting diffusion coefficients from the resulting trajectories, it is common practice to perform a linear fit on mean-squared-displacement curves. However, this strategy is suboptimal and prone to errors. Recently, it was shown that the increments between the observed positions provide a good estimate for the diffusion coefficient, and their statistics are well-suited for likelihood-based analysis methods. Here, we revisit the problem of extracting diffusion coefficients from single-particle tracking experiments subject to static noise and dynamic motion blur using the principle of maximum likelihood. Taking advantage of an efficient real-space formulation, we extend the model to mixtures of subpopulations differing in their diffusion coefficients, which we estimate with the help of the expectation-maximization algorithm. This formulation naturally leads to a probabilistic assignment of trajectories to subpopulations. We employ the theory to analyze experimental tracking data that cannot be explained with a single diffusion coefficient. We test how well a dataset conforms to the assumptions of a diffusion model and determine the optimal number of subpopulations with the help of a quality factor of known analytical distribution. To facilitate use by practitioners, we provide a fast open-source implementation of the theory for the efficient analysis of multiple trajectories in arbitrary dimensions simultaneously.

Arginine multivalency stabilizes protein/RNA condensates

Biomolecular condensates assembled through liquid-liquid phase separation (LLPS) of proteins and RNAs are currently recognized to play an important role in cellular organization. Their assembly depends on the formation of a network of transient, multivalent interactions between flexible scaffold biomolecules. Understanding how protein and RNA sequences determine these interactions and ultimately regulate the phase separation is an open key challenge. Recent in vitro studies have revealed that arginine and lysine residues, which are enriched in most cellular condensates, have markedly distinct propensities to drive the LLPS of protein/RNA mixtures. Here, we employ explicit-solvent atomistic molecular dynamics simulations to shed light on the microscopic origin of this difference by investigating mixtures of polyU oligonucleotides with either polyR/polyK peptides. In agreement with experiments, our simulations indicate that arginine has a higher affinity for polyU than lysine both in highly diluted conditions and in concentrated solutions with a biomolecular density comparable to cellular condensate. The analysis of intermolecular contacts suggests that this differential behavior is due to the propensity of arginine side chains to simultaneously form a higher number of specific interactions with oligonucleotides, including hydrogen bonds and stacking interactions. Our results provide a molecular description of how the multivalency of the guanidinium group enables the coordination of multiple RNA groups by a single arginine residue, thus ultimately stabilizing protein/RNA condensates.

Molecular Details of Protein Condensates Probed by Microsecond Long Atomistic Simulations

The formation of membraneless organelles in cells commonly occurs via liquid-liquid phase separation (LLPS) and is in many cases driven by multivalent interactions between intrinsically disordered proteins (IDPs). Investigating the nature of these interactions, and their effect on dynamics within the condensed phase, is therefore of critical importance but very challenging for either simulation or experiment. Here, we study these interactions and their dynamics by pairing a novel multiscale simulation strategy with microsecond all-atom MD simulations of a condensed, IDP-rich phase. We simulate two IDPs this way, the low complexity domain of FUS and the N-terminal disordered domain of LAF-1, and find good agreement with experimental information about average density, water content, and residue-residue contacts. We go significantly beyond what is known from experiments by showing that ion partitioning within the condensed phase is largely driven by the charge distribution of the proteins and-in the cases considered-shows little evidence of preferential interactions of the ions with the proteins. Furthermore, we can probe the microscopic diffusive dynamics within the condensed phase, showing that water and ions are in dynamic equilibrium between dense and dilute phases, and their diffusion is reduced in the dense phase. Despite their high concentration in the condensate, the protein molecules also remain mobile, explaining the observed liquid-like properties of this phase. We finally show that IDP self-association is driven by a combination of nonspecific hydrophobic interactions as well as hydrogen bonds, salt bridges, and π-π and cation-π interactions. The simulation approach presented here allows the structural and dynamical properties of biomolecular condensates to be studied in microscopic detail and is generally applicable to single- and multicomponent systems of proteins and nucleic acids involved in LLPS.

Estimates of Electrical Conductivity from Molecular Dynamics Simulations: How to Invest the Computational Effort

Although the electrical conductivity of an electrolyte can be estimated from the molecular dynamics trajectory, it is often a challenging task because of the need to obtain a substantial amount of data to ensure sufficient averaging. Here, we present an analysis on the convergence of results with the number of simulated trajectories. A series of molecular dynamics simulations have been performed for a model electrolyte (NaCl in water) and the Einstein relation has been used to calculate the electrical conductivity. The standard deviation of the conductivity estimates is relatively large compared to the mean value, and it has been shown that the off-diagonal contributions to the collective displacement of ions are responsible for large deviations between systems. It has been found that about 40 independent MD simulations may be required to reduce the errors. A procedure to improve the final estimate of the conductivity has been proposed.

Systematic errors in diffusion coefficients from long-time molecular dynamics simulations at constant pressure

In molecular dynamics simulations under periodic boundary conditions, particle positions are typically wrapped into a reference box. For diffusion coefficient calculations using the Einstein relation, the particle positions need to be unwrapped. Here, we show that a widely used heuristic unwrapping scheme is not suitable for long simulations at constant pressure. Improper accounting for box-volume fluctuations creates, at long times, unphysical trajectories and, in turn, grossly exaggerated diffusion coefficients. We propose an alternative unwrapping scheme that resolves this issue. At each time step, we add the minimal displacement vector according to periodic boundary conditions for the instantaneous box geometry. Here and in another paper [J. T. Bullerjahn, S. von Bülow, and G. Hummer, J. Chem. Phys. 153, 024116 (2020)], we apply the new unwrapping scheme to extensive molecular dynamics and Brownian dynamics simulation data. We provide practitioners with a formula to assess if and by how much earlier results might have been affected by the widely used heuristic unwrapping scheme.