Forward sensitivity w.r.t. a switch time encoded as an SBML event (event-time sensitivity dt*/dp != 0)

[wshlavacek]

Summary

A switch time encoded as an SBML event is refused for forward sensitivity, while the same switch, same dynamics, same gradient encoded as piecewise(…, time >= T0, …) is supported once #48 lands. The refusal comes from NetworkModel::event_sensitivity_unsupported_reason (src/model.cpp), which rejects any event whose trigger references a requested sensitivity parameter:

event '…' has a trigger whose crossing time depends on the requested sensitivity parameter 'T0' (the event-time sensitivity dt*/dp is non-zero); event-time sensitivity is not yet supported (GH #212 Phase 2).

A survey of parity_checks/rr_parity/ shows this is not a corner case: 24 models / 88 events have an event whose crossing time is a fitted constant, and the parameter names are exactly the things one fits — treatment_start, ton/toff, tstim, Lockdown_start, tau_1…tau_f.

The encoding sensitivity is arbitrary from the modeller's side

Identical dynamics, identical gradient, opposite outcome:

encoding	forward X(6)	dX/dT0	sensitivity request
piecewise(kin, time >= T0, 0)	5.18799	−0.4060058	works (#48)
event time >= T0 → on := 1, rate kin*on	5.18799	−0.4060081 (FD)	refused

Both are natural ways to write "production starts at T0". Which one you picked — often decided by whichever tool exported the model — decides whether you can fit T0 by gradient.

Evidence: parity_checks/rr_parity

397 SBML files contain events; 236 unique models / 601 events after content-hash dedup (the corpus materializes each model under models/, smoke/, snapshots/ and data/sbml_overrides/, and BioModels COPASI exports reuse ids like model_0000001, so dedup by content, not id).

trigger	events	models	status
fixed-time, literal threshold	366	146	supported (GH #212 Phase 1)
state-dependent	138	54	out of scope
time vs. fitted constant parameter	88	24	refused — this issue
time vs. event-assigned counter	5	5	dose schedule, not a fixed crossing
other	4	2	—

Affected models: BIOMD0000000001, 117, 152, 153, 234, 244, 301, 327, 340, 342, 408, 422, 437, 494, 570, 650, 695, 820, 956, 960, 976, 983, MODEL1710030000, MODEL2310250001.

Threshold parameters seen: treatment_start/treatment_end (Owen1998 tumour treatment), ton/toff (Whitcomb2004), tstim, tGA4on/tGA4off, Lockdown_start/Lockdown_end (Ghanbari2020), tau_1…tau_f (Zongo2020 containment), single_pulse_duration, Treatment_Duration, heparin_infusion_duration_hr, Cycle_Int/N_Cycle.

Two caveats on the count, so it isn't over-read:

• BIOMD152/153 alone contribute 42 of the 88 events (a cAMP spike train), so 24 models is the honest breadth figure, not 88.
• This counts what is encodable, not what anyone is currently fitting.

Notably, Zongo2020 and Ghanbari2020 are the same genre as the #48 exemplar (Lin2021) — COVID containment onset times — just encoded as events rather than if().

Sub-finding: 25 of those events are additionally blocked by a vacuous flag

63 events / 15 models are refused only for the event-time reason. The remaining 25 events / 9 models also trip the delay or persistence checks — but in every one of those cases the delay is a literal 0, and per SBML L3v2 §4.11.3 non-persistence can only matter if the trigger can revert between trigger time and execution time, which a zero-delay event has no window for. So both flags are no-ops there.

Normalizing delay = 0 to "no delay" (and treating persistent=false as vacuous when there is no delay) is worth doing independently of this issue — it would also widen the existing Phase-1 support. Ghanbari2020 and Zongo2020 are blocked by exactly this.

Proposed change

The general forward-sensitivity jump across an event at a parameter-dependent time t*(p), with state reset x⁺ = h(x⁻, p):

s⁺ = ∂h/∂x·(s⁻ + f⁻·∂t*/∂p) + ∂h/∂p − f⁺·∂t*/∂p

Read it as: shift s⁻ along the pre-event flow by how far the event time moves, apply the event Jacobian, shift back along the post-event flow. Both corners are already implemented:

• ∂t*/∂p = 0 (fixed time, state jumps) → s⁺ = ∂h/∂x·s⁻ + ∂h/∂p — the existing apply_event_sensitivity_jump (GH #212).
• h = identity (state continuous, time moves) → s⁺ = s⁻ + (f⁻−f⁺)·∂t*/∂p — the Gradient-based fitting of switch-time parameters: forward sensitivity drops the discontinuity-crossing term at if(t>=σ) steps #48 switch jump.

So this is not new math, it is the two existing halves plus the cross terms. For the common "flag" pattern (on := 1, rate laws read on) it decomposes row-wise: the assigned rows get the GH #212 jump, every other row gets the #48 jump.

Why this is cheaper than #48 was

Neither of #48's two hard parts applies:

• No CVodeSetStopTime. Gradient-based fitting of switch-time parameters: forward sensitivity drops the discontinuity-crossing term at if(t>=σ) steps #48's stall came from a kink in f that CVODES hit while approaching it. An event trigger is not in the RHS — f is smooth right up to the root — and CVODE's root-finding already stops exactly at t*.
• No parameter pinning. ∂f/∂T0 is genuinely zero without help, since T0 never appears in f. Confirmed: with an event root live and sensitivities on for the other parameters, a run completes normally and matches the closed form to 9 digits.

Work items:

1. Compute ∂t*/∂p for an event trigger — the same threshold differentiation bngsim/_switch_sensitivity.py (Gradient-based fitting of switch-time parameters: forward sensitivity drops the discontinuity-crossing term at if(t>=σ) steps #48) already does, pointed at ev.trigger_expr_idx instead of an if() condition.
2. Add the pre-shift f⁻·∂t*/∂p and post-shift −f⁺·∂t*/∂p terms to apply_event_sensitivity_jump in src/cvode_simulator.cpp. ∂h/∂x and ∂h/∂p are already computed there by FD on the assignment expressions.
3. Normalize vacuous delay = 0 / non-persistence (see sub-finding), then narrow event_sensitivity_unsupported_reason to what genuinely remains: state-dependent triggers, real delays, non-persistence with a real delay.

Acceptance criteria

• An event-encoded onset (time >= T0 → on := 1, rate kin*on) returns dX/dT0 matching central finite differences and the closed form, and matching the piecewise encoding of the same dynamics to integration tolerance.
• A parameter-valued reset on a parameter-dependent trigger (both ∂h/∂p ≠ 0 and ∂t*/∂p ≠ 0) matches FD — this is the term neither existing corner exercises.
• The answer does not drift with rtol (the failure mode Gradient-based fitting of switch-time parameters: forward sensitivity drops the discontinuity-crossing term at if(t>=σ) steps #48 hit: a jump read at a perturbed parameter point scales by 1 ∓ √rtol).
• Zongo2020 / Ghanbari2020 / Owen1998 (BIOMD650) run with their onset parameters in sensitivity_params.
• Still refused, with a clear reason: state-dependent triggers, real delays, non-persistence with a real delay.

Relationship to #48

#48 covers the if() / piecewise encoding (BNGL counter clock and SBML literal time). This issue covers the same modelling intent expressed as an event. They share the threshold-differentiation code and the jump concept; #48 should land first.