feat(sensitivity): carry the equilibration's dx/dθ into a parameter scan (#81)

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.

[wshlavacek]

Closes#81.

The gap

A pre-equilibrate → scored dose scan could not be gradient-fit: parameter_scan /
bifurcate refused every sensitivity-configured Simulator. As the issue says, the
refusal is right — each point's initial condition is the equilibrated snapshot x_ss(θ),
so re-seeding it as if it were the model's seed ICs discards the dx/dθ accumulated
during the equilibration. The gap was that there was no correct option instead.

Answering the issue's two questions:

1. Does carry_sensitivities already give the right semantics across the boundary?
   Yes — for a sequential two-phase run. What was missing is that every primitive which
   restores a state dropped the derivative on the way: set_state() /
   restore_concentrations() clear the pending seed, and save_concentrations() cleared it
   and the "carried-over" flag.
2. Does a scan need a per-point clone that inherits both x and dx/dp? Not a clone —
   it needs the restore path to carry the derivative with the state. Then the existing
   per-point loop just passes carry_sensitivities=True.

What changed

• Model.save_concentrations() redefines the IC baseline to the current state; the state
  did not change, so neither did its dx/dθ — the new baseline inherits it, and
  reset() restores both. A baseline saved with nothing carried is θ-independent literal
  ICs, i.e. the pre-No way to carry forward sensitivities through a pre-equilibration into a scored parameter scan (blocks gradient fitting of pre-equilibrated dose-response) #81 fresh start, so reset() on an ordinary model (and the
  every-action-reset backends) is unchanged.
• save_concentrations(label=…) / restore_concentrations(label) capture and restore a
  named snapshot's dx/dθ the same way.
• Simulator.parameter_scan / bifurcate restore the reset target's state and its
  dx/dθ per point, integrate each point with carry_sensitivities=True, and leave the
  model — state, scanned parameter, carried derivative, dirty flag — exactly as found. A
  continuation scan (reset_conc=False) chains each point's dx/dθ from the previous
  point, so the whole sweep is one differentiable protocol.
• NetworkModel gained the write half of the seed accessor,
  set_pending_sensitivity_seed(seed, param_names) — the primitive that lets any protocol
  restore a state together with its θ-derivative (PyBNF drives its own per-point loop, so
  it can use this directly without parameter_scan) — plus
  has_baseline_sensitivity_seed and a writable ic_state_dirty for the restore.

sim=bngsim.Simulator(model, method="ode", sensitivity_params=["kf", "kr", "kdeg"])
sim.run(t_span=(0, 1e6), n_points=2, steady_state=True) # equilibrate ONCEpoints=sim.parameter_scan("L_0", par_min=1e-3, par_max=1e2, n_scan_pts=12,
log_scale=True, t_span=(0, 600), n_points=61,
steady_state=True)
grads= [p.output_sensitivities("observable:pReceptor") forpinpoints]

Measurements

Oracle is a closed form, not a finite difference: preequil_prod_deg.net is
dA/dt = k_prod − (k_deg + dose)·A equilibrated at dose = 0, so every point's
dA/dθ is exact. Four doses off one equilibration:

dose	carried (this PR)	re-seeded per point, t=0	re-seeded, t=3 (k_prod / k_deg)
0.5	1.2e-10	100%	9.5% / 15.3%
1.0	1.7e-10	100%	3.3% / 8.4%
2.0	2.3e-10	100%	0.3% / 1.3%
4.0	1.2e-09	100%	0.0% / 0.0%

The trajectories are identical to 1e-7 in every row — the seed is the only difference. The
same scan scored at each point's own steady state matches dA_ss/dθ = (1/λ, −k_prod/λ²)
to 1e-12. Note the shape of the wrong column: re-seeding is 100% off at t=0 and its
error decays with dose, so it is worst at the low-dose end of a dose-response — the
informative part of the fit — and looks fine at saturating dose.

On the issue's own model family (IGF1R_model_v1, 589 species / 4198 reactions,
equilibrated at a basal ligand dose so 1107 of 1178 carried seed entries are live, then
scanned over three doses applying each with on_point + set_concentration):
d(pY980)/d{kp, kdp} matches central FD of the measured observable over the full
protocol at 1e-8 … 1e-9, stable across h/p = 1e-3, 1e-4, 1e-5 (signal, not FD
noise), and the whole scan runs off ONE equilibration in 2.0 s. The scan is also
bit-identical to the hand-rolled equivalent (restore snapshot → install seed →
run(carry_sensitivities=True)) on every point.

A silent defect this also closes

No scan required: equilibrate → save_concentrations() → run(sensitivities) used to drop
the carry and clear the carried-over flag, so it re-seeded a θ-dependent initial
condition as a fresh start and returned wrong derivatives with no warning (100% off at
t=0 on the model above; reset() after the save had the same hole). That is exactly the
BNG action ordering simulate({steady_state=>1}) → saveConcentrations() → parameter_scan, so a fix that only worked from the live state would not have unblocked
the protocol shape in the issue.

Refusals kept (never a silently re-seeded gradient)

Situation	Why it cannot be answered
the reset target carries no matching dx/dθ	nothing to seed from (no equilibration on this Simulator with these params, or a plain run / set_concentration() / set_state() dropped it)
the scanned parameter is a sensitivity_params entry	each point overwrites it, so the derivative carried into the point was taken at a different value of the same symbol
sensitivity_ic across the boundary	the point starts from a snapshot, not the model's ICs, so ∂y/∂y_k(0) has no meaning
an on_point hook moves a differentiated parameter	same composition problem as scanning one

An on_point hook may apply the usual coupled setConcentration dose override:
assigning species k a literal value makes ∂x_k(0)/∂θ = 0 for that species (what a
literal assignment means, and what set_concentration already documents) while the rest
keep the carried derivative. The one case bngsim cannot infer is an override computed
from a differentiated parameter — such a hook installs the correct rows itself with
set_pending_sensitivity_seed after its concentration writes, and a seed still pending
when the hook returns is used verbatim.

Tests / verification

• python/tests/test_preequilibration_sensitivity.py — new TestScanCarryOver (closed-form
  per point, per-point t=0 seed, observable output sensitivities, scored steady-state
  dose-response, the BNG save_concentrations ordering, reset_to named snapshot,
  continuation == a chain of sequential carry runs, on_point literal-dose zero row and the
  install-your-own-seed escape hatch, "left as we found it", and the five refusals), plus
  baseline-inheritance and seed-write round-trip/validation cases.
• python/tests/test_parameter_scan.py — the old blanket-refusal test split into the two
  refusals that remain.
• python/tests/test_model_clone.py — the Impl clone contract for the pending and
  baseline seeds (new fields ⇒ new cases, per the contract comment).
• Full suite: 2950 passed, 67 skipped (python/tests, in-tree); C++ ctest5/5;
  pre-commit (ruff, ruff-format, clang-format, mypy) clean.
• No codegen / RHS / Jacobian code is touched, so no corpus sweep applies: outside a
  sensitivity workflow the only behavior change is that a plain scan now rewinds the
  carried-over flag it used to leave set.