DOSE_RESPONSE

Overview

The DOSE_RESPONSE function fits a variety of dose-response models to experimental data, commonly used in pharmacology, toxicology, and biochemistry to characterize how biological systems respond to varying concentrations of drugs, chemicals, or other stimuli. A dose-response relationship describes the magnitude of a response as a function of exposure to a stimulus, typically producing a sigmoidal curve when plotted on a logarithmic dose scale.

This implementation uses scipy.optimize.curve_fit from the SciPy library to perform non-linear least squares fitting. The function supports ten different models:

Standard Dose Response Curve: Classic sigmoidal model with asymptotes, midpoint, and slope parameters
Hill Cooperative Binding: Models cooperative binding based on the Hill equation
Hill Binding With Baseline: Hill equation extended with start and end baseline values
Four Parameter Logistic (4PL): Standard logistic model widely used for immunoassay analysis
Five Parameter Logistic (5PL): Adds asymmetry parameter to handle non-symmetric curves
Biphasic Dose Response: Models systems with two distinct response phases
Seven Parameter Biphasic Sigmoid: Complex biphasic model with independent phase parameters
Biphasic Hill Activation Inhibition: Models activation at low doses and inhibition at high doses
Generalized Hyperbola: Flexible hyperbolic response model
Single Site Binding Saturation: Simple Michaelis-Menten-like saturation kinetics

The core Hill equation, which underlies several of these models, is expressed as:

E = \frac{E_{max} \cdot [A]^n}{EC_{50}^n + [A]^n}

where E is the response, [A] is the drug concentration, EC_{50} is the concentration producing 50% maximal response, and n is the Hill coefficient describing cooperativity. The Hill coefficient indicates the steepness of the curve—values greater than 1 indicate positive cooperativity, while values less than 1 suggest negative cooperativity.

The function returns parameter names, fitted values, and standard errors derived from the covariance matrix. These metrics are essential for determining potency measures (EC50, IC50) and efficacy parameters commonly reported in drug development and toxicological assessment.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=DOSE_RESPONSE(xdata, ydata, dose_response_model)

xdata (list[list], required): The xdata value
ydata (list[list], required): The ydata value
dose_response_model (str, required): The dose_response_model value

Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error string.

Examples

Example 1: Demo case 1

Inputs:

dose_response_model	xdata	ydata
standard_dose_response_curve	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=DOSE_RESPONSE("standard_dose_response_curve", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

"non-error"

Example 2: Demo case 2

Inputs:

dose_response_model	xdata	ydata
hill_cooperative_binding	0.1	0.01273095510527576
	1.3250000000000002	0.01
	2.5500000000000003	0.10116243219286802
	3.7750000000000004	0.5533417003379951
	5	1.3344595337609784

Excel formula:

=DOSE_RESPONSE("hill_cooperative_binding", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.01273095510527576;0.01;0.10116243219286802;0.5533417003379951;1.3344595337609784})

Expected output:

"non-error"

Example 3: Demo case 3

Inputs:

dose_response_model	xdata	ydata
hill_binding_with_baseline	0.1	2.7606477079062306
	1.3250000000000002	2.7497078533275006
	2.5500000000000003	2.834608579652217
	3.7750000000000004	3.2127948766463232
	5	3.8660934282364545

Excel formula:

=DOSE_RESPONSE("hill_binding_with_baseline", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.7606477079062306;2.7497078533275006;2.834608579652217;3.2127948766463232;3.8660934282364545})

Expected output:

"non-error"

Example 4: Demo case 4

Inputs:

