Metadynamics for transition paths in irreversible dynamics

Stochastic systems often exhibit multiple viable metastable states that are long-lived. Over very long timescales, fluctuations may push the system to transition between them, drastically changing its macroscopic configuration. In realistic systems, these transitions can happen via multiple physical mechanisms, corresponding to multiple distinct transition channels for a pair of states. In this paper, we use the fact that the transition path ensemble is equivalent to the invariant measure of a gradient flow in pathspace, which can be efficiently sampled via metadynamics. We demonstrate how this pathspace metadynamics, previously restricted to reversible molecular dynamics, is in fact very generally applicable to metastable stochastic systems, including irreversible and time-dependent ones, and allows to estimate rigorously the relative probability of competing transition paths. We showcase this approach on the study of a stochastic partial differential equation describing magnetic field reversal in the presence of advection.


I
Rare events in stochastic systems lie at the core of many applications in the sciences, describing phenomena as wide as regime changes and tipping points in Earth's climate [1][2][3][4], conformation changes of biomolecules and protein folding [5,6], material sciences [7], genetic switches [8,9], etc.These applications have in common that the underlying stochastic system spends very long times close to one macroscopic state, only to eventually switch to a completely different regime due to small and random fluctuations.From the perspective of stochastic analysis, these rare but important transitions correspond to a stochastic bridge, i.e. a realization of the process that starts in one but ends in another basin of attraction.A whole class of rare event algorithms have been developed to obtain information about transition paths, from importance sampling, over genealogical particle algorithms to Freidlin-Wentzell theory and large deviations.They are valid in different regimes, and have different ranges of applicability (see [10] for an overview).
Due to the utmost importance of rare events in applications, several methods have been developed to enhance the probability of observing such events in computer simulations.Transition path sampling (TPS) [11][12][13] and string-based methods [14,15] allow focusing the computational effort in the simulation of the transition events, avoiding spending time in sampling the fluctuations of the system within the metastable states.Umbrella sampling [16][17][18], adaptive force bias [19][20][21][22] and metadynamics [23][24][25], allow reconstructing the free energy landscapes of metastable systems, which are characterized by the presence of at least two deep free energy minima.Free energy methods require a preliminary choice of a small set of collective variables (CVs) which are supposed at least to distinguish all the metastable states.It can be shown that if there are precisely two metastable states then the optimal CV is the committor function [26,27], and many methods have recently been developed using machine learning and data science techniques in order to deduce or approximate it [28][29][30][31].TPS, if used only to simulate the transitions, does not require choosing a CV, but some of its variants aimed at estimating the rate [32,33] also require using a CV.Most TPS-based techniques suffer from metastability in path space: if multiple competing transition mechanisms exist, the trajectories generated by TPS typically remain confined in a single reaction channel.To address this problem it has been suggested to perform metadynamics in path space [12,13,34].The idea, introduced in these works, is driving the transition paths towards unexplored reaction channels by a history-dependent bias defined as a function of a CV capable of distinguishing those channels.
Most of these approaches were developed specifically for applications in molecular and atomistic simulations.In these systems the dynamics (typically) satisfies detailed balance, and samples a stationary and currentless probability measure.However, many stochastic models describing important physical systems do not satisfy these properties and regularly detailed balance is violated [cite].The dynamics does not necessarily sample a stationary probability measure, and even when this is the case this probability measure is often not currentless.Stochastic Partial Differential Equations (SPDEs), which are used for example for describing the mesoscopic dynamics of fluids, active matter system, and atmosphere or ocean dynamics, often include explicitly time-dependent terms.In these condition, defining a free energy is not possible, and most of the methods for studying rare events cannot be applied in their standard form.
In this work we propose a principled approach which allows exploring the bridge ensemble for any kind of dynamics, including irreversible dynamics.The approach is based on the observation that the path probability density defined by the Onsager-Machlup-Jacobi (OMJ) functional defines a currentless probability measure even if the stochastic dynamics from which the action is derived does not [35,36].By performing Langevin dynamics in path space one can generate an ensemble of paths with a probability proportional to the exponential of the OMJ action.This allows defining a meaningful free energy as a function of virtually any collective variable, and use enhanced sampling methods to estimate this free energy.In particular, in order to explore the free energy landscape in path space we use metadynamics, following ref.[34].This approach allows studying rare events generated by irreversible and time-dependent stochastic processes and characterized by multiple competing transition mechanisms.In the context of irreversible dynamics, the advantage of this approach is that we can rely on the rigorously derived equivalence of the bridge measure and the invariant measure of certain stochastic partial differential equations, as proposed in [37][38][39] and formally generalized to arbitrary non-reversible systems [38].
In the following, we will first introduce into bridge path sampling in section 2.1, and subsequently describe the metadynamics algorithm generally in section 2.2.We then show how to apply metadynamics to transition paths to sample multiple competing transition mechanisms in section 2.3.Section 3 then demonstrates the applicability of the algorithm to a multitude of examples, from a temporally oscillating doublewell in section 3.1, to a twodimensional irreversible stochastic differential equation (SDE) in section 3.2, and finally in section 3.3 the stochastic Ginzburg-Landau equation with background advection as a fieldtheoretic example with spatial extent.We conclude with a discussion in section 4.

