Massive Molecular Dynamics Simulations Research Articles

StreaMD: the toolkit for high-throughput molecular dynamics simulations.

Molecular dynamics simulations serve as a prevalent approach for investigating the dynamic behaviour of proteins and protein-ligand complexes. Due to its versatility and speed, GROMACS stands out as a commonly utilized software platform for executing molecular dynamics simulations. However, its effective utilization requires substantial expertise in configuring, executing, and interpreting molecular dynamics trajectories. Existing automation tools are constrained in their capability to conduct simulations for large sets of compounds with minimal user intervention, or in their ability to distribute simulations across multiple servers. To address these challenges, we developed a Python-based tool that streamlines all phases of molecular dynamics simulations, encompassing preparation, execution, and analysis. This tool minimizes the required knowledge for users engaging in molecular dynamics simulations and can efficiently operate across multiple servers within a network or a cluster. Notably, the tool not only automates trajectory simulation but also facilitates the computation of free binding energies for protein-ligand complexes and generates interaction fingerprints across the trajectory. Our study demonstrated the applicability of this tool on several benchmark datasets. Additionally, we provided recommendations for end-users to effectively utilize the tool.Scientific contributionThe developed tool, StreaMD, is applicable to different systems (proteins, ligands and their complexes including co-factors) and requires a little user knowledge to setup and run molecular dynamics simulations. Other features of StreaMD are seamless integration with calculation of MM-GBSA/PBSA binding free energies and protein-ligand interaction fingerprints, and running of simulations within distributed environments. All these will facilitate routine and massive molecular dynamics simulations.

Fluorine Domains Induced Ultrahigh Nitrogen Solubility in Ionic Liquids.

Fluorinated ionic liquids (ILs) are well-known as electrolytes in the nitrogen (N2) electroreduction reaction due to their exceptional gas solubility. However, the influence of fluorinated functional group on N2 solvation and solubility enhancement remains unclear. Massive molecular dynamics simulations and free energy perturbation methods are conducted to investigate the N2 solubility in 11 traditional and 9 fluorinated ILs. It shows that the fluorinated IL of 1-Ethyl-3-methylimidazolium tris(pentafluoroethyl) trifluorophosphate ([Emim]FAP) exhibits ultrahigh solubility, 4.844 × 10-3, approximately 118 times higher than that of traditional IL 1-Ethyl-3-methylimidazolium nitrate ([Emim]NO3). Moreover, fluorinated ILs with more than 10 C-F bonds possess higher N2 solubility than others and show an exothermic nature during solvation. As the C-F bonds number in ILs decreases, the N2 solubility decreases significantly and displays the opposite endothermic behavior. To understand the ultrahigh N2 solubility in fluorinated ILs, we propose a concept of fluorine densification energy (FDE), referring to the average strength of interaction between atoms per unit volume in ILs with fluorine domains, demonstrating a linear relationship with C-F bonds. Physically, lower FDE results in lower N2-anion pair dissociation energy and higher free volume, finally enhancing the N2 solubility. Consequently, medium to long alkyl fluorine tails within a polar environment defines a distinct fluorine domain, emphasizing FDE's role in enhancing N2 solubility. Overall, these quantitative results will not only deepen the understanding of N2 solvation in ILs but may also shed light on the rational design of IL-based high-performance N2 capture and conversion technologies.

Quasiballistic thermal transport in submicron-scale graphene nanoribbons at room-temperature.

Phonon transport in two-dimensional materials has been the subject of intensive studies both theoretically and experimentally. Recently observed unique phenomena such as Poiseuille flow at low temperature in graphene nanoribbons (GNRs) initiated strong interest in similar effects at higher temperatures. Here, we carry out massive molecular dynamics simulations to examine thermal transport in GNRs at room temperature (RT) and demonstrate that non-diffusive behaviors including Poiseuille-like local thermal conductivity and second sound are obtained, indicating quasiballistic thermal transport. For narrow GNRs, a Poiseuille-like thermal conductivity profile develops across the nanoribbon width, and wider GNRs exhibit a mixed nature of phonon transport in that diffusive transport is dominant in the middle region whereas non-uniform behavior is observed near lateral GNR boundaries. In addition, transient heating simulations reveal that the driftless second sound can propagate through GNRs regardless of the GNR width. By decomposing the atomic motion into out-of-plane and in-plane modes, it is further shown that the observed quasiballistic thermal transport is primarily contributed by the out-of-plane motion of C atoms in GNRs.

