Dose-Response Modeling¶

4PL/5PL curve fitting, EC50/IC50 estimation, relative potency, and benchmark dose analysis. GPU-accelerated batch fitting for high-throughput screening (HTS).

Validates against R packages: drc, nplr, BMDS.

Dose-response modeling for preclinical pharmacology.

The workhorse of in vitro assay analysis and toxicology studies. Provides 4PL/5PL curve fitting, EC50/IC50 estimation, relative potency, benchmark dose analysis, and GPU-accelerated batch fitting for HTS campaigns.

Validates against: R packages drc, nplr, BMDS.

class pystatsbio.doseresponse.CurveParams(bottom, top, ec50, hill, asymmetry=None, hormesis=None, model='LL.4')[source]¶

Bases: object

Parameters of a fitted dose-response curve.

For 4PL: bottom + (top - bottom) / (1 + (ec50/x)^hill)

Parameters:

bottom (float)
top (float)
ec50 (float)
hill (float)
asymmetry (float | None)
hormesis (float | None)
model (str)

bottom: float¶

top: float¶

ec50: float¶

hill: float¶

asymmetry: float | None = None¶

hormesis: float | None = None¶

model: str = 'LL.4'¶

predict(dose)[source]¶

Predict response at given dose levels.

Parameters:: dose (ndarray[tuple[Any, ...], dtype[floating]])
Return type:: ndarray[tuple[Any, …], dtype[floating]]

to_array()[source]¶

Return parameter vector in model order (for Jacobian indexing).

Return type:: ndarray[tuple[Any, …], dtype[floating]]

static from_array(params, model)[source]¶

Construct from parameter vector and model name.

Parameters:

params (ndarray[tuple[Any, ...], dtype[floating]])
model (str)

Return type:

CurveParams

class pystatsbio.doseresponse.DoseResponseParams(curve, se, cov, residuals, rss, aic, bic, converged, n_iter, model, dose, response, n_obs, jac)[source]¶

Bases: object

Computed payload of fitting a single dose-response curve.

The fitted-curve value object lives in curve (a CurveParams); the envelope reserves Result.params for this payload, so the field is named curve rather than params. The public DoseResponseSolution re-exposes it as .params to preserve the existing API (fit_result.params.ec50 etc.).

Parameters:

curve (CurveParams)
se (ndarray[tuple[Any, ...], dtype[floating]])
cov (ndarray[tuple[Any, ...], dtype[floating]])
residuals (ndarray[tuple[Any, ...], dtype[floating]])
rss (float)
aic (float)
bic (float)
converged (bool)
n_iter (int)
model (str)
dose (ndarray[tuple[Any, ...], dtype[floating]])
response (ndarray[tuple[Any, ...], dtype[floating]])
n_obs (int)
jac (ndarray[tuple[Any, ...], dtype[floating]])

curve: CurveParams¶

se: ndarray[tuple[Any, ...], dtype[floating]]¶

cov: ndarray[tuple[Any, ...], dtype[floating]]¶

residuals: ndarray[tuple[Any, ...], dtype[floating]]¶

rss: float¶

aic: float¶

bic: float¶

converged: bool¶

n_iter: int¶

model: str¶

dose: ndarray[tuple[Any, ...], dtype[floating]]¶

response: ndarray[tuple[Any, ...], dtype[floating]]¶

n_obs: int¶

jac: ndarray[tuple[Any, ...], dtype[floating]]¶

class pystatsbio.doseresponse.DoseResponseSolution(result)[source]¶

Bases: SolutionReprMixin

Public result of a single dose-response fit — wraps Result[DoseResponseParams].

Exposes every fit output as a read-only property plus the uniform .backend_name / .timing / .warnings / .info metadata and a Jupyter _repr_html_ (via SolutionReprMixin).

.params returns the fitted CurveParams, preserving the public API that callers rely on (fit_result.params.ec50 / .to_array()).

Parameters:: result (Result[DoseResponseParams])

property backend_name: str¶

property timing: dict[str, float] | None¶

property warnings: tuple[str, ...]¶

property info: dict¶

property params: CurveParams¶: The fitted dose-response curve (bottom/top/ec50/hill).

property se: ndarray[tuple[Any, ...], dtype[floating]]¶

property cov: ndarray[tuple[Any, ...], dtype[floating]]¶: Parameter covariance matrix (observed information, matching drc).

property residuals: ndarray[tuple[Any, ...], dtype[floating]]¶

property rss: float¶

property aic: float¶

property bic: float¶

property converged: bool¶

property n_iter: int¶

property model: str¶

property dose: ndarray[tuple[Any, ...], dtype[floating]]¶

property response: ndarray[tuple[Any, ...], dtype[floating]]¶

property n_obs: int¶

property jac: ndarray[tuple[Any, ...], dtype[floating]]¶

predict(dose=None)[source]¶

Predict response. If dose is None, use the fitted dose.

Parameters:: dose (ndarray[tuple[Any, ...], dtype[floating]] | None)
Return type:: ndarray[tuple[Any, …], dtype[floating]]

summary()[source]¶

Human-readable summary, similar to R drc::summary().

Return type:: str

class pystatsbio.doseresponse.BatchDoseResponseParams(ec50, hill, top, bottom, converged, rss, n_compounds)[source]¶

Bases: object

Computed payload of batch-fitting dose-response curves (HTS).

Each array has length n_compounds.

Parameters:

ec50 (ndarray[tuple[Any, ...], dtype[floating]])
hill (ndarray[tuple[Any, ...], dtype[floating]])
top (ndarray[tuple[Any, ...], dtype[floating]])
bottom (ndarray[tuple[Any, ...], dtype[floating]])
converged (ndarray[tuple[Any, ...], dtype[bool]])
rss (ndarray[tuple[Any, ...], dtype[floating]])
n_compounds (int)

ec50: ndarray[tuple[Any, ...], dtype[floating]]¶

hill: ndarray[tuple[Any, ...], dtype[floating]]¶

top: ndarray[tuple[Any, ...], dtype[floating]]¶

bottom: ndarray[tuple[Any, ...], dtype[floating]]¶

converged: ndarray[tuple[Any, ...], dtype[bool]]¶

rss: ndarray[tuple[Any, ...], dtype[floating]]¶

n_compounds: int¶

class pystatsbio.doseresponse.BatchDoseResponseSolution(result)[source]¶

Bases: SolutionReprMixin

Public result of batch dose-response fitting — wraps Result[BatchDoseResponseParams].

Exposes every output as a read-only property plus the uniform .backend_name / .timing / .warnings / .info metadata and a Jupyter _repr_html_ (via SolutionReprMixin).

Parameters:: result (Result[BatchDoseResponseParams])

property backend_name: str¶

property timing: dict[str, float] | None¶

property warnings: tuple[str, ...]¶

property info: dict¶

property ec50: ndarray[tuple[Any, ...], dtype[floating]]¶

property hill: ndarray[tuple[Any, ...], dtype[floating]]¶

property top: ndarray[tuple[Any, ...], dtype[floating]]¶

property bottom: ndarray[tuple[Any, ...], dtype[floating]]¶

property converged: ndarray[tuple[Any, ...], dtype[bool]]¶

property rss: ndarray[tuple[Any, ...], dtype[floating]]¶

property n_compounds: int¶

class pystatsbio.doseresponse.EC50Params(estimate, se, ci_lower, ci_upper, conf_level, method)[source]¶

Bases: object

Computed payload of EC50 (or IC50) estimation with confidence interval.

Parameters:

estimate (float)
se (float)
ci_lower (float)
ci_upper (float)
conf_level (float)
method (str)

estimate: float¶

se: float¶

ci_lower: float¶

ci_upper: float¶

conf_level: float¶

method: str¶

class pystatsbio.doseresponse.EC50Solution(result)[source]¶

Bases: SolutionReprMixin

Public result of EC50 estimation — wraps Result[EC50Params].

Exposes every output as a read-only property plus the uniform .backend_name / .timing / .warnings / .info metadata and a Jupyter _repr_html_ (via SolutionReprMixin).

Parameters:: result (Result[EC50Params])

property backend_name: str¶

property timing: dict[str, float] | None¶

property warnings: tuple[str, ...]¶

property info: dict¶

property estimate: float¶

property se: float¶

property ci_lower: float¶

property ci_upper: float¶

property conf_level: float¶

property method: str¶

class pystatsbio.doseresponse.RelativePotencyParams(ratio, ci_lower, ci_upper, conf_level, method)[source]¶

Bases: object

Computed payload of relative potency (ratio of EC50s) with Fieller’s CI.

Parameters:

ratio (float)
ci_lower (float)
ci_upper (float)
conf_level (float)
method (str)

ratio: float¶

ci_lower: float¶

ci_upper: float¶

conf_level: float¶

method: str¶

class pystatsbio.doseresponse.RelativePotencySolution(result)[source]¶

Bases: SolutionReprMixin

Public result of relative potency — wraps Result[RelativePotencyParams].

Exposes every output as a read-only property plus the uniform .backend_name / .timing / .warnings / .info metadata and a Jupyter _repr_html_ (via SolutionReprMixin).

Parameters:: result (Result[RelativePotencyParams])

property backend_name: str¶

property timing: dict[str, float] | None¶

property warnings: tuple[str, ...]¶

property info: dict¶

property ratio: float¶

property ci_lower: float¶

property ci_upper: float¶

property conf_level: float¶

property method: str¶

class pystatsbio.doseresponse.BMDParams(bmd, bmdl, bmdu, bmr, conf_level, method)[source]¶

Bases: object

Computed payload of benchmark dose analysis.

Parameters:

bmd (float)
bmdl (float)
bmdu (float)
bmr (float)
conf_level (float)
method (str)

bmd: float¶

bmdl: float¶

bmdu: float¶

bmr: float¶

conf_level: float¶

method: str¶

class pystatsbio.doseresponse.BMDSolution(result)[source]¶

Bases: SolutionReprMixin

Public result of benchmark dose analysis — wraps Result[BMDParams].

Exposes every output as a read-only property plus the uniform .backend_name / .timing / .warnings / .info metadata and a Jupyter _repr_html_ (via SolutionReprMixin).

Parameters:: result (Result[BMDParams])

property backend_name: str¶

property timing: dict[str, float] | None¶

property warnings: tuple[str, ...]¶

property info: dict¶

property bmd: float¶

property bmdl: float¶

property bmdu: float¶

property bmr: float¶

property conf_level: float¶

property method: str¶

pystatsbio.doseresponse.ll4(dose, bottom, top, ec50, hill)[source]¶

4-parameter log-logistic (LL.4) model.

\[f(x) = c + \frac{d - c}{1 + \exp\bigl(-b \cdot (\ln x - \ln e)\bigr)}\]

where c = bottom, d = top, e = ec50, b = hill.

Parameters:

dose (array) – Dose (concentration) values. May contain zeros.
bottom (float) – Lower asymptote (response at dose → 0 for hill > 0).
top (float) – Upper asymptote (response at dose → ∞ for hill > 0).
ec50 (float) – Dose producing 50 % of the maximal effect.
hill (float) – Hill slope. Positive for increasing, negative for decreasing.

Returns:

NDArray – Predicted response values.
Validates against (R drc::LL.4())

Return type:

ndarray[tuple[Any, …], dtype[floating]]

pystatsbio.doseresponse.ll5(dose, bottom, top, ec50, hill, asymmetry)[source]¶

5-parameter log-logistic (LL.5) model.

Asymmetric extension of LL.4 with an extra shape parameter f. When asymmetry = 1, this reduces to LL.4.

\[f(x) = c + \frac{d - c}{\bigl(1 + \exp(-b \cdot (\ln x - \ln e))\bigr)^f}\]

Parameters:

dose (array) – Dose (concentration) values.
bottom (float) – Same as ll4().
top (float) – Same as ll4().
ec50 (float) – Same as ll4().
hill (float) – Same as ll4().
asymmetry (float) – Asymmetry parameter (f). 1.0 = symmetric (LL.4).

Returns:

NDArray
Validates against (R drc::LL.5())

