Minimum-Energy FitzHugh-Nagumo Neural Stimulation¶
Background¶
The FitzHugh-Nagumo model compresses excitable-cell dynamics into a fast membrane-voltage variable and a slow recovery variable. It is simple enough for an optimal-control tutorial while retaining the threshold-like transition that makes neural excitation nonlinear.
This example finds a bounded, one-sided current that takes a resting model neuron to a specified voltage at a fixed deadline with minimum normalized stimulus energy. The solution is accepted only if an independent forward integration reaches that voltage while still moving upward. Reaching the target is an upstroke crossing, not a terminal equilibrium or proof of a complete post-deadline action potential.
Problem formulation¶
Let \(v(t)\) be membrane voltage, \(w(t)\) the recovery variable, and \(I(t)\) the injected current. In nondimensional form,
with
The initial condition is the unforced resting equilibrium used by the script:
At the fixed nondimensional horizon \(T=12\), only voltage is prescribed:
The objective is normalized squared current,
Path bounds keep the trajectory in the modeled region:
The optimization does not impose a terminal derivative constraint. After the solve, the independently integrated trajectory must satisfy \(\dot v(T)>0\); otherwise the script rejects it as the wrong crossing direction.
Variables and units¶
| Symbol | Meaning | Unit or value |
|---|---|---|
| \(t\) | Model time | nondimensional |
| \(v\) | Membrane-voltage state | nondimensional |
| \(w\) | Slow recovery state | nondimensional |
| \(I\) | Injected stimulus current | nondimensional |
| \(a\) | Recovery offset | \(0.7\) |
| \(b\) | Recovery gain | \(0.8\) |
| \(\varepsilon\) | Fast/slow time-scale ratio | \(0.08\) |
| \(T\) | Fixed stimulation horizon | \(12\) |
| \(J\) | Current normalized by \(I_{\max}\), squared and integrated | nondimensional |
Every quantity is nondimensional. The plotted current and reported charge or energy therefore cannot be converted to amperes, coulombs, or joules without a separate physical scaling and electrode model.
Modeling choices¶
The nonnegative bound makes this an idealized monophasic stimulus. It is useful for studying nonlinear excitation but does not enforce charge balance. The recovery state is free at the deadline because the task is to reach the voltage on its upstroke, not to settle at a chosen equilibrium.
The script uses an \(80\times4\) Lobatto mesh. A smooth depolarizing state history and moderate current initialize the NLP; they are only a starting guess and do not replace the dynamic constraints.
Independent validation¶
After collocation, the current is reconstructed at 4,001 times. A separate high-accuracy DOP853 integration starts from \((v_r,w_r)\) and uses that sampled current. The script compares the integrated states with the collocation states, checks the forward terminal error, and evaluates
on the independently propagated trajectory. It also checks all dense state and current bounds.
Verified result¶
| Quantity | Verified value |
|---|---|
| Normalized objective \(J\) | \(0.03586644\) |
| Injected charge \(\int I\,\mathrm dt\) | \(0.93127477\) |
| Stimulus energy \(\int I^2\,\mathrm dt\) | \(0.14346575\) |
| Peak current | \(0.21801146\) |
| Forward terminal-voltage error | \(2.272\times10^{-5}\) |
| Forward terminal voltage rate \(\dot v(T)\) | \(0.721925\) |
| Maximum forward-integration error | \(2.289\times10^{-5}\) |
| Maximum dense path-bound violation | \(0\) |
The positive terminal voltage rate verifies an upstroke arrival. Also, \(\int I^2\,\mathrm dt=I_{\max}^2J=4J\) to the reported precision, providing a direct consistency check on the normalized objective.
Scope and safety limitations¶
This is a teaching model, not a clinical waveform or medical recommendation. It omits electrode and tissue dynamics, charge-balancing phases, electrode polarization, spatial activation, safety limits in physical units, parameter variability, noise, and the trajectory after the deadline. A physically deployable stimulation design requires a calibrated biophysical model, charge-balanced hardware constraints, and domain-specific safety review.
Run the example¶
Save the figure without opening a window:
system, phase = build_problem()
guess = initial_guess(phase)
solution = solve_problem(system, guess)
plot_solution(solution)
Source code¶
See the complete runnable example:
examples/neural_stimulation.py.