I have a question about the arrow types in Catalyst models when simulating jump processes. I build the model as a reaction network in Catalyst and simulate it as a JumpProblem. It’s a large model and so I’m trying to optimise the run times. I’ve noticed that using different arrow types can affect the run times (=> vs -->). I understand the difference between the two in that → automatically treats it as a mass action reaction and => just uses the user defined rate. When using the same seed, I expect my solutions to be equivalent no matter what arrows I’m using, as long as the propensities match.
See the MWE below (using a toy model), model 0 and model 3 are equivalent, but models 1 and 2 don’t match them, even though technically the reactions have the same propensities. Can anyone share any information on why this is? Can I treat my models as equivalent if their propensities are the same even though the same seed doesn’t produce the exact same outcome? Thanks in advance!
using Catalyst, JumpProcesses, Random, DataFrames, Statistics
# ============================================================================
# Model Definitions
# ============================================================================
# Model 0: All reactions use --> (mass action)
function model_arrows_0()
rn = @reaction_network begin
k_prod, ∅ --> A
k_bind, A + B --> AB
k_unbind, AB --> A + B
k_deg, A --> ∅
k_deg, B --> ∅
k_deg, AB --> ∅
end
return rn
end
# Model 1: Mixed arrows - unbinding uses explicit =>
function model_arrows_1()
rn = @reaction_network begin
k_prod, ∅ --> A
k_bind, A + B --> AB
k_unbind * AB, AB => A + B
k_deg, A --> ∅
k_deg, B --> ∅
k_deg, AB --> ∅
end
return rn
end
# Model 2: More explicit rates with =>
function model_arrows_2()
rn = @reaction_network begin
k_prod, ∅ --> A
k_bind * A * B, A + B => AB
k_unbind * AB, AB => A + B
k_deg, A --> ∅
k_deg, B --> ∅
k_deg, AB --> ∅
end
return rn
end
# Model 3: All reactions use explicit =>
function model_arrows_3()
rn = @reaction_network begin
k_prod, ∅ => A
k_bind * A * B, A + B => AB
k_unbind * AB, AB => A + B
k_deg * A, A => ∅
k_deg * B, B => ∅
k_deg * AB, AB => ∅
end
return rn
end
# ============================================================================
# Test Function
# ============================================================================
function run_stochastic_sim(model_func; seed=123, tspan=(0.0, 100.0), dt=0.5)
rn = model_func()
# Parameters
pars = Dict(:k_prod => 0.5, :k_bind => 0.01, :k_unbind => 0.1, :k_deg => 0.02)
# Initial conditions
u0 = Dict(:A => 10, :B => 5, :AB => 0)
# Create and solve jump problem
jp = JumpProblem(rn, u0, tspan, pars, aggregator=Direct(), save_positions=(false, false))
Random.seed!(seed)
sol = solve(jp, SSAStepper(), saveat=dt)
return sol
end
# ============================================================================
# Run All Models
# ============================================================================
sol0 = run_stochastic_sim(model_arrows_0)
sol1 = run_stochastic_sim(model_arrows_1)
sol2 = run_stochastic_sim(model_arrows_2)
sol3 = run_stochastic_sim(model_arrows_3)
println("\n" * "="^70)
println("All Models vs Model 0:")
println("="^70)
for i in 1:3
s = [sol0, sol1, sol2, sol3][i+1]
equal = isequal(sol0, s)
status = equal ? "✓ EQUAL" : "✗ DIFFERENT"
println("Model 0 vs Model $i: $status")
end