Merge pull request #7 from weilandtd/patch-1

Higher order reactions

Author: ChrisRackauckas

src/reactions.jl

@@ -59,8 +59,18 @@ ReactionSet(args...) = ReactionSet(tuple(args...))
 function build_jumps_from_reaction(r::Reaction;save_positions=(false,true))
   rate = function (t,u)
     val = r.rate_constant
+    # for higher order interactions i.e. reactants [1,1] vaule = k*u[1]*(u[1]-1)
+    # Since a single entitiy cannot interact with it self!
+    k = 0.0
+    prev_reactant = -1
     @fastmath @inbounds for i in eachindex(r.reactants)
-      val *= u[r.reactants[i]]
+      if prev_reactant == r.reactants[i]
+        k += 1.0
+      else
+        k = 0.0
+      end
+      val *= (u[r.reactants[i]] - k)
+      prev_reactant = r.reactants[i]
     end
     val
   end
@@ -75,8 +85,18 @@ end
 function build_jumps_from_reaction(r::VariableRateReaction;save_positions=(false,true))
   rate = function (t,u)
     val = r.rate_constant
+    # for higher order interactions i.e. reactants [1,1] vaule = k*u[1]*(u[1]-1)
+    # Since a single entitiy cannot interact with it self!
+    k = 0.0
+    prev_reactant = -1
     @fastmath @inbounds for i in eachindex(r.reactants)
-      val *= u[r.reactants[i]]
+      if prev_reactant == r.reactants[i]
+        k += 1.0
+      else
+        k = 0.0
+      end
+      val *= (u[r.reactants[i]] - k)
+      prev_reactant = r.reactants[i]
     end
     val
   end

test/higher_order_reactions.jl

@@ -0,0 +1,40 @@
+using DiffEqBiological, DiffEqJump, DiffEqBase, OrdinaryDiffEq
+
+k = rand()
+# First order
+r1 = Reaction(k,[1],((1,-1),(2,1)))
+jump = build_jumps_from_reaction(r1)
+@test  jump.rate(0.0,[1,0]) == k
+
+
+# Second order
+r2 = Reaction(k,[1,1],((1,-2),(2,1)))
+jump2 = build_jumps_from_reaction(r2)
+@test  jump2.rate(0.0,[1,0]) == 0
+@test  jump2.rate(0.0,[2,0]) == k*2*1
+
+# Thrid order
+r3 = Reaction(k,[1,1,1],((1,-3),(2,1)))
+jump3 = build_jumps_from_reaction(r3)
+@test  jump3.rate(0.0,[1,0]) == 0
+@test  jump3.rate(0.0,[2,0]) == 0
+@test  jump3.rate(0.0,[3,0]) == k*3*2*1
+
+# Mixed Third order
+r21 = Reaction(k ,[1,1,2],((1,-2),(2,-1)))
+jump21 = build_jumps_from_reaction(r21)
+@test  jump21.rate(0.0,[1,1]) == 0
+@test  jump21.rate(0.0,[2,1]) == k*2*1*1
+
+# Mixed Third order different order
+r12 = Reaction(k ,[2,1,1],((1,-2),(2,-1)))
+jump12 = build_jumps_from_reaction(r12)
+  jump12.rate(0.0,[1,1]) == 0
+@test  jump12.rate(0.0,[2,1]) == k*2*1*1
+
+
+# Solve the dimerization problem 
+prob = DiscreteProblem([1,0],(0.0,250.0))
+jump_prob = GillespieProblem(prob,Direct(),r2)
+sol = solve(jump_prob,Discrete())
+@test find( x-> x!=0 ,[u[2] for u in sol.u]) == []

test/runtests.jl

@@ -5,4 +5,5 @@ tic()
 @time @testset "Gillespie Tests" begin include("gillespie.jl") end
 @time @testset "Variable Rate Reaction Tests" begin include("variable_rate_reactions.jl") end
 @time @testset "@reaction_network Tests" begin include("reaction_network.jl") end
+@time @testset "Higher order reaction Tests" begin include("higher_order_reactions.jl") end
 toc()