affinityfit

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

A Python library that fits Kd and related parameters from concentration and signal data, and diagnoses whether the resulting estimate can be trusted.

Any observable that is linear in the fraction bound follows the same model (a saturation curve), regardless of the measurement technique. NMR peak intensity, 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 Kd is otherwise undetermined on its own; sharing a parameter with a better-sampled dataset can make the estimate identifiable.

Usage

uv add affinityfit

from affinityfit import DiagnosticCode, fit, load_csv

conc, signal = load_csv("titration.csv")
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) is how to discover them all 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 input CSV has concentration in the first column and signal in the second. Header rows, comment rows, and blank rows are skipped. Repeated rows at the same concentration are used as replicates as-is.

Missing values (blank cells), nan, and inf raise an error with the row number. Silently dropping a blank row would change the number of points, and with it the degrees of freedom and the diagnostics.

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
R^2      = 0.9998   (n = 8)
AICc     = -71.12   (AIC = -77.12)

No diagnostic issues detected.

If the concentration of the fixed partner (a receptor, a lectin, and so on) 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, notes, warnings_per, and notes_per are structured severity-filtered convenience views; use the scoped collections above when the dataset matters.

Weighting

By default every point counts equally, which assumes the measurement error is the same size everywhere. That holds for a system like SPR response units, 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+Bmax, 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 carries over the relative weight between datasets (one being ten times noisier than another, say).

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 be decided by 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 (measured for Kd=10 with 30% proportional error, 95% CI coverage was 79% unweighted versus 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.

Reporting error from repeated experiments is standard practice in this field. Passing replicate measurements follows that convention.

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)

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 = 4.70e-8 is reported as (4.7 +/- 0.8)e-08 instead.

Model

bmax in langmuir and hill is not restricted in sign. An observable that decreases on binding (fluorescence quenching, a drop in intensity) still has the same meaning for Kd, 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. Fitting the wrong kind of model to the data warns rather than returning a plausible-looking number.

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

Kd, Km and IC50 are optimised on a logarithmic scale. Handled linearly, the fitted result would depend on the unit of concentration, reaching a relative error of 6800% for a picomolar affinity (Kd = 1e-12 M) expressed in molar units. As it stands, Kd is recovered to within 1e-6 relative error across 18 orders of magnitude, from fM to 100 M, giving the same answer whether the same experiment is written in M, nM, or pM.

Dose-Response and Ki

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: a slope near 1 is the ordinary case there 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 IC50 is not a property of the competitor alone: raising the concentration of the labelled ligand being displaced raises the IC50 with it. Dividing by 1 + [tracer]/Kd_tracer removes that dependence, which is what makes Ki 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))
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. Both matter:

The standard form assumes a slope of 1. It is derived for a single site under competition, and a fitted slope whose interval excludes 1 says that 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 choosing between them is yours rather than this library's; what the library does is refuse to stay quiet. Only the FitResult carries the IC50 and the slope together, so only that overload can check it. Passing res.intervals["ic50"] computes the same number without the check.

The tracer constant's own error usually dominates. With r = [T]/Kd*, its relative error reaches Ki damped by r/(1+r), so at [T] = 2.5·Kd* a 20% uncertainty on Kd* puts 14% onto Ki, against about 5% from a well-measured IC50. Treating Kd* as exact here reports ±5% where the experiment supports ±15%. The two are independent, so they combine in quadrature, each side of the interval separately; asymmetry survives, an undetermined limit stays undetermined, and a lower limit driven past zero is reported as undetermined rather than as a Ki of zero. tracer_conc is taken as exact.

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

The relation also assumes the competitor and the tracer exclude each other from one site, and that the free tracer concentration is close to the total. Neither is something a fit can check, since the curve looks the same either way.

Ligand Depletion

Every model in the table except tight_binding assumes the free ligand concentration equals the total one. That stops holding once the receptor is not much more dilute than Kd, because each molecule bound is one fewer left in solution. The curve still looks like a saturation curve, so nothing in the fit statistics gives the problem away: with the receptor at five times Kd, langmuir reports Kd = 3.99 where the truth is 1.0, at R² = 0.9919.

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")

Left free instead, rt is estimated from the data, which 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 is what keeps rt and Kd 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 limits are worth knowing:

• 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 Kd. Depletion reaches the plateau at a lower multiple of Kd than a hyperbola does, so those thresholds are conservative here and may still ask for higher concentrations than this model needs.

Cooperativity Under Depletion

Depletion does not only shift the curve, it steepens it, and hill reads that steepness as an exponent. The same 1:1 data as above, with no cooperativity anywhere in the system, fits hill at n = 1.43 [1.28, 1.59], R² = 0.9989: the interval excludes 1, so the check documented above (intervals["n"].contains(1.0)) reports cooperativity that does not exist.

There is no model here that solves depletion and cooperativity at once, and adding one would mean bolting an exact conservation law onto the Hill equation, which is phenomenological to begin with. So cooperativity has to be measured where depletion is absent: keep the receptor at or below Kd/10 for that question.

A significant n > 1 is therefore always reported with the alternatives that produce the same shape (depletion, self-association, a reading taken before equilibrium). Passing receptor_conc= lets the first of them be checked rather than only listed; without it the warning says so instead of staying silent.

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 parameters, n/k of 2-5).

This has been confirmed empirically as well: with fewer points, the uncorrected AIC tends to favour the model with more parameters.

shared.aicc < free.aicc
fit_global(ds, model=hill).aicc < fit_global(ds, model=langmuir).aicc

When the sample is too small for the correction to be defined (n − k − 1 ≤ 0), aicc is infinite. report() shows AICc first, with AIC in parentheses for reference.

Diagnostics

It is not unusual for Kd to be undetermined even when R² is above 0.99. The following are warned about automatically. A quantity whose name changes by model (Kd, Km) is shown under that model's own label.

There are two kinds of diagnostics: 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, remarks whose premise is removed by sharing or fixing are suppressed. For example, if bmax is shared, "Kd and Bmax cannot be separated" no longer applies, and a note saying "sharing is what made this estimate possible" appears instead. That note is withheld, though, when the fit itself is broken (R² < 0.5, or a collapsed amplitude), since it would otherwise contradict the other advice.

examples/titration_unsaturated.csv is an example where R² = 0.9936 looks like a good fit, but the highest concentration is only a third of Kd, 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 number is invented; 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 fires if either is significant. The false-positive rate is kept under 2%.

Statistics

The warnings are judged at a fixed significance level. Stacking tests across several datasets raises the false-positive rate accordingly, so the underlying statistic and p-value are returned as-is.

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"]. residual_autocorrelation is judged against a fixed threshold rather than a 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 Kd 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.

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

Running this sample on examples/titration_good.csv: