init

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), 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

@@ -1012,3 +1012,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