Catalyst.jl - model reaction network with changing parameters during drug dosing

Specific Domains Modelling & Simulations
1

Hi, I am trying to implement a infusion dosing with my reaction network. The infusion rate will be changing depending on the time. If during the infusion duration, there’s a value for k_inf, otherwise it will be 0. and the dosing is QW. The code need to cover both with treatment and no treatment (where k_inf will be 0 all the time).

Here’s the code I have, and I ran into error on corporating ifelse with the reaction network. Could you provide some suggestions or best strategy moving forward please? Should I use callback instead, if so, how does callback work with reaction network? Thanks so much!!

I tried to implement ifelse as either a function or directly inside reaction_network, both give me errors that I cannot interpret and try to fix


get_ert_rate(total_dosing_time, dose_interval, infusion_duration, inf_rate) = IfElse.ifelse(
    t < total_dosing_time,
    # If within the treatment period, check the periodic window
    IfElse.ifelse(
        t % dose_interval < infusion_duration,
        inf_rate,
        0.0
    ),
    # If outside the treatment period, the rate is permanently zero
    0.0
)   


function model_with_drug()
    
    rn = @reaction_network begin

       @species begin
          plasma_drug(t)
          orange(t)
       end

       @parameters begin
           WT
           k_prod
           
           # treatment related
           dose
           infusion_duration
           dose_interval 
           total_dosing_time
           V
           C
       end

     @equations begin
         inf_rate ~ dose / infusion_duration
         k_inf ~ get_drug_rate(t, total_dosing_time, dose_interval, infusion_duration, inf_rate)

         r_admin ~ k_inf / V
         r_clear ~ C * plasma_drug
 
         r_orange ~ k_prod * orange * plasma_drug
     end

     r_admin, 0 => plasma_drug
     r_clear, plasma_drug => 0
    
     r_orange, 0 => orange
 end 
 return rn
end


#### for running simulations
p_dict = Dict{Symbol, Float64}()
p_dict[:dose] = 5.0
p_dict[:infusion_duration] = 1.5 # 1.5hrs infusion time, k_inf has value only during the 1.5 hrs
p_dict[:dose_interval] = 7 * 24 # QW dosing
p_dict[:total_dosing_time] = 10 * 7 * 24 # 10 cycles of QW dosing


u0 = zeros(2,1)
rn = model_with_drug()
tspan = (0.0, 26*7*24) # 26 weeks
sys = ModelingToolkit.structural_simplify(convert(ODESystem, rn))
prob = ODEProblem(sys, u0, tspan, p_dict)
sol = solve(prob, radau())
2

You probably want to model this with one or two discrete events using @discrete_events to turn it on and off. You can use this in the DSL, but it seems we haven’t merged any documentation on using it yet. You can see how it works in MTK9, which is similar notation, here.

3

Actually we do have a bit on them at: