Improvements to the use of the Trajectory-Adaptive Multilevel Sampling algorithm for the study of rare events

The Trajectory-Adaptive Multilevel Sampling (TAMS) is a promising method to determine probabilities of noise-induced transition in multi-stable high-dimensional dynamical systems. In this paper, we focus on two improvements of the current algorithm related to (i) the choice of the target set and (ii) the formulation of the score function. In particular, we use confidence ellipsoids determined from linearised dynamics in the choice of the target set. Furthermore, we define a score function based on empirical transition paths computed at relatively high noise levels. The suggested new TAMS method is applied to two typical problems illustrating the benefits of the modifications.

1 Introduction

Systems from various areas of physics exhibit multiple stable states. In such multi-stable systems, transitions between states can occur as a result of small-scale processes, usually referred to as noise-induced transitions (Ashwin et al., 2012). Typical elements in the Earth's system which show multistability include the Greenland Ice Sheet (Ridley et al., 2010; Robinson et al., 2012), the Amazon rainforest (Higgins and Scheiter, 2012; Lasslop et al., 2016) and the Atlantic Meridional Overturning Circulation (AMOC). In particular, the latter can undergo transitions to a collapsed state due to fluctuations in the surface freshwater forcing (Castellana et al., 2019).

A central issue in models of these multi-stable systems is the computation of transition probabilities between different states. If we exclude very special classes of systems, analytical results are generally not available. The Eyring–Kramers formula (Eyring, 1935; Kramers, 1940), which allows the computation of transition rates for reversible processes in the zero noise limit, has been recently generalised to non-gradient systems (Bouchet and Reygner, 2016). However, this method involves the calculation of quasi-potentials, which are generally hard to compute from their variational characterisation. From the numerical point of view, the naive method would be to follow a Monte Carlo approach by performing simulations of large ensembles of trajectories and calculate transition probabilities by counting the number of trajectories which actually undergo a transition. However, if the occurrence of a transition is a rare event, such computations are not feasible. For instance, to sample an event of probability $p\sim {\mathrm{10}}^{-\mathrm{8}}$, one would need to compute at least $N>N_{\text{min}}={\mathrm{10}}^{\mathrm{8}}$ trajectories ($N_{\text{min}}\sim \mathrm{1}/p$), which is currently impossible to achieve for large-dimensional dynamical systems, where time integrations are expensive.

In order to sample tails of distributions more effectively, various methods have been developed, generally referred to as rare-event algorithms. One of the promising methods to compute transition probabilities is the Trajectory-Adaptive Multilevel Sampling (TAMS) method (Lestang et al., 2018). Its underlying idea is to perform a selection/mutation process that discards trajectories going away from a certain target set and splits/branches from those that get closer to this set. A very similar algorithm (Adaptive Multilevel Splitting or AMS), based on the same approach, has been used in the study of transitions in Jupiter's turbulent dynamics (Bouchet et al., 2019) and in molecular dynamics to compute the expected dissociation time between a protein and its ligand (Teo et al., 2016). In these studies, AMS proved to be a powerful tool that reduced computational costs by several orders of magnitude. Indeed, the required minimum number of computed trajectories scales like $N_{\text{min}}\sim \mathrm{1}/\log p$ (Cérou et al., 2016), which is exponentially better than that for the Monte Carlo estimation. The aforementioned selection and mutation process of discarding and branching trajectories is carried out according to a score function, which allows us to rank trajectories. Rolland and Simonnet (2015) have shown that the choice of the score function plays an important role in the performance of the algorithm, even for systems with only 2 degrees of freedom. When using non-optimal score functions, especially near phase transitions, the variance of the estimated probability can peak and the convergence of the algorithm can be slow.

The aim of this work is to propose improvements to the use of the TAMS algorithm to be able to compute transitions in multi-stable systems more efficiently. The first type of improvement is the choice of the target set, which is often determined from rather arbitrary thresholds. This choice also raises more broadly the question of a precise definition of what we consider a noise-induced transition between two (stable) states. In the second type of improvement, we propose a more systematic method of defining a score function, i.e. based on empirical transition paths. The modified TAMS method is first applied to an idealised gradient system and then to a system representing a box model of the AMOC (Castellana et al., 2019).

In Sect. 2, we describe the methods developed to improve the TAMS algorithm. In Sect. 3 we show how to incorporate these techniques into the definition of the score function and present the results for idealised dynamical systems and the AMOC model. A discussion follows in Sect. 4, assessing the strengths and the limitations of the new TAMS method.

2 Methods

2.1 Transition probabilities using TAMS

We consider finite-dimensional dynamical systems described by stochastic differential equations (SDEs) of the following form:

$$\begin{array}{}\text{(1)}& \mathrm{d}{\mathit{X}}_t=\mathit{F}\left({\mathit{X}}_t\right)\mathrm{d}t+\mathbf{G}\mathrm{d}{\mathit{W}}_t,\end{array}$$

where X_t∈ℝⁿ and $\mathit{F}:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is the drift field. The noise term W_t∈ℝ^m consists of m independent Wiener processes, with the matrix $\mathbf{G}:{\mathbb{R}}^n\times {\mathbb{R}}^m$ being the noise matrix.

A prominent example of such a system is a model of a free particle moving in a two-dimensional double-well potential (with $n=m=\mathrm{2}$). The drift term in the time-evolution equation for the variables x and y is in this case

$$\begin{array}{}\text{(2)}& \mathit{F}(x,y)\equiv -\mathbf{\nabla }V(x,y)=(x-x^{\mathrm{3}},-y),\end{array}$$

where V(x,y) represents the potential (Fig. 1a). In the deterministic case (i.e. G=0), the stable steady states of the system are ${\mathit{X}}_A=(-\mathrm{1},\mathrm{0})$ and ${\mathit{X}}_B=(+\mathrm{1},\mathrm{0})$, while the unstable steady state is ${\mathit{X}}_C=(\mathrm{0},\mathrm{0})$. Without fluctuations, a particle starting in X_A stays in X_A. However, because of the presence of unresolved processes (such as thermal fluctuations), modelled by the noise term in the SDE, the particle can move away from X_A and in some cases make the transition to the state X_B. An example of such a transition is shown in Fig. 1b.

Figure 1 (a) Iso-potential contours corresponding to the double-well potential (Eq. 2). The two stable steady states ${\mathit{X}}_A=(-\mathrm{1},\mathrm{0})$ and ${\mathit{X}}_B=(+\mathrm{1},\mathrm{0})$ are marked with a filled circle and a cross, respectively. The saddle ${\mathit{X}}_C=(\mathrm{0},\mathrm{0})$ is indicated by a circle. (b) Example of a noise-induced transition from X_A to X_B for a particle moving in the double-well potential in (a). The noise matrix in the general SDE (1) is chosen as G=σI, with σ=0.32. The red circle denotes an arbitrarily defined target set B.

The transition probability that a trajectory starting in X_A reaches a neighbourhood B around X_B before time T is indicated by $P({\mathit{\tau }}_B<T|{\mathit{X}}_{\mathrm{0}}={\mathit{X}}_A)$. Here, τ_B denotes the stopping time associated with reaching B. The TAMS algorithm is based on a selection and mutation process of discarding and branching trajectories, which are ranked according to a score function $\mathit{\varphi }:{\mathbb{R}}^n\times \mathbb{R}\to \mathbb{R}$. For a given state $(\mathit{X},t)\in {\mathbb{R}}^n\times \mathbb{R}$, ϕ(X,t) is supposed to measure how likely it is to start in (X,t) and reach B before time T. As a consequence, if the choice of score function is successful, the probability that a trajectory reaches the target set B keeps increasing at each iteration of the algorithm, which is why this method is more efficient than brute-force techniques. A visual representation of the algorithm is given in Fig. 2 and a step-by-step description is provided in Appendix A (more details can be found in Lestang et al., 2018).

