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, [timesteps], [solve_ivp_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.
timesteps (int, optional, default=10): Number of timesteps to solve for.
solve_ivp_method (string (enum), optional, default=‘RK45’): Integration method. Valid options: 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	timesteps	method
10.0	0.2	5.0	0.0	10.0	10	RK45

Excel formula:


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

Expected output:

t	C
0.000	10.000
1.111	8.007
2.222	6.411
3.333	5.131
4.444	4.109
5.556	3.293
6.667	2.638
7.778	2.112
8.889	1.691
10.000	1.355

Example 2: With Optional Args (RK23 Method)

Inputs:

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

Excel formula:


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

Expected output:

t	C
0.000	5.000
0.556	4.730
1.111	4.474
1.667	4.232
2.222	4.002
2.778	3.785
3.333	3.580
3.889	3.387
4.444	3.204
5.000	3.031

Example 3: All Args Specified (BDF Method)

Inputs:

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

Excel formula:


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

Expected output:

t	C
0.000	20.000
2.222	17.898
4.444	16.014
6.667	14.330
8.889	12.823
11.111	11.475
13.333	10.270
15.556	9.190
17.778	8.225
20.000	7.360

Example 4: Short Interval

Inputs:

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

Excel formula:


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

Expected output:

t	C
0.000	1.000
0.111	0.946
0.222	0.895
0.333	0.847
0.444	0.801
0.556	0.758
0.667	0.717
0.778	0.678
0.889	0.641
1.000	0.607

Python Code


from scipy.integrate import solve_ivp
import numpy as np
 
def compartmental_pk(c_initial, k_el, v_d, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    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.
        timesteps: Number of timesteps to solve for. Default is 10.
        solve_ivp_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)
        ntp = int(timesteps)
    except Exception:
        return "Invalid input: All initial values, parameters, and timesteps must be numbers."
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."
    if k_el <= 0 or v_d <= 0:
        return "Invalid input: k_el and v_d must be positive."
    if solve_ivp_method not in ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']:
        return "Invalid input: solve_ivp_method must be one of 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA'."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)
    
    # ODE system: dC/dt = -k_el * C
    def pk_ode(t, y):
        c = y[0]
        dc_dt = -k_el * c
        return [dc_dt]
    # Integrate with timesteps
    try:
        sol = solve_ivp(
            pk_ode,
            [t0, t1],
            [c0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    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