Merge pull request #330 from TorkelE/reactions_with_multu_stoich_and_empty_sides

Allow reactions with: Multiple copies of one components and either no substrates or no products

Author: ChrisRackauckas

src/systems/reaction/reactionsystem.jl

@@ -3,25 +3,36 @@ struct Reaction{S <: Variable, T <: Number}
     substrates::Vector{Operation}
     products::Vector{Operation}
     substoich::Vector{T}
-    prodstoich::Vector{T}    
+    prodstoich::Vector{T}
     netstoich::Vector{Pair{S,T}}
     only_use_rate::Bool
 end
 
-function Reaction(rate, subs, prods, substoich, prodstoich; 
+function Reaction(rate, subs, prods, substoich, prodstoich;
                   netstoich=nothing, only_use_rate=false, kwargs...)
 
-    subsv  = isnothing(subs) ? Vector{Operation}() : subs
-    prodsv = isnothing(prods) ? Vector{Operation}() : prods
-    ns = isnothing(netstoich) ? get_netstoich(subsv, prodsv, substoich, prodstoich) : netstoich
-    Reaction(rate, subsv, prodsv, substoich, prodstoich, ns, only_use_rate)
+      (isnothing(prods)&&isnothing(subs)) && error("A reaction requires a non-nothing substrate or product vector.")
+      (isnothing(prodstoich)&&isnothing(substoich)) && error("Both substrate and product stochiometry inputs cannot be nothing.")
+      if isnothing(subs)
+        subs = Vector{Operation}()
+        (substoich!=nothing) && error("If substrates are nothing, substrate stiocihometries have to be so too.")
+        substoich = typeof(prodstoich)()
+    end
+    if isnothing(prods)
+        prods = Vector{Operation}()
+        (prodstoich!=nothing) && error("If products are nothing, product stiocihometries have to be so too.")
+        prodstoich = typeof(substoich)()
+    end
+    ns = isnothing(netstoich) ? get_netstoich(subs, prods, substoich, prodstoich) : netstoich
+    Reaction(rate, subs, prods, substoich, prodstoich, ns, only_use_rate)
 end
 
+
 # three argument constructor assumes stoichiometric coefs are one and integers
-function Reaction(rate, subs, prods; kwargs...) 
+function Reaction(rate, subs, prods; kwargs...)
 
-    sstoich = isnothing(subs) ? Int[] : ones(Int,length(subs))
-    pstoich = isnothing(prods) ? Int[] : ones(Int,length(prods))
+    sstoich = isnothing(subs) ? nothing : ones(Int,length(subs))
+    pstoich = isnothing(prods) ? nothing : ones(Int,length(prods))
     Reaction(rate, subs, prods, sstoich, pstoich; kwargs...)
 end
 
@@ -54,18 +65,18 @@ end
 function ReactionSystem(eqs, iv, species, params; systems = ReactionSystem[],
                                                   name = gensym(:ReactionSystem))
 
-    ReactionSystem(eqs, iv, convert.(Variable,species), convert.(Variable,params), 
+    ReactionSystem(eqs, iv, convert.(Variable,species), convert.(Variable,params),
                    name, systems)
 end
 
 # Calculate the ODE rate law
 function oderatelaw(rx)
-    @unpack rate, substrates, substoich, only_use_rate = rx    
+    @unpack rate, substrates, substoich, only_use_rate = rx
     rl = rate
     if !only_use_rate
         coef = one(eltype(substoich))
         for (i,stoich) in enumerate(substoich)
-            coef *= factorial(stoich)        
+            coef *= factorial(stoich)
             rl   *= isone(stoich) ? substrates[i] : substrates[i]^stoich
         end
         (!isone(coef)) && (rl /= coef)
@@ -78,15 +89,15 @@ function assemble_drift(rs)
     eqs = [D(x(rs.iv())) ~ 0 for x in rs.states]
     species_to_idx = Dict((x => i for (i,x) in enumerate(rs.states)))
 
-    for rx in rs.eqs        
+    for rx in rs.eqs
         rl = oderatelaw(rx)
         for (spec,stoich) in rx.netstoich
             i = species_to_idx[spec]
             if iszero(eqs[i].rhs)
                 signedrl = (stoich > zero(stoich)) ? rl : -rl
-                rhs      = isone(abs(stoich)) ? signedrl : stoich * rl                
+                rhs      = isone(abs(stoich)) ? signedrl : stoich * rl
             else
-                Δspec = isone(abs(stoich)) ? rl : abs(stoich) * rl            
+                Δspec = isone(abs(stoich)) ? rl : abs(stoich) * rl
                 rhs   = (stoich > zero(stoich)) ? (eqs[i].rhs + Δspec) : (eqs[i].rhs - Δspec)
             end
             eqs[i] = Equation(eqs[i].lhs, rhs)
@@ -104,7 +115,7 @@ function assemble_diffusion(rs)
         for (spec,stoich) in rx.netstoich
             i            = species_to_idx[spec]
             signedrlsqrt = (stoich > zero(stoich)) ? rlsqrt : -rlsqrt
-            eqs[i,j]     = isone(abs(stoich)) ? signedrlsqrt : stoich * rlsqrt            
+            eqs[i,j]     = isone(abs(stoich)) ? signedrlsqrt : stoich * rlsqrt
         end
     end
     eqs

test/reactionsystem.jl

@@ -1,11 +1,11 @@
 using ModelingToolkit, LinearAlgebra, Test
 
-@parameters t k[1:13]
+@parameters t k[1:15]
 @variables A(t) B(t) C(t) D(t)
 rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
        Reaction(k[2], [B], nothing),            # B -> 0
        Reaction(k[3],[A],[C]),                  # A -> C
-       Reaction(k[4], [C], [A,B]),              # C -> A + B        
+       Reaction(k[4], [C], [A,B]),              # C -> A + B
        Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
        Reaction(k[6], [A,B], [C]),              # A + B -> C
        Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
@@ -14,8 +14,10 @@ rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
        Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
        Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
        Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]), # A+3B+4C -> 2C + 3D
