affinityfit

🚧 This library is currently under development and may change substantially.

A Python library that fits K_d and related parameters from concentration and signal data, and diagnoses whether the measurements constrain the resulting estimate.

Any observable that is linear in the fraction bound follows the same model (a saturation curve), regardless of the measurement technique. Nuclear magnetic resonance (NMR) peak intensity, surface plasmon resonance (SPR) steady-state response, and initial enzyme velocity are all examples.

signal = baseline + Bmax * [L] / (Kd + [L])

Several datasets can be fitted simultaneously, sharing parameters across them or holding some constant. A curve whose measured range does not bracket K_d is undetermined on its own, but sharing a parameter with a better-sampled dataset can render the estimate identifiable.

Usage

uv add affinityfit

import numpy as np
from affinityfit import DiagnosticCode, fit

conc, signal = np.loadtxt(
    "titration.csv",
    delimiter=",",
    skiprows=1,
    unpack=True,
)
res = fit(conc, signal, unit="nM")

print(res.params["kd"], res.intervals["kd"].format("nM"))

# `code` is the stable programmatic contract. `severity` classifies the finding.
for diagnostic in res.diagnostics:
    print(diagnostic.code, diagnostic.severity, diagnostic.message)

if any(diagnostic.code == DiagnosticCode.NOT_SATURATED for diagnostic in res.diagnostics):
    # Extend the concentration range before interpreting Kd or Bmax.
    pass

# Optional, caller-controlled text rendering for a terminal or notebook.
print(res.report())

res.diagnostics is the canonical diagnostics API. Each immutable Diagnostic has a stable code for program logic, a severity of "warning" or "note", and a concise English message for display. Do not branch on message: its prose may be refined without changing the code. res.warnings and res.notes are severity-filtered views containing the same Diagnostic objects.

DiagnosticCode is an enum listing every code this library can emit, so list(DiagnosticCode) enumerates them without reading the source. Each member is a plain string (DiagnosticCode.NOT_SATURATED == "not_saturated"), so comparing diagnostic.code against either the enum member or the equivalent string literal works the same way.

The CSV has one header row, with concentration in the first column and signal in the second.

concentration_nM,signal
0,0.0249
1,0.0963
3,0.2493
10,0.5090
30,0.7657
100,0.9249
300,0.9758
1000,1.0035

Output of the optional report() renderer:

model    : langmuir  (1:1 binding: signal = baseline + Bmax * L / (Kd + L))
interval : profile
Kd       = 10.1 +/- 0.8 nM
Bmax     = 0.995 +/- 0.017
baseline = 0.018 +/- 0.014
AICc     = -63.79   (AIC = -77.12)
R^2      = 0.9998   (n = 8; descriptive only)

No diagnostic issues detected.

If the concentration of the fixed partner (for example, a receptor or a lectin) is known, pass it in to enable the ligand-depletion check.

res = fit(conc, signal, receptor_conc=1.0, unit="nM")

Pass model= to switch models.

from affinityfit import hill

res = fit(conc, signal, model=hill)  # test for cooperativity
print(res.intervals["n"].contains(1.0))

Fitting several datasets simultaneously, sharing or fixing parameters:

from affinityfit import Dataset, fit_global

res = fit_global(
    [Dataset("oxidized", conc, sig_ox), Dataset("reduced", conc, sig_red)],
    shared=["bmax"],  # estimate a single value across all datasets
    fixed={"baseline": 0.0},  # hold at a constant, not estimated
    unit="mM",
)
# Preserve diagnostic scope in a global fit.
for diagnostic in res.fit_diagnostics:
    print("fit", diagnostic.code, diagnostic.severity)
for name, diagnostics in res.diagnostics_per.items():
    for diagnostic in diagnostics:
        print(name, diagnostic.code, diagnostic.severity)

# Pull out one dataset. Its local and fit-wide diagnostics come along with it.
sub = res.result_for("oxidized")
x, y = sub.curve()

# Optional caller-controlled rendering.
print(res.report())

fit_diagnostics contains findings about the complete fit, while diagnostics_per keeps each dataset's findings scoped to its name. result_for() combines those fit-wide and local records into its FitResult. warnings and notes are severity-filtered views over both scopes at once. warnings_per and notes_per provide the corresponding view per dataset and include fit-wide findings in every dataset entry.

Weighting

By default every point counts equally, which assumes the measurement error is the same size everywhere. That holds for a signal measured in resonance units, as in SPR, but not for fluorescence, luminescence, or absorbance, where the error scales with the signal. Since a saturation curve moves the signal from baseline to baseline + B_max, the absolute error near saturation can be tens of times larger than at low concentration.

