Session Y26: Predicting Rare Event Kinetics in Complex Systems with Theory, Simulations and Machine Learning IV

Live

I will describe efforts to understand rare event processes of ions in solution. We employ a range of theoretical techniques comparing full atomistic molecular dynamics simulation for generation of ensembles and dynamic response to reduced models of a Generalized Langevin Equation form. We explore alternatives to conventional distance reaction coordinates such as electric fields, coordination numbers, Google Page Rank, and vibrational frequencies. Extensions of concepts invoked in Marcus theory will be described. Such analysis provides a framework for enhancing understanding of phenomena in response to inhomogeneities, coupling of multiple time scales, and energy transfer. Such techniques show promise for elucidating rare-event kinetics in complex systems.   

Live

Drug-target binding kinetics, as reflected through on-target residence time, is increasingly recognized as an important driver of in vivo drug efficacy and safety. However, due to long time scales and path dependency, computational estimation of residence times is still very challenging and as such has not yet been routinely incorporated into drug discovery project execution. Here we present a novel enhanced sampling based approach to computing residence times which combines automated selection of collective variables and ligand exit pathways with infrequent metadynamics to estimate the off rates. We report preliminary results including predicted residence times for over two dozen ligands across three pharmaceutically relevant protein targets, which are in good agreement with experimental results. The method is scalable, fully automated, and computationally efficient, with turnaround times comparable to free energy calculations routinely used to estimate drug binding affinities - thus allowing it to be readily employed in active drug discovery projects.   

Live

Understanding nucleation of crystals from solution is crucial in calculating nucleation rates and designing crystallization processes to control polymorphism. In this work, we examine the non-classical nucleation mechanism that has been reported in nucleation of NaCl from highly concentrated aqueous solutions [1]. By obtaining the free energy change of nucleation as a function of two structure specific nucleus size coordinates, we demonstrate the thermodynamic preference for nucleation through a composite cluster, where the crystalline nucleus is surrounded by an amorphous layer. The composite cluster model of Iwamatsu [2] describes the free energy change of cluster formation very well. The thermodynamic properties extracted from the fit show that the amorphous phase is less stable than the metastable solution phase, whereas the crystalline phase is the most stable. Finally, we calculate the nucleation rate of the crystalline cluster using our 2D system and compare it with nucleation rates obtained using a 1D free energy profile. This study is preparatory to an investigation of polymorphism in glycine crystallization.

[1] Jiang et al., J. Chem. Phys. 150, no. 12 (2019): 124502.
[2] Iwamatsu, M. J. Chem. Phys., 134 no. 16 (2011): 164508.   

Live

In-depth understanding of the kinetics and mechanisms of rare events in complex systems requires robust sampling of the high-dimensional phase space and the exploration of associated free energy surfaces. In this work, we combine enhanced sampling techniques, such as driven adiabatic free energy dynamics and metadynamics, with a path collective variable defined in a global classifier space. The global classifiers are determined based on local structural environments that are identified using a classification neural network. We demonstrate that the proposed scheme can efficiently sample transformation between different crystalline phases in metallic tungsten and reproduce the free energy landscape.   

Live

Enhanced sampling techniques generate trajectories at a biased potential, such that rare events occur more frequently. Path reweighing techniques recover the transition rates of the unbiased system from the biased trajectories by calculating the path probability ratio. Path reweighing requires that (a) the trajectory has been generated using an integration scheme for stochastic dynamics, and (b) that the formula for the path probability ratio has been tailored for that specific integration scheme. This makes them technically difficult, because a separate reweighing factor for each stochastic integration scheme is needed.

Most published path probability ratios are derived for overdamped Langevin dynamics. Yet, overdamped Langevin dynamics is rarely used in MD simulations. Instead a variety of integration schemes for Langevin dynamics are used. Simply applying the path probability ratio for overdamped Langevin dynamics to a Langevin trajectory introduces a sizeable error.

Here we derive the path probability ratio for the integration scheme of Langevin dynamics implemented in OpenMM. We demonstrate that it to accurately reweights Langevin trajectories.

By comparing this path probability ratio to the path probability ratio for overdamped Langevin dynamics, we then derive an approximate and general path probability ratio for Langevin dynamics. This approximate path probability ratio depends on the random number sequence used to generate the Langevin trajectory and on the bias potential. Because it is independent of the integration scheme, it removes requirement (b) and can be used as multi-purpose probability ratio for any Langevin trajectory. We show that the approximate path probability ratio yields highly accurate results, and discuss the limits of the approximation.

L Donati, BG Keller, J Chem Phys. 149, 072335 (2018)
S Kieninger, L Donati, BG Keller, Curr Opin Struct Biol, 61, 124-131 (2020)   

Live

A variety of enhanced statistical and numerical methods are now routinely used to extract comprehensible and relevant thermodynamic information from the vast amount of complex, high-dimensional data obtained from intensive molecular simulations. The characterization of kinetic properties, such as diffusion coefficients, of molecular systems with significantly high energy barriers, on the other hand, has received less attention. Among others, Markov state models, in which the long-time statistical dynamics of a system is approximated by a Markov chain on a discrete partition of configuration space, have seen widespread use in recent years, with the aim of tackling these fundamental issues. Here, we discuss a general, automatic method to assess multidimensional position-dependent diffusion coefficients within the framework of Markovian stochastic processes and Kramers-Moyal expansion. We apply the formalism to one- and two-dimensional analytic potentials and data from explicit solvent molecular dynamics simulations, including the water-mediated conformations of alanine dipeptide and the transport of drug molecule across three-dimensional heterogeneous porous media.   

