Merge pull request #1222 from SciML/remove_rates_only_bundling

Throw error when bundling reactions where only the rates differe

Author: TorkelE

HISTORY.md

@@ -76,6 +76,21 @@
   functions (as a parameter) in a model. For more details on how to use these,
   please read:
   https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
+- We have introduced a restriction on bundling of reactions in the DSL. Now,
+  bundling is not permitted if multiple rates are provided but only one set each
+  of substrates/products. E.g. this model:
+ ```julia
+  @reaction_network begin
+    (k1,k2), X --> Y
+  end
+ ```
+  will now throw an error. The reason that users attempting to write bi-directional 
+  reactions but typing `-->` instead of `<-->` would get a wrong model. We decided that
+  this kind of bundling was unlikely to be used, and throwing errors for people who
+  made the typo was more important. 
+
+  If you use this type of bundling and it indeed is useful to you, please raise and issue
+  and we will see if we can sort something out.
 - Scoped species/variables/parameters are now treated similar to the latest MTK
   releases (≥ 9.49).
 - A tutorial on making interactive plot displays using Makie has been added.

docs/src/faqs.md

@@ -339,7 +339,7 @@ conserved constant(s), `Γ`, are updated. As an example consider the following
 using Catalyst, NonlinearSolve
 rn = @reaction_network begin
     (k₁,k₂), X₁ <--> X₂
-    (k₃,k₄), X₁ + X₂ --> 2X₃
+    (k₃,k₄), X₁ + X₂ <--> 2X₃
 end
 u0 = [:X₁ => 1.0, :X₂ => 2.0, :X₃ => 3.0]
 ps = [:k₁ => 0.1, :k₂ => 0.2, :k₃ => 0.3, :k₄ => 0.4]

docs/src/model_creation/dsl_basics.md

@@ -188,6 +188,15 @@ end
 ```
 However, like for the above model, bundling reactions too zealously can reduce (rather than improve) a model's readability.
 
+The one exception to reaction bundling is that we do not permit the user to provide multiple rates but only set one set each for the substrates and products. I.e.
+```julia
+rn = @reaction_network begin
+    (k1,k2), X --> Y
+end
+```
+is not permitted (due to this notation's similarity to a bidirectional reaction). However, if multiples are provided for substrates and/or products, like `(k1,k2), (X1,X2) --> Y`, then bundling works.
+
+
 ## [Non-constant reaction rates](@id dsl_description_nonconstant_rates)
 So far we have assumed that all reaction rates are constant (being either a number of a parameter). Non-constant rates that depend on one (or several) species are also possible. More generally, the rate can be any valid expression of parameters and species.
 

src/dsl.jl

@@ -365,6 +365,13 @@ function get_reactions(exprs::Vector{Expr})
         # Reads core reaction information.
         arrow, rate, reaction, metadata = read_reaction_line(line)
 
+        # Currently, reaction bundling where rates (but neither substrates nor products) are
+        # bundled, is disabled. See discussion in https://github.com/SciML/Catalyst.jl/issues/1219.
+        if !in(arrow, double_arrows) && Meta.isexpr(rate, :tuple) &&
+                !Meta.isexpr(reaction.args[2], :tuple) && !Meta.isexpr(reaction.args[3], :tuple)
+            error("Bundling of reactions with multiple rates but singular substrates and product sets is disallowed. This error is potentially due to a bidirectional (`<-->`) reaction being incorrectly typed as `-->`.")
+        end
+
         # Checks which type of line is used, and calls `push_reactions!` on the processed line.
         if in(arrow, double_arrows)
             (typeof(rate) != Expr || rate.head != :tuple) &&

test/dsl/dsl_basic_model_construction.jl

@@ -179,7 +179,7 @@ let
     different_arrow_8 = @reaction_network begin
         p, 2X1 < ∅
         k1, X2 ← X1
-        (k2, k3), X3 ⟻ X2
+        (k2, k3), X3 ⟻ (X2,X2)
         d, ∅ ↼ X3
     end
     push!(identical_networks_1, reaction_networks_standard[8] => different_arrow_8)
@@ -298,10 +298,11 @@ let
 
     no_parameters_10 = @reaction_network begin
         0.01, ∅ ⟶ X1
-        (3.1, 3.2), X1 → X2
-        (0.0, 2.1), X2 → X3
-        (901.0, 63.5), X3 → X4
-        (7, 8), X4 → X5
+        (3.1, 3.2), (X1,X1) → X2
+        (0.0, 2.1), X2 → (X3,X3)
+        (901.0, 63.5), (X3,X3) → (X4,X4)
+        7, X4 → X5
+        8, X4 → X5
         1.0, X5 ⟶ ∅
     end
     push!(identical_networks_3, reaction_networks_standard[10] => no_parameters_10)
@@ -591,3 +592,20 @@ let
         d, X ⇻ 0
     end
 end
+
+# Tests that bundling, but where the rate only have been given multiple alternatives, errors.
+# Checks for two, more than two, and for functional rates. Uses different arrow types and directions.
+let
+    @test_throws Exception @eval @reaction_network begin
+        (k1,k2), X --> 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        (k1,k2,k3), X → 0
+    end
+    @test_throws Exception @eval @reaction_network begin
+        (k1,k1), X => Y
+    end
+    @test_throws Exception @eval @reaction_network begin
+        (mm(X,v1,K1),mm(X,v2,K2)), 0 <-- 0
+    end
+end

test/dsl/dsl_options.jl

@@ -316,7 +316,8 @@ let
     @parameters p1=1.0 p2=2.0 k1=4.0 k2=5.0 v=8.0 K=9.0 n=3 d=10.0
     @species X(t)=4.0 Y(t)=3.0 X2Y(t)=2.0 Z(t)=1.0
         (p1,p2), 0 --> (X,Y)
-        (k1,k2), 2X + Y --> X2Y
+        k1, 2X + Y --> X2Y
+        k2, 2X + Y --> X2Y
         hill(X2Y,v,K,n), 0 --> Z
         d, (X,Y,X2Y,Z) --> 0
     end
@@ -327,7 +328,8 @@ let
     @parameters p1=1.0 p2 k1=4.0 k2 v=8.0 K n=3 d
     @species X(t)=4.0 Y(t) X2Y(t) Z(t)=1.0
         (p1,p2), 0 --> (X,Y)
-        (k1,k2), 2X + Y --> X2Y
+        k1, 2X + Y --> X2Y
+        k2, 2X + Y --> X2Y
         hill(X2Y,v,K,n), 0 --> Z
         d, (X,Y,X2Y,Z) --> 0
     end
@@ -339,7 +341,8 @@ let
     @parameters p1 p2 k1 k2 v K n d
     @species X(t) Y(t) X2Y(t) Z(t)
         (p1,p2), 0 --> (X,Y)
-        (k1,k2), 2X + Y --> X2Y
+        k1, 2X + Y --> X2Y
+        k2, 2X + Y --> X2Y
         hill(X2Y,v,K,n), 0 --> Z
         d, (X,Y,X2Y,Z) --> 0
     end

test/network_analysis/network_properties.jl

@@ -76,7 +76,7 @@ let
         (k₁, k₂), E + S1 <--> ES1
         (k₃, k₄), E + S2 <--> ES2
         (k₅, k₆),  S2 + ES1 <--> ES1S2
-        (k₆, k₇), ES1S2 --> S1 + ES2
+        (k₆, k₇), ES1S2 <--> S1 + ES2
         k₈, ES1S2 --> E+P
         (k₉, k₁₀), S1 <--> 0
         (k₁₀, k₁₁), 0 <--> S2
@@ -148,18 +148,18 @@ let
     @test length(slcs) == 3
     @test length(tslcs) == 2
     @test issubset([[1,2], [3,4,5], [6,7]], slcs)
-    @test issubset([[3,4,5], [6,7]], tslcs) 
+    @test issubset([[3,4,5], [6,7]], tslcs)
 end
 
-# b) Makes the D + E --> G reaction irreversible. Thus, (D+E) becomes a non-terminal linkage class. Checks whether correctly identifies both (A, B+C) and (D+E) as non-terminal 
+# b) Makes the D + E --> G reaction irreversible. Thus, (D+E) becomes a non-terminal linkage class. Checks whether correctly identifies both (A, B+C) and (D+E) as non-terminal
 let
     rn = @reaction_network begin
         (k1, k2), A <--> B + C
         k3, B + C --> D
         k4, D --> E
         (k5, k6), E <--> 2F
         k7, 2F --> D
-        (k8, k9), D + E --> G
+        k8, D + E --> G
     end
 
     rcs, D = reactioncomplexes(rn)
@@ -168,10 +168,10 @@ let
     @test length(slcs) == 4
     @test length(tslcs) == 2
     @test issubset([[1,2], [3,4,5], [6], [7]], slcs)
-    @test issubset([[3,4,5], [7]], tslcs) 
+    @test issubset([[3,4,5], [7]], tslcs)
 end
 
-# From a), makes the B + C <--> D reaction reversible. Thus, the non-terminal (A, B+C) linkage class gets absorbed into the terminal (A, B+C, D, E, 2F) linkage class, and the terminal linkage classes and strong linkage classes coincide. 
+# From a), makes the B + C <--> D reaction reversible. Thus, the non-terminal (A, B+C) linkage class gets absorbed into the terminal (A, B+C, D, E, 2F) linkage class, and the terminal linkage classes and strong linkage classes coincide.
 let
     rn = @reaction_network begin
         (k1, k2), A <--> B + C
@@ -188,13 +188,13 @@ let
     @test length(slcs) == 2
     @test length(tslcs) == 2
     @test issubset([[1,2,3,4,5], [6,7]], slcs)
-    @test issubset([[1,2,3,4,5], [6,7]], tslcs) 
+    @test issubset([[1,2,3,4,5], [6,7]], tslcs)
 end
 
 # Simple test for strong and terminal linkage classes
 let
     rn = @reaction_network begin
-        (k1, k2), A <--> 2B 
+        (k1, k2), A <--> 2B
         k3, A --> C + D
         (k4, k5), C + D <--> E
         k6, 2B --> F
@@ -209,7 +209,7 @@ let
     @test length(slcs) == 3
     @test length(tslcs) == 2
     @test issubset([[1,2], [3,4], [5,6,7]], slcs)
-    @test issubset([[3,4], [5,6,7]], tslcs) 
+    @test issubset([[3,4], [5,6,7]], tslcs)
 end
 
 # Cycle Test: Open Reaction Network
@@ -261,7 +261,7 @@ let
         k2, C + D --> E + F
         k3, C + D --> 2G + H
         k4, 2G + H --> 3I
-        k5, E + F --> J 
+        k5, E + F --> J
         k6, 3I --> K
     end
 
@@ -273,7 +273,7 @@ end
 
 ### Other Network Properties Tests ###
 
-# Tests outgoing complexes matrices (1). 
+# Tests outgoing complexes matrices (1).
 # Checks using dense and sparse representation.
 let
     # Declares network.
@@ -283,7 +283,7 @@ let
         k3, X1 --> X2
         k4, X1 + X2 --> X2
     end
-    
+
     # Compares to manually computed matrix.
     cmplx_out_mat = [
         -1 0 0 0 -1;
@@ -368,17 +368,17 @@ let
                MAPK.KKPP * MAPK.KKPase,
                MAPK.KKPP * MAPK.K,
                MAPK.KKPPK,
-               MAPK.KKPP * MAPK.KP, 
+               MAPK.KKPP * MAPK.KP,
                MAPK.KPKKPP,
                MAPK.KPP * MAPK.KKPP,
                MAPK.KP * MAPK.KPase,
-               MAPK.KPKPase, 
-               MAPK.KKPPKPase, 
+               MAPK.KPKPase,
+               MAPK.KKPPKPase,
                MAPK.K * MAPK.KPase,
                MAPK.KPP * MAPK.KPase,
               ]
     @test isequal(Φ, truevec)
-    
+
     K = Catalyst.fluxmat(MAPK)
     # Construct flux matrix from incidence matrix
     mat = Matrix{Any}(zeros(30, 26))
@@ -391,12 +391,12 @@ let
     @test isequal(K, mat)
     @test isequal(K[1, 1], MAPK.k₁)
     @test all(==(0), K[1, 2:end])
-    @test isequal(K[2, 2], MAPK.k₂) 
+    @test isequal(K[2, 2], MAPK.k₂)
     @test all(==(0), vcat(K[2,1], K[2,3:end]))
     @test isequal(K[3, 2], MAPK.k₃)
     @test all(==(0), vcat(K[3,1], K[3,3:end]))
     @test count(k -> !isequal(k, 0), K) == length(reactions(MAPK))
-    
+
     A_k = Catalyst.laplacianmat(MAPK)
     @test all(col -> sum(col) == 0, eachcol(A_k))
 
@@ -415,7 +415,7 @@ let
     ratetup = Tuple(ratevec)
 
     @test Catalyst.fluxmat(MAPK, ratemap) == Catalyst.fluxmat(MAPK, ratevec) == Catalyst.fluxmat(MAPK, ratetup)
-    
+
     K = Catalyst.fluxmat(MAPK, ratemap)
     A_k = Catalyst.laplacianmat(MAPK, ratemap)
     @test all(col -> sum(col) == 0, eachcol(A_k))
@@ -428,7 +428,7 @@ let
     @test all(iszero, simplify(numeqs - Y*A_k*Φ))
 end
 
-# Test handling for weird complexes and combinatoric rate laws. 
+# Test handling for weird complexes and combinatoric rate laws.
 let
     rn = @reaction_network begin
         k1, 2X + Y + 3Z --> ∅

test/test_networks.jl

@@ -64,7 +64,8 @@ end
 reaction_networks_standard[8] = @reaction_network rns8 begin
     p, ∅ → 2X1
     k1, X1 → X2
-    (k2, k3), X2 → X3
+    k2, X2 → X3
+    k3, X2 → X3
     d, X3 → ∅
 end
 
@@ -78,10 +79,10 @@ end
 
 reaction_networks_standard[10] = @reaction_network rns10 begin
     p, ∅ ⟶ X1
-    (k1, k2), X1 → X2
-    (k3, k4), X2 → X3
-    (k5, k6), X3 → X4
-    (k7, k8), X4 → X5
+    (k1, k2), (X1,X1) → X2
+    (k3, k4), (X2,X2) → X3
+    (k5, k6), (X3,X3) → X4
+    (k7, k8), (X4,X4) → X5
     d, X5 ⟶ ∅
 end