BINDING_MODEL

Overview

The BINDING_MODEL function fits a variety of ligand binding and sigmoidal dose-response models to experimental data. These models are fundamental in biochemistry, pharmacology, and biophysics for characterizing interactions between molecules such as receptor-ligand binding, enzyme kinetics, and drug-response relationships.

The function supports twelve distinct model types including single-site and two-site competitive binding, heterogeneous binding, rectangular hyperbola saturation, Boltzmann sigmoids, and rational function models. Each model is characterized by specific parameters that describe binding affinities, maximum binding capacities, or curve shape characteristics.

At its core, the function uses scipy.optimize.curve_fit from the SciPy library to perform non-linear least squares regression. The algorithm minimizes the sum of squared residuals between observed data and the model prediction:

\min_{\theta} \sum_{i=1}^{n} \left( y_i - f(x_i, \theta) \right)^2

where y_i represents observed values, x_i the independent variable (typically ligand concentration), and \theta the model parameters to be optimized.

For competitive binding models, the function uses variations of the sigmoidal dose-response equation. For example, the single-site competitive binding model follows:

Y = A_2 + \frac{A_1 - A_2}{1 + 10^{(x - \log x_0)}}

where A_1 and A_2 represent the upper and lower asymptotes, and \log x_0 is the inflection point (IC50 in log scale).

The two-site heterogeneous binding model captures systems with two independent binding sites having different affinities:

Y = \frac{B_{\max1} \cdot x}{K_1 + x} + \frac{B_{\max2} \cdot x}{K_2 + x}

The function returns optimal parameter values along with standard errors derived from the covariance matrix, enabling assessment of fit quality and parameter uncertainty. Appropriate bounds are automatically applied for parameters that have physical constraints (e.g., positive binding constants).

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

Excel Usage

=BINDING_MODEL(xdata, ydata, binding_model_model)

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

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

Examples

Example 1: Demo case 1

Inputs:

binding_model_model	xdata	ydata
single_site_competitive_binding	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=BINDING_MODEL("single_site_competitive_binding", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

A1	A2	LOGx0
2.75	2.75	2.55

Example 2: Demo case 2

Inputs:

binding_model_model	xdata	ydata
two_site_heterogeneous_binding	0.1	0.12996017079405636
	1.3250000000000002	1.1366893577378239
	2.5500000000000003	1.8755994176793738
	3.7750000000000004	2.423717070854732
	5	2.7253987644364344

Excel formula:

=BINDING_MODEL("two_site_heterogeneous_binding", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.12996017079405636;1.1366893577378239;1.8755994176793738;2.423717070854732;2.7253987644364344})

Expected output:

Bmax1	Bmax2	k1	k2
0.6185	4.8072	4.8409	4.8414
12611861.8231	12611857.2897	470646.7607	61351.118

Example 3: Demo case 3

Inputs:

binding_model_model	xdata	ydata
two_site_competitive_binding	0.1	2.75
	1.3250000000000002	2.75
	2.5500000000000003	2.75
	3.7750000000000004	2.75
	5	2.75

Excel formula:

=BINDING_MODEL("two_site_competitive_binding", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

A1	A2	logx01	logx02	f
2.75	2.75	1.325	3.775	0.5

Example 4: Demo case 4

Inputs:

binding_model_model	xdata	ydata
rational_linear_linear	0.1	0.9640376137175587
	1.3250000000000002	0.7322857909033379
	2.5500000000000003	0.692559238727459
	3.7750000000000004	0.6682616941483859
	5	0.6666987751234869

Excel formula:

=BINDING_MODEL("rational_linear_linear", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.9640376137175587;0.7322857909033379;0.692559238727459;0.6682616941483859;0.6666987751234869})

Expected output:

a	b	c
2.288	1.04	1.443
0.3121	0.01144	0.2103

Example 5: Demo case 5

Inputs:

binding_model_model	xdata	ydata
rational_normalized_numerator	0.1	0.41997845148249724
	1.3250000000000002	0.8191125738577802
	2.5500000000000003	1.0182805907520838
	3.7750000000000004	1.1214305702261564
	5	1.2246818548206133

Excel formula:

=BINDING_MODEL("rational_normalized_numerator", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.41997845148249724;0.8191125738577802;1.0182805907520838;1.1214305702261564;1.2246818548206133})

Expected output:

a	b	c
2.699	1.13	1.83
0.1065	0.1519	0.1971

Example 6: Demo case 6

Inputs:

binding_model_model	xdata	ydata
rational_normalized_denom_linear	0.1	0.43852598027486356
	1.3250000000000002	0.8450040808154502
	2.5500000000000003	1.0523220604140058
	3.7750000000000004	1.161481225242571
	5	1.2706464006391858

Excel formula:

=BINDING_MODEL("rational_normalized_denom_linear", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.43852598027486356;0.8450040808154502;1.0523220604140058;1.161481225242571;1.2706464006391858})

Expected output:

a	b	c
2.509	0.9765	1.699
0.2731	0.132	0.05509

Example 7: Demo case 7

Inputs:

binding_model_model	xdata	ydata
rational_normalized_num_linear	0.1	0.3945608923626594
	1.3250000000000002	0.4715012541865205
	2.5500000000000003	0.5010986982196675
	3.7750000000000004	0.513787024342421
	5	0.527072218818205

Excel formula:

=BINDING_MODEL("rational_normalized_num_linear", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.3945608923626594;0.4715012541865205;0.5010986982196675;0.513787024342421;0.527072218818205})

Expected output:

a	b	c
2.505	0.9583	1.766
0.2288	0.09172	0.01611

Example 8: Demo case 8

Inputs:

binding_model_model	xdata	ydata
offset_reciprocal_hyperbola	0.1	2.120849355572422
	1.3250000000000002	2.012780422610218
	2.5500000000000003	1.9515255947447598
	3.7750000000000004	1.9076012828530458
	5	1.8871946444394427

Excel formula:

=BINDING_MODEL("offset_reciprocal_hyperbola", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.120849355572422;2.012780422610218;1.9515255947447598;1.9076012828530458;1.8871946444394427})

Expected output:

a	b	c
3.001	1.198	1.735
0.3336	0.1816	0.01873

Example 9: Demo case 9

Inputs:

binding_model_model	xdata	ydata
rational_linear_quadratic	0.1	2.3552484325266896
	1.3250000000000002	0.3937363619329481
	2.5500000000000003	0.17766610281077172
	3.7750000000000004	0.049408836408718515
	5	0.07684997888239337

Excel formula:

=BINDING_MODEL("rational_linear_quadratic", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.3552484325266896;0.3937363619329481;0.17766610281077172;0.049408836408718515;0.07684997888239337})

Expected output:

a	b	c	d
3.124	-0.04454	3.093	1.508
1.929	2.299	7.982	7.402

Example 10: Demo case 10

Inputs:

binding_model_model	xdata	ydata
rectangular_hyperbola_saturation	0.1	0.27985054044043783
	1.3250000000000002	1.5948135745888388
	2.5500000000000003	2.026382357740335
	3.7750000000000004	2.2528257698782923
	5	2.30126094772332

Excel formula:

=BINDING_MODEL("rectangular_hyperbola_saturation", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {0.27985054044043783;1.5948135745888388;2.026382357740335;2.2528257698782923;2.30126094772332})

Expected output:

a	b
2.775	1.05
0.05635	0.08208

Example 11: Demo case 11

Inputs:

binding_model_model	xdata	ydata
boltzmann_sigmoid	0.01	2.75
	2.0075	2.75
	4.005	2.75
	6.0024999999999995	2.75
	8	2.75

Excel formula:

=BINDING_MODEL("boltzmann_sigmoid", {0.01;2.0075;4.005;6.0024999999999995;8}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

A1	A2	x0	dx
2.75	2.75	4.005	1

Example 12: Demo case 12

Inputs:

binding_model_model	xdata	ydata
double_boltzmann_sigmoid	0.01	1.3816540661021661
	2.0075	1.1030136074857504
	4.005	0.8671042453490405
	6.0024999999999995	0.666721135042897
	8	0.46429170454060964

Excel formula:

=BINDING_MODEL("double_boltzmann_sigmoid", {0.01;2.0075;4.005;6.0024999999999995;8}, {1.3816540661021661;1.1030136074857504;0.8671042453490405;0.666721135042897;0.46429170454060964})

Expected output:

y0	A	frac	x01	x02	k1	k2
0.4089	1.039	0.5536	1.704	6.183	0.8241	0.9074

Python Code

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

def binding_model(xdata, ydata, binding_model_model):
    """
    Fits binding_model 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
        binding_model_model (str): The binding_model_model value Valid options: Single Site Competitive Binding, Two Site Heterogeneous Binding, Two Site Competitive Binding, Rational Linear Linear, Rational Normalized Numerator, Rational Normalized Denom Linear, Rational Normalized Num Linear, Offset Reciprocal Hyperbola, Rational Linear Quadratic, Rectangular Hyperbola Saturation, Boltzmann Sigmoid, Double Boltzmann Sigmoid.

    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 = {
        'single_site_competitive_binding': {
            'params': ['A1', 'A2', 'LOGx0'],
            'model': lambda x, A1, A2, LOGx0: A2 + (A1 - A2) / (1.0 + np.power(10.0, x - LOGx0)),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.min(ya)), float(np.median(xa))),
        },
        'two_site_heterogeneous_binding': {
            'params': ['Bmax1', 'Bmax2', 'k1', 'k2'],
            'model': lambda x, Bmax1, Bmax2, k1, k2: Bmax1 * x / (k1 + x) + Bmax2 * x / (k2 + x),
            'guess': lambda xa, ya: (float(np.max(ya) * 0.6), float(np.max(ya) * 0.4), float(np.percentile(xa, 25)), float(np.percentile(xa, 75))),
            'bounds': (0.0, np.inf),
        },
        'two_site_competitive_binding': {
            'params': ['A1', 'A2', 'logx01', 'logx02', 'f'],
            'model': lambda x, A1, A2, logx01, logx02, f: A2 + (A1 - A2) * (f / (1.0 + np.power(10.0, x - logx01)) + (1.0 - f) / (1.0 + np.power(10.0, x - logx02))),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.min(ya)), float(np.percentile(xa, 25)), float(np.percentile(xa, 75)), 0.5),
            'bounds': ([-np.inf, -np.inf, -np.inf, -np.inf, 0.0], [np.inf, np.inf, np.inf, np.inf, 1.0]),
        },
        'rational_linear_linear': {
            'params': ['a', 'b', 'c'],
            'model': lambda x, a, b, c: (b + c * x) / (1.0 + a * x),
            'guess': lambda xa, ya: (float(min(ya)), 1.0, 0.1),
        },
        'rational_normalized_numerator': {
            'params': ['a', 'b', 'c'],
            'model': lambda x, a, b, c: (1.0 + c * x) / (a + b * x),
            'guess': lambda xa, ya: (1.0, 1.0, 0.1),
        },
        'rational_normalized_denom_linear': {
            'params': ['a', 'b', 'c'],
            'model': lambda x, a, b, c: (b + c * x) / (a + x),
            'guess': lambda xa, ya: (1.0, float(np.min(ya)), 0.1),
        },
        'rational_normalized_num_linear': {
            'params': ['a', 'b', 'c'],
            'model': lambda x, a, b, c: (b + x) / (a + c * x),
            'guess': lambda xa, ya: (1.0, float(np.min(ya)), 0.1),
        },
        'offset_reciprocal_hyperbola': {
            'params': ['a', 'b', 'c'],
            'model': lambda x, a, b, c: c + b / (x + a),
            'guess': lambda xa, ya: (1.0, 1.0, float(np.min(ya))),
        },
        'rational_linear_quadratic': {
            'params': ['a', 'b', 'c', 'd'],
            'model': lambda x, a, b, c, d: (a + b * x) / (1.0 + c * x + d * np.square(x)),
            'guess': lambda xa, ya: (float(np.min(ya)), 0.1, 0.0, 0.0),
        },
        'rectangular_hyperbola_saturation': {
            'params': ['a', 'b'],
            'model': lambda x, a, b: a * b * x / (1.0 + b * x),
            'guess': lambda xa, ya: (float(np.max(ya)), 0.1),
            'bounds': ([-np.inf, 0.0], np.inf),
        },
        'boltzmann_sigmoid': {
            'params': ['A1', 'A2', 'x0', 'dx'],
            'model': lambda x, A1, A2, x0, dx: (A1 - A2) / (1.0 + np.exp((x - x0) / dx)) + A2,
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.min(ya)), float(np.median(xa)), 1.0),
        },
        'double_boltzmann_sigmoid': {
            'params': ['y0', 'A', 'frac', 'x01', 'x02', 'k1', 'k2'],
            'model': lambda x, y0, A, frac, x01, x02, k1, k2: y0 + A * (frac / (1.0 + np.exp((x - x01) / k1)) + (1.0 - frac) / (1.0 + np.exp((x - x02) / k2))),
            'guess': lambda xa, ya: (float(np.min(ya)), float(np.ptp(ya) if np.ptp(ya) else 1.0), 0.5, float(np.percentile(xa, 25)), float(np.percentile(xa, 75)), 1.0, 1.0),
            'bounds': ([-np.inf, -np.inf, 0.0, -np.inf, -np.inf, 0.0, 0.0], [np.inf, np.inf, 1.0, np.inf, np.inf, np.inf, np.inf]),
        }
    }

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

    model_info = models[binding_model_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]

Online Calculator