A High-Throughput Steered Molecular Dynamics Study on the Free Energy Profile of Ion Permeation through Gramicidin A

Self-organizing maps of unbiased ligand-target binding pathways and kinetics

The interpretation of ligand-target interactions at atomistic resolution is central to most efforts in computational drug discovery and optimization. However, the highly dynamic nature of protein targets, as well as possible induced fit effects, makes difficult to sample many interactions effectively with docking studies or even with large-scale molecular dynamics (MD) simulations. We propose a novel application of Self-Organizing Maps (SOMs) to address the sampling and dynamic mapping tasks, particularly in cases involving ligand flexibility and induced fit. The SOM approach offers a data-driven strategy to create a map of the interaction process and pathways based on unbiased MD. Furthermore, we show how the preliminary SOM mapping is complementary to kinetic analysis, with the employment of both network-based approaches and Markov state models. We demonstrate the method by comprehensively mapping a large dataset of 640 μs of unbiased trajectories sampling the recognition process between the phosphorylated YEEI peptide and its high-specificity target lck-SH2. The integration of SOM into unbiased simulation protocols significantly advances our understanding of the ligand binding mechanism. This approach serves as a potent tool for mapping intricate ligand-target interactions with unprecedented detail, thereby enhancing the characterization of kinetic properties crucial to drug design.

Mechanism of phosphate release from actin filaments

After ATP-actin monomers assemble filaments, the ATP's [Formula: see text]-phosphate is hydrolyzedwithin seconds and dissociates over minutes. We used all-atom molecular dynamics simulations to sample the release of phosphate from filaments and study residues that gate release. Dissociation of phosphate from Mg2+ is rate limiting and associated with an energy barrier of 20 kcal/mol, consistent with experimental rates of phosphate release. Phosphate then diffuses within an internal cavity toward a gate formed by R177, as suggested in prior computational studies and cryo-EM structures. The gate is closed when R177 hydrogen bonds with N111 and is open when R177 forms a salt bridge with D179. Most of the time, interactions of R177 with other residues occlude the phosphate release pathway. Machine learning analysis reveals that the occluding interactions fluctuate rapidly, underscoring the secondary role of backdoor gate opening in Pi release, in contrast with the previous hypothesis that gate opening is the primary event.

Synthetic K+ Channels Constructed by Rebuilding the Core Modules of Natural K+ Channels in an Artificial System

Different types of natural K⁺ channels share similar core modules and cation permeability characteristics. In this study, we have developed novel artificial K⁺ channels by rebuilding the core modules of natural K⁺ channels in artificial systems. All the channels displayed high selectivity for K⁺ over Na⁺ and exhibited a selectivity sequence of K⁺ ≈Rb⁺ during the transport process, which is highly consistent with the cation permeability characteristics of natural K⁺ channels. More importantly, these artificial channels could be efficiently inserted into cell membranes and mediate the transmembrane transport of K⁺ , disrupting the cellular K⁺ homeostasis and eventually triggering the apoptosis of cells. These findings demonstrate that, by rebuilding the core modules of natural K⁺ channels in artificial systems, the structures, transport behaviors, and physiological functions of natural K⁺ channels can be mimicked in synthetic channels.

Predicting Ion Channel Conductance via Dissipation-Corrected Targeted Molecular Dynamics and Langevin Equation Simulations

Ion channels are important proteins for physiological information transfer and functional control. To predict the microscopic origins of their voltage-conductance characteristics, here we applied dissipation-corrected targeted molecular dynamics in combination with Langevin equation simulations to potassium diffusion through the gramicidin A channel as a test system. Performing a nonequilibrium principal component analysis on backbone dihedral angles, we find coupled protein-ion dynamics to occur during ion transfer. The dissipation-corrected free energy profiles correspond well to predictions from other biased simulation methods. The incorporation of an external electric field in Langevin simulations enables the prediction of macroscopic observables in the form of I-V characteristics.

Metadynamics-Based Approaches for Modeling the Hypoxia-Inducible Factor 2α Ligand Binding Process

Several methods based on enhanced-sampling molecular dynamics have been proposed for studying ligand binding processes. Here, we developed a protocol that combines the advantages of steered molecular dynamics (SMD) and metadynamics. While SMD is proposed for investigating possible unbinding pathways of the ligand and identifying the preferred one, metadynamics, with the path collective variable (PCV) formalism, is suggested to explore the binding processes along the pathway defined on the basis of SMD, by using only two CVs. We applied our approach to the study of binding of two known ligands to the hypoxia-inducible factor 2α, where the buried binding cavity makes simulation of the process a challenging task. Our approach allowed identification of the preferred entrance pathway for each ligand, highlighted the features of the bound and intermediate states in the free-energy surface, and provided a binding affinity scale in agreement with experimental data. Therefore, it seems to be a suitable tool for elucidating ligand binding processes of similar complex systems.

