---
title: Isostable Framework for Nonlinear Dynamics
url: https://www.emergentmind.com/topics/isostable-coordinate-framework
type: topic
---

# Isostable Framework for Nonlinear Dynamics

Searching arXiv for the specified papers to ground the article in current metadata and citations.
arxiv_search.query({"search_query":"id:2507.13997 OR id:2102.04526 OR id:2105.11048 OR id:2310.02725","max_results":10,"sort_by":"relevance","sort_order":"descending"})
arxiv_search.query({"search_query":"ti:\"Identification and Computation of Slow Manifolds Using the Isostable Coordinate System\" OR ti:\"Data-Driven Inference of High-Accuracy Isostable-Based Dynamical Models in Response to External Inputs\" OR ti:\"Isostables for stochastic oscillators\" OR ti:\"Insights into oscillator network dynamics using a phase-isostable framework\"","max_results":10,"sort_by":"submittedDate","sort_order":"descending"})
The isostable coordinate framework is a coordinate description of nonlinear dynamics in which state-space evolution is organized by asymptotic decay rates associated with Koopman eigenfunctions or, for periodic orbits, Floquet modes. For systems with a stable fixed point, the framework focuses on the slowest decaying principal Koopman eigenmodes and defines scalar coordinates whose level sets are isostables; for stochastic oscillators and deterministic limit cycles, closely related constructions produce amplitude-like coordinates complementary to asymptotic phase. In the formulations developed for fixed-point attractors, external-input reductions, stochastic oscillators, and oscillator networks, isostable coordinates provide a common language for invariant geometry, reduced-order modeling, and bifurcation analysis [2507.13997] [2102.04526] [2105.11048] [2310.02725].

## 1. Koopman-theoretic definition and local structure

For a smooth autonomous ODE
$$
\dot x = F(x), \qquad x\in\mathbb{R}^n,
$$
with hyperbolic fixed point \(x_0\), the Koopman operator acts on observables \(g(x)\) by
$$
[K^t g](x) = g\circ \Phi^t(x),
$$
where \(\Phi^t(x)\) is the flow. Although the state-space dynamics are nonlinear, \(K^t\) is linear but generally infinite-dimensional. A Koopman eigenfunction \(\phi_j(x)\) and eigenvalue \(\lambda_j\in\mathbb{C}\) satisfy
$$
\phi_j(\Phi^t(x)) = e^{\lambda_j t}\phi_j(x).
$$
When the Jacobian \(J=DF(x_0)\) has eigenvalues ordered by
$$
|\operatorname{Re}\lambda_1| \le |\operatorname{Re}\lambda_2| \le \cdots \le |\operatorname{Re}\lambda_N|,
$$
the associated \(\phi_1,\dots,\phi_N\) are the principal eigenfunctions, and the real part or modulus of \(\phi_j\) yields an isostable coordinate measuring distance along the \(j\)-th mode [2507.13997].

A constructive limit definition is obtained by fixing left and right eigenvectors \(w_j,v_j\) of \(J\), normalized so that \(w_j^T v_k=\delta_{jk}\), and setting
$$
\psi_j(x) = \lim_{t\to\infty} e^{-\lambda_j t}\, w_j^T[\Phi^t(x)-x_0].
$$
Its level sets \(\{\psi_j(x)=\mathrm{const}\}\) are the isostables associated with decay rate \(\lambda_j\). Differentiation along trajectories gives
$$
\dot \psi_j = \lambda_j \psi_j,
$$
so each principal isostable decays or grows exponentially at rate \(\lambda_j\) [2507.13997].

A closely related normalization appears in isostable reduction for forced systems. There, Koopman eigenfunctions \(\phi_i\) satisfy the spectral PDE
$$
\nabla \phi_i(x)\cdot f(x)=\lambda_i \phi_i(x), \qquad \operatorname{Re}\lambda_i<0,
$$
with normalization
$$
\phi_i(x_0)=0, \qquad \nabla\phi_i(x_0)\cdot v_j=\delta_{ij},
$$
and the real-valued isostable coordinate \(\sigma_i(x)\) is defined by filtering out all but the slowest decaying Koopman mode through a limit of the same form [2102.04526]. Near the fixed point, these definitions recover the linear eigendirections because \(\psi_k(x)\simeq w_k^T(x-x_0)\) locally [2507.13997].

## 2. Slow manifolds as zero sets of fast isostables

