Example: Dosing Protocols

In compartmental models of drug delivery we introduced a simple two-compartment model of drug delivery. The model consists of a central compartment and a peripheral compartment. Transitions between the two compartments and the elimination of the drug from the central compartment are governed by the following equations:

\[ \frac{dq_c}{dt} = \text{Dose}(t) - \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 ). \]

We want to use a hybrid model with stopping conditions and a reset condition to model a three dose protocol where a separate dose level is administered at specific times. To do this, we can encode the stopping conditions as follows:

\[ \mathbf{s}(\mathbf{y}) = \begin{bmatrix} 1 \\ t - t_2 \\ t - t_3 \end{bmatrix} \]

where \(t_1, t_2, t_3\) are the times at which doses are administered. Note that we start the model at \(t_1\). The first stopping condition is never satisfied because we never want to transition to \(N=0\) (i.e. this is the starting value), but the other stopping conditions will be satisfied at the appropriate times. The reset condition can then be defined as follows:

\[ \mathbf{r}(\mathbf{y}) = \begin{bmatrix} q_c + \mathbf{d_N} \\ q_{p1} \end{bmatrix} \]

where \(\mathbf{d}\) is a vector of length 3, and \(\mathbf{d}_N\) is the \(N\)th dose level. This means that when the stopping condition is satisfied, the amount of drug in the central compartment is increased by the appropriate dose level, while the amount of drug in the peripheral compartment remains unchanged.

Finally we give the initial conditions for the model based on the first dose level as follows:

\[ \mathbf{y}(0) = \begin{bmatrix}\mathbf{d}_1 \\ 0 \end{bmatrix} \]

To implement this in diffsol we will use the diffsl language to write down the ODE equations, stop conditions and reset functions. Then we will use the Solution and solve_soln methods to stop through the N segments of the hybrid model.

Solution is a struct for iteratively building up a solution from multiple segments, and solve_soln is the method that (a) runs a solve for each segment, and (b) appends the solve result to the Solution arg. After each segment we move out the internal state of the solver. This drops the solver, which contains a reference to the problem, so that we can then mutate the problem by setting the new model index based on the index of the root found in the previous call to solve_soln. Then we can apply our reset condition to the state, and create another solver with this state in order to continue the iteration.

Note that we could have mutated other aspects of the problem/solver at each stage if we wished, setting parameters, changing solver methods or tolerances, based on how we wanted to implement our hybrid model. In this case our only rule is that N is to be set to the index of the previous found root, before the reset is applied.

use diffsol::{
    CraneliftJitModule, MatrixCommon, OdeBuilder, OdeEquations, OdeSolverMethod, OdeSolverState,
    OdeSolverStopReason, Solution,
};
use plotly::{common::Mode, layout::Axis, layout::Layout, Plot, Scatter};
use std::{fs, path::PathBuf};

type M = diffsol::NalgebraMat<f64>;
type CG = CraneliftJitModule;
type LS = diffsol::NalgebraLU<f64>;

fn main() {
    let model = "
        Vc { 1000.0 } Vp1 { 1000.0 } CL { 100.0 } Qp1 { 50.0 }
        doses_i { 1000.0, 500.0, 2000.0}
        u_i {
            qc = doses_i[0],
            qp1 = 0,
        }
        F_i {
            -qc / Vc * CL - Qp1 * (qc / Vc - qp1 / Vp1),
            Qp1 * (qc / Vc - qp1 / Vp1),
        }
        stop_i {
            1.0,
            t - 6.0,
            t - 12.0,
        }
        reset_i {
            qc + doses_i[N],
            qp1,
        }
    ";
    let mut problem = OdeBuilder::<M>::new()
        .build_from_diffsl::<CG>(model)
        .unwrap();

    let final_time = 24.0;
    let mut soln = Solution::new(final_time);
    let mut state = problem.bdf_state::<LS>().unwrap();

    while !soln.is_complete() {
        state = problem
            .bdf_solver::<LS>(state)
            .unwrap()
            .solve_soln(&mut soln)
            .unwrap()
            .into_state();
        if let Some(OdeSolverStopReason::RootFound(_t_root, root_idx)) = soln.stop_reason {
            problem.eqn_mut().set_model_index(root_idx);
            state.as_mut().apply_reset(problem.eqn()).unwrap();
        }
    }

    let time = soln.ts.clone();
    let q_c: Vec<_> = soln.ys.inner().row(0).into_iter().copied().collect();
    let q_p1: Vec<_> = soln.ys.inner().row(1).into_iter().copied().collect();

    let mut plot = Plot::new();
    plot.add_trace(
        Scatter::new(time.clone(), q_c)
            .mode(Mode::Lines)
            .name("q_c"),
    );
    plot.add_trace(Scatter::new(time, q_p1).mode(Mode::Lines).name("q_p1"));
    plot.set_layout(
        Layout::new()
            .x_axis(Axis::new().title("t [h]"))
            .y_axis(Axis::new().title("amount [ng]")),
    );

    let output_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"))
        .join("../../book/src/primer/images/drug-delivery-hybrid.html");
    let plot_html = plot.to_inline_html(Some("drug-delivery-hybrid"));
    fs::write(output_path, plot_html).expect("Unable to write file");
}