COMPARTMENTAL_PK

Overview

The COMPARTMENTAL_PK function numerically solves the basic compartmental pharmacokinetics system of ordinary differential equations, modeling the concentration of a drug in a single compartment with first-order elimination. This is commonly used in pharmacokinetics to estimate drug concentration over time after administration. The ODE system is:

dC/dt = -k_{el} \cdot C

where $C$ is the drug concentration, $k_{el}$ is the elimination rate constant, and $v_d$ is the volume of distribution. The function uses SciPy’s ODE solver (scipy.integrate.solve_ivp ) to integrate the system from t_start to t_end. Only the basic linear one-compartment model is exposed; multi-compartment, nonlinear, or time-varying dosing models are not supported. The integration method can be selected from several common solvers, but only the most used options are exposed. This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:


=COMPARTMENTAL_PK(c_initial, k_el, v_d, t_start, t_end, [method])

c_initial (float, required): Initial drug concentration.
k_el (float, required): Elimination rate constant.
v_d (float, required): Volume of distribution.
t_start (float, required): Start time for integration.
t_end (float, required): End time for integration.
method (str, optional, default=‘RK45’): Integration method. One of RK45, RK23, DOP853, Radau, BDF, LSODA.

The function returns a 2D array (table) with columns: t, C, representing time and drug concentration at each step. If the input is invalid or integration fails, a string error message is returned.

Examples

Example 1: Basic Case

Inputs:

c_initial	k_el	v_d	t_start	t_end	method
10.0	0.2	5.0	0.0	10.0	RK45

Excel formula:


=COMPARTMENTAL_PK(10.0, 0.2, 5.0, 0.0, 10.0)

Expected output:

t	C
0.000	10.000
10.000	1.353

Example 2: With Optional Args (RK23 Method)

Inputs:

c_initial	k_el	v_d	t_start	t_end	method
5.0	0.1	2.0	0.0	5.0	RK23

Excel formula:


=COMPARTMENTAL_PK(5.0, 0.1, 2.0, 0.0, 5.0, "RK23")

Expected output:

t	C
0.000	5.000
5.000	3.033

Example 3: All Args Specified (BDF Method)

Inputs:

c_initial	k_el	v_d	t_start	t_end	method
20.0	0.05	10.0	0.0	20.0	BDF

Excel formula:


=COMPARTMENTAL_PK(20.0, 0.05, 10.0, 0.0, 20.0, "BDF")

Expected output:

t	C
0.000	20.000
20.000	7.358

Example 4: Short Interval

Inputs:

c_initial	k_el	v_d	t_start	t_end	method
1.0	0.5	1.0	0.0	1.0	RK45

Excel formula:


=COMPARTMENTAL_PK(1.0, 0.5, 1.0, 0.0, 1.0)

Expected output:

t	C
0.000	1.000
1.000	0.607

Python Code


from typing import List, Union
from scipy.integrate import solve_ivp
 
def compartmental_pk(c_initial: float, k_el: float, v_d: float, t_start: float, t_end: float, method: str = 'RK45') -> Union[List[List[float]], str]:
    """
    Numerically solves the basic compartmental pharmacokinetics system of ordinary differential equations.
 
    Args:
        c_initial: Initial drug concentration.
        k_el: Elimination rate constant.
        v_d: Volume of distribution.
        t_start: Start time of integration.
        t_end: End time of integration.
        method: Integration method ('RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'). Default is 'RK45'.
 
    Returns:
        2D list with header row: t, C. Each row contains time and concentration values, or an error message (str) if input is invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
    # Validate input types
    try:
        c0 = float(c_initial)
        k_el = float(k_el)
        v_d = float(v_d)
        t0 = float(t_start)
        t1 = float(t_end)
    except Exception:
        return "Invalid input: All initial values and parameters must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if k_el <= 0 or v_d <= 0:
        return "Invalid input: k_el and v_d must be positive."
    if method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # ODE system: dC/dt = -k_el * C
    def pk_ode(t, y):
        c = y[0]
        dc_dt = -k_el * c
        return [dc_dt]
    # Integrate only at t0 and t1
    try:
        sol = solve_ivp(
            pk_ode,
            [t0, t1],
            [c0],
            method=method,
            dense_output=False,
            t_eval=[t0, t1]
        )
    except Exception as e:
        return f"Integration error: {e}"
    if not sol.success:
        return f"Integration failed: {sol.message}"
    # Format output: header row then data rows
    result = [["t", "C"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        c_val = float(sol.y[0][i])
        # Disallow nan/inf
        if any([
            not isinstance(t_val, float) or not isinstance(c_val, float),
            t_val != t_val or c_val != c_val,  # NaN check
            abs(t_val) == float('inf') or abs(c_val) == float('inf')
        ]):
            return "Invalid output: nan or inf detected."
        result.append([t_val, c_val])
    return result

Example Workbook

Link to Workbook