PySteMoDA: An Open-Source Python Package for the Analysis of Steered Molecular Dynamics Simulations Data

Scorpion venom is a rich source of toxins which have great potential to develop new therapeutic agents. Scorpion chloride channel toxins (ClTxs), such as Chlorotoxin selectively inhibit human Matrix Methaloproteinase-2 (hMMP-2). The inhibitors of hMMP-2 have potential use in cancer therapy. Three new ClTxs, meuCl14, meuCl15 and meuCl16, derived from the venom transcriptome of Iranian scorpion, M. eupeus (Buthidea family), show high sequence identity (71.4%) with Chlorotoxin. Here, 3-D homology model of new ClTxs were constructed. The models were optimized by Molecular Dynamics simulation based on MDFF (molecular dynamics flexible fitting) method. New ClTxs indicate the presence of CSαβ folding of other scorpion toxins. A docking followed by steered molecular dynamics (SMD) simulations to investigate the interactions of meuCl14, meuCl15, and meuCl16 with hMMP-2 was applied. The current study creates a correlation between the unbinding force and the inhibition activities of meuCl14, meuCl15 and meuCl16 to shed some insights as to which toxin may be used as a drug deliverer. To this aim, SMD simulations using Constant Force Pulling method were carried out. The SMD provided useful details related to the changes of electrostatic, van de Waals (vdW), and hydrogen-bonding (H-bonding) interactions between ligands and receptor during the pathway of unbinding. According to SMD results, the interaction of hMMP-2 with meuCl14 is more stable. In addition, Arginine residue was found to contribute significantly in interaction of ClTxs with hMMP-2. All in all, the present study is a dynamical approach whose results are capable of being implemented in structure-based drug design.

This study explores the unbinding dynamics of alpha-ketoglutarate (AKG) from wild-type and mutant IDH1/IDH2 enzymes through steered molecular dynamics (SMD) simulations, examining how mutations influence binding, stability and enzymatic behaviour. Isocitrate dehydrogenase (IDH) enzymes are essential for cellular metabolism, catalyzing the conversion of isocitrate to AKG in the tricarboxylic acid cycle. Mutations in IDH1 and IDH2 lead to the aberrant accumulation of the oncometabolite 2-hydroxyglutarate (2-HG), disrupting normal metabolic processes and contributing to tumorigenesis. SMD simulations were employed to investigate AKG unbinding from both wild-type and mutant IDH1/IDH2. External forces were applied to quantify rupture forces and assess differences in stability among enzyme variants. Wild-type IDH1 exhibited strong and stable AKG interactions, reflected by higher rupture forces and a greater number of hydrogen bonds, consistent with its normal catalytic function. In contrast, the R132H mutation in IDH1 weakened AKG binding, facilitating dissociation and potentially promoting 2-HG formation. Among IDH2 variants, the R140Q mutant demonstrated lower binding stability compared to R172K, while the wild-type enzyme maintained stronger interactions. Mutations in IDH1 and IDH2 disrupt AKG binding and alter the stability, which may contribute to the pathological accumulation of 2-HG. These findings provide molecular insights into the oncogenic effects of IDH mutations and may aid in the development of targeted therapeutic strategies to inhibit mutant enzyme activity in cancer.

A variety of 10 polysaccharide motifs comprising naturally inspired sequences of monosaccharide building blocks are studied to understand the inter-relation between their structural and mechanical behaviors. Equilibrium and steered molecular dynamics (SMD) simulations are employed to investigate the stress-strain relationships and the associated conformational flips of the pyranose moieties along the polysaccharide chains. The presence of a variety of glycosidic linkages connecting the diverse monosaccharide units along with chain-branching in some cases induce wide diversity in the carbohydrate-Ramachandran plots of the glycosidic dihedrals. Similar variations are observed in the Cremer-Pople ring puckering patterns across the polysaccharide variants. The work further provides a comparison between the experimentally obtained atomic force microscopic data of mechanical stretching for some polysaccharides with the stress-strain curves generated from our SMD simulations. Out of all the systems studied, pectin having an axial-axial orientation of the glycosidic linkage showed maximum stretching potential, while acetan-M, with an equatorial-equatorial disposition of the glycosidic bond, stretched the least. The experimental Young's modulus of the corresponding natural polysaccharides could be reasonably compared to the values obtained from our simulation models. Force distribution analysis is done to understand the propagation of punctual stress in the polysaccharides under SMD conditions. Changes in local electrophilicity or nucleophilicity of atomic centers in puckered pyranose rings are estimated through the condensed Fukui functions. All of this information can help understand the physical behavior and chemical reactivity of complex polysaccharides in a complicated milieu of electronic and steric effects experienced by them.

Dynamic information, such as force, structural change, interaction energy, and potential of mean force (PMF), about the desorption of a single cardiotoxin (CTX) protein from a methyl-terminated self-assembled monolayer (SAM) surface was investigated by means of steered molecular dynamics (SMD) simulations. The simulation results indicated that Loop I is the first loop to depart from the SAM surface, which is in good agreement with the results of the nuclear magnetic resonance spectroscopy experiment. The free energy landscape and the thermodynamic force of the CTX desorption process was represented by the PMF and by the derivative of PMF with respect to distance, respectively. By applying Jarzynski's equality, the PMF can be reconstructed from the SMD simulation. The PMFs, calculated by different estimators based upon Jarzynski's equality, were compared with the conventional umbrella sampling method. The best estimation was obtained by using the fluctuation-dissipation estimator with a pulling velocity of v = 0.25 nm/ns for the present study.

The ability of two heteroditopic calix[4]diquinone receptors to transport a KCl ion-pair and a dopamine zwitterion through a water-chloroform interface was investigated via molecular dynamics (MD) simulations. Gas-phase conformational analysis has been carried on KCl and dopamine receptor binding associations and the lowest energy structures found in both cases show that the recognition of KCl and dopamine zwitterion occurs through multiple and cooperative N-H...anion and O...cation bonding interactions, with the receptor adopting equivalent folded conformations stabilized by pi-stacking interactions. The unconstrained MD simulations performed on KCl and dopamine complexes inserted in either the chloroform or water phase revealed that receptors are preferentially located at the interface with the hydrophobic tert-butyl groups of the calix[4]diquinone moiety immersed in the chloroform bulk while the polar anion binding cavity is directed toward the water phase. When the KCl complex is placed in chloroform, the release of the ion-pair occurs only after the first contact with the water interface, being a nonsimultaneous event, with the chloride anion leaving the receptor before the potassium cation. The dopamine, via the -NH(3)(+) binding entity, remains bound to the receptor during the entire time of the MD simulation (10 ns). In contrast, when both complexes were inserted in the water bulk, the full release of KCl and dopamine are fast events. The potentials of mean force (PMFs), associated with the migration of the complexes from chloroform to water through the interface, were calculated from steered molecular dynamics (SMD) simulations. The PMFs for the free KCl and zwitterionic dopamine migrations were also obtained for comparison purposes. The transport of KCl from water to chloroform (the reverse path) mediated by the receptor has a free energy barrier estimated in 6.50 kcal mol(-1), which is 3.0 kcal mol(-1) smaller than that found for the free KCl. The transport of dopamine complex along the reverse path is characterized by downhill energy profile, with a small free energy barrier of 6.56 kcal mol(-1).

Titin, a approximately 1 micron long protein found in striated muscle myofibrils, possesses unique elastic properties. The extensible behavior of titin has been demonstrated in atomic force microscopy and optical tweezer experiments to involve the reversible unfolding of individual immunoglobulin-like (Ig) domains. We have used steered molecular dynamics (SMD), a novel computer simulation method, to investigate the mechanical response of single titin Ig domains upon stress. Simulations of stretching Ig domains I1 and I27 have been performed in a solvent of explicit water molecules. The SMD approach provides a detailed structural and dynamic description of how Ig domains react to external forces. Validation of SMD results includes both qualitative and quantitative agreement with AFM recordings. Furthermore, combining SMD with single molecule experimental data leads to a comprehensive understanding of Ig domains' mechanical properties. A set of backbone hydrogen bonds that link the domains' terminal beta-strands play a key role in the mechanical resistance to external forces. Slight differences in architecture permit a mechanical unfolding intermediate for I27, but not for I1. Refolding simulations of I27 demonstrate a locking mechanism.