Figure 2 Illustration of the TAMS algorithm. First, simulate N trajectories starting in X_A (N=3 blue trajectories in the figure). The trajectories are ranked according to their score ϕ_i, which is the maximum value of the score function ϕ along the trajectory. Then, at each iteration, the trajectory with the lowest score (1 in the figure, with score ϕ₁) is discarded. It is replaced by picking a trajectory uniformly at random from the other ones (2 in the figure), then computing the earliest position at which it reaches a higher score than the discarded trajectory (orange dot) and finally using this position as the branching point for the new trajectory (1^′). This is repeated until all the trajectories reach B or the number of iterations reaches a predefined limit.

2.2 Score functions

Consider a general SDE (1) and two stable states X_A and X_B of the corresponding deterministic system. The TAMS algorithm (see Appendix A) needs a score function to be defined, which allows us to rank trajectories and select the ones to discard. The optimal score function ϕ_com(X,t), i.e. the score function that minimises the variance of the probability estimator, is called the static committor. Its generic expression is given by the following conditional probability (Lestang et al., 2018):

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}_{\text{com}}(\mathit{X},t)=P({\mathit{\tau }}_B<T|\mathit{X},t),\end{array}$$

where T is again the fixed duration of the trajectories and τ_B is the stopping time associated with reaching the target set B. In other words, ϕ_com(X,t) is the probability that a trajectory starting in X at time t reaches B before time T. This expression is quite natural because ideally, the score of (Y,s) should be higher than the score of (X,t), i.e. $\mathit{\varphi }(\mathit{X},t)\le \mathit{\varphi }(\mathit{Y},s)$, if and only if $P({\mathit{\tau }}_B<T|\mathit{X},t)\le P({\mathit{\tau }}_B<T|\mathit{Y},s)$. This condition is clearly satisfied by ϕ_com (and in fact any increasing function of ϕ_com). However, the expression given in Eq. (3) is generally unusable because it is the very quantity that we want to compute. For instance, ϕ_com(X_A,0) is precisely the transition probability that TAMS estimates. As a conditional probability of the form $P(\mathit{Y},s|\mathit{X},t)$, the committor ϕ_com(X,t) satisfies the backward Fokker–Planck equation (Lestang et al., 2018):

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+\mathit{F}\cdot \mathbf{\nabla }\mathit{\varphi }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i,j}{\mathbf{G}}_{ij}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\varphi }}{\partial {\mathit{X}}_i\partial {\mathit{X}}_j}}=\mathrm{0},{\displaystyle }& {\displaystyle }\text{with boundary conditions}\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}\forall \mathit{X}\in \partial B,\forall t\in [\mathrm{0},T],\mathit{\varphi }(\mathit{X},t)=\mathrm{1},\\ \forall \mathit{X}\in {\mathbb{R}}^n\setminus B,\mathit{\varphi }(\mathit{X},T)=\mathrm{0}.\end{array}\right.\end{array}$$

However, solving the backward Fokker–Planck equation in systems with many degrees of freedom is computationally unfeasible. Moreover, even if the committor is available on a grid used for the discretisation of Eq. (4), using interpolation to evaluate it during a TAMS loop can also have a prohibitive computational cost.

Bouchet et al. (2019) made use of a score function based on the distances of the state X from the starting state X_A and the destination equilibrium X_B, respectively. It is defined as

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\text{dist}}\left(\mathit{X}\right)\equiv \left\{\begin{array}{ll}{d}_A/\mathrm{2}{d}_B,& \text{if}{d}_A<d_B,\\ \mathrm{1}-d_B/\mathrm{2}{d}_A,& \text{otherwise,}\end{array}\right.{\displaystyle }& {\displaystyle }\text{with}{d}_A\equiv \parallel \mathit{X}-{\mathit{X}}_A\parallel ,d_B\equiv \parallel \mathit{X}-{\mathit{X}}_B\parallel .\end{array}$$

Alternatively, a Gaussian-shaped score function ϕ_gauss was proposed in Baars (2019). It is defined as

$$\begin{array}{ll}\text{(6)}& {\displaystyle }{\mathit{\varphi }}_{\text{gauss}}\left(\mathit{X}\right)& {\displaystyle }\equiv \mathit{\eta }-\mathit{\eta }{e}^{-\mathit{\beta }\parallel \mathit{X}-{\mathit{X}}_A{\parallel }^{\mathrm{2}}/\parallel {\mathit{X}}_C-{\mathit{X}}_A{\parallel }^{\mathrm{2}}}{\displaystyle }& {\displaystyle }+(\mathrm{1}-\mathit{\eta })e^{-\mathit{\beta }\parallel \mathit{X}-{\mathit{X}}_B{\parallel }^{\mathrm{2}}/\parallel {\mathit{X}}_C-{\mathit{X}}_B{\parallel }^{\mathrm{2}}},\end{array}$$

with $\mathit{\eta }\equiv \parallel {\mathit{X}}_C-{\mathit{X}}_A\parallel /\parallel {\mathit{X}}_A-{\mathit{X}}_B\parallel$. Here β∈ℝ is a parameter controlling the decay and X_C is the saddle state of the system. Due to the general expressions of ϕ_gauss and ϕ_dist, these score functions can be used in systems of any dimension.

2.3 Definition of the target set

Once the score function has been chosen, a threshold needs to be defined for the TAMS algorithm to converge, so that the occurrence of a transition can be detected. In other words, we do not expect each trajectory that undergoes a transition to reach exactly the destination equilibrium X_B, but rather a neighbourhood of it (B). The target set B can then be defined according to a level set ϕ^target of the score function ϕ:

$$\begin{array}{}\text{(7)}& B=\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \mathit{\varphi }(\mathit{X})>{\mathit{\varphi }}^{\text{target}}\mathit{\}}.\end{array}$$

However, different score functions ϕ and different level sets ϕ^target correspond to different target sets B, which can differ in volume and in shape. Moreover, often the level set ϕ^target is defined somewhat arbitrarily. For example, by using ϕ_gauss with target score ${\mathit{\varphi }}^{\text{target}}=\mathrm{0.85},\mathrm{0.9}$ or 0.95, we found that the average of the transition probability estimator for a two-dimensional double-well potential system (Eq. 2) can vary by up to 30 %. This may not be a concern if one only cares about the order of magnitude of the transition probability, but it will be problematic if quantitative comparisons are needed. Moreover, a poor choice of target set can lead to inaccuracies when a trajectory has a score greater than ϕ^target without actually making the transition in all the degrees of freedom. Thus, there is a need to define a canonical choice of a target set B, which needs to fundamentally address what is considered to be a noise-induced transition.

For this purpose, we use the concept of a confidence ellipsoid. This is an ellipsoidal neighbourhood around a stable equilibrium state, inside which a trajectory subject to the locally linearised dynamics stays, within a certain confidence level (Cowan, 1998).

Figure 3 Confidence ellipsoids for the two-dimensional double-well potential system (Eq. 2), with confidence level $\mathrm{1}-\mathit{\alpha }=\mathrm{0.95}$ and noise matrix G=0.25I, where I is the identity matrix in ℝ². Two trajectories of duration T=200 are shown: one (blue) is initialised at the initial state and one (orange) at the target state. 93 % of the points composing a trajectory are inside their corresponding ellipsoid. This is lower than the prescribed confidence level because away from the equilibrium, the first-order dynamics from which the confidence ellipsoids are derived does not hold.

