Steered Molecular Dynamics

Assessing the persistence of chalcogen bonds in solution with neural network potentials

Non-covalent bonding patterns are commonly harvested as a design principle in the field of catalysis, supramolecular chemistry and functional materials to name a few. Yet, their computational description generally neglects finite temperature and environment effects, which promote competing interactions and alter their static gas-phase properties. Recently, neural network potentials (NNPs) trained on Density Functional Theory (DFT) data have become increasingly popular to simulate molecular phenomena in condensed phase with an accuracy comparable to ab initio methods. To date, most applications have centered on solid-state materials or fairly simple molecules made of a limited number of elements. Herein, we focus on the persistence and strength of chalcogen bonds involving a benzotelluradiazole in condensed phase. While the tellurium-containing heteroaromatic molecules are known to exhibit pronounced interactions with anions and lone pairs of different atoms, the relevance of competing intermolecular interactions, notably with the solvent, is complicated to monitor experimentally but also challenging to model at an accurate electronic structure level. Here, we train direct and baselined NNPs to reproduce hybrid DFT energies and forces in order to identify what are the most prevalent non-covalent interactions occurring in a solute-Cl$^-$-THF mixture. The simulations in explicit solvent highlight competition with chalcogen bonds formed with the solvent and the short-range directionality of the interaction with direct consequences for the molecular properties in the solution. The comparison with other potentials (e.g., AMOEBA, direct NNP and continuum solvent model) also demonstrates that baselined NNPs offer a reliable picture of the non-covalent interaction interplay occurring in solution.

In the footsteps of an inhibitor unbinding from the active site of human carbonic anhydrase II

The crystal structure of human carbonic anhydrase (HCA) II bound to an inhibitor molecule, 6-hydroxy-2-thioxocoumarin (FC5), shows FC5 to be located in a hydrophobic pocket at the active site. The present work employs classical molecular dynamics (MD) simulation to follow the FC5 molecule for 1 μs as it unbinds from its binding location, adopts the path of substrate/product diffusion (path 1) to leave the active site at around 75 ns. It is then found to undergo repeated binding and unbinding at different locations on the surface of the enzyme in water. Several transient excursions through different regions of the enzyme are also observed prior to its exit from the active site. These transient paths are combined with functionally relevant cavities/channels to enlist five additional pathways (path 2-6). Pathways 1-6 are subsequently explored using steered MD and umbrella sampling simulations. A free energy barrier of 0.969 kcal mol^-1 is encountered along path 1, while barriers in the range of 0.57-2.84 kcal mol^-1 are obtained along paths 2, 4 and 5. We also analyze in detail the interaction between FC5 and the enzyme along each path as the former leaves the active site of HCA II. Our results indicate path 1 to be the major exit pathway for FC5, although competing contributions may also come from the paths 2, 4 and 5.Communicated by Ramaswamy H. Sarma.

Molecular dynamics simulations suggest possible activation and deactivation pathways in the hERG channel

The elusive activation/deactivation mechanism of hERG is investigated, a voltage-gated potassium channel involved in severe inherited and drug-induced cardiac channelopathies, including the Long QT Syndrome. Firstly, the available structural data are integrated by providing a homology model for the closed state of the channel. Secondly, molecular dynamics combined with a network analysis revealed two distinct pathways coupling the voltage sensor domain with the pore domain. Interestingly, some LQTS-related mutations known to impair the activation/deactivation mechanism are distributed along the identified pathways, which thus suggests a microscopic interpretation of their role. Split channels simulations clarify a surprising feature of this channel, which is still able to gate when a cut is introduced between the voltage sensor domain and the neighboring helix S5. In summary, the presented results suggest possible activation/deactivation mechanisms of non-domain-swapped potassium channels that may aid in biomedical applications.

Costa et al. present the electro-mechanical coupling between the voltage sensor and pore domain of the hERG channel using a combination of molecular dynamics simulations and theoretical network analyses.

Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased “downhill” trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand–protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Computational methods and theory for ion channel research

Ion channels are fundamental biological devices that act as gates in order to ensure selective ion transport across cellular membranes; their operation constitutes the molecular mechanism through which basic biological functions, such as nerve signal transmission and muscle contraction, are carried out. Here, we review recent results in the field of computational research on ion channels, covering theoretical advances, state-of-the-art simulation approaches, and frontline modeling techniques. We also report on few selected applications of continuum and atomistic methods to characterize the mechanisms of permeation, selectivity, and gating in biological and model channels.

Assessing the Stability of Biological Fibrils by Molecular-Scale Simulations

The nanomechanical characterization of several biological fibrils that are the result of protein aggregation via molecular dynamics simulation is nowadays feasible, and together with atomic force microscopy experiments has widened our understanding of the forces in the regime of pN-nN and system sizes of about hundreds of nanometers. Several methodologies have been developed to achieve this target, and they range from the atomistic representation via molecular force fields to coarse-grained strategies that provide comparable results with experiments in a systematic way. In this chapter, we discuss several methodologies for the calculation of mechanical parameters, such as the elastic constants of relevant biological systems. They are presented together with details about parameterization and current limitations. Then, we discuss some of the applications of such methodologies for the description of bacterial filament and β-amyloid systems. Finally, the latest lines of development are discussed.

Correlation in Domain Fluctuations Navigates Target Search of a Viral Peptide along RNA

