Search

We obtain asymptotically sharp error estimates for the consistency error of the Target Measure Diffusion map (TMDmap) (Banisch et al. 2020), a variant of diffusion maps featuring importance sampling and hence allowing input data drawn from an arbitrary density. The derived error estimates include the bias error and the variance error. The resulting convergence rates are consistent with the approximation theory of graph Laplacians. The key novelty of our results lies in the explicit quantification of all the prefactors on leading-order terms. We also prove an error estimate for solutions of Dirichlet BVPs obtained using TMDmap, showing that the solution error is controlled by consistency error. We use these results to study an important application of TMDmap in the analysis of rare events in systems governed by overdamped Langevin dynamics using the framework of transition path theory (TPT). The cornerstone ingredient of TPT is the solution of the committor problem, a boundary value problem for the backward Kolmogorov PDE. Remarkably, we find that the TMDmap algorithm is particularly suited as a meshless solver to the committor problem due to the cancellation of several error terms in the prefactor formula. Furthermore, significant improvements in bias and variance errors occur when using a quasi-uniform sampling density. Our numerical experiments show that these improvements in accuracy are realizable in practice when using $\delta$-nets as spatially uniform inputs to the TMDmap algorithm.

Transition path theory (TPT) is a mathematical framework for quantifying rare transition events between a pair of selected metastable states $A$ and $B$. Central to TPT is the committor function, which describes the probability to hit the metastable state $B$ prior to $A$ from any given starting point of the phase space. Once the committor is computed, the transition channels and the transition rate can be readily found. The committor is the solution to the backward Kolmogorov equation with appropriate boundary conditions. However, solving it is a challenging task in high dimensions due to the need to mesh a whole region of the ambient space. In this work, we explore the finite expression method (FEX, Liang and Yang (2022)) as a tool for computing the committor. FEX approximates the committor by an algebraic expression involving a fixed finite number of nonlinear functions and binary arithmetic operations. The optimal nonlinear functions, the binary operations, and the numerical coefficients in the expression template are found via reinforcement learning. The FEX-based committor solver is tested on several high-dimensional benchmark problems. It gives comparable or better results than neural network-based solvers. Most importantly, FEX is capable of correctly identifying the algebraic structure of the solution which allows one to reduce the committor problem to a low-dimensional one and find the committor with any desired accuracy.

An experimental arrangement and a set of experiments are developed to generate empirical evidence of the effect of noise on a rotating, macro-scale cantilever structure. The experiment is a controlled representation of a rotating machinery blade. Due to the nature of the nonlinear restoring forces acting on the cantilever structure, the structure's response includes regions of multi-stability and hysteresis. Here, a large number of trials are used to show that random perturbations can be used to create a transition between a high amplitude response and a low amplitude response of the cantilever. The observed transition behavior occurs from a high amplitude response to a low amplitude response, but not vice versa. Stochastic modeling of the system, Monte Carlo simulations, and calculations of the stochastic system's quasipotential are used to explain the nearly one-directional transition behavior. These noise-influenced transitions can also occur in other physical systems., Comment: 25 pages, 22 figures

Many processes in nature such as conformal changes in biomolecules and clusters of interacting particles, genetic switches, mechanical or electromechanical oscillators with added noise, and many others are modeled using stochastic differential equations with small white noise. The study of rare transitions between metastable states in such systems is of great interest and importance. The direct simulation of rare transitions is difficult due to long waiting times. Transition path theory is a mathematical framework for the quantitative description of rare events. Its crucial component is the committor function, the solution to a boundary value problem for the backward Kolmogorov equation. The key fact exploited in this work is that the optimal controller constructed from the committor leads to the generation of transition trajectories exclusively. We prove this fact for a broad class of stochastic differential equations. Moreover, we demonstrate that the committor computed for a dimensionally reduced system and then lifted to the original phase space still allows us to construct an effective controller and estimate the transition rate with reasonable accuracy. Furthermore, we propose an all-the-way-through scheme for computing the committor via neural networks, sampling the transition trajectories, and estimating the transition rate without meshing the space. We apply the proposed methodology to four test problems: the overdamped Langevin dynamics with Mueller's potential and the rugged Mueller potential in 10D, the noisy bistable Duffing oscillator, and Lennard-Jones-7 in 2D.

A polymer network fractures by breaking covalent bonds, but the experimentally measured strength of the polymer network is orders of magnitude lower than the strength of covalent bonds. We investigate the effect of statistical variation of the number of links in polymer chains on strength using a parallel chain model. Each polymer chain is represented by a freely-jointed chain, with a characteristic J-shaped force-extension curve. The chain carries entropic forces for most of the extension and carries covalent forces only for a narrow range of extension. The entropic forces are orders of magnitude lower than the covalent forces. Chains with a statistical distribution of the number of links per chain are pulled between two rigid parallel plates. Chains with fewer links attain covalent forces and rupture at smaller extensions, while chains with more links still carry entropic forces. We compute the applied force on the rigid plates as a function of extension and define the strength of the parallel chain model by the maximum force divided by the total number of chains. With the J-shaped force-extension curve of each chain, even a small scatter in the number of links per chain greatly reduces the strength of the parallel chain model. We further show that the strength of the parallel chain model relates to the scatter in the number of links per chain according to a power law., Comment: 17 pages, 5 figures

We present algorithms for solving high-frequency acoustic scattering problems in complex domains. The eikonal and transport partial differential equations from the WKB/geometric optic approximation of the Helmholtz equation are solved recursively to generate boundary conditions for a tree of eikonal/transport equation pairs, describing the phase and amplitude of a geometric optic wave propagating in a complicated domain, including reflection and diffraction. Edge diffraction is modeled using the uniform theory of diffraction. For simplicity, we limit our attention to domains with piecewise linear boundaries and a constant speed of sound. The domain is discretized into a conforming tetrahedron mesh. For the eikonal equation, we extend the jet marching method to tetrahedron meshes. Hermite interpolation enables second order accuracy for the eikonal and its gradient and first order accuracy for its Hessian, computed using cell averaging. To march the eikonal on an unstructured mesh, we introduce a new method of rejecting unphysical updates by considering Lagrange multipliers and local visibility. To handle accuracy degradation near caustics, we introduce several fast Lagrangian initialization algorithms. We store the dynamic programming plan uncovered by the marcher in order to propagate auxiliary quantities along characteristics. We introduce an approximate origin function which is computed using the dynamic programming plan, and whose 1/2-level set approximates the geometric optic shadow and reflection boundaries. We also use it to propagate geometric spreading factors and unit tangent vector fields needed to compute the amplitude and evaluate the high-frequency edge diffraction coefficient. We conduct numerical tests on a semi-infinite planar wedge to evaluate the accuracy of our method. We also show an example with a more realistic building model with challenging architectural features., Comment: 47 pages, 17 figures

The study of phenomena such as protein folding and conformational changes in molecules is a central theme in chemical physics. Molecular dynamics (MD) simulation is the primary tool for the study of transition processes in biomolecules, but it is hampered by a huge timescale gap between the processes of interest and atomic vibrations which dictate the time step size. Therefore, it is imperative to combine MD simulations with other techniques in order to quantify the transition processes taking place on large timescales. In this work, the diffusion map with Mahalanobis kernel, a meshless approach for approximating the Backward Kolmogorov Operator (BKO) in collective variables, is upgraded to incorporate standard enhanced sampling techniques such as metadynamics. The resulting algorithm, which we call the "target measure Mahalanobis diffusion map" (tm-mmap), is suitable for a moderate number of collective variables in which one can approximate the diffusion tensor and free energy. Imposing appropriate boundary conditions allows use of the approximated BKO to solve for the committor function and utilization of transition path theory to find the reactive current delineating the transition channels and the transition rate. The proposed algorithm, tm-mmap, is tested on the two-dimensional Moro-Cardin two-well system with position-dependent diffusion coefficient and on alanine dipeptide in two collective variables where the committor, the reactive current, and the transition rate are compared to those computed by the finite element method (FEM). Finally, tm-mmap is applied to alanine dipeptide in four collective variables where the use of finite elements is infeasible., Comment: Restructured introduction, improved explanation of key algorithms and formulas (Section II.C and III.B,C). Streamlined presentation and proof of Theorem 1

Hydrocarbon pyrolysis is a complex process involving large numbers of chemical species and types of chemical reactions. Its quantitative description is important for planetary sciences, in particular, for understanding the processes occurring in the interior of icy planets, such as Uranus and Neptune, where small hydrocarbons are subjected to high temperature and pressure. We propose a computationally cheap methodology based on an originally developed ten-reaction model, and the configurational model from random graph theory. This methodology yields to accurate predictions for molecule size distributions for a variety of initial chemical compositions and temperatures ranging from 3200K to 5000K. Specifically, we show that the size distribution of small molecules is particularly well predicted, and the size of the largest molecule can be accurately predicted provided that it is not too large., Comment: 27 pages (11 main + 16 supplementary), 15 figures (5 main + 10 supplementary). Submitted to the Journal of Physical Chemistry A

