JuliaDiffEq / DiffEqProblemLibrary.jl / 134

pending completion

Build # 134

Build Type

push

travis-ci

Committed by

ChrisRackauckas

Commit Message

<a href="https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/<a class=hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/<a class="double-link" href="https://git"><a class=hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/2a81b870f383498fdbf78992626c44c065357cae">2a81b870f">&lt;a href=&quot;https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/</a><a class="double-link" href="https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/2a81b870f383498fdbf78992626c44c065357cae">2a81b870f</a><a href="https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/2a81b870f383498fdbf78992626c44c065357cae">&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/JuliaDiffEq/DiffEqProblemLibrary.jl/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https:... (continued)

Run Details

20 of 20 new or added lines in 1 file covered. (100.0%)

102 of 619 relevant lines covered (16.48%)

7.66 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

60.0

/src/jump_premade_problems.jl

""" 
    General structure to hold JumpProblem info. Needed since
    the JumpProblem constructor requires the algorithm, so we 
    don't create the JumpProblem here.
"""
struct JumpProblemNetwork 
    network         # DiffEqBiological network
    rates           # vector of rate constants or nothing
    tstop           # time to end simulation
    u0              # initial values
    discrete_prob   # initialized discrete problem
    prob_data       # additional problem data, stored as a Dict
end

dna_rs = @reaction_network ptype begin
    k1, DNA --> mRNA + DNA
    k2, mRNA --> mRNA + P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
tf = 1000.0
u0 = [1,0,0,0]
prob = DiscreteProblem(u0, (0.0, tf), rates)
Nsims = 8000
expected_avg = 5.926553750000000e+02
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
"""
    DNA negative feedback autoregulatory model. Protein acts as repressor.    
"""
prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob, prob_data)

bd_rs = @reaction_network pdtype begin
    k1, 0 --> A
    k2, A --> 0
end k1 k2
rates = [1000., 10.]
tf = 1.0
u0 = [0]
prob = DiscreteProblem(u0, (0., tf), rates)
Nsims = 16000
expected_avg = t -> rates[1] / rates[2] .* ( 1. - exp.(-rates[2] * t))
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg)
"""
    Simple birth-death process with constant production and degradation.
"""
prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob, prob_data)

nonlin_rs = @reaction_network dtype begin
    k1, 2A --> B
    k2, B --> 2A 
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
rates = [1., 2., .5, .75, .25]
tf = .01
u0 = [200, 100, 150]
prob = DiscreteProblem(u0, (0., tf), rates)
Nsims = 32000
expected_avg = 84.876015624999994
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
"""
    Example with a mix of nonlinear reactions, including third order
"""
prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob, prob_data)


oscil_rs = @reaction_network rnoscType  begin
    0.01, (X,Y,Z) --> 0
    hill(X,3.,100.,-4), 0 --> Y
    hill(Y,3.,100.,-4), 0 --> Z
    hill(Z,4.5,100.,-4), 0 --> X
    hill(X,2.,100.,6), 0 --> R
    hill(Y,15.,100.,4)*0.002, R --> 0
    20, 0 --> S
    R*0.005, S --> SP
    0.01, SP + SP --> SP2
    0.05, SP2 --> 0
end
u0 = [200.,60.,120.,100.,50.,50.,50.]  # Hill equations force use of floats!
tf = 4000.
prob = DiscreteProblem(u0, (0.,tf))
"""
    Oscillatory system, uses a mixture of jump types.
"""
prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, prob, nothing)


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)
"""
    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
"""
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))


# generate the network
N = 10  # number of genes
genenetwork = "@reaction_network twentgtype begin\n"
for i in 1:N
    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"

    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"

    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"
end
genenetwork *= "end"
rs = eval( parse(genenetwork) )
u0 = zeros(Int, length(rs.syms))
for i = 1:(2*N)
    u0[findfirst(rs.syms, Symbol("G$(i)"))] = 1    
end
tf = 2000.0
prob = DiscreteProblem(u0, (0.0, tf))
"""
    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). 
"""
prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)


rn = @reaction_network gnrdtype 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 c1 c2 c3 c4 c5 c6 c7 c8
rnpar = [.09, .05, .001, .0009, .00001, .0005, .005, .9]
varlabels = ["G", "M", "P", "P2","P2G"]
u0 = [1000, 0, 0, 0,0]
tf = 4000.
prob = DiscreteProblem(u0, (0.0, tf), rnpar)
"""
    Negative feedback autoregulatory gene expression model. Dimer is the repressor.
    Taken from Marchetti, Priami and Thanh, 
    "Simulation Algorithms for Comptuational Systems Biology", 
    Springer (2017).
"""
prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob, 
                                                Dict("specs_names" => varlabels))


# diffusion model
function getDiffNetwork(N)
    diffnetwork = "@reaction_network dpldifftype begin\n"
    for i in 1:(N-1)
        diffnetwork *= "\t K, X$(i) --> X$(i+1)\n"
        diffnetwork *= "\t K, X$(i+1) --> X$(i)\n"
    end
    diffnetwork *= "end K"
    rs = eval( parse(diffnetwork) )
    rs
end
params = (1.,)
function getDiffu0(N)
    10*ones(Int64, N)
end
tf = 10.
"""
    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.
"""
prob_jump_diffnetwork = JumpProblemNetwork(getDiffNetwork, params, tf, getDiffu0, nothing, nothing)

1	"""
2	General structure to hold JumpProblem info. Needed since
3	the JumpProblem constructor requires the algorithm, so we
4	don't create the JumpProblem here.
5	"""
6	struct JumpProblemNetwork
7	network # DiffEqBiological network	8✔
8	rates # vector of rate constants or nothing
9	tstop # time to end simulation
10	u0 # initial values
11	discrete_prob # initialized discrete problem
12	prob_data # additional problem data, stored as a Dict
13	end
14
15	dna_rs = @reaction_network ptype begin
16	k1, DNA --> mRNA + DNA
17	k2, mRNA --> mRNA + P
18	k3, mRNA --> 0
19	k4, P --> 0
20	k5, DNA + P --> DNAR
21	k6, DNAR --> DNA + P
22	end k1 k2 k3 k4 k5 k6
23	rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
24	tf = 1000.0
25	u0 = [1,0,0,0]
26	prob = DiscreteProblem(u0, (0.0, tf), rates)
27	Nsims = 8000
28	expected_avg = 5.926553750000000e+02
29	prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
30	"""
31	DNA negative feedback autoregulatory model. Protein acts as repressor.
32	"""
33	prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob, prob_data)
34
35	bd_rs = @reaction_network pdtype begin
36	k1, 0 --> A
37	k2, A --> 0
38	end k1 k2
39	rates = [1000., 10.]
40	tf = 1.0
41	u0 = [0]
42	prob = DiscreteProblem(u0, (0., tf), rates)
43	Nsims = 16000
44	expected_avg = t -> rates[1] / rates[2] .* ( 1. - exp.(-rates[2] * t))
45	prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg)
46	"""
47	Simple birth-death process with constant production and degradation.
48	"""
49	prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob, prob_data)
50
51	nonlin_rs = @reaction_network dtype begin
52	k1, 2A --> B
53	k2, B --> 2A
54	k3, A + B --> C
55	k4, C --> A + B
56	k5, 3C --> 3A
57	end k1 k2 k3 k4 k5
58	rates = [1., 2., .5, .75, .25]
59	tf = .01
60	u0 = [200, 100, 150]
61	prob = DiscreteProblem(u0, (0., tf), rates)
62	Nsims = 32000
63	expected_avg = 84.876015624999994
64	prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
65	"""
66	Example with a mix of nonlinear reactions, including third order
67	"""
68	prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob, prob_data)
69
70
71	oscil_rs = @reaction_network rnoscType begin
72	0.01, (X,Y,Z) --> 0
73	hill(X,3.,100.,-4), 0 --> Y
74	hill(Y,3.,100.,-4), 0 --> Z
75	hill(Z,4.5,100.,-4), 0 --> X
76	hill(X,2.,100.,6), 0 --> R
77	hill(Y,15.,100.,4)*0.002, R --> 0
78	20, 0 --> S
79	R*0.005, S --> SP
80	0.01, SP + SP --> SP2
81	0.05, SP2 --> 0
82	end
83	u0 = [200.,60.,120.,100.,50.,50.,50.] # Hill equations force use of floats!
84	tf = 4000.
85	prob = DiscreteProblem(u0, (0.,tf))
86	"""
87	Oscillatory system, uses a mixture of jump types.
88	"""
89	prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, prob, nothing)
90
91
92	specs_sym_to_name = Dict(
93	:S1 => "R(a,l)",
94	:S2 => "L(r)",
95	:S3 => "A(Y~U,r)",
96	:S4 => "L(r!1).R(a,l!1)",
97	:S5 => "A(Y~U,r!1).R(a!1,l)",
98	:S6 => "A(Y~U,r!1).L(r!2).R(a!1,l!2)",
99	:S7 => "A(Y~P,r!1).L(r!2).R(a!1,l!2)",
100	:S8 => "A(Y~P,r!1).R(a!1,l)",
101	:S9 => "A(Y~P,r)")
102	rates_sym_to_idx = Dict(
103	:R0 => 1, :L0 => 2, :A0 => 3, :kon => 4, :koff => 5,
104	:kAon => 6, :kAoff => 7, :kAp => 8, :kAdp => 9)
105	params = [5360, 1160, 5360, 0.01, 0.1, 0.01, 0.1, 0.01, 0.1]
106	rs = @reaction_network msType begin
107	kon, S1 + S2 --> S4
108	kAon, S1 + S3 --> S5
109	kon, S2 + S5 --> S6
110	koff, S4 --> S1 + S2
111	kAon, S3 + S4 --> S6
112	kAoff, S5 --> S1 + S3
113	koff, S6 --> S2 + S5
114	kAoff, S6 --> S3 + S4
115	kAp, S6 --> S7
116	koff, S7 --> S2 + S8
117	kAoff, S7 --> S4 + S9
118	kAdp, S7 --> S6
119	kon, S2 + S8 --> S7
120	kAon, S1 + S9 --> S8
121	kAon, S4 + S9 --> S7
122	kAoff, S8 --> S1 + S9
123	kAdp, S8 --> S5
124	kAdp, S9 --> S3
125	end kon kAon koff kAoff kAp kAdp
126	rsi = rates_sym_to_idx
127	rates = params[[rsi[:kon], rsi[:kAon], rsi[:koff], rsi[:kAoff], rsi[:kAp], rsi[:kAdp]]]
128	u0 = zeros(Int,9)
129	u0[ findfirst(rs.syms, :S1) ] = params[1]
130	u0[ findfirst(rs.syms, :S2) ] = params[2]
131	u0[ findfirst(rs.syms, :S3) ] = params[3]
132	tf = 100.
133	prob = DiscreteProblem(u0, (0., tf), rates)
134	"""
135	Multistate model from Gupta and Mendes,
136	"An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems",
137	Computation 2018, 6, 9; doi:10.3390/computation6010009
138	Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml
139	"""
140	prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob,
141	Dict("specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rates_sym_to_idx, "params" => params))
142
143
144	# generate the network
145	N = 10 # number of genes
146	genenetwork = "@reaction_network twentgtype begin\n"
147	for i in 1:N
148	genenetwork = "\t 10.0, G$(2i-1) --> G$(2i-1) + M$(2i-1)\n"	10✔
149	genenetwork = "\t 10.0, M$(2i-1) --> M$(2i-1) + P$(2i-1)\n"	10✔
150	genenetwork = "\t 1.0, M$(2i-1) --> 0\n"	10✔
151	genenetwork = "\t 1.0, P$(2i-1) --> 0\n"	10✔
152
153	genenetwork = "\t 5.0, G$(2i) --> G$(2i) + M$(2i)\n"	10✔
154	genenetwork = "\t 5.0, M$(2i) --> M$(2i) + P$(2i)\n"	10✔
155	genenetwork = "\t 1.0, M$(2i) --> 0\n"	10✔
156	genenetwork = "\t 1.0, P$(2i) --> 0\n"	10✔
157
158	genenetwork = "\t 0.0001, G$(2i) + P$(2i-1) --> G$(2i)_ind \n"	10✔
159	genenetwork = "\t 100., G$(2i)_ind --> G$(2i)_ind + M$(2i)\n"	10✔
160	end
161	genenetwork *= "end"
162	rs = eval( parse(genenetwork) )
163	u0 = zeros(Int, length(rs.syms))
164	for i = 1:(2*N)
165	u0[findfirst(rs.syms, Symbol("G$(i)"))] = 1	20✔
166	end
167	tf = 2000.0
168	prob = DiscreteProblem(u0, (0.0, tf))
169	"""
170	Twenty-gene model from McCollum et al,
171	"The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior"
172	Comp. Bio. and Chem., 30, pg. 39-49 (2006).
173	"""
174	prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)
175
176
177	rn = @reaction_network gnrdtype begin
178	c1, G --> G + M
179	c2, M --> M + P
180	c3, M --> 0
181	c4, P --> 0
182	c5, 2P --> P2
183	c6, P2 --> 2P
184	c7, P2 + G --> P2G
185	c8, P2G --> P2 + G
186	end c1 c2 c3 c4 c5 c6 c7 c8
187	rnpar = [.09, .05, .001, .0009, .00001, .0005, .005, .9]
188	varlabels = ["G", "M", "P", "P2","P2G"]
189	u0 = [1000, 0, 0, 0,0]
190	tf = 4000.
191	prob = DiscreteProblem(u0, (0.0, tf), rnpar)
192	"""
193	Negative feedback autoregulatory gene expression model. Dimer is the repressor.
194	Taken from Marchetti, Priami and Thanh,
195	"Simulation Algorithms for Comptuational Systems Biology",
196	Springer (2017).
197	"""
198	prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob,
199	Dict("specs_names" => varlabels))
200
201
202	# diffusion model
203	function getDiffNetwork(N)
204	diffnetwork = "@reaction_network dpldifftype begin\n"	×
205	for i in 1:(N-1)	×
206	diffnetwork *= "\t K, X$(i) --> X$(i+1)\n"	×
207	diffnetwork *= "\t K, X$(i+1) --> X$(i)\n"	×
208	end
209	diffnetwork *= "end K"	×
210	rs = eval( parse(diffnetwork) )	×
211	rs	×
212	end
213	params = (1.,)
214	function getDiffu0(N)
215	10*ones(Int64, N)	×
216	end
217	tf = 10.
218	"""
219	Continuous time random walk (i.e. diffusion approximation) example.
220	Here the network in the JumpProblemNetwork is a function that returns a
221	network given the number of lattice sites.
222	u0 is a similar function that returns the initial condition vector.
223	"""
224	prob_jump_diffnetwork = JumpProblemNetwork(getDiffNetwork, params, tf, getDiffu0, nothing, nothing)

JuliaDiffEq / DiffEqProblemLibrary.jl / 134

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous