From 6c12728f6732b25918dec29216f32873e88349b4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 1 Apr 2024 14:23:35 -0400 Subject: [PATCH] init --- docs/src/api/catalyst_api.md | 4 + .../constraint_equations.md | 2 +- .../catalyst_functionality/dsl_description.md | 62 ++++ src/Catalyst.jl | 1 + src/expression_utils.jl | 7 + src/networkapi.jl | 49 ++++ src/reaction_network.jl | 74 ++++- src/reactionsystem.jl | 20 +- test/dsl/dsl_options.jl | 265 ++++++++++++++++++ 9 files changed, 474 insertions(+), 10 deletions(-) diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md index 5085a60d3f..2a29e3e918 100644 --- a/docs/src/api/catalyst_api.md +++ b/docs/src/api/catalyst_api.md @@ -157,6 +157,10 @@ species nonspecies reactionparams reactions +has_diff_equations +diff_equations +has_alg_equations +alg_equations numspecies numparams numreactions diff --git a/docs/src/catalyst_functionality/constraint_equations.md b/docs/src/catalyst_functionality/constraint_equations.md index b64cf6f8f1..25ef0b8d7b 100644 --- a/docs/src/catalyst_functionality/constraint_equations.md +++ b/docs/src/catalyst_functionality/constraint_equations.md @@ -17,7 +17,7 @@ $\lambda$. For now we'll assume the cell can grow indefinitely. We'll also keep track of one protein $P(t)$, which is produced at a rate proportional $V$ and can be degraded. -## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constratins) +## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constraints) There are several ways we can create our Catalyst model with the two reactions and ODE for $V(t)$. One approach is to use compositional modeling, create diff --git a/docs/src/catalyst_functionality/dsl_description.md b/docs/src/catalyst_functionality/dsl_description.md index 16f31c4588..e53700386c 100644 --- a/docs/src/catalyst_functionality/dsl_description.md +++ b/docs/src/catalyst_functionality/dsl_description.md @@ -701,3 +701,65 @@ rn = @reaction_network begin end nothing # hide ``` + +## Incorporating (differential) equations into reaction network models +Some models cannot be purely described as a reaction network. E.g. consider the growth of a cell, where the rate of change in cell's volume depends on some growth factor. Here, the cell's volume would be described by a normal equation. Such equations can be incorporated into a model using the `@equations` option. Here, we create a model where a growth factor ($G$) is produced and degraded at a linear rates, and the rate of change in cell volume ($V$) os linear to the amount of growth factor: +```@example eqs1 +using Catalyst #hide +rn = @reaction_network begin + @equations begin + D(V) ~ G + end + (p,d), 0 <--> G +end +``` +Here, `D(V)` indicates the (time) derivative with respect to `D`. The differential equation left and right hand sides are separated by a `~`. The left-hand side should contain differential only, the right hand side can contain any algebraic expression. + +We can check the differential equation corresponding to this reaction network using latexify: +```@example eqs1 +using Latexify +latexify(rn; form=:ode) +``` +We can also simulate it using the normal syntax +```@example eqs1 +using DifferentialEquations, Plots # hide +u0 = [:G => 0.0, :V => 0.1] +ps = [:p => 1.0, :d => 0.5] +oprob = ODEProblem(rn, u0, (0.0, 1.0), ps) +sol = solve(oprob) +plot(sol) +``` +Here, growth is indefinite. To improve the model, [a callback](@ref advanced_simulations_callbacks) can be used to half the volume (cell division) once some threshold is reached. + +When creating differential equations this way, the subject of the differential is automatically inferred to be a variable, however, any component on the right-hand side must be declare somewhere in the macro. E.g. to add a scaling parameter ($k$), we must declare that $k$ is a parmaeter using the `@paraemters` option: +```@example eqs1 +rn = @reaction_network begin + @parameters k + @equations begin + D(V) ~ k*G + end + (p,d), 0 <--> G +end +nothing #hide +``` + +It is possible to add several equations to the model. In this case, each have a separate line. E.g. to keep track of a supply of nutrition ($N$) in the growth media, we can use: +```@example eqs1 +rn = @reaction_network begin + @equations begin + D(V) ~ G + D(N) ~ -G + end + (p,d), 0 <--> G +end +nothing #hide +``` + +When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using: +```@example eqs1 +rn = @reaction_network begin + @equations D(V) ~ G + (p,d), 0 <--> G +end +nothing # hide +``` diff --git a/src/Catalyst.jl b/src/Catalyst.jl index cf50cfccc0..5fc467d4f2 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -86,6 +86,7 @@ export mm, mmr, hill, hillr, hillar # functions to query network properties include("networkapi.jl") export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap +export has_diff_equations, diff_equations, has_alg_equations, alg_equations export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap export dependants, dependents, substoichmat, prodstoichmat, netstoichmat diff --git a/src/expression_utils.jl b/src/expression_utils.jl index b13f78b512..ded4d31e2c 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -11,6 +11,13 @@ function get_tup_arg(ex::ExprValues, i::Int) return ex.args[i] end +# Some options takes input on form that is either `@option ...` or `@option begin ... end`. +# This transforms input of the latter form to the former (with only one line in the `begin ... end` block) +function option_block_form(expr) + (expr.head == :block) && return expr + return Expr(:block, expr) +end + # In variable/species/parameters on the forms like: # X # X = 1.0 diff --git a/src/networkapi.jl b/src/networkapi.jl index c1facbaaee..7cd1852db2 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -97,6 +97,55 @@ function reactions(network) [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])] end +""" + diff_equations(network) +Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are differential equations (contains a derivative with respect to any variable). +Notes: +- If `ModelingToolkit.get_systems(network)` is not empty, will allocate. +""" +function diff_equations(network) + eqs = equations(network) + filter!(!isreaction, eqs) + systems = filter_nonrxsys(network) + isempty(systems) && (return rxs) + [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])] +end + +""" + has_diff_equations(network) +Given a [`ReactionSystem`](@ref), check whether it contain any differential equations (i.e. in addition to those generated through reactions). +Notes: +- If `ModelingToolkit.get_systems(network)` is not empty, will allocate. +""" +function has_diff_equations(network) + return !isempty(diff_equations(network)) +end + +""" + alg_equations(network) +Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are algebraic equations (does not contain any derivatives). +Notes: +- If `ModelingToolkit.get_systems(network)` is not empty, will allocate. +""" +function alg_equations(network) + eqs = equations(network) + filter!(!isreaction, eqs) + filter!(!isreaction, eqs) + systems = filter_nonrxsys(network) + isempty(systems) && (return rxs) + [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])] +end + +""" + has_alg_equations(network) +Given a [`ReactionSystem`](@ref), check whether it contain any algebraic equations. +Notes: +- If `ModelingToolkit.get_systems(network)` is not empty, will allocate. +""" +function has_alg_equations(network) + return !isempty(alg_equations(network)) +end + """ speciesmap(network) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 3ca59d4c6a..f314ad44f7 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -82,7 +82,8 @@ const forbidden_variables_error = let end # Declares the keys used for various options. -const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling) +const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling, + :differentials, :equations, :continuous_events, :discrete_events) ### The @species macro, basically a copy of the @variables macro. ### macro species(ex...) @@ -353,7 +354,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) option_lines = Expr[x for x in ex.args if x.head == :macrocall] # Get macro options. - length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) && error("Some options where given multiple times.") + if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) + error("Some options where given multiple times.") + end options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg, option_lines)) @@ -361,12 +364,18 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) default_reaction_metadata = :([]) check_default_noise_scaling!(default_reaction_metadata, options) compound_expr, compound_species = read_compound_options(options) + continuous_events_expr = read_events_option(options, :continuous_events) + discrete_events_expr = read_events_option(options, :discrete_events) # Parses reactions, species, and parameters. reactions = get_reactions(reaction_lines) species_declared = [extract_syms(options, :species); compound_species] parameters_declared = extract_syms(options, :parameters) - variables = extract_syms(options, :variables) + variables_declared = extract_syms(options, :variables) + + # Reads more options. + vars_extracted, add_default_diff, equations = read_equations_options(options, variables_declared) + variables = vcat(variables_declared, vars_extracted) # handle independent variables if haskey(options, :ivs) @@ -391,6 +400,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) species = vcat(species_declared, species_extracted) parameters = vcat(parameters_declared, parameters_extracted) + # Create differential expression. + diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables]) + # Checks for input errors. (sum(length.([reaction_lines, option_lines])) != length(ex.args)) && error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.") @@ -402,7 +414,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) # Creates expressions corresponding to actual code from the internal DSL representation. sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs) - vexprs = haskey(options, :variables) ? options[:variables] : :() + vexprs = get_sexpr(vars_extracted, options, :variables) pexprs = get_pexpr(parameters_extracted, options) ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps") vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars") @@ -412,6 +424,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) for reaction in reactions push!(rxexprs.args, get_rxexprs(reaction)) end + for equation in equations + push!(rxexprs.args, equation) + end quote $ps @@ -420,11 +435,13 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) $sps $observed_vars $comps + $diffexpr Catalyst.remake_ReactionSystem_internal( Catalyst.make_ReactionSystem_internal( $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms), - $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs); + $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs, + continuous_events = $continuous_events_expr, discrete_events = $discrete_events_expr); default_reaction_metadata = $default_reaction_metadata ) end @@ -810,6 +827,53 @@ function check_default_noise_scaling!(default_reaction_metadata, options) end end +# Creates an expression declaring differentials. +function create_differential_expr(options, add_default_diff, used_syms) + # Creates the differential expression. + # If differentials was provided as options, this is used as the initial expression. + # If the default differential (D(...)) was used in equations, this is added to the expression. + diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] : MacroTools.striplines(:(begin end))) + diffexpr = option_block_form(diffexpr) + add_default_diff && push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM)))) + + # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere. + for dexpr in diffexpr.args + (dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.") + (dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.") + in(dexpr.args[1], used_syms) && error("Differential name ($(dexpr.args[1])) is also a species, variable, or parameter. This is ambigious and not allowed.") + in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.") + end + + return diffexpr +end + +# Read the events (continious or discrete) provided as options to the DSL. Returns an expression which evalutes to these. +function read_events_option(options, event_type::Symbol) + # Prepares the events, if required to, converts them to block form. + (event_type in [:continuous_events]) || error("Trying to read an unsupported event type.") + events_input = haskey(options, event_type) ? options[event_type].args[3] : :(begin end) + events_input = option_block_form(events_input) + + # Goes throgh the events, checks for errors, and adds them to the output vector. + events_expr = :([]) + for arg in events_input.args + # Formatting error checks. + # NOTE: Maybe we should move these deeper into the system (rather than the DSL), throwing errors more generally? + if (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3 + error("Events should be on form `condition => affect`, separated by a `=>`. This appears not to be the case for: $(arg).") + end + if (arg.args[2] != :call) && (event_type == :continuous_events) + error("The condition part of continious events (the left-hand side) must be a vector. This is not the case for: $(arg).") + end + if arg.args[3] != :call + error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).") + end + + push!(events_expr.args, arg) + end + return events_expr +end + ### Functionality for expanding function call to custom and specific functions ### #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression. diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 6256fa305d..7eab89973e 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -546,6 +546,11 @@ struct ReactionSystem{V <: NetworkProperties} <: name, systems, defaults, connection_type, nps, cls, cevs, devs, metadata, complete; checks::Bool = true) + # Filters away any potential obervables from `states` and `spcs`. + obs_vars = [obs_eq.lhs for obs_eq in observed] + states = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), states) + spcs = filter(spc -> !any(isequal(spc, obs_var) for obs_var in obs_vars), spcs) + # unit checks are for ODEs and Reactions only currently nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation] if checks && isempty(sivs) @@ -1409,7 +1414,7 @@ end # merge constraint components with the ReactionSystem components # also handles removing BC and constant species -function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false) +function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false, zero_derivatives = false) # if there are BC species, put them after the independent species rssts = get_unknowns(rs) sts = any(isbc, rssts) ? vcat(ists, filter(isbc, rssts)) : ists @@ -1501,7 +1506,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) ists, ispcs = get_indep_sts(fullrs, remove_conserved) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, include_zero_odes) - eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) + eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, zero_derivatives=true) ODESystem(eqs, get_iv(fullrs), sts, ps; observed = obs, @@ -1542,7 +1547,6 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name ists, ispcs = get_indep_sts(fullrs, remove_conserved) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, as_odes = false, include_zero_odes) - error_if_constraint_odes(NonlinearSystem, fullrs) eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) NonlinearSystem(eqs, sts, ps; @@ -1663,12 +1667,20 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, check_length = false, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, remove_conserved = false, - checks = false, kwargs...) + checks = false, structural_simplify=false, kwargs...) u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks, remove_conserved) osys = complete(osys) + + # Handles potential Differential algebraic equations. + if structural_simplify + (osys = MT.structural_simplify(osys)) + # elseif has_alg_equations(rs) + # error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.") + end + return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...) end diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index ad791705cd..5ea24980c7 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -634,4 +634,269 @@ let @observables $X ~ X1 + X2 (k1,k2), X1 <--> X2 end +end + + +### Algebraic Equations ### + +# Checks creation of basic network. +# Check indexing of output solution. +# Check that DAE is solved correctly. +let + rn = @reaction_network rn begin + @parameters k + @variables X(t) Y(t) + @equations begin + X + 5 ~ k*S + 3Y + X ~ S + X*d + end + (p,d), 0 <--> S + end + + @unpack X, Y, S, k, p, d = rn + + # Checks that the internal structures have the correct lengths. + @test length(species(rn)) == 1 + @test length(states(rn)) == 3 + @test length(reactions(rn)) == 2 + @test length(equations(rn)) == 4 + @test !has_diff_equations(rn) + @test isequal(diff_equations(rn), []) + @test has_alg_equations(rn) + @test isequal(alg_equations(rn), [X + 5 ~ k*S, 3Y + X ~ S + X*d]) + + # Checks that the internal structures contain the correct stuff, and are correctly sorted. + @test isspecies(states(rn)[1]) + @test !isspecies(states(rn)[2]) + @test !isspecies(states(rn)[3]) + @test equations(rn)[1] isa Reaction + @test equations(rn)[2] isa Reaction + @test equations(rn)[3] isa Equation + @test equations(rn)[3] isa Equation + @test isequal(equations(rn)[3], X + 5 ~ k*S) + @test isequal(equations(rn)[4], 3Y + X ~ S + X*d) + + # Checks that simulations has the correct output + u0 = Dict([S => 1 + rand(rng), X => 1 + rand(rng), Y => 1 + rand(rng)]) + ps = Dict([p => 1 + rand(rng), d => 1 + rand(rng), k => 1 + rand(rng)]) + oprob = ODEProblem(rn, u0, (0.0, 10000.0), ps; structural_simplify=true) + sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9) + @test sol[S][end] ≈ ps[p]/ps[d] + @test sol[X] .+ 5 ≈ sol[k] .*sol[S] + @test 3*sol[Y] .+ sol[X] ≈ sol[S] .+ sol[X].*sol[d] +end + +# Checks that block form is not required when only a single equation is used. +let + rn1 = @reaction_network rn begin + @parameters k + @variables X(t) + @equations X + 2 ~ k*S + (p,d), 0 <--> S + end + rn2 = @reaction_network rn begin + @parameters k + @variables X(t) + @equations begin + X + 2 ~ k*S + end + (p,d), 0 <--> S + end + @test rn1 == rn2 +end + +# Tries for reaction system without any reactions (just an equation). +# Tries with interpolating a value into an equation. +# Tries using rn.X notation for designating variables. +# Tries for empty parameter vector. +let + c = 6.0 + rn = complete(@reaction_network begin + @variables X(t) + @equations 2X ~ $c - X + end) + + u0 = [rn.X => 0.0] + ps = [] + oprob = ODEProblem(rn, u0, (0.0, 100.0), ps; structural_simplify=true) + sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9) + @test sol[rn.X][end] ≈ 2.0 +end + +# Checks hierarchical model. +let + base_rn = @reaction_network begin + @variables V1(t) + @equations begin + X*3V1 ~ X - 2 + end + (p,d), 0 <--> X + end + @unpack X, V1, p, d = base_rn + + internal_rn = @reaction_network begin + @variables V2(t) + @equations begin + X*4V2 ~ X - 3 + end + (p,d), 0 <--> X + end + + rn = compose(base_rn, [internal_rn]) + + u0 = [V1 => 1.0, X => 3.0, internal_rn.V2 => 2.0, internal_rn.X => 4.0] + ps = [p => 1.0, d => 0.2, internal_rn.p => 2.0, internal_rn.d => 0.5] + oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps; structural_simplify=true) + sol = solve(oprob, Rosenbrock23(); abstol=1e-9, reltol=1e-9) + + @test sol[X][end] ≈ 5.0 + @test sol[X][end]*3*sol[V1][end] ≈ sol[X][end] - 2 + @test sol[internal_rn.X][end] ≈ 4.0 +end + +# Check for combined differential and algebraic equation. +# Check indexing of output solution using Symbols. +let + rn = @reaction_network rn begin + @parameters k + @variables X(t) Y(t) + @equations begin + X + 5 ~ k*S + D(Y) ~ X + S - 5*Y + end + (p,d), 0 <--> S + end + + # Checks that the internal structures have the correct lengths. + @test length(species(rn)) == 1 + @test length(states(rn)) == 3 + @test length(reactions(rn)) == 2 + @test length(equations(rn)) == 4 + @test has_diff_equations(rn) + @test length(diff_equations(rn)) == 1 + @test has_alg_equations(rn) + @test length(alg_equations(rn)) == 1 + + # Checks that the internal structures contain the correct stuff, and are correctly sorted. + @test isspecies(states(rn)[1]) + @test !isspecies(states(rn)[2]) + @test !isspecies(states(rn)[3]) + @test equations(rn)[1] isa Reaction + @test equations(rn)[2] isa Reaction + @test equations(rn)[3] isa Equation + @test equations(rn)[3] isa Equation + + # Checks that simulations has the correct output + u0 = Dict([S => 1 + rand(rng), X => 1 + rand(rng), Y => 1 + rand(rng)]) + ps = Dict([p => 1 + rand(rng), d => 1 + rand(rng), k => 1 + rand(rng)]) + oprob = ODEProblem(rn, u0, (0.0, 10000.0), ps; structural_simplify=true) + sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9) + @test sol[:S][end] ≈ sol[:p]/sol[:d] + @test sol[:X] .+ 5 ≈ sol[:k] .*sol[:S] + @test 5*sol[:Y][end] ≈ sol[:S][end] + sol[:X][end] +end + +# Tests that various erroneous declarations throw errors. +let + # Using = instead of ~ (for equation). + @test_throws Exception @eval @reaction_network begin + @variables X(t) + @equations X = 1 - S + (p,d), 0 <--> S + end + + # Equation with component undeclared elsewhere. + @test_throws Exception @eval @reaction_network begin + @equations X ~ p - S + (P,D), 0 <--> S + end +end + +### Events ### + +# Compares models with complicated events that are created programmatically/with the DSL. +# Checks that simulations are correct. +# Checks that various simulation inputs works. +# Checks continuous, discrete, preset time, and periodic events. +# Tests event affecting non-species components. + +let + # Creates model via DSL. + rn_dsl = @reaction_network rn begin + @parameters thres=1.0 dY_up + @variables Z(t) + @continuous_events begin + [t - 2.5] => [p ~ p + 0.2] + [X - thres, Y - X] => [X ~ X - 0.5, Z ~ Z + 0.1] + end + @discrete_events begin + 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] + [1.0, 5.0] => [p ~ p - 0.1] + [Z > Y, Z > X] => [Z ~ Z - 0.1] + end + + (p, dX), 0 <--> X + (p, dY), 0 <--> Y + end + + # Creates model programmatically. + @variables t Z(t) + @species X(t) Y(t) + @parameters p dX dY thres=1.0 dY_up + rxs = [ + Reaction(p, nothing, [X], nothing, [1]) + Reaction(dX, [X], nothing, [1], nothing) + Reaction(p, nothing, [Y], nothing, [1]) + Reaction(dY, [Y], nothing, [1], nothing) + ] + continuous_events = [ + t - 2.5 => p ~ p + 0.2 + [X - thres, Y - X] => [X ~ X - 0.5, Z ~ Z + 0.1] + ] + discrete_events = [ + 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] + [1.0, 5.0] => [p ~ p - 0.1] + [Z > Y, Z > X] => [Z ~ Z - 0.1] + ] + rn_prog = ReactionSystem([rx1, rx2, eq], t; continuous_events, discrete_events, name=:rn) + + # Tests that approaches yield identical results. + @test isequal(rn_dsl, rn_prog) + + u0 = [X => 1.0, Y => 0.5, Z => 0.25] + tspan = (0.0, 20.0) + ps = [p => 1.0, dX => 0.5, dY => 0.5, dY_up => 0.1] + + sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5()) + sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5()) + @test sol_dsl == sol_prog +end + +# Compares DLS events to those given as callbacks. +# Checks that events works when given to SDEs. +let + # Creates models. + rn = @reaction_network begin + (p, d), 0 <--> X + end + rn_events = @reaction_network begin + @discrete_events begin + [5.0, 10.0] => [X ~ X + 100.0] + end + @continuous_events begin + [X - 90.0] => [X ~ X + 10.0] + end + (p, d), 0 <--> X + end + cb_disc = PresetTimeCallback([5.0, 10.0], int -> (int[:X] += 100.0)) + cb_cont = ContinuousCallback((u, t, int) -> (u[1] - 90.0), int -> (int[:X] += 10.0)) + + # Simulates models,. + u0 = [:X => 100.0] + tspan = (0.0, 50.0) + ps = [:p => 100.0, :d => 1.0] + sol = solve(SDEProblem(rn, u0, tspan, ps), ImplicitEM(); seed = 1234, callback = CallbackSet(cb_disc, cb_cont)) + sol_events = solve(SDEProblem(rn_events, u0, tspan, ps), ImplicitEM(); seed = 1234) + + @test sol == sol_events end \ No newline at end of file