Protein-peptide interactions play essential roles in many cellular processes and their structural characterization is the major focus of current experimental and theoretical research. Two decades ago, it was proposed to employ the steered molecular dynamics (SMD) to assess the strength of protein-peptide interactions. The idea behind using SMD simulations is that the mechanical stability can be used as a promising and an efficient alternative to computationally highly demanding estimation of binding affinity. However, mechanical stability defined as a peak in force-extension profile depends on the choice of the pulling direction. Here we propose an uncommon choice of the pulling direction along resultant dipole moment (RDM) vector, which has not been explored in SMD simulations so far. Using explicit solvent all-atom MD simulations, we apply SMD technique to probe mechanical resistance of ligand-receptor system pulled along two different vectors. A novel pulling direction—when ligand unbinds along the RDM vector—results in stronger forces compared to commonly used ligand unbinding along center of masses vector. Our observation that RDM is one of the factors influencing the mechanical stability of protein-peptide complex can be used to improve the ranking of binding affinities by using mechanical stability as an effective scoring function.

Transthyretin(TTR), as a tetrameric protein, functions as a neuroprotector. The native TTR homotetramer dissociates into dimers and monomers. Dimers and monomers self-assemble into amyloid fibrils, and this process can lead to some diseases. Native TTR homotetramer is a widely accepted model for TTR amyloid formation. In this study, simulations using molecular dynamics(MD) and steered MD(SMD) were performed to explore the mechanisms for glabridin(Glab), a specific inhibitor for TTR binding, for V30A mutant and wild-type(WT) TTR. MD simulation results indicate that, compared with Glab binding to WT and V30A mutant, the WT TTR could lead to the collapse of β-strands from Ser52 to His56 at chain A. This phenomenon facilitated the easy dissociation of chains A and C. Calculations of the binding free energy between the two chains showed that the V30A-Glab TTR complex displayed a lower binding energy than other systems(WT TTR and WT-Glab TTR). Then, SMD simulation was performed to ex-plore the unbinding pathway for Glab through the WT and V30A mutant TTR. The results show that Lys15(chain A) produced a hydrogen bond with Glab at the force peak via the WT TTR tunnel. Meanwhile, in the V30A TTR mutant, the hydrogen bond between Lys15(chain A) and Glab was broken at the force peak. This condition was beneficial for Glab to be taken off from the protein. Our theoretical results will be useful in designing a new specific inhibitor of TTR protein to control the TTR homotetramer dissociation.

Resolving Cadherin Catch Bond Formation at the Single Molecule Level

Self-assembled cyclic peptide nanotubes (CPNs) show a potential use in drug delivery. In this study, the CPN composed of (Trp-d-Leu)4-Gln-d-Leu was synthesized and tested for the transport of the antitumor drug 5-fluorouracil (5-FU). CPN-mediated release of 5-FU from liposomes experimentally tested the transportation function of the synthetic CPNs. To explore the transportation mechanism of CPNs, computational studies have been performed on the CPN models stacked by 8 subuints, including conventional molecular dynamics (CMD) simulations, and steered molecular dynamics (SMD) simulations in the environment of hydrated dimyristoylphosphatidylcholine (DMPC) lipid bilayer. Our CMD simulations demonstrated that the ortho-CPN is the most stable nanotube, in which the Gln residue is in the ortho-position relative to other residues. The calculated diffusion coefficient value for inner water molecules was 1.068 × 10−5 cm2·s−1, almost half that of the bulky water and 24 times faster than that of the typical gramicidin A channel. The CPN conserved its hollow structure along the 10 ns CMD simulations, with a tile angle of 50° relative to the normal of DMPC membrane. Results from SMD simulations showed that the 5-FU molecule was transported by hopping through different potential energy minima distributed along subunits, and finally exited the nanotube by escaping from the kink region at the last two subunits. The hopping of 5-FU was driven by switching from hydrophobic interactions between 5-FU and the interior wall of the nanotube to hydrogen bonding interactions of 5-FU with the backbone carbonyl group and amide group of ortho-CPN. The calculated binding free energy profile of 5-FU interacting with the CPN indicated that there was an energy well near the outer end of the nanotube.

Unbinding pathway energy of glyphosate from the EPSPs enzyme binding site characterized by Steered Molecular Dynamics and Potential of Mean Force