We present a new algorithm, the efficient jet marching method (EJM), for computing the quasipotential and its gradient for two-dimensional SDEs. The quasipotential is a potential-like function for nongradient SDEs that gives asymptotic estimates for the invariant probability measure, expected escape times from basins of attractors, and maximum likelihood escape paths. The quasipotential is a solution to an optimal control problem with an anisotropic cost function which can be solved for numerically via Dijkstra-like label-setting methods. Previous Dijkstra-like quasipotential solvers have displayed in general 1st order accuracy in the mesh spacing. However, by utilizing higher order interpolations of the quasipotential as well as more accurate approximations of the minimum action paths (MAPs), EJM achieves second-order accuracy for the quasipotential and nearly second-order for its gradient. Moreover, by using targeted search neighborhoods for the fastest characteristics following the ideas of Mirebeau, EJM also enjoys a reduction in computation time. This highly accurate solver enables us to compute the prefactor for the WKB approximation for the invariant probability measure and the Bouchet-Reygner sharp estimate for the expected escape time for the Maier-Stein SDE. Our codes are available on GitHub., Comment: 21 figures, 3 tables

The study of rare events in molecular and atomic systems such as conformal changes and cluster rearrangements has been one of the most important research themes in chemical physics. Key challenges are associated with long waiting times rendering molecular simulations inefficient, high dimensionality impeding the use of PDE-based approaches, and the complexity or breadth of transition processes limiting the predictive power of asymptotic methods. Diffusion maps are promising algorithms to avoid or mitigate all these issues. We adapt the diffusion map with Mahalanobis kernel proposed by Singer and Coifman (2008) for the SDE describing molecular dynamics in collective variables in which the diffusion matrix is position-dependent and, unlike the case considered by Singer and Coifman, is not associated with a diffeomorphism. We offer an elementary proof showing that one can approximate the generator for this SDE discretized to a point cloud via the Mahalanobis diffusion map. We use it to calculate the committor functions in collective variables for two benchmark systems: alanine dipeptide, and Lennard-Jones-7 in 2D. For validating our committor results, we compare our committor functions to the finite-difference solution or by conducting a "committor analysis" as used by molecular dynamics practitioners. We contrast the outputs of the Mahalanobis diffusion map with those of the standard diffusion map with isotropic kernel and show that the former gives significantly more accurate estimates for the committors than the latter., Comment: Restructured introduction, additional Theorem 3.1 and Appendix A, B

We develop a family of compact high-order semi-Lagrangian label-setting methods for solving the eikonal equation. These solvers march the total 1-jet of the eikonal, and use Hermite interpolation to approximate the eikonal and parametrize characteristics locally for each semi-Lagrangian update. We describe solvers on unstructured meshes in any dimension, and conduct numerical experiments on regular grids in two dimensions. Our results show that these solvers yield at least second-order convergence, and, in special cases such as a linear speed of sound, third-order of convergence for both the eikonal and its gradient. We additionally show how to march the second partials of the eikonal using cell-based interpolants. Second derivative information computed this way is frequently second-order accurate, suitable for locally solving the transport equation. This provides a means of marching the prefactor coming from the WKB approximation of the Helmholtz equation. These solvers are designed specifically for computing a high-frequency approximation of the Helmholtz equation in a complicated environment with a slowly varying speed of sound, and, to the best of our knowledge, are the first solvers with these properties. We provide a link to a package online providing the solvers, and from which the results of this paper can be reproduced easily., Comment: 26 pages, 11 figures, 2 tables

The study of noise-driven transitions occurring rarely on the time-scale of systems modeled by SDEs is of crucial importance for understanding such phenomena as genetic switches in living organisms and magnetization switches of the Earth. For a gradient SDE, the predictions for transition times and paths between its metastable states are done using the potential function. For a nongradient SDE, one needs to decompose its forcing into a gradient of the so-called quasipotential and a rotational component, which cannot be done analytically in general. We propose a methodology for computing the quasipotential for highly dissipative and chaotic systems built on the example of Lorenz'63 with an added stochastic term. It is based on the ordered line integral method, a Dijkstra-like quasipotential solver, and combines 3D computations in whole regions, a dimensional reduction technique, and 2D computations on radial meshes on manifolds or their unions. Our collection of source codes is available on M. Cameron's web page and on GitHub., Comment: 38 pages, 17 figures, 5 3D figures are linked to movies posted on youtube, software is available on GitHub