Designing artificial ion channels with strict K+/Na+ selectivity toward next-generation electric-eel-mimetic ionic power generation.

A biological potassium channel is >1000 times more permeable to K+ than to Na+ and exhibits a giant permeation rate of ∼108 ions/s. It is a great challenge to construct artificial potassium channels with such high selectivity and ion conduction rate. Herein, we unveil a long-overlooked structural feature that underpins the ultra-high K+/Na+ selectivity. By carrying out massive molecular dynamics simulation for ion transport through carbonyl-oxygen-modified bi-layer graphene nanopores, we find that the twisted carbonyl rings enable strict potassium selectivity with a dynamic K+/Na+ selectivity ratio of 1295 and a K+ conduction rate of 3.5×107 ions/s, approaching those of the biological counterparts. Intriguingly, atomic trajectories of K+ permeation events suggest a dual-ion transport mode, i.e. two like-charged potassium ions are successively captured by the nanopores in the graphene bi-layer and are interconnected by sharing one or two interlayer water molecules. The dual-ion behavior allows rapid release of the exiting potassium ion via a soft knock-on mechanism, which has previously been found only in biological ion channels. As a proof-of-concept utilization of this discovery, we propose a novel way for ionic power generation by mixing KCl and NaCl solutions through the bi-layer graphene nanopores, termed potassium-permselectivity enabled osmotic power generation (PoPee-OPG). Theoretically, the biomimetic device achieves a very high power density of >1000W/m2 with graphene sheets of <1% porosity. This study provides a blueprint for artificial potassium channels and thus paves the way toward next-generation electric-eel-mimetic ionic power generation.

Dynamics of lattice disorder in perovskite materials, polarization nanoclusters and ferroelectric domain wall structures

The nexus between classic ferroelectricity and the structure of perovskite materials hinges on the concept of lattice disorder. Although the ordered perovskites display short-range displacements of the central cations around their equilibrium points, the lattice disorder dynamically unfolds to generate a myriad of distorted rhombohedral lattices characterized by the hopping of the central cations across <111> directions. It is discovered that the lattice disorder correlates with the emergence of minimum configuration energy <100> pathways for the central cations, resulting in spatially modulated ultrafast polarization nanocluster arrangements that are stabilized by the electric charge defects in the material. Through high-resolution phonon dispersion analyses encompassing molecular dynamics (MD) and density functional theory (DFT) simulations, we provide unequivocal evidence linking the hopping of central cations to the development of diffuse soft phonon modes observed throughout the phase transitions of the perovskite. Through massive MD simulations, we unveil the impact of lattice disorder on the structures of domain walls at finite-temperature vis-à-vis collective activation and deactivation of <100> pathways. Furthermore, our simulations demonstrate the development of hierarchical morphotropic phase boundary (MPB) nanostructures under the combined influence of externally applied pressure and stress relaxation, characterized by sudden emergence of zig-zagged monoclinic arrangements that involve dual <111> shifts of the central cations. These findings have implications for tailoring MPBs in thin-film structures and for the light-induced mobilization of DWs. Avenues are finally uncovered to the exploration of lattice disorder through gradual shear strain application.

The thiophene α-terthienylmethanol isolated from Tagetes minuta inhibits angiogenesis by targeting protein kinase C isozymes α and β2.

