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.

Example 1: First-order elimination kinetics over 10 time units

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 8.1873
2 6.70255
3 5.48507
4 4.49026
5 3.67851
6 3.01391
7 2.46728
8 2.01979
9 1.65399
10 1.35453
Example 2: Slow elimination with RK23 integration method

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.5 4.75615
1 4.52401
1.5 4.30292
2 4.09249
2.5 3.89232
3 3.70201
3.5 3.52116
4 3.34937
4.5 3.18601
5 3.03061
Example 3: Long simulation with BDF method for stiff systems

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 18.0977
4 16.3746
6 14.8153
8 13.4055
10 12.1304
12 10.9772
14 9.93315
16 8.98858
18 8.13366
20 7.35977
Example 4: Fast elimination rate over short time interval

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.1 0.951229
0.2 0.904836
0.3 0.860704
0.4 0.818724
0.5 0.778793
0.6 0.740811
0.7 0.704684
0.8 0.670319
0.9 0.63763
1 0.606533

Python Code

Show 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.
    """
    try:
      valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']

      c0 = float(c_initial)
      k = float(k_el)
      vd = float(v_d)
      t0 = float(t_start)
      t1 = float(t_end)
      ntp = int(timesteps)

      if t1 <= t0:
        return "Error: Invalid input: t_end must be greater than t_start."
      if ntp <= 0:
        return "Error: Invalid input: timesteps must be a positive integer."
      if k <= 0 or vd <= 0:
        return "Error: Invalid input: k_el and v_d must be positive."
      if solve_ivp_method not in valid_methods:
        return f"Error: Invalid input: solve_ivp_method must be one of {valid_methods}."

      t_eval = np.linspace(t0, t1, ntp + 1)

      def pk_ode(t, y):
        c = y[0]
        dc_dt = -k * c
        return [dc_dt]

      sol = scipy_solve_ivp(
        pk_ode,
        [t0, t1],
        [c0],
        method=solve_ivp_method,
        dense_output=False,
        t_eval=t_eval
      )
      if not sol.success:
        return f"Error: Integration failed: {sol.message}"

      result = [["t", "C"]]
      for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        c_val = float(sol.y[0][i])
        if math.isnan(t_val) or math.isnan(c_val):
          return "Error: Invalid output: NaN detected in results."
        if math.isinf(t_val) or math.isinf(c_val):
          return "Error: Invalid output: Infinite value detected in results."
        result.append([t_val, c_val])
      return result
    except Exception as e:
      return f"Error: {str(e)}"

Online Calculator

Initial drug concentration (mass/volume).
Elimination rate constant (1/time).
Volume of distribution (volume).
Start time for integration.
End time for integration.
Number of output time points.
Integration method.