dose_response_model	xdata	ydata
four_parameter_logistic	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=DOSE_RESPONSE("four_parameter_logistic", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

"non-error"

Example 5: Demo case 5

Inputs:

dose_response_model	xdata	ydata
five_parameter_logistic_asymmetric	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=DOSE_RESPONSE("five_parameter_logistic_asymmetric", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

"non-error"

Example 6: Demo case 6

Inputs:

dose_response_model	xdata	ydata
biphasic_dose_response_dual_effect	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=DOSE_RESPONSE("biphasic_dose_response_dual_effect", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

"non-error"

Example 7: Demo case 7

Inputs:

dose_response_model	xdata	ydata
seven_param_biphasic_sigmoid	0.1	2.1495589080728332
	0.8	1.9388050696676988
	1.5	1.3007784892461107
	2.2	0.6666853174024305
	2.9	0.4872847952010856
	3.6	0.6468199412679169
	4.3	1.5337915918649858
	5	2.877898193217767
	5.7	3.391432399977618

Excel formula:

=DOSE_RESPONSE("seven_param_biphasic_sigmoid", {0.1;0.8;1.5;2.2;2.9;3.6;4.3;5;5.7}, {2.1495589080728332;1.9388050696676988;1.3007784892461107;0.6666853174024305;0.4872847952010856;0.6468199412679169;1.5337915918649858;2.877898193217767;3.391432399977618})

Expected output:

"non-error"

Example 8: Demo case 8

Inputs:

dose_response_model	xdata	ydata
biphasic_hill_activation_inhibition	0.1	0.027563250575314988
	1.3250000000000002	0.5895242843046955
	2.5500000000000003	2.229390151955633
	3.7750000000000004	2.9772003029928795
	5	2.1565911679591983

Excel formula:

=DOSE_RESPONSE("biphasic_hill_activation_inhibition", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.027563250575314988;0.5895242843046955;2.229390151955633;2.9772003029928795;2.1565911679591983})

Expected output:

"non-error"

Example 9: Demo case 9

Inputs:

dose_response_model	xdata	ydata
generalized_hyperbola	0.1	1.736625907103898
	1.3250000000000002	1.9296383019658383
	2.5500000000000003	2.004903695261842
	3.7750000000000004	2.042460009056909
	5	2.083675785134604

Excel formula:

=DOSE_RESPONSE("generalized_hyperbola", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {1.736625907103898;1.9296383019658383;2.004903695261842;2.042460009056909;2.083675785134604})

Expected output:

"non-error"

Example 10: Demo case 10

Inputs:

dose_response_model	xdata	ydata
single_site_binding_saturation	0.1	0.06498008539702818
	1.3250000000000002	0.5683446788689119
	2.5500000000000003	0.9377997088396869
	3.7750000000000004	1.211858535427366
	5	1.3626993822182172

Excel formula:

=DOSE_RESPONSE("single_site_binding_saturation", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.06498008539702818;0.5683446788689119;0.9377997088396869;1.211858535427366;1.3626993822182172})

Expected output:

"non-error"

Python Code

import numpy as np
from scipy.optimize import curve_fit as scipy_curve_fit
import math

def dose_response(xdata, ydata, dose_response_model):
    """
    Fits dose_response models to data using scipy.optimize.curve_fit. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html for details.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        xdata (list[list]): The xdata value
        ydata (list[list]): The ydata value
        dose_response_model (str): The dose_response_model value Valid options: Standard Dose Response Curve, Hill Cooperative Binding, Hill Binding With Baseline, Four Parameter Logistic, Five Parameter Logistic Asymmetric, Biphasic Dose Response Dual Effect, Seven Param Biphasic Sigmoid, Biphasic Hill Activation Inhibition, Generalized Hyperbola, Single Site Binding Saturation.

    Returns:
        list[list]: 2D list [param_names, fitted_values, std_errors], or error string.
    """
    def _validate_data(xdata, ydata):
        """Validate and convert both xdata and ydata to numpy arrays."""
        for name, arg in [("xdata", xdata), ("ydata", ydata)]:
            if not isinstance(arg, list) or len(arg) < 2:
                raise ValueError(f"{name}: must be a 2D list with at least two rows")
            vals = []
            for i, row in enumerate(arg):
                if not isinstance(row, list) or len(row) == 0:
                    raise ValueError(f"{name} row {i}: must be a non-empty list")
                try:
                    vals.append(float(row[0]))
                except Exception:
                    raise ValueError(f"{name} row {i}: non-numeric value")
            if name == "xdata":
                x_arr = np.asarray(vals, dtype=np.float64)
            else:
                y_arr = np.asarray(vals, dtype=np.float64)

        if x_arr.shape[0] != y_arr.shape[0]:
            raise ValueError("xdata and ydata must have the same number of rows")
        return x_arr, y_arr

    # Model definitions dictionary
    models = {
        'standard_dose_response_curve': {
            'params': ['A1', 'A2', 'LOGx0', 'p'],
            'model': lambda x, A1, A2, LOGx0, p: A1 + (A2 - A1) / (1.0 + np.power(10.0, (LOGx0 - x) * p)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.max(ya)), float(np.median(xa)), 1.0),
        },
        'hill_cooperative_binding': {
            'params': ['Vmax', 'k', 'n'],
            'model': lambda x, Vmax, k, n: Vmax * np.power(x, n) / (np.power(k, n) + np.power(x, n)),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.median(xa)), 1.0),
            'bounds': (0.0, np.inf),
        },
        'hill_binding_with_baseline': {
            'params': ['START', 'END', 'k', 'n'],
            'model': lambda x, START, END, k, n: START + (END - START) * np.power(x, n) / (np.power(k, n) + np.power(x, n)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.max(ya)), float(np.median(xa)), 1.0),
            'bounds': ([-np.inf, -np.inf, 0.0, 0.0], np.inf),
        },
        'four_parameter_logistic': {
            'params': ['A1', 'A2', 'x0', 'p'],
            'model': lambda x, A1, A2, x0, p: (A1 - A2) / (1.0 + np.power(x / x0, p)) + A2,
            'guess': lambda xa, ya: (float(max(ya)), float(min(ya)), float(np.median(xa)), 1.0),
            'bounds': ([-np.inf, -np.inf, -np.inf, 0.0], np.inf),
        },
        'five_parameter_logistic_asymmetric': {
            'params': ['Amin', 'Amax', 'x0', 'h', 's'],
            'model': lambda x, Amin, Amax, x0, h, s: Amin + (Amax - Amin) / np.power(1.0 + np.power(x / x0, -h), s),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.max(ya)), float(np.median(xa)), 1.0, 1.0),
            'bounds': ([-np.inf, -np.inf, -np.inf, -np.inf, 0.0], np.inf),
        },
        'biphasic_dose_response_dual_effect': {
            'params': ['A1', 'A2', 'LOGx01', 'LOGx02', 'h1', 'h2', 'p'],
            'model': lambda x, A1, A2, LOGx01, LOGx02, h1, h2, p: A1 + (A2 - A1) * (p / (1.0 + np.power(10.0, (LOGx01 - x) * h1)) + (1.0 - p) / (1.0 + np.power(10.0, (LOGx02 - x) * h2))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.max(ya)), float(np.median(xa)), float(np.median(xa)), 1.0, 1.0, 0.5),
            'bounds': ([-np.inf, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf, 0.0], [np.inf, np.inf, np.inf, np.inf, np.inf, np.inf, 1.0]),
        },
        'seven_param_biphasic_sigmoid': {
            'params': ['Amin', 'Amax1', 'Amax2', 'x0_1', 'x0_2', 'h1', 'h2'],
            'model': lambda x, Amin, Amax1, Amax2, x0_1, x0_2, h1, h2: Amin + (Amax1 - Amin) / (1.0 + np.power(10.0, (x - x0_1) * h1)) + (Amax2 - Amin) / (1.0 + np.power(10.0, (x0_2 - x) * h2)),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.max(ya)), float(np.max(ya)), float(np.percentile(xa, 25)), float(np.percentile(xa, 75)), 1.0, 1.0),
        },
        'biphasic_hill_activation_inhibition': {
            'params': ['Pm', 'Ka', 'Ki', 'Ha', 'Hi'],
            'model': lambda x, Pm, Ka, Ki, Ha, Hi: Pm / ((1.0 + np.power(Ka / x, Ha)) * (1.0 + np.power(x / Ki, Hi))),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.median(xa)), float(np.median(xa)), 1.0, 1.0),
            'bounds': (0.0, np.inf),
        },
        'generalized_hyperbola': {
            'params': ['a', 'b', 'c', 'd'],
            'model': lambda x, a, b, c, d: a - b / np.power(1.0 + c * x, 1.0 / d),
            'guess': lambda xa, ya: (float(np.max(ya)), 1.0, 0.1, 1.0),
        },
        'single_site_binding_saturation': {
            'params': ['Bmax', 'K1'],
            'model': lambda x, Bmax, K1: Bmax * x / (K1 + x),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.median(xa))),
            'bounds': (0.0, np.inf),
        }
    }

    # Validate model parameter
    if dose_response_model not in models:
        return f"Invalid model: {str(dose_response_model)}. Valid models are: {', '.join(models.keys())}"

    model_info = models[dose_response_model]

    # Validate and convert input data
    try:
        x_arr, y_arr = _validate_data(xdata, ydata)
    except ValueError as e:
        return f"Invalid input: {e}"

    # Perform curve fitting
    try:
        p0 = model_info['guess'](x_arr, y_arr)
        bounds = model_info.get('bounds', (-np.inf, np.inf))
        if bounds == (-np.inf, np.inf):
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, maxfev=10000)
        else:
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, bounds=bounds, maxfev=10000)

        fitted_vals = [float(v) for v in popt]
        for v in fitted_vals:
            if math.isnan(v) or math.isinf(v):
                return "Fitting produced invalid numeric values (NaN or inf)."
    except ValueError as e:
        return f"Initial guess error: {e}"
    except Exception as e:
        return f"curve_fit error: {e}"

    # Calculate standard errors
    std_errors = None
    try:
        if pcov is not None and np.isfinite(pcov).all():
            std_errors = [float(v) for v in np.sqrt(np.diag(pcov))]
    except Exception:
        pass

    return [model_info['params'], fitted_vals, std_errors] if std_errors else [model_info['params'], fitted_vals]

Overview

Excel Usage

Examples

Python Code

Online Calculator