|
10 | 10 | # x_I(t)= 0<=t<=τ ? C*β/(a-γ)*((1-exp(-γ*t))/γ - (1-exp(-a*t))/a) : C*β/a*((1-exp(-γ*τ))/γ + exp(-a*t)*(1-exp(-(a-γ)τ))/(a-γ)) |
11 | 11 | # where a = β + γ |
12 | 12 |
|
13 | | - |
14 | 13 | using Catalyst |
15 | | -using DelaySSAToolkit |
| 14 | +using DelaySSAToolkit |
16 | 15 |
|
17 | 16 | rn = @reaction_network begin |
18 | | - C, 0 --> Xₐ |
19 | | - γ, Xₐ --> 0 |
20 | | - β, Xₐ --> Xᵢ |
21 | | - γ, Xᵢ --> 0 |
| 17 | + C, 0 --> Xₐ |
| 18 | + γ, Xₐ --> 0 |
| 19 | + β, Xₐ --> Xᵢ |
| 20 | + γ, Xᵢ --> 0 |
22 | 21 | end C γ β |
23 | 22 |
|
24 | 23 | jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false) |
25 | 24 |
|
26 | 25 | u0 = [0, 0] |
27 | 26 | tf = 30 |
28 | | -saveat = .1 |
| 27 | +saveat = 0.1 |
29 | 28 | de_chan0 = [[]] |
30 | | -C, γ, β = [2., 0.1, 0.5] |
| 29 | +C, γ, β = [2.0, 0.1, 0.5] |
31 | 30 | p = [C, γ, β] |
32 | | -tspan = (0.,tf) |
| 31 | +tspan = (0.0, tf) |
33 | 32 | # aggregatoralgo = DelayRejection() |
34 | 33 | # aggregatoralgo = DelayMNRM() |
35 | 34 | aggregatoralgo = DelayDirect() |
36 | 35 | # aggregatoralgo = DelayDirectCR() |
37 | 36 | dprob = DiscreteProblem(u0, tspan, p) |
38 | 37 |
|
39 | | -τ = 15. |
| 38 | +τ = 15.0 |
40 | 39 | delay_trigger_affect! = function (integrator, rng) |
41 | | - append!(integrator.de_chan[1], τ) |
| 40 | + append!(integrator.de_chan[1], τ) |
42 | 41 | end |
43 | | -delay_trigger = Dict(3=>delay_trigger_affect!) |
44 | | -delay_complete = Dict(1=>[2=>-1]) |
| 42 | +delay_trigger = Dict(3 => delay_trigger_affect!) |
| 43 | +delay_complete = Dict(1 => [2 => -1]) |
45 | 44 | delay_affect! = function (integrator, rng) |
46 | 45 | i = rand(rng, 1:length(integrator.de_chan[1])) |
47 | | - deleteat!(integrator.de_chan[1],i) |
| 46 | + deleteat!(integrator.de_chan[1], i) |
48 | 47 | end |
49 | | -delay_interrupt = Dict(4=>delay_affect!) |
50 | | -delaysets = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt) |
51 | | -djprob = DelayJumpProblem(jumpsys, dprob, aggregatoralgo, delaysets, de_chan0, save_positions = (false, false), save_delay_channel =false) |
52 | | -sol1 =@time solve(djprob, SSAStepper(), seed = 2, saveat =0.1) |
53 | | - |
54 | | - |
| 48 | +delay_interrupt = Dict(4 => delay_affect!) |
| 49 | +delaysets = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt) |
| 50 | +djprob = DelayJumpProblem(jumpsys, dprob, aggregatoralgo, delaysets, de_chan0, |
| 51 | + save_positions = (false, false), save_delay_channel = false) |
| 52 | +sol1 = @time solve(djprob, SSAStepper(), seed = 2, saveat = 0.1) |
55 | 53 |
|
56 | 54 | ens_prob = EnsembleProblem(djprob) |
57 | 55 | Sample_size = Int(10^4) |
58 | | -@time ens = solve(ens_prob, SSAStepper(),EnsembleThreads(),trajectories = Sample_size, saveat = .1) |
| 56 | +@time ens = solve(ens_prob, SSAStepper(), EnsembleThreads(), trajectories = Sample_size, |
| 57 | + saveat = 0.1) |
59 | 58 |
|
60 | | - |
61 | | -using Plots, DiffEqBase; theme(:vibrant) |
62 | | -plot(ens[1], label = ["X_A" "X_I"], fmt =:svg) |
| 59 | +using Plots, DiffEqBase; |
| 60 | +theme(:vibrant); |
| 61 | +plot(ens[1], label = ["X_A" "X_I"], fmt = :svg) |
63 | 62 | # savefig("docs/src/assets/delay_degradation1.svg") |
64 | 63 |
|
65 | | - |
66 | | -a = β + γ |
67 | | -x_A(t) = C/a*(1-exp(-a*t)) |
68 | | -x_I(t)= 0<=t<=τ ? C*β/(a-γ)*((1-exp(-γ*t))/γ - (1-exp(-a*t))/a) : C*β/a*((1-exp(-γ*τ))/γ + exp(-a*t)*(1-exp((a-γ)τ))/(a-γ)) |
69 | | - |
| 64 | +a = β + γ |
| 65 | +x_A(t) = C / a * (1 - exp(-a * t)) |
| 66 | +function x_I(t) |
| 67 | + 0 <= t <= τ ? C * β / (a - γ) * ((1 - exp(-γ * t)) / γ - (1 - exp(-a * t)) / a) : |
| 68 | + C * β / a * ((1 - exp(-γ * τ)) / γ + exp(-a * t) * (1 - exp((a - γ)τ)) / (a - γ)) |
| 69 | +end |
70 | 70 |
|
71 | 71 | using StatsBase |
72 | 72 | using Catalyst.EnsembleAnalysis |
73 | | -tsm(t) = Catalyst.EnsembleAnalysis.timepoint_mean(ens,t) |
| 73 | +tsm(t) = Catalyst.EnsembleAnalysis.timepoint_mean(ens, t) |
74 | 74 |
|
75 | 75 | mean_A(t) = tsm(t)[1] |
76 | 76 | mean_I(t) = tsm(t)[2] |
77 | 77 | timestamps = 0:0.1:tf |
78 | | -plot(timestamps,x_A.(timestamps),linewidth=3, label = "X_A Exact", xlabel = "Time", ylabel = "# of X_A") |
79 | | -plot!(timestamps,mean_A.(timestamps),linewidth=3,line=:dash, label = "X_A SSA") |
80 | | -plot!(timestamps,x_I.(timestamps),linewidth=3, label = "X_I Exact") |
81 | | -plot!(timestamps,mean_I.(timestamps),linewidth=3,line=:dash, legend = :topleft, label = "X_I SSA") |
| 78 | +plot(timestamps, x_A.(timestamps), linewidth = 3, label = "X_A Exact", xlabel = "Time", |
| 79 | + ylabel = "# of X_A") |
| 80 | +plot!(timestamps, mean_A.(timestamps), linewidth = 3, line = :dash, label = "X_A SSA") |
| 81 | +plot!(timestamps, x_I.(timestamps), linewidth = 3, label = "X_I Exact") |
| 82 | +plot!(timestamps, mean_I.(timestamps), linewidth = 3, line = :dash, legend = :topleft, |
| 83 | + label = "X_I SSA") |
82 | 84 | # savefig("docs/src/assets/delay_degradation2.svg") |
83 | | - |
84 | | - |
0 commit comments