Merge pull request #21 from isaacsas/add-networks-from-mendes-paper

adding multistate network example

Author: ChrisRackauckas

src/DiffEqProblemLibrary.jl

@@ -49,6 +49,8 @@ export prob_dde_1delay, prob_dde_1delay_notinplace, prob_dde_1delay_scalar_notin
 # Jump Example Problems
 export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs,
 # examples mixing mass action and constant rate jumps
-       prob_jump_osc_mixed_jumptypes
+       prob_jump_osc_mixed_jumptypes,
+# examples used in published benchmarks / comparisions
+       prob_jump_multistate
 
 end # module

src/jump_premade_problems.jl

@@ -88,4 +88,56 @@ end
 u0 = [200.,60.,120.,100.,50.,50.,50.]  # Hill equations force use of floats!
 tf = 4000.
 prob = DiscreteProblem(u0, (0.,tf))
-prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)
+prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)
+
+
+"""
+    Multistate model from Gupta and Mendes, 
+    "An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems",
+    Computation 2018, 6, 9; doi:10.3390/computation6010009
+    Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml
+"""
+specs_sym_to_name = Dict(
+    :S1 => "R(a,l)",
+    :S2 => "L(r)",
+    :S3 => "A(Y~U,r)",
+    :S4 => "L(r!1).R(a,l!1)",
+    :S5 => "A(Y~U,r!1).R(a!1,l)",
+    :S6 => "A(Y~U,r!1).L(r!2).R(a!1,l!2)",
+    :S7 => "A(Y~P,r!1).L(r!2).R(a!1,l!2)",
+    :S8 => "A(Y~P,r!1).R(a!1,l)",
+    :S9 => "A(Y~P,r)")
+rates_sym_to_idx = Dict(
+    :R0 => 1,   :L0 => 2,    :A0 => 3,  :kon => 4, :koff => 5,
+    :kAon => 6, :kAoff => 7, :kAp => 8, :kAdp => 9)
+params = [5360, 1160, 5360, 0.01, 0.1, 0.01, 0.1, 0.01, 0.1]
+rs = @reaction_network msType begin
+    kon, S1 + S2 --> S4
+    kAon, S1 + S3 --> S5
+    kon, S2 + S5 --> S6
+    koff, S4 --> S1 + S2
+    kAon, S3 + S4 --> S6
+    kAoff, S5 --> S1 + S3
+    koff, S6 --> S2 + S5
+    kAoff, S6 --> S3 + S4
+    kAp, S6 --> S7
+    koff, S7 --> S2 + S8
+    kAoff, S7 --> S4 + S9
+    kAdp, S7 --> S6
+    kon, S2 + S8 --> S7
+    kAon, S1 + S9 --> S8
+    kAon, S4 + S9 --> S7
+    kAoff, S8 --> S1 + S9
+    kAdp, S8 --> S5
+    kAdp, S9 --> S3
+end kon kAon koff kAoff kAp kAdp
+rsi   = rates_sym_to_idx
+rates = params[[rsi[:kon], rsi[:kAon], rsi[:koff], rsi[:kAoff], rsi[:kAp], rsi[:kAdp]]]
+u0    = zeros(Int,9)
+u0[ findfirst(rs.syms, :S1) ] = params[1]
+u0[ findfirst(rs.syms, :S2) ] = params[2]
+u0[ findfirst(rs.syms, :S3) ] = params[3]
+tf    = 100.
+prob = DiscreteProblem(u0, (0., tf), rates)
+prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob, 
+                Dict("specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rates_sym_to_idx, "params" => params))