Biological macromolecules often exhibit correlations in fluctuations involving distinct domains. This study decodes their functional implications in RNA-protein recognition and target-specific binding. The target search of a peptide along RNA in a viral TAR-Tat complex is closely monitored using atomistic simulations, steered molecular dynamics simulations, free energy calculations, and a machine-learning-based clustering technique. An anticorrelated domain fluctuation is identified between the tetraloop and the bulge region in the apo form of TAR RNA that sets a hierarchy in the domain-specific fluctuations at each binding event and that directs the succeeding binding footsteps. Thus, at each binding footstep, the dynamic partner selects an RNA location for binding where it senses a higher fluctuation, which is conventionally reduced upon binding. This event stimulates an alternate domain fluctuation, which then dictates sequential binding footstep/s and thus the search progresses. Our cross-correlation maps show that the fluctuations relay from one domain to another specific domain until the anticorrelation between those interdomain fluctuations sustains. Artificial attenuation of that hierarchical domain fluctuation inhibits specific RNA binding. The binding is completed with the arrival of a few long-lived water molecules that mediate slightly distant RNA-protein sites and finally stabilize the overall complex. The study underscores the functional importance of naturally designed fluctuating RNA motifs (bulge, tetraloop) and their interplay in dictating the directionality of the search in a highly dynamic environment.

Substrate Access Mechanism in a Novel Membrane-Bound Phospholipase A of Pseudomonas aeruginosa Concordant with Specificity and Regioselectivity

PlaF is a cytoplasmic membrane-bound phospholipase A₁ from Pseudomonas aeruginosa that alters the membrane glycerophospholipid (GPL) composition and fosters the virulence of this human pathogen. PlaF activity is regulated by a dimer-to-monomer transition followed by tilting of the monomer in the membrane. However, how substrates reach the active site and how the characteristics of the active site tunnels determine the activity, specificity, and regioselectivity of PlaF for natural GPL substrates have remained elusive. Here, we combined unbiased and biased all-atom molecular dynamics (MD) simulations and configurational free-energy computations to identify access pathways of GPL substrates to the catalytic center of PlaF. Our results map out a distinct tunnel through which substrates access the catalytic center. PlaF variants with bulky tryptophan residues in this tunnel revealed decreased catalysis rates due to tunnel blockage. The MD simulations suggest that GPLs preferably enter the active site with the sn-1 acyl chain first, which agrees with the experimentally demonstrated PLA₁ activity of PlaF. We propose that the acyl chain-length specificity of PlaF is determined by the structural features of the access tunnel, which results in favorable free energy of binding of medium-chain GPLs. The suggested egress route conveys fatty acid (FA) products to the dimerization interface and, thus, contributes to understanding the product feedback regulation of PlaF by FA-triggered dimerization. These findings open up opportunities for developing potential PlaF inhibitors, which may act as antibiotics against P. aeruginosa.

Application of in vitro PAMPA technique and in silico computational methods for blood-brain barrier permeability prediction of novel CNS drug candidates

Permeability assessment of small molecules through the blood-brain barrier (BBB) plays a significant role in the development of effective central nervous system (CNS) drug candidates. Since in vivo methods for BBB permeability estimation require a lot of time and resources, in silico and in vitro approaches are becoming increasingly popular nowadays for faster and more economical predictions in early phases of drug discovery. In this work, through application of in vitro parallel artificial membrane permeability assay (PAMPA-BBB) and in silico computational methods we aimed to examine the passive permeability of eighteen compounds, which affect serotonin and dopamine levels in the CNS. The data set was consisted of novel six human dopamine transporter (hDAT) substrates that were previously identified as the most promising lead compounds for further optimisation to achieve neuroprotective effect, twelve approved CNS drugs, and their related compounds. Firstly, PAMPA methods was used to experimentally determine effective BBB permeability (P_e) for all studied compounds and obtained results were further submitted for quantitative structure permeability relationship (QSPR) analysis. QSPR models were built by using three different statistical methods: stepwise multiple linear regression (MLR), partial least square (PLS), and support-vector machine (SVM), while their predictive capability was tested through internal and external validation. Obtained statistical parameters (MLR- R²_pred=-0.10; PLS- R²_pred=0.64, r²_m=0.69, r^/2_m=0.44; SVM- R²_pred=0.57, r²_m=0.72, r^/2_m=0.55) indicated that the SVM model is superior over others. The most important molecular descriptors (H₀p and SolvEMt_3D) were identified and used to propose structural modifications of the examined compounds in order to improve their BBB permeability. Moreover, steered molecular dynamics (SMD) simulation was employed to comprehensively investigate the permeability pathway of compounds through a lipid bilayer. Taken together, the created QSPR model could be used as a reliable and fast pre-screening tool for BBB permeability prediction of structurally related CNS compounds, while performed MD simulations provide a good foundation for future in silico examination.

Accelerating Reaction Rates of Biomolecules by Using Shear Stress in Artificial Capillary Systems

Biomimetics is a design principle within chemistry, biology, and engineering, but chemistry biomimetic approaches have been generally limited to emulating nature's chemical toolkit while emulation of nature's physical toolkit has remained largely unexplored. To begin to explore this, we designed biophysically mimetic microfluidic reactors with characteristic length scales and shear stresses observed within capillaries. We modeled the effect of shear with molecular dynamics studies and showed that this induces specific normally buried residues to become solvent accessible. We then showed using kinetics experiments that rates of reaction of these specific residues in fact increase in a shear-dependent fashion. We applied our results in the creation of a new microfluidic approach for the multidimensional study of cysteine biomarkers. Finally, we used our approach to establish dissociation of the therapeutic antibody trastuzumab in a reducing environment. Our results have implications for the efficacy of existing therapeutic antibodies in blood plasma as well as suggesting in general that biophysically mimetic chemistry is exploited in biology and should be explored as a research area.

Targeting N-Terminal Human Maltase-Glucoamylase to Unravel Possible Inhibitors Using Molecular Docking, Molecular Dynamics Simulations, and Adaptive Steered Molecular Dynamics Simulations

