Skip to content

Commit

Permalink
init
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Jun 8, 2024
1 parent bc6d339 commit 00be503
Show file tree
Hide file tree
Showing 7 changed files with 19 additions and 20 deletions.
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# [Introduction to Catalyst and Julia for New Julia users](@id catalyst_for_new_julia_users)
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/).
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/).

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.

Expand Down Expand Up @@ -254,5 +254,5 @@ If you are a new Julia user who has used this tutorial, and there was something

---
## References
[^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)
[^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)
[^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)
[^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)
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# [Parameter Fitting for ODEs using Optimization.jl and DiffEqParamEstim.jl](@id optimization_parameter_fitting)
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.
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.

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/).

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
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](@id programmatic_generative_linear_pathway_intro)
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:
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:
```math
X_{i-1}/\tau,\hspace{0.33cm} ∅ \to X_{i}\\
1/\tau,\hspace{0.33cm} X_{i} \to ∅
Expand All @@ -21,7 +21,7 @@ for some kernel $g(\tau)$. Here, a common kernel is a [gamma distribution](https
```math
g(\tau; \alpha, \beta) = \frac{\beta^{\alpha}\tau^{\alpha-1}}{\Gamma(\alpha)}e^{-\beta\tau}
```
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.
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.

## [Modelling linear pathways using the DSL](@id programmatic_generative_linear_pathway_dsl)
It is known that two linear pathways have similar delays if the following equality holds:
Expand Down Expand Up @@ -132,6 +132,6 @@ plot!(sol_n20; idxs = :Xend, label = "n = 20")

---
## References
[^1]: [J. Metz, O. Diekmann *The Abstract Foundations of Linear Chain Trickery* (1991).](https://ir.cwi.nl/pub/1559/1559D.pdf)
[^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).
[^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)
[^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)
[^2]: [J. Metz, O. Diekmann *The Abstract Foundations of Linear Chain Trickery* (1991).](https://ir.cwi.nl/pub/1559/1559D.pdf)
[^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).
3 changes: 1 addition & 2 deletions docs/src/model_simulation/ode_simulation_performance.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
# [Advice for performant ODE simulations](@id ode_simulation_performance)
We have previously described how to perform ODE simulations of *chemical reaction network* (CRN) models. These simulations are typically fast and require little additional consideration. However, when a model is simulated many times (e.g. as a part of solving an [inverse problem](@ref ref)), or is very large, simulation run
times may become noticeable. Here we will give some advice on how to improve performance for these cases.
We have previously described how to perform ODE simulations of *chemical reaction network* (CRN) models. These simulations are typically fast and require little additional consideration. However, when a model is simulated many times (e.g. as a part of solving an [inverse problem](@ref ref)), or is very large, simulation run times may become noticeable. Here we will give some advice on how to improve performance for these cases [^1].

Generally, there are few good ways to, before a simulation, determine the best options. Hence, while we below provide several options, if you face an application for which reducing run time is critical (e.g. if you need to simulate the same ODE many times), it might be required to manually trial these various options to see which yields the best performance ([BenchmarkTools.jl's](https://github.com/JuliaCI/BenchmarkTools.jl) `@btime` macro is useful for this purpose). It should be noted that the default options typically perform well, and it is primarily for large models where investigating alternative options is worthwhile. All ODE simulations of Catalyst models are performed using the OrdinaryDiffEq.jl package, [which documentation](https://docs.sciml.ai/DiffEqDocs/stable/) provides additional advice on performance.

Expand Down
8 changes: 4 additions & 4 deletions docs/src/steady_state_functionality/dynamical_systems.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# [Analysing model steady state properties with DynamicalSystems.jl](@id dynamical_systems)
The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, however, these are not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems[^1][^2]. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, however, these are not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).

## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will typically be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors.
Expand Down Expand Up @@ -149,6 +149,6 @@ If you want to learn more about analysing dynamical systems, including chaotic b

---
## References
[^1]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7)
[^2]: [S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Westview Press (1994).](http://users.uoa.gr/~pjioannou/nonlin/Strogatz,%20S.%20H.%20-%20Nonlinear%20Dynamics%20And%20Chaos.pdf)
[^3]: [A. M. Lyapunov, *The general problem of the stability of motion*, International Journal of Control (1992).](https://www.tandfonline.com/doi/abs/10.1080/00207179208934253)
[^1]: [S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Westview Press (1994).](http://users.uoa.gr/~pjioannou/nonlin/Strogatz,%20S.%20H.%20-%20Nonlinear%20Dynamics%20And%20Chaos.pdf)
[^2]: [A. M. Lyapunov, *The general problem of the stability of motion*, International Journal of Control (1992).](https://www.tandfonline.com/doi/abs/10.1080/00207179208934253)
[^3]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7)
8 changes: 4 additions & 4 deletions docs/src/steady_state_functionality/homotopy_continuation.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@ method guaranteed to find all steady states for a system that has multiple ones.
However, many chemical reaction networks generate polynomial systems (for
example those which are purely mass action or have only have [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)) with
integer Hill exponents). The roots of these can reliably be found through a
*homotopy continuation* algorithm [^1]. This is implemented in Julia through the
[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/) package.
*homotopy continuation* algorithm[^1][^2]. This is implemented in Julia through the
[HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/) package[^3].

Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, `hc_steady_states`, that can be used to find the steady states of any `ReactionSystem` structure.

Expand Down Expand Up @@ -79,5 +79,5 @@ If you use this functionality in your research, please cite the following paper
---
## References
[^1]: [Andrew J Sommese, Charles W Wampler *The Numerical Solution of Systems of Polynomials Arising in Engineering and Science*, World Scientific (2005).](https://www.worldscientific.com/worldscibooks/10.1142/5763#t=aboutBook)
[^2]: [Paul Breiding, Sascha Timme, *HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia*, International Congress on Mathematical Software (2018).](https://link.springer.com/chapter/10.1007/978-3-319-96418-8_54)
[^43]: [Daniel J. Bates, Paul Breiding, Tianran Chen, Jonathan D. Hauenstein, Anton Leykin, Frank Sottile, *Numerical Nonlinear Algebra*, arXiv (2023).](https://arxiv.org/abs/2302.08585)
[^2]: [Daniel J. Bates, Paul Breiding, Tianran Chen, Jonathan D. Hauenstein, Anton Leykin, Frank Sottile, *Numerical Nonlinear Algebra*, arXiv (2023).](https://arxiv.org/abs/2302.08585)
[^3]: [Paul Breiding, Sascha Timme, *HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia*, International Congress on Mathematical Software (2018).](https://link.springer.com/chapter/10.1007/978-3-319-96418-8_54)
2 changes: 1 addition & 1 deletion docs/src/steady_state_functionality/nonlinear_solve.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# [Finding Steady States using NonlinearSolve.jl and SteadyStateDiffEq.jl](@id steady_state_solving)

Catalyst `ReactionSystem` models can be converted to ODEs (through [the reaction rate equation](@ref introduction_to_catalyst_ratelaws)). We have previously described how these ODEs' steady states can be found through [homotopy continuation](@ref homotopy_continuation). Generally, homotopy continuation (due to its ability to find *all* of a system's steady states) is the preferred approach. However, Catalyst supports two additional approaches for finding steady states:
- Through solving the nonlinear system produced by setting all ODE differentials to 0.
- Through solving the nonlinear system produced by setting all ODE differentials to 0[^1].
- Through forward ODE simulation from an initial condition until a steady state has been reached.

While these approaches only find a single steady state, they offer two advantages as compared to homotopy continuation:
Expand Down

0 comments on commit 00be503

Please sign in to comment.