Potential of mean force and umbrella sampling simulation for the transport of 5-oxazolidinone in heterotetrameric sarcosine oxidase

The structure of heterotetrameric sarcosine oxidase (HSO) contains a highly complex system composed of a large cavity and tunnels, which are essential for the reaction and migration of the reactants, products, and intermediates. Previous geometrical analysis using the CAVER program has predicted that there are three possible tunnels, T1, T2, and T3, for the exit pathway of the iminium intermediate, 5-oxazolidinone (5-OXA), of the enzyme reaction. Previous molecular dynamics (MD) simulation of HSO has identified the regions containing the water channels from the density distribution of water. The simulation indicated that tunnel T3 is the most probable exit pathway of 5-OXA. In the present study, the potential of mean force (PMF) for the transport of 5-OXA through tunnels T1, T2, and T3 was calculated using umbrella sampling (US) MD simulations and the weighted histogram analysis method. The PMF profiles for the three tunnels support the notion that tunnel T3 is the exit pathway of 5-OXA, and that 5-OXA tends to stay at the middle of the tunnel. The maximum errors of the calculated PMF for the predicted exit pathway, tunnel T3, were estimated by repeating the US simulations using different sets of initial positions. The PMF profile was also calculated for the transport of glycine within T3. The PMF profiles from the US simulations were in good agreement with the previous predictions that 5-OXA escape through tunnel T3 and how glycine is released to the outside of HSO was discussed.

Investigation on penetration of saffron components through lipid bilayer bound to spike protein of SARS-CoV-2 using steered molecular dynamics simulation

A coronavirus identified as COVID-19 is the reason for an infection outbreak which is started in December 2019. NO completely effective drugs and treatments are not recognized for this virus. Recently, saffron and its compounds were used to treat different viral diseases. Saffron extract and its major ingredients have shown antiviral effects. In this study, the steered molecular dynamics simulation was used for investigating the effect of four main components of saffron that include: crocin, crocetin, safranal, and picrocrocin as candidate for drug molecules, on COVID-19. The binding energies between drug molecules and spike protein and the main protease of the virus were evaluated. The obtained results based on Lennard-Jones and electrostatic potentials demonstrated that crocetin has a high affinity towards spike protein and also the main protease of the virus. Also, the quantum mechanics calculations elucidated that the crocetin could overcome energy barrier of lipid bilayer with strong dipole moment and polarizability. The pharmacokinetic and ADMET properties proved that crocetin could be a suitable drug candidate. So, crocetin could be a promising drug for treatment of COVID-19.

Efficient Estimation of Absolute Binding Free Energy for a Homeodomain-DNA Complex from Nonequilibrium Pulling Simulations

Estimation of binding free energies is one of the central aims of simulations of biomolecular complexes. We explore the accuracy and efficiency of setups based on nonequilibrium pulling simulations applied to the estimation of binding affinities of DNA-binding proteins. Absolute binding free energies are calculated over a range of temperatures and compared to results obtained previously using an equilibrium method. We show that realistic binding affinities can be obtained with the presented nonequilibrium approach, which also entails lower computational requirements. Errors of the binding free energy estimates are investigated and are shown to be comparable to those observed previously. Bounds are provided on the convergence of the errors with respect to the number of pulling simulations performed and with respect to the applied pull rate.

Atomistic Analysis of ToxN and ToxI Complex Unbinding Mechanism

ToxIN is a triangular structure formed by three protein toxins (ToxNs) and three specific noncoding RNA antitoxins (ToxIs). To respond to stimuli, ToxI is preferentially degraded, releasing the ToxN. Thus, the dynamic character is essential in the normal function interactions between ToxN and ToxI. Here, equilibrated molecular dynamics (MD) simulations were performed to study the stability of ToxN and ToxI. The results indicate that ToxI adjusts the conformation of 3' and 5' termini to bind to ToxN. Steered molecular dynamics (SMD) simulations combined with the recently developed thermodynamic integration in 3nD (TI3nD) method were carried out to investigate ToxN unbinding from the ToxIN complex. The potentials of mean force (PMFs) and atomistic pictures suggest the unbinding mechanism as follows: (1) dissociation of the 5' terminus from ToxN, (2) missing the interactions involved in the 3' terminus of ToxI without three nucleotides (G31, A32, and A33), (3) starting to unfold for ToxI, (4) leaving the binding package of ToxN for three nucleotides of ToxI, (5) unfolding of ToxI. This work provides information on the structure-function relationship at the atomistic level, which is helpful for designing new potent antibacterial drugs in the future.