Nongradient SDEs with small white noise often arise when modeling biological and ecological time-irreversible processes. If the governing SDE were gradient, the maximum likelihood transition paths, transition rates, expected exit times, and the invariant probability distribution would be given in terms of its potential function. The quasipotential plays a similar role for nongradient SDEs. Unfortunately, the quasipotential is the solution of a functional minimization problem that can be obtained analytically only in some special cases. We propose a Dijkstra-like solver for computing the quasipotential on regular rectangular meshes in 3D. This solver results from a promotion and an upgrade of the previously introduced ordered line integral method with the midpoint quadrature rule for 2D SDEs. The key innovations that have allowed us to keep the CPU times reasonable while maintaining good accuracy are $(i)$ a new hierarchical update strategy, $(ii)$ the use of Karush-Kuhn-Tucker theory for rejecting unnecessary simplex updates, and $(iii)$ pruning the number of admissible simplexes and a fast search for them. An extensive numerical study is conducted on a series of linear and nonlinear examples where the quasipotential is analytically available or can be found at transition states by other methods. In particular, the proposed solver is applied to Tao's examples where the transition states are hyperbolic periodic orbits, and to a genetic switch model by Lv et al. (2014). The C source code implementing the proposed algorithm is available at M. Cameron's web page., Comment: 11 figures, 3 tables

Nongradient stochastic differential equations (SDEs) with position-dependent and anisotropic diffusion are often used in biological modeling. The quasi-potential is a crucial function in the Large Deviation Theory that allows one to estimate transition rates between attractors of the corresponding ordinary differential equation and find the maximum likelihood transition paths. Unfortunately, the quasi-potential can rarely be found analytically. It is defined as the solution to a certain action minimization problem. In this work, the recently introduced Ordered Line Integral Method (OLIM) is extended for computing the quasi-potential for 2D SDEs with anisotropic and position-dependent diffusion scaled by a small parameter on a regular rectangular mesh. The presented solver employs the dynamical programming principle. At each step, a local action minimization problem is solved using straight line path segments and the midpoint quadrature rule. The solver is tested on two examples where analytic formulas for the quasi-potential are available. The dependence of the computational error on the mesh size, the update factor K (a key parameter of OLIMs), as well as the degree and the orientation of anisotropy is established. The effect of anisotropy on the quasi-potential and the maximum likelihood paths is demonstrated on the Maier-Stein model. The proposed solver is applied to find the quasi-potential and the maximum likelihood transition paths in a model of the genetic switch in Lambda Phage between the lysogenic state where the phage reproduces inside the infected cell without killing it, and the lytic state where the phage destroys the infected cell., Comment: 26 pages, 6 figures, 3 tables

The quasi-potential is a key function in the Large Deviation Theory. It characterizes the difficulty of the escape from the neighborhood of an attractor of a stochastic non-gradient dynamical system due to the influence of small white noise. It also gives an estimate of the invariant probability distribution in the neighborhood of the attractor up { to} the exponential order. We present a new family of methods for computing the quasi-potential on a regular mesh named the Ordered Line Integral Methods (OLIMs). In comparison with the first proposed quasi-potential finder based on the Ordered Upwind Method (OUM) (Cameron, 2012), the new methods are 1.5 to 4 times faster, can produce error two to three orders of magnitude smaller, and may exhibit faster convergence. Similar to the OUM, OLIMs employ the dynamical programming principle. Contrary to it, they (i) have an optimized strategy for the use of computationally expensive { triangle} updates leading to a notable speed-up, and (ii) directly solve local minimization problems using quadrature rules instead of solving the corresponding Hamilton-Jacobi-type equation by the first order finite difference upwind scheme. The OLIM with the right-hand quadrature rule is equivalent to OUM. The use of higher order quadrature rules in local minimization problems dramatically boosts up the accuracy of OLIMs. We offer a detailed discussion on the origin of numerical errors in OLIMs and propose rules-of-thumb for the choice of the important parameter, the update factor, in the OUM and OLIMs. Our results are supported by extensive numerical tests on two challenging 2D examples., Comment: 13 figures, 4 tables

