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:

\[\begin{split}\begin{aligned} \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) \end{aligned}\end{split}\]

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

This Python version solves the model over a 24 hour window with multiple bolus dosing events applied via hybrid stop/reset callbacks.

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,
        }
        tdose_i { -1.0, 6.0, 12.0, 18.0 }
        dose_i { 0.0, 1000.0, 2000.0, 1000.0 }
        u_i {
            qc = 1000.0,
            qp1 = 0.0,
        }
        F_i {
            - qc / Vc * CL - Qp1 * (qc / Vc - qp1 / Vp1),
            Qp1 * (qc / Vc - qp1 / Vp1),
        }
        stop_i {
            tdose_i - t,
        }
        reset_i {
            qc + dose_i[N],
            qp1,
        }
        """,
        matrix_type=ds.nalgebra_dense,
        ode_solver=ds.tsit45,
        linear_solver=ds.lu,
    )

    params = np.array([1000.0, 1000.0, 100.0, 50.0])
    solution = ode.solve(params, 24.0)

    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()