When there is a pronounced spectral gap at index \(\beta\),
$$
|\operatorname{Re}\lambda_\beta| \ll |\operatorname{Re}\lambda_{\beta+1}|,
$$
the framework defines a \(\beta\)-dimensional slow manifold by setting all faster isostable amplitudes to zero:
$$
W^s = \{x\in\mathbb{R}^n \mid \psi_k(x)=0 \ \text{for}\ k=\beta+1,\dots,N\}.
$$
This identifies the slow manifold as the set on which the fastest decaying principal Koopman coordinates vanish [2507.13997].

The definition immediately implies invariance. On \(W^s\), the coordinates \(\psi_{\beta+1},\dots,\psi_N\) are identically zero, so their derivatives vanish under \(\dot \psi_j=\lambda_j\psi_j\), and the flow remains on \(W^s\). Near \(x_0\), the manifold reduces to \(\operatorname{span}\{v_1,\dots,v_\beta\}\), the span of the slow eigenvectors. Away from the linear regime, \(W^s\) is a nonlinear manifold of codimension \(N-\beta\), locally transverse to the fast eigendirections [2507.13997].

This formulation is not merely a local tangent-space statement. The motivating point of the construction is that the condition \(\psi_{k> \beta}=0\) remains meaningful beyond the linear neighborhood of the attractor, so the same coordinate language used for spectral decomposition also gives a geometric definition of the slow manifold in the nonlinear regime [2507.13997]. This suggests a direct bridge between Koopman spectral objects and model order reduction: the reduced dynamics are obtained by retaining the slowest coordinates and eliminating fast amplitudes by an invariant constraint.

## 3. Backward-time computation and numerical strategies

Direct backward integration of the original ODE from a point on \(W^s\) is numerically unstable once the trajectory leaves a small neighborhood of \(x_0\), because any small error in a fast isostable \(\psi_{j>\beta}\) blows up like \(e^{|\operatorname{Re}\lambda_j|t}\). To avoid this, the framework rewrites the backward dynamics in isostable coordinates. With backward-time variable \(\tau=-t\), one has
$$
\frac{d\psi_j}{d\tau} = \nabla\psi_j(x)^T \frac{dx}{d\tau} = I_j^T \frac{dx}{d\tau},
$$
where \(I_j(x)=\partial \psi_j/\partial x\). Stacking \(I_1,\dots,I_N\) into a matrix and enforcing \(\psi_{\beta+1}=\cdots=\psi_N\equiv 0\) on \(W^s\) yields
$$
\frac{dx}{d\tau} = [I_1\ \cdots\ I_N]^{-1}\,[ -\lambda_1\psi_1\ \cdots\ -\lambda_\beta\psi_\beta,\ 0,\dots,0]^T.
$$
If \(I_j\) and \(\psi_j\) are known along \(W^s\), this evolution marches backward without exciting fast modes [2507.13997].

For the slow coordinates \(j\le \beta\), the isostable gradients satisfy the adjoint variational equation
$$
\dot I_j = -[J(x(t))^T-\lambda_j I]\,I_j,
$$
with initial condition \(I_j(x_0)=w_j\). Because these are slow modes, backward integration remains accurate over substantial intervals. The fast gradients \(I_{j>\beta}\), however, are not computed directly. Instead, the computation uses their duals \(g_k\), which solve
$$
\dot g_k = [J(x(t))-\lambda_k I]\,g_k,
$$
with bi-orthonormality \(I_j^T g_k=\delta_{jk}\). Since the backward formula only needs the subspace spanned by \(I_{\beta+1},\dots,I_N\), it is enough to approximate the orthogonal complement of \(\operatorname{span}\{g_1,\dots,g_\beta\}\) [2507.13997].

Two strategies are developed. The asymptotic-expansion approach expands \(x\) and \(g_k\) in Taylor series in the isostables near \(x_0\), for example
$$
x = x_0 + \sum \psi_j v_j + \sum \psi_j\psi_k h^{jk}+\cdots,
$$
and then matches with \(dx/dt=\sum \lambda_j\psi_j g_j\) to solve for \(g_j(\psi)\) to arbitrary order. The method works well in low dimension or mild nonlinearity but becomes unwieldy beyond \(4\)–\(8\)-th order. The predictor–corrector approach instead uses the approximation \(\operatorname{span}\{g_1,\dots,g_\beta\}\approx \operatorname{span}\{v_1,\dots,v_\beta\}\), integrates backward using the slow gradients, then corrects the resulting point by solving for a displacement \(\Delta x\) in the fast eigendirections so that the exact isostable-gradient condition is restored [2507.13997].