Consider a general SDE system given by Eq. (1). Because the drift vanishes at an equilibrium state F(X_B)=0, its first-order approximation around the equilibrium state X_B via Taylor expansion is

$$\begin{array}{}\text{(8)}& \mathit{F}\left(\mathit{X}\right)=\mathbf{A}\left({\mathit{X}}_B\right)(\mathit{X}-{\mathit{X}}_B)+\mathcal{O}(\parallel \mathit{X}-{\mathit{X}}_B{\parallel }^{\mathrm{2}}),\end{array}$$

where $\mathbf{A}\left({\mathit{X}}_B\right)=\mathrm{\nabla }F\left({\mathit{X}}_B\right)\in {\mathbb{R}}^{n\times n}$ is the Jacobian matrix of F at X_B. Using a translation $\stackrel{\mathrm{̃}}{\mathit{X}}=\mathit{X}-{\mathit{X}}_B$, the first-order approximation of the SDE is then

$$\begin{array}{}\text{(9)}& \mathrm{d}\stackrel{\mathrm{̃}}{{\mathit{X}}_t}=\mathbf{A}\left({\mathit{X}}_B\right)\stackrel{\mathrm{̃}}{{\mathit{X}}_t}\mathrm{d}t+\mathbf{G}\mathrm{d}{\mathbf{W}}_t.\end{array}$$

Because the drift term has been linearised, this is the equation for an n-dimensional Ornstein–Uhlenbeck process. The stationary probability density function (PDF) f of the approximating process is Gaussian (Cowan, 1998) and given by

$$\begin{array}{}\text{(10)}& f\left(\mathit{X}\right)={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{2}\mathit{\pi }{)}^{\frac{n}{\mathrm{2}}}\right|{\mathbf{C}}_B{|}^{\frac{\mathrm{1}}{\mathrm{2}}}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}\parallel \mathit{X}-{\mathit{X}}_B{\parallel }_{{\mathbf{C}}_B^{-\mathrm{1}}}^{\mathrm{2}}},\end{array}$$

where ${\mathbf{C}}_B\in {\mathbb{R}}^{n\times n}$ is the covariance matrix of the system calculated in X_B and $\parallel .{\parallel }_{{\mathbf{C}}_B^{-\mathrm{1}}}$ the norm induced by its inverse ${\mathbf{C}}_B^{-\mathrm{1}}$, defined by $\parallel \mathit{X}{\parallel }_{{\mathbf{C}}_B^{-\mathrm{1}}}^{\mathrm{2}}\equiv {\mathit{X}}^{\top }{\mathbf{C}}_B^{-\mathrm{1}}\mathit{X}$. The covariance matrix C_B can be thought of heuristically as the matrix containing the correlations 𝔼(x_ix_j) (with $\mathit{X}=(x_{\mathrm{1}},\mathrm{\dots },x_n)$), which generalises the notion of variance in n dimensions. C_B is obtained by solving the Lyapunov equation (see Kuehn, 2012, for the full derivation):

$$\begin{array}{}\text{(11)}& \mathbf{A}\left({\mathit{X}}_B\right){\mathbf{C}}_B+{\mathbf{C}}_B\mathbf{A}({\mathit{X}}_B{)}^T+{\mathbf{GG}}^T=\mathrm{0}.\end{array}$$

We then define the confidence ellipsoid, which has ${\mathbf{C}}_B^{-\mathrm{1}}$ as a shape matrix, as follows:

$$\begin{array}{ll}{\displaystyle }\mathcal{E}& {\displaystyle }=\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \parallel \mathit{X}-{\mathit{X}}_B{\parallel }_{{\mathbf{C}}_B^{-\mathrm{1}}}^{\mathrm{2}}\\ \text{(12)}& {\displaystyle }& {\displaystyle }\equiv (\mathit{X}-{\mathit{X}}_B{)}^T{\mathbf{C}}_B^{-\mathrm{1}}(\mathit{X}-{\mathit{X}}_B)<Q_{\mathit{\alpha }}\mathit{\}},\end{array}$$

where Q_α is the quantile of confidence level 1−α of the n-dimensional χ² distribution (Cowan, 1998), usually $\mathrm{1}-\mathit{\alpha }=\mathrm{0.95}$. The directions of symmetry of the ellipsoids are given by the eigenvectors of the covariance matrix C_B and the radii are given by the corresponding eigenvalues and the confidence level 1−α. Intuitively, the greater the eigenvalue, the more a trajectory fluctuates in the given direction.

The (1−α)-covariance ellipsoid represents the n-dimensional volume where a trajectory is confined with confidence level 1−α, provided its dynamics is well approximated by the first-order expansion at the equilibrium point. An illustration of the 0.95-confidence ellipsoid for the double-well potential is shown in Fig. 3. The confidence ellipsoid ℰ constitutes a way to meaningfully define the target set B with minimal arbitrary parameters. As shown in Fig. 3, when initialising a trajectory at X_A or X_B in the two-dimensional double-well system, it stays inside the correspondent ellipsoid with a certain confidence $\mathrm{1}-\mathit{\alpha }=\mathrm{0.95}$. In the next section, we show how to incorporate this choice of target set into the score function ϕ.

2.4 Estimating the typical transition path using histograms

The second line of improvement of the score function concerns the estimation of typical transition paths of the dynamical system. In the zero noise limit, the Freidlin–Wentzell theory of large deviations predicts that transition paths will cluster around the most probable transition path, called the instanton (Freidlin and Wentzell, 1984). On the other hand, in the finite noise regime, transition paths may deviate from the instanton. Moreover, instantons may be computationally inaccessible for systems with many degrees of freedom. Therefore, it can prove more relevant to estimate empirically the typical transition path that the system follows at a given finite noise level, which is the approach we follow here.

The idea is to first accumulate transition paths at a noise level where transitions are frequent enough (typically $p>{\mathrm{10}}^{-\mathrm{3}}$) so that any sampling method (direct Monte Carlo or TAMS with naive score functions) can be used. Then, we compute the spatial histogram of the transition paths over a discretised phase space using n-dimensional boxes. This provides the spatial distribution of the transition paths, which is concentrated around a typical transition path, reminiscent of an instanton phenomenology, which was also observed in more complex systems (Bouchet et al., 2019). From the spatial histogram, we extract a typical transition path.

The main steps of the path-finding algorithm are listed below:

• i.

  the trajectory of the typical transition path starts in the box of the histogram containing the initial state X_A;

• ii.

  the next box in this trajectory corresponds to the neighbour which has the highest nonzero histogram value but which has not already been visited by the typical transition path;

• iii.

  the algorithm stops if it reaches the box containing the target state X_B.

In addition, the full typical path estimation algorithm uses a self-correcting method to avoid dead ends when there are no valid neighbours to be the next point in the trajectory. The spirit of the path-finding algorithm is similar to the depth-first search algorithm (Cormen et al., 2009). We found that, as long as the histogram is not fragmented, i.e. there is a sufficient number of accumulated trajectories or large enough histogram bins, the algorithm converges. Possible artefacts created by this estimation include spiralling near the initial equilibrium (because of the concentric shell structure of the local probability density function, Eq. 10) and zigzagging at the histogram box size. They can both be addressed by a clean-up algorithm: starting from the first box, at each box X_j, if the trajectory goes back to one of its neighbours at a later time, with X_j+k being the latest neighbour visit, we erase the points ${\mathit{X}}_{j+\mathrm{1}},\mathrm{\dots },{\mathit{X}}_{j+k}$ from the trajectory.

