Fixes extension support, add more test, fix conversion bug

Author: TorkelE

Project.toml

@@ -6,6 +6,7 @@ version = "13.5.1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"

ext/CatalystHomotopyContinuationExtension.jl

@@ -2,6 +2,7 @@ module CatalystHomotopyContinuationExtension
 
 # Fetch packages.
 using Catalyst
+import DynamicPolynomials
 import ModelingToolkit as MT
 import HomotopyContinuation as HC
 import Setfield: @set

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -53,7 +53,8 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
     eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
     eqs_intexp = make_int_exps.(eqs)
-    return Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+    ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+    return poly_type_convert(ss_poly)
 end
 
 # If u0s are not given while conservation laws are present, throws an error.
@@ -94,7 +95,7 @@ end
 function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
     var_names_extended = String.(Symbol.(HC.variables(ss_poly)))
     var_names = [Symbol(s[1:prevind(s, findlast('_', s))]) for s in var_names_extended]
-    sort_pattern = indexin(MT.getname.(species(rs)), var_names)
+    sort_pattern = indexin(MT.getname.(unknowns(rs)), var_names)
     foreach(sol -> permute!(sol, sort_pattern), sols)
 end
 
@@ -104,4 +105,13 @@ function filter_negative_f(sols; neg_thres=-1e-20)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end
     return filter(sol -> all(>=(0), sol), sols)
+end
+
+# Sometimes (when polynomials are created from hybrid CRN/DAEs), the steady state polynomial have the wrong type.
+# This converts it to the correct type, which homotopy continuation can handle.
+const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
+const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
+function poly_type_convert(ss_poly)
+    (typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
+    return ss_poly
 end

src/reactionsystem.jl

@@ -993,7 +993,7 @@ function get_indep_sts(rs::ReactionSystem, remove_conserved = false)
     sts = get_unknowns(rs)
     nps = get_networkproperties(rs)
     indepsts = if remove_conserved
-        filter(s -> (s ∈ nps.indepspecs) && (!isbc(s)), sts)
+        filter(s -> ((s ∈ nps.indepspecs) || (!isspecies(s))) && (!isbc(s)), sts)
     else
         filter(s -> !isbc(s), sts)
     end
@@ -1573,9 +1573,9 @@ 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, zero_derivatives=true)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, zero_derivatives=true)
 
-    ODESystem(eqs, get_iv(fullrs), sts, ps;
+    ODESystem(eqs, get_iv(fullrs), us, ps;
               observed = obs,
               name,
               defaults = _merge(defaults,defs),
@@ -1614,10 +1614,10 @@ 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)
-    eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
-    NonlinearSystem(eqs, sts, ps;
+    NonlinearSystem(eqs, us, ps;
                     name,
                     observed = obs,
                     defaults = _merge(defaults,defs),
@@ -1669,13 +1669,13 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     eqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, include_zero_odes,
                          remove_conserved)
     noiseeqs = assemble_diffusion(flatrs, ists, ispcs; combinatoric_ratelaws, remove_conserved)
-    eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved)
 
     if any(isbc, get_unknowns(flatrs))
         @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
     end
 
-    SDESystem(eqs, noiseeqs, get_iv(flatrs), sts, ps;
+    SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
               observed = obs,
               name,
               defaults = defs,

test/extensions/bifurcation_kit.jl

@@ -180,6 +180,29 @@ end
 
 ### Other Tests ###
 
+# Checks that bifurcation diagrams can be computed for hybrid CRN/DAE systems.
+let
+    # Prepares the model (production/degradation of X, with equations for volume and X concentration).
+    rs = @reaction_network begin
+        @parameters k
+        @variables C(t)
+        @equations begin
+            D(V) ~ k*X - V
+            C ~ X/V
+        end
+        (p/V,d/V), 0 <--> X
+    end
+    u0_guess = [:X => 1.0, :V => 1.0, :C => 1.0]
+    p_start = [:p => 2.0, :d => 1.0, :k => 5.0]
+
+    # Computes bifurcation diagram.
+    bprob = BifurcationProblem(rs, u0_guess, p_start, :p; plot_var = :C)
+    p_span = (0.1, 6.0)
+    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 all(getfield.(bif_dia.γ.branch, :x) .≈ 0.2)
+end
+
 # Checks that `BifurcationProblem`s cannot be generated from non-complete `ReactionSystems`s.
 let 
     # Create model.

test/extensions/homotopy_continuation.jl

@@ -109,6 +109,25 @@ end
 
 ### Other Tests ###
 
+# Tests that homotopy continuation can be applied to hybrid DAE/CRN systems.
+let
+    # Prepares the model (production/degradation of X, with equations for volume and X concentration).
+    rs = @reaction_network begin
+        @parameters k
+        @variables C(t)
+        @equations begin
+            D(V) ~ k*X - V
+            C ~ X/V
+        end
+        (p/V,d/V), 0 <--> X
+    end
+
+    # Checks that homotopy continuation correctly find the system's single steady state.
+    ps = [:p => 2.0, :d => 1.0, :k => 5.0]
+    hc_ss = hc_steady_states(rs, ps)
+    @test hc_ss ≈ [[2.0, 0.2, 10.0]]
+end
+
 # Checks that `hc_steady_states` cannot be applied to non-complete `ReactionSystems`s.
 let 
     # Create model.

test/extensions/structural_identifiability.jl

@@ -312,6 +312,39 @@ end
 
 ### Other Tests ###
 
+# Checks that identifiability can be assessed for hybrid CRN/DAE systems.
+let
+    rs = @reaction_network begin
+        @parameters k c1 c2
+        @variables C(t)
+        @equations begin
+            D(V) ~ k*X - V
+            C ~ (c1 + c2) * X/V
+        end
+        (p/V,d/V), 0 <--> X
+    end
+    @unpack p, d, k, c1, c2 = rs
+    
+    # Tests identifiability assessment when all unknowns are measured.
+    gi_1 = assess_identifiability(rs; measured_quantities=[:X, :V, :C])
+    li_1 = assess_local_identifiability(rs; measured_quantities=[:X, :V, :C])
+    ifs_1 = find_identifiable_functions(rs; measured_quantities=[:X, :V, :C])
+    @test sym_dict(gi_1) == Dict([:X => :globally, :C => :globally, :V => :globally, :k => :globally, 
+                      :c1 => :nonidentifiable, :c2 => :nonidentifiable, :p => :globally, :d => :globally])
+    @test sym_dict(li_1) == Dict([:X => 1, :C => 1, :V => 1, :k => 1, :c1 => 0, :c2 => 0, :p => 1, :d => 1])
+    @test issetequal(ifs_1, [d, p, k, c1 + c2])
+    
+    # Tests identifiability assessment when only variables are measured. 
+    # Checks that a parameter in an equation can be set as known.
+    gi_2 = assess_identifiability(rs; measured_quantities=[:V, :C], known_p = [:c1])
+    li_2 = assess_local_identifiability(rs; measured_quantities=[:V, :C], known_p = [:c1])
+    ifs_2 = find_identifiable_functions(rs; measured_quantities=[:V, :C], known_p = [:c1])
+    @test sym_dict(gi_2) == Dict([:X => :nonidentifiable, :C => :globally, :V => :globally, :k => :nonidentifiable, 
+                      :c1 => :globally, :c2 => :nonidentifiable, :p => :nonidentifiable, :d => :globally])
+    @test sym_dict(li_2) == Dict([:X => 0, :C => 1, :V => 1, :k => 0, :c1 => 1, :c2 => 0, :p => 0, :d => 1])
+    @test issetequal(ifs_2, [d, c1, k*p, c1*p + c2*p])
+end
+
 # Checks that identifiability functions cannot be applied to non-complete `ReactionSystems`s.
 let 
     # Create model.