Ajuste de curvas dosis-respuesta: cálculo automático de EC50 e intervalos de confianza con numpy

Al finalizar este tema

Con numpy, matplotlib y la optimización voraz aprendidos en los libros de texto, puedes crear directamente una herramienta que calcule automáticamente el EC50 y sus intervalos de confianza a partir de datos de dosis-respuesta en el cribado de fármacos. El tiempo que antes se dedicaba a calcular manualmente uno por uno los resultados de las placas de 96 pocillos se reduce a 30 segundos.

Este texto es un ejemplo genérico con fines educativos. En la práctica, se utilizan herramientas como GraphPad Prism o el paquete drc de R. Aquí se abordan con precisión los principios internos de dichas herramientas.

«¿Esto todo a mano?» — La trampa de los cálculos repetitivos

Imaginemos que están realizando un cribado de nuevos fármacos. De una sola placa de 96 pocillos, se obtienen los siguientes datos:

• Se coloca en cada pocillo diferentes concentraciones de compuestos candidatos (por ejemplo, 0.001, 0.01, 0.1, 1, 10, 100 μM)
• Se mide la viabilidad celular (%) en cada concentración

Lo que quieren saber: EC50 (la concentración que produce la mitad de la respuesta máxima) y la respuesta máxima.

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

Parámetros:

• bottom: Respuesta mínima (en ausencia de fármaco)
• top: Respuesta máxima (con una dosis suficiente)
• EC50: Concentración a la que la respuesta es (top+bottom)/2
• hill: Pendiente (la inclinación de la curva)

En 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)

En una escala logarítmica, la visualización es más sencilla, por lo que, en general, se aplica una transformación logarítmica a la dosis antes de graficarla.

Parte 2: Optimización voraz

¿Cómo se determinan los 4 parámetros? Se busca minimizar el error (residual) entre los datos y la predicción. Método de los mínimos cuadrados.

minimize sum_i (y_i - y_hat_i)^2

No lineal de mínimos cuadrados no tiene una solución analítica y se resuelve mediante algoritmos iterativos, como Levenberg-Marquardt o el método de la región de confianza. Estos algoritmos funcionan con un enfoque voraz: en la posición actual, se mueven ligeramente en la dirección que reduce más el error y, a continuación, repiten el proceso desde esa nueva posición.

scipy.optimize.curve_fit es una implementación de este algoritmo.

from scipy.optimize import curve_fit


def fit_hill(doses: np.ndarray, responses: np.ndarray) -> dict:
    # Estimación de valores iniciales
    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))  # Error estándar

    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]
    }

Componente 3: Visualización de Matplotlib

Gráfico estándar para verificar los resultados del ajuste.

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()

Ejemplo práctico

Comprobación con datos simulados:

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)

Salida esperada:

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

Desvanecimiento: tres espacios en blanco para completar.

Espacio en blanco 1: intervalo de confianza del 95%

La desviación estándar por sí sola no es suficiente. Calcule el intervalo de confianza del 95% utilizando la distribución t.

from scipy import stats


def compute_ci(fit: dict, n_data: int, alpha: float = 0.05) -> dict:
    """
    Calcula los intervalos de confianza del 95 %.
    """
    dof = n_data - 4  # Cuatro parámetros
    t_val = stats.t.ppf(1 - alpha / 2, dof)

    # TODO: para cada parámetro, popt ± t_val * se
    # Devuelve: {"ec50_ci": (low, high), "hill_ci": (low, high), ...}
    pass

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

Espacio en blanco 2: Procesamiento por lotes de 96 pocillos

Ajuste automático para cada compuesto a partir de un archivo de placa (CSV).

def batch_fit_plate(csv_path: str) -> "pd.DataFrame":
    """
    CSV: columns = [compound, dose, response]
    Llama a fit_hill para cada compound y devuelve los resultados como DataFrame.
    """
    import pandas as pd
    df = pd.read_csv(csv_path)

    results = []
    for compound, group in df.groupby("compound"):
        # TODO: extraer doses y responses, y llamar a fit_hill
        # Añadir al resultado junto con el nombre de compound
        pass

    return pd.DataFrame(results)

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

Espacio en blanco 3: Gráfico de comparación de múltiples compuestos

Compara las curvas de varios compuestos en un solo eje.

def plot_multi_compounds(fits: dict, doses_dict: dict) -> None:
    """
    fits: {compound_name: fit_result_dict}
    doses_dict: {compound_name: (doses, responses)}
    Superpone las curvas ajustadas de cada compuesto con colores distintos.
    """
    fig, ax = plt.subplots(figsize=(9, 6))

    # TODO: para cada compuesto, dibujar el diagrama de dispersión y la curva ajustada con un color distinto
    # Mostrar la EC50 de cada compuesto en la leyenda
    pass

Pista: 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}").

Reflexión: diferencias con las herramientas de ajuste de curvas profesionales

Regresión robusta: las herramientas profesionales utilizan métodos de ajuste robustos frente a valores atípicos. Su función curve_fit trata todos los datos por igual. Las alternativas profesionales son la pérdida de Huber o RANSAC.

Ajuste ponderado: las herramientas profesionales tienen en cuenta que cada punto de datos tiene una precisión de medición diferente. Por ejemplo, la incertidumbre relativa es mayor en las respuestas bajas. Puede asignar pesos a cada punto mediante el parámetro sigma.

Selección de modelos: además de la ecuación de Hill, existen varios modelos (logística de 4 parámetros, bifásica, sigmoide Emax, etc.). Las herramientas profesionales comparan varios modelos utilizando AIC/BIC para seleccionar el mejor.

GraphPad Prism: es una herramienta estándar en el campo de la farmacología clínica. Su enfoque de Python es bueno para la automatización y la reproducibilidad, pero Prism tiene una mejor accesibilidad a través de su interfaz gráfica.

Intervalos de confianza bootstrap: los intervalos de confianza basados en la distribución t asumen que los errores siguen una distribución normal. Cuando esta suposición no se cumple, se pueden obtener intervalos de confianza mediante el remuestreo bootstrap. Las herramientas profesionales también utilizan este enfoque con frecuencia.

Proyectos de ampliación

1. Aplicación Streamlit: los usuarios cargan un archivo CSV, se realiza un ajuste y una visualización automáticos, y se genera un informe en PDF.

2. Intervalos de confianza bootstrap: se remuestrean los datos y se realiza el ajuste varias veces para obtener la distribución de EC50 y calcular los intervalos de confianza.

3. Base de datos de respuesta a dosis: se almacenan los resultados de varios compuestos en SQLite y se realiza un análisis de la relación estructura-actividad (SAR) a lo largo del tiempo.

4. Interacción entre fármacos: se implementa el modelo de Bliss/Loewe para determinar si la combinación de dos fármacos es aditiva, sinérgica o antagónica.

Mapa de componentes de este fragmento

• [F] numpy: manipulación de matrices, transformación de escala logarítmica, cálculo de la desviación estándar.
• [F] matplotlib: diagrama de dispersión en escala logarítmica + curva de ajuste + línea de referencia.
• [F] optimización voraz: comprender el algoritmo de Levenberg-Marquardt dentro de scipy.optimize.curve_fit.
• [W] E/S de archivos: análisis de CSV (se proporciona un script completo).

[F] = usted lo implementa / [W] = se proporciona como código completo.