Solution of SteadyStateProblem is extremely slow

[SmearingMap]

I switched some code over from DiffEqBiological to Catalyst and noticed that Catalyst was much, much slower.
The problem here is to determine steady states for a MAPK cascade for values of the input E1 that range over a few orders of magnitude (this is the same problem and with similar solved in this 1996 paper, so we're not talking about anything that should cause major problems at the solver level). Using DiffEqBiological I was able to solve this problem for 10,000 points spanning 5 orders of magnitude in under a minute (I never bothered timing, because it was fast enough not to be an issue). When I switched over to Catalyst, trying to do the same thing was not feasible at all.

This code shows the issue. When I ran it, I noticed a couple of things:

1. Making multiple problems is very slow. Using btime tells me that it takes about 200 milliseconds to make each problem, and that doubling the number of problems doubles the amount of time taken.
2. It only takes 35 milliseconds to solve 50 problems, but a lot more than that to solve 100. I started running this several minutes ago, and btime still hasn't told me how long it takes to solve 100 problems.

Again, this problem was trivial before switching over from DiffEqBiological to Catalyst, so it seems like it has to be an issue with Catalyst, since I was using the same version of DifferentialEquations along with DiffEqBiological.

using Catalyst
using DifferentialEquations
using BenchmarkTools

MAPK_cascade = @reaction_network begin
    
    # Raf phosphorylation
    a_1, Raf0 + E1 --> Raf0E1
    d_1, Raf0E1 --> Raf0 + E1
    k_1, Raf0E1 --> Raf1 + E1
    
    # Raf de-phosphorylation
    a_2, Raf1 + E2 --> Raf1E2
    d_2, Raf1E2 --> Raf1 + E2
    k_2, Raf1E2 --> Raf0 + E2
    
    # First MEK phosphorylation
    a_3, Raf1 + Mek0 --> Raf1Mek0
    d_3, Raf1Mek0 --> Raf1 + Mek0
    k_3, Raf1Mek0 --> Raf1 + Mek1
    
    # First MEK de-phosphorylation
    a_4, Mek1 + MekF --> Mek1MekF
    d_4, Mek1MekF --> Mek1 + MekF
    k_4, Mek1MekF --> Mek0 + MekF
    
    # Second MEK phosphorylation
    a_5, Raf1 + Mek1 --> Raf1Mek1
    d_5, Raf1Mek1 --> Raf1 + Mek1
    k_5, Raf1Mek0 --> Raf1 + Mek2
    
    # Second MEK de-phosphorylation
    a_6, Mek2 + MekF --> Mek2MekF
    d_6, Mek2MekF --> Mek2 + MekF
    k_6, Mek2MekF --> Mek1 + MekF
    
    # First ERK phosphorylation
    a_7, Mek2 + Erk0 --> Mek2Erk0
    d_7, Mek2Erk0 --> Mek2 + Erk0
    k_7, Mek2Erk0 --> Mek2 + Erk1
    
    # First ERK de-phosphorylation
    a_8, Erk1 + ErkF --> Erk1ErkF
    d_8, Erk1ErkF --> Erk1 + ErkF
    k_8, Erk1ErkF --> Erk0 + ErkF
    
    # Second ERK phosphorylation
    a_9, Mek2 + Erk1 --> Mek2Erk1
    d_9, Mek2Erk1 --> Mek2 + Erk1
    k_9, Mek2Erk1 --> Mek2 + Erk2
    
    # Second ERK de-phosphorylation
    a_10, Erk2 + ErkF --> Erk2ErkF
    d_10, Erk2ErkF --> Erk2 + ErkF
    k_10, Erk2ErkF --> Erk1 + ErkF
end a_1 d_1 k_1 a_2 d_2 k_2 a_3 d_3 k_3 a_4 d_4 k_4 a_5 d_5 k_5 a_6 d_6 k_6 a_7 d_7 k_7 a_8 d_8 k_8 a_9 d_9 k_9 a_10 d_10 k_10


function unbound_initial_condition(E1, E2,
                                   Raf,
                                   Mek, MekF,
                                   Erk, ErkF
                                    )
    #=
    Set all complex concentrations to zero initially.
    Setting entries is hardcoded to sidestep dictionary issue.
    =#
    x0 = zeros(length(MAPK_cascade.states))
    x0[1] = Raf
    x0[2] = E1
    x0[5] = E2
    x0[7] = Mek
    x0[10] = MekF
    x0[15] = Erk
    x0[18] = ErkF
    x0
end

function unbound_initial_condition(E1::Number)
    #=
    Set only E1 and fill in the rest with default values
    =#
    unbound_initial_condition(E1, 
        DEFAULT_CONCENTRATIONS[2:length(DEFAULT_CONCENTRATIONS)]...
    )
end

#= Setting default rates to obtain K_m = 100 nanomolar
   and rate constants to be order-of-magnitude consistent
   with Table 1 of the review by Piala 2014.
   Note: Ferrell 1996 uses K_m = 300 nanomolar for all
   Michaelis constants.
=#

DEFAULT_RATES = repeat([100.0, 1.0, 9.0], 10)


DEFAULT_CONCENTRATIONS = [
                              0.001, #E1
                              0.001, #E2
                              0.003, #Raf
                              1.2,   #Mek
                              0.001, #MekF
                              1.2,   #Erk
                              0.1    #ErkF
                            ];

xs = collect(10.0^i for i in -7:0.1:-2)

x0s = unbound_initial_condition.(xs)

#=
Define functions needed for the steady state response function
=#
odesys = convert(ODESystem, MAPK_cascade)
make_problem(x0) = SteadyStateProblem(odesys, x0, DEFAULT_RATES)
solve_problem(p) = solve(p, DynamicSS(Rodas5()); abstol=1e-4, reltol=1e-3)
response_function(x) = solve_problem(make_problem(unbound_initial_condition(x)))


println("Time to make 50 problems:")
@btime make_problem.(x0s)
problems = make_problem.(x0s)
println("Time to solve 50 problems:")
@btime solve_problem.(problems)


#= Increasing the number of problems seems to increase the time faster than linearly. This could be a stiffness problem, but I was able to solve this problem on 10,000 points in just a few seconds with DiffEqBiological.
=#
 
xs = collect(10.0^i for i in -7:0.05:2)
x0s = unbound_initial_condition.(xs);

println("Time to make 100 problems:")
@btime make_problem.(x0s)
problems = make_problem.(x0s)
println("Time to solve 100 problems:")
@btime solve_problem.(problems)

[ChrisRackauckas]

@shashi could you look at the timings here?

[SmearingMap]

Now that the script is finished, btime tells me that the solution of 100 problems only took 224 milliseconds, but the script ran for about 20 minutes to get that results. I don't really know how btime works, but I do know that trying to solve a number of problems that DiffEqBiological handled in about a minute had my laptop at work for about 30 minutes before I stopped it.

Let me know if there's anything else you need.

[ChrisRackauckas]

Sounds like something had an absurd compile time.

[isaacsas]

@mtrenfield I can't seem to run this code, getting an error at

x0s = unbound_initial_condition.(xs)

that

ERROR: LoadError: KeyError: key "Raf0" not found

Also, what version of the various packages are you on?

[SmearingMap]

@isaacsas I'm using Julia 1.4.1, Catalyst v5.1.1, and DifferentialEquations v6.16.0. I did just check and I can run this code copy-pasted from a text file.

I'll edit this example and hard-code the initial condition so it's easier to work with.

[SmearingMap]

@isaacsas I edited the code in the opening issue to remove dependence on the dictionary, which is what was causing problems. Let me know if you're able to run the code above as it is now.

[isaacsas]

Those are quite old versions of Julia and Catalyst at this point. In particular, being on that old of a Catalyst version means you won't have the latest versions of ModelingToolkit and Symbolics, which are what handle generating the compiled ODE derivative function (and should have had a lot of performance improvements since Catalyst v5, which is when we first switched to using them so not well optimized). I'd also strong recommend updating to Julia 1.6.2 as there are many improvements that go along with it in compile times and time to first plot.

Can you update to the latest Catalyst (v9) and DifferentialEquations.jl (v6.19). If you don't want to alter your default environment you can create a blank temporary environment by hitting "]" in the Julia REPL to enter package mode, and then entering

1. activate --temp
2. add Catalyst, DifferentialEquations, BenchmarkTools
3. type "st" to check you have the versions I mentioned above installed.

and then hitting "delete" to leave package mode. Your REPL will then be in a new environment with just the three packages for your script available.

At that point please can you see if your code still has the issues of your original post?

Thanks!

[isaacsas]

I just ran your script on Catalyst v9. The code

println("Time to make 50 problems:")
@btime make_problem.($x0s)
problems = make_problem.(x0s)
println("Time to solve 50 problems:")
@btime solve_problem.($problems)

gives

Time to make 50 problems:
  164.916 ms (1947548 allocations: 78.54 MiB)
Time to solve 50 problems:
  14.650 ms (13167 allocations: 2.22 MiB)

Note, @btime will run the command you give it a bunch of times to collect statistics on how long the command takes. It took maybe a minute to run these two commands on my laptop.

[isaacsas]

Also, note that the latest Catalyst requires Julia 1.5 at a minimum, but I'd suggest just using the latest Julia.

[SmearingMap]

The problem does not seem to persist after updating packages, so I'm marking this as closed.

[shashi]

I'm so glad.