Unbiased Molecular Dynamics Simulations Research Articles

Molecular Simulation to Investigate Open-Close Motion of a Flagellar Export Apparatus Protein FlhAC.

Molecular dynamics (MD) simulation and parallel cascade selection molecular dynamics (PaCS-MD) are widely used to investigate large-amplitude motions of proteins. PaCS-MD is an enhanced conformational sampling method consisting of cycles of parallel unbiased MD simulations combined with a selection of MD snapshots as the initial structures for the next cycle. In addition, free energy calculation can be achieved by the combination of PaCS-MD and the Markov state model (MSM). In this chapter, the protocols to investigate the open-close motion of a flagellar export apparatus protein, FlhAC, by MD and the combination of PaCS-MD and MSM are described.

Conformational Heterogeneity of RNA Stem-Loop Hairpins Bound to FUS-RNA Recognition Motif with Disordered RGG Tail Revealed by Unbiased Molecular Dynamics Simulations.

RNA-protein complexes use diverse binding strategies, ranging from structurally well-defined interfaces to completely disordered regions. Experimental characterization of flexible segments is challenging and can be aided by atomistic molecular dynamics (MD) simulations. Here, we used an extended set of microsecond-scale MD trajectories (400 μs in total) to study two FUS-RNA constructs previously characterized by nuclear magnetic resonance (NMR) spectroscopy. The FUS protein contains a well-structured RNA recognition motif domain followed by a presumably disordered RGG tail that binds RNA stem-loop hairpins. Our simulations not only provide several suggestions complementing the experiments but also reveal major methodological difficulties in studies of such complex RNA-protein interfaces. Despite efforts to stabilize the binding via system-specific force-field adjustments, we have observed progressive distortions of the RNA-protein interface inconsistent with experimental data. We propose that the dynamics is so rich that its converged description is not achievable even upon stabilizing the system. Still, after careful analysis of the trajectories, we have made several suggestions regarding the binding. We identify substates in the RNA loops, which can explain the NMR data. The RGG tail localized in the minor groove remains disordered, sampling countless transient interactions with the RNA. There are long-range couplings among the different elements contributing to the recognition, which can lead to allosteric communication throughout the system. Overall, the RNA-FUS systems form dynamical ensembles that cannot be fully represented by single static structures. Thus, albeit imperfect, MD simulations represent a viable tool to investigate dynamic RNA-protein complexes.

CHARMM-GUI Enhanced Sampler for various collective variables and enhanced sampling methods.

Enhanced sampling methodologies modifying underlying Hamiltonians can be used for the systems with a rugged potential energy surface that makes it hard to observe convergence using conventional unbiased molecular dynamics (MD) simulations. We present CHARMM-GUI Enhanced Sampler, a web-based tool to prepare various enhanced sampling simulations inputs with user-selected collective variables (CVs). Enhanced Sampler provides inputs for the following nine methods: accelerated MD, Gaussian accelerated MD, conformational flooding, metadynamics, adaptive biasing force, steered MD, temperature replica exchange MD, replica exchange solute tempering 2, and replica exchange umbrella sampling for the method-implemented MD packages including AMBER, CHARMM, GENESIS, GROMACS, NAMD, and OpenMM. Users only need to select a group of atoms via intuitive web-implementation in order to define commonly used nine CVs of interest: center of mass based distance, angle, dihedral, root-mean-square-distance, radius of gyration, distance projected on axis, two types of angles projected on axis, and coordination numbers. The enhanced sampling methods are tested with several biological systems to illustrate their efficiency over conventional MD. Enhanced Sampler with carefully optimized system-dependent parameters will help users to get meaningful results from their enhanced sampling simulations.

Function-Related Asymmetry of the Interactions between Matrix Loops and Conserved Sequence Motifs in the Mitochondrial ADP/ATP Carrier.

The ADP/ATP carrier (AAC) plays a central role in oxidative metabolism by exchanging ATP and ADP across the inner mitochondrial membrane. Previous experiments have shown the involvement of the matrix loops of AAC in its function, yet potential mechanisms remain largely elusive. One obstacle is the limited information on the structural dynamics of the matrix loops. In the current work, unbiased all-atom molecular dynamics (MD) simulations were carried out on c-state wild-type AAC and mutants. Our results reveal that: (1) two ends of a matrix loop are tethered through interactions between the residue of triplet 38 (Q38, D143 and Q240) located at the C-end of the odd-numbered helix and residues of the [YF]xG motif located before the N-end of the short matrix helix in the same domain; (2) the initial progression direction of a matrix loop is determined by interactions between the negatively charged residue of the [DE]G motif located at the C-end of the short matrix helix and the capping arginine (R30, R139 and R236) in the previous domain; (3) the two chemically similar residues D and E in the highly conserved [DE]G motif are actually quite different; (4) the N-end of the M3 loop is clamped by the [DE]G motif and the capping arginine of domain 2 from the two sides, which strengthens interactions between domain 2 and domain 3; and (5) a highly asymmetric stable core exists within domains 2 and 3 at the m-gate level. Moreover, our results help explain almost all extremely conserved residues within the matrix loops of the ADP/ATP carriers from a structural point of view. Taken together, the current work highlights asymmetry in the three matrix loops and implies a close relationship between asymmetry and ADP/ATP transport.

Partial Dissociation of Antigenic Peptides from MHC I -- Linking NMR Data to Microstates Observed in MD Simulations

Although several crystal structures of Major Histocompatibility Complex I (MHC I) binding an antigenic peptide (pMHC I) have been solved, the structural dynamics of pMHC I at the cell surface remain largely elusive. According to a recent NMR study on HLA-B*35:01 (Yanaka et al., J. Biol. Chem., 2014), pMHC I complexes can adopt a minor state in which the antigenic peptide is bound tightly to MHC I and a major state in which the peptide is bound more loosely. In unbiased molecular dynamics (MD) simulations of the pMHC I studied by Yanaka et al., the peptide N-terminus dissociated from the MHC I binding groove during a few hundred nanoseconds. This finding suggests that pMHC I with completely bound peptide may correspond to the proposed minor state, whereas pMHC I with partially dissociated peptide may constitute the major state. To underpin this, the potential of mean force (PMF) along the distance between the anchor residue of the peptide N-terminus and its binding partner in the MHC I binding groove has been calculated. The resulting free energy differences can be compared to values derived from the published NMR data. Preliminary results from umbrella sampling (1 microsecond of MD simulation per umbrella window) suggest that pMHC I in which the peptide N-terminus has dissociated from MHC I are lowest in free energy. Contrarily, Hamiltonian replica exchange simulations (2 ×2 microseconds), in which the effective temperature of the MHC I binding groove and the antigenic peptide is increased, predict pMHC I with completely bound peptide to be the global minimum of the free energy. This seemingly contradicting difference is discussed and further elucidated.

Multiagent Reinforcement Learning-Based Adaptive Sampling for Conformational Dynamics of Proteins.

Machine learning is increasingly applied to improve the efficiency and accuracy of molecular dynamics (MD) simulations. Although the growth of distributed computer clusters has allowed researchers to obtain higher amounts of data, unbiased MD simulations have difficulty sampling rare states, even under massively parallel adaptive sampling schemes. To address this issue, several algorithms inspired by reinforcement learning (RL) have arisen to promote exploration of the slow collective variables (CVs) of complex systems. Nonetheless, most of these algorithms are not well-suited to leverage the information gained by simultaneously sampling a system from different initial states (e.g., a protein in different conformations associated with distinct functional states). To fill this gap, we propose two algorithms inspired by multiagent RL that extend the functionality of closely related techniques (REAP and TSLC) to situations where the sampling can be accelerated by learning from different regions of the energy landscape through coordinated agents. Essentially, the algorithms work by remembering which agent discovered each conformation and sharing this information with others at the action-space discretization step. A stakes function is introduced to modulate how different agents sense rewards from discovered states of the system. The consequences are three-fold: (i) agents learn to prioritize CVs using only relevant data, (ii) redundant exploration is reduced, and (iii) agents that obtain higher stakes are assigned more actions. We compare our algorithm with other adaptive sampling techniques (least counts, REAP, TSLC, and AdaptiveBandit) to show and rationalize the gain in performance.

Molecular simulations of fluoxetine in hydrated lipid bilayers, as well as in aqueous solutions containing β-cyclodextrin

Fluoxetine, which is a well-known antidepressant drug, is studied in hydrated cholesterol-free and cholesterol-containing lipid bilayers through unbiased and biased atomistic molecular dynamics simulations. The latter are conducted for the calculation of the potential of mean force (PMF) of fluoxetine along an axis perpendicular to the two leaflets of the bilayer. The PMF indicates that the drug prefers to reside inside the lipid phase and allows us to calculate important thermodynamic properties, such as the Gibbs energy difference of partitioning from the water to the lipid phase and the Gibbs energy barrier for hopping events between the two leaflets of the bilayer. The results from the biased simulations are in accord with the mass density profiles calculated from the unbiased simulations. Moreover, we estimate the effect of fluoxetine mole fraction on the order parameters of the lipid alkyl chains and on the area per lipid. It is also found that fluoxetine forms a hydrogen bond network with lipids and water molecules penetrating into the lipid phase. In addition, fluoxoetine is studied in detail in aqueous solutions containing β-cyclodextrin. It is observed from unbiased molecular dynamics simulations that the two aforementioned molecules form a noncovalent complex spontaneously and the calculated binding free energy is in agreement with the literature.

Interacting binding insights and conformational consequences of the differential activity of cannabidiol with two endocannabinoid-activated G-protein-coupled receptors

Cannabidiol (CBD), the major non-psychoactive phytocannabinoid present in the plant Cannabis sativa, has displayed beneficial pharmacological effects in the treatment of several neurological disorders including, epilepsy, Parkinson’s disease, and Alzheimer’s disease. In particular, CBD is able to modulate different receptors in the endocannabinoid system, some of which belong to the family of G-protein-coupled receptors (GPCRs). Notably, while CBD is able to antagonize some GPCRs in the endocannabinoid system, it also seems to activate others. The details of this dual contrasting functional feature of CBD, that is, displaying antagonistic and (possible) agonistic ligand properties in related receptors, remain unknown. Here, using computational methods, we investigate the interacting determinants of CBD in two closely related endocannabinoid-activated GPCRs, the G-protein-coupled receptor 55 (GPR55) and the cannabinoid type 1 receptor (CB1). While in the former, CBD has been demonstrated to function as an antagonist, the way by which CBD modulates the CB1 receptor remains unclear. Namely, CBD has been suggested to directly trigger receptor’s activation, stabilize CB1 inactive conformations or function as an allosteric modulator. From microsecond-length unbiased molecular dynamics simulations, we found that the presence of the CBD ligand in the GPR55 receptor elicit conformational changes associated with antagonist-bound GPCRs. In contrast, when the GPR55 receptor is simulated in complex with the selective agonist ML186, agonist-like conformations are sampled. These results are in agreement with the proposed modulatory function of each ligand, showing that the computational techniques utilized to characterize the GPR55 complexes correctly differentiate the agonist-bound and antagonist-bound systems. Prompted by these results, we investigated the role of the CBD compound on the CB1 receptor using similar computational approaches. The all-atom MD simulations reveal that CBD induces conformational changes linked with agonist-bound GPCRs. To contextualize the results we looked into the CB1 receptor in complex with a well-established antagonist. In contrast to the CBD/CB1 complex, when the CB1 receptor is simulated in complex with the ligand antagonist AM251, inactive conformations are explored, showing that the computational techniques utilized to characterize the CB1 complexes correctly differentiate the agonist-bound and antagonist-bound systems. In addition, our results suggest a previously unknown sodium-binding site located in the extracellular domain of the CB1 receptor. From our detailed characterization, we found particular interacting loci in the binding sites of the GPR55 and the CB1 receptors that seem to be responsible for the differential functional features of CBD. Our work will pave the way for understanding the CBD pharmacology at a molecular level and aid in harnessing its potential therapeutic use.

Counteractive Effects of Choline Geranate (CAGE) ILs and Ethanol on Insulin’s Stability—A Leap Forward towards Oral Insulin Formulation

Choline geranate (CAGE) ionic liquids (ILs) stabilize insulin, thereby aiding its oral delivery, whereas ethanol (EtOH) affects its stability by disrupting the hydrophobic interactions. In this study, cognizance of the stabilization mechanism of insulin dimer in the presence of both CAGE ILs and EtOH mixtures is achieved through biased and unbiased molecular dynamics (MD) simulations. Here, two order parameters are employed to study the insulin dimer dissociation using well-tempered metadynamics (WT-MetaD). The stability of insulin is found to be strongly maintained until a 0.20 mole fraction of EtOH. Besides, higher concentrations of EtOH marginally affect the insulin stability. Moreover, geranate anions form a higher number of H-bonding interactions with water molecules, which aids insulin stabilization. Conversely, the addition of EtOH minimizes the water-mediated H-bonding interactions of geranate. Additionally, geranate traps the EtOH molecules, thereby preventing the interactions between insulin and EtOH. Furthermore, the free energy landscape (FEL) reveals the absence of dimer dissociation along with noticeable deviations in the distances R and the number of contacts Q. The dimerization free energy of insulin was calculated to be −16.1 kcal/mol at a 0.20 mole fraction of EtOH. Moreover, increments in mole fractions of EtOH effectuate a decrease in the insulin stability. Thus, the present study represents CAGE ILs as efficient insulin dimer stabilizes at low concentrations of EtOH.

Comparative study of the unbinding process of some HTLV-1 protease inhibitors using unbiased molecular dynamics simulations.

The HTLV-1 protease is one of the major antiviral targets to overwhelm this virus. Several research groups have developed protease inhibitors, but none has been successful. In this regard, developing new HTLV-1 protease inhibitors to fix the defects in previous inhibitors may overcome the lack of curative treatment for this oncovirus. Thus, we decided to study the unbinding pathways of the most potent (compound 10, PDB ID 4YDF, Ki = 15 nM) and one of the weakest (compound 9, PDB ID 4YDG, Ki = 7900 nM) protease inhibitors, which are very structurally similar. We conducted 12 successful short and long simulations (totaling 14.8 μs) to unbind the compounds from two monoprotonated (mp) forms of protease using the Supervised Molecular Dynamics (SuMD) without applying any biasing force. The results revealed that Asp32 or Asp32' in the two forms of mp state similarly exert powerful effects on maintaining both potent and weak inhibitors in the binding pocket of HTLV-1 protease. In the potent inhibitor's unbinding process, His66' was a great supporter that was absent in the weak inhibitor's unbinding pathway. In contrast, in the weak inhibitor's unbinding process, Trp98/Trp98' by pi-pi stacking interactions were unfavorable for the stability of the inhibitor in the binding site. In our opinion, these results will assist in designing more potent and effective inhibitors for the HTLV-1 protease.

Transient Hoogsteen Base Pairs Observed in Unbiased Molecular Dynamics Simulations of DNA.

Duplex DNA is modeled as canonical B-DNA displaying the characteristic Watson-Crick base pairs. A less common and short-lived pairing of the nucleobases is the Hoogsteen (HG) conformation. The low population of the HG base pairs (<1%) necessitates extended sampling times in order to analyze through unbiased molecular dynamics (MD) simulations. We have discovered that with extended sampling times using multiple independent copies of an 18-mer sequence, the MD trajectories reproduce the expected and transient HG base pairing. Consistent with experiment, the percentage of the HG events are within the range of ∼0.1-1.0% over the combined aggregate sampling time of more than 3.6 ms. We present the reliability of the current AMBER set of nucleic acid force fields and tools to accurately simulate naturally occurring base-pairing and opening events without any bias or restraints. The mechanism consists of base pair fraying, flipping of the purine, and reformation with HG base pairs.

Enhanced Expansion and Reduced Kiss-and-Run Events in Fusion Pores Steered by Synaptotagmin-1 C2B Domains

The fusion pore controls the release of exocytotic vesicle contents through a precise orchestration of lipids from the fusing membranes and proteins. There is a major lipid reorganization during the different stages in life of the fusion pore (membrane fusion, nucleation, and expansion) that can be scrutinized thermodynamically. In this work, using umbrella sampling simulations we describe the expansion of the fusion pore. We have calculated free energy profiles to drive a nascent, just nucleated, fusion pore to its expanded configuration. We have quantified the effects on the free energy of one and two Synaptotagmin-1 C2B domains in the cytosolic space. We show that C2B domains cumulatively reduce the cost for expansion, favoring the system to evolve toward full fusion. Finally, by conducting thousands of unbiased molecular dynamics simulations, we show that C2B domains significantly decrease the probability of kiss-and-run events.