res = fit(conc, signal, sigma=0.01 + 0.10 * signal)  # per-point standard deviation

sigma is treated as relative; the overall scale is still estimated from the residuals, so the result is unchanged under a constant rescaling. In a global fit, Dataset(..., sigma=...) can be given per dataset, which also sets the relative weight between datasets (for example, one dataset ten times noisier than another).

sigma must be given for either every dataset or none of them. Giving it to only some would leave the others with an implicit weight of 1, so the relative weight between datasets would depend on the absolute scale of sigma (whether it is written as a fraction, a percentage, or ppm). Mixing the two raises an error.

Fitting heteroscedastic data without weights costs more than precision: it also narrows the confidence interval. Over 200 simulated 12-point titrations with 30% proportional error and K_d near the top of the measured range, 95% interval coverage was 79% unweighted against 94% with sigma supplied.

The diagnostics warn when the size of the residuals grows with the fitted value.

Confidence intervals

ci= selects one of three methods.

replicates has shape (n_replicates, n_points).

res = fit(conc, signal, ci="bootstrap", replicates=reps, n_boot=2000)

Dataset also accepts replicates on their own (signal is then their mean).

Dataset("oxidized", conc, replicates=reps)

Under ci="bootstrap", n_boot must be at or above the minimum needed for a percentile interval (100), or it raises an error; below that, no interval can be formed.

The reported precision follows the uncertainty. A measurement with a 17% relative error does not warrant three significant figures, so Kd = 3.80e-8 is reported as (3.8 +/- 0.6)e-08 instead.

Model

bmax in langmuir and hill is not restricted in sign. For an observable that decreases on binding (fluorescence quenching, a drop in intensity), K_d retains its meaning and the decrease is expressed as a negative response coefficient in the same equation. vmax in michaelis is a rate and stays non-negative, so use langmuir or hill for decreasing data. A model of the wrong sign yields a diagnostic rather than a plausible-looking number.

Km = (k_off + k_cat)/k_on, and only equals K_d when k_cat is negligible.

K_d, K_m, and IC₅₀ are optimised on a logarithmic scale. On a linear scale the fitted result depends on the unit of concentration: a picomolar affinity (K_d = 1 × 10⁻¹² M) expressed in molar units comes out about 6700% too high on a 12-point titration. On the logarithmic scale K_d is recovered to within 10⁻⁶ relative error across 18 orders of magnitude, from fM to 100 M, and the same experiment written in M, nM, or pM gives the same answer.

Dose–response and K_i

ic50 is the four-parameter logistic. It is the same equation as hill, under the names used to report a dose–response curve, and without the cooperativity checks, since a slope near 1 is the ordinary case in a dose–response experiment rather than a finding.

from affinityfit import ic50

res = fit(conc, response, model=ic50, unit="nM")
print(res.params["ic50"], res.params["hillslope"])

A displacement IC₅₀ is not a property of the competitor alone: raising the concentration of the labelled ligand being displaced raises the IC₅₀ with it. Dividing that dependence out makes K_i comparable between assays run at different tracer concentrations.

from affinityfit import Interval, ki_from_ic50

ki = ki_from_ic50(
    res,
    tracer_conc=5.0,
    tracer_kd=Interval(point=2.0, lower=1.6, upper=2.4),
    receptor_conc=3.45,  # selects the exact correction
)
print(ki.format("nM"))

Pass the whole result rather than one interval out of it, and give the tracer constant as an Interval when its uncertainty is known.

Give receptor_conc whenever it is known. Without it the standard Cheng–Prusoff form IC50 / (1 + [T]/Kd) is used, which assumes the receptor is dilute enough that it neither depletes the tracer nor takes up an appreciable share of the competitor. Where that fails the form is biased, and the bias does not decay for a weak competitor: it settles at an offset set by how much tracer the receptor holds. Measured for a receptor at 6.5 times the tracer constant, the approximation came out 96% high at K_i = 0.05 and 9% low at K_i = 100, in the same units. Choosing a different tracer concentration to plug in does not fix it: the exact denominator falls between the ones the total and the free tracer concentration give. With receptor_conc the exact correction is used instead, which agrees with the Munson–Rodbard result. fit(..., receptor_conc=...) is carried on the FitResult, so passing the result is enough. Without it, a UserWarning reports that the bias could not be assessed.

