up

Author: TorkelE

docs/src/api/catalyst_api.md

@@ -180,7 +180,7 @@ netstoichmat
 reactionrates
 ```
 
-## Functions to extend or modify a network
+## (Functions to extend or modify a network)(@id api_network_extension_and_modification)
 `ReactionSystem`s can be programmatically extended using
 [`@add_reactions`](@ref), [`addspecies!`](@ref), [`addparam!`](@ref) and
 [`addreaction!`](@ref), or using [`ModelingToolkit.extend`](@ref) and

docs/src/catalyst_applications/bifurcation_diagrams.md

@@ -110,6 +110,10 @@ This bifurcation diagram does not contain any interesting features (such as bifu
 - When providing the concentrations for computing the conserved quantities (in `u0`), we only have to designate the concentrations of species that are actually involved in conservation laws. For larger systems, determining which one are may, however, be difficult. In this case, it might still be wise to provide concentrations for all species.
 - The steady state guess in `u_guess` does not actually have to fulfil the conserved concentrations provided in `u0`.
 
+!!! note
+    Computing bifurcation diagrams for [hierarchical models created using compostability](@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)
 If you use this functionality in your research, please cite the following paper to support the author of the BifurcationKit package:

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -10,7 +10,7 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
     (u0 isa Vector{Pair{Symbol, Float64}}) && (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.

src/miscellaneous.jl

@@ -2,8 +2,10 @@
 
 
 # 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)) || 

test/extensions/bifurcation_kit.jl

@@ -80,4 +80,62 @@ let
         p = rand(rng, 4)
         @test bprob_BK.VF.F(u0, p) == bprob.VF.F(u0, p)
     end
+end
+
+# Creates a system where rn is composed of 4, somewhat nested, networks.
+# Tests with defaults within nested networks.
+let
+    rn1 = @reaction_network rn1 begin
+        @parameters p=1.0
+        (p, d), 0 <--> X
+    end
+    rn2 = @reaction_network rn2 begin
+        @parameters p=2.0
+        (p, d), 0 <--> X
+    end
+    rn3 = @reaction_network rn3 begin
+        @parameters p=3.0
+        (p, d), 0 <--> X
+    end
+    rn4 = @reaction_network rn4 begin
+        @parameters p=4.0
+        (p, d), 0 <--> X
+    end
+    @named rn3 =compose(rn3, [rn4])
+    @named rn = compose(rn1, [rn2, rn3])
+
+    # Declares parameter values and initial u guess.
+    @unpack X, d = rn
+    u0_guess = [X => 1.0, rn2.X => 1.0, rn3.X => 1.0, rn3.rn4.X => 1.0]
+    p_start = [d => 1.0, rn2.d => 1.0, rn3.d => 1.0, rn3.rn4.d => 1.0]
+
+
+    # Computes bifurcation diagrams and check them (the final point at [p_value = 6] should be =[p/d]).
+    p_span = (0.1, 6.0)
+
+    let
+        # Checks top layer.
+        bprob = BifurcationProblem(rn, u0_guess, p_start, d; plot_var = X)
+        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)
+        @test bif_dia.γ.branch[end].x ≈ 1.0/6
+
+        # Checks second layer (1).
+        bprob = BifurcationProblem(rn, u0_guess, p_start, rn2.d; plot_var = rn2.X)
+        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)
+        @test bif_dia.γ.branch[end].x ≈ 2.0/6
+
+        # Checks second layer (2).
+        bprob = BifurcationProblem(rn, u0_guess, p_start, rn3.d; plot_var = rn3.X)
+        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)
+        @test bif_dia.γ.branch[end].x ≈ 3.0/6
+
+        # Checks third layer.
+        bprob = BifurcationProblem(rn, u0_guess, p_start, rn3.rn4.d; plot_var = rn3.rn4.X)
+        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)
+        @test bif_dia.γ.branch[end].x ≈ 4.0/6
+    end
 end

test/programmatic_model_creation/component_based_model_creation.jl

@@ -451,3 +451,33 @@ let
     @test numspecies(rs123) == 6
     @test issetequal(nonspecies(rs123), [V1, rs23.V2, rs23.rs3.V3])
 end
+
+# Tests that conversion with defaults works for a composed model.
+let 
+    rn1 = @reaction_network rn1 begin
+        @parameters p=1.0 r=2.0
+        @species X(t) = 3.0 Y(t) = 4.0
+        (p1, d), 0 <--> X
+        (p2, r), 0 <--> Z
+    end
+    rn2 = @reaction_network rn1 begin
+        @parameters p=10. q=20.0
+        @species X(t) = 30.0 Z(t) = 40.0
+        (p1, d), 0 <--> X
+        (p2, q), 0 <--> Z
+    end
+    composed_reaction_system = compose(rn1, [rn2])
+    osys = convert(ODESystem, composed_reaction_system)
+    parameters(osys)[1].metadata
+    
+    osys.defaults
+    @unpack p, r, X, Y = rn1
+    osys.defaults[p] == 1.0
+    osys.defaults[r] == 2.0
+    osys.defaults[X] == 3.0
+    osys.defaults[Y] == 4.0
+    osys.defaults[rn2.p] == 10.0
+    osys.defaults[rn2.q] == 20.0
+    osys.defaults[rn2.X] == 30.0
+    osys.defaults[rn2.Z] == 40.0
+end