Calcium-dependent lectin I from Pseudomonas aeruginosa (PA-IL) binds specifically to oligosaccharides presenting an alpha-galactose residue at their nonreducing end, such as the disaccharides alphaGal1-2betaGalOMe, alphaGal1-3betaGalOMe, and alphaGal1-4betaGalOMe. This provides a unique model for studying the effect of the glycosidic linkage of the ligands on structure and thermodynamics of the complexes by means of experimental and theoretical tools. The structural features of PA-IL in complex with the three disaccharides were established by docking and molecular dynamics simulations and compared with those observed in available crystal structures, including PA-IL.alphaGal1-2betaGalOMe complex, which was solved at 2.4 A resolution and reported herein. The role of a structural bridge water molecule in the binding site of PA-IL was also elucidated through molecular dynamics simulations and free energy calculations. This water molecule establishes three very stable hydrogen bonds with O6 of nonreducing galactose, oxygen from Pro-51 main chain, and nitrogen from Gln-53 main chain of the lectin binding site. Binding free energies for PA-IL in complex with the three disaccharides were investigated, and the results were compared with the experimental data determined by titration microcalorimetry. When the bridge water molecule was included in the free energy calculations, the simulations predicted the correct binding affinity trends with the 1-2-linked disaccharide presenting three times stronger affinity ligand than the other two. These results highlight the role of the water molecule in the binding site of PA-IL and indicate that it should be taken into account when designing glycoderivatives active against P. aeruginosa adhesion.

Cellobiohydrolase A from Ruminiclostridium thermocellum (Cbh9A) is a processive exoglucanase from family 9 and is an important cellobiohydrolase that hydrolyzes cello-oligosaccharide into cellobiose. Residues Tyr555 and Trp678 considerably affect catalytic activity, but their mechanisms are still unknown. To investigate how the Tyr555 and Trp678 affect the processivity of Cbh9A, conventional molecular dynamics, steered molecular dynamics, and free energy calculation were performed to simulate the processive process of wild type (WT)-Cbh9A, Y555S mutant, and W678G mutant. Analysis of simulation results suggests that the binding free energies between the substrate and WT-Cbh9A are lower than those of Y555S and W678G mutants. The pull forces and energy barrier in Y555S and W678G mutants also reduced significantly during the steered molecular dynamics (SMD) simulation compared with that of the WT-Cbh9A. And the potential mean force calculations showed that the pulling energy barrier of Y555S and W678G mutants is much lower than that of WT-Cbh9A.

Simulations of lipid bilayers at different levels of detail is of major interest in computational biophysics. With a recently proposed algorithm, by mapping the information from a Molecular Dynamics (MD) simulation of a lipid bilayer on a 2D lattice model, it is possible to gain insight into the enthalpic and entropic contributions, governing the interaction of adjacent lipids and the individual chain entropy, respectively. The contributions are obtained as a function of the lipid chain order parameter. Here we use this approach for the case of pure DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphocholine) to compare the all-atom CHARMM36 DPPC lipid with the respective MARTINI coarse-grained DPPC lipid. We individually compare the enthalpy and entropy functions as well as the resulting free energies. This allows us to gain new insight into the comparison of both levels of description. When numerically solving the 2D lattice model via Monte Carlo (MC) simulations, the 2D model displays the gel/liquid transition at the same temperature as the respective MD simulation. The dramatic increase in efficiency of the 2D lattice model as compared to the MD simulation allows us to estimate the equilibrium order parameters of the gel phase, inaccessible by the MD simulations even after 3 μs.

The understanding of interactions of poly(ethylene glycol) (PEG) or poly(ethylene oxide) (PEO) with biological interfaces has important technological application in industry and in medicine. In this paper, structural and dynamical properties of PEO at the dimyristoylphospatidylcholine (DMPC) bilayer/water interface have been investigated by molecular dynamics (MD) and steered molecular dynamics (SMD) simulations. The structural properties of a PEO chain in bulk water, at the water/vacuum interface, and in the presence of the membrane were compared with available experimental data. The presence of a barrier for the PEO penetration into the DMPC bilayer has been found. A qualitative estimation of the barrier provided a value equal to approximately 19 kJ/mol, that is, 7 times the value of kT at 310 K.