Dose-Response Curve Fitting — Automatically Calculate EC50 and Confidence Intervals with NumPy

After Completing This Topic

You will be able to create your own tool that combines the numpy, matplotlib, and greedy optimization concepts learned in the textbook to automatically calculate EC50 and confidence intervals from dose-response data in drug screening. This will reduce the time spent manually calculating the results of a 96-well plate from one by one to just 30 seconds.

This article is an educational, general example. In practice, tools like GraphPad Prism or the drc package in R are used. This section focuses on accurately explaining the internal principles of those tools.

"Are you really doing all of this manually?" – The Pitfalls of Repetitive Calculations

Let's say you're performing a drug screening. You obtain the following data from a 96-well plate:

• Each well contains a different concentration of a candidate compound (e.g., 0.001, 0.01, 0.1, 1, 10, 100 μM)
• Cell viability (%) is measured at each concentration

What you want to know: EC50 (the concentration that causes half of the response) and maximum response.

The most straightforward approach: Paste the data into Excel, create a scatter plot, and visually identify the midpoint.

The practical problems with this approach:

Problem 1: Accuracy. The curve drawn by eye will vary from person to person. No reproducibility.

Problem 2: Scale. 100 candidate compounds × 8 concentrations each = 800 data points. It would take all day to click through them in Excel.

Problem 3: Confidence intervals. When the EC50 is, for example, 2.3 μM, you can't know "how accurate is this value, really?".

The real approach is non-linear curve fitting. Find the parameters that best fit the Hill equation (sigmoid function) to the data, and then calculate the statistical confidence interval for those parameters.

From Black Box to Components

Component 1: Hill Equation

The standard mathematical model for a dose-response curve is the Hill equation.

response(dose) = bottom + (top - bottom) / (1 + (EC50 / dose)^hill)

It has four parameters:

• bottom: The minimum response (when no drug is present)
• top: The maximum response (when a large amount of drug is present)
• EC50: The concentration at which the response is (top+bottom)/2
• hill: The slope (steepness of the curve)

In Python:

import numpy as np


def hill_equation(dose: np.ndarray, bottom: float, top: float, ec50: float, hill: float) -> np.ndarray:
    return bottom + (top - bottom) / (1 + (ec50 / dose) ** hill)

Because it is conveniently visualized on a log scale, the dose is generally log-transformed and then plotted.

Component 2: Greedy Optimization

How do we find the four parameters? By minimizing the error (residual) between the data and the prediction. This is least squares.

minimize sum_i (y_i - y_hat_i)^2

Non-linear least squares does not have a closed-form solution and is solved using an iterative algorithm, such as Levenberg-Marquardt or trust-region. These algorithms operate using a greedy approach: they move a small distance in the direction that most reduces the error from the current position, and then repeat from that position.

scipy.optimize.curve_fit is an implementation of this algorithm.

from scipy.optimize import curve_fit


def fit_hill(doses: np.ndarray, responses: np.ndarray) -> dict:
    # Guess initial values
    p0 = [
        responses.min(),  # bottom
        responses.max(),  # top
        np.median(doses), # EC50
        1.0               # hill
    ]
    
    popt, pcov = curve_fit(hill_equation, doses, responses, p0=p0, maxfev=5000)
    perr = np.sqrt(np.diag(pcov))  # Standard error
    
    return {
        "bottom": popt[0],
        "top": popt[1],
        "ec50": popt[2],
        "hill": popt[3],
        "bottom_se": perr[0],
        "top_se": perr[1],
        "ec50_se": perr[2],
        "hill_se": perr[3]
    }

Component 3: Matplotlib Visualization

A standard graph for checking the fitting results.

import matplotlib.pyplot as plt


def plot_dose_response(doses: np.ndarray, responses: np.ndarray, fit: dict) -> None:
    fig, ax = plt.subplots(figsize=(7, 5))
    
    ax.scatter(doses, responses, color="steelblue", s=50, alpha=0.7, label="Data")
    
    dose_grid = np.logspace(np.log10(doses.min() * 0.5), np.log10(doses.max() * 2), 200)
    fit_curve = hill_equation(dose_grid, fit["bottom"], fit["top"], fit["ec50"], fit["hill"])
    ax.plot(dose_grid, fit_curve, color="darkred", linewidth=2, label=f"Fit (EC50={fit['ec50']:.3g})")
    
    ax.axvline(fit["ec50"], linestyle="--", color="gray", alpha=0.5)
    ax.axhline((fit["top"] + fit["bottom"]) / 2, linestyle="--", color="gray", alpha=0.5)
    
    ax.set_xscale("log")
    ax.set_xlabel("Dose (μM)")
    ax.set_ylabel("Response (%)")
    ax.legend()
    plt.tight_layout()
    plt.show()

Practical Example

Let's verify with virtual data:

doses = np.array([0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100])
responses = np.array([2, 3, 5, 12, 28, 55, 78, 90, 96, 98, 99])

fit = fit_hill(doses, responses)
print(f"EC50: {fit['ec50']:.3g} ± {fit['ec50_se']:.3g}")
print(f"Hill: {fit['hill']:.2f}")
print(f"Range: {fit['bottom']:.1f} → {fit['top']:.1f}")

plot_dose_response(doses, responses, fit)

Expected output:

EC50: 0.271 ± 0.0084
Hill: 0.98
Range: 1.2 → 99.5

Fading — Three Blanks for You to Fill

Blank 1: 95% Confidence Interval

Standard error alone is not enough. Calculate the 95% confidence interval using the t-distribution.

from scipy import stats


def compute_ci(fit: dict, n_data: int, alpha: float = 0.05) -> dict:
    """
    Calculate the 95% confidence interval.
    """
    dof = n_data - 4  # 4 parameters
    t_val = stats.t.ppf(1 - alpha / 2, dof)
    
    # TODO: For each parameter, calculate popt ± t_val * se
    # Return: {"ec50_ci": (low, high), "hill_ci": (low, high), ...}
    pass

Hint: ec50_low = fit["ec50"] - t_val * fit["ec50_se"], ec50_high = fit["ec50"] + t_val * fit["ec50_se"].

Blank 2: 96-Well Plate Processing

Automatically fit each compound from a plate file (CSV).

def batch_fit_plate(csv_path: str) -> "pd.DataFrame":
    """
    CSV: columns = [compound, dose, response]
    Call fit_hill for each compound, return the results as a DataFrame.
    """
    import pandas as pd
    df = pd.read_csv(csv_path)
    
    results = []
    for compound, group in df.groupby("compound"):
        # TODO: Extract doses and responses, call fit_hill
        # Append the results to the list, along with the compound name
        pass
    
    return pd.DataFrame(results)

Hint: doses = group["dose"].values; responses = group["response"].values; fit = fit_hill(doses, responses); results.append({"compound": compound, **fit}).

Blank 3: Multi-Compound Comparison Plot

Compare the curves of multiple compounds on the same axis.

def plot_multi_compounds(fits: dict, doses_dict: dict) -> None:
    """
    fits: {compound_name: fit_result_dict}
    doses_dict: {compound_name: (doses, responses)}
    Plot the fitted curves of each compound in different colors, overlapping them.
    """
    fig, ax = plt.subplots(figsize=(9, 6))
    
    # TODO: For each compound, plot the scatter plot and fitted curve in different colors
    # Display the EC50 of each compound in the legend
    pass

Hint: colors = plt.cm.tab10.colors; for i, (name, fit) in enumerate(fits.items()): ax.scatter(..., color=colors[i]); ax.plot(..., color=colors[i], label=f"{name} EC50={fit['ec50']:.3g}").

Reflections — Differences from a Practical Curve Fitting Tool

Robust regression: Practical tools employ robust fitting methods that are resistant to outliers. Your curve_fit treats all data points equally. Huber loss or RANSAC are practical alternatives.

Weighted fitting: Practical tools account for varying measurement errors in each data point. For example, the relative error is larger at low responses. You can assign weights to each point using the sigma parameter.

Model selection: In addition to the Hill equation, there are several other models (4-parameter logistic, biphasic, sigmoid Emax, etc.). Practical tools compare multiple models using AIC/BIC to select the optimal one.

GraphPad Prism: A standard tool in the field of clinical pharmacology. While your Python approach excels in automation and reproducibility, Prism has an advantage in GUI accessibility.

Bootstrap confidence intervals: Confidence intervals based on the t-distribution assume that the errors are normally distributed. When this assumption is violated, confidence intervals can be obtained using bootstrap resampling. Practical tools often use this approach.

Expansion Project

1. Streamlit App: Upon CSV upload by the user, automatically perform fitting, visualization, and generate a PDF report.

2. Bootstrap Confidence Interval: Resample the data and perform multiple fittings to obtain the EC50 distribution and calculate the confidence interval.

3. Dose-Response DB: Accumulate results for multiple compounds in SQLite and perform time-based SAR (structure-activity relationship) analysis.

4. Drug Combination Interaction: Implement the Bliss/Loewe model to determine addition, synergy, or antagonism from combinations of two drugs.

Component Guide for This Module

• [F] numpy: Array manipulation, log-scale transformation, standard error calculation.
• [F] matplotlib: Log-scale scatter plot + fitting curve + reference line.
• [F] greedy optimization: Understanding the Levenberg-Marquardt algorithm within scipy.optimize.curve_fit.
• [W] File I/O: CSV parsing (complete script provided).

[F] = You will implement this yourself / [W] = Complete code will be provided.