We model an isothermal aggregation process of particles/atoms interacting according to the Lennard-Jones pair potential by mapping the energy landscapes of each cluster size $N$ onto stochastic networks, computing transition probabilities {from} the network for an $N$-particle cluster to the one for $N+1$, and connecting these networks into a single joint network. The attachment rate is a control parameter. The resulting network representing the aggregation of up to 14 particles contains {6427} vertices. It is not only time-irreversible but also reducible. To analyze its transient dynamics, we introduce the sequence of the expected initial and pre-attachment distributions and compute them for a wide range of attachment rates and three values of temperature. As a result, we find the {configurations most likely to be observed} in the process of aggregation for each cluster size. We examine the attachment process and conduct a structural analysis of the sets of local energy minima for every cluster size. We show that both processes taking place in the network, attachment and relaxation, lead to the dominance of icosahedral packing in small (up to 14 atom) clusters., Comment: 28 pages, 13 figures, 6 tables

Large continuous-time Markov chains with exponentially small transition rates arise in modeling complex systems in physics, chemistry and biology. We propose a constructive graph-algorithmic approach to determine the sequence of critical timescales at which the qualitative behavior of a given Markov chain changes, and give an effective description of the dynamics on each of them. This approach is valid for both time-reversible and time-irreversible Markov processes, with or without symmetry. Central to this approach are two graph algorithms, Algorithm 1 and Algorithm 2, for obtaining the sequences of the critical timescales and the hierarchies of Typical Transition Graphs or T-graphs indicating the most likely transitions in the system {without and with} symmetry respectively. The sequence of {critical} timescales includes the subsequence of the reciprocals { of the real parts } of eigenvalues. Under a certain assumption, we prove sharp asymptotic estimates for eigenvalues (including prefactors) and show how one can extract them from the output of Algorithm 1. We discuss the relationship between Algorithms 1 and 2, and explain how one needs to interpret the output of Algorithm 1 if it is applied in the case with symmetry instead of Algorithm 2. Finally, we analyze an example motivated by R. D. Astumian's model of the dynamics of kinesin, a molecular motor, by means of Algorithm 2., Comment: 46 pages, 17 figures

We consider stochastic networks with pairwise transition rates of the exponential form where the temperature T is a small parameter. Such networks arise in physics and chemistry and serve as mathematically tractable models of complex systems. Typically, such networks contain large numbers of states and widely varying pairwise transition rates. We present a methodology for spectral analysis and clustering of such networks that takes advance of the small parameter T and consists of two steps: (1) computing zero-temperature asymptotics for eigenvalues and the collection of quasi-invariant sets, and (2) finite temperature continuation. Step (1) is re- ducible to a sequence of optimization problems on graphs. A novel single-sweep algorithm for solving them is introduced. Its mathematical justification is provided. This algorithm is valid for both time-reversible and time-irreversible networks. For time-reversible networks, a finite temperature continuation technique combining lumping and truncation with Rayleigh quotient iteration is developed. The proposed methodology is applied to the network representing the energy landscape of the Lennard-Jones-75 cluster containing 169,523 states and 226,377 edges. The transition process between its two major funnels, is analyzed. The corresponding eigenvalue is shown to have a kink at the solid-solid phase transition temperature., Comment: 33 pages, 16 figures

QPot is an R package for analyzing two-dimensional systems of stochastic differential equations. It provides users with a wide range of tools to simulate, analyze, and visualize the dynamics of these systems. One of QPot's key features is the computation of the quasi-potential, an important tool for studying stochastic systems. Quasi-potentials are particularly useful for comparing the relative stabilities of equilibria in systems with alternative stable states. This paper describes QPot's primary functions, and explains how quasi-potentials can yield insights about the dynamics of stochastic systems. Three worked examples guide users through the application of QPot's functions.

The concept of metastability has caused a lot of interest in recent years. The spectral decomposition of the generator matrix of a stochastic network exposes all of the transition processes in the system. The assumption of the existence of a low lying group of eigenvalues separated by a spectral gap, leading to factorization of the dynamics, has become a popular theme. We consider stochastic networks representing potential energy landscapes where the states and the edges correspond to local minima and transition states respectively, and the pairwise transition rates are given by the Arrhenuis formula. Using the minimal spanning tree, we construct the asymptotics for eigenvalues and eigenvectors of the generator matrix starting from the low lying group. This construction gives rise to an efficient algorithm for computing the asymptotic spectrum suitable for large and complex networks. We apply it to Wales's Lennard-Jones-38 network with 71887 states and 119853 edges where the underlying potential energy landscape has a double-funnel structure. Our results demonstrate that the concept of metastability should be applied with care to this system. In particular, for the full network, there is no significant spectral gap separating the eigenvalue corresponding to the exit from the wider and shallower icosahedral funnel at any reasonable temperature range., Comment: Submitted to Journal Networks and Heterogeneous Media on Feb. 25, 2014, 36 pages, 14 figures