Unveiling the Releasing Processes of Pt(II)-Based Anticancer Drugs from Oxidized Carbon Nanohorn: An In Silico Study.

About half of all cancer chemotherapies currently applied involve medication with the three worldwide approved Pt(II)-based drugs, cisplatin (cddp), carboplatin (cpx), and oxaliplatin (oxa), due to their notable antitumor activity for several cancers. However, this wide application is accompanied by severe side effects, such as nephrotoxicity, myelosuppression, and neurotoxicity, as a result of their low bioavailability and selectivity for cancer cells. To mitigate these drawbacks, the use of chemically functionalized carbon nanohorns (CNH) as nanocarriers represents a potential formulation since CNH has been noted for their biodegradability, biocompatibility, low toxicity, and cavities dimensionally compatible with small drugs. This work reports energetic and dynamic analyses of complexes formed by oxidized CNH (CNHox) and the cddp, cpx, and oxa drugs. Using unbiased molecular dynamics (MD) simulations, we show that the encapsulated formulations (cddp@CNHox, cpx@CNHox, and oxa@CNHox) were more stable by ∼11.0 kcal mol-1 than the adsorbed ones (cddp > CNHox, cpx > CNHox, and oxa > CNHox). This high stability, mainly governed by van der Waals interactions, was responsible for the drug confinement during the entire simulation time (200 ns). The biased MD simulations of the inclusion complexes confirmed the nonspontaneity of the drug release since the potentials of mean force (PMF) indicated the endergonic character of this process. Additionally, the releasing energy profiles pointed out that the free energy barrier (ΔΔG≠) for the escape from CNHox cavity follows the order oxa > cpx ∼ cddp, with the value for the oxa complex (21-26 kcal mol-1) found to be about 36 and 30% larger than those for cpx and cddp, respectively. While the approximate residence time (tres) of the oxa drug inside the CNHox cavity was 5.45 × 108 s, the same measure for the cddp and cpx drugs was 5.3 × 105 and 1.60 × 103 s. Simulations also revealed that the escape of oxa with the oxalate group facing the nanowindow was the most unfavorable process, giving tres = 1.09 × 109 s. Besides reinforcing and extending the nanovectorization of cddp, cpx, and oxa in CNHox for cancer chemotherapies, all features considered may provide interpretations for experimental data and encourage new investigations aiming to propose less aggressive treatments for oncological diseases.

Expanding the Accessible Chemical Space of SIRT2 Inhibitors through Exploration of Binding Pocket Dynamics.

Considerations of binding pocket dynamics are one of the crucial aspects of the rational design of binders. Identification of alternative conformational states or cryptic subpockets could lead to the discovery of completely novel groups of the ligands. However, experimental characterization of pocket dynamics, besides being expensive, may not be able to elucidate all of the conformational states relevant for drug discovery projects. In this study, we propose the protocol for computational simulations of sirtuin 2 (SIRT2) binding pocket dynamics and its integration into the structure-based virtual screening (SBVS) pipeline. Initially, unbiased molecular dynamics simulations of SIRT2:inhibitor complexes were performed using optimized force field parameters of SIRT2 inhibitors. Time-lagged independent component analysis (tICA) was used to design pocket-related collective variables (CVs) for enhanced sampling of SIRT2 pocket dynamics. Metadynamics simulations in the tICA eigenvector space revealed alternative conformational states of the SIRT2 binding pocket and the existence of a cryptic subpocket. Newly identified SIRT2 conformational states outperformed experimentally resolved states in retrospective SBVS validation. After performing prospective SBVS, compounds from the under-represented portions of the SIRT2 inhibitor chemical space were selected for in vitro evaluation. Two compounds, NDJ18 and NDJ85, were identified as potent and selective SIRT2 inhibitors, which validated the in silico protocol and opened up the possibility for generalization and broadening of its application. The anticancer effects of the most potent compound NDJ18 were examined on the triple-negative breast cancer cell line. Results indicated that NDJ18 represents a promising structure suitable for further evaluation.

Solvent's Role in Cavity-Ligand Recognition Would Depend on the Mode of Ligand Diffusion.

Solvent is known to play crucial roles in dictating the thermodynamics and kinetics of the biomolecular recognition process. Here, we show that the extent of significance of water in modulating the ligand recognition process is critically contingent on the ligand diffusion and on the constraints introduced on it. Toward the end, we use a well-known prototypical system of spherical ligand diffusing freely toward a hydrophobic concave cavity in explicit water. We analyze a large series of adaptively sampled unbiased molecular dynamics simulation trajectories within the framework of time-structured independent component analysis (TICA). Our quantitative investigations reveal that water would play a significant role in the ligand recognition process, provided that the ligand is constricted to diffuse along a centro-symmetric fashion. On the contrary, water's contribution in the ligand recognition process would diminish to a negligible value if the ligand freely diffuses toward the pocket. A Markov state model (MSM) constructed using the simulated trajectories identifies a set of transiently populated metastable states comprising partially ligand-unbound macro states, alongside ligand-bound and ligand-unbound pose and gives rise to multiple transition paths of ligand in its way to the hydrophobic cavity. Lifting the restriction on ligand movement changes its binding pathway, time scales, and the extent of the role of solvent in modulating the recognition process.

Ligand-Dependent Modulation of the Dynamics of Intracellular Loops Dictates Functional Selectivity of 5-HT2AR

The serotonin 2A receptor (5-HT2AR) subtype of the G protein-coupled receptor (GPCR) family is involved in a plethora of neuromodulatory functions (e.g., neurogenesis, sleep, and cognitive processes). 5-HT2AR is the target of pharmacologically distinct classes of ligands, binding of which either activate or inactivate the receptor. Although high-resolution structures of 5-HT2AR as well as several other 5-HT GPCRs provided snapshots of both active and inactive conformational states, these structures, representing a truncated form of the receptor, cannot fully explain the mechanism of conformational transitions during their function. Importantly, biochemical studies have suggested the importance of intracellular loops in receptor functions. In our previous study, a model of the ligand-free form of 5-HT2AR with the third intracellular loop (ICL3) has been meticulously built. Here, we have investigated the functional regulation of 5-HT2AR with intact intracellular loops in ligand-free and five distinct ligand-bound configurations using unbiased atomistic molecular dynamics (MD) simulations. The selected ligands belong to either of the full, partial, or inverse agonist classes, which exert distinct pharmacological responses. We have observed significant structural, dynamic, and thermodynamic differences within ligand-bound complexes. Our results revealed, for the first time, that either activation or inactivation of the receptor upon specific ligand binding is primarily achieved through conformational transitions of its second and third intracellular loops (ICL2 and ICL3). A remarkable allosteric cross-talk between the ligand-binding site and the distal intracellular parts of the receptor, where binding of a specific ligand thermodynamically controls (either stabilizes or destabilizes) the intracellular region, consisting of crucial dynamic elements ICL2 and ICL3, and differential conformational transitions of these loops determine ligand-dependent functional selectivity.

Signatures of sluggish dynamics and local structural ordering during ice nucleation.

We investigate the microscopic pathway of spontaneous crystallization in the ST2 model of water under deeply supercooled conditions via unbiased classical molecular dynamics simulations. After quenching below the liquid-liquid critical point, the ST2 model spontaneously separates into low-density liquid (LDL) and high-density liquid phases, respectively. The LDL phase, which is characterized by lower molecular mobility and enhanced structural order, fosters the formation of a sub-critical ice nucleus that, after a stabilization time, develops into the critical nucleus and grows. Polymorphic selection coincides with the development of the sub-critical nucleus and favors the formation of cubic (Ic) over hexagonal (Ih) ice. We rationalize polymorphic selection in terms of geometric arguments based on differences in the symmetry of second neighbor shells of ice Ic and Ih, which are posited to favor formation of the former. The rapidly growing critical nucleus absorbs both Ic and Ih crystallites dispersed in the liquid phase, a crystal with stacking faults. Our results are consistent with, and expand upon, recent observations of non-classical nucleation pathways in several systems.

Reconstruction of the unbinding pathways of noncovalent SARS-CoV and SARS-CoV-2 3CLpro inhibitors using unbiased molecular dynamics simulations.

