diff --git a/algebra/derive_four_state_r_equivalent_coefficients.ipynb b/algebra/derive_four_state_r_equivalent_coefficients.ipynb index cb0d4eb..88810e0 100644 --- a/algebra/derive_four_state_r_equivalent_coefficients.ipynb +++ b/algebra/derive_four_state_r_equivalent_coefficients.ipynb @@ -5,365 +5,506 @@ "id": "overview", "metadata": {}, "source": [ - "# How `bindcurve` uses the four-state models\n", + "# Four-state incomplete-competition model in true free receptor\n", "\n", - "This notebook is intentionally implementation-oriented. Instead of re-deriving the full Roehrl algebra symbolically, it shows the runtime workflow used by `bindcurve` for the incomplete-competition models `comp_4st_specific` and `comp_4st_total`.\n", + "This notebook derives and validates the quintic used by `bindcurve` for the four-state competitive-binding models.\n", "\n", - "**Reference**\n", + "The filename still contains `r_equivalent` for historical reasons. The implementation no longer solves the four-state model in the artificial coordinate\n", "\n", - "Roehrl, M. H. A.; Wang, J. Y.; Wagner, G. *Biochemistry* **2004**, 43(51), 16056-16066. DOI: [10.1021/bi048233g](https://doi.org/10.1021/bi048233g)\n", + "$$\n", + "r_{\\mathrm{equiv}} = K_d^*\\frac{F_b^*}{1-F_b^*}.\n", + "$$\n", + "\n", + "Instead, the runtime model solves the quintic directly in literal free receptor concentration\n", + "\n", + "$$\n", + "R \\equiv [R].\n", + "$$\n", + "\n", + "The response is then computed from the actual four-state bound tracer population, `RS + RLS`.\n", + "\n", + "Reference:\n", + "\n", + "Roehrl, M. H. A.; Wang, J. Y.; Wagner, G. *Biochemistry* **2004**, 43(51), 16056-16066. DOI: [10.1021/bi048233g](https://doi.org/10.1021/bi048233g)\n" + ] + }, + { + "cell_type": "markdown", + "id": "species-model", + "metadata": {}, + "source": [ + "## Species model\n", + "\n", + "The four-state incomplete-competition model uses free species\n", + "\n", + "$$\n", + "R,\\quad S \\equiv L^*,\\quad L\n", + "$$\n", + "\n", + "and complexes\n", + "\n", + "$$\n", + "RS,\\quad RL,\\quad RLS.\n", + "$$\n", "\n", - "The key detail is that `bindcurve` does **not** solve the four-state model directly in literal free receptor concentration. It solves a quintic in the transformed coordinate\n", + "The equilibrium concentration expressions are\n", "\n", "$$\n", - "r_{\\mathrm{equiv}} = K_{ds}\\frac{FSB}{1-FSB},\n", + "[RS] = \\frac{[R][S]}{K_d^*},\n", "$$\n", "\n", - "and then converts the selected physical root back to tracer-bound fraction with\n", + "$$\n", + "[RL] = \\frac{[R][L]}{K_d},\n", + "$$\n", + "\n", + "$$\n", + "[RLS] = \\frac{[R][L][S]}{K_dK_{d3}}.\n", + "$$\n", + "\n", + "Here `Kds` in the source code means $K_d^*$, the tracer dissociation constant.\n", + "\n", + "The mass balances are\n", + "\n", + "$$\n", + "R_T = R + RS + RL + RLS,\n", + "$$\n", "\n", "$$\n", - "FSB = \\frac{r_{\\mathrm{equiv}}}{K_{ds} + r_{\\mathrm{equiv}}}.\n", + "S_T = S + RS + RLS,\n", "$$\n", "\n", - "That recovered `FSB` is what gets mapped to the experimental signal.\n" + "$$\n", + "L_T = L + RL + RLS.\n", + "$$\n", + "\n", + "The experimental tracer-bound fraction is\n", + "\n", + "$$\n", + "F_b^* = \\frac{RS + RLS}{S_T}.\n", + "$$\n" ] }, { "cell_type": "markdown", - "id": "when-to-use", + "id": "symbolic-intro", "metadata": {}, "source": [ - "## When the four-state models are used\n", + "## Symbolic derivation of the quintic in free receptor\n", + "\n", + "The following cell eliminates free tracer `S` and free competitor `L` from the three mass-balance equations, leaving one polynomial in true free receptor `R`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "symbolic-derive", + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sp\n", + "\n", + "R, S, L = sp.symbols(\"R S L\")\n", + "RT, ST, LT = sp.symbols(\"RT ST LT\")\n", + "Kds, Kd, Kd3 = sp.symbols(\"Kds Kd Kd3\")\n", "\n", - "The four-state models are for **incomplete competition**. In addition to `RLs` and `RL`, the system allows the ternary complex `RLLs`, so the tracer can still bind after the competitor is already present.\n", + "RS = R * S / Kds\n", + "RL = R * L / Kd\n", + "RLS = R * L * S / (Kd * Kd3)\n", "\n", - "- `Kds`: tracer affinity for free receptor\n", - "- `Kd`: competitor affinity for free receptor\n", - "- `Kd3`: tracer affinity for the receptor-competitor complex\n", - "- `comp_4st_specific`: solves the specific-binding model\n", - "- `comp_4st_total`: uses the same machinery, but replaces `Kd` by `(1 + N) * Kd`\n", + "eq_receptor = RT - R - RS - RL - RLS\n", + "eq_tracer = ST - S - RS - RLS\n", + "eq_competitor = LT - L - RL - RLS\n", "\n", - "The implementation lives in `bindcurve/modeling/competitive.py`.\n" + "common_denominator = Kds * Kd * Kd3\n", + "p_receptor = sp.expand(eq_receptor * common_denominator)\n", + "p_tracer = sp.expand(eq_tracer * common_denominator)\n", + "p_competitor = sp.expand(eq_competitor * common_denominator)\n", + "\n", + "relation_ligand = sp.resultant(p_tracer, p_competitor, S)\n", + "relation_receptor = sp.resultant(p_receptor, p_tracer, S)\n", + "\n", + "relation_ligand = sp.factor(relation_ligand / (Kd3 * Kds))\n", + "relation_receptor = sp.factor(-relation_receptor / (Kd3 * Kds))\n", + "\n", + "resultant = sp.resultant(relation_ligand, relation_receptor, L)\n", + "factored = sp.factor(resultant)\n", + "\n", + "factored\n" + ] + }, + { + "cell_type": "markdown", + "id": "extract-intro", + "metadata": {}, + "source": [ + "The resultant contains extraneous factors from the elimination. The physical receptor polynomial is the quintic factor.\n" ] }, { "cell_type": "code", - "execution_count": 4, - "id": "imports", + "execution_count": null, + "id": "extract-quintic", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Example parameter set used below:\n", - "{'Kd': 1.6,\n", - " 'Kd3': 0.5,\n", - " 'Kds': 0.02,\n", - " 'LsT': 0.005,\n", - " 'N': 0.35,\n", - " 'RT': 0.05,\n", - " 'ymax': 95.0,\n", - " 'ymin': 5.0}\n", - "\n", - "Effective Kd for comp_4st_total: 2.160000 uM\n" - ] - } - ], + "outputs": [], "source": [ - "from pprint import pprint\n", + "extraneous = -Kd3 * Kd**4 * Kds**2 * R**2 * ST\n", + "quintic = sp.factor(factored / extraneous)\n", + "poly = sp.Poly(quintic, R)\n", "\n", - "import numpy as np\n", - "import pandas as pd\n", + "assert poly.degree() == 5\n", "\n", - "from bindcurve.modeling.competitive import (\n", - " CompetitiveFourStateSpecificKdModel,\n", - " CompetitiveFourStateTotalKdModel,\n", - " _competitive_four_state_coefficients,\n", - " _direct_bound_fraction,\n", - " _equivalent_receptor_bounds,\n", - " _equivalent_receptor_from_bound_fraction,\n", - " _select_physical_root,\n", - ")\n", + "coefficients = poly.all_coeffs()\n", + "coefficients\n" + ] + }, + { + "cell_type": "markdown", + "id": "coefficients-intro", + "metadata": {}, + "source": [ + "Written as\n", "\n", - "np.set_printoptions(precision=6, suppress=True)\n", + "$$\n", + "aR^5 + bR^4 + cR^3 + dR^2 + eR + f = 0,\n", + "$$\n", "\n", - "params = {\n", - " \"RT\": 0.05,\n", - " \"LsT\": 0.005,\n", - " \"Kds\": 0.02,\n", - " \"Kd\": 1.6,\n", - " \"Kd3\": 0.5,\n", - " \"N\": 0.35,\n", - " \"ymin\": 5.0,\n", - " \"ymax\": 95.0,\n", - "}\n", + "the coefficients are the expressions below.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "print-coefficients", + "metadata": {}, + "outputs": [], + "source": [ + "for symbol, coefficient in zip(\"abcdef\", coefficients, strict=True):\n", + " print(f\"{symbol} =\")\n", + " print(sp.factor(coefficient))\n", + " print()\n" + ] + }, + { + "cell_type": "markdown", + "id": "runtime-helper-intro", + "metadata": {}, + "source": [ + "## Runtime coefficient helper\n", "\n", - "print(\"Example parameter set used below:\")\n", - "pprint(params)\n", - "print(f\"\\nEffective Kd for comp_4st_total: {(1 + params['N']) * params['Kd']:.6f} uM\")\n" + "The implementation in `bindcurve.modeling.competitive._competitive_four_state_coefficients` should match the symbolic coefficients above. This cell verifies that correspondence symbolically.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "runtime-symbolic-check", + "metadata": {}, + "outputs": [], + "source": [ + "runtime_coefficients = [\n", + " Kds - Kd3,\n", + " (\n", + " -Kd3 * Kd\n", + " - Kd3 * Kds\n", + " - Kd3 * LT\n", + " + Kd3 * RT\n", + " - Kd3 * ST\n", + " + Kd * Kds\n", + " + Kds**2\n", + " + Kds * LT\n", + " - 2 * Kds * RT\n", + " + Kds * ST\n", + " ),\n", + " (\n", + " Kd3 * Kd * RT\n", + " - Kd3 * Kd * ST\n", + " - Kd3 * Kds * LT\n", + " + Kd3 * Kds * RT\n", + " + Kd * Kds**2\n", + " - 2 * Kd * Kds * RT\n", + " + 2 * Kd * Kds * ST\n", + " + 2 * Kds**2 * LT\n", + " - 2 * Kds**2 * RT\n", + " - Kds * LT * RT\n", + " + Kds * LT * ST\n", + " + Kds * RT**2\n", + " - Kds * RT * ST\n", + " ),\n", + " (\n", + " Kd3 * Kd**2 * Kds\n", + " + Kd3 * Kd * Kds**2\n", + " + Kd3 * Kd * Kds * LT\n", + " + Kd3 * Kd * Kds * ST\n", + " + Kd * Kds**2 * LT\n", + " - 2 * Kd * Kds**2 * RT\n", + " + Kd * Kds**2 * ST\n", + " + Kd * Kds * RT**2\n", + " - 2 * Kd * Kds * RT * ST\n", + " + Kd * Kds * ST**2\n", + " + Kds**2 * LT**2\n", + " - 2 * Kds**2 * LT * RT\n", + " + Kds**2 * RT**2\n", + " ),\n", + " (\n", + " Kd3 * Kd**2 * Kds**2\n", + " - Kd3 * Kd**2 * Kds * RT\n", + " + Kd3 * Kd**2 * Kds * ST\n", + " + Kd3 * Kd * Kds**2 * LT\n", + " - Kd3 * Kd * Kds**2 * RT\n", + " - Kd * Kds**2 * LT * RT\n", + " + Kd * Kds**2 * LT * ST\n", + " + Kd * Kds**2 * RT**2\n", + " - Kd * Kds**2 * RT * ST\n", + " ),\n", + " -Kd3 * Kd**2 * Kds**2 * RT,\n", + "]\n", + "\n", + "assert all(\n", + " sp.simplify(symbolic - runtime) == 0\n", + " for symbolic, runtime in zip(coefficients, runtime_coefficients, strict=True)\n", + ")\n" ] }, { "cell_type": "markdown", - "id": "bounds-explanation", + "id": "root-selection", "metadata": {}, "source": [ - "## The transformed coordinate and its physical interval\n", + "## Numeric root selection\n", "\n", - "In the three-state model, the equivalent coordinate collapses to literal free receptor. In the four-state model it generally does **not**. The physically valid interval for `r_equiv` is determined by two endpoint bound fractions:\n", + "For each competitor concentration, `bindcurve` evaluates the coefficient helper, finds all numeric roots, and keeps the effectively real root in the physical free-receptor interval\n", "\n", - "- the direct-binding limit at `LT = 0`, governed by `Kds`\n", - "- the asymptotic high-competitor limit, governed by `Kd3`\n", + "$$\n", + "0 \\le R \\le R_T.\n", + "$$\n", "\n", - "Those two endpoint fractions are converted into `r_equiv`, and that interval is used to filter the numeric roots of the quintic.\n" + "If `Kds == Kd3`, the leading coefficient is zero and the quintic degenerates to a quartic. The runtime helper trims leading near-zero coefficients before calling `numpy.roots`.\n" ] }, { "cell_type": "code", - "execution_count": 5, - "id": "bounds-demo", + "execution_count": null, + "id": "numeric-root", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Tracer-bound fraction at LT = 0: 0.699265\n", - "Asymptotic tracer-bound fraction at very large LT: 0.090163\n", - "r_equiv at LT = 0: 0.046504\n", - "Asymptotic r_equiv: 0.001982\n", - "Feasible interval used for root selection: [0.001982, 0.046504]\n", - "RT for comparison: 0.050000\n" - ] - } - ], + "outputs": [], "source": [ - "initial_fraction = _direct_bound_fraction(\n", - " Kd=params[\"Kds\"],\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - ")\n", - "asymptotic_fraction = _direct_bound_fraction(\n", - " Kd=params[\"Kd3\"],\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - ")\n", + "import numpy as np\n", "\n", - "initial_r_equiv = _equivalent_receptor_from_bound_fraction(\n", - " initial_fraction,\n", - " Kds=params[\"Kds\"],\n", + "from bindcurve.modeling.competitive import (\n", + " _competitive_four_state_coefficients,\n", + " _select_physical_root,\n", ")\n", - "asymptotic_r_equiv = _equivalent_receptor_from_bound_fraction(\n", - " asymptotic_fraction,\n", - " Kds=params[\"Kds\"],\n", + "\n", + "params = {\n", + " \"RT\": 0.05,\n", + " \"LsT\": 0.005,\n", + " \"Kds\": 0.02,\n", + " \"Kd\": 1.6,\n", + " \"Kd3\": 0.5,\n", + "}\n", + "\n", + "ligand_total = 0.1\n", + "numeric_coefficients = _competitive_four_state_coefficients(\n", + " ligand_total,\n", + " **params,\n", ")\n", - "lower_bound, upper_bound = _equivalent_receptor_bounds(\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - " Kds=params[\"Kds\"],\n", - " Kd3=params[\"Kd3\"],\n", + "\n", + "R_free = _select_physical_root(\n", + " numeric_coefficients,\n", + " lower_bound=0.0,\n", + " upper_bound=params[\"RT\"],\n", ")\n", "\n", - "print(f\"Tracer-bound fraction at LT = 0: {initial_fraction:.6f}\")\n", - "print(f\"Asymptotic tracer-bound fraction at very large LT: {asymptotic_fraction:.6f}\")\n", - "print(f\"r_equiv at LT = 0: {initial_r_equiv:.6f}\")\n", - "print(f\"Asymptotic r_equiv: {asymptotic_r_equiv:.6f}\")\n", - "print(f\"Feasible interval used for root selection: [{lower_bound:.6f}, {upper_bound:.6f}]\")\n", - "print(f\"RT for comparison: {params['RT']:.6f}\")\n" + "print(\"coefficients:\", numeric_coefficients)\n", + "print(\"roots:\", np.roots(numeric_coefficients))\n", + "print(\"selected R:\", R_free)\n", + "print(\"residual:\", np.polyval(numeric_coefficients, R_free))\n", + "assert 0.0 <= R_free <= params[\"RT\"]\n", + "assert np.isclose(np.polyval(numeric_coefficients, R_free), 0.0)\n" ] }, { "cell_type": "markdown", - "id": "quintic-explanation", + "id": "fraction-derivation", "metadata": {}, "source": [ - "## One competitor concentration: build the quintic and pick the physical root\n", + "## Bound tracer fraction from true free receptor\n", + "\n", + "For a fixed free receptor concentration `R`, the remaining free concentrations can be recovered without another nonlinear solve.\n", "\n", - "For each competitor concentration `LT`, `bindcurve` builds a quintic in `r_equiv`:\n", + "Define\n", "\n", "$$\n", - "a r_{\\mathrm{equiv}}^5 + b r_{\\mathrm{equiv}}^4 + c r_{\\mathrm{equiv}}^3 + d r_{\\mathrm{equiv}}^2 + e r_{\\mathrm{equiv}} + f = 0.\n", + "A = 1 + \\frac{R}{K_d^*},\n", "$$\n", "\n", - "The coefficients are the closed-form expressions stored in `competitive.py`. The solver finds all numeric roots, rejects nonphysical ones, and keeps the feasible root with the smallest scaled residual.\n" + "$$\n", + "B = 1 + \\frac{R}{K_d},\n", + "$$\n", + "\n", + "$$\n", + "C = \\frac{R}{K_dK_{d3}}.\n", + "$$\n", + "\n", + "The tracer balance gives\n", + "\n", + "$$\n", + "S_T = S(A + CL),\n", + "$$\n", + "\n", + "and the competitor balance gives\n", + "\n", + "$$\n", + "L_T = L(B + CS).\n", + "$$\n", + "\n", + "Substituting the first into the second gives a quadratic in free competitor `L`:\n", + "\n", + "$$\n", + "BC L^2 + (AB + CS_T - CL_T)L - AL_T = 0.\n", + "$$\n", + "\n", + "After selecting the non-negative root, the tracer-bound fraction can be evaluated as\n", + "\n", + "$$\n", + "F_b^* = \\frac{RS + RLS}{S_T}\n", + " = \\frac{R/K_d^* + CL}{1 + R/K_d^* + CL}.\n", + "$$\n", + "\n", + "In source-code names, this is implemented as `R/Kds + C*L`, where `C = R/(Kd*Kd3)`.\n" ] }, { "cell_type": "code", - "execution_count": 6, - "id": "single-point-demo", + "execution_count": null, + "id": "species-validation", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Quintic coefficients at LT = 0.1 uM:\n", - " a = -0.0001\n", - " b = -0.00016454\n", - " c = -2.699484e-06\n", - " d = 2.400551e-07\n", - " e = 7.656216e-09\n", - " f = 6.15128e-11\n", - "\n", - "All roots:\n", - " -1.62792952837 + 0j\n", - " 0.0439994824223 + 0j\n", - " -0.0214754533465 + 0j\n", - " -0.0200203685679 + 0j\n", - " -0.0199741321382 + 0j\n", - "\n", - "Selected physical root: r_equiv = 0.0439994824223\n", - "Recovered tracer-bound fraction: FSB = 0.687497472745\n", - "Polynomial residual at selected root: 1.965e-24\n" - ] - } - ], + "outputs": [], "source": [ - "LT = 0.1\n", - "coefficients = _competitive_four_state_coefficients(\n", - " LT,\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - " Kds=params[\"Kds\"],\n", - " Kd=params[\"Kd\"],\n", - " Kd3=params[\"Kd3\"],\n", - ")\n", + "def species_from_free_receptor(R_free, ligand_total, *, RT, LsT, Kds, Kd, Kd3):\n", + " a = 1.0 + R_free / Kds\n", + " b = 1.0 + R_free / Kd\n", + " c = R_free / (Kd * Kd3)\n", + "\n", + " quadratic_a = b * c\n", + " quadratic_b = a * b + c * LsT - c * ligand_total\n", + " quadratic_c = -a * ligand_total\n", + "\n", + " discriminant = quadratic_b**2 - 4.0 * quadratic_a * quadratic_c\n", + " if abs(quadratic_a) > np.finfo(float).tiny:\n", + " L_free = (-quadratic_b + np.sqrt(discriminant)) / (2.0 * quadratic_a)\n", + " else:\n", + " L_free = ligand_total / b\n", + "\n", + " S_free = LsT / (a + c * L_free)\n", + "\n", + " RS = R_free * S_free / Kds\n", + " RL = R_free * L_free / Kd\n", + " RLS = R_free * L_free * S_free / (Kd * Kd3)\n", + "\n", + " return {\n", + " \"R\": R_free,\n", + " \"S\": S_free,\n", + " \"L\": L_free,\n", + " \"RS\": RS,\n", + " \"RL\": RL,\n", + " \"RLS\": RLS,\n", + " }\n", "\n", - "print(f\"Quintic coefficients at LT = {LT} uM:\")\n", - "for name, coefficient in zip(\"abcdef\", coefficients, strict=True):\n", - " print(f\" {name} = {coefficient:.12g}\")\n", "\n", - "roots = np.roots(coefficients)\n", - "print(\"\\nAll roots:\")\n", - "for root in roots:\n", - " sign = \"+\" if np.imag(root) >= 0 else \"-\"\n", - " print(f\" {np.real(root): .12g} {sign} {abs(np.imag(root)):.12g}j\")\n", + "species = species_from_free_receptor(R_free, ligand_total, **params)\n", + "fraction_bound = (species[\"RS\"] + species[\"RLS\"]) / params[\"LsT\"]\n", "\n", - "selected_root = _select_physical_root(\n", - " coefficients,\n", - " lower_bound=lower_bound,\n", - " upper_bound=upper_bound,\n", - ")\n", - "fraction_bound = selected_root / (params[\"Kds\"] + selected_root)\n", + "print(species)\n", + "print(\"FSB:\", fraction_bound)\n", "\n", - "print(f\"\\nSelected physical root: r_equiv = {selected_root:.12g}\")\n", - "print(f\"Recovered tracer-bound fraction: FSB = {fraction_bound:.12g}\")\n", - "print(f\"Polynomial residual at selected root: {np.polyval(coefficients, selected_root):.3e}\")\n" + "assert np.isclose(\n", + " species[\"R\"] + species[\"RS\"] + species[\"RL\"] + species[\"RLS\"],\n", + " params[\"RT\"],\n", + ")\n", + "assert np.isclose(\n", + " species[\"S\"] + species[\"RS\"] + species[\"RLS\"],\n", + " params[\"LsT\"],\n", + ")\n", + "assert np.isclose(\n", + " species[\"L\"] + species[\"RL\"] + species[\"RLS\"],\n", + " ligand_total,\n", + ")\n" ] }, { "cell_type": "markdown", - "id": "model-evaluation-explanation", + "id": "model-evaluation", "metadata": {}, "source": [ "## End-to-end model evaluation\n", "\n", - "Once the physical `r_equiv` root is selected, `bindcurve` converts it to `FSB` and then to signal:\n", + "The model response is\n", "\n", "$$\n", - "y = ymin + (ymax - ymin)FSB.\n", + "y = ymin + (ymax - ymin)F_b^*.\n", "$$\n", "\n", - "The total-binding model does not use a separate polynomial form. It simply evaluates the same four-state machinery with `Kd_eff = (1 + N) * Kd`.\n" + "The total/nonspecific variant does not need a separate algebraic derivation. It uses the same free-receptor quintic, but with\n", + "\n", + "$$\n", + "K_{d,\\mathrm{eff}} = (1 + N)K_d.\n", + "$$\n" ] }, { "cell_type": "code", - "execution_count": 7, - "id": "end-to-end-demo", + "execution_count": null, + "id": "end-to-end", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Response table (percent signal):\n", - " LT_uM specific total_N0 total_with_N\n", - " 0.001000 67.923039 67.923039 67.925776\n", - " 0.005179 67.877994 67.877994 67.892154\n", - " 0.026827 67.645811 67.645811 67.718647\n", - " 0.138950 66.472860 66.472860 66.836925\n", - " 0.719686 61.110011 61.110011 62.686485\n", - " 3.727594 44.690967 44.690967 48.545368\n", - " 19.306977 24.571622 24.571622 27.505069\n", - "100.000000 15.785939 15.785939 16.658413\n", - "\n", - "max |specific - total(N=0)| = 0.000e+00\n", - "Effective Kd used by comp_4st_total = 2.160000 uM\n", - "max |total(N=0.35) - specific(Kd_eff)| = 0.000e+00\n" - ] - } - ], + "outputs": [], "source": [ + "from bindcurve.modeling.competitive import (\n", + " CompetitiveFourStateSpecificKdModel,\n", + " CompetitiveFourStateTotalKdModel,\n", + ")\n", + "\n", "concentration = np.logspace(-3, 2, 8)\n", "specific_model = CompetitiveFourStateSpecificKdModel()\n", "total_model = CompetitiveFourStateTotalKdModel()\n", "\n", "specific = specific_model.evaluate(\n", " concentration,\n", - " ymin=params[\"ymin\"],\n", - " ymax=params[\"ymax\"],\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - " Kds=params[\"Kds\"],\n", - " Kd=params[\"Kd\"],\n", - " Kd3=params[\"Kd3\"],\n", + " ymin=5.0,\n", + " ymax=95.0,\n", + " **params,\n", ")\n", - "total_zero_n = total_model.evaluate(\n", + "total_n0 = total_model.evaluate(\n", " concentration,\n", - " ymin=params[\"ymin\"],\n", - " ymax=params[\"ymax\"],\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - " Kds=params[\"Kds\"],\n", - " Kd=params[\"Kd\"],\n", - " Kd3=params[\"Kd3\"],\n", + " ymin=5.0,\n", + " ymax=95.0,\n", " N=0.0,\n", + " **params,\n", ")\n", + "\n", + "N = 0.35\n", "total = total_model.evaluate(\n", " concentration,\n", - " ymin=params[\"ymin\"],\n", - " ymax=params[\"ymax\"],\n", - " RT=params[\"RT\"],\n", - " LsT=params[\"LsT\"],\n", - " Kds=params[\"Kds\"],\n", - " Kd=params[\"Kd\"],\n", - " Kd3=params[\"Kd3\"],\n", - " N=params[\"N\"],\n", + " ymin=5.0,\n", + " ymax=95.0,\n", + " N=N,\n", + " **params,\n", ")\n", - "\n", - "effective_kd = (1 + params[\"N\"]) * params[\"Kd\"]\n", "specific_with_effective_kd = specific_model.evaluate(\n", " concentration,\n", - " ymin=params[\"ymin\"],\n", - " ymax=params[\"ymax\"],\n", + " ymin=5.0,\n", + " ymax=95.0,\n", " RT=params[\"RT\"],\n", " LsT=params[\"LsT\"],\n", " Kds=params[\"Kds\"],\n", - " Kd=effective_kd,\n", + " Kd=(1.0 + N) * params[\"Kd\"],\n", " Kd3=params[\"Kd3\"],\n", ")\n", "\n", - "table = pd.DataFrame(\n", - " {\n", - " \"LT_uM\": concentration,\n", - " \"specific\": specific,\n", - " \"total_N0\": total_zero_n,\n", - " \"total_with_N\": total,\n", - " }\n", - ")\n", + "assert np.allclose(total_n0, specific)\n", + "assert np.allclose(total, specific_with_effective_kd)\n", "\n", - "print(\"Response table (percent signal):\")\n", - "print(table.to_string(index=False, float_format=lambda value: f\"{value:0.6f}\"))\n", - "print()\n", - "print(f\"max |specific - total(N=0)| = {np.max(np.abs(specific - total_zero_n)):.3e}\")\n", - "print(f\"Effective Kd used by comp_4st_total = {effective_kd:.6f} uM\")\n", - "print(\n", - " f\"max |total(N={params['N']}) - specific(Kd_eff)| = \"\n", - " f\"{np.max(np.abs(total - specific_with_effective_kd)):.3e}\"\n", - ")\n" + "np.column_stack([concentration, specific, total])\n" ] }, { @@ -373,16 +514,14 @@ "source": [ "## Takeaway\n", "\n", - "The practical four-state workflow in `bindcurve` is:\n", + "The clean four-state implementation is now:\n", "\n", - "1. Build the quintic coefficients in `r_equiv` for each competitor concentration.\n", - "2. Solve the quintic numerically.\n", - "3. Keep the effectively real root inside the feasible transformed-coordinate interval.\n", - "4. Convert that root back to `FSB` with `FSB = r_equiv / (Kds + r_equiv)`.\n", - "5. Convert `FSB` to the observed response using `ymin` and `ymax`.\n", - "6. For `comp_4st_total`, use the same machinery with `Kd_eff = (1 + N) * Kd`.\n", + "1. Derive and solve the quintic in true free receptor concentration `R`.\n", + "2. Select the physical root in `0 <= R <= RT`.\n", + "3. Compute the observable tracer-bound fraction from the actual species `RS + RLS`.\n", + "4. Use the same machinery for `comp_4st_total`, with `Kd` replaced by `(1 + N) * Kd`.\n", "\n", - "So the important implementation choice is not the raw symbolic derivation itself, but the transformed coordinate, the physical root-selection rule, and the `Kd` scaling used by the total-binding variant.\n" + "This keeps the four-state implementation consistent with the three-state implementation, where the primary equilibrium coordinate is also true free receptor concentration.\n" ] } ], diff --git a/bindcurve/modeling/competitive.py b/bindcurve/modeling/competitive.py index 76506bf..f871206 100644 --- a/bindcurve/modeling/competitive.py +++ b/bindcurve/modeling/competitive.py @@ -19,52 +19,6 @@ def _competition_guess(compound: CompoundData) -> dict[str, float]: return {"ymin": ymin, "ymax": ymax, "Kd": kd_guess} -def _direct_bound_fraction(*, Kd: float, RT: float, LsT: float) -> float: - """Return the direct-binding bound tracer fraction for finite totals.""" - discriminant = (Kd + LsT + RT) ** 2 - 4.0 * LsT * RT - return float((Kd + LsT + RT - np.sqrt(discriminant)) / (2.0 * LsT)) - - -def _equivalent_receptor_from_bound_fraction( - fraction_bound: float, - *, - Kds: float, -) -> float: - """Transform FSB to the legacy four-state polynomial coordinate.""" - return float(Kds * fraction_bound / (1.0 - fraction_bound)) - - -def _equivalent_receptor_bounds( - *, - RT: float, - LsT: float, - Kds: float, - Kd3: float, -) -> tuple[float, float]: - """Return the physically feasible range for the transformed coordinate. - - The four-state polynomial is expressed in the coordinate - ``r_equiv = Kds * FSB / (1 - FSB)``. This is equal to free receptor only in - three-state/complete competition. In incomplete competition, the physical - interval is determined by the initial direct-binding fraction, governed by - ``Kds``, and the asymptotic bound fraction, governed by ``Kd3``. - """ - initial_fraction = _direct_bound_fraction(Kd=Kds, RT=RT, LsT=LsT) - asymptotic_fraction = _direct_bound_fraction(Kd=Kd3, RT=RT, LsT=LsT) - initial_coordinate = _equivalent_receptor_from_bound_fraction( - initial_fraction, - Kds=Kds, - ) - asymptotic_coordinate = _equivalent_receptor_from_bound_fraction( - asymptotic_fraction, - Kds=Kds, - ) - return ( - min(initial_coordinate, asymptotic_coordinate), - max(initial_coordinate, asymptotic_coordinate), - ) - - def _competitive_four_state_coefficients( ligand_total: float, *, @@ -74,12 +28,15 @@ def _competitive_four_state_coefficients( Kd: float, Kd3: float, ) -> np.ndarray: - """Return quintic coefficients for the four-state competitive model. + """Return quintic coefficients for true free receptor in the four-state model. + + The polynomial variable is literal free receptor concentration ``R``: + + ``a*R**5 + b*R**4 + c*R**3 + d*R**2 + e*R + f = 0``. - The polynomial variable is not literal free receptor in the four-state - model. It is the transformed coordinate ``r_equiv = Kds * FSB / (1 - FSB)``. - This coordinate lets the model recover ``FSB`` as - ``FSB = r_equiv / (Kds + r_equiv)`` after solving the quintic. + After the physical root is selected in ``0 <= R <= RT``, the observable + tracer-bound fraction is computed from the actual four-state species + ``RS + RLS`` rather than from a transformed receptor-like coordinate. The total/nonspecific model should call this with an effective ``Kd`` of ``(1 + N) * Kd`` rather than maintaining a duplicated coefficient @@ -87,17 +44,86 @@ def _competitive_four_state_coefficients( """ LT = float(ligand_total) - a = -Kds**2 * Kd3**2 - b = Kds**2 * Kd3 * (Kds * Kd - 3 * Kds * Kd3 + Kds * LT - 2 * Kds * LsT + Kds * RT - Kd * Kd3 - Kd3 * LT - Kd3 * LsT + Kd3 * RT) # noqa: E501 - c = Kds**2 * (3 * Kds**2 * Kd * Kd3 - 3 * Kds**2 * Kd3**2 + 3 * Kds**2 * Kd3 * LT - 4 * Kds**2 * Kd3 * LsT + 3 * Kds**2 * Kd3 * RT + Kds**2 * LT * LsT - Kds**2 * LT * RT - Kds**2 * LsT**2 + Kds**2 * LsT * RT - 3 * Kds * Kd * Kd3**2 + Kds * Kd * Kd3 * LsT - Kds * Kd * Kd3 * RT - 3 * Kds * Kd3**2 * LT - 2 * Kds * Kd3**2 * LsT + 3 * Kds * Kd3**2 * RT - Kds * Kd3 * LT * LsT + Kds * Kd3 * LT * RT - 2 * Kds * Kd3 * LsT**2 + 3 * Kds * Kd3 * LsT * RT - Kds * Kd3 * RT**2 - Kd * Kd3**2 * LsT + Kd * Kd3**2 * RT) # noqa: E501 - d = Kds**3 * (3 * Kds**2 * Kd * Kd3 - Kds**2 * Kd3**2 + 3 * Kds**2 * Kd3 * LT - 2 * Kds**2 * Kd3 * LsT + 3 * Kds**2 * Kd3 * RT + 2 * Kds**2 * LT * LsT - 3 * Kds**2 * LT * RT - Kds**2 * LsT**2 + 2 * Kds**2 * LsT * RT - 3 * Kds * Kd * Kd3**2 + 2 * Kds * Kd * Kd3 * LsT - 3 * Kds * Kd * Kd3 * RT - 3 * Kds * Kd3**2 * LT - Kds * Kd3**2 * LsT + 3 * Kds * Kd3**2 * RT - 2 * Kds * Kd3 * LT * LsT + 3 * Kds * Kd3 * LT * RT - 2 * Kds * Kd3 * LsT**2 + 6 * Kds * Kd3 * LsT * RT - 3 * Kds * Kd3 * RT**2 - Kds * LsT**3 + 2 * Kds * LsT**2 * RT - Kds * LsT * RT**2 - 2 * Kd * Kd3**2 * LsT + 3 * Kd * Kd3**2 * RT) # noqa: E501 - e = Kds**4 * (Kds**2 * Kd * Kd3 + Kds**2 * Kd3 * LT + Kds**2 * Kd3 * RT + Kds**2 * LT * LsT - 3 * Kds**2 * LT * RT + Kds**2 * LsT * RT - Kds * Kd * Kd3**2 + Kds * Kd * Kd3 * LsT - 3 * Kds * Kd * Kd3 * RT - Kds * Kd3**2 * LT + Kds * Kd3**2 * RT - Kds * Kd3 * LT * LsT + 3 * Kds * Kd3 * LT * RT + 3 * Kds * Kd3 * LsT * RT - 3 * Kds * Kd3 * RT**2 + 2 * Kds * LsT**2 * RT - 2 * Kds * LsT * RT**2 - Kd * Kd3**2 * LsT + 3 * Kd * Kd3**2 * RT) # noqa: E501 - f = Kds**5 * RT * (-Kds**2 * LT - Kds * Kd * Kd3 + Kds * Kd3 * LT - Kds * Kd3 * RT - Kds * LsT * RT + Kd * Kd3**2) # noqa: E501 + a = Kds - Kd3 + b = ( + -Kd3 * Kd + - Kd3 * Kds + - Kd3 * LT + + Kd3 * RT + - Kd3 * LsT + + Kd * Kds + + Kds**2 + + Kds * LT + - 2.0 * Kds * RT + + Kds * LsT + ) + c = ( + Kd3 * Kd * RT + - Kd3 * Kd * LsT + - Kd3 * Kds * LT + + Kd3 * Kds * RT + + Kd * Kds**2 + - 2.0 * Kd * Kds * RT + + 2.0 * Kd * Kds * LsT + + 2.0 * Kds**2 * LT + - 2.0 * Kds**2 * RT + - Kds * LT * RT + + Kds * LT * LsT + + Kds * RT**2 + - Kds * RT * LsT + ) + d = ( + Kd3 * Kd**2 * Kds + + Kd3 * Kd * Kds**2 + + Kd3 * Kd * Kds * LT + + Kd3 * Kd * Kds * LsT + + Kd * Kds**2 * LT + - 2.0 * Kd * Kds**2 * RT + + Kd * Kds**2 * LsT + + Kd * Kds * RT**2 + - 2.0 * Kd * Kds * RT * LsT + + Kd * Kds * LsT**2 + + Kds**2 * LT**2 + - 2.0 * Kds**2 * LT * RT + + Kds**2 * RT**2 + ) + e = ( + Kd3 * Kd**2 * Kds**2 + - Kd3 * Kd**2 * Kds * RT + + Kd3 * Kd**2 * Kds * LsT + + Kd3 * Kd * Kds**2 * LT + - Kd3 * Kd * Kds**2 * RT + - Kd * Kds**2 * LT * RT + + Kd * Kds**2 * LT * LsT + + Kd * Kds**2 * RT**2 + - Kd * Kds**2 * RT * LsT + ) + f = -Kd3 * Kd**2 * Kds**2 * RT return np.array([a, b, c, d, e, f], dtype=float) +def _trim_leading_near_zero( + coefficients: np.ndarray, + *, + relative_tolerance: float = 1.0e-14, +) -> np.ndarray: + """Remove leading coefficients that are zero at working precision.""" + coefficients = np.asarray(coefficients, dtype=float) + scale = float(np.max(np.abs(coefficients))) if coefficients.size else 0.0 + if scale == 0.0: + return np.array([0.0], dtype=float) + + threshold = relative_tolerance * scale + for index, coefficient in enumerate(coefficients): + if abs(float(coefficient)) > threshold: + return coefficients[index:] + + return np.array([0.0], dtype=float) + + def _scaled_polynomial_residual(coefficients: np.ndarray, root: float) -> float: + coefficients = _trim_leading_near_zero(coefficients) degree = len(coefficients) - 1 root_scale = max(1.0, abs(root)) denominator = sum( @@ -117,13 +143,14 @@ def _select_physical_root( imaginary_tolerance: float = 1.0e-7, interval_tolerance: float = 1.0e-8, ) -> float: - """Select the physical four-state transformed-coordinate root. + """Select the physical four-state free-receptor root. The physical root must be effectively real and lie in the feasible interval - for ``r_equiv = Kds * FSB / (1 - FSB)``. This interval is not generally - ``0 <= r_equiv <= RT``. Among feasible candidates, the root with the - smallest scaled polynomial residual is selected. + for literal free receptor concentration. For the four-state receptor + polynomial this interval is ``0 <= R <= RT``. Among feasible candidates, the + root with the smallest scaled polynomial residual is selected. """ + coefficients = _trim_leading_near_zero(coefficients) roots = np.roots(coefficients) interval_scale = max(1.0, abs(lower_bound), abs(upper_bound)) lower = lower_bound - interval_tolerance * interval_scale @@ -140,9 +167,9 @@ def _select_physical_root( if not candidates: raise ValueError( - "No physical four-state root found in the feasible transformed " - f"coordinate interval {lower_bound} <= r_equiv <= {upper_bound}. " - f"Roots were: {roots!r}" + "No physical four-state root found in the feasible free-receptor " + f"interval {lower_bound} <= R <= {upper_bound}. Roots were: " + f"{roots!r}" ) return min( @@ -151,7 +178,7 @@ def _select_physical_root( ) -def _competitive_four_state_equivalent_receptor( +def _competitive_four_state_receptor_free( ligand_total: np.ndarray, *, RT: float, @@ -161,14 +188,11 @@ def _competitive_four_state_equivalent_receptor( Kd3: float, ) -> np.ndarray: ligand_total = np.asarray(ligand_total, dtype=float) + if RT == 0.0: + return np.zeros_like(ligand_total, dtype=float) + flat_ligand_total = ligand_total.ravel() - lower_bound, upper_bound = _equivalent_receptor_bounds( - RT=RT, - LsT=LsT, - Kds=Kds, - Kd3=Kd3, - ) - equivalent_receptor = [] + receptor_free = [] for concentration in flat_ligand_total: coefficients = _competitive_four_state_coefficients( @@ -179,15 +203,65 @@ def _competitive_four_state_equivalent_receptor( Kd=Kd, Kd3=Kd3, ) - equivalent_receptor.append( + receptor_free.append( _select_physical_root( coefficients, - lower_bound=lower_bound, - upper_bound=upper_bound, + lower_bound=0.0, + upper_bound=RT, ) ) - return np.asarray(equivalent_receptor, dtype=float).reshape(ligand_total.shape) + return np.asarray(receptor_free, dtype=float).reshape(ligand_total.shape) + + +def _competitive_four_state_fraction_tracer_bound( + receptor_free: np.ndarray, + ligand_total: np.ndarray, + *, + LsT: float, + Kds: float, + Kd: float, + Kd3: float, +) -> np.ndarray: + """Return ``(RS + RLS) / LsT`` from true free receptor concentration. + + For fixed free receptor ``R``, the tracer and competitor mass balances can + be reduced to a quadratic equation in free competitor concentration ``L``. + The tracer-bound fraction is then obtained from the actual four-state + species without explicitly constructing free tracer concentration. + """ + receptor_free = np.asarray(receptor_free, dtype=float) + ligand_total = np.asarray(ligand_total, dtype=float) + + a = 1.0 + receptor_free / Kds + b = 1.0 + receptor_free / Kd + c = receptor_free / (Kd * Kd3) + + quadratic_a = b * c + quadratic_b = a * b + c * LsT - c * ligand_total + quadratic_c = -a * ligand_total + + discriminant = np.maximum( + quadratic_b**2 - 4.0 * quadratic_a * quadratic_c, + 0.0, + ) + denominator = 2.0 * quadratic_a + fallback = np.divide( + ligand_total, + b, + out=np.zeros_like(ligand_total, dtype=float), + where=b != 0.0, + ) + ligand_free = np.array(fallback, copy=True, dtype=float) + np.divide( + -quadratic_b + np.sqrt(discriminant), + denominator, + out=ligand_free, + where=np.abs(denominator) > np.finfo(float).tiny, + ) + + bound_tracer_ratio = receptor_free / Kds + c * ligand_free + return bound_tracer_ratio / (1.0 + bound_tracer_ratio) class CompetitiveFourStateSpecificKdModel(BaseDoseResponseModel): @@ -219,7 +293,7 @@ def evaluate( Kd3: float, Kd: float, ) -> np.ndarray: - equivalent_receptor = _competitive_four_state_equivalent_receptor( + receptor_free = _competitive_four_state_receptor_free( x, RT=RT, LsT=LsT, @@ -227,7 +301,14 @@ def evaluate( Kd=Kd, Kd3=Kd3, ) - fraction_tracer_bound = equivalent_receptor / (Kds + equivalent_receptor) + fraction_tracer_bound = _competitive_four_state_fraction_tracer_bound( + receptor_free, + x, + LsT=LsT, + Kds=Kds, + Kd=Kd, + Kd3=Kd3, + ) return ymin + (ymax - ymin) * fraction_tracer_bound def guess(self, compound: CompoundData) -> dict[str, float]: @@ -265,15 +346,23 @@ def evaluate( N: float, Kd: float, ) -> np.ndarray: - equivalent_receptor = _competitive_four_state_equivalent_receptor( + effective_kd = (1.0 + N) * Kd + receptor_free = _competitive_four_state_receptor_free( x, RT=RT, LsT=LsT, Kds=Kds, - Kd=(1.0 + N) * Kd, + Kd=effective_kd, + Kd3=Kd3, + ) + fraction_tracer_bound = _competitive_four_state_fraction_tracer_bound( + receptor_free, + x, + LsT=LsT, + Kds=Kds, + Kd=effective_kd, Kd3=Kd3, ) - fraction_tracer_bound = equivalent_receptor / (Kds + equivalent_receptor) return ymin + (ymax - ymin) * fraction_tracer_bound def guess(self, compound: CompoundData) -> dict[str, float]: diff --git a/pyproject.toml b/pyproject.toml index 2758566..65a9387 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,7 @@ lint = [ notebook = [ "ipykernel", "jupyterlab", + "sympy", ] dev = [ {include-group = "test"}, diff --git a/tests/test_competitive_four_state_models.py b/tests/test_competitive_four_state_models.py index 9f6bbde..d88dabb 100644 --- a/tests/test_competitive_four_state_models.py +++ b/tests/test_competitive_four_state_models.py @@ -7,7 +7,6 @@ import bindcurve as bc from bindcurve.modeling.competitive import ( _competitive_four_state_coefficients, - _equivalent_receptor_bounds, _select_physical_root, ) @@ -97,7 +96,7 @@ def test_registry_contains_competitive_four_state_models(): ) -def test_four_state_root_selector_returns_physical_transformed_root(): +def test_four_state_root_selector_returns_physical_free_receptor_root(): kwargs = { "RT": 0.05, "LsT": 0.005, @@ -106,31 +105,36 @@ def test_four_state_root_selector_returns_physical_transformed_root(): "Kd3": 0.5, } coefficients = _competitive_four_state_coefficients(0.1, **kwargs) - lower_bound, upper_bound = _equivalent_receptor_bounds( - RT=kwargs["RT"], - LsT=kwargs["LsT"], - Kds=kwargs["Kds"], - Kd3=kwargs["Kd3"], - ) root = _select_physical_root( coefficients, - lower_bound=lower_bound, - upper_bound=upper_bound, + lower_bound=0.0, + upper_bound=kwargs["RT"], ) - assert lower_bound <= root <= upper_bound + assert 0.0 <= root <= kwargs["RT"] assert np.isclose(np.polyval(coefficients, root), 0.0, atol=1.0e-10) -def test_four_state_root_bounds_can_extend_beyond_receptor_total(): - lower_bound, upper_bound = _equivalent_receptor_bounds( - RT=0.05, - LsT=0.005, - Kds=0.02, - Kd3=0.001, +def test_four_state_root_selector_handles_degenerate_quartic_case(): + kwargs = { + "RT": 0.05, + "LsT": 0.005, + "Kds": 0.02, + "Kd": 0.8, + "Kd3": 0.02, + } + coefficients = _competitive_four_state_coefficients(0.1, **kwargs) + + assert coefficients[0] == 0.0 + + root = _select_physical_root( + coefficients, + lower_bound=0.0, + upper_bound=kwargs["RT"], ) - assert lower_bound < 0.05 < upper_bound + assert 0.0 <= root <= kwargs["RT"] + assert np.isclose(np.polyval(coefficients, root), 0.0, atol=1.0e-10) def test_four_state_root_selector_rejects_polynomial_with_no_physical_root():