Live

In drug manufacturing, when solvent-based methods are used for the crystallization of active pharmaceutical ingredients (APIs), often, the solvent can weakly interact with the API resulting in the formation of a new solid form, the so-called solvate. Thus, when desolvation occurs upon heating it can result in either a known stable form or the recrystallization to a new solid form. In this work, we researched the desolvation kinetics of the fluorobenzene (FB) solvates of Ibrutinib (IBR) by combining thermogravimetric analysis (TGA), in-situ powder X-ray diffraction (PXRD), all-atom molecular dynamics (MD) simulations, and macroscopic modelling of the desolvation kinetics. Using model-fitting and iso-conversional methods, we accurately predicted the desolvation kinetics by validating both TGA and XRPD data as well as calculated the activation energy of desolvation. Performing a large set of MD simulations, we traced individual FB molecules and then calculated the activation energy for their diffusion. Our results show that the desolvation kinetics is rather affected by the diffusion of FB in the crystal lattice than by polymorphic transformations.   

Live

Simulating the dynamics of complex systems requires the precise identification of atomic diffusions, including rare-events, through an exhaustive exploration of the energy landscape. To reduce the associated computational cost we propose to use the Static Mode (SM) method [Ref1] together with ab initio calculations (QM). SM are first calculated from the Hessian matrix associated to an initial relaxed configuration. Each SM represents the strain field of a molecular system submitted to the stress (displacement) of a given atom in a given direction. Using this information, we screen the perturbation of all atoms in all directions and evaluate the local response of the system. The most relevant deformations are then selected in regards to a criterion chosen by the user to guide the evolution of the system. Finally, from the deformed states, ab initio calculations are used to reach new minima. The QMSM method that we propose can be applied to any type of system. Its validation will be presented, based on two different applications: how DNA can graft on alumina surface and how defect can migrate in silicon bulk.

Ref1 : M. Brut et al. European Physical Journal E vol. 28,17-25 (2009)   

Live

We hereby report on the use of Restrained Molecular Dynamics [1] to investigate the effect of the presence of small concentrations of hydrophobic gases on the phase behaviour of water confined in hydrophobic nanopores [2].

Metastability, and thus rare events, are a typical feature of water in extreme hydrophobic confinement: nanopores can either be wet or occupied by vapor with the relative stability of the two states depending on the size and hydrophobicity of the porous material [3].

Our results show that the presence of even a single hydrophobic gas atom within the pore can dramatically accelerate evaporation and alter the relative stability of the wet and vapor-filled states.

These findings might shed a light on a possible unspecific mechanism of action of volatile anaesthetic gases on certain ion channels, whose ionic currents may be block by vapour bubbles, and on the control of undesired bubble formation in artificial nanopores for sensing and chromatography.

[1] Maragliano, L., & Vanden-Eijnden, E., (2006). Chem.Phys.Letters 426.1-3
[2] Camisasca, G., Tinti, A., & Giacomello, A. (2020). JPCL 11
[3] Tinti, A., Giacomello, A., Grosu, Y., & Casciola C.M. (2017). PNAS 114.48   

Live

We study first-passage time (FPT) distributions of target finding events through complex environments with elongated obstacles distributed with different anisotropies and volume occupation fractions. For isotropic systems and for low densities of aligned obstacles, FPTs are exponentially distributed. At large enough densities of aligned obstructions, elongated channels emerge, gradually reducing the dynamics dimensionality from 3D, to 1D in the case of narrow structures. We analyze how the local structure of the channels, such as geometry and size, modifies the FPT distribution. We find that channel size and geometry determines the shape of the FPT distribution and its mean first-passage time. Moreover, we develop an exactly solvable model for synthetic rectangular channels that captures the effects of the tortuous local structure of the channels that naturally emerge in our system. For arbitrary values of the nematic order parameter of fiber orientations, we develop a mapping to the simpler situation of fully aligned fibers at some other effective volume occupation fraction. Because of the complex nature of fibrous biological environments in tissues, we suggest that our results shed light on the understanding of the molecular transport occurring between cells.   

Live

In the assembly of nanoscopic and mesoscopic materials, clusters form as important precursors to larger aggregates. Long lived stable, and metastable states within these assemblies determine the structure and dynamics of subsequent assembly. In metallic nanoparticles, specific cluster geometries are sought to control the particle's catalytic properties. Predicting long-lived aggregates is a complex problem, and conventional structural analyses based on spectroscopy or diffraction provide only averaged and not instantaneous structures. Molecular simulations offer an opportunity to examine the formation and fluctuation of metallic clusters in a highly controlled environment where the preferred conformations can be determined. In this work, we use advanced sampling algorithms coupled with ab initio molecular dynamics to estimate cluster conformations' free energy. We explore the conformational free energy landscape of the small metallic cluster Au₄, a simplified system enabling comprehensive study in neutral and charged configurations. We analyze the thermodynamics of conformational isomerization and predict transition rates with transition state theory. These simulations offer a quantitative understanding of the fluxionality of cluster structures.