up

Author: TorkelE

docs/src/catalyst_applications/bifurcation_diagrams.md

@@ -112,9 +112,6 @@ Finally, for additional clarity, we reiterate the purpose of the two `u` argumen
 - `u_guess`: A guess of the initial steady states (which BifurcationKit uses to find its starting point). Typically, most trivial guesses work (e.g. setting all species concentrations to `1.0`). `u_guess` *does not* have to fulfill the conserved concentrations provided in `u0`.
 - `u0`: Used to compute the concentrations of any conserved quantities (e.g. in our example $X + Xp = 1.0$). Technically, values are only required for species that are involved in conservation laws (in our case we do not need to provide a value for $K$). However, sometimes determining which species are actually involved in conservation laws can be difficult, and it might be easier to simply provide concentrations for all species.
 
-!!! note
-    Computing bifurcation diagrams for [hierarchical models created by composing multiple systems](@ref compositional_modeling), that also contain conservations laws, is currently not supported. For these, please [flatten](@ref api_network_extension_and_modification) the system first.
-
 
 ---
 ## [Citation](@id bifurcation_kit_citation)

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -1,18 +1,21 @@
 ### Dispatch for BifurcationKit BifurcationProblems ###
 
 # Creates a BifurcationProblem, using a ReactionSystem as an input.
-function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...; plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)  
+function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...; 
+                               plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)  
+
     # Converts symbols to symbolics.
     (bif_par isa Symbol) && (bif_par = rs.var_to_name[bif_par])
     (plot_var isa Symbol) && (plot_var = rs.var_to_name[plot_var])
-    (u0_bif isa Vector{Pair{Symbol, Float64}}) && (u0_bif = symmap_to_varmap(rs, u0_bif))
-    (ps isa Vector{Pair{Symbol, Float64}}) && (ps = symmap_to_varmap(rs, ps))
-    (u0 isa Vector{Pair{Symbol, Float64}}) && (u0 = symmap_to_varmap(rs, u0))
+    ((u0_bif isa Vector{<:Pair{Symbol,<:Any}}) || (u0_bif isa Dict{Symbol, <:Any})) && (u0_bif = symmap_to_varmap(rs, u0_bif))
+    ((ps isa Vector{<:Pair{Symbol,<:Any}}) || (ps isa Dict{Symbol, <:Any})) && (ps = symmap_to_varmap(rs, ps))
+    ((u0 isa Vector{<:Pair{Symbol,<:Any}}) || (u0 isa Dict{Symbol, <:Any})) && (u0 = symmap_to_varmap(rs, u0))
 
     # Creates NonlinearSystem.
-    #Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
+    Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
     nsys = convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0))
 
-    # Makes BifurcationProblem.
-    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var, record_from_solution=record_from_solution, jac=jac, kwargs...)
+    # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).
+    return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var, 
+                                 record_from_solution=record_from_solution, jac=jac, kwargs...)
 end

src/Catalyst.jl

@@ -106,9 +106,6 @@ export @compound
 export components, iscompound, coefficients
 export balance_reaction
 
-# for functions I am unsure where to best place them.
-include("miscellaneous.jl")
-
 ### Extensions ###
 
 # HomotopyContinuation

src/miscellaneous.jl

@@ -1279,13 +1279,11 @@ conservedquantities(state, cons_laws) = cons_laws * state
 
 
 # If u0s are not given while conservation laws are present, throws an error.
-# If the system is nested, the function is skipped (conservation laws cannot be computed for nested systems).
 # Used in HomotopyContinuation and BifurcationKit extensions.
 function conservationlaw_errorcheck(rs, pre_varmap)
-    isempty(get_systems(rs)) || return
     vars_with_vals = Set(p[1] for p in pre_varmap)
     any(s -> s in vars_with_vals, species(rs)) && return
-    isempty(conservedequations(rs)) || 
+    isempty(conservedequations(Catalyst.flatten(rs))) || 
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
 

src/networkapi.jl

@@ -138,4 +138,36 @@ let
         bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
         @test bif_dia.γ.branch[end].x ≈ 4.0/6
     end
+end
+
+# Tests for nested model with conservation laws.
+let
+    # Creates model.
+    rn1 = @reaction_network rn1 begin
+        (k1, k2), X1 <--> X2
+    end
+    rn2 = @reaction_network rn2 begin
+        (l1, l2), Y1 <--> Y2
+    end
+    @named rn = compose(rn1, [rn2])
+
+    # Creates input parameter and species vectors.
+    @unpack X1, X2, k1, k2 = rn1
+    u_guess = [X1 => 1.0, X2 => 1.0, rn2.Y1 => 1.0, rn2.Y2 => 1.0]
+    p_start = [k1 => 1.0, k2 => 1.0, rn2.l1 => 1.0, rn2.l2 => 1.0]
+    u0 = [X1 => 1.0, X2 => 0.0, rn2.Y1 => 1.0, rn2.Y2 => 0.0]
+
+    # Computes bifurcation diagram.
+    p_span = (0.2, 5.0)
+    bprob = BifurcationProblem(rn, u_guess, p_start, k1; plot_var = X1, u0=u0)
+    opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.001, ds = 0.0001, max_steps = 10000, p_min = p_span[1], p_max = p_span[2], n_inversion = 4)
+    bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+
+    # Checks that the bifurcation diagram is correct.
+    xs = getfield.(bif_dia.γ.branch, :x)
+    k1s = getfield.(bif_dia.γ.branch, :param)
+    all(1 ./ k1s .* (1 .- xs) .≈ xs)
+
+    # Checks that there is an error if information for conserved quantities computation is not provided.
+    @test_throws Exception bprob = BifurcationProblem(rn, u_guess, p_start, k1; plot_var = X1)
 end