HeunDeterministic: predictor-step bounds clamping degrades integration from O(dt²) to O(dt) near state variable boundaries

[Ziaeemehr]

Summary

HeunDeterministic.scheme() applies integration_bound_and_clamp() to the predictor
intermediate state before evaluating k2. This causes Heun to silently degrade
from 2nd-order to 1st-order accuracy at every integration step where a state variable
boundary is active.

Location

tvb/simulator/integrators.py — HeunDeterministic.scheme():

m_dx_tn = dfun(X, coupling, local_coupling)
inter = X + self.dt * (m_dx_tn + stimulus)
self.integration_bound_and_clamp(inter)           # ← problem here

dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0

X_next = X + dX + self.dt * stimulus
self.integration_bound_and_clamp(X_next)

Why this is incorrect

Heun achieves 2nd-order accuracy by using two slope estimates:

k1 = f(Xₙ)
k2 = f(Xₙ + dt·k1)          ← must be the genuine Euler step
X_{n+1} = Xₙ + dt/2·(k1+k2)

When inter is clamped before evaluating k2, k2 is no longer evaluated at the
genuine predictor point Xₙ + dt·f(Xₙ) but at a boundary-projected point. The average
(k1+k2)/2 is no longer a consistent 2nd-order approximation of the true average slope.
The global truncation error at every step where clamping activates becomes O(dt)
instead of O(dt²) — effectively 1st-order (Euler) accuracy at those steps.

Affected models

Any model with state_variable_boundaries where the state can transiently violate a
bound during the predictor sub-step:

• GastSchmidtKnosche_SD / _SF — r ≥ 0
• CoombesByrne2D — r ≥ 0
• MontbrioPazoRoxin — r ≥ 0
• WilsonCowan — E, I ∈ [0, 1]
• ReducedWongWang — S ∈ [0, 1]

For small dt or mild dynamics the predictor rarely violates bounds and the effect is
invisible. For larger dt or aggressive parameter regimes the 1st-order error
accumulates silently.

Impact on reproducibility

Third-party simulators that implement standard Heun (post-corrector clamp only) produce
trajectories that diverge from TVB even with identical model equations, parameters,
connectivity, and random seeds. This makes cross-validation against TVB and reproduction
of TVB results with other tools difficult.

We encountered this while cross-validating a Numba/JAX reimplementation against TVB:
trajectories matched exactly for all models where bounds were never activated during
integration, but diverged for GastSchmidtKnosche_SD/SF where r transiently goes
below 0 during the predictor step.

Suggested fix

Apply integration_bound_and_clamp only after the corrector — not to the
predictor. This preserves 2nd-order accuracy while guaranteeing the returned state
satisfies all bounds:

def scheme(self, X, dfun, coupling, local_coupling, stimulus):
    m_dx_tn = dfun(X, coupling, local_coupling)
    inter = X + self.dt * (m_dx_tn + stimulus)
    # No clamping here — k2 evaluated at the genuine predictor point

    dX = (m_dx_tn + dfun(inter, coupling, local_coupling)) * self.dt / 2.0

    X_next = X + dX + self.dt * stimulus
    self.integration_bound_and_clamp(X_next)   # clamp once, post-corrector
    return X_next

For QIF-family models the ODE itself provides a natural restoring force at r = 0:
dr/dt|_{r=0} = Δ/(πτ) > 0 for Δ > 0, so small transient predictor violations are
self-correcting and post-corrector clamping is sufficient for both safety and accuracy.

Note on HeunStochastic

The same pattern appears in HeunStochastic.scheme() and should be addressed
consistently.

References

• Hairer, Nørsett, Wanner. Solving Ordinary Differential Equations I, §II.1
  (Runge-Kutta methods, order conditions)
• Montbrio, Pazo, Roxin. Macroscopic Description for Networks of Spiking Neurons.
  Physical Review X, 2015.
• Gast, Schmidt, Knösche. A mean-field description of bursting dynamics in spiking
  neural networks with short-term adaptation. Neural Computation, 2020.