Capstan-like mechanism in hyaluronan-phospholipid systems

Functionality of articular cartilage results from complex interactions between its molecular components. Among many biomolecules, two are of prime importance for lubrication: hyaluronic acid (HA) and phospholipids (PL). The purpose of this study is to discuss a mechanism of interaction between these two components and how their synergies contribute to nanobiolubrication of articular cartilage. Preliminary molecular dynamics simulations have been performed to investigate these interactions by adopting a capstan-like mechanism of action. By applying a constant pulling force to both ends of a HA molecule, wrapped around a PL micelle, we viewed the rotation of the PL micelle. The simulations were performed upon two physicochemical constraints: force- and solvent-dependency. The results show the efficiency of rotation from intermolecular bond creation and annihilation. We found a direct relation between the available surface of the micelle and the magnitude of the force, which varies significantly through the unwinding. The movement of the attached molecules is characterized by a slide-to-roll relation, which is affected by the viscosity of the surrounding medium. As a consequence, two solvents were studied for specific force conditions and the molecular dynamics simulation exhibited double the slide-to-roll coefficient for the viscous solvent as compared to its low-viscosity limit.

A Polarizable Atomic Multipole-Based Force Field for Molecular Dynamics Simulations of Anionic Lipids

In all of the classical force fields, electrostatic interaction is simply treated and explicit electronic polarizability is neglected. The condensed-phase polarization, relative to the gas-phase charge distributions, is commonly accounted for in an average way by increasing the atomic charges, which remain fixed throughout simulations. Based on the lipid polarizable force field DMPC and following the same framework as Atomic Multipole Optimized Energetics for BiomoleculAr (AMOEBA) simulation, the present effort expands the force field to new anionic lipid models, in which the new lipids contain DMPG and POPS. The parameters are compatible with the AMOEBA force field, which includes water, ions, proteins, etc. The charge distribution of each atom is represented by the permanent atomic monopole, dipole and quadrupole moments, which are derived from the ab initio gas phase calculations. Many-body polarization including the inter- and intramolecular polarization is modeled in a consistent manner with distributed atomic polarizabilities. Molecular dynamics simulations of the two aqueous DMPG and POPS membrane bilayer systems, consisting of 72 lipids with water molecules, were then carried out to validate the force field parameters. Membrane width, area per lipid, volume per lipid, deuterium order parameters, electron density profile, electrostatic potential difference between the center of the bilayer and water are all calculated, and compared with limited experimental data.

Computational studies of G protein-coupled receptor complexes: Structure and dynamics

G protein-coupled receptors (GPCRs) are ubiquitously expressed transmembrane proteins associated with a wide range of diseases such as Alzheimer's, Parkinson, schizophrenia, and also implicated in in several abnormal heart conditions. As such, this family of receptors is regarded as excellent drug targets. However, due to the high number of intracellular signaling partners, these receptors have a complex interaction networks and it becomes challenging to modulate their function. Experimentally determined structures give detailed information on the salient structural properties of these signaling complexes but they are far away from providing mechanistic insights into the underlying process. This chapter presents some of the computational tools, namely molecular dynamics, molecular docking, and molecular modeling and related analyses methods that have been used to complement experimental findings.

Modeling the human Nav1.5 sodium channel: structural and mechanistic insights of ion permeation and drug blockade

Abnormalities in the human Nav1.5 (hNav1.5) voltage-gated sodium ion channel (VGSC) are associated with a wide range of cardiac problems and diseases in humans. Current structural models of hNav1.5 are still far from complete and, consequently, their ability to study atomistic interactions of this channel is very limited. Here, we report a comprehensive atomistic model of the hNav1.5 ion channel, constructed using homology modeling technique and refined through long molecular dynamics simulations (680 ns) in the lipid membrane bilayer. Our model was comprehensively validated by using reported mutagenesis data, comparisons with previous models, and binding to a panel of known hNav1.5 blockers. The relatively long classical MD simulation was sufficient to observe a natural sodium permeation event across the channel's selectivity filters to reach the channel's central cavity, together with the identification of a unique role of the lysine residue. Electrostatic potential calculations revealed the existence of two potential binding sites for the sodium ion at the outer selectivity filters. To obtain further mechanistic insight into the permeation event from the central cavity to the intracellular region of the channel, we further employed "state-of-the-art" steered molecular dynamics (SMD) simulations. Our SMD simulations revealed two different pathways through which a sodium ion can be expelled from the channel. Further, the SMD simulations identified the key residues that are likely to control these processes. Finally, we discuss the potential binding modes of a panel of known hNav1.5 blockers to our structural model of hNav1.5. We believe that the data presented here will enhance our understanding of the structure-property relationships of the hNav1.5 ion channel and the underlying molecular mechanisms in sodium ion permeation and drug interactions. The results presented here could be useful for designing safer drugs that do not block the hNav1.5 channel.

