# How to model a chemical reaction network with changing parameters? (Catalyst.jl)

**URL:** <https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489>\
**Category:** Chemistry\
**Tags:** question, catalyst\
**Created:** [April 14, 2022, 1:50pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489 "2022-04-14T13:50:06Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![govissers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/govissers/32/35435_2.png) [@govissers](https://discourse.julialang.org/u/govissers)\
**Post date:** [April 14, 2022, 1:50pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/1 "2022-04-14T13:50:06Z")

</div>

Hello,

I am hoping to model a chemical reaction network within a cell, where an organism may change the rate of a reaction given an external stimulus at time _t_. Is it possible to model a CRN wherea reaction constant, specified in **rates** changes at time _t_? Code pasted below.

```julia
using Catalyst
using Plots
using DifferentialEquations

rn = @reaction_network begin
    a, K --> RP + K
    b, RP --> P + RP
    (c, d), P + K <--> KP
    e, P --> 0
    f, RP --> 0
end a b c d e f

rates = (
    :a => 0.1,
    :b => 0.001,
    :c => 0.01,
    :d => 0.1,
    :e => 0.1,
    :f => 0.1
)
tspan = (0., 1000.)
u0 = [
    :K => 10,
    :RP => 0,
    :P => 0,
    :KP => 0
]

oprob = ODEProblem(rn, u0, tspan, rates)
solution = solve(oprob, Tsit5(), saveat=10.)

plot(solution)

```

Concretely, is it possible to generate a piecewise-like function in which the rate constant _e_ changes values at a given time _t_?

Thank you!  
Graeme

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [April 14, 2022, 2:33pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/2 "2022-04-14T14:33:34Z")

</div>

If there is a single change of rate at a know moment you could do something like this, where you change the value of a at t=1000 to 1.0  
But probably there are more elegant solutions.

```julia
rates = (
    :a => 0.1,
    :b => 0.001,
    :c => 0.01,
    :d => 0.1,
    :e => 0.1,
    :f => 0.1
)
tspan = (0., 1000.)
u0 = [
    :K => 10,
    :RP => 0,
    :P => 0,
    :KP => 0
]

oprob = ODEProblem(rn, u0, tspan, rates)
solution = solve(oprob, Tsit5(), saveat=10.)

rates = (
    :a => 1,
    :b => 0.001,
    :c => 0.01,
    :d => 0.1,
    :e => 0.1,
    :f => 0.1
)
tspan = (1000., 2000.)
u0 = [
    :K => solution.u[end][1],
    :RP => solution.u[end][2],
    :P => solution.u[end][3],
    :KP => solution.u[end][4]
]

oprob = ODEProblem(rn, u0, tspan, rates)
solutionpart2 = solve(oprob, Tsit5(), saveat=10.)

```

and then combine the two solution parts

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [April 14, 2022, 3:05pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/3 "2022-04-14T15:05:52Z")

</div>

There are a bunch of ways to achieve this. Catalyst doesn’t support events yet, but you could explicitly convert `rn` to a ModelingToolkit ODESystem and pass in a symbolic event. You could also just add an event via the DifferentialEquations.jl interface, see [here](https://diffeq.sciml.ai/latest/features/callback_functions/).

Finally, you could use a simple function for the rate expression to manually flip the rate at a specified time (note I then also set `tstops` in the call to `solve` to ensure the ODE solver steps to the switching time exactly):

```julia
using Catalyst, OrdinaryDiffEq, Plots
rate(k1,k2,tswitch,t) = (t < tswitch) ? k1 : k2

# this ensures our custom function works with ModelingToolkit/Symbolics.jl
@register rate(k1,k2,tswitch,t)   

rn = @reaction_network begin
       rate(k1,k2,tswitch,t), A --> 0
       end k1 k2 tswitch
p = (:k1 => 0.0, :k2 => 1.0, :tswitch => 1.0)
u0 = [:A => 10.0]
tspan = (0.0,3.0)
oprob = ODEProblem(rn,u0,tspan,p)
sol = solve(oprob, Tsit5(), tstops=[1.0])
plot(sol)

```

giving  
 ![plot](https://global.discourse-cdn.com/julialang/original/3X/6/f/6fbeb81348bd3dcd3c981219a650aa51ef46c6d8.png)

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [April 14, 2022, 3:50pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/4 "2022-04-14T15:50:57Z")

</div>

And here is a version that uses `IfElse.ifelse`, and so probably works better in applications that might need AD:

```julia
using Catalyst, OrdinaryDiffEq, Plots, IfElse

rate(k1,k2,tswitch,t) = IfElse.ifelse(t < tswitch, k1, k2)

rn = @reaction_network begin
       rate(k1,k2,tswitch,t), A --> 0
       end k1 k2 tswitch
p = (:k1 => 0.0, :k2 => 1.0, :tswitch => 1.0)
u0 = [:A => 10.0]
tspan = (0.0,3.0)
oprob = ODEProblem(rn,u0,tspan,p)
sol = solve(oprob, Tsit5(), tstops=[1.0])
plot(sol)

```

---

<div class="post-metadata">

**Author:** ![govissers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/govissers/32/35435_2.png) [@govissers](https://discourse.julialang.org/u/govissers)\
**Post date:** [April 14, 2022, 4:19pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/5 "2022-04-14T16:19:17Z")

</div>

Fantastic. I ended up using IfElse, and it did what I needed. Thank you!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [April 14, 2022, 5:38pm UTC](https://discourse.julialang.org/t/how-to-model-a-chemical-reaction-network-with-changing-parameters-catalyst-jl/79489/6 "2022-04-14T17:38:12Z")

</div>

👍