High: [Kinetics] Reversible-reaction analytical Jacobian ignores the equilibrium reverse term

[bernalde]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:kinetics.
• Tackle batch: P4 - high-impact numerical/model correctness.

The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.

---

Original audit details

Part of the model-correctness audit — see #17.

Severity: high | Category: kinetics | Location: PharmaPy/Kinetics.py lines 564-598, 603-609

What's wrong

get_rxn_rates computes net rates with the reverse term when keq_params is not None: f_terms = equilibrium_model(conc, temp, delta_hrxn) which returns f_term - r_term (lines 606-607, 529). However, when jac=True (line 596-598) it calls derivatives(conc, temp), and derivatives ALWAYS uses self.kinetic_model (= elem_f_model, forward only) and self.df_dstates (= elem_df_dstates, forward only) at lines 566-571. The reverse contribution -d/dC[ prod(C_prod^beta)/Keq ] is never added to the Jacobian.

Why it's incorrect

The net rate is r = kfprod(C_react^alpha) - (kf/Keq)prod(C_prod^beta). Its derivative wrt concentrations must include the reverse term: dr/dC_k = kfalpha_kf/C_k - (kf/Keq)beta_kg/C_k. The code returns only the first (forward) part, so the analytical state Jacobian disagrees with the actual rate model for any reversible (keq) reaction. Any implicit/Newton integrator or sensitivity analysis using this Jacobian gets an inconsistent (wrong) derivative, degrading convergence and giving incorrect sensitivities.

Suggested regression test

For a reversible reaction A<->B with finite Keq, compare derivatives(conc,temp) against a finite-difference Jacobian of get_rxn_rates(conc,temp,overall_rates=True). They will differ by the omitted reverse-term derivative whenever the reverse rate is nonzero.

---

Verification: 2/2 adversarial reviewers confirmed (refuted: 0).

[bernalde]

Independent Codex verification: CONFIRMED

severity_assessment: high
evidence: PharmaPy/Kinetics.py:505-529 computes reversible net rates as forward minus reverse, but PharmaPy/Kinetics.py:564-571 differentiates only self.kinetic_model, the forward term. For A <-> B with Keq=2, finite difference gives Jacobian columns [[-1, 0.5], [1, -0.5]], while the code misses the B column.
fix_note: Add the reverse-term derivative to derivatives when keq_params is set.

[parkyr]

Triage: reproduced ✅ (failing regression test attached)

Root cause is exactly as diagnosed in the audit and the independent Codex verification above: derivatives differentiates only self.df_dstates (→ elem_f_model, forward-only), so the reverse-term derivative is never added to the state Jacobian for keq_params-reactions. No restatement needed — confirmed against current base source.

Base commit: 1a3ab31 (origin/master)
Test command: python -m pytest tests/test_reversible_jacobian_reverse_term.py -q (pure NumPy; no assimulo marker needed)

Observed vs expected

Reversible A <-> B, Keq = 2, kf = 1, conc = [1, 1, 1, 1]. Analytical get_rxn_rates(..., jac=True) vs finite-difference Jacobian of the net rate — the (A, B) block:

entry	analytical (derivatives)	expected (finite diff / hand-calc)
∂ṙ_A/∂C_B	0.0	+0.5
∂ṙ_B/∂C_B	0.0	-0.5

The entire C_B column (the reverse-term contribution −β_k·(g/Keq)/C_k) is dropped.

Failing regression test

Place at tests/test_reversible_jacobian_reverse_term.py:

# -*- coding: utf-8 -*-
"""Regression test for issue #23.

For a reversible reaction the analytical state Jacobian returned by
``RxnKinetics.get_rxn_rates(..., jac=True)`` (``derivatives``) must equal the
finite-difference Jacobian of the net rate ``get_rxn_rates(..., jac=False)``.

``derivatives`` differentiates only the forward term (``self.df_dstates`` ->
``elem_df_dstates``), so for any reversible (``keq_params``) reaction the
reverse-term contribution ``-d/dC[prod(C_prod**beta)/Keq]`` is missing from the
Jacobian. This is pure NumPy (no solver backend), so no assimulo marker.
"""

