Source-linked AI summary
Fast Surrogate Modeling of Excitable and Oscillatory FitzHugh-Nagumo Dynamics with Parametric Neural Operators
Andrew Franck, Justin Li
TL;DR
Classical solvers make repeated FHN parameter-space solves expensive, despite their relevance to connecting voltage recordings with model parameters. The paper trains parameter-conditioned Fourier Neural Operators for voltage and recovery fields, reporting low in-range error and substantial speedup over finite differences.
Problem
Classical finite-difference or spectral solvers require fine discretization and make the many fast, differentiable forward solves needed for parameter sweeps, inverse inference, and stimulus optimization expensive.
Source: §1 · p. 2Method
The paper trains a parameter-conditioned Fourier Neural Operator that uses feature-wise linear modulation to incorporate λ into every Fourier layer.
Source: Contributions · p. 2, Parameter Conditioning · p. 9Results
Sub-0.1% relative L2 error on both fields and nearly three orders of magnitude faster than the finite-difference baseline within the training parameter range.
Source: Contributions · p. 2Takeaways & Limitations
The surrogate is positioned as a fast, differentiable forward model for downstream inverse parameter inference and closed-loop stimulation, which the paper leaves for future work.
Source: Contributions · p. 2, Limitations and Future Work · p. 16All experiments are one-dimensional; extension to 2D/3D geometries is left for future work.
Source: Limitations and Future Work · p. 15
Abstract
from arXiv · showhide
The FitzHugh-Nagumo (FHN) system serves as a simplified model of neuronal voltage dynamics, capturing the activator-inhibitor structure behind both isolated action potentials and the rhythmic spiking seen across the brain. Exploring its 5D physiological parameter space is important for neuromodulation and mapping voltage recordings back to biophysics, yet classical finite-difference solvers make rapid parameter sweeps expensive. We train parameter-conditioned Fourier Neural Operators (FNOs) as fast, differentiable surrogates for the FHN voltage and recovery fields on a one-dimensional spatial domain, conditioning each Fourier layer on the parameter vector $\lambda = (D_u, D_v, a, b, \tau)$ via feature-wise linear modulation (FiLM). We apply a single bifurcation analysis that delimits the two distinct regimes the model spans, oscillatory (tonic firing) and excitable (action-potential propagation), and we train one operator in each. In the oscillatory regime the surrogate attains sub-$0.1\%$ relative $L^2$ error on both fields, runs nearly three orders of magnitude faster than the finite-difference baseline, generalizes uniformly across the parameter space, and extrapolates to low single-digit percentage errors outside of the training bounds. In the excitable regime the same operator accurately reproduces the firing threshold and the $c \propto \sqrt{D_u}$ conduction-velocity law and replicates full traveling pulses, fully capturing the excitable bifurcation structure rather than just smoothly interpolating fields.
1 Introduction
The paper targets fast, differentiable FHN forward solves for neural decoding, inverse modeling, and neuromodulation, where classical solvers make parameter sweeps expensive. It introduces a parameter-conditioned operator and evaluates it across oscillatory and excitable regimes.
Contributions:
- Feature-wise linear modulation injects λ = (Du, Dv, a, b, τ) into every Fourier layer, allowing one network to cover the 5D parameter space.The cited prior reaction-diffusion work fixed the physical parameters.Source: Contributions · p. 2
- A bifurcation analysis places the sampled parameter ranges in the oscillatory and excitable regimes of firing neurons.Source: Contributions · p. 2
- All parameter-error correlations satisfy |ρ| < 0.36, indicating relatively uniform generalization across the parameter space.The paper also reports extrapolation errors and identifies diffusion magnitude as the primary failure mode, while reaction and time-scale axes extrapolate accurately.Source: Contributions · p. 2
- Sub-0.1% relative L2 error on both fields and nearly three orders of magnitude faster than the finite-difference baseline in the training range.The paper positions the surrogate as a fast, differentiable forward model for downstream inverse-inference and closed-loop stimulation workloads, which it leaves for future work.Source: Contributions · p. 2
- In the excitable regime, the operator reproduces the all-or-none firing threshold, c ∝ √Du conduction-velocity scaling, and full traveling action potentials over 100-step rollouts.Source: Contributions · p. 2
2 Related Work and Mathematical Background
The paper frames FHN as a two-field reaction-diffusion model and situates its operator-learning approach against classical solvers and prior FNO applications. It learns a parameter-dependent map from initial fields to future fields using an FNO and a mean-squared-error training objective.
2.1 The FitzHugh-Nagumo System
- FHN describes coupled voltage and recovery dynamics, and different parameter combinations produce oscillatory or excitable regimes.Source: §2.1 · pp. 2–3
2.2 Learning Parametric Solution Operators with FNOs
- Classical numerical methods make large parameter sweeps expensive, while prior FNO work fixed physical parameters and evaluated generalization across initial conditions.This work instead considers stiff dynamics, sharper activations, and a five-dimensional parameter space governing distinct ODE bifurcations.Source: §2.1 · p. 3, §2.2 · p. 3
- The learned operator maps initial fields to fields one time step later, and repeated application provides longer-time predictions.Source: §2.2 · pp. 3–4
2.3 Fourier Neural Operator Architecture
- The FNO lifts input channels to a wider representation, applies Fourier layers, then projects back to the two output fields.Fourier layers perform convolutions in the frequency domain using learnable weights.Source: §2.3 · p. 4
2.4 Training Objective
- Training minimizes mean squared error over the dataset, with θ optimized using AdamW and learning-rate scheduling.Source: §2.4 · p. 4
3 Implementation
The implementation first verifies that the modeled parameter boxes occupy their intended dynamical regimes, then generates trajectories with a semi-implicit finite-difference solver.
3.1 Dynamical Regimes and Parameter-Box Placement
- The homogeneous-system fixed points are obtained from a cubic, and the trace and determinant of their Jacobian classify local flow as excitable or oscillatory; three real fixed points indicate bistability.Source: §3.1 · p. 5
- The oscillatory parameter cube lies entirely in the oscillatory regime, while the excitable box lies in the stable-rest region.A finer sweep over τ ∈ [1, 20] confirms that 100% of the oscillatory cube is oscillatory.Source: §3.1 · p. 5
3.2 Data Generation
- Training and validation trajectories are generated on a uniform periodic grid with second-order central differences and implicit treatment of the stiff diffusion term.The implicit operator's LU factorization is precomputed once per trajectory and reused at each step.Source: §3.2 · p. 5
Initial Conditions
Initial fields use Gaussian Random Fields, with Fourier coefficients transformed back to space and normalized. The parameter boxes are also placed within their respective dynamical regimes.
Initial Conditions
- Gaussian Random Fields are the standard initial-condition type used in one dimension.Source: Initial Conditions · p. 5
- Fourier coefficients use a standard normal draw with spectral decay controlled by α = 2.0, followed by inverse Fourier transform and field normalization.The initial v field is manually scaled by 0.5 to give the inhibitor a smaller amplitude than the activator.Source: Initial Conditions · p. 5
- The regime map distinguishes excitable, oscillatory, and bistable regions and places the two training boxes wholly within their respective regimes.Source: Figure 1 · p. 6
Parameter Sampling
The oscillatory training cube spans physiologically motivated values of five FHN parameters, while a separate dataset supports excitable-regime experiments.
Parameter Sampling
- The sampled ranges are Du ∈ U(0.01, 0.1), Dv ∈ U(0.005, 0.05), a ∈ U(−0.1, 0.1), b ∈ U(0.1, 0.5), and τ ∈ U(1.0, 20.0).The entire sampled parameter cube lies in the oscillatory regime of tonically firing neurons.Source: Parameter Sampling · p. 6
Dataset Specifications
- The primary dataset has 8,000 trajectories on a periodic domain with 256 points, evolved to T = 1.0 with 50 saved snapshots.A separate, disjoint held-out test set contains 1,600 trajectories.Source: Dataset Specifications · pp. 6–7
- Excitable-regime experiments use 1,250 trajectories on a longer domain, with localized super-threshold Gaussian-bump initial conditions.The dataset is split into 1,000 training and 250 validation trajectories and evolved to T = 40.Source: Dataset Specifications · p. 7
3.3 Model Architecture
The FNO processes the coupled voltage and recovery fields through Fourier and pointwise pathways, with lifting, projection, and a learnable residual connection framing the layers.
3.3 Model Architecture
- The architecture uses 16 Fourier modes, hidden width 64, six Fourier layers, and two input and output channels for the coupled fields.GELU is used for its smooth gradient properties, despite higher computational cost than ReLU.Source: §3.3 · p. 7
Lifting and Projection Networks
- Lifting expands the two-channel input to 128 channels, applies GELU, and projects to width 64; the projection network maps the hidden representation back to (u, v).The projection mirrors the lifting structure, using Conv1D layers with GELU between them.Source: Lifting and Projection Networks · p. 7
- A learnable global residual scale initialized at 0.1 lets the network learn small perturbations to the input state rather than reconstructing the output from scratch.Source: Lifting and Projection Networks · p. 7
3.3 Model Architecture
- Each Fourier layer sums a truncated spectral pathway and a pointwise Conv1D pathway before applying GELU.The spectral pathway retains the first 16 modes and zeroes the rest.Source: Fourier Layer Design · p. 7
Fourier Layer Design
The Fourier-mode ablation identifies kmax = 16 as an accuracy-and-size trade-off, with diffusion providing a rationale for retaining low-frequency modes.
Fourier Layer Design
- ∼4× validation rel-L2 improvement from kmax = 4 to kmax = 16; kmax = 64 is statistically indistinguishable from 16 but uses 3.4× more parameters.Figure 2 plots validation rel-L2 for u and v and parameter count across kmax settings.Source: Fourier Layer Design · p. 8, Figure 2 · p. 8
- kmax = 32 has a slightly better error score than 16, but kmax = 16 is more stable during autoregressive rollouts.Rollout stability is the stated ultimate use case for the surrogate.Source: Fourier Layer Design · p. 8
- Diffusion damps Fourier mode k at a rate proportional to k^2, so even sharp FHN spikes concentrate energy in the lowest few wavenumbers.Source: Fourier Layer Design · p. 8
- The ablation varies kmax while holding width 64 and depth six fixed, leaving a fuller width/depth sweep for future work.Source: Fourier Layer Design · p. 8
Parameter Conditioning
The model conditions a shared FNO on the FHN parameter vector so it can represent parameter-specific dynamics across the parameter space.
Parameter Conditioning
- Because the unconditioned lifting, Fourier, and projection blocks receive only (ut, vt), they cannot distinguish FHN systems with different parameters.Source: Parameter Conditioning · p. 8
- A three-layer MLP encodes λ, and independent per-layer heads use that encoding to produce channel-wise scales and shifts for Fourier-layer outputs.The encoder has widths 5 → 2w → w → w; modulation acts identically at every spatial location.Source: Parameter Conditioning · pp. 8–9
- Training the encoder and per-layer heads with the network turns the shared operator into a λ-indexed family of operators.Source: Parameter Conditioning · p. 9
3.4 Training and Evaluation
Training pairs are drawn from trajectories with independently normalized voltage and recovery fields, alongside each trajectory’s parameter vector. The section describes optimization, accuracy metrics, and repeated-seed evaluation.
3.4 Training and Evaluation
- Training pairs use independently z-normalized u and v fields, while the parameter vector λ is standardized before entering the FiLM encoder.Normalization statistics are computed from the training set.Source: §3.4 · p. 9
- Optimization uses per-element MSE over both fields with AdamW, learning rate 10−3, weight decay 10−4, batch size 32, and gradient clipping 1.0.Training runs for up to 1000 epochs on an A100 GPU; Table 1 covers data-generation, architecture, and training hyperparameters.Source: §3.4 · p. 9, Table 1 · p. 9
- Accuracy is assessed with per-field relative L2 error and MSE, plus autoregressive rollout error over up to 50 consecutive steps.Five-seed runs underlie Table 2 and the ±std values in Table 3; other tables use seed 42.Source: §3.4 · p. 9
4 Results
The results evaluate the oscillatory-regime surrogate on a disjoint held-out test set and assess single-step accuracy, parameter-wise error, autoregressive rollouts, extrapolation, and speed. A separately trained operator is tested in the excitable regime.
4.1 Single-Step Prediction Accuracy
- Relative L2 errors are 6.66 × 10−4 (0.067%) for u and 8.01 × 10−4 (0.080%) for v; FiLM also outperforms both reported baselines.Against channel-broadcast FNO, error is roughly 30% lower with 26% fewer parameters; DeepONet errors are roughly 58× higher on u and 50× on v.Source: §4.1 · pp. 10–12, Table 2 · p. 10, Table 3 · p. 12
4.2 Parameter Error Correlation
- Across the five parameter axes, the largest error correlation is |ρ| = 0.357 for τ–v; error dependence is moderate and concentrated in the recovery field.The activator’s strongest correlation is ρ = +0.143 with Du, while reaction parameters a and b have |ρ| ≤ 0.05 for both fields.Source: §4.2 · pp. 12–13, Table 4 · p. 12
4.3 Long-Duration Autoregressive Rollout
- After 50 autoregressive steps, relative L2 error is ∼0.5% for u and ∼0.7% for v, with limit-cycle geometry preserved and gradual phase drift.The longer-horizon error could grow to concerning levels; the authors suggest fixed-horizon rollout training as future work.Source: §4.3 · p. 13, Figure 4 · p. 13
4.4 Generalization to Unseen Parameters
- Extrapolation is weakest when diffusion exceeds training bounds: errors reach 3.4% for u and 1.6% for v, while reaction- and time-scale-axis cases remain accurate.The authors attribute diffusion-axis degradation to sharper spatial gradients being truncated by spectral layers.Source: §4.4 · p. 14, Table 5 · p. 14
4.5 Computational Efficiency
- The FNO is 111× faster at batch 1 and 920× faster at batch 32 than the finite-difference solver, with memory use of ∼150 MB versus 8 MB.The comparison is for one 50-step trajectory rollout.Source: §4.5 · p. 14, Table 6 · p. 15
4.6 Excitable Regime: Action-Potential Propagation
- In the excitable regime, the operator reproduces the sharp threshold at A⋆≈0.7 and c ∝ √Du, with 100-step relative L2 error below 0.9%.Across the full pulse, error is confined to a thin band at the moving wavefront.Source: §4.6 · p. 15
Limitations and Future Work
- The study is limited to one-dimensional experiments; diffusion-axis extrapolation degrades, and the linear parameter-sensitivity screen could miss nonlinear dependencies.Actual recovered-parameter or optimized-stimulus results are left for future work.Source: Limitations and Future Work · pp. 15–17
Ethical Statement
The paper flags the environmental cost of neural-operator training and warns against over-reliance on AI approximations and unrepresentative training data, particularly in biomedical use.
Ethical Statement
- Training neural operators is computationally intensive and carries a notable environmental footprint.The authors plan to limit training runs and use mixed precision and tuned batch sizes to improve hardware utilization.Source: Ethical Statement · p. 17
- Blind trust in FHN surrogate predictions could lead to incorrect features or research based on false data, so the models should augment rather than replace decision making.The authors emphasize this concern for high-stakes domains, especially medicine.Source: Ethical Statement · p. 17
- Unrepresentative simulated training scenarios could yield systematically biased outputs on underrepresented conditions, with potential for health disparities in biomedical applications.The authors recommend cross-checking the model for differing accuracy across patient groups.Source: Ethical Statement · p. 17