change format and compat info

Project.toml

@@ -22,13 +22,13 @@ DataStructures = "0.17, 0.18"
 DiffEqBase = "6.122"
 DocStringExtensions = "0.8, 0.9"
 FunctionWrappers = "1.0"
-JumpProcesses = "9.0"
-ModelingToolkit = "6.3, 7.0, 8.0"
+JumpProcesses = "9"
+ModelingToolkit = "6.3, 7.0, 8"
 RandomNumbers = "1.3"
 Reexport = "0.2, 1.0"
-SimpleUnPack = "1"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-julia = "1.6, 1.8"
+SimpleUnPack = "1"
+julia = "1.9"
 
 [extras]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"

README.md

@@ -4,9 +4,9 @@
 |:-----------------:|:----------------:|:----------------:|
 | [![doc dev badge]][doc dev link] | [![ci badge]][ci link] [![cov badge]][cov link] | [![download badge]][download link]| -->
 
-| **Documentation** | **Build Status** | **Code Style** |
-|:-----------------:|:----------------:|:----------------:|
-| [![doc dev badge]][doc dev link] | [![ci badge]][ci link] [![cov badge]][cov link] | [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
+| **Documentation** | **Build Status** |
+|:-----------------:|:----------------:|
+| [![doc dev badge]][doc dev link] | [![ci badge]][ci link] [![cov badge]][cov link] | 
 
 [doc dev badge]: https://img.shields.io/badge/docs-dev-blue.svg
 [doc dev link]: https://palmtree2013.github.io/DelaySSAToolkit.jl/dev/

docs/make.jl

@@ -1,39 +1,40 @@
 using DelaySSAToolkit
 using Documenter
 
-DocMeta.setdocmeta!(DelaySSAToolkit, :DocTestSetup, :(using DelaySSAToolkit);
-                    recursive = true)
+DocMeta.setdocmeta!(
+    DelaySSAToolkit, :DocTestSetup, :(using DelaySSAToolkit); recursive=true
+)
 
 makedocs(;
-         modules = [DelaySSAToolkit],
-         authors = "Xiaoming Fu",
-         repo = "https://github.com/palmtree2013/DelaySSAToolkit.jl/blob/{commit}{path}#{line}",
-         sitename = "DelaySSAToolkit.jl",
-         format = Documenter.HTML(;
-                                  mathengine = Documenter.Writers.HTMLWriter.MathJax2(),
-                                  prettyurls = get(ENV, "CI", "false") == "true",
-                                  canonical = "https://palmtree2013.github.io/DelaySSAToolkit.jl",
-                                  assets = String[]),
-         pages = [
-             "Home" => "index.md",
-             "Tutorials" => [
-                 "tutorials/tutorials.md",
-                 "tutorials/bursty.md",
-                 "tutorials/delay_degradation.md",
-                 "tutorials/heterogeneous_delay.md",
-                 "tutorials/delay_oscillator.md",
-                 "tutorials/stochastic_delay.md",
-             ],
-             "Algorithm" => [
-                 "algorithms/notations.md",
-                 "algorithms/delayrejection.md",
-                 "algorithms/delaydirect.md",
-                 "algorithms/delaymnrm.md",
-             ],
-             "Theory" => "theory.md",
-             "API" => "api.md",
-         ])
+    modules=[DelaySSAToolkit],
+    authors="Xiaoming Fu",
+    repo="https://github.com/palmtree2013/DelaySSAToolkit.jl/blob/{commit}{path}#{line}",
+    sitename="DelaySSAToolkit.jl",
+    format=Documenter.HTML(;
+        mathengine=Documenter.Writers.HTMLWriter.MathJax2(),
+        prettyurls=get(ENV, "CI", "false") == "true",
+        canonical="https://palmtree2013.github.io/DelaySSAToolkit.jl",
+        assets=String[],
+    ),
+    pages=[
+        "Home" => "index.md",
+        "Tutorials" => [
+            "tutorials/tutorials.md",
+            "tutorials/bursty.md",
+            "tutorials/delay_degradation.md",
+            "tutorials/heterogeneous_delay.md",
+            "tutorials/delay_oscillator.md",
+            "tutorials/stochastic_delay.md",
+        ],
+        "Algorithm" => [
+            "algorithms/notations.md",
+            "algorithms/delayrejection.md",
+            "algorithms/delaydirect.md",
+            "algorithms/delaymnrm.md",
+        ],
+        "Theory" => "theory.md",
+        "API" => "api.md",
+    ],
+)
 
-deploydocs(;
-           repo = "github.com/palmtree2013/DelaySSAToolkit.jl",
-           devbranch = "main")
+deploydocs(; repo="github.com/palmtree2013/DelaySSAToolkit.jl", devbranch="main")

docs/src/algorithms/delaydirect.md

@@ -1,42 +1,49 @@
-# Delay Direct Method Algorithm 
+# Delay Direct Method Algorithm
 
 The number of discarded $\Delta$ will be approximately equal to the number of delayed reactions that initiate. This follows because, other than the stored completions at the time the code terminates, every delayed completion will cause one computed $\Delta$ to be discarded. Thus, Cai [1] developed an algorithm, called the Direct Method for systems with delays, in which no random variables are discarded.
 
-The principle of Direct Method is the same as that of the original Gillespie Algorithm and the Rejection Method above: use one random variable to calculate when the next reaction initiates and use another random variable to calculate which reaction occurs at that future time. However, Direct Method updates the state of the system and propensity functions due to stored delayed reactions during the search for the next initiation time. In this way it is ensured that no random variables are discarded as in the Rejection Method. 
+The principle of Direct Method is the same as that of the original Gillespie Algorithm and the Rejection Method above: use one random variable to calculate when the next reaction initiates and use another random variable to calculate which reaction occurs at that future time. However, Direct Method updates the state of the system and propensity functions due to stored delayed reactions during the search for the next initiation time. In this way it is ensured that no random variables are discarded as in the Rejection Method.
 
 ## Algorithm
 
 Suppose that at time $t$ there are ongoing delayed reactions set to complete at times $t+T_1, t+T_2, \ldots, t+T_d$. Define $T_0=0$ and $T_{d+1}=\infty$.
 
 Define *Tstruct*, whose $i$th $(i=1,\dots,d)$ row stores $T_i$ and the index, $\mu_i$, of the reaction that $T_i$ is associated with.
-1. Initialize. Set the initial number of molecules of each species and set  $t=0$. Clear *Tstruct*.
-2. Calculate the propensity of function $a_k$, for each reaction $k \in 1,\ldots, M$.
-3. Set $a_0=\sum_{k=1}^M{a_k}$.
-4. Generate  $\Delta$.
-   - Input the time $t$ and $a_0=\sum_{k=1}^M{a_k}$.
-   - Generate an independent uniform $(0,1)$ random number $r_1$.
-   - If *Tstruct* is empty 
-     - This means there is no ongoing delay reactions, set $\Delta = 1/a_0\ln(1/r_1)$.
-   - Else
-     - Set $i=1$, $a_t = a_0T_1$ and  $F=1-e^{-a_t}$.
-     -  While $F < r_1$
-         - Update the state vector $\mathbf{x}$ due to the finish of the delayed reaction $t+T_i$.
-         - Calculate propensity $a_k(t+T_{i+1})$ due to the finish of the delayed reaction at $t+T_{i+1}$ and calculate $a_0(t+T_{i+1})$.
-         - Update $a_t=a_t+a_0(t+T_{i+1})(T_{i+1}-T_i)$.
-         - Update $F=1-e^{-a_t},i=i+1$.
-     - EndWhile
-     - Calculate propensity $a_k(t+T_i)$ due to the finish of the delayed reaction at $t+T_i$ and calculate $a_0(t+T_i)$.
-     - Set $\Delta=T_i-(\ln(1-r_1)+a_t-a_0(t+T_i)(T_{i+1}-T_i))/a_0(t+T_i)$.
-   - EndIf
-5. If $\Delta\in[T_i,T_{i+1})$, delete the columns $1,\ldots,i$ of $T_i$ and set $T_j=T_j-\Delta$.
-6. Generate an independent uniform $(0,1)$ random number $r_2$.
-7. Find $\mu\in[1,\dots,m]$ such that
-   ```math
-   \sum_{k=1}^{\mu-1} a_k < r_2 \leq \sum_{k=1}^{\mu}a_k
-   ```
-   where the $a_k$ and $a_0$ are generated in step 4.
-8. If $\mu\in \text{ND}$ , update the number of each molecular species according to the reaction $\mu$.
-9. If $\mu\in \text{CD}$, update *Tstruct* by adding the row $[\tau_\mu,\mu]$ so that $Tstruct(i,1)<Tstruct(i+1,1)$ still holds for all **i**.
+
+ 1. Initialize. Set the initial number of molecules of each species and set  $t=0$. Clear *Tstruct*.
+ 2. Calculate the propensity of function $a_k$, for each reaction $k \in 1,\ldots, M$.
+ 3. Set $a_0=\sum_{k=1}^M{a_k}$.
+ 4. Generate  $\Delta$.
+
+      + Input the time $t$ and $a_0=\sum_{k=1}^M{a_k}$.
+      + Generate an independent uniform $(0,1)$ random number $r_1$.
+      + If *Tstruct* is empty
+
+          * This means there is no ongoing delay reactions, set $\Delta = 1/a_0\ln(1/r_1)$.
+      + Else
+
+          * Set $i=1$, $a_t = a_0T_1$ and  $F=1-e^{-a_t}$.
+          * While $F < r_1$
+
+              - Update the state vector $\mathbf{x}$ due to the finish of the delayed reaction $t+T_i$.
+              - Calculate propensity $a_k(t+T_{i+1})$ due to the finish of the delayed reaction at $t+T_{i+1}$ and calculate $a_0(t+T_{i+1})$.
+              - Update $a_t=a_t+a_0(t+T_{i+1})(T_{i+1}-T_i)$.
+              - Update $F=1-e^{-a_t},i=i+1$.
+          * EndWhile
+          * Calculate propensity $a_k(t+T_i)$ due to the finish of the delayed reaction at $t+T_i$ and calculate $a_0(t+T_i)$.
+          * Set $\Delta=T_i-(\ln(1-r_1)+a_t-a_0(t+T_i)(T_{i+1}-T_i))/a_0(t+T_i)$.
+      + EndIf
+ 5. If $\Delta\in[T_i,T_{i+1})$, delete the columns $1,\ldots,i$ of $T_i$ and set $T_j=T_j-\Delta$.
+ 6. Generate an independent uniform $(0,1)$ random number $r_2$.
+ 7. Find $\mu\in[1,\dots,m]$ such that
+
+    ```math
+    \sum_{k=1}^{\mu-1} a_k < r_2 \leq \sum_{k=1}^{\mu}a_k
+    ```
+
+    where the $a_k$ and $a_0$ are generated in step 4.
+ 8. If $\mu\in \text{ND}$ , update the number of each molecular species according to the reaction $\mu$.
+ 9. If $\mu\in \text{CD}$, update *Tstruct* by adding the row $[\tau_\mu,\mu]$ so that $Tstruct(i,1)<Tstruct(i+1,1)$ still holds for all **i**.
 10. If $\mu\in \text{ICD}$, update the system according to the initiation of $\mu$ and update *Tstruct* by adding the row $[\tau_\mu,\mu]$ so that $Tstruct(i,1)<Tstruct(i+1,1)$ still holds for all $i$.
 11. Set $t=t+\Delta$.
 12. Return to Step 2 or quit.

docs/src/algorithms/delaymnrm.md

@@ -1,57 +1,73 @@
 # Delay Modified Next Reaction Method Algorithm
 
-## Representation using Poisson processes 
+## Representation using Poisson processes
 
 Let $v_k$, $v'_k\in Z^N_{\geq 0}$ be the vectors representing the number of each species consumed and created in the $k$th reaction, respectively. Then, if $N_k(t)$ is the number of initiations of reaction $k$ by time $t$, the state of the system at time $t$ is
+
 ```math
 X(t)=X(0)+\sum_{k=1}^M{N_k(t)(v^{'}_k-v_k)}.
 ```
+
 However, based upon the fundamental premise of stochastic chemical kinetics, $N_k(t)$ is a counting process with intensity $a_k(X(t))$ such that $\text{P}(N_k(t+\Delta t)-N_k(t)=1|X(s),s\leqslant t)=a_k(X(t))\Delta t$ for small $\Delta t$. Therefore, based upon the [counting process interpretation](https://en.wikipedia.org/wiki/Poisson_point_process), we have
+
 ```math
 N_k(t)=Y_k\Big(\int^t_0a_k(X(s))ds\Big),\tag{1} 
 ```
+
 where the $Y_k$ are independent, unit rate Poisson processes. Thus, $X(t)$ can be represented as the solution to the following equation:
+
 ```math
 X(t)=X(0)+\sum_{k=1}^M{Y_k\Big(\int^t_0a_k(X(s))ds\Big)(v'_k-v_k)}.
 ```
+
 All the randomness in the system is encapsulated in the $Y_k$'s and has therefore been separated from the state of the system. Thus, the system only changes when one of the $Y_k$'s changes. There are actually $M + 1$ relevant time frames in the system. The first time frame is the actual, or absolute time, $t$. However, each Poisson process $Y_k$ brings its own time frame. Thus, if we define $T_k(t)=\int^t_0a_k(X(s))ds$ for each $k$, then it is relevant for us to consider $Y_k(T_k(t))$. We will call $T_k(t)$ the "**internal time**" for reaction $k$.
 
 We denote by $P_k$ the first firing time of $Y_k$, in the time frame of $Y_k$, which is strictly larger than $T_k$. That is, $P_k=\min \left\{s>T_k:Y_k(s)>Y(T_k)\right\}$. The main idea of the following algorithm is that by Eq.(1) the value
+
 ```math
 \Delta t_k=\frac{P_k-T_k}{a_k}
 ```
+
 gives the amount of absolute time needed until the Poisson process $Y_k$ fires assuming that $a_k$ remains constant. $a_k$ does remain constant until the next reaction takes place. Therefore, a minimum of the different $\Delta t_k$ gives the time until the next reaction takes place.
 
 Now we extend the above notions to the delay system. No matter whether a reaction is contained in $\text{ND}$, $\text{CD}$, or $\text{ICD}$, the number of initiations at absolute time $t$ will be given by
+
 ```math
 \text{number of initiations of reaction } k \text{ by time } t = Y_k\Big(\int_{0}^{t} a_k(X(s)\, \mathrm{d}s\Big)
 ```
+
 where the $Y_k$ are independent, unit rate Poisson processes.
 
 Because the initiations are still given by the firing times of independent Poisson processes. Therefore, if $T_k$ is the current internal time of $Y_k$, $P_k$ the first internal time after $T_k$ at which $Y_k$ fires, and the propensity function for the $k$th reaction channel is given by $a_k$, then the time until the next initiation of reaction $k$(assuming no other reactions initiate or complete) is still given by $\Delta t_k= (P_k−T_k)/a_k$. The only change to the algorithm will be in keeping track and storing the delayed completions. To each delayed reaction channel we therefore assign a vector, $s_l$, that stores the completion times of that reaction in ascending order for $l=1,\ldots,L$, where $L$ is the total number of delay channels. Thus, the time until there is a change in the state of the system, be it an initiation or a completion, will be given by:
+
 ```math
 \Delta = \min\{\min_{k\in \{1,\ldots,M\}}\{\Delta t_k\}, \min_{l\in\{1,\ldots,L\}}\{s_l(1) − t\}\},
 ```
+
 where $t$ is the current time of the system. These ideas form the heart of the Next Reaction Method [1] for systems with delays.
 
 ## Algorithm
 
-1. Initialize. Set the initial number of molecules of each species and set $t = 0$. For each $k \leq M$, set $P_k = 0$ and $T_k = 0$, and for each delayed reaction channel set $s_l = [\infty]$.
-2. Calculate the propensity function, $a_k$, for each reaction.
-3. Generate $M$ independent, uniform $(0,1)$ random numbers, $r_k$, and set $P_k = \ln(1/r_k)$.
-4. Set $\Delta t_k = (P_k − T_k)/a_k$.
-5. Set $\Delta = \min\{\min_{k\in \{1,\ldots,M\}}\{\Delta t_k\}, \min_{l\in\{1,\ldots,L\}}\{s_l(1) − t\}\}$.
-6. Set $t = t + \Delta$.
-7. If we chose the completion of the delayed reaction $\mu$:
-   - Update the system based upon the completion of the reaction $\mu$.
-   - Delete the first row of $S_\mu$.
-8. Elseif reaction $\mu$ initiated and $\mu\in \text{ND}$
-   - Update the system according to reaction $\mu$.
-9. Elseif reaction $\mu$ initiated and $\mu\in \text{CD}$
-   - Update $s_\mu$ by inserting $t + \tau_\mu$ into $s_\mu$ in the second to last position.
+ 1. Initialize. Set the initial number of molecules of each species and set $t = 0$. For each $k \leq M$, set $P_k = 0$ and $T_k = 0$, and for each delayed reaction channel set $s_l = [\infty]$.
+ 2. Calculate the propensity function, $a_k$, for each reaction.
+ 3. Generate $M$ independent, uniform $(0,1)$ random numbers, $r_k$, and set $P_k = \ln(1/r_k)$.
+ 4. Set $\Delta t_k = (P_k − T_k)/a_k$.
+ 5. Set $\Delta = \min\{\min_{k\in \{1,\ldots,M\}}\{\Delta t_k\}, \min_{l\in\{1,\ldots,L\}}\{s_l(1) − t\}\}$.
+ 6. Set $t = t + \Delta$.
+ 7. If we chose the completion of the delayed reaction $\mu$:
+
+      + Update the system based upon the completion of the reaction $\mu$.
+      + Delete the first row of $S_\mu$.
+ 8. Elseif reaction $\mu$ initiated and $\mu\in \text{ND}$
+
+      + Update the system according to reaction $\mu$.
+ 9. Elseif reaction $\mu$ initiated and $\mu\in \text{CD}$
+
+      + Update $s_\mu$ by inserting $t + \tau_\mu$ into $s_\mu$ in the second to last position.
 10. Elseif reaction $\mu$ initiated and $\mu\in \text{ICD}$
-    - Update the system based upon the initiation of reaction $\mu$.
-    - Update $s_\mu$ by inserting $t + \tau_\mu$ into $s_\mu$ in the second to last position.
+
+      + Update the system based upon the initiation of reaction $\mu$.
+      + Update $s_\mu$ by inserting $t + \tau_\mu$ into $s_\mu$ in the second to last position.
 11. For each k, set $T_k = T_k + a_k \Delta$.
 12. If reaction $\mu$ initiated, let $r$ be uniform$(0,1)$ and set $P_\mu = P_\mu + \ln(1/r)$.
 13. Recalculate the propensity functions, $a_k$.

docs/src/algorithms/delayrejection.md

@@ -4,26 +4,30 @@ Simulation methods for systems with delays need to calculate when reactions init
 
 ## Algorithm
 
-1. Initialize. Set the initial number of molecules of each species and set $t = 0$.
-2. Calculate the propensity function, $a_k$, for each reaction.
-3. Set $a_0 = \sum_{k=1}^M a_k$.
-4. Generate an independent uniform $(0,1)$ random number, $r_1$, and set $\Delta = 1/a_0 \ln(1/r_1)$.
-5. If there is a delayed reaction set to finish in $[t, t + \Delta)$
-   - Discard $\Delta$.
-   - Update $t$ to be the time of the next delayed reaction,$\mu$.
-   - Update $\mathbf{x}$ according to the stored reaction $\mu$.
-6. Else
-   - Generate an independent uniform $(0,1)$ random number $r_2$.
-   - Find $\mu\in[1,\ldots, m]$ such that
-   ```math
-   \sum_{k=1}^{\mu-1} a_k(t) < r_2 a_0 \leq \sum_{k=1}^\mu a_k(t)
-   ```
-   - If $\mu\in \text{ND}$, update the number of each molecular species according to reaction $\mu$.
-   - If $\mu\in \text{CD}$, store the information that at time $t+\tau_\mu$ the system must be updated according to reaction $\mu$.
-   - If $\mu\in \text{ICD}$, update the system according to the initiation of $\mu$ and store that at time $t+\tau_\mu$ the system must be updated according to the completion of reaction $\mu$.
-   - Set $t = t +\Delta$
-7. Endif
-8. Return to step 2 or quit.
+ 1. Initialize. Set the initial number of molecules of each species and set $t = 0$.
+ 2. Calculate the propensity function, $a_k$, for each reaction.
+ 3. Set $a_0 = \sum_{k=1}^M a_k$.
+ 4. Generate an independent uniform $(0,1)$ random number, $r_1$, and set $\Delta = 1/a_0 \ln(1/r_1)$.
+ 5. If there is a delayed reaction set to finish in $[t, t + \Delta)$
+
+      + Discard $\Delta$.
+      + Update $t$ to be the time of the next delayed reaction,$\mu$.
+      + Update $\mathbf{x}$ according to the stored reaction $\mu$.
+ 6. Else
+
+      + Generate an independent uniform $(0,1)$ random number $r_2$.
+      + Find $\mu\in[1,\ldots, m]$ such that
+
+    ```math
+    \sum_{k=1}^{\mu-1} a_k(t) < r_2 a_0 \leq \sum_{k=1}^\mu a_k(t)
+    ```
+
+      + If $\mu\in \text{ND}$, update the number of each molecular species according to reaction $\mu$.
+      + If $\mu\in \text{CD}$, store the information that at time $t+\tau_\mu$ the system must be updated according to reaction $\mu$.
+      + If $\mu\in \text{ICD}$, update the system according to the initiation of $\mu$ and store that at time $t+\tau_\mu$ the system must be updated according to the completion of reaction $\mu$.
+      + Set $t = t +\Delta$
+ 7. Endif
+ 8. Return to step 2 or quit.
 
 ## References
 
@@ -32,4 +36,3 @@ Simulation methods for systems with delays need to calculate when reactions init
 
 [2]  Manuel Barrio, Kevin Burrage, André Leier, Tianhai Tian. "Oscillatory Regulation of Hes1: Discrete Stochastic Delay Modelling and Simulation", PLoS Computational Biology, 10.1371(2006).
 [https://doi.org/10.1371/journal.pcbi.0020117](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.0020117)
-

docs/src/algorithms/notations.md

@@ -23,9 +23,8 @@ Because the assumption above only pertains to the initiation times of reactions,
 
 **Case 3**: If reaction $k$ is in $\text{ICD}$ and initiates at time $t$, then the system is updated by losing the reactant species at the time of initiation, $t$, and is updated by gaining the product species at the time of completion, $t + \tau_k$.
 
-
-
 ## References
+
 [1] Dmitri A. Bratsun, Dmitri N. Volfson, Jeff Hasty, and Lev S. Tsimring "Non-Markovian processes in gene regulation (Keynote Address)", Proc. SPIE 5845, Noise in Complex Systems and Stochastic Dynamics III, (23 May 2005).
 [https://doi.org/10.1117/12.609707](https://doi.org/10.1117/12.609707)
 

docs/src/api.md

@@ -10,10 +10,11 @@ DelayJumpProblem
 ```
 
 ## Types and Algorithms
+
 ```@docs
 AbstractDelayAggregatorAlgorithm
 DelayDirect
 DelayRejection
 DelayMNRM
 DelayDirectCR
-```
+```