From 012ca2c4eed9e47622cae853a09edfee5b01a66c Mon Sep 17 00:00:00 2001 From: paulflang Date: Mon, 2 Dec 2024 12:57:11 +0100 Subject: [PATCH 1/9] add test for is_event_assignment --- test/simresults.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/simresults.jl b/test/simresults.jl index 16a0f2c..aad1a7d 100644 --- a/test/simresults.jl +++ b/test/simresults.jl @@ -8,7 +8,8 @@ const case_ids = [7, # boundary_condition 140, # compartment size overridden with assignmentRule 170, # Model using parameters and rules only 325, # One reactions and two rate rules with four species in a 2D compartment - 679 # Initial value calculated by assignmentRule in compartment of non-unit size # 1208, # Non-hOSU species with rateRule in variable compartment -> require MTK fix. + 679, # Initial value calculated by assignmentRule in compartment of non-unit size # 1208, # Non-hOSU species with rateRule in variable compartment -> require MTK fix. + 944 # Non-constant parameter as target of event_assignment ] const logdir = joinpath(@__DIR__, "logs") From f7e0f82ae8593fe249257be9fc3b7eefdcdfa056 Mon Sep 17 00:00:00 2001 From: paulflang Date: Mon, 2 Dec 2024 14:18:39 +0100 Subject: [PATCH 2/9] fix is_event_assignment --- src/drafts.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/drafts.jl b/src/drafts.jl index 2d313b2..4a0c3c2 100644 --- a/src/drafts.jl +++ b/src/drafts.jl @@ -1,5 +1,5 @@ function is_event_assignment(k, model) - for ev in values(model.events) + for ev in last.(model.events) for as in ev.event_assignments if as.variable == k return true From 27c1def456d533114ddd6319566be0b3cc2e521e Mon Sep 17 00:00:00 2001 From: paulflang Date: Mon, 2 Dec 2024 14:22:06 +0100 Subject: [PATCH 3/9] bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3fb013c..3a0203c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SBMLToolkit" uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e" authors = ["paulflang", "anandijain"] -version = "0.1.29" +version = "0.1.30" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" From 087083ff266fe2f161d4cfed931c01f0fa7201a4 Mon Sep 17 00:00:00 2001 From: paulflang Date: Mon, 2 Dec 2024 16:51:25 +0100 Subject: [PATCH 4/9] define zero ODE for parameters/compartments with event assignments --- src/systems.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/systems.jl b/src/systems.jl index ffc6e2f..5dc53d9 100644 --- a/src/systems.jl +++ b/src/systems.jl @@ -87,11 +87,20 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...) # Todo: requires # Dict(Catalyst.DEFAULT_IV.val => 0))) push!(raterules_subs, rhs ~ o.lhs) end + zero_rates = [] + for (k, v) in merge(model.parameters, model.compartments) + if is_event_assignment(k, model) + if v.constant + ErrorException("$k appears in event assignment but should be constant.") + end + push!(zero_rates, D(create_var(k, IV; isbcspecies = true)) ~ 0) + end + end if haskey(kwargs, :defaults) defs = ModelingToolkit._merge(defs, kwargs[:defaults]) kwargs = filter(x -> !isequal(first(x), :defaults), kwargs) end - rs = ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged...], + rs = ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged..., zero_rates...], IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML), continuous_events = get_events(model), From 20a54303e9e3848acdfbde5f8051093cf26f7cb5 Mon Sep 17 00:00:00 2001 From: paulflang Date: Mon, 2 Dec 2024 17:52:17 +0100 Subject: [PATCH 5/9] format --- src/systems.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems.jl b/src/systems.jl index 5dc53d9..e20a7ef 100644 --- a/src/systems.jl +++ b/src/systems.jl @@ -100,7 +100,8 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...) # Todo: requires defs = ModelingToolkit._merge(defs, kwargs[:defaults]) kwargs = filter(x -> !isequal(first(x), :defaults), kwargs) end - rs = ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged..., zero_rates...], + rs = ReactionSystem( + [rxs..., algrules..., raterules_subs..., obsrules_rearranged..., zero_rates...], IV, first.(u0map), first.(parammap); defaults = defs, name = gensym(:SBML), continuous_events = get_events(model), From 79686c5866135a9539660379979eae02a2b93515 Mon Sep 17 00:00:00 2001 From: paulflang Date: Tue, 10 Dec 2024 21:34:36 +0100 Subject: [PATCH 6/9] handle ifelse in get_massaction; handle occurences of reaction IDs in math elements --- src/reactions.jl | 1 + src/rules.jl | 46 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/src/reactions.jl b/src/reactions.jl index 8bf819b..93b581d 100644 --- a/src/reactions.jl +++ b/src/reactions.jl @@ -177,6 +177,7 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing}, check_args(::Val{false}, x::SymbolicUtils.BasicSymbolic{Real}) = isspecies(x) ? NaN : 0 # Species vs Parameter leaf node check_args(::Real) = 0 # Real leaf node + check_args(::SymbolicUtils.BasicSymbolic{Bool}) = NaN check_args(x) = throw(ErrorException("Cannot handle $(typeof(x)) types.")) # Unknow leaf node if isnothing(reactants) && isnothing(stoich) diff --git a/src/rules.jl b/src/rules.jl index d6f99a7..6f609fe 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -21,12 +21,56 @@ function get_rules(model) algeqs, obseqs, raterules end + +extensive_kinetic_math(m::SBML.Model, formula::SBML.Math) = SBML.interpret_math( # TODO: move this to SBML.jl + formula, + map_apply = (x, rec) -> SBML.MathApply(x.fn, rec.(x.args)), + map_const = identity, + map_ident = (x::SBML.MathIdent) -> begin + haskey(m.reactions, x.id) && return m.reactions[x.id].kinetic_math + haskey(m.species, x.id) || return x + sp = m.species[x.id] + sp.only_substance_units && return x + if isnothing(m.compartments[sp.compartment].size) && + !seemsdefined(sp.compartment, m) + if m.compartments[sp.compartment].spatial_dimensions == 0 + # If the comparment ID doesn't seem directly defined anywhere + # and it is a zero-dimensional unsized compartment, just avoid + # any sizing questions. + return x + else + # In case the compartment is expected to be defined, complain. + throw( + DomainError( + sp.compartment, + "compartment size is insufficiently defined", + ), + ) + end + else + # Now we are sure that the model either has the compartment with + # constant size, or the definition is easily reachable. So just use + # the compartment ID as a variable to compute the concentration (or + # area-centration etc, with different dimensionalities) by dividing + # it. + return SBML.MathApply("/", [x, SBML.MathIdent(sp.compartment)]) + end + end, + map_lambda = (x, _) -> error( + ErrorException("converting lambdas to extensive kinetic math is not supported"), + ), + map_time = identity, + map_avogadro = identity, + map_value = identity, +) + + function get_var_and_assignment(model, rule) if !haskey(merge(model.species, model.compartments, model.parameters), rule.variable) error("Cannot find target for rule with ID `$(rule.variable)`") end var = create_var(rule.variable, IV) - math = SBML.extensive_kinetic_math(model, rule.math) + math = extensive_kinetic_math(model, rule.math) vc = get_volume_correction(model, rule.variable) if !isnothing(vc) math = SBML.MathApply("*", [SBML.MathIdent(vc), math]) From 858a2b5e4a6470616edcd66fe6d7e6749e9add11 Mon Sep 17 00:00:00 2001 From: paulflang Date: Tue, 10 Dec 2024 21:35:20 +0100 Subject: [PATCH 7/9] test handling of ifelse in get_massaction --- test/reactions.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/reactions.jl b/test/reactions.jl index 0d30fbc..2223b05 100644 --- a/test/reactions.jl +++ b/test/reactions.jl @@ -139,4 +139,5 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r)) @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1])) @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1])) +@test isnan(SBMLToolkit.get_massaction(s1 + ifelse(s1 < k1, s1, k1), [s1], [1])) @test_throws ErrorException SBMLToolkit.get_massaction(k1, nothing, [1]) From 5b52c815a00b527981fc00a35f850c752b4d1127 Mon Sep 17 00:00:00 2001 From: paulflang Date: Tue, 10 Dec 2024 23:01:53 +0100 Subject: [PATCH 8/9] fix typo --- src/rules.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rules.jl b/src/rules.jl index 6f609fe..f942a66 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -34,7 +34,7 @@ extensive_kinetic_math(m::SBML.Model, formula::SBML.Math) = SBML.interpret_math( if isnothing(m.compartments[sp.compartment].size) && !seemsdefined(sp.compartment, m) if m.compartments[sp.compartment].spatial_dimensions == 0 - # If the comparment ID doesn't seem directly defined anywhere + # If the compartment ID doesn't seem directly defined anywhere # and it is a zero-dimensional unsized compartment, just avoid # any sizing questions. return x From 751c5c55f15a7fd51b447537ea06245537e93eaf Mon Sep 17 00:00:00 2001 From: paulflang Date: Tue, 10 Dec 2024 23:04:51 +0100 Subject: [PATCH 9/9] format --- src/rules.jl | 80 ++++++++++++++++++++++++++-------------------------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/src/rules.jl b/src/rules.jl index f942a66..9b506e7 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -21,49 +21,49 @@ function get_rules(model) algeqs, obseqs, raterules end - -extensive_kinetic_math(m::SBML.Model, formula::SBML.Math) = SBML.interpret_math( # TODO: move this to SBML.jl - formula, - map_apply = (x, rec) -> SBML.MathApply(x.fn, rec.(x.args)), - map_const = identity, - map_ident = (x::SBML.MathIdent) -> begin - haskey(m.reactions, x.id) && return m.reactions[x.id].kinetic_math - haskey(m.species, x.id) || return x - sp = m.species[x.id] - sp.only_substance_units && return x - if isnothing(m.compartments[sp.compartment].size) && - !seemsdefined(sp.compartment, m) - if m.compartments[sp.compartment].spatial_dimensions == 0 - # If the compartment ID doesn't seem directly defined anywhere - # and it is a zero-dimensional unsized compartment, just avoid - # any sizing questions. - return x - else - # In case the compartment is expected to be defined, complain. - throw( - DomainError( +function extensive_kinetic_math(m::SBML.Model, formula::SBML.Math) + SBML.interpret_math( # TODO: move this to SBML.jl + formula, + map_apply = (x, rec) -> SBML.MathApply(x.fn, rec.(x.args)), + map_const = identity, + map_ident = (x::SBML.MathIdent) -> begin + haskey(m.reactions, x.id) && return m.reactions[x.id].kinetic_math + haskey(m.species, x.id) || return x + sp = m.species[x.id] + sp.only_substance_units && return x + if isnothing(m.compartments[sp.compartment].size) && + !seemsdefined(sp.compartment, m) + if m.compartments[sp.compartment].spatial_dimensions == 0 + # If the compartment ID doesn't seem directly defined anywhere + # and it is a zero-dimensional unsized compartment, just avoid + # any sizing questions. + return x + else + # In case the compartment is expected to be defined, complain. + throw( + DomainError( sp.compartment, - "compartment size is insufficiently defined", + "compartment size is insufficiently defined" ), - ) + ) + end + else + # Now we are sure that the model either has the compartment with + # constant size, or the definition is easily reachable. So just use + # the compartment ID as a variable to compute the concentration (or + # area-centration etc, with different dimensionalities) by dividing + # it. + return SBML.MathApply("/", [x, SBML.MathIdent(sp.compartment)]) end - else - # Now we are sure that the model either has the compartment with - # constant size, or the definition is easily reachable. So just use - # the compartment ID as a variable to compute the concentration (or - # area-centration etc, with different dimensionalities) by dividing - # it. - return SBML.MathApply("/", [x, SBML.MathIdent(sp.compartment)]) - end - end, - map_lambda = (x, _) -> error( - ErrorException("converting lambdas to extensive kinetic math is not supported"), - ), - map_time = identity, - map_avogadro = identity, - map_value = identity, -) - + end, + map_lambda = (x, _) -> error( + ErrorException("converting lambdas to extensive kinetic math is not supported"), + ), + map_time = identity, + map_avogadro = identity, + map_value = identity + ) +end function get_var_and_assignment(model, rule) if !haskey(merge(model.species, model.compartments, model.parameters), rule.variable)