Revealing the binding modes and the unbinding of 14-3-3σ proteins and inhibitors by computational methods

The 14-3-3σ proteins are a family of ubiquitous conserved eukaryotic regulatory molecules involved in the regulation of mitogenic signal transduction, apoptotic cell death, and cell cycle control. A lot of small-molecule inhibitors have been identified for 14-3-3 protein-protein interactions (PPIs). In this work, we carried out molecular dynamics (MD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method to study the binding mechanism between a 14-3-3σ protein and its eight inhibitors. The ranking order of our calculated binding free energies is in agreement with the experimental results. We found that the binding free energies are mainly from interactions between the phosphate group of the inhibitors and the hydrophilic residues. To improve the binding free energy of Rx group, we designed the inhibitor R9 with group R9 = 4-hydroxypheny. However, we also found that the binding free energy of inhibitor R9 is smaller than that of inhibitor R1. By further using the steer molecular dynamics (SMD) simulations, we identified a new hydrogen bond between the inhibitor R8 and residue Arg64 in the pulling paths. The information obtained from this study may be valuable for future rational design of novel inhibitors, and provide better structural understanding of inhibitor binding to 14-3-3σ proteins.

Characterizing the Free-Energy Landscape of MDM2 Protein-Ligand Interactions by Steered Molecular Dynamics Simulations

Inhibition of p53-MDM2 interaction by small molecules is considered to be a promising approach to re-activate wild-type p53 for tumor suppression. Several inhibitors of the MDM2-p53 interaction were designed and studied by the experimental methods and the molecular dynamics simulation. However, the unbinding mechanism was still unclear. The steered molecular dynamics simulations combined with Brownian dynamics fluctuation-dissipation theorem were employed to obtain the free-energy landscape of unbinding between MDM2 and their four ligands. It was shown that compounds 4 and 8 dissociate faster than compounds 5 and 7. The absolute binding free energies for these four ligands are in close agreement with experimental results. The open movement of helix II and helix IV in the MDM2 protein-binding pocket upon unbinding is also consistent with experimental MDM2-unbound conformation. We further found that different binding mechanisms among different ligands are associated with H-bond with Lys51 and Glu25. These mechanistic results may be useful for improving ligand design.

Recognition of 2',5'-linked oligoadenylates by human ribonuclease L: molecular dynamics study

The capability of current MD simulations to be used as a tool in rational design of agonists of medically interesting enzyme RNase L was tested. Dimerization and enzymatic activity of RNase L is stimulated by 2',5'-linked oligoadenylates (pA₂₅A₂₅A; 2-5A). First, it was necessary to ensure that a complex of monomeric human RNase L and 25A was stable in MD simulations. It turned out that Glu131 had to be protonated. The non-protonated Glu131 caused dissociation of 2-5A from RNase L. Because of the atypical 2'-5' internucleotide linkages and a specific spatial arrangement of the 25A trimer, when a single molecule carries all possible conformers of the glycosidic torsion angle, several versions of the AMBER force field were tested. One that best maintained functionally important interactions of 25A and RNase L was selected for subsequent MD simulations. Furthermore, we wonder whether powerful GPUs are able to produce MD trajectories long enough to convincingly demonstrate effects of subtle perturbations of interactions between 25A and RNase L. Detrimental impacts of various point mutations of RNase L (R155A, F126A, W60A, K89A) on 2-5A binding were observed on a time scale of 200 ns. Finally, 2-5A analogues with a bridged 3'--O,4'--C-alkylene linkage (B) introduced into the adenosine units (A) were used to assess ability of MD simulations to distinguish on the time scale of hundreds of nanoseconds between agonists of RNase L (pA₂₅A₂₅B, pB₂₅A₂₅A, pB₂₅A₂₅B) and inactive analogs (pA₂₅B₂₅A, pA₂₅B₂₅B, pB₂₅B₂₅A, pB₂₅B₂₅B). Agonists were potently bound to RNase L during 200 ns MD runs. For inactive 2-5A analogs, by contrast, significant disruptions of their interactions with RNase L already within 100 ns MD runs were found.

Bringing Clarity to the Prediction of Protein-Ligand Binding Free Energies via "Blurring"

We present a method to evaluate the free energies of ligand binding utilizing a Monte Carlo estimation of the configuration integrals concomitant with uncertainty quantification. Ensembles for integration are built through systematically perturbing an initial ligand conformation in a rigid binding pocket, which is optimized separately prior to incorporation of the ligand. We call the procedure producing the ensembles "blurring", and it is carried out using an in-house developed code. The Boltzmann factor contribution of each pose to the configuration integral is computed and from there the free energy is obtained. Potential function uncertainties are estimated using a fragment-based error propagation method. This method has been applied to a set of small aromatic ligands complexed with T4 Lysozyme L99A mutant. Microstate energies have been determined with the force fields ff99SB and ff94, and the semiempirical method PM6DH2 in conjunction with continuum solvation models including Generalized Born (GB), the Conductor-like Screening Model (COSMO), and SMD. Of the methods studied, PM6DH2-based scoring gave binding free energy estimates, which yielded a good correlation to the experimental binding affinities (R2 = 0.7). All methods overestimated the calculated binding affinities. We trace this to insufficient sampling, the single static protein structure, and inaccuracies in the solvent models we have used in this study.

Analyzing machupo virus-receptor binding by molecular dynamics simulations

In many biological applications, we would like to be able to computationally predict mutational effects on affinity in protein-protein interactions. However, many commonly used methods to predict these effects perform poorly in important test cases. In particular, the effects of multiple mutations, non alanine substitutions, and flexible loops are difficult to predict with available tools and protocols. We present here an existing method applied in a novel way to a new test case; we interrogate affinity differences resulting from mutations in a host-virus protein-protein interface. We use steered molecular dynamics (SMD) to computationally pull the machupo virus (MACV) spike glycoprotein (GP1) away from the human transferrin receptor (hTfR1). We then approximate affinity using the maximum applied force of separation and the area under the force-versus-distance curve. We find, even without the rigor and planning required for free energy calculations, that these quantities can provide novel biophysical insight into the GP1/hTfR1 interaction. First, with no prior knowledge of the system we can differentiate among wild type and mutant complexes. Moreover, we show that this simple SMD scheme correlates well with relative free energy differences computed via free energy perturbation. Second, although the static co-crystal structure shows two large hydrogen-bonding networks in the GP1/hTfR1 interface, our simulations indicate that one of them may not be important for tight binding. Third, one viral site known to be critical for infection may mark an important evolutionary suppressor site for infection-resistant hTfR1 mutants. Finally, our approach provides a framework to compare the effects of multiple mutations, individually and jointly, on protein-protein interactions.

Steered molecular dynamics simulations for studying protein-ligand interaction in cyclin-dependent kinase 5

In this study, we applied steered molecular dynamics (SMD) simulations to investigate the unbinding mechanism of nine inhibitors of the enzyme cyclin-dependent kinase 5 (CDK5). The study had two major objectives: (i) to create a correlation between the unbinding force profiles and the inhibition activities of these compounds expressed as IC50 values; (ii) to investigate the unbinding mechanism and to reveal atomistic insights, which could help identify accessory binding sites and transient interactions. Overall, we carried out 1.35 μs of cumulative SMD simulations. We showed that SMD could qualitatively discriminate binders from nonbinders, while it failed to properly rank series of inhibitors, particularly when IC50 values were too similar. From a mechanistic standpoint, SMD provided useful insights related to transient and dynamical interactions, which could complement static description obtained by X-ray crystallography experiments. In conclusion, the present study represents a further step toward a systematic exploitation of SMD and other dynamical approaches in structure-based drug design and computational medicinal chemistry.

Substrate recognition by norovirus polymerase: microsecond molecular dynamics study

Molecular dynamics simulations of complexes between Norwalk virus RNA dependent RNA polymerase and its natural CTP and 2dCTP (both containing the O5'-C5'-C4'-O4' sequence of atoms bridging the triphosphate and sugar moiety) or modified coCTP (C5'-O5'-C4'-O4'), cocCTP (C5'-O5'-C4'-C4'') substrates were produced by means of CUDA programmable graphical processing units and the ACEMD software package. It enabled us to gain microsecond MD trajectories clearly showing that similar nucleoside triphosphates can bind surprisingly differently into the active site of the Norwalk virus RNA dependent RNA polymerase. It corresponds to their different modes of action (CTP-substrate, 2dCTP-poor substrate, coCTP-chain terminator, cocCTP-inhibitor). Moreover, extremely rare events-as repetitive pervasion of Arg182 into a potentially reaction promoting arrangement-were captured.