0
votes

I am new to Julia (version 1.0.2) and currently trying the @reaction_network from the package DiffEqBiological (also current version, I can't find the version number here):

tspan = (0.0, 50.0);
y0 = [100.0 50.0 0.0 0.0] #[substrate enzyme complex product]
S, E = y0[1], y0[2]
for r in 0.1:0.1:1.0
    println("Creating Michaelis-Menten reaction model...")
    r1, r2, r3 = r, r, r;
    michaelismenten = @reaction_network rType begin
        r1, S + E ⟶ C
        r2, C ⟶ S + E
        r3, C ⟶ P + E
    end

    y = ODEProblem(michaelismenten, y0, tspan)
    sol = solve(y, CVODE_BDF(), reltol=1e-8, abstol=1e-8)
end

When I try to compile this, I get the syntax error: unsupported 'const' declaration on local variable around C:\Users\...\.julia\packages\DiffEqBiological\nujlA\src\reaction_network.jl:447. I tried really hard to find the error and maybe my knowledge is to low to actually understand similar questions, as I tried the answers, but still get the same error message. I tried the usual

r = 0.1
michaelismenten = @reaction_network rType begin
            r, S + E ⟶ C
            r, C ⟶ S + E
            r, C ⟶ P + E
        end
y = ODEProblem(michaelismenten, y0, tspan)
sol = solve(y, CVODE_BDF(), reltol=1e-8, abstol=1e-8)

and this works perfectly fine. As soon as I add the for-loop to vary my reaction rate, I get this error, meaning this code

for i = 1:10
    r = i/10 # here I thought that Julia maybe did not like the 0.1:0.1:1.0 from above
    michaelismenten = @reaction_network rType begin
                    r, S + E ⟶ C
                    r, C ⟶ S + E
                    r, C ⟶ P + E
                end
end

already throws the same error. That's why, I think that I am doing something wrong with the for loop, that I wanted to use to vary the reaction rate r so that I am able to animate my plots via

anim = @animate for r = 0.1:0.1:1.0
     # create michaelismenten and solve in this for-loop before plotting (see above)
     plot(sol)
end
gif(anim, "$S_$E_RRE.gif", fps = 15)

I hope, that someone can please help me find the solution for my problem. And please help me to make my asked question more clear with formatting.

To summarize the comments a bit:

  • declaring r as global does not work and throws the same error
  • wrapping the michaelismenten reaction network in a function with r as input throws the same error of unsupported 'const' declaration on local variable
1
I don't have the immediate answer to your problem, but it would likely help me, and others, to provide more details. I.e., where are y0 and tspan defined, as well as more details on what versions of Julia and the DiffEqBiological package you are running. Does the first example on github.com/JuliaDiffEq/DiffEqBiological.jl work?jvz
I added your suggested improvements (definitions of y0, tspan, etc) and tested the first example and this works fine.Lissa
I guess it is due to the scoping behaviour. What happens if in your third code block, you put global r on the line after starting the for-loop? (sorry, don't have access to Julia here to try myself). See also the section "Local scope" in the manual: docs.julialang.org/en/v1/manual/variables-and-scoping/…jvz
It still throws the same error message although I declared the variable as global. I looked into line 447 in reaction_network.jl where the error occurs & it seems to be this function (I don't know how broad I have to look when the error says "around reaction_network.jl:447" #Turns an array of expressions to a expression block with corresponding expressions. function expr_arr_to_block(exprs) block = :(begin end) foreach(expr -> push!(block.args, expr), exprs) return block end I am sorry for the weird code formatting, as I tried my best for it to look nice, but it didn't work.Lissa
What happens if you wrap the michaelismenten = @reaction_network rType begin r, S + E ⟶ C r, C ⟶ S + E r, C ⟶ P + E end in a function that takes in r and returns michaelismenten? If that works, you could just call that function.jvz

1 Answers

0
votes

I struggled quite a bit, but this code now works (for those who are interested in a solution ;) ). Thanks, zundertj, for your help! You triggered some ideas in my brain, which lead to this solution!

using Plots
using DifferentialEquations
using Distributions
using DiffEqBiological

function solveRRE(michaelismenten, y0, tspan, p)
    # create and solve problem
    println("Creating ODEProblem...")
    y = ODEProblem(michaelismenten, y0, tspan, p)
    println("Solving Reaction Rate Equation...")
    sol = solve(y, CVODE_BDF(), reltol=1e-8, abstol=1e-8)
    return sol
end

tspan = (0.0, 50.0);
y0 = [100.0 50.0 0.0 0.0] #[substrate enzyme complex product]
S, E = y0[1], y0[2]

michaelismenten = @reaction_network rType begin
    r, S + E ⟶ C
    r, C ⟶ S + E
    r, C ⟶ P + E
end r
anim = @animate for i in 0.1:0.1:1.0
    p = i
    sol = solveRRE(michaelismenten, y0, tspan, p)
    plot(sol)
end
gif(anim, "RRE.gif", fps = 15)

Another possible occuring error is that ffmpeg is not installed on Windows/Linux. For further instructions for that, this website helped me lots: https://video.stackexchange.com/questions/20495/how-do-i-set-up-and-use-ffmpeg-in-windows

Can someone please tell me, how to close this topic as solved? Thanks in advance!