There are multiple drugs for the treatment of type 2 diabetes, including traditional sulfonylureas biguanides, glinides, thiazolidinediones, α-glucosidase inhibitors, glucagon-like peptide-1 (GLP-1) receptor agonists, dipeptidyl peptidase IV (DPP-4) inhibitors, and sodium-glucose cotransporter 2 (SGLT2) inhibitors. α-Glucosidase inhibitors have been used to control postprandial glucose levels caused by type 2 diabetes since 1990. α-Glucosidases are rather crucial in the human metabolic system and are principally found in families 13 and 31. Maltase-glucoamylase (MGAM) belongs to glycoside hydrolase family 31. The main function of MGAM is to digest terminal starch products left after the enzymatic action of α-amylase; hence, MGAM becomes an efficient drug target for insulin resistance. In order to explore the conformational changes in the active pocket and unbinding pathway for NtMGAM, molecular dynamics (MD) simulations and adaptive steered molecular dynamics (ASMD) simulations were performed for two NtMGAM-inhibitor [de-O-sulfonated kotalanol (DSK) and acarbose] complexes. MD simulations indicated that DSK bound to NtMGAM may influence two domains (inserted loop 1 and inserted loop 2) by interfering with the spiralization of residue 497–499. The flexibility of inserted loop 1 and inserted loop 2 can influence the volume of the active pocket of NtMGAM, which can affect the binding progress for DSK to NtMGAM. ASMD simulations showed that compared to acarbose, DSK escaped from NtMGAM easily with lower energy. Asp542 is an important residue on the bottleneck of the active pocket of NtMGAM and could generate hydrogen bonds with DSK continuously. Our theoretical results may provide some useful clues for designing new α-glucosidase inhibitors to treat type 2 diabetes.

Machine Learning and Enhanced Sampling Simulations for Computing the Potential of Mean Force and Standard Binding Free Energy

Computational capabilities are rapidly increasing, primarily because of the availability of GPU-based architectures. This creates unprecedented simulative possibilities for the systematic and robust computation of thermodynamic observables, including the free energy of a drug binding to a target. In contrast to calculations of relative binding free energy, which are nowadays widely exploited for drug discovery, we here push the boundary of computing the binding free energy and the potential of mean force. We introduce a novel protocol that leverages enhanced sampling, machine learning, and ad hoc algorithms to limit human intervention, computing time, and free parameters in free energy calculations. We first validate the method on a host-guest system, and then we apply the protocol to glycogen synthase kinase 3 beta, a protein kinase of pharmacological interest. Overall, we obtain a good correlation with experimental values in relative and absolute terms. While we focus on protein-ligand binding, the strategy is of broad applicability to any complex event that can be described with a path collective variable. We systematically discuss key details that influence the final result. The parameters and simulation settings are available at PLUMED-NEST to allow full reproducibility.

Free Energy Landscapes from SARS-CoV-2 Spike Glycoprotein Simulations Suggest that RBD Opening Can Be Modulated via Interactions in an Allosteric Pocket

The SARS-CoV-2 coronavirus is an enveloped, positive-sense single-stranded RNA virus that is responsible for the COVID-19 pandemic. The spike is a class I viral fusion glycoprotein that extends from the viral surface and is responsible for viral entry into the host cell and is the primary target of neutralizing antibodies. The receptor binding domain (RBD) of the spike samples multiple conformations in a compromise between evading immune recognition and searching for the host-cell surface receptor. Using atomistic simulations of the glycosylated wild-type spike in the closed and 1-up RBD conformations, we map the free energy landscape for RBD opening and identify interactions in an allosteric pocket that influence RBD dynamics. The results provide an explanation for experimental observation of increased antibody binding for a clinical variant with a substitution in this pocket. Our results also suggest the possibility of allosteric targeting of the RBD equilibrium to favor open states via binding of small molecules to the hinge pocket. In addition to potential value as experimental probes to quantify RBD conformational heterogeneity, small molecules that modulate the RBD equilibrium could help explore the relationship between RBD opening and S1 shedding.

Protein Dynamics and Substrate Protonation States Mediate the Catalytic Action of trans-4-Hydroxy-l-Proline Dehydratase

The enzyme trans-4-hydroxy-l-proline (Hyp) dehydratase (HypD) is among the most abundant glycyl radical enzymes (GREs) in the healthy human gut microbiome and is considered a promising antibiotic target for the prominent antibiotic-resistant pathogen Clostridium difficile. Although an enzymatic mechanism has been proposed, the role of the greater HypD protein environment in mediating radical reactivity is not well understood. To fill this gap in understanding, we investigate HypD across multiple time- and length-scales using electronic structure modeling and classical molecular dynamics. We observe that the Hyp substrate protonation state significantly alters both its enzyme-free reactivity and its dynamics within the enzyme active site. Accurate coupled-cluster modeling suggests the deprotonated form of Hyp to be the most reactive protonation state for C5-H_pro-S activation. In the protein environment, hydrophobic interactions modulate the positioning of the Cys434 radical to enhance the reactivity of C5-H_pro-S abstraction. Long-time dynamics reveal that changing Hyp protonation states triggers the switching of a Leu643-gated water tunnel, a functional feature that has not yet been observed for members of the GRE superfamily.

Electrochemical adsorption of hydrogen on mixed Pd2Pt nanostructures

In the present contribution we have focused on the electrochemical adsorption of a proton from the solution-the Volmer reaction-on a variety of systems based on bimetallic nanostructures-clusters and wires-of Pd and Pt deposited on a surface of Au(111). We have calculated the free energy surface for the electron transfer step by a combination of DFT calculations, MD simulations and the theory of electrocatalysis. We analyze in detail the interaction of the metal d band with the valence orbital of the hydrogen and its effect on the catalytic activity as well as several aspects that influence the electrode reactivity such as spatial arrangements of the nanostructures, the solvation shell and chemical factors. We found that the mixed Pd₂Pt wire interacts strongly with hydrogen, and retains an almost complete solvation shell, which is reflected in a substantially reduced activation energy for the Volmer step. Thus, Pd₂Pt wires on Au(111) are predicted to be efficient electrocatalysts for the reaction.

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.

Adaptative Steered Molecular Dynamics Study of Mutagenesis Effects on Calcium Affinity in the Regulatory Domain of Cardiac Troponin C

Calcium-dependent cardiac muscle contraction is regulated by the protein complex troponin (cTn) and specifically by the regulatory N-terminal domain (N-cTnC) which contains one active Ca²⁺ binding site (site II). It has been previously shown that cardiac muscle contractility and functionality is affected by mutations in N-cTnC which alter calcium binding affinity. Here, we describe the application of adaptive steered molecular dynamics to characterize the influence of N-cTnC mutations on site II calcium binding affinity. We observed the correct trends for all of the studied calcium sensitizing and desensitizing mutants, in conjunction with loop II perturbations. Additionally, the potential of mean force accuracy was shown to increase substantially with increasingly slower speeds and using fewer trajectories. This study presents a novel approach to computationally estimate the Ca²⁺ binding affinity of N-cTnC structures and is a valuable potential tool to support the design and characterization of novel mutations with potential therapeutic benefits.

Molecular conformations and dynamics of nucleotide repeats associated with neurodegenerative diseases: double helices and CAG hairpin loops

Pathogenic DNA secondary structures have been identified as a common and causative factor for expansion in trinucleotide, hexanucleotide, and other simple sequence repeats. These expansions underlie about fifty neurological and neuromuscular disorders known as “anticipation diseases”. Cell toxicity and death have been linked to the pathogenic conformations and functional changes of the RNA transcripts, of DNA itself and, when trinucleotides are present in exons, of the translated proteins. We review some of our results for the conformations and dynamics of pathogenic structures for both RNA and DNA, which include mismatched homoduplexes formed by trinucleotide repeats CAG and GAC; CCG and CGG; CTG(CUG) and GTC(GUC); the dynamics of DNA CAG hairpins; mismatched homoduplexes formed by hexanucleotide repeats (GGGGCC) and (GGCCCC); and G-quadruplexes formed by (GGGGCC) and (GGGCCT). We also discuss the dynamics of strand slippage in DNA hairpins formed by CAG repeats as observed with single-molecule Fluorescence Resonance Energy Transfer. This review focuses on the rich behavior exhibited by the mismatches associated with these simple sequence repeat noncanonical structures.

Bias free multiobjective active learning for materials design and discovery

The design rules for materials are clear for applications with a single objective. For most applications, however, there are often multiple, sometimes competing objectives where there is no single best material and the design rules change to finding the set of Pareto optimal materials. In this work, we leverage an active learning algorithm that directly uses the Pareto dominance relation to compute the set of Pareto optimal materials with desirable accuracy. We apply our algorithm to de novo polymer design with a prohibitively large search space. Using molecular simulations, we compute key descriptors for dispersant applications and drastically reduce the number of materials that need to be evaluated to reconstruct the Pareto front with a desired confidence. This work showcases how simulation and machine learning techniques can be coupled to discover materials within a design space that would be intractable using conventional screening approaches.

Insight Into Seeded Tau Fibril Growth From Molecular Dynamics Simulation of the Alzheimer's Disease Protofibril Core

Aggregates of the microtubule associated tau protein are a major constituent of neurofibrillary lesions that define Alzheimer’s disease (AD) pathology. Increasing experimental evidence suggests that the spread of tau neurofibrillary tangles results from a prion-like seeding mechanism in which small oligomeric tau fibrils template the conversion of native, intrinsically disordered, tau proteins into their pathological form. By using atomistic molecular dynamics (MD) simulations, we investigate the stability and dissociation thermodynamics of high-resolution cryo-electron microscopy (cryo-EM) structures of both the AD paired-helical filament (PHF) and straight filament (SF). Non-equilibrium steered MD (SMD) center-of-mass pulling simulations are used to probe the stability of the protofibril structure and identify intermolecular contacts that must be broken before a single tau peptide can dissociate from the protofibril end. Using a combination of exploratory metadynamics and umbrella sampling, we investigate the complete dissociation pathway and compute a free energy profile for the dissociation of a single tau peptide from the fibril end. Different features of the free energy surface between the PHF and SF protofibril result from a different mechanism of tau unfolding. Comparison of wild-type tau PHF and post-translationally modified pSer356 tau shows that phosphorylation at this site changes the dissociation free energy surface of the terminal peptide. These results demonstrate how different protofibril morphologies template the folding of endogenous tau in distinct ways, and how post-translational modification can perturb the folding mechanism.

Serine 477 plays a crucial role in the interaction of the SARS-CoV-2 spike protein with the human receptor ACE2

Since the worldwide outbreak of the infectious disease COVID-19, several studies have been published to understand the structural mechanism of the novel coronavirus SARS-CoV-2. During the infection process, the SARS-CoV-2 spike (S) protein plays a crucial role in the receptor recognition and cell membrane fusion process by interacting with the human angiotensin-converting enzyme 2 (hACE2) receptor. However, new variants of these spike proteins emerge as the virus passes through the disease reservoir. This poses a major challenge for designing a potent antigen for an effective immune response against the spike protein. Through a normal mode analysis (NMA) we identified the highly flexible region in the receptor binding domain (RBD) of SARS-CoV-2, starting from residue 475 up to residue 485. Structurally, the position S477 shows the highest flexibility among them. At the same time, S477 is hitherto the most frequently exchanged amino acid residue in the RBDs of SARS-CoV-2 mutants. Therefore, using MD simulations, we have investigated the role of S477 and its two frequent mutations (S477G and S477N) at the RBD during the binding to hACE2. We found that the amino acid exchanges S477G and S477N strengthen the binding of the SARS-COV-2 spike with the hACE2 receptor.

Recent advances in computational methods for biosensor design

Biosensors are analytical tools with a great application in healthcare, food quality control, and environmental monitoring. They are of considerable interest to be designed by using cost-effective and efficient approaches. Designing biosensors with improved functionality or application in new target detection has been converted to a fast-growing field of biomedicine and biotechnology branches. Experimental efforts have led to valuable successes in the field of biosensor design; however, some deficiencies restrict their utilization for this purpose. Computational design of biosensors is introduced as a promising key to eliminate the gap. A set of reliable structure prediction of the biosensor segments, their stability, and accurate descriptors of molecular interactions are required to computationally design biosensors. In this review, we provide a comprehensive insight into the progress of computational methods to guide the design and development of biosensors, including molecular dynamics simulation, quantum mechanics calculations, molecular docking, virtual screening, and a combination of them as the hybrid methodologies. By relying on the recent advances in the computational methods, an opportunity emerged for them to be complementary or an alternative to the experimental methods in the field of biosensor design.

High-throughput virtual screening of drug databanks for potential inhibitors of SARS-CoV-2 spike glycoprotein

COVID-19, which is caused by a novel coronavirus known as SARS-CoV-2, has spread rapidly around the world, and it has infected more than 29 million individuals as recorded on 16 September 2020. Much effort has been made to stop the virus from spreading, and there are currently no approved pharmaceutical products to treat COVID-19. Here, we apply an in silico approach to investigate more than 3800 FDA approved drugs on the viral RBD S₁-ACE2 interface as a target. The compounds were investigated through flexible ligand docking, ADME property calculations and protein–ligand interaction maps. Molecular dynamics (MD) simulations were also performed on eleven compounds to study the stability and the interactions of the protein–ligand complexes. The MD simulations show that bagrosin, chidamide, ebastine, indacaterol, regorafenib, salazosulfadimidine, silodosin and tasosartan are relatively stable near the C terminal domain (CTD1) of the S₁ subunit of the viral S protein. The relative MMGBSA binding energies show that silodosin has the best binding to the target. The constant velocity steered molecular dynamics (SMD) simulations show that silodosin preferentially interacts with the RBD S₁ and has potential to act as an interfering compound between viral spike–host ACE2 interactions.

Communicated by Ramaswamy H. Sarma

Nicotinic receptor pharmacology in silico: Insights and challenges

Nicotinic acetylcholine receptors (nAChR) are homo- or hetero-pentameric ligand-gated ion channels of the Cys-loop superfamily and play important roles in the nervous system and muscles. Studies on nAChR benefit from in silico modeling due to the lack of high-resolution structures for most receptor subtypes and challenges in experiments addressing the complex mechanism of activation involving allosteric sites. Although there is myriad of computational modeling studies on nAChR, the multitude of the methods and parameters used in these studies makes modeling nAChR a daunting task, particularly for the non-experts in the field. To address this problem, the modeling literature on Torpedo nAChR and α7 nAChR were focused on as examples of heteromeric and homomeric nAChR, and the key in silico modeling studies between the years 1995-2019 were concisely reviewed. This was followed by a critical analysis of these studies by comparing the findings with each other and with the emerging experimental and computational data on nAChR. Based on these critical analyses, suggestions were made to guide the future researchers in the field of in silico modeling of nAChR. This article is part of the special issue on 'Contemporary Advances in Nicotine Neuropharmacology'.

How Good is Jarzynski's Equality for Computer-Aided Drug Design?

Accurate determination of the binding affinity of the ligand to the receptor remains a difficult problem in computer-aided drug design. Here, we study and compare the efficiency of Jarzynski’s equality (JE) combined with steered molecular dynamics and the linear interaction energy (LIE) method by assessing the binding affinity of 23 small compounds to six receptors, including β-lactamase, thrombin, factor Xa, HIV-1 protease (HIV), myeloid cell leukemia-1, and cyclin-dependent kinase 2 proteins. It was shown that Jarzynski’s nonequilibrium binding free energy ΔG_neq^Jar correlates with the available experimental data with the correlation levels R = 0.89, 0.86, 0.83, 0.80, 0.83, and 0.81 for six data sets, while for the binding free energy ΔG_LIE obtained by the LIE method, we have R = 0.73, 0.80, 0.42, 0.23, 0.85, and 0.01. Therefore, JE is recommended to be used for ranking binding affinities as it provides accurate and robust results. In contrast, LIE is not as reliable as JE, and it should be used with caution, especially when it comes to new systems.

Slippery and Wear-Resistant Surfaces Enabled by Interface Engineered Graphene

Friction and wear remain the primary cause of mechanical energy dissipation and system failure. Recent studies reveal graphene as a powerful solid lubricant to combat friction and wear. Most of these studies have focused on nanoscale tribology and have been limited to a few specific surfaces. Here, we uncover many unknown aspects of graphene's contact-sliding at micro- and macroscopic tribo-scales over a broader range of surfaces. We discover that graphene's performance reduces for surfaces with increasing roughness. To overcome this, we introduce a new type of graphene/silicon nitride (SiN_x, 3 nm) bilayer overcoats that exhibit superior performance compared to native graphene sheets (mono and bilayer), that is, display the lowest microscale friction and wear on a range of tribologically poor flat surfaces. More importantly, two-layer graphene/SiN_x bilayer lubricant (<4 nm in total thickness) shows the highest macroscale wear durability on tape-head (topologically variant surface) that exceeds most previous thicker (∼7-100 nm) overcoats. Detailed nanoscale characterization and atomistic simulations explain the origin of the reduced friction and wear arising from these nanoscale coatings. Overall, this study demonstrates that engineered graphene-based coatings can outperform conventional coatings in a number of technologies.

In Silico Modeling of the α7 Nicotinic Acetylcholine Receptor: New Pharmacological Challenges Associated with Multiple Modes of Signaling