In a general system, transitions induced at high noise do not necessarily follow similar transition paths at lower noise levels. However, high-noise estimates are robust as long as there is no drastic change in the behaviour of the system at an intermediate noise level, which would signal a physical phase transition in the system (Rolland and Simonnet, 2015). Moreover, we can check a posteriori if such a transition occurs by starting at high noise, gradually decrease the noise level and apply the empirical estimation of the typical path or simply compare histograms each time the noise level is decreased.

3 Results

In this section, we apply both modifications to TAMS (ellipsoids in the score function and typical path estimation) to different problems.

3.1 Incorporating ellipsoids in the score function

First, we apply the modified TAMS method to the two-dimensional double-well system defined by Fig. 1. As outlined in Sect. 2.3, given any score function ϕ, we want the target set $B=\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \mathit{\varphi }(\mathit{X})>{\mathit{\varphi }}^{\text{target}}\mathit{\}}$ to coincide with the 0.95-confidence ellipsoid ℰ associated with the equilibrium state X_B. However, because the contour levels of ϕ, in general, do not have the shape of an ellipsoid, there is no suitable choice of ϕ^target such that B coincides with ℰ. Here, we propose a general method to modify any score function ϕ so that we are able to choose the target set B to be exactly the 0.95-confidence ellipsoid ℰ of X_B.

Figure 4 (a) Contour levels of the modified score function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{gauss}}$, for the two-dimensional double-well system represented in Fig. 1, according to the procedure described in Eq. (14). The level set ϕ^target (dotted blue line) of the former score function ϕ_gauss is tangent to the ellipsoid. The level set ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}$ (dotted red) of the modified score function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{gauss}}$ coincides with the ellipsoid. (b) Same plot for the score function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{dist}}$.

We first compute the level ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}$, defined as the minimum of the score function ϕ on the confidence ellipsoid ℰ:

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\equiv \underset{\mathit{X}\in \mathcal{E}}{min}\mathit{\varphi }\left(\mathit{X}\right),\end{array}$$

such that the set $\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \mathit{\varphi }(\mathit{X})>{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\mathit{\}}$ contains the ellipsoid ℰ and is tangent to ℰ. This can be done numerically by generating a mesh of points around X_B, then selecting the points inside ℰ by comparing their norm $\parallel \mathit{X}-{\mathit{X}}_B{\parallel }_{C_B^{-\mathrm{1}}}^{\mathrm{2}}$ with the quantile Q_α and finally computing the minimum ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}$ of ϕ on these points. Then, define the modified score function $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$ in the following way:

$$\begin{array}{}\text{(14)}& \stackrel{\mathrm{̃}}{\mathit{\varphi }}\left(\mathit{X}\right)\equiv \left\{\begin{array}{ll}\mathrm{1}& \mathrm{\text{if}}\mathit{X}\in \mathcal{E},\text{i.e.}\parallel \mathit{X}-{\mathit{X}}_B{\parallel }_{C_B^{-\mathrm{1}}}^{\mathrm{2}}<Q_{\mathit{\alpha }},\\ {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}& \text{if}\mathit{\varphi }\left(\mathit{X}\right)>{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\text{and}\mathit{X}\notin \mathcal{E},\\ \mathit{\varphi }\left(\mathit{X}\right)& \text{otherwise.}\end{array}\right.\end{array}$$

The target set $B=\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \stackrel{\mathrm{̃}}{\mathit{\varphi }}(\mathit{X})>{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\mathit{\}}$ for the modified score function $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$ turns out to coincide with ℰ. We apply this procedure on both the score functions ϕ_gauss and ϕ_dist, and the results for the improved score functions ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{dist}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{gauss}}$ are shown in Fig. 4.

3.2 Designing a score function based on a typical transition path

In order to show how to design a score function based on a typical transition path, we consider a two-dimensional system slightly less trivial than the double-well system, i.e. a two-dimensional system with the following potential,

$$\begin{array}{ll}{\displaystyle }V(x,y)& {\displaystyle }=\underset{\text{global confinement}}{\underbrace{\mathrm{0.1}{x}^{\mathrm{2}}+\mathrm{0.05}{y}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+\underset{\text{potential barrier at}y<\mathrm{15},x=\mathrm{0}}{\underbrace{\mathrm{30}{e}^{-(x/\mathrm{2}{)}^{\mathrm{2}}}\left(\mathrm{1}+\tanh (\mathrm{15}-y)\right)}}\\ {\displaystyle }& {\displaystyle }-\underset{\text{left potential well}}{\underbrace{\mathrm{10}{e}^{\left(\right(x+\mathrm{6})/\mathrm{2}{)}^{\mathrm{2}}-(y/\mathrm{2}{)}^{\mathrm{2}}}}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }-\underset{\text{right potential well}}{\underbrace{\mathrm{10}{e}^{-\left(\right(x-\mathrm{6})/\mathrm{2}{)}^{\mathrm{2}}-(y/\mathrm{2}{)}^{\mathrm{2}}}}},\end{array}$$

depicted in Fig. 5a. It consists of two energy minima at ${\mathit{X}}_A\approx (-\mathrm{5.77},\mathrm{0})$ and ${\mathit{X}}_B\approx (\mathrm{5.77},\mathrm{0})$ and a potential barrier spanning y<15 and at x=0. The dynamics is then given by the SDE (1), with drift $\mathit{F}\left({\mathit{X}}_t\right)=-\mathbf{\nabla }V\left({\mathit{X}}_t\right)$.

The dynamics of the system is quite interesting, as it exhibits two distinct regimes for transition paths, depending on the noise level σ (assuming G=σI). At high noise (σ²≫ΔV, where the potential barrier height is $\mathrm{\Delta }V=V(\mathrm{0},\mathrm{0})-V_A$), trajectories are likely to cross the potential barrier. At low noise (σ²≪ΔV), trajectories are not likely to cross the barrier, and the trajectories which undergo the transition instead go through the upper channel at y≈15. Typical examples of such trajectories are shown in Fig. 5b. Rolland and Simonnet (2015) investigated the convergence properties of AMS using a triple-well potential system, which also exhibits two regimes of preferred transition paths depending on the noise. They found a strong dependency of the statistics of the rare-event algorithm (e.g. the number of iterations) and the duration of reactive trajectories on the choice of score function. We expect the same behaviour when applying TAMS to the system with potential (Eq. 15) and that this system will reveal differences in performance between various score functions.

Figure 5 (a) Potential landscape of the two-dimensional gradient system defined by Eq. (15). Two energy minima are located at ${\mathit{X}}_A\approx (-\mathrm{5.77},\mathrm{0})$ and ${\mathit{X}}_B\approx (\mathrm{5.77},\mathrm{0})$. They are separated by a potential barrier spanning y<15 and at x=0. (b) Typical transition paths from X_A to X_B for two noise levels σ. At high noise, trajectories can cross the potential barrier (σ=10, blue). At low noise, trajectories go through the upper channel (σ=3, red).

Figure 6a shows a histogram computed with 300 transition paths for the system with the potential given by Eq. (15), with noise level set at σ=3. We need enough trajectories so that the histogram gives an accurate representation of all possible transition paths and that fluctuations along the typical path can be effectively averaged out by the path estimation algorithm. One way to verify that there are enough trajectories is to check that the histogram near the equilibria is similar to the Gaussian stationary distribution that arises from locally linearised Ornstein–Uhlenbeck dynamics. For example, starting from an equilibrium point and picking any direction, then if for the first few histogram cells the histogram value is not decreasing (as would be the case for a Gaussian distribution), this would mean that too few trajectories were used. The typical transition path was estimated using the algorithm sketched in Sect. 2.4. As already mentioned, some artefacts created by the estimation (such as spiraling or zigzagging) can be corrected using a clean-up algorithm. The result is shown in Fig. 6b: the empirical estimation of the typical path resembles the instanton around which the transition paths are clustered at lower noise (Fig. 6b, blue path). The instanton was calculated by implementing the geometric minimum action method (Heymann and Vanden-Eijnden, 2008). Note that if, instead, we accumulated trajectories at high noise σ>10, we would obtain trajectories going from X_A to X_B in a straight line, which is the typical path at high noise similar to Rolland and Simonnet (2015). This typical transition path is then substantially different from the instanton, which goes through the upper channel. Thus, our method can be advantageous in multistable systems where the typical transition path depends on the noise level. We can start at high noise and reapply the empirical estimation of the typical path each time the noise level is decreased.

Figure 6 (a) Histogram of N=300 transition paths at noise level σ=3, using the modified score function ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{gauss}}$ with β=1.5, defined by Eq. (6) and implemented with Eq. (14). The corresponding transition probability ($p=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$) is high enough so that direct Monte Carlo sampling could have been used to produce a similar histogram. The bin resolution (Δx=1.4, Δy=1.75) is coarse for illustration purposes. The histogram is used as input in the path-finding algorithm which produces the transition path in red. Grid-scale spiralling occurs near the initial state X_A because of the concentric shell structure of the local probability density function given by Eq. (10). (b) Same histogram as panel (a). The estimated typical transition path (red) has been cleaned up from its grid-scale spiraling and zigzags with the clean-up algorithm and has then been smoothed. It strongly resembles the instanton (blue) which was computed by implementing the geometric action minimum method (Heymann and Vanden-Eijnden, 2008).

Given a typical transition path 𝒞, we present the design of a score function ϕ_𝒞 which encourages trajectories to follow the transition path 𝒞 such that it gives a reasonable approximation of the static committor. Let us consider a trajectory 𝒞(s) in ℝⁿ parameterised by arclength $s\in [\mathrm{0},\mathrm{1}]$. Then we define the score function ϕ_𝒞, called the path-based score function, such that it grows from 0 to 1 along the trajectory from X_A to X_B and decays exponentially along the direction transverse to the trajectory:

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_{\mathcal{C}}:{\mathbb{R}}^n& {\displaystyle }\to [\mathrm{0},\mathrm{1}]\\ \text{(16)}& {\displaystyle }\mathit{X}& {\displaystyle }\mapsto s(\mathit{X},\mathcal{C})\exp \left(-{\displaystyle \frac{d(\mathit{X},\mathcal{C}{)}^{\mathrm{2}}}{d_{\mathrm{0}}^{\mathrm{2}}}}\right),\end{array}$$

where $d(\mathit{X},\mathcal{C})={inf}_{s\in [\mathrm{0},\mathrm{1}]}\Vert \mathit{X}-\mathcal{C}\left(s\right){\Vert }_d$ is the distance between X and the trajectory 𝒞(s), s(X,𝒞) is the curvilinear coordinate of the position on the trajectory that satisfies the infimum in the definition of d(X,𝒞) and d₀ is the characteristic decay length (free parameter). The score function ϕ_𝒞 is shown in Fig. 7 for the estimated transition path 𝒞 shown in Fig. 6 and two values of decay length d₀=20 and d₀=200.

Figure 7 (a) Contour levels of the score function ϕ_𝒞 defined in Eq. (16) associated with the estimated transition path 𝒞 shown in Fig. 6 and with decay length d₀=20. (b) Same figure for d₀=200.

The score function increases from 0 to 1 along the trajectory. Thus, it encodes the preferred direction that the system has to follow. This contrasts with the generic score functions ϕ_gauss and ϕ_dist, which are symmetrical in y: they do not contain the information that the system has to increase in y in order to make the transition to X_B. Note that the method we developed here can be applied, in principle, to systems of any dimension. As shown in Fig. 7b, the score function ϕ_𝒞 is discontinuous because the trajectory has positive curvature. The discontinuity is located near the axis x=0. Indeed, when crossing the axis x=0, the closest point on 𝒞 changes and s(X,𝒞) is discontinuous.

Next, we applied TAMS with the path-based score function ϕ_𝒞 to the two-dimensional system in Eq. (15). We compare its performance with the previously defined score functions ϕ_dist and ϕ_gauss. In fact, we use the associated modified score functions, such that the target set B matches the 0.95-confidence ellipsoid (we drop the tildes for readability). We use the following parameters: duration of a trajectory T=20, time step dt=0.01, and N=100 trajectories in the ensemble, which we found was enough to ensure that the interquartile range of the independent realisations spans less than 1 order of magnitude around the mean for most noise levels.

We show in Fig. 8a the transition probability estimates using the score functions ϕ_dist, ϕ_gauss (with β=1.5) and ϕ_𝒞 (with decay parameters $d_{\mathrm{0}}=\mathrm{2},\mathrm{20},\mathrm{200}$) averaged over 10 instances of the algorithm. The probability estimates are in good agreement with each other and with a Monte Carlo estimation for σ>2.5. The score function ϕ_𝒞 is robust with respect to the choice of the decay length d₀.

The performances of the numerical methods are next measured using the work-normalised relative error ϵ, which combines the variance of the algorithm and its computational cost (Glynn et al., 2009):

$$\begin{array}{}\text{(17)}& \mathit{ϵ}\equiv {\displaystyle \frac{{\mathit{\sigma }}_{\widehat{p}}}{{\mathit{\mu }}_{\widehat{p}}}}\times \mathit{\omega },\end{array}$$

where ${\mathit{\mu }}_{\widehat{p}}$ and ${\mathit{\sigma }}_{\widehat{p}}$ are the mean and standard deviation of the probability estimate over the different instances, and ω is the average number of time steps calculated in one realisation. In short, ϵ measures how precise the numerical method is at equal computational cost. The smaller the ϵ, the better the algorithm performs.

The results are shown in Fig. 8b. For the system in Eq. (15), the score functions ϕ_dist, ϕ_gauss and ϕ_𝒞 have little difference in performance. For the lowest noise values, the path-based score function has at most a 30 % smaller error ϵ than the score function ϕ_gauss. Changing the decay length d₀ hardly changes the error ϵ. When adjusting the parameters of the potential V or applying the same method to the triple-well system used in Rolland and Simonnet (2015), the performance gain, while often present, never systematically exceeded 30 %. All in all, in this category of two-dimensional systems, using the path-based score function approach yields little improvement.

Figure 8 (a) Transition probability p as a function of noise parameter σ. Mean estimates and interquartile range (error bars) over N_samples=10 independent realisations of the TAMS algorithm using the score functions ϕ_dist (green down triangle), ϕ_gauss (with β=1.5, blue circle) and ϕ_𝒞 with decay parameters $d_{\mathrm{0}}=\mathrm{2},\mathrm{20},\mathrm{200}$ (yellow square, pink diamond, purple up triangle) are shown. The Monte Carlo estimate (black cross) has been run with $N=\mathrm{5}\times {\mathrm{10}}^{\mathrm{5}}$ trajectories with the target sets coinciding with the confidence ellipsoids around X_B, with associated standard deviation ${\mathit{\sigma }}_p=\sqrt{p(\mathrm{1}-p)/N}$ (error bars). There is an overall good agreement between the numerical estimations. The path-based score function is robust to the choice of decay parameter. (b) Performance of the score functions measured by the work-normalised error ϵ as a function of noise σ (same markers as panel a). At most, there is a 30 % decrease in error when using the path-based score function ϕ_𝒞 at low noise. Note that at high noise σ=3, using a short decay length d₀=2 with ϕ_𝒞 leads to poor performance because the greater values of ϕ_𝒞 are tightly concentrated around the estimated instanton, whereas typical transition paths are not necessarily clustered around it. Otherwise, the performances of the score functions are roughly similar.

3.3 Transition probabilities in a box model of the AMOC

Finally, as a main application of the modification of the target set B into a confidence ellipsoid ℰ, we consider the dynamical system in Castellana et al. (2019), which represents a box model of the AMOC. The system consists of a set of stochastic differential equations, plus one algebraic constraint, the latter representing the salt conservation in the model (as shown in Appendix B). This model can be formulated as

$$\begin{array}{ll}{\displaystyle }\mathrm{d}{\mathit{Y}}_t& {\displaystyle }={\mathit{F}}_{\mathbf{1}}({\mathit{Y}}_t,Z_t)\mathrm{d}t+\mathbf{G}\mathrm{d}{\mathbf{W}}_t,\\ \text{(18)}& {\displaystyle }\mathrm{0}& {\displaystyle }={\mathit{F}}_{\mathbf{2}}({\mathit{Y}}_t,Z_t).\end{array}$$

In the equations above, we split the state of the system X_t into two parts: Y_t, which includes salinities of four of the boxes plus the depth of the pycnocline D, and Z_t, which represents the salinity of the deep box (S_d). As the noise is applied only to the asymmetric component of the atmospheric freshwater flux (E_a), it directly affects only two of the variables: S_n and S_s. Moreover, the stochastic increments associated with the two variables are identical to make sure that each decrease in freshwater forcing in the southern box results in the same increase in it in the northern box in the model. Therefore, the noise is not spatially independent, and the noise matrix G is no longer diagonal: it consists of a (5×1) row vector, with only two elements different from zero. For a reasonable choice of the parameters, the deterministic system is in a bistable regime (Castellana et al., 2019), which means that there are two possible equilibrium states under the same forcing conditions. In general, we are interested in studying transitions between the present-day AMOC (X_A) and the collapsed state (X_B).

For a differential-algebraic system of equations (DAEs), such as the system in Eq. (18), we need to be particularly careful while computing the covariance ellipsoids. First of all, we make use of the Schur complement of the Jacobian of the system, which allows us to calculate the covariance matrix when an algebraic constraint is present (Baars et al., 2017). Nevertheless, the resulting matrix C_B is singular, with two eigenvalues being equal to zero. One of the corresponding eigenvectors is a vector pointing in the direction of the variable D (depth of the pycnocline). The reason behind it is that the differential equation governing the evolution of D does not contain any of the other variables (see Castellana et al., 2019) when the system is in the collapsed state (X_B). This results in D not being affected by the noise, as this is imposed only on two of the salinities. Hence, we compute the covariance matrix ${\mathbf{C}}_B^S$ relative to the salinities, removing 1 degree of freedom from the original matrix. Unfortunately, such a matrix still has one zero eigenvalue, which is due to the salinity conservation (the algebraic equation in the system Eq. 18). To overcome this problem, we compute the Moore–Penrose inverse (or pseudo-inverse) of the covariance matrix, ${\mathbf{C}}_B^{S+}$, by performing a singular value decomposition of ${\mathbf{C}}_B^S$ and removing the zero eigenvalue, together with the corresponding eigenvector (Ben-Israel and Greville, 2003). A two-dimensional projection of the ellipsoid constructed for the box model is shown in Fig. 9.

For the system in Eq. (18), the modified score function is more complicated, as the covariance matrix used to construct the ellipsoid contains only the degrees of freedom related to the salinities of the model, leaving the variable D (depth of the pycnocline) out. From a geometric point of view, that means that the covariance ellipsoid around X_B is degenerated along the D direction. Figure 9 shows a projection of the ellipsoid – once the noise amplitude is fixed – onto the plane identified by the variables S_n and S_s (respectively, the salinity of the northern box and the one of the southern box in the model). The projection was obtained by calculating the conditional covariance matrix of the two variables under consideration, given that the other variables are set on their mean value (Wasserman, 2013). The two-dimensional ellipse contains the large majority of the projected time points (on the same plane) of a trajectory that wanders around the equilibrium. By construction, the confidence level of the confinement is higher than the one prescribed for the full-dimensional ellipsoid (in this case 0.95). Clearly the level sets of the score function ϕ_gauss do not coincide with the one of the ellipsoid. Hence, the importance of modifying the score function appears evident.

Figure 9 Level sets of the score function ϕ_gauss in proximity to the collapsed equilibrium state X_B of the system in Castellana et al. (2019), projected onto the plane identified by the variables S_n and S_s (respectively, the salinity of the northern box and the one of the southern box in the model). A trajectory of the system, initiated in the collapsed state, has been projected onto the same plane (blue circles). The red ellipse is the two-dimensional projection of the covariance ellipsoid constructed by using the matrix ${\mathbf{C}}_B^S$ and a confidence level of 0.95; 98 % of the points comprising the trajectory are inside the ellipsoid.

When constructing an improved score function, in order to evaluate whether a state belongs to the neighbourhood of X_B, we need to check two conditions: (i) whether the state of the system is inside the salinity covariance ellipsoid drawn around the destination equilibrium and (ii) whether the variable D of the state is the same as the one of X_B. Hence, the improved score function for the box model reads as

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{box}}\left(\mathit{X}\right)\equiv \left\{\begin{array}{ll}\mathrm{1}& \text{if}{\mathit{X}}^S\in {\mathcal{E}}^S\text{and}{\mathit{X}}^D={\mathit{X}}_B^D,\\ {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}& \text{if}{\mathit{\varphi }}_{\text{gauss}}\left(\mathit{X}\right)>{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\text{and}\\ & ({\mathit{X}}^S\notin {\mathcal{E}}^S\text{or}{\mathit{X}}^D\ne {\mathit{X}}_B^D),\\ {\mathit{\varphi }}_{\text{gauss}}\left(\mathit{X}\right)& \text{otherwise},\end{array}\right.\end{array}$$

where X^S indicates the part of the state vector representing the set of the salinities, whereas X^D represents the variable D. As already mentioned, to check whether a certain state belongs to the salinity ellipsoid, we made use of the pseudo-inverse of ${\mathbf{C}}_B^S$ in the definition of Eq. (12).

To be able to assess the relevance of a proper definition of the target set in the TAMS algorithm and hence the importance of using the improved version of the score function, we computed transition probabilities of the AMOC from the present-climate state to the collapsed state for reasonable values of the atmospheric forcing and noise. In particular, we chose ${\overline{E}}_{\mathrm{a}}=\mathrm{0.20}\mathrm{Sv}$ and f_σ=0.16. This last value represents the ratio between the standard deviation of the noise in the atmospheric forcing and its mean value (Castellana et al., 2019). The number of trajectories used in the algorithm was set to 100, and the timescales at which the probabilities were evaluated were chosen to be 1000 and 10 000 years, respectively. For each transition probability, we used three different versions of the algorithm, based on three different settings for the score function. The first two versions were implemented by using ϕ_gauss in Eq. (6), with two different choices for the threshold ϕ^target. In the third version, we used ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{box}}$ as defined in Eq. (19), where the starting score function was ϕ_gauss. Each probability was calculated by running 15 instances of the algorithm and then computing the median and the interquartile range (IQR). The results are shown in Table 1 below.

Table 1 Results of the transition probabilities for the AMOC model using different score functions.

When running TAMS to compute transition probabilities between two states of the Atlantic circulation in this model, with different versions of the score function, we found a considerable discrepancy between the obtained values. In particular, it appears that setting a very high threshold in the Gaussian score function makes the algorithm detect no transitions (we set up the algorithm so that it stops when the probabilities involved are smaller than 10⁻⁹): the reason behind this is that, because of the presence of the noise, we do not expect the state of the system to stay indefinitely close to the destination equilibrium but rather to wander around it. Therefore, the score function, which assigns the maximum score only to a very small neighbourhood of the equilibrium, is not able to properly recognise transitions. Moreover, it is not surprising that, when using ϕ_gauss with a smaller value of ϕ^target (0.99) or ${\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\text{box}}$, we obtain different results: as the shape of the covariance ellipsoid is not spherical (see Fig. 9), we expect ϕ_gauss to detect transitions even though the state is actually still far from the destination equilibrium, at least in certain directions. As a general rule, we expect ϕ_gauss to give results different from ϕ_box as long as the ellipsoid of the system is not spherical, regardless of the choice of the threshold ϕ^target.

4 Summary and discussion

We presented and applied several improvements to the TAMS rare-event algorithm, when used to compute transitions in multistable systems. The first improvement was based on a more rigorous criterion to define noise-induced transitions involving confidence ellipsoids ℰ to formalise this criterion. In turn, this led to the rigorous choice of the target set B=ℰ_B, which was traditionally set by rather arbitrary thresholds. We then showed how to incorporate this definition of B into the score function $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$. For certain classes of systems, like the ones containing algebraic constraints in addition to differential equations or when the noise does not affect one or more directions in the variable space (i.e. the associated covariance matrix is singular), this method requires some precautions. In particular, for the box model in Castellana et al. (2019), we had to adapt the definition of the improved score function in Eq. (19) as well as calculate the pseudo-inverse matrix of the covariance matrix in order to compute the ellipsoid.

This method, while quite general, has several limitations. While the modified score function $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$ is continuous outside B, it is constant in the domain $M=\mathit{\{}\mathit{X}\in {\mathbb{R}}^n\mid \mathit{\varphi }(\mathit{X})>{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}^{\text{target}}\mathit{\}}\setminus \mathcal{E}$ (see Fig. 4). This means that in the TAMS algorithm, branching will never occur inside M, but at the boundary ∂M. This can have an influence on the convergence of TAMS if the level sets of the initial score function ϕ have a pathological shape near X_B and the spatial extension of D is not negligible. Nevertheless, we expect this to have little impact because this phenomenon is localised near the target state X_B. Therefore, the trajectories will naturally converge towards X_B as a result of the dynamics, even without the help of the branching process of TAMS. However, to ensure that the confidence ellipsoid ℰ defines a meaningful target set, one needs to be sure that ℰ is contained inside the basin of attraction of the target state X_B. While this is the case in the limit of small noise σ→0, it might not be the case for finite noise. A solution to this issue would be to compute the basin of attraction 𝒱 of X_B and define the target set B as the intersection $B=\mathcal{E}\cap \mathcal{V}$. However, we expect that, in the generic case, this occurs when the noise level σ is high enough so that transitions are less rare and that a direct Monte Carlo estimation is sufficient to estimate transition probabilities.

Next, we proposed a systematic method of defining a score function, designed to approximate the static committor, based on empirical transition paths. We proposed an algorithm to estimate the typical transition path under a high noise level, which is then used to define a family of score functions with a single decay parameter d₀. We applied our method to a two-dimensional well with a potential barrier. We found that our typical path estimation gave satisfactory results and that the associated score function ϕ_𝒞, while discontinuous, remained unbiased and relatively insensitive to the change in the decay parameter d₀. While we did not find significant performance improvements over existing non-trivial score functions, we think that differences will become apparent if applied to higher-dimensional systems, where there are more directions to fluctuate in. We did not show any results of the application of this method to the box model of the AMOC (Castellana et al., 2019), because for this rather complicated case, we were not able to efficiently compute typical transition paths from the trajectory histograms of the system.

One key limitation of our approach of constructing the path-based score function ϕ_𝒞 is the computer memory needed to store the trajectory histogram, which becomes prohibitively huge for high-dimensional systems such as discretised partial differential equation (PDE) systems. As an example, a 50×50 two-dimensional grid storing four variables in each cell (e.g. two velocity components, pressure and a tracer) with 10 bins of resolution in each degree of freedom would require more than 10 PB of memory, which is unfeasible. However, this limitation can be easily bypassed by defining the objects needed to run the TAMS algorithm, namely the score function ϕ and the target set B, in a reduced space of much fewer dimensions. For instance, Bouchet et al. (2019) studied the dynamics of the barotropic beta-plane quasi-geostrophic equations describing Jupiter's turbulent atmosphere. While the PDE system was evolved in the full phase space, their rare-event algorithm was run in a reduced three-dimensional phase space defined by three Fourier coefficients. The target set B and the score function were defined on this reduced space. Moreover, they accumulated transition trajectories in a three-dimensional histogram and showed their concentration around instantons. By applying the path-finding algorithm, an empirical estimation of the instanton can be made. This offers a viable alternative to solving a minimisation problem in the full space to compute the instanton and then project it in the reduced space, which is next to intractable for this system.

Another way to define the reduced space 𝒱 in which to run TAMS is to consider the principal components, also called empirical orthogonal functions (EOFs), which are the eigenvectors of the covariance matrix. One idea, suggested by Baars (2019), is to retain the EOFs $\mathit{\{}{\mathit{Y}}_{\mathrm{1}}^A,\mathrm{\dots },{\mathit{Y}}_k^A,{\mathit{Y}}_{\mathrm{1}}^B,\mathrm{\dots },{\mathit{Y}}_{k^{\prime }}^B\mathit{\}}$ with the largest variance (i.e. the eigenvalue of the covariance matrix) at the initial state X_A and target state X_B. EOFs represent the directions in which the system fluctuates the most. They are then assumed to be the directions which capture best the noise-driven dynamics. When studying transitions in a PDE model of the Atlantic Meridional Overturning Circulation, den Toom et al. (2011) and Baars et al. (2021) projected the dynamics in a reduced space $\mathcal{W}\equiv \text{Span}\mathit{\{}{\mathit{X}}_A,{\mathit{X}}_B,{\mathit{Y}}_{\mathrm{1}}^A,\mathrm{\dots },{\mathit{Y}}_k^A,{\mathit{Y}}_{\mathrm{1}}^B,\mathrm{\dots },{\mathit{Y}}_{k^{\prime }}^B\mathit{\}}$ of a dimension (≈500) still too large to apply the histogram method directly. However, one idea is that the TAMS algorithm could be run in an even smaller space – $\mathcal{V}\equiv \text{Span}\mathit{\{}{\mathit{X}}_A,{\mathit{X}}_B,{\mathit{Y}}_{\mathrm{1}}^A,\mathrm{\dots },{\mathit{Y}}_d^A,{\mathit{Y}}_{\mathrm{1}}^B,\mathrm{\dots },{\mathit{Y}}_{d^{\prime }}^B\mathit{\}}$ (of dimension<10) – while still computing the dynamics in the space 𝒲. Then, the memory required to store a histogram becomes reasonable and our histogram method can be applied.

Another potential issue of our modified TAMS method is the fact that the score function ϕ_𝒞 is discontinuous because the trajectory has positive curvature, as shown in Fig. 7b, blue path. The discontinuity is located near the axis x=0. Indeed, when crossing the axis x=0, the closest point on 𝒞 changes and s(X,𝒞) is discontinuous. In fact, in the mathematical proofs about the statistical and convergence properties of the probability estimator (Cérou et al., 2016), the score functions are assumed to be continuous. Nevertheless, in our applications, we did not detect any statistically significant bias in the probability estimator due to the discontinuity. Moreover, some meaning can be attributed to the discontinuity: it is located at the boundary between the attraction basins of X_A and X_B, and it thus reflects a qualitative change in behaviour in the system. Crossing this boundary means that the trajectory converges to X_B instead of X_A if σ=0. In addition, the remnant of a discontinuity is observed in the static committor of the similar triple-well system used in Rolland and Simonnet (2015). Indeed, their Fig. 4c shows the contour plots of the static committor in the low-noise regime. A steep gradient is present at the x=0 boundary, which gives further evidence that the discontinuity of ϕ_𝒞 may not be problematic.

Further testing of the ideas presented in this work in high-dimensional systems such as discretised PDEs would give more insight as to the effectiveness of our approach compared to the more generic score functions used up to now. Moreover, incorporating some form of time dependence into the score function ϕ to specifically optimise TAMS would constitute an interesting project.

Appendix A: Trajectory-Adaptive Multilevel Sampling (TAMS) algorithm

Appendix B: Box model for the Atlantic Meridional Ocean Circulation

The equations determining the evolution of the AMOC in this model are the salinity budgets of the different boxes, together with the variation of the volume of the pycnocline, and the salt and volume conservation equations (Castellana et al., 2019):

$$\begin{array}{llll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}\left(V_{\mathrm{t}}{S}_{\mathrm{t}}\right)}{\mathrm{d}t}}& {\displaystyle }& {\displaystyle }=q_{\mathrm{S}}\left(\mathit{\theta }\right(q_{\mathrm{S}})S_{\text{ts}}+\mathit{\theta }(-q_{\mathrm{S}}\left)S_{\mathrm{t}}\right)\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }+q_{\mathrm{U}}{S}_{\mathrm{d}}-\mathit{\theta }\left(q_{\mathrm{N}}\right)q_{\mathrm{N}}{S}_{\mathrm{t}}\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }+r_{\mathrm{S}}(S_{\text{ts}}-S_{\mathrm{t}})\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }+r_{\mathrm{N}}(S_{\mathrm{n}}-S_{\mathrm{t}})+\mathrm{2}{E}_{\mathrm{s}}{S}_{\mathrm{0}},\\ {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left(V_{\text{ts}}{S}_{\text{ts}}\right)}{\mathrm{d}t}}& {\displaystyle }& {\displaystyle }=q_{\text{Ek}}{S}_{\mathrm{s}}-q_{\mathrm{e}}{S}_{\text{ts}}-q_{\mathrm{S}}\left(\mathit{\theta }\right(q_{\mathrm{S}})S_{\text{ts}}\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }+\mathit{\theta }(-q_{\mathrm{S}})S_{\mathrm{t}})+r_{\mathrm{S}}(S_{\mathrm{t}}-S_{\text{ts}}),\\ {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left(V_{\mathrm{n}}{S}_{\mathrm{n}}\right)}{\mathrm{d}t}}& {\displaystyle }& {\displaystyle }=\mathit{\theta }\left(q_{\mathrm{N}}\right)q_{\mathrm{N}}(S_{\mathrm{t}}-S_{\mathrm{n}})+r_{\mathrm{N}}(S_{\mathrm{t}}-S_{\mathrm{n}})\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }-(E_{\mathrm{s}}+E_{\mathrm{a}})S_{\mathrm{0}},\\ {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left(V_{\mathrm{s}}{S}_{\mathrm{s}}\right)}{\mathrm{d}t}}& {\displaystyle }& {\displaystyle }=q_{\mathrm{S}}\left(\mathit{\theta }\right(q_{\mathrm{S}})S_{\mathrm{d}}+\mathit{\theta }(-q_{\mathrm{S}}\left)S_{\mathrm{s}}\right)+q_{\mathrm{e}}{S}_{\text{ts}}\\ {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }-q_{\text{Ek}}{S}_{\mathrm{s}}-(E_{\mathrm{s}}-E_{\mathrm{a}})S_{\mathrm{0}},\\ {\displaystyle }& {\displaystyle }\left(A+{\displaystyle \frac{L_{\text{xA}}{L}_y}{\mathrm{2}}}\right){\displaystyle \frac{\mathrm{d}D}{\mathrm{d}t}}& {\displaystyle }& {\displaystyle }=q_{\mathrm{U}}+q_{\text{Ek}}-q_{\mathrm{e}}-\mathit{\theta }\left(q_{\mathrm{N}}\right)q_{\mathrm{N}},\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{0}}{V}_{\mathrm{0}}& {\displaystyle }& {\displaystyle }=V_{\mathrm{n}}{S}_{\mathrm{n}}+V_{\mathrm{d}}{S}_{\mathrm{d}}+V_{\mathrm{t}}{S}_{\mathrm{t}}\\ \text{(B1)}& {\displaystyle }& {\displaystyle }& {\displaystyle }& {\displaystyle }+V_{\text{ts}}{S}_{\text{ts}}+V_{\mathrm{s}}{S}_{\mathrm{s}},\end{array}$$