Return type:

ndarray[tuple[Any, …], dtype[floating]]

pystatsbio.doseresponse.weibull1(dose, bottom, top, ec50, hill)[source]¶

Weibull type 1 (W1.4) model.

Asymmetric dose-response, left-skewed.

\[f(x) = c + (d - c) \exp\bigl(-\exp(-b \cdot (\ln x - \ln e))\bigr)\]

Parameters:

dose (array) – Dose (concentration) values.
bottom (float) – Same as ll4().
top (float) – Same as ll4().
ec50 (float) – Same as ll4().
hill (float) – Same as ll4().

Returns:

NDArray
Validates against (R drc::W1.4())

Return type:

ndarray[tuple[Any, …], dtype[floating]]

pystatsbio.doseresponse.weibull2(dose, bottom, top, ec50, hill)[source]¶

Weibull type 2 (W2.4) model.

Asymmetric dose-response, right-skewed.

\[f(x) = c + (d - c) \bigl(1 - \exp(-\exp(-b \cdot (\ln x - \ln e)))\bigr)\]

Parameters:

dose (array) – Dose (concentration) values.
bottom (float) – Same as ll4().
top (float) – Same as ll4().
ec50 (float) – Same as ll4().
hill (float) – Same as ll4().

Returns:

NDArray
Validates against (R drc::W2.4())

Return type:

ndarray[tuple[Any, …], dtype[floating]]

pystatsbio.doseresponse.brain_cousens(dose, bottom, top, ec50, hill, hormesis)[source]¶

Brain-Cousens hormesis model (BC.5).

Biphasic dose-response with low-dose stimulation. The hormesis parameter adds a linear term to the numerator that can push the response above the upper asymptote at low doses.

\[f(x) = c + \frac{d - c + f \cdot x}{1 + \exp(-b \cdot (\ln x - \ln e))}\]

Parameters:

dose (array) – Dose (concentration) values.
bottom (float) – Same as ll4().
top (float) – Same as ll4().
ec50 (float) – Same as ll4().
hill (float) – Same as ll4().
hormesis (float) – Hormesis coefficient (f). 0.0 reduces to LL.4.

Returns:

NDArray
Validates against (R drc::BC.5())

Return type:

ndarray[tuple[Any, …], dtype[floating]]

pystatsbio.doseresponse.fit_drm(dose, response, *, model='LL.4', weights=None, start=None, lower=None, upper=None)[source]¶

Fit a dose-response model to a single curve.

Uses Trust Region Reflective nonlinear least squares (scipy.optimize.least_squares).

Parameters:

dose (array) – Dose (concentration) values.
response (array) – Response values.
model (str) – Model name: 'LL.4', 'LL.5', 'W1.4', 'W2.4', 'BC.5'.
weights (array or None) – Optional observation weights.
start (dict or None) – Starting values for parameters. If None, uses self-starting estimates derived from the data.
lower (dict or None) – Box constraints on parameters.
upper (dict or None) – Box constraints on parameters.

Return type:

DoseResponseSolution

Notes

Parameter standard errors (.se) and the covariance (.cov) use the observed-information covariance s²·(½H)⁻¹ — the inverse of the scaled Hessian of the residual sum of squares, s² = RSS/(n-p) — which is what R drc::drm() reports. This differs from the Gauss-Newton approximation s²·(JᵀJ)⁻¹ used by scipy.optimize.curve_fit / R nls by the second-order curvature term; the two are asymptotically equivalent but can differ by ~10% on a poorly-conditioned coefficient (e.g. the Hill slope).

Examples

>>> import numpy as np
>>> dose = np.array([0, 0.01, 0.1, 1, 10, 100])
>>> response = np.array([10, 12, 30, 55, 85, 92])
>>> result = fit_drm(dose, response, model='LL.4')
>>> round(result.params.ec50, 1)
1.0

Validates against: R drc::drm()

pystatsbio.doseresponse.fit_drm_batch(dose_matrix, response_matrix, *, model='LL.4', backend='auto', max_iter=100, tol=1e-08)[source]¶

Batch-fit dose-response curves across many compounds.

