Merge pull request #1237 from SciML/species_diffeq_error

Error when having differentials w.r.t. species

Author: TorkelE

src/reactionsystem.jl

@@ -367,6 +367,13 @@ struct ReactionSystem{V <: NetworkProperties} <:
             end
         end
 
+        # Checks that no (non-reaction) equation contains a differential w.r.t. a species.
+        for eq in eqs
+            (eq isa Reaction) && continue
+            (hasnode(is_species_diff, eq.lhs) || hasnode(is_species_diff, eq.rhs)) &&
+                error("An equation ($eq) contains a differential with respect to a species. This is currently not supported. If this is a functionality you require, please raise an issue on the Catalyst GitHub page and we can consider the best way to implement it.")
+        end
+
         rs = new{typeof(nps)}(
             eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed,
             name, systems, defaults, connection_type, nps, cls, cevs,
@@ -376,6 +383,13 @@ struct ReactionSystem{V <: NetworkProperties} <:
     end
 end
 
+# Checks if a symbolic expression constains a differential with respect to a species (either directly
+# or somehwere within the differential expression).
+function is_species_diff(expr)
+    Symbolics.is_derivative(expr) || return false
+    return hasnode(ex -> (ex isa Symbolics.BasicSymbolic) && isspecies(ex) && !isbc(ex), expr)
+end
+
 # Four-argument constructor. Permits additional inputs as optional arguments.
 # Calls the full constructor.
 function ReactionSystem(eqs, iv, unknowns, ps;
@@ -482,7 +496,7 @@ function ReactionSystem(rxs::Vector, iv = Catalyst.DEFAULT_IV; kwargs...)
     make_ReactionSystem_internal(rxs, iv, [], []; kwargs...)
 end
 
-# One-argument constructor. Creates an emtoy `ReactionSystem` from a time independent variable only.
+# One-argument constructor. Creates an empty `ReactionSystem` from a time independent variable only.
 function ReactionSystem(iv; kwargs...)
     ReactionSystem(Reaction[], iv, [], []; kwargs...)
 end

test/reactionsystem_core/coupled_equation_crn_systems.jl

@@ -978,17 +978,6 @@ let
     @variables V1(t)
     @species S1(t) S2(t)
 
-    # Coupled system with additional differential equation for species.
-    eqs = [
-        Reaction(p1, [S1], [S2]),
-        D(S1) ~ p2 - S1
-    ]
-    @named rs = ReactionSystem(eqs, t)
-    rs = complete(rs)
-    u0 = [S1 => 1.0, S2 => 2.0]
-    ps = [p1 => 2.0, p2 => 3.0]
-    @test_throws Exception ODEProblem(rs, u0, (0.0, 1.0), ps; structural_simplify = true)
-
     # Coupled system overconstrained due to additional algebraic equations (without variables).
     eqs = [
         Reaction(p1, [S1], [S2]),
@@ -1012,3 +1001,32 @@ let
     ps = [p1 => 2.0, p2 => 3.0]
     @test_throws Exception ODEProblem(rs, u0, (0.0, 1.0), ps; structural_simplify = true)
 end
+
+# Checks that equations cannot contain differentials with respect to species.
+let
+    # Basic case.
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(X) ~ 1.0
+        d, X --> 0
+    end
+
+    # Complicated differential.
+    @test_throws Exception @eval @reaction_network begin        
+        @equations V + X^2 ~ 1.0 - D(V + log(1 + X))
+        d, X --> 0
+    end
+
+    # Case where the species is declared, but not part of a reaction.
+    @test_throws Exception @eval @reaction_network begin
+        @species Y(t)
+        @equations D(Y) ~ 1.0
+        d, X --> 0
+    end
+
+    # Case where the equation also declares a new non-species variable.
+    # At some point something like this could be supported, however, not right now.
+    @test_throws Exception @eval @reaction_network begin
+        @equations D(V) ~ 1.0 + D(X)
+        d, X --> 0
+    end
+end

test/reactionsystem_core/reactionsystem.jl

@@ -770,21 +770,6 @@ let
 end
 
 # Fix for SBML test 305.
-let
-    @parameters k1 k2 S2 [isconstantspecies = true]
-    @species S1(t) S3(t)
-    rx = Reaction(k2, [S1], nothing)
-    ∂ₜ = default_time_deriv()
-    eq = ∂ₜ(S3) ~ k1 * S2
-    @mtkbuild osys = ODESystem([eq], t)
-    @named rs = ReactionSystem([rx, eq], t)
-    rs = complete(rs)
-    @test issetequal(unknowns(rs), [S1, S3])
-    @test issetequal(parameters(rs), [S2, k1, k2])
-    osys = convert(ODESystem, rs)
-    @test issetequal(unknowns(osys), [S1, S3])
-    @test issetequal(parameters(osys), [S2, k1, k2])
-end
 let
     @parameters k1 k2 S2 [isconstantspecies = true]
     @species S1(t) S3(t) [isbcspecies = true]