The main protease (3CLpro) is one of the essential components of the SARS-CoVs viral life cycle, which makes it an interesting target for overpowering these viruses. Although many covalent and noncovalent inhibitors have been designed to inhibit this molecular target, none have gained FDA approval as a drug. Because of the high rate of COVID-19 pandemic development, in addition to laboratory research, we require in silico methods to accelerate rational drug design. The unbinding pathways of two SARS-CoV and SARS-CoV-2 3CLpro noncovalent inhibitors with the PDB IDs: 3V3M, 4MDS, 6W63, 5RF7 were explored from a comparative perspective using unbiased molecular dynamics (UMD) simulations. We uncovered common weak points for selected inhibitors that could not interact significantly with a binding pocket at specific residues by all their fragments. So water molecules entered the free binding S regions and weakened protein-inhibitor fundamental interactions gradually. N142, G143, and H163 are the essential residues, which cause key protein-ligand interactions in the binding pocket. We believe that these results will help design new potent inhibitors against SARS-CoV-2.

A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls.

The string method with swarms of trajectories (SMwST) is an algorithm that identifies a physically meaningful transition pathway─a one-dimensional curve, embedded within a high-dimensional space of selected collective variables. The SMwST algorithm leans on a series of short, unbiased molecular dynamics simulations spawned at different locations of the discretized path, from whence an average dynamic drift is determined to evolve the string toward an optimal pathway. However conceptually simple in both its theoretical formulation and practical implementation, the SMwST algorithm is computationally intensive and requires a careful choice of parameters for optimal cost-effectiveness in applications to challenging problems in chemistry and biology. In this contribution, the SMwST algorithm is presented in a self-contained manner, discussing with a critical eye its theoretical underpinnings, applicability, inherent limitations, and use in the context of path-following free-energy calculations and their possible extension to kinetics modeling. Through multiple simulations of a prototypical polypeptide, combining the search of the transition pathway and the computation of the potential of mean force along it, several practical aspects of the methodology are examined with the objective of optimizing the computational effort, yet without sacrificing accuracy. In light of the results reported here, we propose some general guidelines aimed at improving the efficiency and reliability of the computed pathways and free-energy profiles underlying the conformational transitions at hand.

Stem-loop II motif (s2m) structural and dynamical differences between SARS-CoV-2 and b.1.617.2 (Delta) variant

The stem-loop II motif (s2m) present in the 3'-untranslated region of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) genome differs from the B.1.617.2 (Delta) variant by one mutation (G15U). Preliminary bioinformatic and experimental results suggest a relationship between the Delta variant and the s2m G15U mutation. In addition, the G15U s2m secondary structure is unclear compared to known SARS-CoV-2 s2m base pairing, which could influence suspected viral lifecycle events. To define structural and dynamic differences, unbiased molecular dynamics simulations of the G15U s2m were carried out for 3.5 µs, using starting coordinates from our previously reported SARS-CoV-2 s2m model based on NMR NOE data. Simulations were carried out at 283 K, to match to NMR data, and 310 K, to match physiological conditions. Results from the 3.5 µs trajectories show significant differences in three-dimensional structure and secondary structure resulting in key tertiary motif changes through hydrogen bonding and base-stacking interactions in the terminal loop, stem, and bulges. In this work, structural and dynamical differences are defined between the wildtype SARS-CoV-2 and G15U s2m, enabling further analysis of the effects on possible lifecycle events and identification of new pharmacological antiviral targets. The stem-loop II motif (s2m) present in the 3'-untranslated region of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) genome differs from the B.1.617.2 (Delta) variant by one mutation (G15U). Preliminary bioinformatic and experimental results suggest a relationship between the Delta variant and the s2m G15U mutation. In addition, the G15U s2m secondary structure is unclear compared to known SARS-CoV-2 s2m base pairing, which could influence suspected viral lifecycle events. To define structural and dynamic differences, unbiased molecular dynamics simulations of the G15U s2m were carried out for 3.5 µs, using starting coordinates from our previously reported SARS-CoV-2 s2m model based on NMR NOE data. Simulations were carried out at 283 K, to match to NMR data, and 310 K, to match physiological conditions. Results from the 3.5 µs trajectories show significant differences in three-dimensional structure and secondary structure resulting in key tertiary motif changes through hydrogen bonding and base-stacking interactions in the terminal loop, stem, and bulges. In this work, structural and dynamical differences are defined between the wildtype SARS-CoV-2 and G15U s2m, enabling further analysis of the effects on possible lifecycle events and identification of new pharmacological antiviral targets.