The standard form assumes a slope of 1. It is derived for a single site under competition. A bounded, non-zero-width fitted-slope interval that excludes 1 indicates that the derivation does not apply. The modified forms that cover a slope away from 1 raise the terms to powers and do not agree with one another, so the choice between them is left to the caller; the library reports the violation rather than applying one silently. Ligand depletion steepens a displacement curve, so it is the first cause to rule out before adopting a modified form. Only the FitResult carries the IC₅₀, slope, and receptor concentration together, so only that overload can check the slope and inherit receptor_conc. Passing res.intervals["ic50"] skips the slope check and uses the standard approximation unless receptor_conc is passed explicitly.

The tracer constant's own error can dominate. Under the standard Cheng–Prusoff form, with r = [T]/Kd*, the relative error of Kd* propagates to K_i damped by r/(1+r). At [T] = 2.5·Kd*, a 20% uncertainty on Kd* therefore puts 14% onto K_i, against about 5% from a well-measured IC₅₀. Treating Kd* as exact here reports ±5% where the experiment supports ±15%. The exact correction instead evaluates the sensitivity to Kd* through the depleted-tracer relation.

When the mapped IC₅₀ interval brackets its point estimate, the two independent uncertainties combine in quadrature on each side. A percentile interval that excludes its point estimate is widened directly at its outer endpoints instead. In either case, asymmetry is preserved, an undetermined limit stays undetermined, and a lower limit driven past zero is reported as undetermined rather than as a K_i of zero. tracer_conc and receptor_conc are taken as exact.

For competitive enzyme inhibition the same expression applies with [S] and Km.

Both forms also assume the competitor and the tracer exclude each other from one site. A fit cannot check that, because the curve has the same shape either way.

Ligand depletion

Every model in the table except tight_binding assumes the free ligand concentration equals the total one. This assumption fails once the receptor is not much more dilute than K_d, because each molecule bound is one fewer left in solution. The curve still looks like a saturation curve, and R² alone does not reveal the problem. On noiseless 1:1 data with the receptor at five times K_d (15 points: 14 positive concentrations spanning 10⁻² to 10² times K_d and a blank, baseline held at 0), langmuir reports K_d = 3.99 against a true 1.0, at R² = 0.9919; RESIDUAL_STRUCTURE and HETEROSCEDASTIC are raised on the same fit and point at the mismatch that R² alone misses.

tight_binding solves the 1:1 equilibrium without that assumption. Pass the total receptor concentration as rt, normally as a constant.

from affinityfit import tight_binding

res = fit(conc, signal, model=tight_binding, fixed={"rt": 5.0}, unit="uM")

When rt is left free instead, it is estimated from the data and measures the active concentration: the fraction of the immobilised or pipetted receptor that is actually binding. Depletion changes the shape of the curve and not only its midpoint, which prevents rt and K_d from trading off against each other.

res = fit(conc, signal, model=tight_binding, unit="uM")
print(res.intervals["rt"].format("uM"))

The diagnostic that recommends this model (see below) is suppressed once the model is in use. Two limitations apply:

• fixed= applies one value to every dataset, so in a global fit spanning several receptor concentrations rt cannot be fixed per dataset. Fit those separately.
• The checks on where the measurements sit compare the range against K_d. Depletion reaches the plateau at a lower multiple of K_d than a hyperbola does, so those thresholds are conservative here and may still call for higher concentrations than this model requires.

Cooperativity under depletion

Depletion does not only shift the curve, it steepens it, and hill reads that steepness as an exponent. The same data as above, with no cooperativity anywhere in the system, fits hill at n = 1.43 [1.32, 1.54]. hill reduces to langmuir exactly at n = 1, so under unweighted least squares its R² cannot be lower and does not itself argue for cooperativity. The displayed unweighted R² does not have this guarantee when sigma is supplied. The interval on n excluding 1 is the relevant check, and intervals["n"].contains(1.0) accordingly reports cooperativity that does not exist.

No model here solves depletion and cooperativity at once. Adding one would impose an exact conservation law on the Hill equation, which is phenomenological to begin with. Cooperativity must therefore be measured where depletion is absent: keep the receptor at or below K_d/10 for that question.

HILL_N_ABOVE_ONE reports an n significantly above 1 together with the alternatives that produce the same shape: depletion, self-association, and a reading taken before equilibrium. Its message does not change with receptor_conc. Supplying receptor_conc enables the separate LIGAND_DEPLETION check, which is emitted when the receptor concentration exceeds K_d/10.

Model selection

Use aicc (the corrected Akaike information criterion) to compare models or parameter-sharing schemes. aic is kept for reference, but the uncorrected AIC is only asymptotically valid and favours the model with more parameters at the scale of a typical titration (6–15 points, 3–6 fitted coefficients, n/k of 2–5 once the residual-variance parameter below is counted). Lower AICc is preferred; the ordering is determined by the data rather than by the model name or sharing scheme.