M P
2.1.Bridge Path Sampling and Transition Path Sampling.Consider for X t ∈ R n the stochastic differential equation where b ∶ R n ×[0, T ] → R n is the deterministic drift vector field, and W t is R n -dimensional Brownian motion.We are interested in drawing sample paths from (1), but conditional on a (possibly rare) final condition X T = x + .Such trajectories are called stochastic bridges, as they connect x − to x + via the stochastic process (1).Sampling such bridges, i.e. drawing paths from the bridge ensemble, is non-trivial, as simply integrating forward-in-time the stochastic equation ( 1) has a very low probability of arriving anywhere near x + , in particular in high dimension and/or if the final point is difficult to reach for the dynamics (i.e. if hitting near x + is rare).
Traditional examples for bridge sampling problems include molecular dynamics and chemical reactions, where we are interested in transitions between a reactant and product state of chemical molecules.Here, a well-developed theory of the form of transition state theory (TST) and transition path sampling (TPS) exists.Many methods employed in TPS, though, assume reversibility of the underlying stochastic process (with some notable exceptions, such as forward flux sampling [33]), which is indeed fulfilled in the classical molecular dynamics case.For equation (1), reversibility corresponds to demanding that the drift is of the particular form b(x) = −∇U (x) for a potential U ∶ R n → R. In that case, the stationary density is explicitly known to be the Gibbs distribution ρ(x) ∝ exp(−ε −1 U (x)).
For general drifts b, i.e. for the irreversible setting, bridge path sampling can be attempted as well, for example in the form of pathspace Langevin Markov Chain Monte Carlo (MCMC) [37][38][39], which at least formally allows a generalization to non-reversible systems [38].This is the path we will take in the following.
2.1.1.The Onsager-Machlup functional and large deviation theory.In the vanishing noise limit, ε → 0, transitions from x − to x + of (1) are described by the minimizers, also termed instantons, of the Freidlin-Wentzell action functional [40], ( This is because the probability of observing any given path φ(t) is given by P (φ(t)) ∼ exp(−ε −1 I[φ]), as can be intuited by "solving" the SDE (1) for the noise, which has a Gaussian distribution.Thus, considering all possible transition pathways φ(t) between x + and x − , the minimizer will exponentially dominate all other scenarios.In other words, for vanishing noise, the bulk of the transition trajectories between x − and x + will lie in the vicinity of the instanton trajectory.These minimizers can be computed efficiently by numerical means [10,[41][42][43][44][45].
For finite noise, this functional must be amended by an additional term [46] to obtain the Onsager-Machlup functional While finite noise minimizers of the Onsager-Machlup function (3) can be found, for example by performing a gradient descent in trajectory space, a key realization is that, similar to above, paths φ(t) of the process (1) for finite noise can formally be thought of as drawn from a Gibbs-type distribution As a consequence, one can sample the bridge ensemble by instead integrating the stochastic partial differential equation where φ(t, τ ) ∶ [0, T ] × [0, ∞] → R n is a field in physical time t and virtual time τ , with boundary conditions φ(0, τ ) = x − , φ(T, τ ) = x + .Here, η(t, τ ) is white in virtual and physical time, i.e.
Effectively, the physical time t is treated as spatial variable, and the evolution is happening in the newly introduced virtual time variable.Sampling the SPDE (5) amounts to Langevin Markov chain Monte Carlo (MCMC) in path space, i.e. the SPDE corresponds to a Markov process (namely the functional Langevin equation) with an invariant measure equivalent to the original transition path ensemble.
If the SDE ( 1) is in one dimension, n = 1, then combining (3) and ( 5) yields the SPDE which is essentially a gradient flow with noise, with the gradient computed via the Euler-Lagrange equation of ( 3).This result has been rigorously derived in [37][38][39].
If n > 1 and the system is non-gradient, the Euler-Lagrange equation of the the action (3) instead yields where the additional term proportional to ∂ t φ vanishes in the case of detailed balance, where ∇b(φ) is self-adjoint.This SPDE has been conjectured in [38, sec. 9], and also used in [47].Indeed, ( 5) can be used to effectively sample the bridge ensemble, since every generated path automatically fulfills the initial and final conditions of the bridge and further the SPDE has the bridge ensemble as invariant measure.The problem that motivates this work is that the Onsager-Machlup action (3) itself might have multiple local minima in trajectory space, each corresponding to a locally effective transition trajectory.This is the case, for example, when there are multiple competing reaction channels, or multiple physical mechanisms for the transition to occur.In that case, bridge sampling via (5) will be inefficient, since the SPDE will spend exponentially long times in one of the minima, and only very rarely explore others.We here propose to overcome this limitation by introducing metadynamics to the sampling of (5).

Metadynamics.
Metadynamics is a rare event simulation technique that allows effective estimation of the free energy even in the presence of multiple local minima with large barriers.It was originally developed, and is mostly used, to sample high-dimensional atomistic systems as present in molecular dynamics, with applications in e.g.material science or biophysics.
Its main idea is to avoid a Markov Chain Monte Carlo sampler to get stuck in a single maximum of the Gibbs measure, by gradually 'filling' the current free energy minimum until saturation, so that the system eventually starts exploring all other states.This filling is done in a controlled way that later allows to reconstruct the original probability landscape.A multitude of variants of metadynamics exist [25], but they are restricted in application to reversible systems, i.e. those with an underlying free energy landscape that can be filled.
The key realization of this paper is the fact that even for irreversible systems of the form (1), the corresponding transition path sampling problem is a gradient flow in path space and thus amenable to metadynamics.In short, while the dynamics might not have an underlying potential landscape, its trajectories do, in the form of the Onsager-Machlup action (3).

2.2.1.
Well-tempered metadynamics.In order to implement metadynamics, consider a reversible stochastic process q t ∈ R n evolving according to implying that the stationary distribution of q t is given by the Gibbs distribution q ∼ ρ(q) ∝ exp(−ε −1 U (q)).In metadynamics, the gradient dynamics (10) are augmented by a biasing potential where the biasing potential is large (and thus repelling) in regions where the system spent a lot of time in the past.That way, exploration is increasingly encouraged over time.
The time-dependent bias potential is assumed to depend on the coordinates q only via a smooth function f ∶ R n → R m with m ≪ n (in typical applications m < 4).The function f should be chosen in such a way that the marginal distribution with respect to f , is multimodal, namely characterized by the presence of at least two distinct maxima, separated by a region in which ρ f (s) is small.The function f is typically called Collective Variable (CV), as it is assumed to be able to recapitulate the most salient features of the probability density.
A typical protocol for building up the biasing potential is well-tempered metadynamics, in which the biasing potential and process co-evolve according to The above intuitively deposits droplets of shape ψ δ (typically Gaussian profiles of width δ) at the current position of the CV, f (q t ), with a weight w > 0. The weight is exponentially decreasing for increasing V (s, t), with a scale given by κ > 0 (see [24,25]).In the longtime limit, the biasing potential V ∞ (s) yields an estimate of ρ f (s) via As a consequence, κ can be seen as an effective sampling temperature or noise amplitude, different from the physical noise amplitude given by ε, and (14) translates between the two.

Pathspace Metadynamics.
Applying now well-tempered metadynamics (13) in collective variables to the Onsager-Machlup stochastic gradient flow (5) yields an effective method for sampling transitions in the presence of multiple possible transition channels.This technique was presented in [34] in the context of equilibrium systems from molecular dynamics.We will derive this technique in our general notation in order to sample transition paths in nonequilibrium systems.
, reducing from continuous trajectories on R n × [0, T ] to m collective variables.For example, one could consider the location of the trajectory at t = T 2, which corresponds to f [φ] = φ(T 2) and thus δf δφ = δ(t − T 2).
Then, for a trajectory φ(t), integrate, in virtual time τ , We term the above pathspace metadynamics, and will use it in the following in several examples of transitions in irreversible processes that exhibit multiple transition channels.

N E
3.1.Two competing transition paths.The simplest case of coexistence of multiple transition paths is 2-dimensional gradient diffusion in a landscape that has two separate valleys connecting the local minima.To this end, consider for Z = (X, Y ) ∈ R 2 the system

F
. The 2-dimensional test model, as defined by equation ( 16).For γ = 0 (left), the drift is a gradient flow in the potential (17).For nonzero γ, such as γ = 0.4 (right), an additional swirl is added.There are two fixed points of the system (white dots), and notably two distinct most likely transition pathways between them (upper and lower).Nonvanishing swirl makes one transition path preferred over the other.In general it is hard to know the relative likelihood of the two transition channels.

F
. Trajectories sampled by transition path sampling (left and center) versus trajectories explored by pathspace metadynamics (right), for the slightly asymmetric 2d example (γ = 0.2).Traditional transition path sampling can only explore one of the two available transition mechanisms, depending on the initial transition trajectory fed to the algorithm.Metadynamics can explore the whole space of possible transitions.and we choose This system is not reversible for γ ≠ 0 and corresponds to a diffusion in a potential with two minima,

F
. Probability of taking the upper or lower transition path as a function of γ.Transitioning from left to right, the upper transition path (i.e.z > 0) is preferred for positive γ: The probability to observe a positive y-value at the center of the time interval, y(T 2) > 0, is increasing, but both upper and lower channel remain present.Shown is the scaled logarithm of the probability, ε log P (z) + C, with arbitrary normalization constant C.
pathway is preferred for a left-to-right transition.For a general system of this type, it is hard to compute all possible transition pathways and their relative likelihood.
Applying pathspace metadynamics to this example corresponds to integrating equation (15) in virtual time for trajectories φ(t, τ ) = (φ x (t, τ ), φ y (t, τ )) connecting the left to the right fixed point.At every instance in virtual time, the solution to the SPDE yields a transition path, and integrating for long times samples from the transition path ensemble.Without employing metadynamics, the SPDE will quickly converge to one of the two transition scenarios, and remain stuck there for exponentially long (virtual) times.Employing metadynamics, instead, allows us to gradually fill up one of the transition channels and eventually explore the other as well.In particular in the scenario γ ≠ 0, where the relative likelihood of the two channels is not obvious, this strategy helps identify their relative weight.
As collective variable, we choose m = 1 and f [φ] = φ y (T 2), i.e. the y-value of the trajectory at the center of the temporal interval.We expect this value to be positive if the upper channel is chosen, and negative for the lower.Note though that f [φ] = 0 not necessarily corresponds to a trajectory crossing the domain at the center.Instead, it might correspond to a trajectory that waits at the initial point until t = T 2 and only afterwards initiates a transition through either channel.
Figure 2 shows a comparison of naive pathspace MCMC (left, center) to pathspace metadynamics (right) for γ = 0.2.While MCMC results in the exploration of only one of the two transition mechanisms (either the upper transition or the lower transition, depending on initialization of the sampler), metadynamics can efficiently explore both these mechanisms, and compare their relative likelihood.As expected for positive γ the lower pathway is less likely for a left-to-right transition than the upper pathway.This picture is confirmed when looking at the probability distribution of the CV, for different values of γ, in figure 3.While for γ = 0, both upper and lower pathway are of equal likelihood, for increasing values of γ the upper pathway becomes more likely.Note that the relative likelihood of the two, for the noise of ε = 10 −2 , makes the upper pathway exponentially more likely even for rather small positive γ. Figure 3 shows the rescaled log-likelihood.

F
. Left: Path density plot obtained from pathspace MCMC sampling for transitions in the oscillating double well potential from x = −1 to x = 1, shown in the x-t-plane.Even for a very large number of samples, the trajectory never leaves the current local minimum determined by the initial condition of the sampling procedure.Right: Path density plot obtained from pathspace metadynamics for the same system.As clearly visible, the process explores all relevant transition trajectories, irrespective of the initial condition of initialization of the sampling procedure.For both plots, ε = 10 −2 and ∆T = 2.
. Left: Limiting biasing potential V ∞ (s) in the space of collective coordinates (s 1 , s 2 ) = (x(T 3), x(2T 3)) for the moving double well experiment.There are three dominant states, A, B, and C, which correspond to late transition, early transition, and dwelling around the saddle, respectively.Right: Sketch of these three transition scenarios.in time is done via second order finite differences.Integration in virtual time is done with Euler-Maruyama with stepsize ∆τ = 10 −3 .

Oscillating doublewell.
Even in one dimension there can be multiple competing reaction channels, for example when the system has explicit time dependence and there are multiple points in time at which a transition is advantageous.As an example, consider the SDE given in equation (1) This is simply a diffusion in a double-well potential that is periodically tilted left or right.In fact, within the time period T the potential tilts exactly twice, allowing for two "optimal" opportunities to jump from the left to the right basin.Therefore, the ensemble of transition paths contains at least two different distinct trajectories, and sampling both of them efficiently is a hard task.Concretely, we want to sample trajectories that transition within the interval [0, T ] from the left basin, X − = −1, to the right basin, X + = 1.We could consider employing classical transition path sampling by integrating φ(t, τ ) via equation ( 8) with φ(0, τ ) = X − and φ(T, τ ) = X + .Then, depending on the initial choice of transition path, we would sample almost certainly only one of the multiple possible transition pathways.This is depicted in figure 4 (left), which shows as streamlines the deterministic, time-dependent drift b(x, t), and in color a weighted histogram of observed samples.If instead we perform transition path sampling with the help of pathspace metadynamics then we see how the whole space of feasible transition paths is explored.In particular, it becomes clear that there are at least three noteworthy classes of transition trajectories: those that transition at the earliest opportunity, those at the latest, and those that spend one oscillation in close proximity to the saddle point (which, despite its instability, turns out to be a very relevant transition trajectory).This is shown in figure 4 (right), where again we show a histogram of transition paths {x(t)} T t=0 .As CVs, we choose i.e. the location of the sampled trajectory at 1 3 and 2 3 of the time interval.We obtain, after running the algorithm ( 15) sufficiently long, a clear separation of the available options, as visible in the converged biasing potential V ∞ (s) in figure 5 (left): Either s 1 remains small, corresponding to a late transition and remaining in the original basin for the first oscillation, labeled as region A. Alternatively, s 2 is large, corresponding to an early transition where the trajectory is in the second basin already during the second oscillation, labeled as region B. Lastly, both s 1 and s 2 are of moderate value, corresponding to an early transition to the saddle, but then dwelling at the saddle for a full oscillation, labeled as region C.The biasing potential gives an estimate of the relative likelihood of these options.These alternatives are sketched in figure 5 (right).

Magnetic field reversal in the presence of advection.
To demonstrate the feasibility of applying pathspace metadynamics to systems with a large number of variables, we will lastly consider a system with infinitely many degrees of freedom, a stochastic partial differential equation.Concretely, consider the Ginzburg-Landau (or equivalently φ 4 -or Allen-Cahn) equation for a field u(x, t) Here, at each point x in space, u(x, t) is traveling in a double-well potential, driving locally towards ±1.The diffusion constant ν couples neighboring locations, and η(x, t) is white in space and time stochastic noise.The system prefers to be in one of the two states u + (x) or u − (x), corresponding to almost all of the domain being at 1 or −1, respectively, except for a small boundary layer, since the boundary is kept at 0 via Dirichlet boundary conditions.The natural question is to ask by what mechanism the model switches from u − to u + .This is a classical example of nucleation theory, with applications in phase separation or magnetic field reversal, and has found particular attention in the small-noise regime and large deviation theory [48][49][50][51][52]. Intuitively, a critical nucleus must be formed to drive parts of the system over the potential barrier.After this, front propagation ensures prevailing of

F
. Snapshots of samples of the four possible transition pathways from the negative configuration u − (x) to the positive configuration u + (x): nucleation from left wall (L), nucleation from right wall (R), nucleation from both walls (B), or nucleation in the center (C).

F
. Left: Limiting biasing potential V ∞ (s) in CV space (s 1 , s 2 ), corresponding to the sin and cos-components at t = T 2. Intuitively, s 1 corresponds to the left-right asymmetry, and s 2 to the center-boundary asymmetry of the configuration in space.Four separate transition channels are identified (R, L, C, B).Right: The same, but in the presence of background flow, γ = 0.03, as in equation (27).Since fronts are propagated predominantly left to right, the nucleation from the left boundary (L) is now clearly preferred over nucleation from the right (R).
the new phase throughout the whole domain.Since the boundaries are kept at 0, nucleating from one of the boundaries is generally easier, so that the system can flip either via a traveling wave from the left or from the right.Depending on the size of the domain L (or equivalently the diffusion constant ν), a nucleus might also form within the domain.
Note that the O(ε)-term, which is the gradient of a divergence in the finite dimensional case, now includes the (functional) trace of the operator δb δφ.Unfortunately, this operator is not necessarily trace class, and the corresponding term might diverge.This is indeed the case for the Ginzburg-Landau-type cubic nonlinearity in (24), and is problem well-known and discussed in the field theory literature [53][54][55].We follow common practice [53] here and drop the diverging term, but remark in passing that this does not rest on any rigorous footing, and in fact a renormalized term of O(ε) might be more appropriate.
Concretely, for our choice of b(φ) = ν∂ 2 x φ − φ 3 + φ and given boundary conditions, this amounts to with boundary conditions As can be seen by the missing ∂ t φ-term, this system is actually reversible, but we will add a term breaking reversibility below.
Integrating (24) corresponds to a Markov chain Monte Carlo for the Ginzburg-Landau equation, and indeed this samples from the transition ensemble from u − to u + .Again, though, starting the system close to one of transition channels (for example nucleating from the left boundary as in figure 6, top left) means that the sampler will effectively never leave its vicinity.Employing pathspace metadynamics, instead, results in exploring all possible transition channels.
As collective variables, we suggest m = 2 with taking the t = T 2-configuration and projecting it onto − sin and cos, respectively.This yields a value for the left-right asymmetry as s 1 , and for the center-boundary asymmetry as s 2 .For example, nucleating from the left boundary means that around t = T 2 most mass is concentrated in the left half-space, yielding a positive value for s 1 (and roughly 0 for s 2 ).
The converged biasing potential in the CV-plane is depicted in figure 7 (left) and reveals the possible transition channels: The landscape clearly has 4 distinct peaks, corresponding to a nucleation from the left boundary (L), from the right boundary (R), from both boundaries at the same time (B), and from the center (C).L and R have identical weight by symmetry, and are clearly preferred over C and B. Figure 6 shows samples of the four individual transition scenarios as individual realizations of the sampler within the four regions, respectively.
The situation becomes more complicated when adding additionally an advection term to the Ginzburg-Landau equation, namely with velocity field v(x) chosen here to be v(x) = γ exp(− ).Since for γ > 0 this velocity field is positive within the domain, and only drops to 0 at the boundaries, one expects it to favor nucleation from the left boundary: The rightward advection helps pushing the phase front across the domain.Nucleation from the right boundary, on the other hand, is suppressed, since the nucleus must expand against the background flow.Nucleation in the center, or from both boundaries, will be distorted by the advection in one direction.Note further that the advection term is not L 2 self-adjoint, and thus the resulting system is no longer reversible.This is intuitively clear since one expects the most likely forward transition from u − to u + to nucleate from the left boundary, but backward from u + to u − to nucleate from the left boundary as well, instead of the time-reversed forward transition.
This intuition is confirmed by the biasing potential, as visible in figure 7 (right).Here, the nucleation from the left boundary (L) is clearly preferred over R for a value of γ = 0.03.Increasing γ further leads to a disappearance of R, while C and B merge with L into the only surviving transition channel.

D
We demonstrated a novel algorithm to sample transition trajectories in complex stochastic systems.The method makes use of the fact that the transition path ensemble can be sampled efficiently via Langevin MCMC and accelerated with metadynamics.This makes sure that the algorithm samples from all possible transition paths even in the presence of multiple competing transition channels.The flexibility comes at the price of choosing the right collective coordinates: Without correct CVs, metadynamics is ineffective or yields misleading results.Choosing suitable CVs necessitates physical intuition about the system at hand, and it is not always easy to find a good set of CVs.The problem of finding CVs for generic sampling problems is classical, and essential not only for our algorithm here, but for all rare event algorithms.It can be shown that the optimal collective variable is the committor function, which is usually not available for any problem of practical relevance.As mentioned in the introduction, good heuristics for CVs can come from physical intuition, and machine learning and data science techniques might offer more systematic insight [28][29][30][31].
The presented approach requires storing in memory the whole pathway, posing an intrinsic limitation to the size of the systems which can be studied with it.However, because of the increasing speed of modern computers, we demonstrate that even stochastic partial differential equations are feasible to be treated without any special effort or hardware.Since the underlying pathspace Langevin equation is the same as the gradient decent equation for action minimization (for example for instanton computation in the Freidlin-Wentzell sense), similar optimizations can be employed to sample transition paths in even higher dimensional PDEs, as demonstrated for example for the 3D Navier-Stokes equation in [57].
Metadynamics has been suggested for use in pathspace before [12,13,34].Generally in those works, the O(ε)-term of the Onsager-Machlup action (3) is omitted, even though it has been pointed out in other contexts [46,58] that the term is crucial for the correct behavior of transition paths at non-zero temperature.Here, we chose to discretize the continuous limiting SPDE, which must include this term, but conceivably one might alternatively discretize the original stochastic process first and then derive the bridge sampler fully discretely.This might allow to sidestep issues of the renormalizability of the bridge sampling SPDE, in particular for higher dimensional continuous problems.
All previous applications of metadynamics to pathspace restrict the procedure to the setup of equilibrium systems and gradient flows, where b(X) = −∇V (X).As shown in our examples above, the method presented in this work is generally applicable to a wide range of situations, including irreversible dynamics without underlying potential landscape, explicit time dependence of the drift (or noise) terms, etc. Due to the efficiency of metadynamics, it works even for high-or infinite-dimensional systems such as SPDEs.Additionally, our method is readily generalized to situations such as computing expectations over transition paths, sampling time-periodic paths, sampling trajectories that visit a sequence of points, sampling trajectories with a time-integrated constraint, and many other scenarios encountered in different problem domains.
), but with an additional non-gradient force l(z) of strength γ added.The potential, and the flowlines of the full drift, are depicted in figure1.Notably, there are two distinct most likely transition pathways between the left and the right fixed point, through the upper or the one lower channel.If γ = 0, both of these are equally likely due to symmetry.If γ > 0, as visible in figure1(right), the upper