The reported numerical issues are threefold. Timescale separation makes direct backward integration exponentially unstable in fast modes; the reformulated backward equation avoids this by keeping \(\psi_{j>\beta}\equiv 0\). Stiffness arises in the adjoint equation when \(\lambda_j\) is very negative, motivating a stiff integrator or small steps. Non-uniqueness requires choosing \(\beta\) so that there are no resonances \(\lambda_j+\lambda_k=\lambda_\ell\) for slow indices; this guarantees convergence of the asymptotic expansions [2507.13997].

## 4. Reduced-order models with forcing and data-driven inference

For systems with a stable fixed point driven by an input \(u(t)\) entering through \(g(x)\), an \(M\)-mode isostable reduction keeps only the slowest coordinates \(\sigma=(\sigma_1,\dots,\sigma_M)^T\) and represents the reduced dynamics and output as
$$
\dot \sigma_i = \lambda_i \sigma_i + I_i(\sigma)\,u(t), \qquad y(t)=y_0+G(\sigma).
$$
Here \(I_i(\sigma)\) and \(G(\sigma)\) are expanded in multivariate Taylor series in \(\sigma\), and truncating at total degree \(j\) yields a \(j\)-th-order accurate reduced model [2102.04526].

When the underlying equations are unknown, the coefficients can be inferred from steady-state responses to sinusoidal probing. With rank-one input
$$
u(t)=\epsilon \sin(\omega t),
$$
the steady-state output admits an \(\epsilon\)-expansion and Fourier series
$$
y_{ss}(\omega,t)-y_0
=
\sum_{m=1}^\infty \epsilon^m\Bigl[T_m(\omega)\sin(m\omega t)+S_m(\omega)\cos(m\omega t)\Bigr]
+O(\epsilon^{M+1}).
$$
The key structural fact is that the \(m\)-th harmonic first appears at \(O(\epsilon^m)\), and its amplitude is a linear combination of the \(m\)-th-order Taylor coefficients in \(G\) and \(I_i\). At first order,
$$
\frac{1}{\epsilon}\Gamma_1=\pi\,\Xi_1\,\Upsilon_1+O(\epsilon),
$$
with \(\Gamma_1\) collecting measured first harmonics, \(\Xi_1\) known analytically in terms of \(\lambda_i\), and \(\Upsilon_1\) containing the unknown coefficient products. Pseudoinverse solution gives
$$
\Upsilon_1=\frac{1}{\pi}\Xi_1^\dagger \Gamma_1.
$$
Second and higher orders have the same structure, with remainder terms depending only on lower-order quantities. The eigenvalues \(\lambda_i\) can also be refined by Newton iteration on the first-order relation [2102.04526].

The algorithmic summary consists of a preliminary data-driven estimate of \(\lambda_i\) from delay-embedded outputs and POD, symbolic assembly of response matrices \(\Xi_1,\dots,\Xi_j\), harmonic extraction at probe frequencies \(\omega_1,\dots,\omega_q\), pseudoinverse recovery of the coefficient vectors \(\Upsilon_m\), and optional Newton refinement before assembling the final reduced model [2102.04526].

The examples emphasize accuracy under finite forcing amplitudes. In a two-dimensional test model with \(\mu=-0.05\) and \(\lambda=-1\), one isostable suffices, the third-order reduction matches the full dynamics perfectly when \(f\) and \(g\) are known, and data-driven inference from noisy output recovers \(\lambda_1\approx -0.046\) together with coefficients up to cubic order. For a population of \(1000\) synaptically coupled Morris–Lecar-type neurons, two complex-conjugate isostables suffice, and the second-order inferred model outperforms the linear one by an order of magnitude in mean absolute error under composite transient inputs. For the one-dimensional Burgers’ equation projected onto five POD modes, a second-order three-isostable model reduces \(L^2\)-error by two orders of magnitude compared to the first-order model under composite multi-frequency boundary inputs [2102.04526].

A related input-driven reduction appears in the Goodwin oscillator example used for slow-manifold computation. After adding an external input \(u(t)\) to the \(B\)-equation and restricting to the manifold \(W^s\) with \(\psi_3=0\), the reduced model becomes
$$
\dot \psi_1 = \lambda_1\psi_1 + i(\psi_1)\,u(t),
$$
where \(i(\psi_1)=\nabla\psi_1\cdot[1,0,0]^T\) is tabulated along the computed manifold. In that example, the reduced one-complex-dimensional ODE reproduces steady-state periodic and period-doubling behavior far beyond the linear approximation; the full and isostable-reduced models both show a period-doubling bifurcation at input amplitude \(a\approx 0.02\)–\(0.023\), whereas a naive linearization fails to capture the bifurcation, and the maximum amplitude of \(B\) versus \(a\) matches to within a few percent [2507.13997].