Parameters:

dose_matrix (array, shape (n_compounds, n_doses)) – Dose values for each compound.
response_matrix (array, shape (n_compounds, n_doses)) – Response values for each compound.
model (str) – Model name (currently only 'LL.4' for batch fitting).
backend (str) – Execution target (device and precision), per the pystatistics convention: 'cpu' (float64), 'gpu' (float32, CUDA or Apple Silicon/MPS), 'gpu_fp64' (float64, CUDA only), or 'auto' (CUDA float32 if present, else CPU). GPU uses batched Levenberg-Marquardt via PyTorch for massive parallelism.
max_iter (int) – Maximum LM iterations per compound (default 100).
tol (float) – Convergence tolerance on relative RSS change (default 1e-8).

Return type:

BatchDoseResponseSolution

Notes

GPU backend requires pip install pystatsbio[gpu] (PyTorch). On CPU, curves are fit sequentially using scipy.optimize. On GPU, all curves are fit simultaneously using batched Jacobian computation and batched normal equations.

pystatsbio.doseresponse.ec50(fit_result, *, conf_level=0.95, method='delta')[source]¶

Extract the EC50 (ED50) with confidence interval from a fitted model.

The EC50 is the dose producing a response half-way between the fitted lower and upper asymptotes, obtained by solving the fitted curve — matching R drc::ED(type="relative"). This equals the model’s e location parameter for the symmetric LL.4 and differs for the asymmetric models (LL.5, W1.4, W2.4, BC.5). The confidence interval is a raw-scale symmetric Wald interval (estimate ± t·SE) with the SE from the delta method applied to the solved ED50, matching drc::ED(interval="delta").

Parameters:

fit_result (DoseResponseSolution) – A fitted dose-response model.
conf_level (float) – Confidence level (default 0.95).
method (str) – 'delta' (delta method).

Returns:

EC50Solution
Validates against (R drc::ED())

Return type:

EC50Solution

pystatsbio.doseresponse.relative_potency(fit1, fit2, *, conf_level=0.95)[source]¶

Relative potency: ratio of EC50s (ED50s) between two curves, Fieller’s CI.

Computes rho = ED50_2 / ED50_1 with a confidence interval based on Fieller’s theorem for the ratio of two independent estimates.

Both ED50s are obtained by solving each fitted curve for its half-maximal response (the same quantity ec50() returns), not by ratioing the raw e location parameters. For the symmetric LL.4 — and for any two parallel curves, where the e-to-ED50 factor is identical and cancels — the two are the same. They differ for non-parallel asymmetric fits (LL.5/W1.4/W2.4/BC.5), where ratioing e would not be the potency ratio at all.

Note this function does not test parallelism; relative potency is only meaningful for curves of the same shape, and assessing that is the caller’s responsibility.

Parameters:

fit1 (DoseResponseSolution) – First fitted model (reference).
fit2 (DoseResponseSolution) – Second fitted model (test).
conf_level (float) – Confidence level (default 0.95).

Returns:

RelativePotencySolution
Validates against (R drc::compParm(), drc::EDcomp())

Return type:

RelativePotencySolution

pystatsbio.doseresponse.bmd(fit_result, *, bmr=0.1, bmr_type='extra', conf_level=0.95, method='delta')[source]¶

Compute benchmark dose (BMD) with BMDL/BMDU.

Parameters:

fit_result (DoseResponseSolution) – A fitted dose-response model.
bmr (float) – Benchmark response level (default 10 % = 0.10).
bmr_type (str) – 'extra' (extra risk) or 'additional' (additional risk).
conf_level (float) – Confidence level (default 0.95).
method (str) – 'delta' (delta method).

Return type:

BMDSolution

Raises:

ValidationError – If inputs are invalid or the BMR target is outside the curve range.
NumericalError – If numerical BMD computation or CI estimation fails.

Notes

For extra risk: the target response is top - bmr * (top - bottom) (i.e. a bmr fraction of the full range from the upper asymptote).

For additional risk: the target is top - bmr * |top - bottom|.

Validates against: EPA BMDS software, R BMDL packages