I'm new to using Julia and I'm attempting to simulate a reaction network while a disturbance changes the concentration of one species over time. I've learned that I can use a callback to alter the value of a specific species at a given time to simulate this external disturbance. Although I've read the documentation, but I couldn't find how I can do it. Here's the code in which I've tried to implement this simulation. At first, I defined the reaction network and its parameters:
using Catalyst
using DifferentialEquations
rn = @reaction_network begin
α, A + B --> 2B
β, B --> A
end α,β
p = [:α => 1, :β => 2]
tspan = (0.0,20.0)
u0 = [:A => 5.0, :B=> 5.0]
op = ODEProblem(rn, u0, tspan, p)
Next, for the callback that triggers at t=10
condition(u, t, integrator) = t == 10
to increase the value of species A by two
affect!(integrator) = integrator.??? += 2
But this is where I'm unsure about which parameters to use to change the value of A.
And the rest of my code:
cb = DiscreteCallback(condition, affect!)
sol = solve(prob, Tsit5(), callback = cb, tstops = [10.0])
using Plots
plot(sol)
I'd appreciate any help to improve my understanding in this area.
In this case, you want to say affect!(integrator) = integrator.u[1] += 2. You need to access an indexed state of the integrator. To figure out which index you need to use, you can check states(rn) to see the order of the unknowns that you need to use to make sure you're accessing the right one. If you had written the first reaction as B + A instead, then B would be the first state listed by states(rn) and you would need to increment integrator.u[2].
If you know the time that you want to create the disturbance exactly, I'd also recommend using a PresetTimeCallback instead (see here). DiscreteCallback needs to evaluate the condition after every integrator step, which can add unnecessary overhead.