COMPARTMENTAL_PK

Overview

The COMPARTMENTAL_PK function simulates drug concentration over time using a one-compartment pharmacokinetic model with first-order elimination kinetics. This model represents the simplest pharmacokinetic system where the drug distributes instantaneously throughout a single compartment (the body) and is eliminated at a rate proportional to its concentration. The governing differential equation is:

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

where C is the drug concentration, t is time, and k_{el} is the elimination rate constant. The analytical solution to this equation is C(t) = C_0 \cdot e^{-k_{el} \cdot t}, where C_0 is the initial concentration. However, this function uses numerical integration to solve the system, which provides a foundation for extending to more complex multi-compartment models or non-linear elimination processes.

The implementation uses scipy.integrate.solve_ivp, a versatile ordinary differential equation (ODE) solver from the SciPy library. The function supports multiple integration methods including explicit Runge-Kutta methods (RK45, RK23, DOP853) suitable for non-stiff problems, and implicit methods (Radau, BDF, LSODA) for stiff systems. The default RK45 method uses a fifth-order accurate Runge-Kutta formula with adaptive step sizing and error control.

One-compartment models are widely used in early-stage drug development, clinical pharmacology, and therapeutic drug monitoring. The volume of distribution (V_d) parameter relates the administered dose to the resulting plasma concentration, while the elimination rate constant determines the drug’s half-life (t_{1/2} = \ln(2)/k_{el}). For more information on compartmental modeling in pharmacokinetics, see Compartmental models in pharmacokinetics on Wikipedia.

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

Excel Usage

=COMPARTMENTAL_PK(c_initial, k_el, v_d, t_start, t_end, timesteps, solve_ivp_method)

c_initial (float, required): Initial drug concentration (mass/volume).
k_el (float, required): Elimination rate constant (1/time).
v_d (float, required): Volume of distribution (volume).
t_start (float, required): Start time for integration.
t_end (float, required): End time for integration.
timesteps (int, optional, default: 10): Number of output time points.
solve_ivp_method (str, optional, default: “RK45”): Integration method.

Returns (list[list]): 2D list with header [t, C], or error message string.

Examples

Example 1: Demo case 1

Inputs:

c_initial	k_el	v_d	t_start	t_end	timesteps
10	0.2	5	0	10	10

Excel formula:

=COMPARTMENTAL_PK(10, 0.2, 5, 0, 10, 10)

Expected output:

t	C
0	10
1.1111	8.0074
2.2222	6.4106
3.3333	5.1308
4.4444	4.1091
5.5556	3.2933
6.6667	2.6376
7.7778	2.1116
8.8889	1.6911
10	1.3545

Example 2: Demo case 2

Inputs:

c_initial	k_el	v_d	t_start	t_end	timesteps	solve_ivp_method
5	0.1	2	0	5	10	RK23

Excel formula:

=COMPARTMENTAL_PK(5, 0.1, 2, 0, 5, 10, "RK23")

Expected output:

t	C
0	5
0.5556	4.7298
1.1111	4.4739
1.6667	4.2316
2.2222	4.0023
2.7778	3.7854
3.3333	3.5804
3.8889	3.3868
4.4444	3.2038
5	3.0306

Example 3: Demo case 3

Inputs:

c_initial	k_el	v_d	t_start	t_end	timesteps	solve_ivp_method
20	0.05	10	0	20	10	BDF

Excel formula:

=COMPARTMENTAL_PK(20, 0.05, 10, 0, 20, 10, "BDF")

Expected output:

t	C
0	20
2.2222	17.8976
4.4444	16.0144
6.6667	14.3295
8.8889	12.8231
11.1111	11.4754
13.3333	10.2698
15.5556	9.1904
17.7778	8.2245
20	7.3598

Example 4: Demo case 4

Inputs:

c_initial	k_el	v_d	t_start	t_end	timesteps
1	0.5	1	0	1	10

Excel formula:

=COMPARTMENTAL_PK(1, 0.5, 1, 0, 1, 10)

Expected output:

t	C
0	1
0.1111	0.946
0.2222	0.8948
0.3333	0.8465
0.4444	0.8007
0.5556	0.7575
0.6667	0.7165
0.7778	0.6778
0.8889	0.6412
1	0.6065

Python Code

from scipy.integrate import solve_ivp as scipy_solve_ivp
import numpy as np
import math

def compartmental_pk(c_initial, k_el, v_d, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the basic one-compartment pharmacokinetics ODE using scipy.integrate.solve_ivp.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

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

    Args:
        c_initial (float): Initial drug concentration (mass/volume).
        k_el (float): Elimination rate constant (1/time).
        v_d (float): Volume of distribution (volume).
        t_start (float): Start time for integration.
        t_end (float): End time for integration.
        timesteps (int, optional): Number of output time points. Default is 10.
        solve_ivp_method (str, optional): Integration method. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header [t, C], or error message string.
    """
    valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']

    # Validate input types
    try:
        c0 = float(c_initial)
        k = float(k_el)
        vd = 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 <= 0 or vd <= 0:
        return "Invalid input: k_el and v_d must be positive."
    if solve_ivp_method not in valid_methods:
        return f"Invalid input: solve_ivp_method must be one of {valid_methods}."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)

    # ODE system: dC/dt = -k * C (first-order elimination)
    def pk_ode(t, y):
        c = y[0]
        dc_dt = -k * c
        return [dc_dt]

    # Integrate using scipy solve_ivp
    try:
        sol = scipy_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])
        # Validate output values
        if math.isnan(t_val) or math.isnan(c_val):
            return "Invalid output: NaN detected in results."
        if math.isinf(t_val) or math.isinf(c_val):
            return "Invalid output: Infinite value detected in results."
        result.append([t_val, c_val])
    return result

Overview

Excel Usage

Examples

Python Code

Online Calculator