-       Reaction(k[13]*A/(2+A), [A], nothing; only_use_rate=true)  # A -> 0 with custom rate 
-       ]  
+       Reaction(k[13], [A,B], nothing, [3,1], nothing), # 3A+B -> 0
+       Reaction(k[14], nothing, [A], nothing, [2]), # 0 -> 2A
+       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true)  # A -> 0 with custom rate
+       ]
 rs = ReactionSystem(rxs,t,[A,B,C,D],k)
 odesys = convert(ODESystem,rs)
 sdesys = convert(SDESystem,rs)
@@ -24,8 +26,8 @@ sdesys = convert(SDESystem,rs)
 function oderhs(u,k,t)
        A = u[1]; B = u[2]; C = u[3]; D = u[4];
        du = zeros(eltype(u),4)
-       du[1] = k[1] - k[3]*A + k[4]*C + 2*k[5]*C - k[6]*A*B + k[7]*B^2/2 - k[9]*A*B - k[10]*A^2 - k[11]*A^2/2 - k[12]*A*B^3*C^4/144 - k[13]*A/(2+A)
-       du[2] = -k[2]*B + k[4]*C - k[6]*A*B - k[7]*B^2 - k[8]*A*B - k[9]*A*B + k[11]*A^2/2 - 3*k[12]*A*B^3*C^4/144
+       du[1] = k[1] - k[3]*A + k[4]*C + 2*k[5]*C - k[6]*A*B + k[7]*B^2/2 - k[9]*A*B - k[10]*A^2 - k[11]*A^2/2 - k[12]*A*B^3*C^4/144 - 3*k[13]*A^3*B/6 + 2*k[14] - k[15]*A/(2+A)
+       du[2] = -k[2]*B + k[4]*C - k[6]*A*B - k[7]*B^2 - k[8]*A*B - k[9]*A*B + k[11]*A^2/2 - 3*k[12]*A*B^3*C^4/144 - k[13]*A^3*B/6
        du[3] = k[3]*A - k[4]*C - k[5]*C  + k[6]*A*B + k[8]*A*B + k[9]*A*B + k[10]*A^2/2 - 2*k[12]*A*B^3*C^4/144
        du[4] = k[9]*A*B + k[10]*A^2/2 + 3*k[12]*A*B^3*C^4/144
        du
@@ -49,8 +51,10 @@ function sdenoise(u,k,t)
             -2*sqrt(k[10]*A^2/2) z sqrt(k[10]*A^2/2) sqrt(k[10]*A^2/2);
             -sqrt(k[11]*A^2/2) sqrt(k[11]*A^2/2) z z;
             -sqrt(k[12]*A*B^3*C^4/144) -3*sqrt(k[12]*A*B^3*C^4/144) -2*sqrt(k[12]*A*B^3*C^4/144) 3*sqrt(k[12]*A*B^3*C^4/144);
-            -sqrt(k[13]*A/(2+A)) z z z]'
-       
+            -3*sqrt(k[13]*A^3*B/6) -sqrt(k[13]*A^3*B/6) z z
+            2*sqrt(k[14]) z z z;
+            -sqrt(k[15]*A/(2+A)) z z z]'
+
        return G
 end