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 = 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:
Return type:

CurveParams

class pystatsbio.doseresponse.DoseResponseResult(params, se, residuals, rss, aic, bic, converged, n_iter, model, dose, response, n_obs, jac)[source]

Bases: object

Result of fitting a single dose-response curve.

Parameters:
params: CurveParams
se: 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]]
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.BatchDoseResponseResult(ec50, hill, top, bottom, converged, rss, n_compounds)[source]

Bases: object

Result 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
class pystatsbio.doseresponse.EC50Result(estimate, se, ci_lower, ci_upper, conf_level, method)[source]

Bases: object

EC50 (or IC50) with confidence interval.

Parameters:
estimate: float
se: float
ci_lower: float
ci_upper: float
conf_level: float
method: str
class pystatsbio.doseresponse.RelativePotencyResult(ratio, ci_lower, ci_upper, conf_level, method)[source]

Bases: object

Relative potency (ratio of EC50s) with Fieller’s CI.

Parameters:
ratio: float
ci_lower: float
ci_upper: float
conf_level: float
method: str
class pystatsbio.doseresponse.BMDResult(bmd, bmdl, bmdu, bmr, conf_level, method)[source]

Bases: object

Benchmark dose result.

Parameters:
bmd: float
bmdl: float
bmdu: float
bmr: float
conf_level: float
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:
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:
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:

DoseResponseResult

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) – 'cpu', 'gpu', or 'auto'. 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:

BatchDoseResponseResult

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 EC50 with confidence interval from a fitted model.

For models where EC50 is a direct parameter (LL.4, LL.5, W1.4, W2.4, BC.5), the standard error comes from the parameter covariance matrix. The confidence interval is constructed on the log scale (since EC50 is positive) and back-transformed.

Parameters:

  • fit_result (DoseResponseResult) – A fitted dose-response model.

  • conf_level (float) – Confidence level (default 0.95).

  • method (str) – 'delta' (delta method).

Returns:

  • EC50Result

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

Return type:

EC50Result

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

Relative potency: ratio of EC50s between two curves with Fieller’s CI.

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

Parameters:
Returns:

  • RelativePotencyResult

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

Return type:

RelativePotencyResult

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 (DoseResponseResult) – 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:

BMDResult

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