diff --git a/lib/JumpProblemLibrary/Project.toml b/lib/JumpProblemLibrary/Project.toml index 178dc05e..6e7ce779 100644 --- a/lib/JumpProblemLibrary/Project.toml +++ b/lib/JumpProblemLibrary/Project.toml @@ -1,18 +1,16 @@ name = "JumpProblemLibrary" uuid = "faf0f6d7-8cee-47cb-b27c-1eb80cef534e" -version = "2.0.3" +version = "3.0.0" [deps] -Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" -RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" [compat] -Catalyst = "16" Pkg = "1.10" -RuntimeGeneratedFunctions = "0.5" SafeTestsets = "0.1" -Test = "1.10" +SciMLBase = "3" SciMLTesting = "2.10" +Test = "1.10" julia = "1.10" [extras] diff --git a/lib/JumpProblemLibrary/docs/Project.toml b/lib/JumpProblemLibrary/docs/Project.toml index 6c9a7a3e..c8ce874b 100644 --- a/lib/JumpProblemLibrary/docs/Project.toml +++ b/lib/JumpProblemLibrary/docs/Project.toml @@ -4,5 +4,5 @@ JumpProblemLibrary = "faf0f6d7-8cee-47cb-b27c-1eb80cef534e" [compat] Documenter = "1" -JumpProblemLibrary = "2" +JumpProblemLibrary = "3" julia = "1.10" diff --git a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl index 9263e67e..c2f2b4f1 100644 --- a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl +++ b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl @@ -1,11 +1,6 @@ module JumpProblemLibrary -using Catalyst: Catalyst, @named, @parameters, @reaction_network, @species, Num, - Reaction, ReactionSystem, complete, default_t, hill, make_empty_network, - species, unknowns - -import RuntimeGeneratedFunctions -RuntimeGeneratedFunctions.init(@__MODULE__) +using SciMLBase: DiscreteProblem # Jump Example Problems export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs, @@ -19,103 +14,257 @@ export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs, """ JumpProblemNetwork -Container for the inputs needed to construct a `JumpProblem` from a -Catalyst `ReactionSystem`. Consumers materialize the actual `JumpProblem` -via the Catalyst API: +Container for the inputs needed to construct a `JumpProblem`. It does not store a +`JumpProblem` because that requires choosing an aggregator. + +Fields: + + - `rates`: rate constants of the mass action reactions, one per reaction. + - `tstop`: final simulation time. + - `u0`: initial species counts, ordered as `prob_data["species"]`. + - `discrete_prob`: `DiscreteProblem(u0, (0.0, tstop), rates)`. + - `prob_data`: a `Dict` holding the reaction network and any reference results. + +Every network stores these `prob_data` entries, with jumps ordered mass action +reactions first, then constant rate jumps: + + - `"species"`: species names, in state-vector order. + - `"reactant_stoich"`, `"net_stoich"`: stoichiometry of the mass action reactions as + vectors of `species index => coefficient` pairs (the `MassActionJump` format). + - `"constant_rate_jumps"`: vector of `(rate, net_stoich)` tuples for the reactions + whose propensity is not of mass action form, with `rate(u, p, t)` the propensity. + - `"jumps"`: `jumps(u, p, t)` returns the vector of all propensities, with + mass action propensities `k * ∏ binomial(u[s], c)`. + - `"nu"`: dense state change vector of each jump. + - `"dep_graph"`: `dep_graph[j]` lists the jumps whose propensity changes when jump + `j` fires, including `j`. + - `"vartojumps_map"`: `vartojumps_map[s]` lists the jumps whose propensity + depends on species `s`. + - `"jumptovars_map"`: `jumptovars_map[j]` lists the species changed by jump `j`. + +Example, using JumpProcesses: ```julia using JumpProcesses -JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) +jpn = prob_jump_dnarepressor +d = jpn.prob_data +majumps = MassActionJump(jpn.rates, d["reactant_stoich"], d["net_stoich"]) +jprob = JumpProblem(jpn.discrete_prob, NRM(), majumps; dep_graph = d["dep_graph"]) ``` + +`prob_jump_diffnetwork` instead stores functions of the lattice size, see its docstring. """ struct JumpProblemNetwork - network::Any # Catalyst ReactionSystem - rates::Any # vector of rate constants or nothing - tstop::Any # time to end simulation - u0::Any # initial values - prob_data::Any # additional problem data, stored as a Dict + rates::Any + tstop::Any + u0::Any + discrete_prob::Any + prob_data::Any +end + +""" +Convert reactions given as `reactants => products`, with each side a vector of +`species => coefficient` pairs, into `MassActionJump` reactant and net stoichiometry. +""" +function stoichiometry(species, reactions) + idx = Dict(s => i for (i, s) in enumerate(species)) + toidx(side) = sort!(Pair{Int, Int}[idx[s] => c for (s, c) in side]) + reactant_stoich = [toidx(r) for (r, _) in reactions] + net_stoich = map(reactions) do (r, p) + net = Dict{Int, Int}() + for (s, c) in r + net[idx[s]] = get(net, idx[s], 0) - c + end + for (s, c) in p + net[idx[s]] = get(net, idx[s], 0) + c + end + sort!(Pair{Int, Int}[s => c for (s, c) in net if c != 0]) + end + return reactant_stoich, net_stoich end -dna_rs = @reaction_network begin - k1, DNA --> mRNA + DNA - k2, mRNA --> mRNA + P - k3, mRNA --> 0 - k4, P --> 0 - k5, DNA + P --> DNAR - k6, DNAR --> DNA + P +function massaction_propensity(u, k, reactants) + a = float(k) + for (s, c) in reactants + a *= prod(u[s] - i for i in 0:(c - 1)) / factorial(c) + end + return a end -rates = [ - :k1 => 0.5, - :k2 => (20 * log(2.0) / 120.0), - :k3 => (log(2.0) / 120.0), - :k4 => (log(2.0) / 600.0), - :k5 => 0.025, - :k6 => 1.0, + +""" +Assemble the `prob_data` entries of a network. Each constant +rate jump is a `(rate, net_stoich, rate_species)` tuple, with `rate_species` the +indices of the species its propensity depends on. +""" +function network_data( + species, reactant_stoich, net_stoich; + constant_rate_jumps = Tuple{Function, Vector{Pair{Int, Int}}, Vector{Int}}[] + ) + nspecies = length(species) + crj_rates = [f for (f, _, _) in constant_rate_jumps] + all_net_stoich = vcat(net_stoich, [ns for (_, ns, _) in constant_rate_jumps]) + rate_species = vcat( + [first.(rs) for rs in reactant_stoich], [deps for (_, _, deps) in constant_rate_jumps] + ) + + function jumps(u, p, t) + props = [massaction_propensity(u, p[j], rs) for (j, rs) in enumerate(reactant_stoich)] + return vcat(props, [f(u, p, t) for f in crj_rates]) + end + + nu = map(all_net_stoich) do ns + v = zeros(Int, nspecies) + for (s, c) in ns + v[s] = c + end + v + end + + vartojumps_map = [Int[] for _ in 1:nspecies] + for (j, specs) in enumerate(rate_species), s in specs + push!(vartojumps_map[s], j) + end + jumptovars_map = [first.(ns) for ns in all_net_stoich] + dep_graph = [ + sort!(union([j], (vartojumps_map[s] for s in jumptovars_map[j])...)) + for j in eachindex(all_net_stoich) + ] + + return Dict{String, Any}( + "species" => species, + "reactant_stoich" => reactant_stoich, "net_stoich" => net_stoich, + "constant_rate_jumps" => [(f, ns) for (f, ns, _) in constant_rate_jumps], + "jumps" => jumps, "nu" => nu, "dep_graph" => dep_graph, + "vartojumps_map" => vartojumps_map, "jumptovars_map" => jumptovars_map + ) +end + +function build_network(rates, tstop, u0, species, reactions, extra_data = Dict(); kwargs...) + data = network_data(species, stoichiometry(species, reactions)...; kwargs...) + merge!(data, extra_data) + prob = DiscreteProblem(u0, (0.0, tstop), rates) + return JumpProblemNetwork(rates, tstop, u0, prob, data) +end + +const ∅ = Pair{Symbol, Int}[] + +dna_species = [:DNA, :mRNA, :P, :DNAR] +dna_reactions = [ + [:DNA => 1] => [:mRNA => 1, :DNA => 1], + [:mRNA => 1] => [:mRNA => 1, :P => 1], + [:mRNA => 1] => ∅, + [:P => 1] => ∅, + [:DNA => 1, :P => 1] => [:DNAR => 1], + [:DNAR => 1] => [:DNA => 1, :P => 1], ] -tf = 1000.0 -u0 = [:DNA => 1, :mRNA => 0, :P => 0, :DNAR => 0] -Nsims = 8000 -expected_avg = 5.92655375e+2 -prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) +rates = [0.5, (20 * log(2.0) / 120.0), (log(2.0) / 120.0), (log(2.0) / 600.0), 0.025, 1.0] """ DNA negative feedback autoregulatory model. Protein acts as repressor. + +Species `[DNA, mRNA, P, DNAR]`, reactions: + + 1. `k1, DNA --> mRNA + DNA` + 2. `k2, mRNA --> mRNA + P` + 3. `k3, mRNA --> 0` + 4. `k4, P --> 0` + 5. `k5, DNA + P --> DNAR` + 6. `k6, DNAR --> DNA + P` """ -prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob_data) +prob_jump_dnarepressor = build_network( + rates, 1000.0, [1, 0, 0, 0], dna_species, dna_reactions, + Dict("num_sims_for_mean" => 8000, "expected_mean" => 5.92655375e+2) +) -bd_rs = @reaction_network begin - k1, 0 --> A - k2, A --> 0 -end -rates = [:k1 => 1000.0, :k2 => 10.0] -tf = 1.0 -u0 = [:A => 0] -Nsims = 16000 -expected_avg = t -> rates[1] / rates[2] .* (1.0 - exp.(-rates[2] * t)) -prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg) +bd_rates = [1000.0, 10.0] +expected_avg = t -> bd_rates[1] / bd_rates[2] .* (1.0 .- exp.(-bd_rates[2] * t)) """ Simple birth-death process with constant production and degradation. + +Species `[A]`, reactions: + + 1. `k1, 0 --> A` + 2. `k2, A --> 0` """ -prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob_data) - -nonlin_rs = @reaction_network begin - k1, 2A --> B - k2, B --> 2A - k3, A + B --> C - k4, C --> A + B - k5, 3C --> 3A -end -rates = [:k1 => 1.0, :k2 => 2.0, :k3 => 0.5, :k4 => 0.75, :k5 => 0.25] -tf = 0.01 -u0 = [:A => 200, :B => 100, :C => 150] -Nsims = 32000 -expected_avg = 84.876015624999994 -prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) +prob_jump_constproduct = build_network( + bd_rates, 1.0, [0], [:A], [∅ => [:A => 1], [:A => 1] => ∅], + Dict("num_sims_for_mean" => 16000, "expected_mean_at_t" => expected_avg) +) + +nonlin_reactions = [ + [:A => 2] => [:B => 1], + [:B => 1] => [:A => 2], + [:A => 1, :B => 1] => [:C => 1], + [:C => 1] => [:A => 1, :B => 1], + [:C => 3] => [:A => 3], +] """ -Example with a mix of nonlinear reactions, including third order +Example with a mix of nonlinear reactions, including third order. + +Species `[A, B, C]`, reactions: + + 1. `k1, 2A --> B` + 2. `k2, B --> 2A` + 3. `k3, A + B --> C` + 4. `k4, C --> A + B` + 5. `k5, 3C --> 3A` """ -prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob_data) - -oscil_rs = @reaction_network begin - 0.01, (X, Y, Z) --> 0 - hill(X, 3.0, 100.0, -4), 0 --> Y - hill(Y, 3.0, 100.0, -4), 0 --> Z - hill(Z, 4.5, 100.0, -4), 0 --> X - hill(X, 2.0, 100.0, 6), 0 --> R - hill(Y, 15.0, 100.0, 4) * 0.002, R --> 0 - 20, 0 --> S - R * 0.005, S --> SP - 0.01, SP + SP --> SP2 - 0.05, SP2 --> 0 -end -u0 = [ - :X => 200.0, :Y => 60.0, :Z => 120.0, :R => 100.0, :S => 50.0, :SP => 50.0, - :SP2 => 50.0, -] # Hill equations force use of floats! -tf = 4000.0 +prob_jump_nonlinrxs = build_network( + [1.0, 2.0, 0.5, 0.75, 0.25], 0.01, [200, 100, 150], [:A, :B, :C], nonlin_reactions, + Dict("num_sims_for_mean" => 32000, "expected_mean" => 84.876015624999994) +) + +""" + hill(x, v, K, n) + +Hill function `v * x^n / (K^n + x^n)`, written in a form that is also finite at `x = 0`. +""" +hill(x, v, K, n) = v / (1 + (K / x)^n) + +oscil_species = [:X, :Y, :Z, :R, :S, :SP, :SP2] +oscil_reactions = [ + [:X => 1] => ∅, + [:Y => 1] => ∅, + [:Z => 1] => ∅, + ∅ => [:S => 1], + [:SP => 2] => [:SP2 => 1], + [:SP2 => 1] => ∅, +] +oscil_crjs = [ + ((u, p, t) -> hill(u[1], 3.0, 100.0, -4), [2 => 1], [1]), + ((u, p, t) -> hill(u[2], 3.0, 100.0, -4), [3 => 1], [2]), + ((u, p, t) -> hill(u[3], 4.5, 100.0, -4), [1 => 1], [3]), + ((u, p, t) -> hill(u[1], 2.0, 100.0, 6), [4 => 1], [1]), + ((u, p, t) -> hill(u[2], 15.0, 100.0, 4) * 0.002 * u[4], [4 => -1], [2, 4]), + ((u, p, t) -> 0.005 * u[4] * u[5], [5 => -1, 6 => 1], [4, 5]), +] """ Oscillatory system, uses a mixture of jump types. + +Species `[X, Y, Z, R, S, SP, SP2]`. Mass action reactions: + + 1. `0.01, X --> 0` + 2. `0.01, Y --> 0` + 3. `0.01, Z --> 0` + 4. `20, 0 --> S` + 5. `0.01, SP + SP --> SP2` + 6. `0.05, SP2 --> 0` + +Constant rate jumps, with `hill(x, v, K, n) = v * x^n / (K^n + x^n)`: + + 7. `hill(X, 3.0, 100.0, -4), 0 --> Y` + 8. `hill(Y, 3.0, 100.0, -4), 0 --> Z` + 9. `hill(Z, 4.5, 100.0, -4), 0 --> X` + 10. `hill(X, 2.0, 100.0, 6), 0 --> R` + 11. `hill(Y, 15.0, 100.0, 4) * 0.002, R --> 0` + 12. `R * 0.005, S --> SP` + +Species counts are floats since the Hill functions are evaluated on them. """ -prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, nothing) +prob_jump_osc_mixed_jumptypes = build_network( + [0.01, 0.01, 0.01, 20.0, 0.01, 0.05], 4000.0, + [200.0, 60.0, 120.0, 100.0, 50.0, 50.0, 50.0], oscil_species, oscil_reactions; + constant_rate_jumps = oscil_crjs +) specs_sym_to_name = Dict( :S1 => "R(a,l)", @@ -133,161 +282,169 @@ rsi = Dict( :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 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 -rates = [ - :kon, - :kAon, - :koff, - :kAoff, - :kAp, - :kAdp, -] .=> params[ - [ - rsi[:kon], rsi[:kAon], rsi[:koff], rsi[:kAoff], rsi[:kAp], rsi[:kAdp], - ], +multistate_reactions = [ + :kon => [:S1 => 1, :S2 => 1] => [:S4 => 1], + :kAon => [:S1 => 1, :S3 => 1] => [:S5 => 1], + :kon => [:S2 => 1, :S5 => 1] => [:S6 => 1], + :koff => [:S4 => 1] => [:S1 => 1, :S2 => 1], + :kAon => [:S3 => 1, :S4 => 1] => [:S6 => 1], + :kAoff => [:S5 => 1] => [:S1 => 1, :S3 => 1], + :koff => [:S6 => 1] => [:S2 => 1, :S5 => 1], + :kAoff => [:S6 => 1] => [:S3 => 1, :S4 => 1], + :kAp => [:S6 => 1] => [:S7 => 1], + :koff => [:S7 => 1] => [:S2 => 1, :S8 => 1], + :kAoff => [:S7 => 1] => [:S4 => 1, :S9 => 1], + :kAdp => [:S7 => 1] => [:S6 => 1], + :kon => [:S2 => 1, :S8 => 1] => [:S7 => 1], + :kAon => [:S1 => 1, :S9 => 1] => [:S8 => 1], + :kAon => [:S4 => 1, :S9 => 1] => [:S7 => 1], + :kAoff => [:S8 => 1] => [:S1 => 1, :S9 => 1], + :kAdp => [:S8 => 1] => [:S5 => 1], + :kAdp => [:S9 => 1] => [:S3 => 1], ] -u0 = [ - :S1 => params[1], :S2 => params[2], :S3 => params[3], :S4 => 0, :S5 => 0, - :S6 => 0, :S7 => 0, :S8 => 0, :S9 => 0, -] -tf = 100.0 """ 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 + +Species `[S1, ..., S9]`, reactions: + + 1. `kon, S1 + S2 --> S4` + 2. `kAon, S1 + S3 --> S5` + 3. `kon, S2 + S5 --> S6` + 4. `koff, S4 --> S1 + S2` + 5. `kAon, S3 + S4 --> S6` + 6. `kAoff, S5 --> S1 + S3` + 7. `koff, S6 --> S2 + S5` + 8. `kAoff, S6 --> S3 + S4` + 9. `kAp, S6 --> S7` + 10. `koff, S7 --> S2 + S8` + 11. `kAoff, S7 --> S4 + S9` + 12. `kAdp, S7 --> S6` + 13. `kon, S2 + S8 --> S7` + 14. `kAon, S1 + S9 --> S8` + 15. `kAon, S4 + S9 --> S7` + 16. `kAoff, S8 --> S1 + S9` + 17. `kAdp, S8 --> S5` + 18. `kAdp, S9 --> S3` + +`prob_data["params"]` holds the initial counts and rate constants, indexed by +`prob_data["rates_sym_to_idx"]`. """ -prob_jump_multistate = JumpProblemNetwork( - rs, rates, tf, u0, +prob_jump_multistate = build_network( + [params[rsi[k]] for (k, _) in multistate_reactions], 100.0, + [params[1], params[2], params[3], 0, 0, 0, 0, 0, 0], Symbol.(:S, 1:9), + last.(multistate_reactions), Dict( "specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rsi, "params" => params ) ) -# generate the network -N = 10 # number of genes -t = Catalyst.default_t() -@species (G(t))[1:(2N)] (M(t))[1:(2N)] (P(t))[1:(2N)] (G_ind(t))[1:(2N)] - -function construct_genenetwork(N) - rxs = Reaction[] +function twentygenes_reactions(N) + G(i) = Symbol(:G, i) + M(i) = Symbol(:M, i) + P(i) = Symbol(:P, i) + G_ind(i) = Symbol(:G_ind, i) + rates = Float64[] + reactions = Pair{Vector{Pair{Symbol, Int}}, Vector{Pair{Symbol, Int}}}[] for i in 1:N - push!(rxs, Reaction(10.0, [G[2 * i - i]], [G[2 * i - i], M[2 * i - i]])) - push!(rxs, Reaction(10.0, [M[2 * i - i]], [M[2 * i - i], P[2 * i - i]])) - push!(rxs, Reaction(1.0, [M[2 * i - i]], nothing)) - push!(rxs, Reaction(1.0, [P[2 * i - i]], nothing)) - # genenetwork *= "\t 10.0, G$(2*i-1) --> G$(2*i-1) + M$(2*i-1)\n" - # genenetwork *= "\t 10.0, M$(2*i-1) --> M$(2*i-1) + P$(2*i-1)\n" - # genenetwork *= "\t 1.0, M$(2*i-1) --> 0\n" - # genenetwork *= "\t 1.0, P$(2*i-1) --> 0\n" - - push!(rxs, Reaction(5.0, [G[2 * i]], [G[2 * i], M[2 * i]])) - push!(rxs, Reaction(5.0, [M[2 * i]], [M[2 * i], P[2 * i]])) - push!(rxs, Reaction(1.0, [M[2 * i]], nothing)) - push!(rxs, Reaction(1.0, [P[2 * i]], nothing)) - # genenetwork *= "\t 5.0, G$(2*i) --> G$(2*i) + M$(2*i)\n" - # genenetwork *= "\t 5.0, M$(2*i) --> M$(2*i) + P$(2*i)\n" - # genenetwork *= "\t 1.0, M$(2*i) --> 0\n" - # genenetwork *= "\t 1.0, P$(2*i) --> 0\n" - - push!(rxs, Reaction(0.0001, [G[2 * i], P[2 * i - i]], [G_ind[2 * i]])) - push!(rxs, Reaction(100.0, [G_ind[2 * i]], [G_ind[2 * i], M[2 * i]])) - # genenetwork *= "\t 0.0001, G$(2*i) + P$(2*i-1) --> G$(2*i)_ind \n" - # genenetwork *= "\t 100., G$(2*i)_ind --> G$(2*i)_ind + M$(2*i)\n" + for (g, k) in ((2i - 1, 10.0), (2i, 5.0)) + append!(rates, [k, k, 1.0, 1.0]) + push!(reactions, [G(g) => 1] => [G(g) => 1, M(g) => 1]) + push!(reactions, [M(g) => 1] => [M(g) => 1, P(g) => 1]) + push!(reactions, [M(g) => 1] => ∅) + push!(reactions, [P(g) => 1] => ∅) + end + append!(rates, [0.0001, 100.0]) + push!(reactions, [G(2i) => 1, P(2i - 1) => 1] => [G_ind(2i) => 1]) + push!(reactions, [G_ind(2i) => 1] => [G_ind(2i) => 1, M(2i) => 1]) end - spcs = reduce(vcat, collect.((G, M, P, G_ind))) - @named genenetwork = ReactionSystem(rxs, t, spcs, []) - return complete(genenetwork) + species = [f(i) for f in (G, M, P, G_ind) for i in 1:(2N)] + return rates, species, reactions end -rs = construct_genenetwork(N) -u0 = Num.(unknowns(rs)) .=> zeros(Int, length(unknowns(rs))) -for i in 1:(2 * N) - u0[findfirst(isequal(G[i]), unknowns(rs))] = (G[i] => 1) -end -tf = 2000.0 - +rates, species, reactions = twentygenes_reactions(10) """ Twenty-gene model from McCollum et al, "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior" Comp. Bio. and Chem., 30, pg. 39-49 (2006). + +Species `[G1, ..., G20, M1, ..., M20, P1, ..., P20, G_ind1, ..., G_ind20]`. For each +gene `g` with rate `k = 10` for odd and `k = 5` for even `g`: + + - `k, G_g --> G_g + M_g` + - `k, M_g --> M_g + P_g` + - `1, M_g --> 0` + - `1, P_g --> 0` + +and for each gene pair `i = 1, ..., 10`: + + - `0.0001, G_(2i) + P_(2i-1) --> G_ind_(2i)` + - `100, G_ind_(2i) --> G_ind_(2i) + M_(2i)` """ -prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, nothing) - -rn = @reaction_network begin - c1, G --> G + M - c2, M --> M + P - c3, M --> 0 - c4, P --> 0 - c5, 2P --> P2 - c6, P2 --> 2P - c7, P2 + G --> P2G - c8, P2G --> P2 + G -end -rnpar = [ - :c1 => 0.09, :c2 => 0.05, :c3 => 0.001, :c4 => 0.0009, :c5 => 0.00001, - :c6 => 0.0005, :c7 => 0.005, :c8 => 0.9, +prob_jump_twentygenes = build_network( + rates, 2000.0, [i <= 20 ? 1 : 0 for i in eachindex(species)], species, reactions +) + +dimer_species = [:G, :M, :P, :P2, :P2G] +dimer_reactions = [ + [:G => 1] => [:G => 1, :M => 1], + [:M => 1] => [:M => 1, :P => 1], + [:M => 1] => ∅, + [:P => 1] => ∅, + [:P => 2] => [:P2 => 1], + [:P2 => 1] => [:P => 2], + [:P2 => 1, :G => 1] => [:P2G => 1], + [:P2G => 1] => [:P2 => 1, :G => 1], ] -varlabels = ["G", "M", "P", "P2", "P2G"] -u0 = [:G => 1000, :M => 0, :P => 0, :P2 => 0, :P2G => 0] -tf = 4000.0 """ Negative feedback autoregulatory gene expression model. Dimer is the repressor. Taken from Marchetti, Priami and Thanh, "Simulation Algorithms for Computational Systems Biology", Springer (2017). + +Species `[G, M, P, P2, P2G]`, reactions: + + 1. `c1, G --> G + M` + 2. `c2, M --> M + P` + 3. `c3, M --> 0` + 4. `c4, P --> 0` + 5. `c5, 2P --> P2` + 6. `c6, P2 --> 2P` + 7. `c7, P2 + G --> P2G` + 8. `c8, P2G --> P2 + G` """ -prob_jump_dnadimer_repressor = JumpProblemNetwork( - rn, rnpar, tf, u0, - Dict("specs_names" => varlabels) +prob_jump_dnadimer_repressor = build_network( + [0.09, 0.05, 0.001, 0.0009, 0.00001, 0.0005, 0.005, 0.9], 4000.0, + [1000, 0, 0, 0, 0], dimer_species, dimer_reactions, + Dict("specs_names" => string.(dimer_species)) ) -# diffusion model -function getDiffNetwork(N) - diffnetwork = make_empty_network() - @parameters K - t = default_t() - @species (X(t))[1:N] - rxs = Reaction[] +function getDiffNetwork(N, K = 1.0) + species = Symbol.(:X, 1:N) + reactions = Pair{Vector{Pair{Symbol, Int}}, Vector{Pair{Symbol, Int}}}[] for i in 1:(N - 1) - push!(rxs, Reaction(K, [X[i]], [X[i + 1]])) - push!(rxs, Reaction(K, [X[i + 1]], [X[i]])) + push!(reactions, [species[i] => 1] => [species[i + 1] => 1]) + push!(reactions, [species[i + 1] => 1] => [species[i] => 1]) end - @named diffnetwork = ReactionSystem(rxs, t, collect(X), [K]) - return complete(diffnetwork) -end -params = [:K => 1.0] -function getDiffu0(diffnetwork, N) - return species(diffnetwork) .=> (10 * ones(Int64, N)) + rates = fill(float(K), length(reactions)) + return build_network(rates, 10.0, getDiffu0(N), species, reactions) end -tf = 10.0 +getDiffu0(N) = fill(10, N) """ -Continuous time random walk (i.e. diffusion approximation) example. -Here the network in the JumpProblemNetwork is a function that returns a -network given the number of lattice sites. -u0 is a similar function that returns the initial condition vector. +Continuous time random walk (i.e. diffusion approximation) example on a 1D lattice +of `N` sites with hopping rate `K` between neighboring sites: + + - `K, X_i --> X_(i+1)` + - `K, X_(i+1) --> X_i` + +for `i = 1, ..., N - 1`, starting with 10 particles per site. As the network depends +on `N`, `prob_data["network_func"](N, K)` returns the full `JumpProblemNetwork` for a +lattice of `N` sites, and `u0(N)` returns its initial condition. """ prob_jump_diffnetwork = JumpProblemNetwork( - getDiffNetwork, params, tf, getDiffu0, nothing + [1.0], 10.0, getDiffu0, nothing, Dict("network_func" => getDiffNetwork) ) end # module diff --git a/lib/JumpProblemLibrary/test/network_tests.jl b/lib/JumpProblemLibrary/test/network_tests.jl new file mode 100644 index 00000000..4a11a010 --- /dev/null +++ b/lib/JumpProblemLibrary/test/network_tests.jl @@ -0,0 +1,122 @@ +using JumpProblemLibrary, Test + +networks = [ + "dnarepressor" => prob_jump_dnarepressor, + "constproduct" => prob_jump_constproduct, + "nonlinrxs" => prob_jump_nonlinrxs, + "osc_mixed_jumptypes" => prob_jump_osc_mixed_jumptypes, + "multistate" => prob_jump_multistate, + "twentygenes" => prob_jump_twentygenes, + "dnadimer_repressor" => prob_jump_dnadimer_repressor, + "diffnetwork N = 6" => prob_jump_diffnetwork.prob_data["network_func"](6, 2.0), +] + +propensities(jpn, u) = jpn.prob_data["jumps"](u, jpn.rates, 0.0) + +@testset "$name: network data" for (name, jpn) in networks + d = jpn.prob_data + nspecies = length(d["species"]) + nma = length(d["reactant_stoich"]) + njumps = nma + length(d["constant_rate_jumps"]) + @test length(jpn.u0) == nspecies + @test length(jpn.rates) == nma == length(d["net_stoich"]) + @test jpn.discrete_prob.u0 == jpn.u0 + @test jpn.discrete_prob.tspan == (0.0, jpn.tstop) + @test length(d["nu"]) == length(d["dep_graph"]) == length(d["jumptovars_map"]) == njumps + @test length(d["vartojumps_map"]) == nspecies + + all_net_stoich = vcat(d["net_stoich"], last.(d["constant_rate_jumps"])) + for (nu, ns) in zip(d["nu"], all_net_stoich) + @test length(nu) == nspecies + @test nu == [sum((c for (s, c) in ns if s == i); init = 0) for i in 1:nspecies] + end + @test d["jumptovars_map"] == [findall(!iszero, nu) for nu in d["nu"]] + + # The dependency graphs are checked against the propensity functions themselves: + # perturbing species s changes only propensities in vartojumps_map[s], and firing + # jump j changes only those in dep_graph[j]. Taking the union over several states + # guards against a perturbation that leaves a propensity unchanged by coincidence. + states = [convert.(eltype(jpn.u0), 50 .+ mod.(37 .* (1:nspecies) .+ 11m^2, 101)) for m in 1:5] + for s in 1:nspecies + changed = map(states) do u + up = copy(u) + up[s] += 1 + findall(propensities(jpn, up) .!= propensities(jpn, u)) + end + @test all(c -> c ⊆ d["vartojumps_map"][s], changed) + @test sort!(union(changed...)) == d["vartojumps_map"][s] + end + for j in 1:njumps + changed = map(states) do u + findall(propensities(jpn, u .+ d["nu"][j]) .!= propensities(jpn, u)) + end + @test all(c -> c ⊆ d["dep_graph"][j], changed) + @test sort!(union([j], changed...)) == d["dep_graph"][j] + end + @test all(u -> length(propensities(jpn, u)) == njumps && all(>(0), propensities(jpn, u)), states) +end + +@testset "Propensities" begin + k = prob_jump_dnarepressor.rates + u = [1, 5, 7, 2] + @test propensities(prob_jump_dnarepressor, u) ≈ + [k[1] * 1, k[2] * 5, k[3] * 5, k[4] * 7, k[5] * 1 * 7, k[6] * 2] + + @test propensities(prob_jump_constproduct, [3]) ≈ [1000.0, 10.0 * 3] + @test prob_jump_constproduct.prob_data["expected_mean_at_t"](1.0) ≈ 100 * (1 - exp(-10)) + + k = prob_jump_nonlinrxs.rates + @test propensities(prob_jump_nonlinrxs, [200, 100, 150]) ≈ [ + k[1] * 200 * 199 / 2, k[2] * 100, k[3] * 200 * 100, k[4] * 150, + k[5] * 150 * 149 * 148 / 6, + ] + + hill(x, v, K, n) = v * x^n / (K^n + x^n) + X, Y, Z, R, S, SP, SP2 = prob_jump_osc_mixed_jumptypes.u0 + @test propensities(prob_jump_osc_mixed_jumptypes, prob_jump_osc_mixed_jumptypes.u0) ≈ [ + 0.01 * X, 0.01 * Y, 0.01 * Z, 20.0, 0.01 * SP * (SP - 1) / 2, 0.05 * SP2, + hill(X, 3.0, 100.0, -4), hill(Y, 3.0, 100.0, -4), hill(Z, 4.5, 100.0, -4), + hill(X, 2.0, 100.0, 6), hill(Y, 15.0, 100.0, 4) * 0.002 * R, 0.005 * R * S, + ] + @test propensities(prob_jump_osc_mixed_jumptypes, [0.0, 0.0, 0.0, R, S, SP, SP2])[7:10] == + [3.0, 3.0, 4.5, 0.0] + + k = prob_jump_dnadimer_repressor.rates + G, M, P, P2, P2G = 1000, 10, 20, 5, 3 + @test propensities(prob_jump_dnadimer_repressor, [G, M, P, P2, P2G]) ≈ [ + k[1] * G, k[2] * M, k[3] * M, k[4] * P, k[5] * P * (P - 1) / 2, k[6] * P2, + k[7] * P2 * G, k[8] * P2G, + ] + + kon, koff, kAon, kAoff, kAp, kAdp = 0.01, 0.1, 0.01, 0.1, 0.01, 0.1 + S = collect(1:9) .+ 10 + @test propensities(prob_jump_multistate, S) ≈ [ + kon * S[1] * S[2], kAon * S[1] * S[3], kon * S[2] * S[5], koff * S[4], + kAon * S[3] * S[4], kAoff * S[5], koff * S[6], kAoff * S[6], kAp * S[6], + koff * S[7], kAoff * S[7], kAdp * S[7], kon * S[2] * S[8], kAon * S[1] * S[9], + kAon * S[4] * S[9], kAoff * S[8], kAdp * S[8], kAdp * S[9], + ] + @test prob_jump_multistate.u0 == [5360, 1160, 5360, 0, 0, 0, 0, 0, 0] + + d = prob_jump_twentygenes.prob_data + @test length(d["species"]) == 80 + @test prob_jump_twentygenes.u0 == [fill(1, 20); fill(0, 60)] + sp = Dict(s => i for (i, s) in enumerate(d["species"])) + u = collect(1:80) + a = propensities(prob_jump_twentygenes, u) + @test length(a) == 100 + for i in 1:10, (g, kg) in ((2i - 1, 10.0), (2i, 5.0)) + @test kg * u[sp[Symbol(:G, g)]] in a + end + for i in 1:10 + cross = 0.0001 * u[sp[Symbol(:G, 2i)]] * u[sp[Symbol(:P, 2i - 1)]] + @test a[10i - 1] ≈ cross + @test d["net_stoich"][10i - 1] == + sort([sp[Symbol(:G, 2i)] => -1, sp[Symbol(:P, 2i - 1)] => -1, sp[Symbol(:G_ind, 2i)] => 1]) + end + + diff3 = prob_jump_diffnetwork.prob_data["network_func"](3, 2.0) + @test diff3.u0 == prob_jump_diffnetwork.u0(3) == [10, 10, 10] + @test propensities(diff3, [1, 2, 3]) == 2.0 .* [1, 2, 2, 3] + @test diff3.prob_data["nu"] == [[-1, 1, 0], [1, -1, 0], [0, -1, 1], [0, 1, -1]] +end diff --git a/lib/JumpProblemLibrary/test/qa/Project.toml b/lib/JumpProblemLibrary/test/qa/Project.toml index 1f3821d5..9b6c181a 100644 --- a/lib/JumpProblemLibrary/test/qa/Project.toml +++ b/lib/JumpProblemLibrary/test/qa/Project.toml @@ -6,7 +6,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" -JumpProblemLibrary = "2" +JumpProblemLibrary = "3" SciMLTesting = "2.4" Test = "1.10" julia = "1.10" diff --git a/lib/JumpProblemLibrary/test/runtests.jl b/lib/JumpProblemLibrary/test/runtests.jl index d84b237f..7e294d87 100644 --- a/lib/JumpProblemLibrary/test/runtests.jl +++ b/lib/JumpProblemLibrary/test/runtests.jl @@ -10,11 +10,11 @@ function activate_qa_env() return Pkg.instantiate() end -# The Core test is simply that all of the examples build (load the module). if TEST_GROUP == "Core" || TEST_GROUP == "All" @time @testset "Load Tests" begin @test JumpProblemLibrary isa Module end + @time @safetestset "Network data" include("network_tests.jl") end # Quality assurance: no undefined exports, stale dependencies, etc.