DelaySSAToolkit

| Documentation | Build Status | |:-----------------:|:----------------:| | | |

DelaySSAToolkit.jl is a tool developed on top of JumpProcesses.jl in Julia which solves the stochastic simulation [1] coupled with delays. A portion of this library’s code is taken from the MIT licensed JumpProcesses.jl library.

For R implementation by Ziyan Jin, Xinyi Zhou, Zhaoyuan Fang, please check out https://github.com/Zoey-JIN/DelaySSA.

That code is copyright (c) 2017: Chris Rackauckas. This package contains the following features:

Features

• Various delay stochastic simulation algorithms are provided [2-6];
• Stochastic delay type is supported;
• Multiple delay channels and simultaneous delay reactions are supported;
• A cascade of delay reactions is supported (a delay reaction that causes other delay reactions);
• Priority queue and dependency graph are integrated for high computational performance;
• Ecosystem with Catalyst.jl, JumpProcesses.jl, DifferentialEquations.jl and more...

More information is available in the documentation. Please feel free to open issues and submit pull requests!

Installation

DelaySSAToolkit can be installed through the Julia package manager:

julia using Pkg Pkg.add("DelaySSAToolkit")

Examples

To run the following two examples, Catalyst.jl has to be installed by

julia using Pkg Pkg.add("Catalyst")

SEIR model

Check this example for more details.

```julia using Catalyst using DelaySSAToolkit

Model: Markovian part

rn = @reaction_network begin ρ, S + I --> E + I r, I --> R end u0 = [999, 1, 0, 0] # S, I, E, R tf = 400.0 tspan = (0, tf) ps = [1e-4, 1e-2]

Model: non-Markovian part (delay reactions)

delaytriggeraffect! = function (integrator, rng) τ = 20.0 append!(integrator.dechan[1], τ) # add a delay time τ to the first delay channel end delaytrigger = Dict(1 => delaytriggeraffect!) # the first reaction S+I -> E+I will trigger a delay reaction: E --> I after τ time.
delaycomplete = Dict(1 => [2 => 1, 3 => -1]) # E --> I after τ time: transfer from E (minus 1) to I (plus 1) after the completed delay reaction delayinterrupt = Dict() dechan0 = [[]] # initial condition for delay channel: no on-going delay reactions delayjumpset = DelayJumpSet(delaytrigger, delaycomplete, delayinterrupt)

convert the ReactionSystem to a JumpSystem

jumpsys = convert(JumpSystem, rn; combinatoricratelaws=false) dprob = DiscreteProblem(jumpsys, u0, tspan, ps) djprob = DelayJumpProblem( jumpsys, dprob, DelayRejection(), delayjumpset, dechan0; save_positions=(true, true) ) sol = solve(djprob, SSAStepper(); seed=1234)

] add Plots

using Plots; theme(:vibrant); plot( sol; label=["S" "I" "E" "R"], linewidth=3, legend=:top, ylabel="# of individuals", xlabel="Time", fmt=:png, ) ```

A bursty model [7]

Check this example for more details.

```julia using DelaySSAToolkit using Catalyst

Model: Markovian part

@parameters a b @variables t @species X(t) burstsup = 30 rxs = [Reaction(a * b^i / (1 + b)^(i + 1), nothing, [X], nothing, [i]) for i in 1:burstsup] rxs = vcat(rxs) @named rs = ReactionSystem(rxs, t, [X], [a, b]) u0 = [0] tf = 200.0 tspan = (0, tf) ps = [0.0282, 3.46]

Model: non-Markovian part

delaytriggeraffect! = [] for i in 1:burstsup push!(delaytriggeraffect!, function (integrator, rng) τ = 130.0 append!(integrator.dechan[1], fill(τ, i)) end) end delaytrigger = Dict([Pair(i, delaytriggeraffect![i]) for i in 1:burstsup]) delaycomplete = Dict(1 => [1 => -1]) delayinterrupt = Dict() dechan0 = [[]] delayjumpset = DelayJumpSet(delaytrigger, delaycomplete, delayinterrupt)

convert the ReactionSystem to a JumpSystem

jumpsys = convert(JumpSystem, rs; combinatoricratelaws=false) dprob = DiscreteProblem(jumpsys, u0, tspan, ps) djprob = DelayJumpProblem( jumpsys, dprob, DelayRejection(), delayjumpset, dechan0; save_positions=(false, false) ) ensprob = EnsembleProblem(djprob) ens = solve(ensprob, SSAStepper(), EnsembleThreads(); trajectories=10^5) ```

Recommendations

To solve a DelayJumpProblem, here are few recommendations for good performance:

• Use Catalyst.jl to build your Markovian model (model without delays). For certain algorithms that need dependency graph, it will be auto-generated.

• For a small number of jumps, DelayRejection and DelayDirect will often perform better than other aggregators.

• For large numbers of jumps with sparse chain like structures and similar jump rates, for example continuous time random walks, DelayDirectCR and DelayMNRM often have the best performance.

Other related packages

• FiniteStateProjection.jl: Finite State Projection algorithms for chemical reaction networks.
• MomentClosure.jl: Tools to generate and study moment equations for any chemical reaction network using various moment closure approximations.

Citation

@article{fuDelaySSAToolkitJlStochastic2022, title = {{{DelaySSAToolkit}}.Jl: Stochastic Simulation of Reaction Systems with Time Delays in {{Julia}}}, shorttitle = {{{DelaySSAToolkit}}.Jl}, author = {Fu, Xiaoming and Zhou, Xinyi and Gu, Dongyang and Cao, Zhixing and Grima, Ramon}, editor = {Valencia, Alfonso}, year = {2022}, month = sep, journal = {Bioinformatics}, volume = {38}, number = {17}, pages = {4243--4245}, issn = {1367-4803, 1460-2059}, doi = {10.1093/bioinformatics/btac472}, copyright = {All rights reserved}, langid = {english} }

References

[1] Daniel T. Gillespie, "Exact stochastic simulation of coupled chemical reactions", The Journal of Physical Chemistry 1977 81 (25), 2340-2361. https://doi.org/10.1021/j100540a008

[2] Barrio, Manuel, Kevin Burrage, André Leier, and Tianhai Tian. "Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation." PLoS computational biology 2, no. 9 (2006): e117. https://doi.org/10.1371/journal.pcbi.0020117

[3] Xiaodong Cai, "Exact stochastic simulation of coupled chemical reactions with delays", The Journal of Chemical Physics 126, 124108(2007). https://doi/10.1063/1.2710253

[4] David F. Anderson, "A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays", The Journal of Chemical Physics 128, 109903(2008). https://doi/10.1063/1.2799998

[5] Slepoy, Alexander, Aidan P. Thompson, and Steven J. Plimpton. "A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks." The journal of chemical physics 128, no. 20 (2008): 05B618. https://doi.org/10.1063/1.2919546

[6] Mauch, Sean, and Mark Stalzer. "Efficient formulations for exact stochastic simulation of chemical systems." IEEE/ACM Transactions on Computational Biology and Bioinformatics 8, no. 1 (2009): 27-35. https://doi.org/10.1109/TCBB.2009.47

[7] Qingchao Jiang, Xiaoming Fu, Shifu Yan, Runlai Li, Wenli Du, Zhixing Cao, Feng Qian, Ramon Grima, "Neural network aided approximation and parameter inference of non-Markovian models of gene expression". Nature communications, (2021) 12(1), 1-12. https://doi.org/10.1038/s41467-021-22919-1