import numpy as np
from affinityfit import Dataset, fit_global, hill, langmuir

conc = np.concatenate(([0.0], np.logspace(-1, 2, 11)))
rng = np.random.default_rng(7)


def simulate(kd):
    clean = langmuir(conc, kd, 1.0, 0.02)
    return clean + rng.normal(0.0, 0.01, conc.size)


datasets = [
    Dataset("sample_a", conc, simulate(3.0)),
    Dataset("sample_b", conc, simulate(30.0)),
]
free_langmuir = fit_global(datasets, model=langmuir)
shared_langmuir = fit_global(datasets, model=langmuir, shared=["bmax"])
free_hill = fit_global(datasets, model=hill)

print({"free": free_langmuir.aicc, "shared": shared_langmuir.aicc})
print({"langmuir": free_langmuir.aicc, "hill": free_hill.aicc})

report() shows AICc first, with AIC in parentheses for reference. Least squares also estimates the residual variance, so k here is one more than the number of fitted coefficients. When the sample is too small for the correction to be defined (n − k − 1 ≤ 0), aicc is infinite, and FitResult.report() drops the line instead of printing an infinity.

Diagnostics

K_d is frequently undetermined even when R² exceeds 0.99. The conditions below are detected automatically. A quantity whose name changes by model (K_d, K_m) is shown under that model's own label.

The diagnostics fall into two groups: whether the fit itself is sound, and whether the measurements are placed well.

Diagnostic.code (a DiagnosticCode member) is the value to branch on; see DiagnosticCode for how to list every code from Python.

ki_from_ic50 raises a UserWarning rather than a diagnostic, since it is called after the fit and has no FitResult of its own to attach one to.

In a global fit, sharing the amplitude suppresses NOT_SATURATED and WEAKLY_SATURATED, because the plateau is constrained across datasets. An unsaturated dataset receives SHARED_AMPLITUDE_IDENTIFIES_LOCATION instead. Fixing the baseline suppresses NO_LOW_CONC. The shared-amplitude note is withheld when the fit itself is broken (NO_FIT or AMPLITUDE_COLLAPSED), since it would otherwise contradict the other findings.

In examples/titration_unsaturated.csv, R² = 0.9936 suggests a good fit, but the highest concentration is only a third of K_d, so the upper limit is undetermined and the result can only be reported as Kd > 76 nM.

When a confidence interval cannot be computed (zero degrees of freedom, a rank-deficient Jacobian, an inner refit that fails to converge, or too few bootstrap resamples succeeding), no value is substituted; the interval is reported as undetermined, with the reason given. The point estimate itself is still returned.

Residuals are tested with both a runs test (Wald–Wolfowitz) and lag-1 autocorrelation, and a warning is emitted when either crosses its threshold. Over 2000 simulated 15-point titrations fitted with the correct model, the false-positive rate was 1.1%.

Statistics

The thresholded statistical checks do not share one criterion. model_vs_constant and heteroscedasticity use 0.01, while residual_runs uses 0.025. When every non-zero residual has the same sign, residual_sign_test records the exact two-sided P value, and RESIDUAL_STRUCTURE is emitted without comparing it with a threshold. residual_autocorrelation is compared with a fixed correlation threshold of 0.3 rather than a P value. Stacking tests across several datasets raises the false-positive rate, so the underlying statistics and P values are returned for a caller-selected multiple-comparison correction.

res = fit(conc, signal)
for s in res.statistics:
    print(s.name, s.statistic, s.p_value)

For several datasets, these are available per dataset through GlobalFitResult.statistics_per["dataset name"]. The names are model_vs_constant, residual_runs, residual_autocorrelation, and heteroscedasticity, plus residual_sign_test, which stands in for residual_runs when every residual shares a sign and the runs count is degenerate. residual_autocorrelation has no P value, so its p_value is None.

Plotting

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(*res.curve(), "-")
ax.plot(conc, signal, "o")
ax.set_xscale("log")

Keep the concentration axis logarithmic. On a linear axis the low-concentration region collapses, hiding both the curvature near K_d and design issues such as the measurements being skewed toward saturation.

examples/plot_fit.py is a sample that draws the measured points, the fitted curve, and a residual panel in a single figure. It is not part of the distribution. Run it from a source checkout with the development dependencies installed; Matplotlib is not installed with the library runtime dependencies.

uv run python examples/plot_fit.py examples/titration_good.csv --unit nM --out examples/fit_good.png

Running this sample on examples/titration_good.csv: