Search

We introduce a novel enhanced sampling approach named OPES flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the On-the-fly-Probability-Enhanced-Sampling (OPES) approach [Invernizzi and Parrinello, JPC Lett. 2020], which has been recently developed for calculating converged free energy surfaces for complex systems. In this paper, we describe the theoretical details of the OPES flooding technique and demonstrate the application on three systems of increasing complexity: barrier crossing in a two-dimensional double well potential, conformational transition in the alanine dipeptide in gas phase, and the folding and unfolding of the chignolin polypeptide in aqueous environment. From extensive tests, we show that the calculation of accurate kinetics not only requires the transition state to be bias-free, but the amount of bias deposited should also not exceed the effective barrier height measured along the chosen collective variables. In this vein, the possibility of computing rates from biasing suboptimal order parameters has also been explored. Furthermore, we describe the choice of optimum parameter combinations for obtaining accurate results from limited computational effort.

Monoclonal antibodies are emerging as a viable treatment for the coronavirus disease 19 (COVID-19). However, newly evolved variants of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) can reduce the efficacy of currently available antibodies and can diminish vaccine-induced immunity. Here, we demonstrate that the microscopic dynamics of neutralizing monoclonal antibodies can be profoundly modified by the mutations present in the spike proteins of the SARS-COV-2 variants currently circulating in the world population. The dynamical perturbations within the antibody structure, which alter the thermodynamics of antigen recognition, are diverse and can depend both on the nature of the antibody and on the spatial location of the spike mutation. The correlation between the motion of the antibody and that of the spike receptor binding domain (RBD) can also be changed, modulating binding affinity. Using protein-graph-connectivity networks, we delineated the mutant-induced modifications in the information-flow along allosteric pathway throughout the antibody. Changes in the collective dynamics were spatially distributed both locally and across long-range distances within the antibody. On the receptor side, we identified an anchor-like structural element that prevents the detachment of the antibodies; individual mutations there can significantly affect the antibody binding propensity. Our study provides insight into how virus neutralization by monoclonal antibodies can be impacted by local mutations in the epitope via a change in dynamics. This realization adds a new layer of sophistication to the efforts for rational design of monoclonal antibodies against new variants of SARS-CoV2, taking the allostery in the antibody into consideration.

In concurrence with the pursuit of clean and green medium, recent years have witnessed an unprecedented rise in the usage of ionic liquids (ILs).

The study of the rare transitions that take place between long lived metastable states is a major challenge in molecular dynamics simulations. Many of the methods suggested to address this problem rely on the identification of the slow modes of the system, which are referred to as collective variables. Recently, machine learning methods have been used to learn the collective variables as functions of a large number of physical descriptors. Among many such methods, Deep Targeted Discriminant Analysis has proven to be useful. This collective variable is built from data harvested from short unbiased simulations in the metastable basins. Here, we enrich the set of data on which the Deep Targeted Discriminant Analysis collective variable is built by adding data from the transition path ensemble. These are collected from a number of reactive trajectories obtained using the On-the-fly Probability Enhanced Sampling flooding method. The collective variables thus trained lead to more accurate sampling and faster convergence. The performance of these new collective variables is tested on a number of representative examples.

We present six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for enabling applications of the weighted ensemble (WE) path sampling strategy to even larger systems and/or slower processes. The tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of “binless” schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling at the system biology level. Applications of the tutorials range from atomistic to residue-level to non-spatial models, and include complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics simulations and completed the previous suite of WESTPA tutorials.

Allostery in proteins involves, broadly speaking, ligand-induced conformational transitions that modulate function at active sites distal to where the ligand binds. In contrast, the concept of cooperativity (in the sense used in phase transition theory) is often invoked to understand protein folding and, therefore, function. The modern view on allostery is one based on dynamics and hinges on the time-dependent interactions between key residues in a complex network, interactions that determine the free-energy profile for the reaction at the distal site. Here, we merge allostery and cooperativity, and we discuss a joint model with features of both. In our model, the active-site reaction is replaced by the reaction pathway that leads to protein folding, and the presence or absence of the effector is replaced by mutant-vs-wild type changes in key residues. To this end, we employ our recently introduced time-lagged independent component analysis (tICA) correlation approach [Ray et al. Proc. Natl. Acad. Sci. 118(43) (2021), e2100943118] to identify the allosteric role of distant residues in the folded-state dynamics of a large protein. In this work, we apply the technique to identify key residues that have a significant role in the folding of a small, fast folding-protein, chignolin. Using extensive enhanced sampling simulations, we critically evaluate the accuracy of the predictions by mutating each residue one at a time and studying how the mutations change the underlying free energy landscape of the folding process. We observe that mutations in those residues whose associated backbone torsion angles have a high correlation score can indeed lead to loss of stability of the folded configuration. We also provide a rationale based on interaction energies between individual residues with the rest of the protein to explain this effect. From these observations, we conclude that the tICA correlation score metric is a useful tool for predicting the role of individual residues in the correlated dynamics of proteins and can find application to the problem of identifying regions of protein that are either most vulnerable to mutations or— mutatis mutandis—to binding events that affect their functionality.

