Example: Compartmental Models of Drug Delivery

Pharmacokinetic (PK) models describe how a drug moves through the body using compartments and rate processes (absorption, distribution, metabolism, and excretion). In a two-compartment PK model we track drug amount in:

  • a central compartment, \(q_c\)

  • a peripheral compartment, \(q_{p1}\)

For linear clearance from the central compartment and linear exchange between central and peripheral compartments, the model is:

\[\frac{dq_c}{dt} = - \frac{q_c}{V_c} CL - Q_{p1} \left( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right)\]
\[\frac{dq_{p1}}{dt} = Q_{p1} \left( \frac{q_c}{V_c} - \frac{q_{p1}}{V_{p1}} \right)\]

with initial conditions \(q_c(0)=1000\) ng and \(q_{p1}(0)=0\) ng.

We use:

  • \(V_c = 1000\) mL

  • \(V_{p1} = 1000\) mL

  • \(CL = 100\) mL/h

  • \(Q_{p1} = 50\) mL/h

Bolus doses of 1000 ng are applied at 0, 6, 12, and 18 hours. In this Python version, each bolus event is implemented by updating solution.current_state for the central compartment before continuing the solve. The example uses ode.solve(...) segment-by-segment between dose times.

import numpy as np
import matplotlib.pyplot as plt
import pydiffsol as ds

def solve():
    ode = ds.Ode(
        """
        in_i {
            Vc = 1000.0,
            Vp1 = 1000.0,
            CL = 100.0,
            Qp1 = 50.0,
            dose0 = 1000.0,
        }
        u_i {
            qc = dose0,
            qp1 = 0.0,
        }
        F_i {
            - qc / Vc * CL - Qp1 * (qc / Vc - qp1 / Vp1),
            Qp1 * (qc / Vc - qp1 / Vp1),
        }
        """,
        matrix_type=ds.nalgebra_dense,
        method=ds.tsit45,
        linear_solver=ds.lu,
    )

    # [Vc, Vp1, CL, Qp1, dose0]
    params = np.array([1000.0, 1000.0, 100.0, 50.0, 1000.0])
    # (time [h], bolus amount [ng])
    doses = [(0.0, 1000.0), (6.0, 1000.0), (12.0, 1000.0), (18.0, 1000.0)]
    final_time = 24.0

    solution = None
    for i, (_t_dose, dose) in enumerate(doses):
        if solution:
            # Apply bolus directly to the central compartment.
            y = solution.current_state
            y[0] += dose
            solution.current_state = y

        next_t = doses[i + 1][0] if i + 1 < len(doses) else final_time
        solution = ode.solve(params, next_t, solution=solution)

    ts = solution.ts
    q_c, q_p1 = solution.ys

    fig, ax = plt.subplots()
    ax.plot(ts, q_c, label="q_c (central)")
    ax.plot(ts, q_p1, label="q_p1 (peripheral)")
    ax.set_xlabel("t [h]")
    ax.set_ylabel("amount [ng]")
    ax.legend()
    fig.savefig("docs/images/compartmental_drug_delivery.svg")

if __name__ == "__main__":
    solve()
compartmental_drug_delivery.svg