Background: Tumor angiogenesis is considered as a crucial pathologic feature of cancer with a key role in multidrug resistance (MDR). Adverse effects of the currently available drugs and the development of resistance to these remain as the hardest obstacles to defeat. Objetive: This work explores flora from Argentina as a source of new chemical entities with antiangiogenic activity. Methods: Tube formation assay using bovine aortic endothelial cells (BAECs) was the experiment of choice to assess antiangiogenic activity. The effect of the pure compound in cell invasiveness was investigated through the trans-well migration assay. The inhibitory effect of the pure compound on VEGFR-2 and PKC isozymes α and β2 activation was studied by molecular and massive dynamic simulations. Cytotoxicity on peripheral blood mononuclear cells and erythrocyte cells was evaluated by means of MTT and hemolysis assay, respectively. In silico prediction of pharmacological properties (ADME) and evaluation of drug-likeness features were performed using the SwissADME online tool. Results: Among the plants screened, T. minuta, showed an outstanding effect with an IC50 of 33.6 ± 3.4 μg/ml. Bio-guided isolation yielded the terthiophene α-terthienylmethanol as its active metabolite. This compound inhibited VEGF-induced tube formation with an IC50 of 2.7 ± 0.4 μM and significantly impaired the invasiveness of bovine aortic endothelial cells (BAECs) as well as of the highly aggressive breast cancer cells, MDA-MB-231, when tested at 10 μM. Direct VEGFR-2 and PKC inhibition were both explored by means of massive molecular dynamics simulations. The results obtained validated the inhibitory effect on protein kinase C (PKC) isozymes α and β2 as the main mechanism underlying its antiangiogenic activity. α-terthienylmethanol showed no evidence of toxicity against peripheral blood mononuclear and erythrocyte cells. Conclusion: These findings support this thiophene as a promising antiangiogenic phytochemical to fight against several types of cancer mainly those with MDR phenotype.

Entropy driving highly selective CO2 separation in nanoconfined ionic liquids

Nowadays, the global greenhouse effect has led to the imminent development of CO2 capture, separation, and storage technologies. Hybrid membranes with nanoconfined ionic liquids (ILs) show great potential for CO2 separation, but the intrinsic mechanism is still obscure. Herein, the thermodynamical properties and solvating processes of CO2 and CH4 in ILs confined in graphene oxide were studied via performing massive molecular dynamics simulations. It was first identified that selectivity rises from 25.01 to 149.20 as the interlayer distance decreases from 3.00 to 1.50 nm, showing an ultrahigh separating selectivity. Interestingly, the solubility of CO2 in confined ILs increases by almost two orders of magnitude compared with that in bulk ILs, which is far larger than CH4 in confined ILs. The high solubility mainly originates from the fact that the confined ILs can induce the structure rearrangement and provide abundant CO2 adsorbing sites, raising the configurational entropy of CO2 in the confined ILs, and further driving the high separation selectivity of CO2 over CH4. Finally, quantitative relations between solubility, diffusion capacity, permeability, selectivity, and structural entropy of gas in confined ILs are constructed, which are meaningful for the theoretical understanding, rational design, and applications of highly efficient and low-cost separation of CO2.

Dynamic and equilibrium properties of finite-size polymer models of chromosome folding.

Novel technologies are revealing that chromosomes have a complex three-dimensional organization within the cell nucleus that serves functional purposes. Models from polymer physics have been developed to quantitively understand the molecular principles controlling their structure and folding mechanisms. Here, by using massive molecular-dynamics simulations we show that classical scaling laws combined with finite-size effects of a simple polymer model can effectively explain the scaling behavior that chromatin exhibits at the topologically associating domains level, as revealed by experimental observations. Model results are then validated against recently published high-resolution in situ Hi-C data.

Understanding imprint formation, plastic instabilities and hardness evolutions in FCC, BCC and HCP metal surfaces

Nanoindentation experiments in metal surfaces are characterized by the onset of plastic instabilities along with the development of permanent nanoimprints and dense defect networks. This investigation concerns massive molecular dynamics simulations of nanoindentation experiments in FCC, BCC and HCP metals using blunted (spherical) tips of realistic size, and the detailed comparison of the results with experimental measurements. Our findings shed light on the defect processes which dictate the contact resistance to plastic deformation, the development of a transitional stage with abrupt plastic instabilities, and the evolution towards a self-similar steady-state characterized by the plateauing hardness p p at constant dislocation density ρ p . The onset of permanent nanoimprints is governed by stacking fault and nanotwin interlocking, the buildup of nanostructured regions and crystallites throughout the imprint, the cross-slip and cross-kinking of surfaced screw dislocations, and the occurrence of defect remobilization events within the plastic zone. As a result of these mechanisms, the ratio between the hardness p p and the Young's modulus E becomes higher in BCC Ta and Fe, followed by FCC Al, HCP Mg and large stacking fault width FCC Ni and Cu. Finally, when nanoimprint formation is correlated with the uniaxial response of the indented minuscule material volume, the hardness to yield strength ratio, p p / σ ys , varies from ≈ 7 to ≈ 10, which largely exceeds the continuum plasticity bound of ≈ 2.8. Our results have general implications to the understanding of indentation size-effects, where the onset of extreme nanoscale hardness values is associated with the occurrence of unique imprint-forming processes under large strain gradients.

Phonon dispersion relation of metallic glasses

Experimental data on the phase sound speed of metallic glasses show anomalies in the terahertz range, reflecting an underlying complex behavior of their phonon dispersion spectrum not yet explained. We determine the phonon dispersion curve of metallic glasses by means of massive molecular dynamics simulations, allowing us to obtain the $\text{low-}q$ region behavior with unprecedented detail. Results confirm that the sound speed is constant below the THz range, down to the macroscopic limit. On the contrary, a hardening of the sound speed, more notable in the transverse case, is found in the THz range. This behavior is modeled in terms of a relaxation model. The model gives quantitative agreement and allows us to determine a new threshold frequency ${\ensuremath{\omega}}_{h}$, at the end of the boson-peak region. Above ${\ensuremath{\omega}}_{h}$ the shear modulus increases dramatically, reflecting the end of the amorphous-like acoustic propagation region characterized by the excess density of vibrational states.

Heat dissipation guides activation in signaling proteins

Life is fundamentally a nonequilibrium phenomenon. At the expense of dissipated energy, living things perform irreversible processes that allow them to propagate and reproduce. Within cells, evolution has designed nanoscale machines to do meaningful work with energy harnessed from a continuous flux of heat and particles. As dictated by the Second Law of Thermodynamics and its fluctuation theorem corollaries, irreversibility in nonequilibrium processes can be quantified in terms of how much entropy such dynamics produce. In this work, we seek to address a fundamental question linking biology and nonequilibrium physics: can the evolved dissipative pathways that facilitate biomolecular function be identified by their extent of entropy production in general relaxation processes? We here synthesize massive molecular dynamics simulations, Markov state models (MSMs), and nonequilibrium statistical mechanical theory to probe dissipation in two key classes of signaling proteins: kinases and G-protein-coupled receptors (GPCRs). Applying machinery from large deviation theory, we use MSMs constructed from protein simulations to generate dynamics conforming to positive levels of entropy production. We note the emergence of an array of peaks in the dynamical response (transient analogs of phase transitions) that draw the proteins between distinct levels of dissipation, and we see that the binding of ATP and agonist molecules modifies the observed dissipative landscapes. Overall, we find that dissipation is tightly coupled to activation in these signaling systems: dominant entropy-producing trajectories become localized near important barriers along known biological activation pathways. We go on to classify an array of equilibrium and nonequilibrium molecular switches that harmonize to promote functional dynamics.

Molecular Dynamics Simulations of Ibuprofen Release from pH-Gated Silica Nanochannels.

The iboprufen delivery process from cylindrical silica pores of 3 nm diameter, with polyamine chains anchored at the pore outlets, was investigated by means of massive molecular dynamics simulations. Effects from pH were introduced by considering polyamine chains with different degrees of protonation. High, low, and intermediate pH environments were investigated. The increment of the acidity of the environment leads to a significant decrease of the pore aperture, yielding an effective diameter, for the lowest pH case, that is 3.5 times smaller than the one associated with the highest pH one. Using a biased sampling procedure, Gibbs free energy profiles for the ibuprofen delivery process were obtained. The joint analysis of the corresponding profiles, time evolution of the ibuprofen position within the channel, orientation of the molecule, and instantaneous effective diameter of the gate suggest a three-step mechanism for ibuprofen delivery. A complementary analysis of the translational mobility of ibuprofen along the axial direction of the channel revealed a subdiffusive dynamics in the low and intermediate pH cases. Deviations from Brownian diffusive dynamics are discussed and compared with direct experimental results.

Effects of oncogenic mutations on the conformational free-energy landscape of EGFR kinase

Activating mutations in the epidermal growth factor receptor (EGFR) tyrosine kinase are frequently found in many cancers. It has been suggested that changes in the equilibrium between its active and inactive conformations are linked to its oncogenic potential. Here, we quantify the effects of some of the most common single (L858R and T790M) and double (T790M-L858R) oncogenic mutations on the conformational free-energy landscape of the EGFR kinase domain by using massive molecular dynamics simulations together with parallel tempering, metadynamics, and one of the best force-fields available. Whereas the wild-type EGFR catalytic domain monomer is mostly found in an inactive conformation, our results show a clear shift toward the active conformation for all of the mutants. The L858R mutation stabilizes the active conformation at the expense of the inactive conformation and rigidifies the αC-helix. The T790M gatekeeper mutant favors activation by stabilizing a hydrophobic cluster. Finally, T790M with L858R shows a significant positive epistasis effect. This combination not only stabilizes the active conformation, but in nontrivial ways changes the free-energy landscape lowering the transition barriers.

From Quantity To Quality: Massive Molecular Dynamics Simulation of Nanostructures under Plastic Ddeformation in Desktop and Service Grid Distributed Computing Infrastructure

The distributed computing infrastructure (DCI) on the basis of BOINC and EDGeS-bridge technologies for high-performance distributed computing is used for porting the sequential molecular dynamics (MD) application to its parallel version for DCI with Desktop Grids (DGs) and Service Grids (SGs). The actual metrics of the working DG-SG DCI were measured, and the normal distribution of host performances, and signs of log-normal distributions of other characteristics (CPUs, RAM, and HDD per host) were found. The practical feasibility and high efficiency of the MD simulations on the basis of DG-SG DCI were demonstrated during the experiment with the massive MD simulations for the large quantity of aluminum nanocrystals ($\sim10^2$-$10^3$). Statistical analysis (Kolmogorov-Smirnov test, moment analysis, and bootstrapping analysis) of the defect density distribution over the ensemble of nanocrystals had shown that change of plastic deformation mode is followed by the qualitative change of defect density distribution type over ensemble of nanocrystals. Some limitations (fluctuating performance, unpredictable availability of resources, etc.) of the typical DG-SG DCI were outlined, and some advantages (high efficiency, high speedup, and low cost) were demonstrated. Deploying on DG DCI allows to get new scientific $\it{quality}$ from the simulated $\it{quantity}$ of numerous configurations by harnessing sufficient computational power to undertake MD simulations in a wider range of physical parameters (configurations) in a much shorter timeframe.

The Different Flexibility of c-Src and c-Abl Kinases Regulates the Accessibility of a Druggable Inactive Conformation

c-Src and c-Abl are two closely related protein kinases that constitute important anticancer targets. Despite their high sequence identity, they show different sensitivities to the anticancer drug imatinib, which binds specifically to a particular inactive conformation in which the Asp of the conserved DFG motif points outward (DFG-out). We have analyzed the DFG conformational transition of the two kinases using massive molecular dynamics simulations, free energy calculations, and isothermal titration calorimetry. On the basis of the reconstruction of the free energy surfaces for the DFG-in to DFG-out conformational changes of c-Src and c-Abl, we propose that the different flexibility of the two kinases results in a different stability of the DFG-out conformation and might be the main determinant of imatinib selectivity.

Relating Molecular Interactions with N-BAR Domains to the Mesoscopic Nature of Membrane Remodeling

Theoretical approaches to studying biological problems have seen appreciable advancements in the past decade, allowing detailed structural and thermodynamic description of fundamental cellular processes, such as protein folding, lipid self-assembly, large-scale structural rearrangements of macromolecules, etc. Nevertheless, a full atomic description of even simple microscopic biological systems requires billions of atoms, beyond the reach of current computational methods. Membrane remodeling induced by members of the BAR domain protein family, is an innately multiscale problem, in which molecular interactions between the BAR protein and the lipids induces local curvature, and ultimately leads to large-scale reticulations of liposomes. In our study, we start with a discretized field-theoretical description of a liposome and carry out continuum mechanics simulations that generate reticulated topologies, quite similar to the ones seen in experiments. The obtained configurations are used as templates for mapping coarse-grained lipids and BARs over the continuum model. Subsequently, we use the coarse-grained representation of the liposome to run massive molecular dynamics simulations, with the aim of capturing the underlying molecular mechanisms that direct the liposome to undergo such large-scale biological restructuring and detecting the preferred localization of BARs at high resolution.View Large Image | View Hi-Res Image | Download PowerPoint Slide

VIRTUAL LABORATORY MOLDYNGRID AS A PART OF SCIENTIFIC INFRASTRUCTURE FOR BIOMOLECULAR SIMULATIONS

Virtual laboratory MolDynGrid has been established for computational researches in structural biology and bioinformatics, especially for molecular dynamics (MD) simulations of biological macromolecules and their complexes. Biomolecular simulation usually requires a very high processing power and huge storage space for MD trajectories. In this study we describe the development of virtual laboratory infrastructure to run massive MD simulations on a Grid infrastructure in a distributed manner. The aim of MolDynGrid development was to provide an efficient infrastructure for automation of MD job processing in Grid. It will assist researchers to use grid technology without having them to be familiar with highly complicated internals of the Grid. Web-portal of MolDynGrid virtual laboratory has been developed (http://moldyngrid.org), which consists of computational, analytical and educational blocks аnd data base of MD simulations.

High-Performance Drug Discovery: Computational Screening by Combining Docking and Molecular Dynamics Simulations

Virtual compound screening using molecular docking is widely used in the discovery of new lead compounds for drug design. However, this method is not completely reliable and therefore unsatisfactory. In this study, we used massive molecular dynamics simulations of protein-ligand conformations obtained by molecular docking in order to improve the enrichment performance of molecular docking. Our screening approach employed the molecular mechanics/Poisson-Boltzmann and surface area method to estimate the binding free energies. For the top-ranking 1,000 compounds obtained by docking to a target protein, approximately 6,000 molecular dynamics simulations were performed using multiple docking poses in about a week. As a result, the enrichment performance of the top 100 compounds by our approach was improved by 1.6–4.0 times that of the enrichment performance of molecular dockings. This result indicates that the application of molecular dynamics simulations to virtual screening for lead discovery is both effective and practical. However, further optimization of the computational protocols is required for screening various target proteins.

Uncommon features in the dynamic structure factor of a molten transition metal

The collective SC(Q,ω), single-particle SS(Q,ω) dynamic structure factors as well as the spectra for transverse-current correlations CT(Q,ω) for molten nickel close to melting have been calculated from a massive molecular dynamics simulation using an Embedded Atom Model potential. A number of features are found in the structure factors that are absent in the spectra of simpler liquids such as molten alkali metals. Estimates for the shear component of the generalized, wavevector-dependent viscosity are obtained from analysis of CT(Q,ω). Such data together with those concerning the strength of the thermal conduction channels for the excitation decay enable us to subtract such contributions to form the linewidths of Brillouin peaks and thus to evaluate the Q-dependent bulk viscosity coefficient. The result shows that contrary to inferences based upon geophysical data, both bulk and shear-viscosity coefficients are of the same order of magnitude.

Switch interactions control energy frustration and multiple flagellar filament structures

Bacterial flagellar filament is a macromolecular assembly consisting of a single protein, flagellin. Bacterial swimming is controlled by the conformational transitions of this filament between left- and right-handed supercoils induced by the flagellar motor torque. We present a massive molecular dynamics simulation that was successful in constructing the atomic-level supercoil structures consistent with various experimental data and further in elucidating the detailed underlying molecular mechanisms of the polymorphic supercoiling. We have found that the following three types of interactions are keys to understanding the supercoiling mechanism. "Permanent" interactions are always maintained between subunits in the various supercoil structures. "Sliding" interactions are formed between variable hydrophilic or hydrophobic residue pairs, allowing intersubunit shear without large change in energy. The formation and breakage of "switch" interactions stabilize inter- and intrasubunit interactions, respectively. We conclude that polymorphic supercoiling is due to the energy frustration between them. The transition between supercoils is achieved by a "transform and relax" mechanism: the filament structure is geometrically transformed rapidly and then slowly relaxes to energetically metastable states by rearranging interactions.