GitHub

Transition Path Theory

Transition Path Theory (TPT) studies rare transitions between metastable states of an SDE (Vanden-Eijnden, 2006; E and Vanden-Eijnden, 2010; E. and Vanden-Eijnden, 2006). Given two disjoint sets $A$ and $B$, it tells you how, where, and how often a noise-driven trajectory crosses from one to the other — without simulating the (typically exponentially long) waiting times.

TPT is built on top of the discretised generator and the standard boundary-value problems on it (committors, invariant density). Read the diffusion operator manual first if you have not already.

Reactive trajectories and observables

The key objects are the reactive trajectories: pieces of a long path that have just left $A$ and will hit $B$ before returning. TPT characterises this sub-ensemble through three observables, all built from the forward committor $q^{+}$, the backward committor $q^{-}$, and the invariant density $\rho$:

  • Reactive density $\rho_R = \rho\,q^{+}\,q^{-}$ — probability density of being on a reactive piece at a given time.
  • Reactive transition rate $k_{AB}$ — events per unit time, the flux of probability leaving $A$ along reactive trajectories.
  • Reactive current $J_R$ — vector field whose streamlines reveal the dominant transition channels.

For reversible systems $q^{-} = 1 - q^{+}$; in general the two are independent and both must be computed.

The reactive current splits into a reversible (gradient-flow) part and an irreversible (cyclic) part — the latter vanishes for systems satisfying detailed balance and exposes the circulating component of the flow for non-equilibrium systems. The pieces are accessed via reactive_current_reversible and reactive_current_irreversible; they sum to the full reactive_current.

When TPT is the right tool

Use TPT when you have a time-homogeneous SDE with metastable wells, and you want spatial information (where reactive paths concentrate, the committor surface, transition channels) — not just a scalar rate. Practical sweet spot: dimensionless barrier $\Delta U / \sigma^2$ between $\sim 1$ and $10$, state-space dimension $\lesssim 3$.

Look elsewhere when:

  • Very high barriers ($\Delta U/\sigma^2 \gg 10$): use sgMAM / string method (minimize_action, string_method) — they compute the most likely path and its rate directly.
  • Very shallow barriers: just simulate (transitions).
  • Only need a rate: use mean_first_passage_time or Kramers asymptotics.
  • Not a diffusion (jumps, fractional noise, deterministic chaos): outside the SDE framework on which TPT is built.
  • High-dimensional ($D \gtrsim 4$) with no reaction coordinate: a uniform grid is infeasible. The standard remedy is to first learn a 1–3D collective variable and then run TPT in CV space — outside this package.
  • Non-autonomous drift: not implemented here.

A useful sanity check is to compare the TPT rate against an independent estimate (Kramers, simulation, or large-deviation action).

The API

ReactiveTransition runs the full pipeline ($\rho$, both committors, the rate) in one call and caches the results.

result = ReactiveTransition(sys, grid, A, B)         # or
result = ReactiveTransition(gen, A, B)               # reuse a generator
result = ReactiveTransition(sys, grid, A, B; reverse = sys_rev)

reactive_rate(result)                # cached
reactive_density(result)
J_nodes, J_faces = reactive_current(result)
probability_reactive(result)
probability_last_A(result)

forward_committor, backward_committor, and stationary_distribution also dispatch on a ReactiveTransition and return the cached fields. Predicates / BitVector / linear-index masks are accepted for A and B, just as on a DiffusionGenerator.

CriticalTransitions.ReactiveTransition — Type
struct ReactiveTransition{D, BC, T<:AbstractFloat}

Transition-path-theory result for transitions from set A to set B.

ReactiveTransition stores the generator, set masks, invariant density, forward and backward committors, discrete adjoint generator, and reactive rate. It represents the stationary ensemble of trajectories that most recently left A and will hit B before returning to A.

Fields

  • generator::DiffusionGenerator{D, BC} where {D, BC}: The diffusion generator the analysis was built on.

  • A_mask::BitVector: Boolean mask selecting cells in set A (length ncells(grid)).

  • B_mask::BitVector: Boolean mask selecting cells in set B (length ncells(grid)).

  • ρ::Vector{T} where T<:AbstractFloat: Invariant probability density (dot(ρ, weights) = 1).

  • qplus::Vector{T} where T<:AbstractFloat: Forward committor q⁺ (probability of reaching B before A).

  • qminus::Vector{T} where T<:AbstractFloat: Backward committor q⁻.

  • Qadj::SparseArrays.SparseMatrixCSC{T, Int64} where T<:AbstractFloat: Discrete adjoint generator w.r.t. ρ; used by reactive_current_{reversible,irreversible}.

  • rate::AbstractFloat: A → B reactive transition rate (cached).

  • physical_reverse::Bool: True if qminus came from a user-supplied physical reverse generator; false if from the discrete adjoint.

Construction

ReactiveTransition(gen::DiffusionGenerator, A, B; alg=UMFPACKFactorization(), reverse=nothing)
ReactiveTransition(sys::CoupledSDEs, grid, A, B; bc=Reflecting(), reverse=nothing)

The first form takes a pre-built DiffusionGenerator — useful for sweeping multiple (A, B) pairs without rediscretising. The second is a convenience that builds the generator from sys internally.

A and B accept:

  • a predicate x::SVector{D} -> Bool evaluated at each cell center,
  • a BitVector of length ncells(grid), or
  • a vector of linear cell indices.

If reverse is provided (a CoupledSDEs for the high-level form, a DiffusionGenerator for the generator form), the backward committor is computed on that physical reverse-drift generator instead of the discrete adjoint.

res = ReactiveTransition(gen, x -> x[1] < -1, x -> x[1] > 1)
k = reactive_rate(res)
J_nodes, J_faces = reactive_current(res)

Default alg is UMFPACKFactorization() (sparse direct LU) for all internal solves. Pass any SciMLLinearSolveAlgorithm (e.g. KrylovJL_GMRES()) for an iterative solver — useful for large grids where direct LU runs out of memory. Further kwargs flow to LinearSolve.solve.

See also forward_committor, backward_committor, reactive_density, and reactive_current.

source
CriticalTransitions.reactive_rate — Function
reactive_rate(r::ReactiveTransition) -> AbstractFloat

A → B reactive transition rate $k_{AB}$: events per unit time, the rate at which probability leaves A along reactive trajectories.

In CTMC form,

k_{AB} = ∑_{i ∈ A, j ∉ A} ρ[i] · v · q⁻[i] · Q[i, j] · q⁺[j] ,

with v = cell_volume(grid).

source
CriticalTransitions.reactive_current — Function
reactive_current(r::ReactiveTransition) -> Tuple{Any, Any}

Reactive probability current as (J_nodes, J_faces):

  • J_faces :: NTuple{D, Array{Float64, D}} — net signed flux across each axis-aligned face, divided by the corresponding axis spacing. Element k has shape nbox along all axes except axis k, where it has nbox[k] - 1 entries for reflecting / absorbing boundary conditions (one per interior face), or nbox[k] entries for periodic boundary conditions (the extra slot is the wraparound face from cell nbox[k] back to cell 1).
  • J_nodes :: NTuple{D, Array{Float64, D}} — face fluxes averaged onto cell centers along each axis; element k has shape nbox.

The face flux from n to m is μ[n] · q⁻[n] · Q[n, m] · q⁺[m], where μ = ρ · v and Q is the rate matrix.

For non-equilibrium systems the current splits into a reversible (gradient-flow) part and an irreversible (cyclic) part — see reactive_current_reversible and reactive_current_irreversible to split gradient-flow and cyclic components.

source
CriticalTransitions.reactive_current_reversible — Function
reactive_current_reversible(
    r::ReactiveTransition
) -> Tuple{Any, Any}

Reversible (gradient-flow) part of the reactive current. Built from the symmetric part of the generator, Q_sym = (Q + Q̃)/2, where Q̃ is the discrete adjoint with respect to the invariant density. For a reversible system this equals the full reactive_current.

source
CriticalTransitions.reactive_current_irreversible — Function
reactive_current_irreversible(
    r::ReactiveTransition
) -> Tuple{Any, Any}

Irreversible (cyclic / time-asymmetric) part of the reactive current. Built from the antisymmetric part of the generator, Q_antisym = (Q − Q̃)/2. Vanishes for systems satisfying detailed balance; for non-equilibrium systems it exposes the circulating component of the reactive flow.

Returns (J_nodes, J_faces) with the same shapes as reactive_current.

source
CriticalTransitions.probability_last_A — Function
probability_last_A(r::ReactiveTransition{D, BC, T}) -> Any

Total probability that the most recent metastable visit was set A.

The complementary probability is the chance that the most recent visit was B.

source

Complete example

A 2D Maier–Stein-like system with stable wells at $(\pm 1, 0)$:

using CriticalTransitions
using LinearAlgebra: norm

b(u, p, t) = SA[u[1] - u[1]^3 - 10*u[1]*u[2]^2, -(1 + u[1]^2)*u[2]]
sys = CoupledSDEs(b, [0.0, 0.0]; noise_strength = 0.3)

grid = CartesianGrid((-1.6, 1.6, 121), (-1.0, 1.0, 81))
A = x -> norm(x .- SA[-1.0, 0.0]) < 0.25
B = x -> norm(x .- SA[ 1.0, 0.0]) < 0.25

result = ReactiveTransition(sys, grid, A, B)
@show reactive_rate(result)
J_nodes, _ = reactive_current(result)

To sweep over $(A, B)$ choices without re-discretising:

gen = DiffusionGenerator(sys, grid)
for r in (0.2, 0.3, 0.4)
    A = x -> norm(x .- SA[-1.0, 0.0]) < r
    B = x -> norm(x .- SA[ 1.0, 0.0]) < r
    println("r = $r → k_AB = ", reactive_rate(ReactiveTransition(gen, A, B)))
end

A fully worked, plotted version lives in examples/transition_path_theory_double_well.jl.

Tips

  • Choose $A$ and $B$ generously around the wells. Too small → noisy; overlapping the saddle → distorted current. A ball of radius $\sim \sqrt{\sigma^2/\kappa}$ (with $\kappa$ the well curvature) is a reasonable default.
  • The rate is in SDE time units. Dividing by an appropriate occupation probability gives the mean transition time; see (E and Vanden-Eijnden, 2010) for the precise relation.
  • Sweeping $(A, B)$? Build DiffusionGenerator once and reuse it.
  • Small-noise regime. As $\sigma \to 0$ the committor develops a thin transition layer near the saddle — uniform grids must be globally fine, which is wasteful. sgMAM / string method are the better tool there.

References