The α7 nicotinic acetylcholine receptor is a homopentameric ion-channel of the Cys-loop superfamily characterized by its low probability of opening, high calcium permeability, and rapid desensitization. The α7 receptor has been targeted for the treatment of the cognitive symptoms of schizophrenia, depression, and Alzheimer's disease, but it is also involved in inflammatory modulation as a part of the cholinergic anti-inflammatory pathway. Despite its functional importance, in silico studies of the α7 receptor cannot produce a general model explaining the structural features of receptor activation, nor predict the mode of action for various ligand classes. Two particular problems in modeling the α7 nAChR are the absence of a high-resolution structure and the presence of five potentially nonequivalent orthosteric ligand binding sites. There is wide variability regarding the templates used for homology modeling, types of ligands investigated, simulation methods, and simulation times. However, a systematic survey focusing on the methodological similarities and differences in modeling α7 has not been done. In this work, we make a critical analysis of the modeling literature of α7 nAChR by comparing the findings of computational studies with each other and with experimental studies under the main topics of structural studies, ligand binding studies, and comparisons with other nAChR. In light of our findings, we also summarize current problems in the field and make suggestions for future studies concerning modeling of the α7 receptor.

A comparative study for the intermediate states of myelin oligodendrocyte glycoprotein in the absence and presence of glycan - A computational approach

Myelin Oligodendrocyte glycoprotein (MOG) is found to play an important role in providing structural integrity to myelin sheath at the same time it acts as an auto-antigen which might lead to Multiple Sclerosis (MS). What causes this specific property of being an auto-antigen is still not known. Here we present molecular dynamics simulation studies of unfolding and folding of the protein MOG in both the absence and presence of N-glycan in order to understand the role of glycosylation in the stability and flexibility of the protein. The main results from these studies show that the glycosylation increases the stability of the protein MOG and inhibits the complete unfolding of MOG in the SMD. From the folding studies using TMD, it was observed that the glycan helps the protein to attain the near-native folded conformation. However, it was also observed from the direct TMD studies that the pathway of protein folding was enhanced by the trace-back of intermediate states in the presence of glycan.

Catalytic Mechanism and Covalent Inhibition of UDP-N-Acetylglucosamine Enolpyruvyl Transferase (MurA): Implications to the Design of Novel Antibacterials

UDP-N-acetylglucosamine enolpyruvyl transferase (MurA) catalyzes the first step in the biosynthesis of the bacterial cell wall. This pathway is essential for the growth of bacteria but missing in mammals, that nominates MurA as an attractive antibacterial target. MurA has a flexible loop whose conformational change is known to be part of the activation mechanism of the enzyme. We have shown that the loop closed conformation makes the proton transfer from Cys115 to His394 possible by a low barrier exothermic process. QM/MM MD simulations revealed that the activated thiolate is able to react with phosphoenolpyruvate (PEP), the natural substrate of MurA. The binding free energy profile of several covalent inhibitors with various warheads reacting with the activated Cys115 was calculated by QM/MM MD simulations and confirmed that reaction barrier heights tend to separate active from inactive compounds. Our results give new insight into the catalytic mechanism and covalent inhibition of MurA and suggest that QM/MM MD simulations are able to support ligand design by providing sensible relative free energy barriers for covalent inhibitors with various warheads reacting with thiolate nucleophiles.

Interplay of Nanoparticle Rigidity and Its Translocation Ability through Cell Membrane

Understanding the endocytic process of nanoparticles (NPs) with different mechanical rigidities is critical to develop effective drug delivery vectors. Here, we perform experiments, coarse-grained molecular dynamics simulations, and theoretical analyses to investigate the role of NPs' mechanical rigidity in the cellular endocytic process. Experiments based on two types of engineered Au NPs that have similar properties but different rigidities are performed in order to investigate their cellular uptake efficiencies, and it has been found that the more rigid NPs can achieve a higher cellular uptake efficiency. Simulation results confirm that rigid NPs can achieve full internalization by forming a complete double-layer endosome coating, while relatively soft NPs can only reach 40% surface coverage by membrane lipids. Simulation results capture an intriguing translocation of multiple NPs with different rigidities in a cooperative manner where the NPs' mechanical rigidities regulate their translocation efficiencies. We find that theoretically rigid NPs require less energy to overcome the energy barrier for membrane internalization than soft NPs do, which is in good agreement with experiment and simulation results. This synergetic study offers useful insight into the design principle of a general NP-based drug delivery vector as well as the promising biomedical application of NP-based medicine.

ART-RRT: As-Rigid-As-Possible search for protein conformational transition paths

The possible functions of a protein are strongly related to its structural rearrangements in the presence of other molecules or environmental changes. Hence, the evaluation of transition paths of proteins, which encodes conformational changes between stable states, is important since it may reveal the underlying mechanisms of the biochemical processes related to these motions. During the last few decades, different geometry-based methods have been proposed to predict such transition paths. However, in the cases where the solution requires complex motions, these methods, which typically constrain only locally the molecular structures, could produce physically irrelevant solutions involving self-intersection. Recently, we have proposed ART-RRT, an efficient method for finding ligand-unbinding pathways. It relies on the exploration of energy valleys in low-dimensional spaces, taking advantage of some mechanisms inspired from computer graphics to ensure the consistency of molecular structures. This article extends ART-RRT to the problem of finding probable conformational transition between two stable states for proteins. It relies on a bidirectional exploration rooted on the two end states and introduces an original strategy to attempt connections between the explored regions. The resulting method is able to produce at low computational cost biologically realistic paths free from self-intersection. These paths can serve as valuable input to other advanced methods for the study of proteins. A better understanding of conformational changes of proteins is important since it may reveal the underlying mechanisms of the biochemical processes related to such motions. Recently, the ART-RRT method has been introduced for finding ligand-unbinding pathways. This article presents an adaptation of the method for finding probable conformational transition between two stable states of a protein. The method is not only computationally cost-effective but also able to produce biologically realistic paths which are free from self-intersection.

A Quantitative Metric for the Design of Selective Supercritical CO2 Extraction of Lithium from Geothermal Brine

As demand grows for lithium, its recovery from geothermal brines provides an attractive alternative to slow mining. One promising extraction method uses crown ethers as extractants in supercritical carbon dioxide with cation exchangers to facilitate extraction from brine. Molecular dynamics modeling is used to understand the mechanism of binding between lithium (or sodium) and combinations of 14-crown-4 ethers and cation exchangers, and the predictive capability of computational modeling to test lithium selectivity is established for four combinations of crown ethers [methylene-14-crown-4 (M14C4) and a fluorinated 14-crown-4 (F14C4)] and cation exchangers [di(2-ethyl-hexyl)phosphoric acid (HDEHP) and tetraethylammonium perfluoro-1-octanesulfonate (TPFOS)]. Binding free energies (given in kcal mol^-1 ) of lithium and sodium, respectively, to crown ether-cation exchangers are 85 and 71 for M14C4-HDEHP, 90 and 71 for F14C4-HDEHP, 93 and 80 for M14C4-TPFOS, and 104 and 93 for F14C4-TPFOS. Good agreement is found between computational predictions and supercritical carbon dioxide extraction experiments at 60 °C and 250 bar. Binding free energy gives a suitable metric to describe extraction efficiency. Differences in the binding free energies of sodium and lithium to crown ethers determine the extraction selectivity. Fluorine groups are found to exert a positive influence to optimize extraction efficiency. Of the systems studied, F14C4 with TPFOS offers the most selective and efficient extraction system.

Surface Topography Effects of Globular Biomolecules on Hydration Water

Hydration water serves as a microscopic manifestation of structural stability and functions of biomolecules. To develop bio-nanomaterials in applications, it is important to study how the surface topography and heterogeneity of biomolecules result in their diversity of the hydration dynamics and energetics. We here performed molecular dynamics simulations combined with the steered molecular dynamics and umbrella sampling to investigate the dynamics and escape process associated with the free energy change of water molecules close to a globular biomolecule, i.e., hemoglobin (Hb) and G-quadruplex DNA (GDNA). The residence time, power of long-time tail, and dipole relaxation time were found to display drastic changes within the averaged hydration shell of 3.0-5.0 Å. Compared with bulk water, in the inner hydration shell, the water dipole moment displays a slower relaxation process and is more oriented toward GDNA than toward Hb, forming a hedgehog-like structure when it surrounds GDNA. In particular, a spine water structure is observed in the GDNA narrow groove. The water isotope effect not only prolongs the dynamic time scales of libration motion in the inner hydration shell and the dipole relaxation processes in the bulk but also strengthens the DNA spine water structure. The potential of the mean force profile reflects the integrity of the hydration shell structure and enables us to obtain detailed insights into the structures formed by water, such as the caged H-bond network and the edge bridge structures; it also reveals that local hydration shell free energy (LHSFE) depends on H-bond rupture processes and ranges from 0.2 to 4.2 kcal/mol. Our results demonstrate that the surface topography of a biomolecule influences the integrity of the hydration shell structure and LHSFE. Our studies are able to identify various further applications in the areas of microfluid devices and nano-dewetting on bioinspired surfaces.

Quantum Mechanical Description of Electrostatics Provides a Unified Picture of Catalytic Action Across Methyltransferases

Methyl transferases (MTases) are a well-studied class of enzymes for which competing enzymatic enhancement mechanisms have been suggested, ranging from structural methyl group CH···X hydrogen bonds (HBs) to electrostatic- and charge-transfer-driven stabilization of the transition state (TS). We identified all Class I MTases for which reasonable resolution (<2.0 Å) crystal structures could be used to form catalytically competent ternary complexes for multiscale (i.e., quantum-mechanical/molecular-mechanical or QM/MM) simulation of the S_N2 methyl transfer reaction coordinate. The four Class I MTases studied have both distinct functions (e.g., protein repair or biosynthesis) and substrate nucleophiles (i.e., C, N, or O). While CH···X HBs stabilize all reactant complexes, no universal TS stabilization role is found for these interactions in MTases. A consistent picture is instead obtained through analysis of charge transfer and electrostatics, wherein much of cofactor-substrate charge separation is maintained in the TS region, and electrostatic potential is correlated with substrate nucleophilicity (i.e., intrinsic reactivity).

Molecular Modeling of Complex Cross-Linked Networks of PEGDA Nanogels

Poly(ethylene glycol) (PEG)-based nanogels are attractive for biomedical applications due to their biocompatibility, versatile end group chemistry, and ability to sterically shield encapsulated drug molecules. The characteristics of a hydrogel network govern the encapsulation and efficient delivery of drug molecules for a target application. A molecular-level description of network topology can complement experimental investigations to understand its effects on the structural properties of these nanogels. In this work, atomistic molecular simulations of heterogeneous, nonideal PEG-diacrylate (PEGDA) nanogels are presented. The effects of cross-linking density and topological features on the structural properties of PEGDA nanogels were studied. The average functionality was controlled to systematically study the effect of cross-linking density on the radius of gyration, shape, and mesh size of the nanogels. For a given average functionality, the impact of distinct network topologies on the structural properties was also studied. The aspect ratios, based on the gyration tensor, were calculated to characterize the shapes of these nanogels for different topologies. Nanogel structures with higher cross-linking densities showed a globular shape, while structures with lower cross-linking density showed shape anisotropy. The distribution and connectivity of the cross-linked junctions played a key role in determining the size and shape anisotropy of PEGDA nanogels; the number of unreacted chain ends and their connectivity directly affected the anisotropy. The mesh size, denoted by the limiting "free volume element" present in the nanogel samples, does not show a significant change with increasing average functionality. This work provides insight into the structural properties of heterogeneous hydrogels that aid the design of nonideal nanogel networks for a targeted drug delivery application.

Disclosing the Template-Induced Misfolding Mechanism of Tau Protein by Studying the Dissociation of the Boundary Chain from the Formed Tau Fibril Based on a Steered Molecular Dynamics Simulation

