You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
+3-3Lines changed: 3 additions & 3 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -1,5 +1,5 @@
1
1
# [Introduction to Catalyst and Julia for New Julia users](@id catalyst_for_new_julia_users)
2
-
The Catalyst tool for the modelling of chemical reaction networks is based in the Julia programming language. While experience in Julia programming is advantageous for using Catalyst, it is not necessary for accessing most of its basic features. This tutorial serves as an introduction to Catalyst for those unfamiliar with Julia, while also introducing some basic Julia concepts. Anyone who plans on using Catalyst extensively is recommended to familiarise oneself more thoroughly with the Julia programming language. A collection of resources for learning Julia can be found [here](https://julialang.org/learning/), and a full documentation is available [here](https://docs.julialang.org/en/v1/). A more practical (but also extensive) guide to Julia programming can be found [here](https://modernjuliaworkflows.github.io/writing/).
2
+
The Catalyst tool for the modelling of chemical reaction networks is based in the Julia programming language[^1][^2]. While experience in Julia programming is advantageous for using Catalyst, it is not necessary for accessing most of its basic features. This tutorial serves as an introduction to Catalyst for those unfamiliar with Julia, while also introducing some basic Julia concepts. Anyone who plans on using Catalyst extensively is recommended to familiarise oneself more thoroughly with the Julia programming language. A collection of resources for learning Julia can be found [here](https://julialang.org/learning/), and a full documentation is available [here](https://docs.julialang.org/en/v1/). A more practical (but also extensive) guide to Julia programming can be found [here](https://modernjuliaworkflows.github.io/writing/).
3
3
4
4
Julia can be downloaded [here](https://julialang.org/downloads/). Generally, it is recommended to use the [*juliaup*](https://github.com/JuliaLang/juliaup) tool to install and update Julia. Furthermore, *Visual Studio Code* is a good IDE with [extensive Julia support](https://code.visualstudio.com/docs/languages/julia), and a good default choice.
5
5
@@ -254,5 +254,5 @@ If you are a new Julia user who has used this tutorial, and there was something
254
254
255
255
---
256
256
## References
257
-
[^1]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)
258
-
[^2]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
257
+
[^1]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
258
+
[^2]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)
Copy file name to clipboardExpand all lines: docs/src/inverse_problems/optimization_ode_param_fitting.md
+1-1Lines changed: 1 addition & 1 deletion
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -1,5 +1,5 @@
1
1
# [Parameter Fitting for ODEs using Optimization.jl and DiffEqParamEstim.jl](@id optimization_parameter_fitting)
2
-
Fitting parameters to data involves solving an optimisation problem (that is, finding the parameter set that optimally fits your model to your data, typically by minimising a cost function). The SciML ecosystem's primary package for solving optimisation problems is [Optimization.jl](https://github.com/SciML/Optimization.jl). It provides access to a variety of solvers via a single common interface by wrapping a large number of optimisation libraries that have been implemented in Julia.
2
+
Fitting parameters to data involves solving an optimisation problem (that is, finding the parameter set that optimally fits your model to your data, typically by minimising a cost function)[^1]. The SciML ecosystem's primary package for solving optimisation problems is [Optimization.jl](https://github.com/SciML/Optimization.jl). It provides access to a variety of solvers via a single common interface by wrapping a large number of optimisation libraries that have been implemented in Julia.
3
3
4
4
This tutorial demonstrates both how to create parameter fitting cost functions using the [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl) package, and how to use Optimization.jl to minimise these. Optimization.jl can also be used in other contexts, such as [finding parameter sets that maximise the magnitude of some system behaviour](@ref behaviour_optimisation). More details on how to use these packages can be found in their [respective](https://docs.sciml.ai/Optimization/stable/)[documentations](https://docs.sciml.ai/DiffEqParamEstim/stable/).
Copy file name to clipboardExpand all lines: docs/src/model_creation/examples/programmatic_generative_linear_pathway.md
+5-5Lines changed: 5 additions & 5 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -2,7 +2,7 @@
2
2
This example will show how to use programmatic, generative, modelling to model a system implicitly. I.e. rather than listing all system reactions explicitly, the reactions are implicitly generated from a simple set of rules. This example is specifically designed to show how [programmatic modelling](@ref programmatic_CRN_construction) enables *generative workflows* (demonstrating one of its advantages as compared to [DSL-based modelling](@ref dsl_description)). In our example, we will model linear pathways, so we will first introduce these. Next, we will model them first using the DSL, and then using a generative programmatic workflow.
Linear pathways consists of a series of species ($X_0$, $X_1$, $X_2$, ..., $X_n$) where each activates the subsequent one. These are often modelled through the following reaction system:
5
+
Linear pathways consists of a series of species ($X_0$, $X_1$, $X_2$, ..., $X_n$) where each activates the subsequent one[^1]. These are often modelled through the following reaction system:
6
6
```math
7
7
X_{i-1}/\tau,\hspace{0.33cm} ∅ \to X_{i}\\
8
8
1/\tau,\hspace{0.33cm} X_{i} \to ∅
@@ -21,7 +21,7 @@ for some kernel $g(\tau)$. Here, a common kernel is a [gamma distribution](https
When this is converted to an ODE, this generates an integro-differential equation. These (as well as the simpler delay differential equations) can be difficult to solve and analyse (especially when SDE or jump simulations are desired). Here, *the linear chain trick* can be used to instead model the delay as a linear pathway of the form described above[^1]. A result by Fargue shows that this is equivalent to a gamma-distributed delay, where $\alpha$ is equivalent to $n$ (the number of species in our linear pathway) and $\beta$ to %\tau$ (the delay length term)[^2]. While modelling time delays using the linear chain trick introduces additional system species, it is often advantageous as it enables simulations using standard ODE, SDE, and Jump methods.
24
+
When this is converted to an ODE, this generates an integro-differential equation. These (as well as the simpler delay differential equations) can be difficult to solve and analyse (especially when SDE or jump simulations are desired). Here, *the linear chain trick* can be used to instead model the delay as a linear pathway of the form described above[^2]. A result by Fargue shows that this is equivalent to a gamma-distributed delay, where $\alpha$ is equivalent to $n$ (the number of species in our linear pathway) and $\beta$ to %\tau$ (the delay length term)[^3]. While modelling time delays using the linear chain trick introduces additional system species, it is often advantageous as it enables simulations using standard ODE, SDE, and Jump methods.
25
25
26
26
## [Modelling linear pathways using the DSL](@id programmatic_generative_linear_pathway_dsl)
27
27
It is known that two linear pathways have similar delays if the following equality holds:
[^1]: [J. Metz, O. Diekmann *The Abstract Foundations of Linear Chain Trickery* (1991).](https://ir.cwi.nl/pub/1559/1559D.pdf)
136
-
[^2]: D. Fargue *Reductibilite des systemes hereditaires a des systemes dynamiques (regis par des equations differentielles aux derivees partielles)*, Comptes rendus de l'Académie des Sciences (1973).
137
-
[^3]: [N. Korsbo, H. Jönsson *It’s about time: Analysing simplifying assumptions for modelling multi-step pathways in systems biology*, PLoS Computational Biology (2020).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007982)
135
+
[^1]: [N. Korsbo, H. Jönsson *It’s about time: Analysing simplifying assumptions for modelling multi-step pathways in systems biology*, PLoS Computational Biology (2020).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007982)
136
+
[^2]: [J. Metz, O. Diekmann *The Abstract Foundations of Linear Chain Trickery* (1991).](https://ir.cwi.nl/pub/1559/1559D.pdf)
137
+
[^3]: D. Fargue *Reductibilite des systemes hereditaires a des systemes dynamiques (regis par des equations differentielles aux derivees partielles)*, Comptes rendus de l'Académie des Sciences (1973).
Copy file name to clipboardExpand all lines: docs/src/model_creation/examples/smoluchowski_coagulation_equation.md
+66-50Lines changed: 66 additions & 50 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -4,72 +4,89 @@ This tutorial shows how to programmatically construct a [`ReactionSystem`](@ref)
4
4
The Smoluchowski coagulation equation describes a system of reactions in which monomers may collide to form dimers, monomers and dimers may collide to form trimers, and so on. This models a variety of chemical/physical processes, including polymerization and flocculation.
5
5
6
6
We begin by importing some necessary packages.
7
-
```julia
7
+
```@example smcoag1
8
8
using ModelingToolkit, Catalyst, LinearAlgebra
9
-
usingDiffEqBase, JumpProcesses
9
+
using JumpProcesses
10
10
using Plots, SpecialFunctions
11
11
```
12
12
Suppose the maximum cluster size is `N`. We assume an initial concentration of monomers, `Nₒ`, and let `uₒ` denote the initial number of monomers in the system. We have `nr` total reactions, and label by `V` the bulk volume of the system (which plays an important role in the calculation of rate laws since we have bimolecular reactions). Our basic parameters are then
13
-
```julia
14
-
## Parameter
15
-
N =10# maximum cluster size
16
-
Vₒ = (4π/3)*(10e-06*100)^3# volume of a monomers in cm³
The [Smoluchowski coagulation equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation) Wikipedia page illustrates the set of possible reactions that can occur. We can easily enumerate the `pair`s of multimer reactants that can combine when allowing a maximal cluster size of `N` monomers. We initialize the volumes of the reactant multimers as `volᵢ` and `volⱼ`
26
35
27
-
```julia
36
+
```@example smcoag1
28
37
# possible pairs of reactant multimers
29
38
pair = []
30
39
for i = 2:N
31
-
push!(pair, [1:integ(i/2) i .- (1:integ(i/2))])
40
+
halfi = floor(Int, i/2)
41
+
push!(pair, [(1:halfi) (i .- (1:halfi))])
32
42
end
33
43
pair = vcat(pair...)
34
-
vᵢ =@view pair[:,1] # Reactant 1 indices
35
-
vⱼ =@view pair[:,2] # Reactant 2 indices
36
-
volᵢ = Vₒ*vᵢ # cm⁻³
37
-
volⱼ = Vₒ*vⱼ # cm⁻³
44
+
vᵢ = @view pair[:, 1] # Reactant 1 indices
45
+
vⱼ = @view pair[:, 2] # Reactant 2 indices
46
+
volᵢ = Vₒ * vᵢ # cm⁻³
47
+
volⱼ = Vₒ * vⱼ # cm⁻³
38
48
sum_vᵢvⱼ = @. vᵢ + vⱼ # Product index
49
+
nothing #hide
39
50
```
40
51
We next specify the rates (i.e. kernel) at which reactants collide to form products. For simplicity, we allow a user-selected additive kernel or constant kernel. The constants(`B` and `C`) are adopted from Scott's paper [2](https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469_1968_025_0054_asocdc_2_0_co_2.xml)
41
-
```julia
52
+
```@example smcoag1
42
53
# set i to 1 for additive kernel, 2 for constant
43
54
i = 1
44
-
if i==1
45
-
B =1.53e03# s⁻¹
46
-
kv =@. B*(volᵢ + volⱼ)/V # dividing by volume as its a bi-molecular reaction chain
55
+
if i == 1
56
+
B = 1.53e03 # s⁻¹
57
+
58
+
# dividing by volume as it is a bimolecular reaction chain
59
+
kv = @. B * (volᵢ + volⱼ) / V
47
60
elseif i==2
48
-
C =1.84e-04# cm³ s⁻¹
49
-
kv =fill(C/V, nr)
61
+
C = 1.84e-04 # cm³ s⁻¹
62
+
kv = fill(C / V, nr)
50
63
end
64
+
nothing #hide
51
65
```
52
-
We'll store the reaction rates in `pars` as `Pair`s, and set the initial condition that only monomers are present at ``t=0`` in `u₀map`.
53
-
```julia
54
-
# unknown variables are X, pars stores rate parameters for each rx
66
+
We'll set the parameters and the initial condition that only monomers are present at ``t=0`` in `u₀map`.
67
+
```@example smcoag1
68
+
# k is a vector of the parameters, with values given by the vector kv
69
+
@parameters k[1:nr] = kv
70
+
71
+
# create the vector of species X_1,...,X_N
55
72
t = default_t()
56
-
@species k[1:nr] (X(t))[1:N]
57
-
pars =Pair.(collect(k), kv)
73
+
@species (X(t))[1:N]
58
74
59
75
# time-span
60
76
if i == 1
61
-
tspan = (0. ,2000.)
77
+
tspan = (0.0, 2000.0)
62
78
elseif i == 2
63
-
tspan = (0. ,350.)
79
+
tspan = (0.0, 350.0)
64
80
end
65
81
66
82
# initial condition of monomers
67
83
u₀ = zeros(Int64, N)
68
84
u₀[1] = uₒ
69
-
u₀map =Pair.(collect(X), u₀) # map variable to its initial value
85
+
u₀map = Pair.(collect(X), u₀) # map species to its initial value
86
+
nothing #hide
70
87
```
71
88
Here we generate the reactions programmatically. We systematically create Catalyst `Reaction`s for each possible reaction shown in the figure on [Wikipedia](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation). When `vᵢ[n] == vⱼ[n]`, we set the stoichiometric coefficient of the reactant multimer to two.
We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/jump_types/) documentation
0 commit comments