The molecular features that dictate interactions between functionalized nanoparticles and biomolecules are not well understood. This is in part because for highly charged nanoparticles in solution, establishing a clear connection between the molecular features of surface ligands and common experimental observables such as ζ potential requires going beyond the classical models based on continuum and mean field models. Motivated by these considerations, molecular dynamics simulations are used to probe the electrostatic properties of functionalized gold nanoparticles and their interaction with a charged peptide in salt solutions. Counterions are observed to screen the bare ligand charge to a significant degree even at the moderate salt concentration of 50 mM. As a result, the apparent charge density and ζ potential are largely insensitive to the bare ligand charge densities, which fall in the range of ligand densities typically measured experimentally for gold nanoparticles. While this screening effect was predicted by classical models such as the Manning condensation theory, the magnitudes of the apparent surface charge from microscopic simulations and mean-field models are significantly different. Moreover, our simulations found that the chemical features of the surface ligand (e.g., primary vs. quaternary amines, heterogeneous ligand lengths) modulate the interfacial ion and water distributions and therefore the interfacial potential. The importance of interfacial water is further highlighted by the observation that introducing a fraction of hydrophobic ligands enhances the strength of electrostatic binding of the charged peptide. Finally, the simulations highlight that the electric double layer is perturbed upon binding interactions. As a result, it is the bare charge density rather than the apparent charge density or ζ potential that better correlates with binding affinity of the nanoparticle to a charged peptide. Overall, our study highlights the importance of molecular features of the nanoparticle/water interface and underscores a set of design rules for the modulation of electrostatic driven interactions at nano/bio interfaces.

The aggregation phenomena of two ionic liquid surfactants—1-decyl-3-methylimidazolium chloride (C10MeImCl) and 1-hexadecyl-3-methylimidazolium chloride (C16MeImCl) have been investigated in aqueous solutions in the presence of a polyelectrolyte, sodium polystyrenesulfonate (NaPSS), using electrical conductivity, surface tension, fluorescence and vapor pressure measurements. The counterion-condensation behavior of aqueous NaPSS has also been studied and approximately 62% of the counterions have been found to remain uncondensed in these solutions. The characteristic concentrations signifying different kinds of interactions, namely the critical aggregation concentration and the polymer saturation concentration, in the premicellar regime have been identified. Surface-active complexes with bound surfactant monomers onto the polyion backbone, non-surface active aggregates with small micellar units wrapped by polyion chains, and non-surface active micelles of ionic liquid surfactant monomers, were found to form in aqueous ionic liquid surfactant solutions in the presence of NaPSS. The surfactant concentration where polyelectrolyte–surfactant monomer complexes begin to form (cac), the concentration where the polyion chains get saturated with the small micellar units (psc), and the concentration where micelles of ionic liquid surfactant molecules begin to appear (cmc*) were identified in aqueous C16MeImCl solutions in the presence of NaPSS. Only two critical concentrations, namely the psc, and cmc* were, however, detected in aqueous C10MeImCl solutions in the presence of NaPSS. Effects of alkyl chain length of the ionic liquid surfactants, temperature, NaPSS concentration, and the charge parameter of NaPSS on the self-aggregation of the ionic liquid surfactants have been considered to elucidate the interactions in these mixed systems. The thermodynamics of micellization in these systems have also been studied and the spontaneity of the polyelectrolyte-induced micellization processes have been rationalized for the systems investigated, with the entropy terms superseding the enthalpy terms.

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-timescale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond timescale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy, but it also allows for a scheme to reconstruct a multidimensional free energy landscape along additional degrees of freedom which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to four orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential application in kinetics and free-energy based computational drug design.

We consider the recently developed weighted ensemble milestoning (WEM) scheme [J. Chem. Phys. 152, 234114 (2020)], and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na$^+$/Cl$^-$ ion pair and 4-hydroxy-2-butanone (BUT) ligand with FK506 binding protein (FKBP). As proof or principle, we show that the WEM formalism reproduces the Na$^+$/Cl$^-$ ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constant and the committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant ($k_{\text{on}}$) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Genetic information is encoded in the DNA double helix which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleobase pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen base pairing pattern in naked DNA duplexes and estimated its relative stability and lifetime. In contrast, HG transitions in RNA were not observed. Understanding Hoogsteen (HG) base pairing is important because the underlying "breathing" can modulate significantly DNA/RNA recognition by proteins. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force bias method) and Markov State modeling to obtain accurate free energy, kinetics and the intermediates in the transition pathway between WC and HG base pair for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extra-helical intermediates in this seemingly simple process of base pair conformation switching in B-DNA. Extending our calculation to A-RNA, for which HG base pair is not observed experimentally, resulted in relatively unstable single hydrogen bonded distorted Hoogsteen like base pair. Unlike B-DNA the transition pathway primarily involved base paired and intra-helical intermediates with transition timescales much higher than that of B-DNA. The seemingly obvious flip-over reaction coordinate, i.e., the glycosidic torsion angle is unable to resolve the intermediates; so a multidimensional picture, involving backbone dihedral angles and distance between atoms participating in hydrogen bonds, is required to gain insight into the molecular mechanism.SIGNIFICANCEFormation of unconventional Hoogsteen (HG) base pairing is an important problem in DNA biophysics owing to its key role in facilitating the binding of DNA repairing enzymes, proteins and drugs to damaged DNA. X-ray crystallography and NMR relaxation experiments revealed the presence of HG base pair in naked DNA duplex and protein-DNA complex but no HG base pair was observed in RNA. Molecular dynamics simulations could reproduce the experimental free energy cost of HG base pairing in DNA although a detailed mechanistic insight is still missing. We performed enhanced sampling simulation and Markov state modeling to obtain accurate free energy, kinetics and the intermediates in the transition pathway between WC and HG base pair for both B-DNA and A-RNA.

The aggregation behaviour of two imidazolium-based surface active ionic liquids (SAILs), e.g. , 1-decyl-3-methylimidazolium chloride (C 10 MeImCl) and 1-hexadecyl- 3-methylimidazolium chloride (C 16 MeImCl) in aqueous solutions was investigated in presence of varying amount of poly(ethylene oxide) (PEO) to explore the possible SAIL-PEO interactions. Electrical conductivity measurements were performed in aqueous solutions of PEO having two different molar mass over a selected temperature range. The critical micellization concentrations of the investigated SAILs in PEO-water media were always found to be higher than those in pure water. PEO was found to play a dual role in the SAIL-PEO-water ternaries. It has been concluded that the PEO molecules added to water modify the property of the medium besides binding a very small fraction of the surfactant ions prior to the commencement of the micellization of SAILs. The thermodynamic parameters of micellization such as Gibbs free energy ( Δ G m 0 ), standard enthalpy ( Δ H m 0 ), and standard entropy ( Δ S m 0 ) of aggregation were estimated from the temperature dependence of the critical micellization concentration values. The effects of concentrations of PEO, and molar mass of PEO on the self-aggregation of the C 10 MeImCl and C 16 MeImCl have been discussed in terms of the phenomena prevailing in these solutions. The observed enthalpy-entropy compensation effect in these SAIL-PEO-water ternaries provided important information as to the behaviour of the PEO-water media with respect to the ions present in solution.

To directly simulate rare events using atomistic molecular dynamics is a significant challenge in computational biophysics. Well-established enhanced-sampling techniques do exist to obtain the thermodynamic functions for such systems. However, developing methods for obtaining the kinetics of long timescale processes from simulation at atomic detail is comparatively less developed an area. Milestoning and the weighted ensemble (WE) method are two different stratification strategies; both have shown promise for computing long timescales of complex biomolecular processes. Nevertheless, both require a significant investment of computational resources. We have combined WE and milestoning to calculate observables in orders-of-magnitude less central processing unit and wall-clock time. Our weighted ensemble milestoning method (WEM) uses WE simulation to converge the transition probability and first passage times between milestones, followed by the utilization of the theoretical framework of milestoning to extract thermodynamic and kinetic properties of the entire process. We tested our method for a simple one-dimensional double-well potential, for an eleven-dimensional potential energy surface with energy barrier, and on the biomolecular model system alanine dipeptide. We were able to recover the free energy profiles, time correlation functions, and mean first passage times for barrier crossing events at a significantly small computational cost. WEM promises to extend the applicability of molecular dynamics simulation to slow dynamics of large systems that are well beyond the scope of present day brute-force computations.

Genetic information is encoded in the DNA double helix, which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleo-base pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen (HG) base pairing pattern in naked DNA duplexes, and estimated its relative stability and lifetime. In contrast with DNA, such structures were not observed in RNA duplexes. Understanding HG base pairing is important because the underlying "breathing" motion between the two conformations can significantly modulate protein binding. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force-bias method) and Markov state modeling to obtain accurate free energy, kinetics, and the intermediates in the transition pathway between Watson-Crick and HG base pairs for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extrahelical intermediates in the seemingly simple process of base flipping in B-DNA. Extending our calculation to A-RNA, for which HG base pairing is not observed experimentally, resulted in relatively unstable, single-hydrogen-bonded, distorted Hoogsteen-like bases. Unlike B-DNA, the transition pathway primarily involved base paired and intrahelical intermediates with transition timescales much longer than that of B-DNA. The seemingly obvious flip-over reaction coordinate (i.e., the glycosidic torsion angle) is unable to resolve the intermediates. Instead, a multidimensional picture involving backbone dihedral angles and distance between hydrogen bond donor and acceptor atoms is required to gain insight into the molecular mechanism.

A theoretical study on the ionization dynamics of carbon atom irradiated with a few-cycle, intense laser field is performed within a quasiclassical model to get mechanistic insights into an earlier reported carrier-envelope phase dependency of ionization probabilities of an atom [ Phys. Rev. Lett. 2013, 110, 083602]. The carrier-envelope phase of the laser pulse is found to govern the overall dynamics, reflecting its importance in controlling electronic motion. To understand the origin of this effect, individual trajectories were analyzed at a particular laser intensity. We found that a variation in the carrier-envelope phase affects the angle of ejection of the electrons and subsequently the attainment of the desired final state.

The influence of alloying on mode-selectivity in H

Electrical conductivities of aqueous solutions of an anionic polyelectrolyte sodium carboxymethylcellulose with different degrees of substitution have been reported as a function of concentration and temperature. The conductance-concentration data have been analyzed with the modified Manning counterion-condensation model appropriate for semidilute polyelectrolyte solutions, considering the de Gennes scaling description of the polyion configuration. 36–55% of the counterions were found to remain condensed and the condensation process was shown to be spontaneous. Although the level of counterion-condensation depends appreciably on the degree of substitution, it was found to remain invariant with the polyelectrolyte concentration and the temperature. The free energies of counterion-condensation, the ratios of the mobility of the polyion to that of the counterion, the polyion transference numbers, and the coefficients of monomer-solvent friction have been evaluated. The effects of the polyion charge, its size, and the bulk viscosity of the media on the polyion-counterion interactions have also been discussed.

Oxidation of the three isomeric pentanols by potassium dichromate in aqueous acidic medium was investigated spectrophotometrically at 25 °C. The oxidation leads to the formation of corresponding carbonyl products as identified by 1 H NMR measurement and other techniques. The effects of cationic (cetylpyridinium chloride, CPCl), anionic (sodium dodecylsulphate, NaDS) and nonionic surfactant (polyoxyethylene octyl phenyl ether, TX-100) on reaction kinetics were studied. Moreover the catalytic effect of surfactant in presence of three heteroaromatic nitrogen bases: 2-picolinic acid (2-PA), 2,2′-bipyridyl (2,2′-bipy) and 1,10-phenanthroline (1,10-phen) on the oxidation were also investigated. The value of the rate constant in promoted and unpromoted reaction has been found to be higher in NaDS than in TX-100 and CPCl micellar medium. Conductivity and tensiometry measurements were carried out to determine critical micellization concentration and free energy of micellization (Δ G m °) of the used surfactant in presence of three pentanols. Among the different surfactant and promoter the combination of 2,2′-bipy and NaDS is efficient to exhibit ∼1000 fold rate enhancement for oxidation of 1-pentanol and it is ‘selective’ over the rest of two pentanols oxidation.

Copper-Nickel bimetallic alloys are emerging heterogeneous catalysts for water dissociation which is the rate-determining step of the industrially important water gas shift (WGS) reaction. Yet, the detailed quantum dynamics studies of water-surface scattering in the literature are limited to pure metal surfaces. We present here a three-dimensional wave packet dynamics study of water dissociation on Cu-Ni alloy surfaces, using a pseudodiatomic model of water on a London-Eyring-Polanyi-Sato (LEPS) potential energy surface in order to study the effect of initial vibration, rotation, and orientation of the water molecule on the reactivity. For all the chosen surfaces, reactivity increases significantly with vibrational excitation. In general, for lower vibrational states the reactivity increases with increasing rotational excitation but it decreases in higher vibrational states. Molecular orientation strongly affects reactivity by helping the molecule to align along the reaction path at higher vibrational states. For different alloys, the reaction probability follows the trend of barrier heights and the surfaces having all Ni atoms in the uppermost layer are much more reactive than the ones with Cu atoms. Hence the nature of the alloy surface and initial quantum state of the incoming molecule significantly influence the reactivity in surface catalyzed water dissociation.

Aggregation behavior of a surface active ionic liquid 1-decyl-3-methylimidazolium chloride (C 10 MeImCl) was studied in aqueous solutions in absence and in presence of sodium carboxymethylcellulose (NaCMC) by electrical conductivity, surface tension, vapor pressure, and fluorescence measurements. Ion-association behavior of C 10 MeImCl (aq) in the premicellar regime has also been investigated. Two characteristic concentrations, namely the critical aggregation concentration and polymer saturation concentration, before free C 10 MeImCl micelles appear in C 10 MeImCl–NaCMC solutions were identified. Effects of temperature, NaCMC concentration, and the bulk solution structural property on the self-aggregation of C 10 MeImCl have been discussed to elucidate C 10 MeImCl–NaCMC interactions. Thermodynamics of the micellization processes provided important insight regarding the (a) release of water molecules from the hydration layer around the hydrophilic domain, and from the water cage around the hydrophobic moiety of the SAIL, and (b) transfer of the hydrocarbon chains into the micelle and restoration of the H-bonding structure of the water around the micelle.

The influence of alloying on mode-selectivity in H2O dissociation on Cu/Ni bimetallic surfaces has been studied using a fully quantum approach based on reaction path Hamiltonian. Both the metal alloy catalyst surface and the normal modes of H2O impact the chemical reactivity of H2O dissociation. A combination of these two different factors will enhance their influence reasonably. Among all the bimetallic surfaces, one monolayer (Ni4_Cu(111)) and 12 monolayer of Ni on Cu surface (Ni2_Cu(111)) show lowest barrier to the dissociation. Excitation of bending mode and symmetric stretching mode enhances the reactivity remarkably due to a significant decrease in their frequencies near the transition state in the vibrational adiabatic approximation. In the presence of non-adiabatic coupling between the modes, asymmetric stretching also shows similar enhancement in reactivity to that of symmetric stretching for all the systems. Inclusion of lattice motion using a sudden model enhances the dissociation probability at surface temperature 300 K and at lower incident energy, compared to that of the static surface approximation. The mode selective behaviour of H2O molecules is almost similar on all the Cu- and Ni-based surfaces. The excitation of symmetric stretching vibration by one quantum is shown to have largest efficacy for promoting reactions for all the systems. Overall, the dissociation probabilities for all the systems are enhanced by vibrational excitation of normal modes and become more significant with the non-adiabatic coupling effect.

The influence of sodium carboxymethylcellulose (NaCMC) on the aggregation phenomena of a surface active ionic liquid 1-hexadecyl-3-methylimidazolium chloride (C 16 MeImCl) was studied in aqueous solutions using electrical conductivity and surface tension measurements. The counterion condensation behavior of NaCMC (aq) and the premicellar ion-association behavior of C 16 MeImCl (aq) were also investigated. Two characteristic concentrations, namely the critical aggregation concentration and polymer saturation concentration, before free C 16 MeImCl micelles appear in C 16 MeImCl–NaCMC solutions have been identified. Effects of temperature, NaCMC concentration, and the charge density parameter of NaCMC on the self-aggregation of the C 16 MeImCl have been discussed to elucidate C 16 MeImCl–NaCMC interactions. The thermodynamic parameters for micellization of C 16 MeImCl were estimated both in absence and in the presence of NaCMC. The observed enthalpy–entropy compensation effect in C 16 MeImCl and C 16 MeImCl–NaCMC systems provided important insight as to how micellization processes are governed by the bulk structural property of the solution with respect to that of the water.