The level of tau aggregation into neurofibrillary tangles, including paired helical filament (PHF) and straight filament (SF), is closely associated with Alzheimer's disease. Despite the pathological importance of misfolding and aggregation of tau, the corresponding mechanism remains unclear. Therefore, to uncover the misfolding mechanism of the tau monomer upon induction of formed PHF and SF, in this study, a conventional molecular dynamics simulation combined with a steered molecular dynamics simulation was performed to study the dissociation of the boundary chain. Interestingly, our results show that the dissociation mechanisms of the boundary chain in PHF and SF are different. In PHF, the boundary chain begins to dissociate from regions β2 and β3 and ends at β8. However, in SF, it is simultaneously dissociated from β1 and β8 and ends at β5. The dissociation of the boundary chain is the reverse of template-induced misfolding of the monomer. Therefore, we can deduce the misfolding mechanism of the monomer upon induction of the template. For PHF, β8 first interacts with the template by hydrophobic interaction. Then β7, β6, β5, β4, and β1 sequentially bind to the template by electrostatic and hydrophobic interactions. After β1 binds to the template, β2 and β3 very quickly bind to the template through hydrophobic interaction. For SF, β5 of the monomer first interacts with the template by electrostatic attraction. Then β4 and β6, β3 and β7, and β2 and β8 bind to the template in turn. Finally, β1 and β8 are fully bound to the template by hydrophobic interaction. The obtained results will be vital for understanding the earlier events during misfolding and aggregation of tau.

Simulated (Un)Binding of Arachidonic Acid in the Cyclooxygenase site of Prostaglandin H2 Synthase-1

Molecular dynamics simulations with external forces are employed to study the unbinding and binding of arachidonic acid (AA) in the cyclooxygenase (COX) site of prostaglandin H2 synthase-1. Simulations with AA inside the COX binding channel reveal sequences of concerted bond rotations in the fatty acid alkyl chain which obviate the need for gross conformational changes in the protein and substrate during unbinding and binding. The all-cis structure of AA, with double bonds separated by two single bonds, facilitates easy access to the COX channel and correct positioning inside the active site for the COX chemistry to occur. Two derivatives of AA, one with a cis double bond changed to a trans configuration and the other with a double bond reduced to a single bond, are also studied. In both cases the concertedness of bond rotations in the fatty acid chain is diminished and larger forces are required to move the fatty acid inside the COX channel. Important motions of residues near the mouth of the COX channel are found and analyzed. In particular, a conformational “switch” involving Arg83, Glu524 and Arg120 is seen to mediate the movement of the substrate from the membrane to the channel.

Force Sensing in Nanoscale: Integration of Virtual Spring for Sensing the Interactions of β-Amyloid Grabbed by Prefoldin

Analysis of atomic forces in molecular scale including the reciprocal internal force directions, quantity, and Young's modulus, and more detailed analysis via diagrams, is one of the most current research topics. This research aims to analyze the bilateral forces applied between a nano-actuator and a cargo from a new perspective. While virtual springs have been used as a tool to convey forces in the previous steered molecular dynamics simulations, this study uses the spring as a sensor for measuring the internal forces. This study introduces atomic sensing via the steered molecular dynamics method. Following the previous studies, the protein employed in this study is the mutated Archetype Prefoldin being used to control the pathogenic cargo beta-amyloid (Alzheimer's). The powerful Gromacs software carries out the simulations for the calculation of the total force and force for each branch. The simulation results illustrate the total force between the cargo and nano-actuator is ~ 2.8 nN, while each branch needs a force of ~ 1.2-1.5 nN to release the cargo. The results demonstrate the validity of the method and applicability of the virtual sensor for assessing the microscopic forces. This investigation is a pioneer study for the advent of the sensor as an assessment tool for the mechanical analysis and precise atomic force studies.

Dissociation Process of a MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and the Markov State Model

Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to a complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within an average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained in combination with one-dimensional-, three-dimensional (3D)- or C_α-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. C_α-MSM based on all C_α coordinates of TAD-p53 reproduced the experimentally measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.

Extended Efficacy of Low-Dose Valganciclovir for Prevention of Cytomegalovirus Disease in Intermediate-Risk Kidney Transplant Recipients: Two-Year Follow-Up

OBJECTIVES

In a previous study, we evaluated 1-year outcomes of using low-dose valganciclovir prophylaxis for cytomegalovirus infection in intermediate-risk kidney transplant recipients. Whether this effect persists in the long term is unknown. We aimed to evaluate the 2-year follow up of such adopted prophylaxis.

MATERIALS AND METHODS

We randomized 2 matched groups of kidney transplant recipients (1:1) to receive valganciclovir as 450 mg daily (group 1) or 900 mg daily (group 2) for the first 6 months after kidney transplant. The final analysis included 196 patients as intermediate-risk patients (98 in each treatment group) after exclusion of 5 high-risk patients. Serologically, all patients were at moderate risk for cytomegalovirus infection. Long-term outcomes including cytomegalovirus disease, acute rejection, new-onset diabetes after transplant, graft loss, and patient survival were assessed.

RESULTS

Through year 2 of follow-up, cytomegalovirus infection was reported in only 1 patient in group 1 (at month 13) and 1 patient in group 2 (at month 19) (not significant). Biopsy-proven acute rejection episodes were not statistically different between the groups (2 episodes in group 1 and 6 in group 2; P = .431). New-onset diabetes posttransplant was reported in 8.1% in group 1 and 13.2% in group 2 (P = .535). Graft failure was equal in both groups (1 in each group) at 2 years of follow up (not significant). Patient survival was comparable in both groups (100% in group 1 versus 97.9% in group 2; P = .661). The total number of cytomegalovirus infections at 2 years was numerically less in group 1 (P = .128).

CONCLUSIONS

Low-dose valganciclovir prophylaxis for 6 months was associated with sustained reduction of cytomegalovirus infection up to 2 years after kidney transplant without significant impact on the acute rejection, new-onset diabetes posttransplant, or patient and graft outcomes.