where the function θ(x) is the Heaviside step function. The transports depend on the variables via the following relations:

$$\begin{array}{ll}{\displaystyle }{q}_{\text{Ek}}& {\displaystyle }={\displaystyle \frac{\mathit{\tau }{L}_{\text{xS}}}{{\mathit{\rho }}_{\mathrm{0}}\left|f_{\mathrm{S}}\right|}},\\ {\displaystyle }{q}_{\mathrm{e}}& {\displaystyle }=A_{\text{GM}}{\displaystyle \frac{L_{\text{xA}}}{L_y}}D,\\ {\displaystyle }{q}_{\mathrm{U}}& {\displaystyle }={\displaystyle \frac{\mathit{\kappa }A}{D}},\\ {\displaystyle }{q}_{\mathrm{N}}& {\displaystyle }=\mathit{\eta }{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{n}}-{\mathit{\rho }}_{\text{ts}}}{{\mathit{\rho }}_{\mathrm{0}}}}{D}^{\mathrm{2}},\\ \text{(B2)}& {\displaystyle }{q}_{\mathrm{S}}& {\displaystyle }=q_{\text{Ek}}-q_{\mathrm{e}},\end{array}$$

where the density of the generic box i is defined as

$$\begin{array}{}\text{(B3)}& {\mathit{\rho }}_i={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-\mathit{\alpha }(T_i-T_{\mathrm{0}})+\mathit{\beta }(S_i-S_{\mathrm{0}})\right).\end{array}$$

The volumes depend, in turn, on the variable D:

$$\begin{array}{ll}{\displaystyle }{V}_{\mathrm{t}}& {\displaystyle }=AD,\\ {\displaystyle }{V}_{\text{ts}}& {\displaystyle }={\displaystyle \frac{L_{\text{xA}}{L}_y}{\mathrm{2}}}D,\\ \text{(B4)}& {\displaystyle }{V}_{\mathrm{d}}& {\displaystyle }=V_{\mathrm{0}}-V_{\mathrm{n}}-V_{\mathrm{s}}-V_{\mathrm{t}}-V_{\text{ts}}.\end{array}$$

The reference parameter values are shown in Table B1.

Table B1 Reference parameters used in Eqs. (B1)–(B4) from Castellana et al. (2019).