## 5. Stochastic and phase–amplitude extensions

For stochastic oscillators governed by the Itô SDE
$$
dX_t=f(X_t)\,dt+g(X_t)\,dW_t,
$$
with diffusion tensor \(D(x)=\tfrac12 g(x)g(x)^T\), the relevant operator is the backward Kolmogorov, or stochastic Koopman, operator
$$
\mathcal{L}^\dagger \phi(x)=f(x)\cdot \nabla \phi(x)+D_{ij}(x)\,\partial^2_{x_i x_j}\phi(x).
$$
Under standard ellipticity and suitable boundary conditions, \(\mathcal{L}\) and \(\mathcal{L}^\dagger\) admit discrete spectra and biorthogonal eigenbases. A system is termed robustly oscillatory when the nonzero eigenvalue of \(\mathcal{L}^\dagger\) with maximal real part is a complex conjugate pair \(\lambda_\pm=\mu\pm i\omega\), with \(\mu<0\), \(\omega>0\), and all other eigenvalues satisfy \(\operatorname{Re}\lambda'\le 2\mu\). The asymptotic phase eigenfunction is then \(Q_+(x)\), and the asymptotic phase is
$$
\theta(x)=\operatorname{Arg}[Q_+(x)].
$$
If there is also a unique real eigenvalue \(\lambda_{\mathrm{iso}}=\alpha<0\) with maximal real part among real eigenvalues, the corresponding real eigenfunction \(Q_{\mathrm{iso}}(x)\) is interpreted as the stochastic isostable coordinate [2105.11048].

Normalization is fixed by biorthogonality, and \(Q_{\mathrm{iso}}\) is typically shifted so that its zero-level set defines an effective limit cycle
$$
\Sigma_0=\{x:Q_{\mathrm{iso}}(x)=0\}.
$$
In the phase–amplitude variables
$$
\theta=\operatorname{Arg}\psi_{\mathrm{phase}}(x), \qquad r=Q_{\mathrm{iso}}(x),
$$
Itô’s formula yields
$$
d\theta_t=[\omega+A_\theta(x)]\,dt+B_\theta(x)\,dW_t,
$$
$$
dr_t=\alpha r_t\,dt+B_r(x)\,dW_t.
$$
Thus, in the mean sense,
$$
E[d\theta_t]/dt=\omega, \qquad E[dr_t]/dt=\alpha\,E[r_t].
$$
The noisy linear focus provides a closed-form example: for the planar Ornstein–Uhlenbeck process with
$$
A=\begin{bmatrix}\mu&-\omega\\ \omega&\mu\end{bmatrix}, \qquad BB^T=2\epsilon I,
$$
one obtains
$$
Q_\pm(x)=x_1\pm i x_2,\qquad \lambda_\pm=\mu\pm i\omega,
$$
$$
Q_{\mathrm{iso}}(x)=2+(\mu/\epsilon)(x_1^2+x_2^2),\qquad \lambda_{\mathrm{iso}}=2\mu,
$$
so the effective limit cycle is the circle \(x_1^2+x_2^2=2\epsilon/|\mu|\) [2105.11048].

For deterministic limit cycles, the analogous construction uses Floquet theory rather than Koopman eigenfunctions at a fixed point. If the uncoupled ODE admits a stable \(T\)-periodic orbit \(\gamma=\{x^\gamma(t)\}\) with dominant nontrivial Floquet exponent \(\kappa<0\), the asymptotic phase \(\theta\) is defined by \(\dot \theta=\omega=2\pi/T\), and the slowest-decaying transverse coordinate is an isostable \(\Psi(x)\) satisfying
$$
\Psi(x)=\lim_{k\to\infty} e^{-\kappa t_k} w^T[\phi(t_k,x)-x^\gamma(0)], \qquad \dot \Psi=\kappa \Psi.
$$
Near the cycle, the phase–isostable normal form under weak perturbation \(P(t)\) is
$$
\dot \theta = \omega + Z^{(0)}(\theta)\cdot P(t) + O(\Psi),
$$
$$
\dot \Psi = \kappa \Psi + I^{(0)}(\theta)\cdot P(t) + O(\Psi),
$$
where \(Z^{(0)}\) is the iPRC and \(I^{(0)}\) is the iIRC [2310.02725]. The phase and isostable coordinates therefore separate tangential timing from the slowest transverse amplitude decay.

## 6. Network formulations, comparative scope, and limitations

For \(N\) identical coupled oscillators
$$
\dot x_i=F(x_i)+\epsilon\sum_j w_{ij}G(x_i,x_j),
$$
retaining the phase \(\theta_i\) and a single slow isostable \(\Psi_i\) at each node gives, after Taylor expansion and first-order averaging in \(\epsilon\),
$$
\dot \theta_i
=
\omega
+
\epsilon\sum_j w_{ij}\bigl[
H_1(\theta_j-\theta_i)
+\Psi_i H_2(\theta_j-\theta_i)
+\Psi_j H_3(\theta_j-\theta_i)
\bigr],
$$
$$
\dot \Psi_i
=
\kappa \Psi_i
+
\epsilon\sum_j w_{ij}\bigl[
H_4(\theta_j-\theta_i)
+\Psi_i H_5(\theta_j-\theta_i)
+\Psi_j H_6(\theta_j-\theta_i)
\bigr].
$$
The six \(2\pi\)-periodic coupling functions \(H_k\) are defined by averaging \(h_k\) over one period and depend on \(Z^{(0)},Z^{(1)},I^{(0)},I^{(1)}\), the coupling Jacobians, and the Floquet eigenfunction \(g^{(1)}\) [2310.02725].

The reduced network equations support explicit existence and stability conditions for phase-locked states. For a \(1{:}1\) phase-locked solution \(\theta_i=\phi_i+\Omega t\) with constant \(\Psi_i\), the locked amplitudes satisfy \(\Psi=Q^{-1}q\) and the collective frequency shift satisfies
$$
\Omega-\omega=\epsilon[\,p-PQ^{-1}q\,]\cdot 1_N.
$$
Linearization produces a \(2N\times 2N\) Jacobian in block form \(J=[H^1\ H^2; H^3\ H^4]\), and stability requires all eigenvalues except the rotational zero mode to have negative real part. In the globally coupled case, synchrony reduces to a \(2\times 2\) matrix \(M\), so stability requires
$$
\operatorname{Trace}M<0,\qquad \det M>0,\qquad \kappa+\epsilon[H_5(0)+H_6(0)]<0.
$$
Taking \(\Psi\to 0\) recovers the classical phase-only criterion \(H_1'(0)<0\) [2310.02725].

The mean-field complex Ginzburg–Landau equation provides an analytic benchmark. There, the phase-isostable reduction reproduces the synchrony boundary exactly, yields compact analytic expressions for the splay and antisynchrony boundaries, and tracks the true bifurcation loci qualitatively over a wide parameter range, whereas second- and third-order phase reductions agree only locally near \(\epsilon=0\). It also predicts bistability regions between synchrony and splay in agreement with the full system. In globally coupled Morris–Lecar networks, the \(2N\)-dimensional phase-isostable equations predict loss and restabilization of synchrony, a narrow window of stable antisynchrony, off-invariant-manifold phase-locked states, quasiperiodic solutions born in Hopf bifurcations, and cluster states for \(N=200\), while the first-order phase model predicts only trivial stable antisynchrony and unstable synchrony [2310.02725].

Across the fixed-point, forced, stochastic, and network settings, several limitations are explicit. Higher-order Taylor inference suffers from combinatorial growth in the number of coefficients and from noise sensitivity because fitting divides by \(\epsilon^j\); careful choice of \(\epsilon_m\) can mitigate this. The data-driven fixed-point framework requires a stable fixed point and forcing experiments that remain inside its basin of attraction. In slow-manifold computation, stiffness and fast-mode instability constrain backward marching, and asymptotic expansions require a nonresonance condition on slow indices. In network reduction, higher-order pure phase reductions become cumbersome for large \(N\), whereas the phase-isostable model remains compact by retaining one extra coordinate per node and six coupling functions \(H_k\) [2102.04526] [2507.13997] [2310.02725].

Taken together, these formulations define the isostable coordinate framework as a family of reductions in which the dominant decay structure of a nonlinear system is represented explicitly. At fixed points, the coordinates identify invariant slow manifolds by the condition that fast isostables vanish; under forcing, they support high-order reduced models and purely data-driven inference; in stochastic oscillators, they supply an amplitude variable dual to asymptotic phase; and in oscillator networks, they extend phase reduction by tracking slow transverse deviations from the attracting cycle [2507.13997] [2102.04526] [2105.11048] [2310.02725].

Source: https://www.emergentmind.com/topics/isostable-coordinate-framework