import numpy as np

from PharmaPy.Kinetics import RxnKinetics


# A <-> B (reversible); C and solv are inert padding species in the database.
STOICH = [[-1, 1, 0, 0]]
PARTIC = ["A", "B", "C", "solv"]


def _db(data_path):
    return str(data_path["integration"] / "pfr_test_pure_comp.json")


def _reversible_kinetics(data_path):
    return RxnKinetics(
        path=_db(data_path),
        k_params=[1.0],        # kf = 1 at temp_ref = inf
        ea_params=[0.0],
        stoich_matrix=STOICH,
        partic_species=PARTIC,
        keq_params=[2.0],
        delta_hrxn=[0.0],      # Keq constant = 2, independent of temperature
        tref_hrxn=298.15,
    )


def test_reversible_jacobian_includes_reverse_term(data_path):
    kin = _reversible_kinetics(data_path)

    temp = 298.15
    conc = np.array([1.0, 1.0, 1.0, 1.0])   # reverse rate nonzero at C_B = 1
    deltah = np.array([0.0])

    analytical = np.asarray(
        kin.get_rxn_rates(conc, temp, overall_rates=True, jac=True)
    )

    # Central finite-difference Jacobian of the net species rates.
    def net_rates(c):
        return np.asarray(
            kin.get_rxn_rates(c, temp, overall_rates=True, delta_hrxn=deltah)
        )

    n = conc.size
    fd = np.zeros((n, n))
    h = 1e-6
    for j in range(n):
        cp = conc.copy(); cp[j] += h
        cm = conc.copy(); cm[j] -= h
        fd[:, j] = (net_rates(cp) - net_rates(cm)) / (2 * h)

    # Hand-computed truth for the participating (A, B) block: Keq = 2, kf = 1.
    # net = [-(C_A - C_B/2), +(C_A - C_B/2), 0, 0]
    #   d/dC_A = [-1, 1];  d/dC_B = [+0.5, -0.5]
    np.testing.assert_allclose(fd[0, 1], 0.5, atol=1e-5)
    np.testing.assert_allclose(fd[1, 1], -0.5, atol=1e-5)

    # The analytical Jacobian must reproduce the reverse-term B column.
    np.testing.assert_allclose(analytical, fd, atol=1e-5)

Red output against 1a3ab31:

tests/test_reversible_jacobian_reverse_term.py:73: in test_reversible_jacobian_includes_reverse_term
    np.testing.assert_allclose(analytical, fd, atol=1e-5)
E   AssertionError:
E   Not equal to tolerance rtol=1e-07, atol=1e-05
E   Mismatched elements: 2 / 16 (12.5%)
E   Max absolute difference: 0.5
E    x: array([[-1.,  0.,  0.,  0.],
E          [ 1.,  0.,  0.,  0.],
E          [ 0.,  0.,  0.,  0.],
E          [ 0.,  0.,  0.,  0.]])
E    y: array([[-1. ,  0.5,  0. ,  0. ],
E          [ 1. , -0.5,  0. ,  0. ],
E          [ 0. ,  0. ,  0. ,  0. ],
E          [ 0. ,  0. ,  0. ,  0. ]])

Test is not vacuously red: adding the reverse-term derivative to derivatives (df_dstates -= β_k·r_term/C_k when keq_params is set) turns it green.

Acceptance criteria

• derivatives (and get_rxn_rates(..., jac=True)) equals the finite-difference Jacobian of the net rate for reversible reactions, within FD tolerance.
• The forward-only Jacobian is unchanged when keq_params is None.
• Correct across a batch/time axis (2-D conc/temp), consistent with the vectorized equilibrium_model net-rate path.
• Existing kinetics tests still pass.

Next: gh-issue-to-pr 23.