diff --git a/docs/Project.toml b/docs/Project.toml index 5abb2891e3..35c1fa384c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" @@ -10,8 +11,10 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" GlobalSensitivity = "af5da776-676b-467e-8baf-acd8249e4f0f"# HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" +IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" @@ -43,12 +46,13 @@ Distributions = "0.25" Documenter = "1.4.1" GlobalSensitivity = "2.4.0" HomotopyContinuation = "2.6" +JumpProcesses = "9" Latexify = "0.15, 0.16" ModelingToolkit = "9.5" NonlinearSolve = "3.4.0" Optim = "1" Optimization = "3.19" -OptimizationBBO = "0.1.5" +OptimizationBBO = "0.1.5, 0.2" OptimizationNLopt = "0.1.8" OptimizationOptimJL = "0.1.14" OptimizationOptimisers = "0.1.1" diff --git a/docs/pages.jl b/docs/pages.jl index 1c42d369ba..96e6587bf3 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -6,8 +6,8 @@ pages = Any[ # Advanced introduction. ], "Model Creation and Properties" => Any[ - "model_creation/dsl_description.md", - # DSL - Advanced + "model_creation/dsl_basics.md", + "model_creation/dsl_advanced.md", "model_creation/programmatic_CRN_construction.md", "model_creation/compositional_modeling.md", "model_creation/constraint_equations.md", @@ -18,18 +18,20 @@ pages = Any[ "model_creation/network_analysis.md", "model_creation/chemistry_related_functionality.md", "Model creation examples" => Any[ - "model_creation/examples/basic_CRN_examples.md", + "model_creation/examples/basic_CRN_library.md", "model_creation/examples/programmatic_generative_linear_pathway.md", "model_creation/examples/hodgkin_huxley_equation.md", "model_creation/examples/smoluchowski_coagulation_equation.md" ] ], "Model simulation" => Any[ + "model_simulation/simulation_introduction.md", # Simulation introduction. - # Simulation plotting. + "model_simulation/simulation_plotting.md", "model_simulation/simulation_structure_interfacing.md", "model_simulation/ensemble_simulations.md", # Stochastic simulation statistical analysis. + "model_simulation/ode_simulation_performance.md", # ODE Performance considerations/advice. # SDE Performance considerations/advice. # Jump Performance considerations/advice. diff --git a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md index 5e23f8235f..82c297b33d 100644 --- a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md +++ b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md @@ -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 some of its basic features. This tutorial serves as an introduction to Catalyst for those unfamiliar with Julia, 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/). +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/). 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. @@ -19,12 +19,13 @@ area = length * width ```@example ex1 min(1.0, 3.0) ``` +Julia has a specific *help mode*, which can be [queried for information about any function](https://docs.julialang.org/en/v1/stdlib/REPL/#Help-mode) (including those defined by Catalyst). -Each Julia variable has a specific *type*, designating what type of value it is. While not directly required to use Catalyst, this is useful to be aware of. To learn the type of a specific variable, use the `typeof` function. More information about types can be [found here](https://docs.julialang.org/en/v1/manual/types/). +Each Julia variable has a specific *type*, designating what type of value it contains. While not directly required to use Catalyst, this is useful to be aware of. To learn the type of a specific variable, use the `typeof` function. More information about types can be [found here](https://docs.julialang.org/en/v1/manual/types/). ```@example ex1 typeof(1.0) ``` -Here, `Float64` denotes decimal-valued numbers. Integer-valued numbers instead are of the `Int64` type. +Here, `Float64` denotes decimal-valued numbers. Integer-valued numbers instead have the `Int64` type. ```@example ex1 typeof(1) ``` @@ -65,7 +66,7 @@ using Catalyst using DifferentialEquations using Plots ``` -Here, if we restart Julia, these commands *do need to be rerun. +Here, if we restart Julia, these `using` commands *must be rerun*. A more comprehensive (but still short) introduction to package management in Julia is provided at [the end of this documentation page](@ref catalyst_for_new_julia_users_packages). It contains some useful information and is hence highly recommended reading. For a more detailed introduction to Julia package management, please read [the Pkg documentation](https://docs.julialang.org/en/v1/stdlib/Pkg/). @@ -79,18 +80,18 @@ The `@reaction_network` command is followed by the `begin` keyword, which is fol Here, we create a simple *birth-death* model, where a single species ($X$) is created at rate $b$, and degraded at rate $d$. The model is stored in the variable `rn`. ```@example ex2 rn = @reaction_network begin - b, 0 --> X - d, X --> 0 + b, 0 --> X + d, X --> 0 end ``` -For more information on how to use the Catalyst model creator (also known as the Catalyst DSL), please read [the corresponding documentation](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/). +For more information on how to use the Catalyst model creator (also known as *the Catalyst DSL*), please read [the corresponding documentation](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/). Next, we wish to simulate our model. To do this, we need to provide some additional information to the simulator. This is * The initial condition. That is, the concentration (or copy numbers) of each species at the start of the simulation. -* The timespan. That is, the timeframe over which we wish to run the simulation. +* The time span. That is, the time frame over which we wish to run the simulation. * The parameter values. That is, the values of the model's parameters for this simulation. -The initial condition is given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specific within brackets, `[]`, and separated by `,`. Since we only have one species, the vector holds a single element. In this element, we set the value of $X$ using the `:X => 1.0` syntax. Here, we first denote the name of the species (with a `:` pre-appended, i.e. creating a `Symbol`), next follows a `=>` and then the value of $X$. Since we wish to simulate the *concentration* of X over time, we will let the initial condition be decimal valued. +The initial condition is given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specific within brackets, `[]`, and separated by `,`. Since we only have one species, the vector holds a single element. In this element, we set the value of $X$ using the `:X => 1.0` syntax. Here, we first denote the name of the species (with a `:` pre-appended, which creates a `Symbol`), next follows a `=>` and then the value of $X$. Since we wish to simulate the *concentration* of X over time, we will let the initial condition be decimal valued. ```@example ex2 u0 = [:X => 1.0] ``` @@ -107,7 +108,7 @@ params = [:b => 1.0, :d => 0.2] Please read here for more information on [vectors](https://docs.julialang.org/en/v1/manual/arrays/) and [tuples](https://docs.julialang.org/en/v1/manual/types/#Tuple-Types). -Next, before we can simulate our model, we bundle all the required information together in a so-called `ODEProblem`. Note that the order in which the input (the model, the initial condition, the timespan, and the parameter values) is provided to the ODEProblem matters. E.g. the parameter values cannot be provided as the first argument, but have to be the fourth argument. Here, we save our `ODEProblem` in the `oprob` variable. +Next, before we can simulate our model, we bundle all the required information together in a so-called `ODEProblem`. Note that the order in which the input (the model, the initial condition, the timespan, and the parameter values) is provided to the `ODEProblem` matters. E.g. the parameter values cannot be provided as the first argument, but have to be the fourth argument. Here, we save our `ODEProblem` in the `oprob` variable. ```@example ex2 oprob = ODEProblem(rn, u0, tspan, params) ``` @@ -129,7 +130,7 @@ For more information about the numerical simulation package, please see the [Dif ## Additional modelling example To make this introduction more comprehensive, we here provide another example, using a more complicated model. Instead of simulating our model as concentrations evolve over time, we will now simulate the individual reaction events through the [Gillespie algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm) (a common approach for adding *noise* to models). -Remember, unless we have restarted Julia, we do not need to activate our packages (through the `using` command) again. +Remember (unless we have restarted Julia) we do not need to activate our packages (through the `using` command) again. This time, we will declare a so-called [SIR model for an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model). Note that even if this model does not describe a set of chemical reactions, it can be modelled using the same framework. The model consists of 3 species: * $S$, the amount of *susceptible* individuals. @@ -147,13 +148,13 @@ Each reaction is also associated with a specific rate (corresponding to a parame We declare the model using the `@reaction_network` macro, and store it in the `sir_model` variable. ```@example ex2 sir_model = @reaction_network begin - b, S + I --> 2I - k, I --> R + b, S + I --> 2I + k, I --> R end ``` Note that the first reaction contains two different substrates (separated by a `+` sign). While there is only a single product (*I*), two copies of *I* are produced. The *2* in front of the product *I* denotes this. -Next, we declare our initial condition, time span, and parameter values. Since we want to simulate the individual reaction events that discretely change the state of our model, we want our initial conditions to be integer-valued. We will start with a mostly susceptible population, but where a single individual has been infected through some means. +Next, we declare our initial condition, time span, and parameter values. Since we want to simulate the individual reaction events that discretely change the state of our model, we want our initial conditions to be integer-valued. We will start with a mostly susceptible population, but to which a single infected individual has been introduced. ```@example ex2 u0 = [:S => 50, :I => 1, :R => 0] tspan = (0.0, 10.0) @@ -166,10 +167,9 @@ Previously we have bundled this information into an `ODEProblem` (denoting a det dprob = DiscreteProblem(sir_model, u0, tspan, params) jprob = JumpProblem(sir_model, dprob, Direct()) ``` +Again, the order in which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. The last argument to the `JumpProblem` (`Direct()`) denotes which simulation method we wish to use. For now, we recommend that users simply use the `Direct()` option, and then consider alternative ones (see the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/)) when they are more familiar with modelling in Catalyst and Julia. -Again, the order with which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. The last argument to the `JumpProblem` (`Direct()`) denotes which simulation method we wish to use. For now, we recommend that users simply use the `Direct()` option, and then consider alternative ones (see the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/)) when they are more familiar with modelling in Catalyst and Julia. - -Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. Here, the `solve` function also has a second argument (`SSAStepper()`). This is a time stepping algorithm that calls the `Direct` solver to advance a simulation. Again, we recommend at this stage you simply use this option, and then explore exactly what this means at a later stage. +Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. Here, the `solve` function also has a second argument (`SSAStepper()`). This is a time-stepping algorithm that calls the `Direct` solver to advance a simulation. Again, we recommend at this stage you simply use this option, and then explore exactly what this means at a later stage. ```@example ex2 sol = solve(jprob, SSAStepper()) sol = solve(jprob, SSAStepper(); seed=1234) # hide @@ -194,7 +194,7 @@ This will: 2. Switch your current Julia session to use the current folder's environment. !!! note - If you check any folder which has been designated as a Julia environment, it contains a Project.toml and a Manifest.toml file. These store all information regarding the corresponding environment. For non-advanced users, it is recommended to never touch these files directly (and instead do so using various functions from the Pkg package, the important ones which are described in the next two subsections). + If you check any folder which has been designated as a Julia environment, it contains a Project.toml and a Manifest.toml file. These store all information regarding the corresponding environment. For non-advanced users, it is recommended to never touch these files directly (and instead do so using various functions from the Pkg package, the important ones which are described in the next two subsections). ### [Installing and importing packages in Julia](@id catalyst_for_new_julia_users_packages_installing) Package installation and import have been described [previously](@ref catalyst_for_new_julia_users_packages_intro). However, for the sake of this extended tutorial, let us repeat the description by demonstrating how to install the [Latexify.jl](https://github.com/korsbo/Latexify.jl) package (which enables e.g. displaying Catalyst models in Latex format). First, we import the Julia Package manager ([Pkg](https://github.com/JuliaLang/Pkg.jl)) (which is required to install Julia packages): @@ -228,22 +228,23 @@ Pkg.status() So, why is this required, and why cannot we simply import any package installed on our computer? The reason is that most packages depend on other packages, and these dependencies may be restricted to only specific versions of these packages. This creates complicated dependency graphs that restrict what versions of what packages are compatible with each other. When you use `Pkg.add("PackageName")`, only a specific version of that package is actually added (the latest possible version as permitted by the dependency graph). Here, Julia environments both define what packages are available *and* their respective versions (these versions are also displayed by the `Pkg.status()` command). By doing this, Julia can guarantee that the packages (and their versions) specified in an environment are compatible with each other. -The reason why all this is important is that it is *highly recommended* to, for each project, define a separate environment. To these, only add the required packages. General-purpose environments with a large number of packages often produce package-incompatibility issues. While these might not prevent you from installing all desired package, they often mean that you are unable to use the latest version of some packages. +The reason why all this is important is that it is *highly recommended* to, for each project, define a separate environment. To these, only add the required packages. General-purpose environments with a large number of packages often, in the long term, produce package incompatibility issues. While these might not prevent you from installing all desired package, they often mean that you are unable to use the latest version of some packages. !!! note - A not-infrequent cause for reported errors with Catalyst (typically the inability to replicate code in tutorials) is package incompatibilities in large environments preventing the latest version of Catalyst from being installed. Hence, whenever an issue is encountered, it is useful to run `Pkg.status()` to check whenever the latest version of Catalyst is being used. + A not-infrequent cause for reported errors with Catalyst (typically the inability to replicate code in tutorials) is package incompatibilities in large environments preventing the latest version of Catalyst from being installed. Hence, whenever an issue is encountered, it is useful to run `Pkg.status()` to check whenever the latest version of Catalyst is being used. Some additional useful Pkg commands are: - `Pk.rm("PackageName")` removes a package from the current environment. -- `Pkg.update(PackageName")`: updates the designated package. +- `Pkg.update("PackageName")`: updates the designated package. +- `Pkg.update()`: updates all packages. !!! note - A useful feature of Julia's environment system is that enables the exact definition of what packages and versions were used to execute a script. This supports e.g. reproducibility in academic research. Here, by providing the corresponding Project.toml and Manifest.toml files, you can enable someone to reproduce the exact program just to perform some set of analyses. + A useful feature of Julia's environment system is that enables the exact definition of what packages and versions were used to execute a script. This supports e.g. reproducibility in academic research. Here, by providing the corresponding Project.toml and Manifest.toml files, you can enable someone to reproduce the exact program used to perform some set of analyses. --- ## Feedback -If you are a new Julia user who has used this tutorial, and there was something you struggled with or would have liked to have explained better, please [raise an issue](https://github.com/SciML/Catalyst.jl/issues). That way, we can continue improving this tutorial. +If you are a new Julia user who has used this tutorial, and there was something you struggled with or would have liked to have explained better, please [raise an issue](https://github.com/SciML/Catalyst.jl/issues). That way, we can continue improving this tutorial. The same goes for any part of the Catalyst documentation: It is written to help new users understand how to use the package, and if it is not doing so successfully we would like to know! --- ## References diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md index 8f0d530e0c..c579a68aa7 100644 --- a/docs/src/inverse_problems/optimization_ode_param_fitting.md +++ b/docs/src/inverse_problems/optimization_ode_param_fitting.md @@ -1,11 +1,11 @@ # [Parameter Fitting for ODEs using SciML/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. -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. 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/). +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 ref). 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/). ## [Basic example](@id optimization_parameter_fitting_basics) -Let us consider a simple catalysis network, where an enzyme ($E$) turns a substrate ($S$) into a product ($P$): +Let us consider a [Michaelis-Menten enzyme kinetics model](@ref ref), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$): ```@example diffeq_param_estim_1 using Catalyst rn = @reaction_network begin @@ -18,17 +18,17 @@ From some known initial condition, and a true parameter set (which we later want ```@example diffeq_param_estim_1 # Define initial conditions and parameters. u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] -p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] +ps_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] # Generate synthetic data. using OrdinaryDiffEq -oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true) -true_sol = solve(oprob_true, Tsit5()) -data_sol = solve(oprob_true, Tsit5(); saveat=1.0) +oprob_true = ODEProblem(rn, u0, (0.0, 10.0), ps_true) +true_sol = solve(oprob_true) +data_sol = solve(oprob_true; saveat=1.0) data_ts = data_sol.t[2:end] data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] -# Plots the true solutions and the (synthetic data) measurements. +# Plots the true solutions and the (synthetic) data measurements. using Plots plot(true_sol; idxs=:P, label="True solution", lw=8) plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue) @@ -37,22 +37,22 @@ plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color Next, we will use DiffEqParamEstim to build a loss function to measure how well our model's solutions fit the data. ```@example diffeq_param_estim_1 using DiffEqParamEstim, Optimization -p_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0] -oprob = ODEProblem(rn, u0, (0.0, 10.0), p_dummy) +ps_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0] +oprob = ODEProblem(rn, u0, (0.0, 10.0), ps_dummy) loss_function = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); - maxiters=10000, verbose=false, save_idxs=4) + maxiters = 10000, verbose = false, save_idxs = 4) nothing # hide ``` To `build_loss_objective` we provide the following arguments: - `oprob`: The `ODEProblem` with which we simulate our model (using some dummy parameter values, since we do not know these). -- `Tsit5()`: The numeric integrator we wish to simulate our model with. +- `Tsit5()`: The [numeric integrator](@ref ref) we wish to simulate our model with. - `L2Loss(data_ts, data_vals)`: Defines the loss function. While [other alternatives](https://docs.sciml.ai/DiffEqParamEstim/stable/getting_started/#Alternative-Cost-Functions-for-Increased-Robustness) are available, `L2Loss` is the simplest one (measuring the sum of squared distances between model simulations and data measurements). Its first argument is the time points at which the data is collected, and the second is the data's values. - `Optimization.AutoForwardDiff()`: Our choice of [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) framework. Furthermore, we can pass any number of additional optional arguments, these are then passed to the internal `solve()` function (which is used to solve our ODE). Here we provide the following additional arguments: -- `maxiters=10000`: If the ODE integrator takes a very large number of steps, that can be a sign of a very poor fit (or stiffness in the ODEs, but that is not a concern for our current example). Reducing the `maxiters` threshold reduces the time we waste on evaluating such models. -- `verbose=false`: The simulation of models with highly unsuitable parameter sets typically generate various warnings (such as premature simulation termination due to reaching `maxiters` timesteps). To avoid an overflow of such (here unnecessary) warnings, as we evaluate a large number of parameter sets, we turn warnings off. -- `save_idxs=4`: The measured species ($P$) is the 4th species in our species vector (`species(rn)`). Since we only assume data is available for $P(t)$ there is no reason to save any other species. +- `maxiters = 10000`: If the ODE integrator takes a very large number of steps, that can be a sign of a very poor fit (or stiffness in the ODEs, but that is not a concern for our current example). Reducing the `maxiters` threshold reduces the time we waste on evaluating such models. +- `verbose = false`: The simulation of models with highly unsuitable parameter sets typically generate various warnings (such as premature simulation termination due to reaching `maxiters` time steps). To avoid an overflow of such (here unnecessary) warnings, as we evaluate a large number of parameter sets, we turn warnings off. +- `save_idxs = 4`: The measured species ($P$) is the 4th species in our species vector (`species(rn)`). Since data is available for $P(t)$, we will only save the value of this species. Now we can create an `OptimizationProblem` using our `loss_function` and some initial guess of parameter values from which the optimiser will start: ```@example diffeq_param_estim_1 @@ -63,7 +63,7 @@ nothing # hide !!! note `OptimizationProblem` cannot currently accept parameter values in the form of a map (e.g. `[:kB => 1.0, :kD => 1.0, :kP => 1.0]`). These must be provided as individual values (using the same order as the parameters occur in in the `parameters(rs)` vector). Similarly, `build_loss_objective`'s `save_idxs` uses the species' indexes, rather than the species directly. These inconsistencies should be remedied in future DiffEqParamEstim releases. -Finally, we can optimise `optprob` to find the parameter set that best fits our data. Optimization.jl only provides a few optimisation methods natively. However, for each supported optimisation package, it provides a corresponding wrapper-package to import that optimisation package for use with Optimization.jl. E.g., if we wish to use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationOptimJL package. A summary of all, by Optimization.jl, supported optimisation packages can be found [here](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers). Here, we import the Optim.jl package and uses it to minimise our cost function (thus finding a parameter set that fits the data): +Finally, we can optimise `optprob` to find the parameter set that best fits our data. Optimization.jl only provides a few optimisation methods natively. However, for each supported optimisation package, it provides a corresponding wrapper-package to import that optimisation package for use with Optimization.jl. E.g., if we wish to use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationOptimJL package. A summary of all, by Optimization.jl supported, optimisation packages can be found [here](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers). Here, we import the Optim.jl package and uses it to minimise our cost function (thus finding a parameter set that fits the data): ```@example diffeq_param_estim_1 using OptimizationOptimJL optsol = solve(optprob, Optim.NelderMead()) @@ -87,7 +87,7 @@ sol = solve(optprob, NLopt.LD_LBFGS()) nothing # hide ``` -## Optimisation problems with data for multiple species +## [Optimisation problems with data for multiple species](@id optimization_parameter_fitting_multiple_species) Imagine that, in our previous example, we had measurements of the concentration of both *S* and *P*: ```@example diffeq_param_estim_1 data_vals_S = (0.8 .+ 0.4*rand(10)) .* data_sol[:S][2:end] @@ -98,32 +98,32 @@ plot!(data_ts, data_vals_S; label="Measured S", seriestype=:scatter, ms=6, color plot!(data_ts, data_vals_P; label="Measured P", seriestype=:scatter, ms=6, color=:red) ``` -In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1,4]` arguments in `loss_function`: +In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1, 4]` arguments in `loss_function`: ```@example diffeq_param_estim_1 -loss_function_S_P = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, Array(hcat(data_vals_S, data_vals_P)')), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=[1,4]) +loss_function_S_P = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, Array(hcat(data_vals_S, data_vals_P)')), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=[1, 4]) nothing # hide ``` Here, `Array(hcat(data_vals_S, data_vals_P)')` is required to put the data in the right form (in this case, a 2x10 matrix). We can now fit our model to data and plot the results: ```@example diffeq_param_estim_1 -optprob_S_P = OptimizationProblem(loss_function_S_P, [1.0,1.0, 1.0]) +optprob_S_P = OptimizationProblem(loss_function_S_P, [1.0, 1.0, 1.0]) optsol_S_P = solve(optprob_S_P, Optim.NelderMead()) -oprob_fitted_S_P = remake(oprob; p=optsol_S_P.u) -fitted_sol_S_P = solve(oprob_fitted_S_P, Tsit5()) -plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle=:dash, lw=6, color=[:lightblue :pink]) +oprob_fitted_S_P = remake(oprob; p = optsol_S_P.u) +fitted_sol_S_P = solve(oprob_fitted_S_P) +plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink]) ``` -## Setting parameter constraints and boundaries -Sometimes, it is desired to set boundaries on parameter values. Indeed, this can speed up the optimisation process (by preventing searching through unfeasible parts of parameter space), and can also be a requirement for some optimisation methods. This can be done by passing the `lb` (lower bounds) and `up` (upper bounds) arguments to `OptimizationProblem`. These are vectors (of the same length as the number of parameters), with each argument corresponding to the boundary value of the parameter with the same index (as used in the `parameters(rn)` vector). If we wish to constrain each parameter to the interval $(0.1, 10.0)$ this can be done through: +## [Setting parameter constraints and boundaries](@id optimization_parameter_fitting_constraints) +Sometimes, it is desirable to set boundaries on parameter values. Indeed, this can speed up the optimisation process (by preventing searching through unfeasible parts of parameter space), and can also be a requirement for some optimisation methods. This can be done by passing the `lb` (lower bounds) and `up` (upper bounds) arguments to `OptimizationProblem`. These are vectors (of the same length as the number of parameters), with each argument corresponding to the boundary value of the parameter with the same index (as used in the `parameters(rn)` vector). If we wish to constrain each parameter to the interval $(0.1, 10.0)$ this can be done through: ```@example diffeq_param_estim_1 -optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0]; lb = [0.1, 0.1, 0.1], ub = [10.0, 10.0, 10.0]) +optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0]; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1]) nothing # hide ``` In addition to boundaries, Optimization.jl also supports setting [linear and non-linear constraints](https://docs.sciml.ai/Optimization/stable/tutorials/constraints/#constraints) on its output solution for some optimizers. -## Parameter fitting with known parameters +## [Parameter fitting with known parameters](@id optimization_parameter_fitting_known_parameters) If we from previous knowledge know that $kD = 0.1$, and only want to fit the values of $kB$ and $kP$, this can be achieved through `build_loss_objective`'s `prob_generator` argument. First, we create a function (`fixed_p_prob_generator`) that modifies our `ODEProblem` to incorporate this knowledge: ```@example diffeq_param_estim_1 fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2])) @@ -142,8 +142,8 @@ optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead()) nothing # hide ``` -## [Fitting parameters on the logarithmic scale](@id optimization_parameter_fitting_logarithmic_scale) -Often it can be advantageous to fit parameters on a [logarithmic, rather than linear, scale](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646). The best way to proceed is to simply replace each parameter in the model definition by its logarithmic version: +## [Fitting parameters on the logarithmic scale](@id optimization_parameter_fitting_log_scale) +Often it can be advantageous to fit parameters on a [logarithmic scale (rather than on a linear scale)](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646). The best way to do this is to simply replace each parameter in the model definition by its logarithmic version: ```@example diffeq_param_estim_2 using Catalyst rn = @reaction_network begin @@ -152,23 +152,39 @@ rn = @reaction_network begin 10^kP, SE --> P + E end ``` -And then proceeding, by keeping in mind that parameter values are logarithmic. Here, setting +And then going forward, by keeping in mind that parameter values are logarithmic. Here, setting ```@example diffeq_param_estim_2 p_true = [:kB => 0.0, :kD => -1.0, :kP => 10^(0.5)] nothing # hide ``` corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`). -## Parameter fitting to multiple experiments +## [Parameter fitting to multiple experiments](@id optimization_parameter_fitting_multiple_experiments) Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/). -## Optimisation solver options +## [Optimisation solver options](@id optimization_parameter_fitting_solver_options) Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument: ```@example diffeq_param_estim_1 -optsol_fixed_kD = solve(optprob, Optim.NelderMead(); maxtime=100) +optsol_fixed_kD = solve(optprob, Optim.NelderMead(); maxtime = 100) nothing # hide ``` +--- +## [Citation](@id structural_identifiability_citation) +If you use this functionality in your research, please cite the following paper to support the authors of the Optimization.jl package: +``` +@software{vaibhav_kumar_dixit_2023_7738525, + author = {Vaibhav Kumar Dixit and Christopher Rackauckas}, + month = mar, + publisher = {Zenodo}, + title = {Optimization.jl: A Unified Optimization Package}, + version = {v3.12.1}, + doi = {10.5281/zenodo.7738525}, + url = {https://doi.org/10.5281/zenodo.7738525}, + year = 2023 +} +``` + --- ## References [^1]: [Alejandro F. Villaverde, Dilan Pathirana, Fabian Fröhlich, Jan Hasenauer, Julio R. Banga, *A protocol for dynamic model calibration*, Briefings in Bioinformatics (2023).](https://academic.oup.com/bib/article/23/1/bbab387/6383562?login=false) \ No newline at end of file diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md new file mode 100644 index 0000000000..8a1026bbae --- /dev/null +++ b/docs/src/model_creation/dsl_advanced.md @@ -0,0 +1,467 @@ +# [The Catalyst DSL - Advanced Features and Options](@id dsl_advanced_options) +Within the Catalyst DSL, each line can represent either *a reaction* or *an option*. The [previous DSL tutorial](@ref dsl_description) described how to create reactions. This one will focus on options. These are typically used to supply a model with additional information. Examples include the declaration of initial condition/parameter default values, or the creation of observables. + +All option designations begin with a declaration starting with `@`, followed by its input. E.g. the `@observables` option allows for the generation of observables. Each option can only be used once within each use of `@reaction_network`. A full list of options can be found [here](@ref ref), with most (but not all) being described in more detail below. This tutorial will also describe some additional advanced DSL features that do not involve using an option. + +As a first step, we import Catalyst (which is required to run the tutorial): +```@example dsl_advanced_explicit_definitions +using Catalyst +``` + +## [Explicit specification of network species and parameters](@id dsl_advanced_options_declaring_species_and_parameters) +[Previously](@ref ref), we mentioned that the DSL automatically determines which symbols correspond to species and which to parameters. This is done by designating everything that appears as either a substrate or a product as a species, and all remaining quantities as parameters (i.e. those only appearing within rates or [stoichiometric constants](@ref ref)). Sometimes, one might want to manually override this default behaviour for a given symbol. I.e. consider the following model, where the conversion of a protein `P` from its inactive form (`Pᵢ`) to its active form (`Pₐ`) is catalysed by an enzyme (`E`). Using the most natural description: +```@example dsl_advanced_explicit_definitions +catalysis_sys = @reaction_network begin + k*E, Pᵢ --> Pₐ +end +``` +`E` (as well as `k`) will be considered a parameter, something we can confirm directly: +```@example dsl_advanced_explicit_definitions +parameters(catalysis_sys) +``` +If we want `E` to be considered a species, we can designate this using the `@species` option: +```@example dsl_advanced_explicit_definitions +catalysis_sys = @reaction_network begin + @species E(t) + k*E, Pᵢ --> Pₐ +end +parameters(catalysis_sys) +``` +!!! note + When declaring species using the `@species` option, the species symbol must be followed by `(t)`. The reason is that species are time-dependent variables, and this time-dependency must be explicitly specified ([designation of non-`t` dependant species is also possible](@ref ref)). + +Similarly, the `@parameters` option can be used to explicitly designate something as a parameter: +```@example dsl_advanced_explicit_definitions +catalysis_sys = @reaction_network begin + @parameters k + k*E, Pᵢ --> Pₐ +end +``` +Here, while `k` is explicitly defined as a parameter, no information is provided about `E`. Hence, the default case will be used (setting `E` to a parameter). The `@species` and `@parameter` options can be used simultaneously (although a quantity cannot be declared *both* as a species and a parameter). They may be followed by a full list of all species/parameters, or just a subset. + +While designating something which would default to a parameter as a species is straightforward, the reverse (creating a parameter which occurs as a substrate or product) is more involved. This is, however, possible, and described [here](@ref dsl_advanced_options_constant_species). + +Rather than listing all species/parameters on a single line after the options, a `begin ... end` block can be used (listing one species/parameter on each line). E.g. in the following example we use this notation to explicitly designate all species and parameters of the system: +```@example dsl_advanced_explicit_definitions +catalysis_sys = @reaction_network begin + @species begin + E(t) + Pᵢ(t) + Pₐ(t) + end + @parameters begin + k + end + k*E, Pᵢ --> Pₐ +end +``` + +A side-effect of using the `@species` and `@parameter` options is that they specify *the order in which the species and parameters are stored*. I.e. lets check the order of the parameters in the parameters in the following dimerisation model: +```@example dsl_advanced_explicit_definitions +dimerisation = @reaction_network begin + (p,d), 0 <--> X + (kB,kD), 2X <--> X2 +end +parameters(dimerisation) +``` +The default order is typically equal to the order with which the parameters (or species) are encountered in the DSL (this is, however, not guaranteed). If we specify the parameters using `@parameters`, the order used within the option is used instead: +```@example dsl_advanced_explicit_definitions +dimerisation = @reaction_network begin + @parameters kB kD p d + (p,d), 0 <--> X + (kB,kD), 2X <--> X2 +end +parameters(dimerisation) +``` +!!! danger + Generally, Catalyst and the SciML ecosystem *do not* guarantee that parameter and species order are preserved throughout various operations on a model. Writing programs that depend on these orders is *strongly discouraged*. There are, however, some legacy packages which still depend on order (one example is provided [here](@ref ref)). In these situations, this might be useful. However, in these cases, it is recommended that the user is extra wary, and also checks the order manually. + +!!! note + The syntax of the `@species` and `@parameters` options is identical to that used by the `@species` and `@parameters` macros [used in programmatic modelling in Catalyst](@ref programmatic_CRN_construction) (for e.g. designating metadata or initial conditions). Hence, if one has learnt how to specify species/parameters using either approach, that knowledge can be transferred to the other one. + +Generally, there are four main reasons for specifying species/parameters using the `@species` and `@parameters` options: +1. To designate a quantity, that would otherwise have defaulted to a parameter, as a species. +2. To designate default values for parameters/species initial conditions (described [here](@ref dsl_advanced_options_default_vals)). +3. To designate metadata for species/parameters (described [here](@ref dsl_advanced_options_species_and_parameters_metadata)). +4. To designate a species or parameters that do not occur in reactions, but are still part of the model (e.g a [parametric initial condition](@ref dsl_advanced_options_parametric_initial_conditions)) + +!!!! warn + Catalyst's DSL automatically infer species and parameters from the input. However, it only does so for *quantities that appear in reactions*. Until now this has not been relevant. However, this tutorial will demonstrate cases where species/parameters that are not part of reactions are used. These *must* be designated using either the `@species` or `@parameters` options (or the `@variables` option, which is described [later](@ref ref)). + +### [Setting default values for species and parameters](@id dsl_advanced_options_default_vals) +When declaring species/parameters using the `@species` and `@parameters` options, one can also assign them default values (by appending them with `=` followed by the desired default value). E.g here we set `X`'s default initial condition value to $1.0$, and `p` and `d`'s default values to $1.0$ and $0.2$, respectively: +```@example dsl_advanced_defaults +using Catalyst # hide +rn = @reaction_network begin + @species X(t)=1.0 + @parameters p=1.0 d=0.1 + (p,d), 0 <--> X +end +``` +Next, if we simulate the model, we do not need to provide values for species or parameters that have default values. In this case all have default values, so both `u0` and `ps` can be empty vectors: +```@example dsl_advanced_defaults +using OrdinaryDiffEq, Plots +u0 = [] +tspan = (0.0, 10.0) +p = [] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob) +plot(sol) +``` +It is still possible to provide values for some (or all) initial conditions/parameters in `u0`/`ps` (in which case these overrides the default values): +```@example dsl_advanced_defaults +u0 = [:X => 4.0] +p = [:d => 0.5] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob) +plot(sol) +``` +It is also possible to declare a model with default values for only some initial conditions/parameters: +```@example dsl_advanced_defaults +using Catalyst # hide +rn = @reaction_network begin + @species X(t)=1.0 + (p,d), 0 <--> X +end + +tspan = (0.0, 10.0) +p = [:p => 1.0, :D => 0.2] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob) +plot(sol) +``` +API for checking the default values of species and parameters can be found [here](@ref ref). + +### [Setting parametric initial conditions](@id dsl_advanced_options_parametric_initial_conditions) +In the previous section, we designated default values for initial conditions and parameters. However, the right-hand side of the designation accepts any valid expression (not only numeric values). While this can be used to set up some advanced default values, the most common use case is to designate a species's initial condition as a parameter. E.g. in the following example we represent the initial condition of `X` using the parameter `X₀`. +```@example dsl_advanced_defaults +rn = @reaction_network begin + @species X(t)=X₀ + @parameters X₀ + (p,d), 0 <--> X +end +``` +Please note that as the parameter `X₀` does not occur as part of any reactions, Catalyst's DSL cannot infer whether it is a species or a parameter. This must hence be explicitly declared. We can now simulate our model while providing `X`'s value through the `X₀` parameter: +```@example dsl_advanced_defaults +u0 = [] +p = [:X₀ => 1.0, :p => 1.0, :d => 0.5] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob, Tsit5()) +plot(sol) +``` +It is still possible to designate $X$'s value in `u0`, in which case this overrides the default value. +```@example dsl_advanced_defaults +u0 = [:X => 0.5] +p = [:X₀ => 1.0, :p => 1.0, :d => 0.5] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob, Tsit5()) +plot(sol) +``` +Please note that `X₀` is still a parameter of the system, and as such its value must still be designated to simulate the model (even if it is not actually used). + +### [Designating metadata for species and parameters](@id dsl_advanced_options_species_and_parameters_metadata) +Catalyst permits the user to define *metadata* for species and parameters. This permits the user to assign additional information to these, which can be used for a variety of purposes. Some Catalyst features depend on using metadata (with each such case describing specifically how this is done). Here we will introduce how to set metadata, and describe some common metadata types. + +Whenever a species/parameter is declared using the `@species`/`@parameters` options, it can be followed by a `[]` within which the metadata is given. Each metadata entry consists of the metadata's name, followed by a `=`, followed by its value. E.g. the `description` metadata allows you to attach a [`String`](https://docs.julialang.org/en/v1/base/strings/) to a species/parameter. Here we create a simple model where we add descriptions to all species and parameters. +```@example dsl_advanced_metadata +using Catalyst # hide +two_state_system = @reaction_network begin + @species Xi(t) [description="The X's inactive form"] Xa(t) [description="The X's active form"] + @parameters kA [description="X's activation rate"] kD [description="X's deactivation rate"] + (ka,kD), Xi <--> Xa +end +``` +A metadata can be given to only a subset of a system's species/parameters, and a quantity can be given several metadata entries. To give several metadata, separate each by a `,`. Here we only provide a description for `kA`, for which we also provide a [bounds metadata](@ref https://docs.sciml.ai/ModelingToolkit/dev/basics/Variable_metadata/#Bounds), +```@example dsl_advanced_metadata +two_state_system = @reaction_network begin + @parameters kA [description="X's activation rate", bound=(0.01,10.0)] + (ka,kD), Xi <--> Xa +end +``` + +It is possible to add both default values and metadata to a parameter/species. In this case, first provide the default value, and next the metadata. I.e. to in the above example set $kA$'s default value to $1.0$ we use +```@example dsl_advanced_metadata +two_state_system = @reaction_network begin + @parameters kA=1.0 [description="X's activation rate", bound=(0.01,10.0)] + (ka,kD), Xi <--> Xa +end +``` + +When designating metadata for species/parameters in `begin ... end` blocks the syntax changes slightly. Here, a `,` must be inserted before the metadata (but after any potential default value). I.e. a version of the previous example can be written as +```@example dsl_advanced_metadata +two_state_system = @reaction_network begin + @parameters begin + kA, [description="X's activation rate", bound=(0.01,10.0)] + kD = 1.0, [description="X's deactivation rate"] + end + (ka,kD), Xi <--> Xa +end +``` + +Each metadata has its own getter functions. E.g. we can get the description of the parameter `kA` using `getdescription` (here we use [system indexing](@ref ref) to access the parameter): +```@example dsl_advanced_metadata +getdescription(two_state_system.kA) +``` + +It is not possible for the user to directly designate their own metadata. These have to first be added to Catalyst. Doing so is somewhat involved, and described in detail [here](@ref ref). A full list of metadata that can be used for species and/or parameters can be found [here](@ref ref). + +### [Designating constant-valued/fixed species parameters](@id dsl_advanced_options_constant_species) + +Catalyst enables the designation of parameters as `constantspecies`. These parameters can be used as species in reactions, however, their values are not changed by the reaction and remain constant throughout the simulation (unless changed by e.g. the [occurrence of an event]@ref ref). Practically, this is done by setting the parameter's `isconstantspecies` metadata to `true`. Here, we create a simple reaction where the species `X` is converted to `Xᴾ` at rate `k`. By designating `X` as a constant species parameter, we ensure that its quantity is unchanged by the occurrence of the reaction. +```@example dsl_advanced_constant_species +using Catalyst # hide +rn = @reaction_network begin + @parameters X [isconstantspecies=true] + k, X --> Xᴾ +end +``` +We can confirm that $Xᴾ$ is the only species of the system: +```@example dsl_advanced_constant_species +species(rn) +``` +Here, the produced model is actually identical to if $X$ had simply been a parameter in the reaction's rate: +```@example dsl_advanced_constant_species +rn = @reaction_network begin + k*X, 0 --> Xᴾ +end +``` + +A common use-case for constant species is when modelling systems where some species are present in such surplus that their amounts the reactions' effect on it is negligible. A system which is commonly modelled this way is the [Brusselator](https://en.wikipedia.org/wiki/Brusselator). + +### [Designating parameter types](@id dsl_advanced_options_parameter_types) +Sometimes it is desired to designate that a parameter should have a specific [type](@ref ref). When supplying this parameter's value to e.g. an `ODEProblem`, that parameter will then be restricted to that specific type. Designating a type is done by appending the parameter with `::` followed by its type. E.g. in the following example we specify that the parameter `n` (the number of `X` molecules in the `Xn` polymer) must be an integer (`Int64`) +```@example dsl_advanced_parameter_types +using Catalyst # hide +polymerisation_network = @reaction_network begin + @parameters n::Int64 + (kB,kD), n*X <--> Xn +end +nothing # hide +``` +Generally, when simulating models with mixed parameter types, it is recommended to [declare parameter values as tuples, rather than vectors](@ref ref), e.g.: +```@example dsl_advanced_parameter_types +ps = (:kB => 0.2, :kD => 1.0, :n => 2) +nothing # hide +``` + +If a parameter has a type, metadata, and a default value, they are designated in the following order: +```@example dsl_advanced_parameter_types +polymerisation_network = @reaction_network begin + @parameters n::Int64 = 2 [description="Parameter n, which is an integer and defaults to the value 2."] + (kB,kD), n*X <--> Xn +end +nothing # hide +``` + +### [Vector-valued species or parameters](@id dsl_advanced_options_vector_variables) +Sometimes, one wishes to declare a large number of similar parameters or species. This can be done by *creating them as vectors*. E.g. below we create a [two-state system](@ref ref). However, instead of declaring `X1` and `X2` (and `k1` and `k2`) as separate entities, we declare them as vectors: +```@example dsl_advanced_vector_variables +using Catalyst # hide +two_state_model = @reaction_network begin + @parameters k[1:2] + @species X(t)[1:2] + (k[1],k[2]), X[1] <--> X[2] +end +``` +Now, we can also declare our initial conditions and parameter values as vectors as well: +```@example dsl_advanced_vector_variables +using OrdinaryDiffEq, Plots # hide +u0 = [:X => [0.0, 2.0]] +tspan = (0.0, 1.0) +ps = [:k => [1.0, 2.0]] +oprob = ODEProblem(two_state_model, u0, tspan, ps) +sol = solve(oprob) +plot(sol) +``` + +## [Naming reaction networks](@id dsl_advanced_options_naming) +Each reaction network model has a name. It can be accessed using the `nameof` function. By default, some generic name is used: +```@example dsl_advanced_names +using Catalyst # hide +rn = @reaction_network begin + (p,d), 0 <--> X +end +nameof(rn) +``` +A specific name can be given as an argument between the `@reaction_network` and the `begin`. E.g. to name a network `my_network` we can use: +```@example dsl_advanced_names +rn = @reaction_network my_network begin + (p,d), 0 <--> X +end +nameof(rn) +``` + +A consequence of generic names being used by default is that networks, even if seemingly identical, by default are not. E.g. +```@example dsl_advanced_names +rn1 = @reaction_network begin + (p,d), 0 <--> X +end +rn2 = @reaction_network begin + (p,d), 0 <--> X +end +rn1 == rn2 +``` +The reason can be confirmed by checking that their respective (randomly generated) names are different: +```@example dsl_advanced_names +nameof(rn1) == nameof(rn2) +``` +By designating the networks to have the same name, however, identity is achieved. +```@example dsl_advanced_names +rn1 = @reaction_network my_network begin + (p,d), 0 <--> X +end +rn2 = @reaction_network my_network begin + (p,d), 0 <--> X +end +rn1 == rn2 +``` + +Setting model names is primarily useful for [hierarchical modelling](@ref ref), where network names are appended to the display names of subnetworks' species and parameters. + +## [Creating observables](@id dsl_advanced_options_observables) +Sometimes one might want to use observable variables. These are variables with values that can be computed directly from a system's state (rather than having their values implicitly given by reactions or equations). Observables can be designated using the `@observables` option. Here, the `@observables` option is followed by a `begin ... end` block with one line for each observable. Each line first gives the observable, followed by a `~` (*not* a `=`!), followed by an expression describing how to compute it. + +Let us consider a model where two species (`X` and `Y`) can bind to form a complex (`XY`, which also can dissociate back into `X` and `Y`). If we wish to create a representation for the total amount of `X` and `Y` in the system, we can do this by creating observables `Xtot` and `Ytot`: +```@example dsl_advanced_observables +using Catalyst # hide +rn = @reaction_network begin + @observables begin + Xtot ~ X + XY + Ytot ~ Y + XY + end + (kB,kD), X + Y <--> XY +end +``` +We can now simulate our model using normal syntax (initial condition values for observables should not, and can not, be provided): +```@example dsl_advanced_observables +using OrdinaryDiffEq +u0 = [:X => 1.0, :Y => 2.0, :XY => 0.0] +tspan = (0.0, 10.0) +ps = [:kB => 1.0, :kD => 1.5] +oprob = ODEProblem(rn, u0, tspan, ps) +sol = solve(oprob, Tsit5()) +nothing # hide +``` + +Next, we can use [symbolic indexing](@ref simulation_structure_interfacing) of our solution object, but with the observable as input. E.g. we can use +```@example dsl_advanced_observables +sol[:Xtot] +``` +to get a vector with `Xtot`'s value throughout the simulation. We can also use +```@example dsl_advanced_observables +using Plots +plot(sol; idxs = [:Xtot, :Ytot]) +``` +to plot the observables (rather than the species). + +Observables can be defined using complicated expressions containing species, parameters, and [variables](@ref ref) (but not other observables). In the following example (which uses a [parametric stoichiometry](@ref ref)) `X` polymerises to form a complex `Xn` containing `n` copies of `X`. Here, we create an observable describing the total number of `X` molecules in the system: +```@example dsl_advanced_observables +rn = @reaction_network begin + @observables Xtot ~ X + n*Xn + (kB,kD), n*X <--> Xn +end +nothing # hide +``` +!!! + If only a single observable is declared, the `begin .. end` block is not required and the observable can be declared directly after the `@observables` option. + +[Metadata](@ref dsl_advanced_options_species_and_parameters_metadata) can be supplied to an observable directly after its declaration (but before its formula). If so, the metadata must be separated from the observable with a `,`, and the observable plus the metadata encapsulated by `()`. E.g. to add a [description metadata](@ref dsl_advanced_options_species_and_parameters_metadata) to our observable we can use +```@example dsl_advanced_observables +rn = @reaction_network begin + @observables (Xtot, [description="The total amount of X in the system."]) ~ X + n*Xn + (kB,kD), n*X <--> Xn +end +nothing # hide +``` + +Observables are by default considered [variables](@ref ref) (not species). To designate them as a species, they can be pre-declared using the `@species` option. I.e. Here `Xtot` becomes a species: +```@example dsl_advanced_observables +rn = @reaction_network begin + @species Xtot(t) + @observables Xtot ~ X + n*XnXY + (kB,kD), n*X <--> Xn +end +nothing # hide +``` + +Some final notes regarding observables: +- The left-hand side of the observable declaration must contain a single symbol only (with the exception of metadata, which can also be supplied). +- All quantities appearing on the right-hand side must be declared elsewhere within the `@reaction_network` call (either by being part of a reaction, or through the `@species`, `@parameters`, or `@variables` options). +- Observables may not depend on other observables. +- Observables have their [dependent variable(s)](@ref ref) automatically assigned as the union of the dependent variables of the species and variables on which it depends. + +## [Specifying non-time independent variables](@id dsl_advanced_options_ivs) + +As [described elsewhere](@ref ref), Catalyst's `ReactionSystem` models depend on a *time independent variable*, and potentially one or more *spatial independent variables*. By default, the independent variable `t` is used. We can declare another independent variable (which is automatically used as the default one) using the `@ivs` option. E.g. to use `τ` instead of `t` we can use +```@example dsl_advanced_ivs +using Catalyst # hide +rn = @reaction_network begin + @ivs τ + (ka,kD), Xi <--> Xa +end +nothing # hide +``` +We can confirm that `Xi` and `Xa` depend on `τ` (and not `t`): +```@example dsl_advanced_ivs +species(rn) +``` + +It is possible to designate several independent variables using `@ivs`. If so, the first one is considered the default (time) independent variable, while the following one(s) are considered spatial independent variable(s). If we want some species to depend on a non-default independent variable, this has to be explicitly declared: +```@example dsl_advanced_ivs +rn = @reaction_network begin + @ivs τ x + @species X(τ) Y(x) + (p1,d1), 0 <--> X + (p2,d2), 0 <--> Y +end +species(rn) +``` +It is also possible to have species which depends on several independent variables: +```@example dsl_advanced_ivs +rn = @reaction_network begin + @ivs t x + @species Xi(t,x) Xa(t,x) + (ka,kD), Xi <--> Xa +end +species(rn) +``` + +!!! note + Setting spatial independent variables is primarily intended for the modelling of spatial systems on continuous domains. Catalyst's support for this is currently under development. Hence, the utility of specifying spatial independent variables is limited. + + +## [Setting reaction metadata](@id dsl_advanced_options_reaction_metadata) +It is possible to supply reactions with *metadata*, containing some additional information of the reaction. A reaction's metadata follows after its declaration (first using the metadata's name, then a `=`, then its value) and is encapsulated by `[]` (where individual entries are separated by `,`). Here, we add a `description` metadata to the reactions of a birth-death process: +```@example dsl_advanced_reaction_metadata +using Catalyst # hide +bd_model = @reaction_network begin + p, 0 --> X, [description="A production reaction"] + d, X --> 0, [description="A degradation reaction"] +end +nothing # hide +``` + +When [bundling reactions](@ref dsl_description_reaction_bundling), reaction metadata can be bundled using the same rules as rates. Bellow we re-declare our birth-death process, but on a single line: +```@example dsl_advanced_reaction_metadata +bd_model = @reaction_network begin + (p,d), 0 --> X, ([description="A production reaction"], [description="A degradation reaction"]) +end +nothing # hide +``` + +Here we declare a model where we also provide a `misc` metadata (which can hold any quantity we require) to our birth reaction: +```@example dsl_advanced_reaction_metadata +bd_model = @reaction_network begin + p, 0 --> X, [description="A production reaction", misc=:value] + d, X --> 0, [description="A degradation reaction"] +end +nothing # hide +``` + +A reaction's metadata can be accessed using specific functions, e.g. `Catalyst.hasdescription` and `Catalyst.getdescription` can be used to check if a reaction have a description metadata, and to retrieve it, respectively: +```@example dsl_advanced_reaction_metadata +rx = @reaction p, 0 --> X, [description="A production reaction"] +Catalyst.getdescription(rx) +``` + +A list of all available reaction metadata can be found [here](@ref ref). \ No newline at end of file diff --git a/docs/src/model_creation/dsl_basics.md b/docs/src/model_creation/dsl_basics.md new file mode 100644 index 0000000000..19057d2d65 --- /dev/null +++ b/docs/src/model_creation/dsl_basics.md @@ -0,0 +1,372 @@ +# [The Catalyst DSL - Introduction](@id dsl_description) +In the [introduction to Catalyst](@ref introduction_to_catalyst) we described how the `@reaction_network` [macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros) can be used to create chemical reaction network (CRN) models. This macro enables a so-called [domain-specific language](https://en.wikipedia.org/wiki/Domain-specific_language) (DSL) for creating CRN models. This tutorial will give a basic introduction on how to create Catalyst models using this macro (from now onwards called "*the Catalyst DSL*"). A [follow-up tutorial](@ref dsl_advanced_options) will describe some of the DSL's more advanced features. + +The Catalyst DSL generates a [`ReactionSystem`](@ref) (the [julia structure](https://docs.julialang.org/en/v1/manual/types/#Composite-Types) Catalyst uses to represent CRN models). These can be created through alternative methods (e.g. [programmatically](@ref programmatic_CRN_construction) or [compositionally](@ref compositional_modeling)). A summary of the various ways to create `ReactionSystems`s can be found [here](@ref ref). [Previous](@ref ref) and [following](@ref ref) tutorials describe how to simulate models once they have been created using the DSL. This tutorial will solely focus on model creation. + +Before we begin, we will first load the Catalyst package (which is required to run the code). +```@example dsl_basics_intro +using Catalyst +``` + +### [Quick-start summary](@id dsl_description_quick_start) +The DSL is initiated through the `@reaction_network` macro, which is followed by one line for each reaction. Each reaction consists of a *rate*, followed lists first of the substrates and next of the products. E.g. a [Michaelis-Menten enzyme kinetics system](@ref ref) can be written as +```@example dsl_basics_intro +rn = @reaction_network begin + (kB,kD), S + E <--> SE + kP, SE --> P + E +end +``` +Here, `<-->` is used to create a bi-directional reaction (with forward rate `kP` and backward rate `kD`). Next, the model (stored in the variable `rn`) can be used as input to various types of [simulations](@ref ref). + +## [Basic syntax](@id dsl_description_basic_syntax) +The basic syntax of the DSL is +```@example dsl_basics +rn = @reaction_network begin + 2.0, X --> Y + 1.0, Y --> X +end +``` +Here, you start with `@reaction_network begin`, next list all of the model's reactions, and finish with `end`. Each reaction consists of +- A *rate*. +- A (potentially empty) set of *substrates*. +- A (potentially empty) set of *products*. + +Each reaction line declares, in order, the rate, the substrate(s), and the product(s). The rate is separated from the substrate(s) by a `,`, and the substrate(s) from the production by a `-->` (other arrows, however, are [also possible](@ref dsl_description_symbols_arrows)). In the above example, our model consists of two reactions. In the first one, `X` (the single substrate) becomes `Y` (the single product) at rate `2.0`. In the second reaction, `Y` becomes `X` at rate `1.0`. + +Finally, `rn = ` is used to store the model in the variable `rn` (a normal Julia variable, which does not need to be called `rn`). + +## [Defining parameters and species in the DSL](@id dsl_description_parameters_basics) +Typically, the rates are not constants, but rather parameters (which values can be set e.g. at [the beginning of each simulation](@ref ref)). To set parametric rates, simply use whichever symbol you wish to represent your parameter with. E.g. to set the above rates to `a` and `b`, we use: +```@example dsl_basics +rn1 = @reaction_network begin + a, X --> Y + b, Y --> X +end +``` + +Here we have used single-character symbols to designate all species and parameters. Multi-character symbols, however, are also permitted. E.g. we could call the rates `kX` and `kY`: +```@example dsl_basics +rn1 = @reaction_network begin + kX, X --> Y + kY, Y --> X +end +nothing # hide +``` +Generally, anything that is a [permitted Julia variable name](@id https://docs.julialang.org/en/v1/manual/variables/#man-allowed-variable-names) can be used to designate a species or parameter in Catalyst. + +## [Different types of reactions](@id dsl_description_reactions) + +### [Reactions with multiple substrates or products](@id dsl_description_reactions_multiples) +Previously, our reactions have had a single substrate and a single product. However, reactions with multiple substrates and/or products are possible. Here, all the substrates (or products) are listed and separated by a `+`. E.g. to create a model where `X` and `Y` bind (at rate `kB`) to form `XY` (which then can dissociate, at rate `kD`, to form `XY`) we use: +```@example dsl_basics +rn2 = @reaction_network begin + kB, X + Y --> XY + kD, XY --> X + Y +end +``` +Reactions can have any number of substrates and products, and their names do not need to have any relationship to each other, as demonstrated by the following mock model: +```@example dsl_basics +rn3 = @reaction_network begin + k, X + Y + Z --> A + B + C + D +end +``` + +### [Reactions with degradation or production](@id dsl_description_reactions_degradation_and_production) +Some reactions have no products, in which case the substrate(s) are degraded (i.e. removed from the system). To denote this, set the reaction's right-hand side to `0`. Similarly, some reactions have no substrates, in which case the product(s) are produced (i.e. added to the system). This is denoted by setting the left-hand side to `0`. E.g. to create a model where a single species `X` is both created (in the first reaction) and degraded (in a second reaction), we use: +```@example dsl_basics +rn4 = @reaction_network begin + p, 0 --> X + d, X --> 0 +end +``` + +### [Reactions with non-unitary stoichiometries](@id dsl_description_reactions_stoichiometries) +Reactions may include multiple copies of the same reactant (i.e. a substrate or a product). To specify this, the reactant is preceded by a number indicating its number of copies (also called the reactant's *stoichiometry*). E.g. to create a model where two copies of `X` dimerise to form `X2` (which then dissociate back to two `X` copies) we use: +```@example dsl_basics +rn5 = @reaction_network begin + kB, 2X --> X2 + kD, X2 --> 2X +end +``` +Reactants whose stoichiometries are not defined are assumed to have stoichiometry `1`. Any integer number can be used, furthermore, [decimal numbers and parameters can also be used as stoichiometries](@ref ref). A discussion of non-unitary (i.e. not equal to `1`) stoichiometries affecting the created model can be found [here](@ref ref). + +Stoichiometries can be combined with `()` to define them for multiple reactants. Here, the following (mock) model declares the same reaction twice, both with and without this notation: +```@example dsl_basics +rn6 = @reaction_network begin + k, 2X + 3(Y + 2Z) --> 5(V + W) + k, 2X + 3Y + 6Z --> 5V + 5W +end +nothing # hide +``` + +## [Bundling of similar reactions](@id dsl_description_reaction_bundling) + +### [Bi-directional (or reversible) reactions](@id dsl_description_reaction_bundling_reversible) +As is the case for the following two-state model: +```@example dsl_basics +rn7 = @reaction_network begin + k1, X1 --> X2 + k2, X2 --> X1 +end +``` +it is common that reactions occur in both directions (so-called *bi-directional* reactions). Here, it is possible to bundle the reactions into a single line by using the `<-->` arrow. When we do this, the rate term must include two separate rates (one for each direction, these are enclosed by a `()` and separated by a `,`). I.e. the two-state model can be declared using: +```@example dsl_basics +rn7 = @reaction_network begin + (k1,k2), X1 <--> X2 +end +``` +Here, the first rate (`k1`) denotes the *forward rate* and the second rate (`k2`) the *backwards rate*. + +Catalyst also permits writing pure backwards reactions. These use identical syntax to forward reactions, but with the `<--` arrow: +```@example dsl_basics +rn8 = @reaction_network begin + k, X <-- Y +end +``` +Here, the substrate(s) are on the right-hand side and the product(s) are on the left-hand side. Hence, the above model can be written identically using: +```@example dsl_basics +rn8 = @reaction_network begin + k, Y --> X +end +``` +Generally, using forward reactions is clearer than backwards ones, with the latter typically never being used. + +### [Bundling similar reactions on a single line](@id dsl_description_reaction_bundling_similar) +There exist additional situations where models contain similar reactions (e.g. systems where all system components degrade at identical rates). Reactions which share either rates, substrates, or products can be bundled into a single line. Here, the parts which are different for the reactions are written using `(,)` (containing one separate expression for each reaction). E.g., let us consider the following model where species `X` and `Y` both degrade at the rate `d`: +```@example dsl_basics +rn8 = @reaction_network begin + d, X --> 0 + d, Y --> 0 +end +``` +These share both their rates (`d`) and products (`0`), however, the substrates are different (`X` and `Y`). Hence, the reactions can be bundled into a single line using the common rate and product expression while providing separate substrate expressions: +```@example dsl_basics +rn8 = @reaction_network begin + d, (X,Y) --> 0 +end +``` +This declaration of the model is identical to the previous one. Reactions can share any subset of the rate, substrate, and product expression (the cases where they share all or none, however, do not make sense to use). I.e. if the two reactions also have different degradation rates: +```@example dsl_basics +rn9 = @reaction_network begin + dX, X --> 0 + dY, Y --> 0 +end +``` +This can be represented using: +```@example dsl_basics +rn9 = @reaction_network begin + (dX,dY), (X,Y) --> 0 +end +``` + +It is possible to use bundling for any number of reactions. E.g. in the following model we bundle the conversion of a species $X$ between its various forms (where all reactions use the same rate $k$): +```@example dsl_basics +rn10 = @reaction_network begin + k, (X0,X1,X2,X3) --> (X1,X2,X3,X4) +end +``` + +It is possible to combine bundling with bi-directional reactions. In this case, the rate is first split into the forward and backwards rates. These may then (or may not) indicate several rates. We exemplify this using the two following two (identical) networks, created with and without bundling. +```@example dsl_basics +rn11 = @reaction_network begin + kf, S --> P1 + kf, S --> P2 + kb_1, P1 --> S + kb_2, P2 --> S +end +``` +```@example dsl_basics +rn11 = @reaction_network begin + (kf, (kb_1, kb_2)), S <--> (P1,P2) +end +``` + +Like when we designated stoichiometries, reaction bundling can be applied very generally to create some truly complicated reactions: +```@example dsl_basics +rn12 = @reaction_network begin + ((pX, pY, pZ),d), (0, Y0, Z0) <--> (X, Y, Z1+Z2) +end +``` +However, like for the above model, bundling reactions too zealously can reduce (rather than improve) a model's readability. + +## [Non-constant reaction rates](@id dsl_description_nonconstant_rates) +So far we have assumed that all reaction rates are constant (being either a number of a parameter). Non-constant rates that depend on one (or several) species are also possible. More generally, the rate can be any valid expression of parameters and species. + +Let us consider a model with an activator (`A`, which degraded at a constant rate) and a protein (`P`). The production rate of `P` depends both on `A` and a parameter (`kP`). We model this through: +```@example dsl_basics +rn_13 = @reaction_network begin + d, A --> 0 + kP*A, 0 --> P +end +``` +Here, `P`'s production rate will be reduced as `A` decays. We can [print the ODE this model produces with `Latexify`](@ref ref): +```@example dsl_basics +using Latexify +latexify(rn_13; form=:ode) +``` + +In this case, we can generate an equivalent model by instead adding `A` as both a substrate and a product to `P`'s production reaction: +```@example dsl_basics +rn_13_alt = @reaction_network begin + d, A --> 0 + kp, A --> A + P +end +``` +We can confirm that this generates the same ODE: +```@example dsl_basics +latexify(rn_13_alt; form=:ode) +``` +Here, while these models will generate identical ODE, SDE, and jump simulations, the chemical reaction network models themselves are not equivalent. Generally, as pointed out in the two notes below, using the second form is preferable. +!!! warn + While `rn_13` and `rn_13_alt` will generate equivalent simulations, for jump simulations, the first model will have [reduced performance](@ref ref) (which generally are more performant when rates are constant). + +!!! danger + Catalyst automatically infers whether quantities appearing in the DSL are species or parameters (as described [here](@ref dsl_advanced_options_declaring_species_and_parameters)). Generally, anything that does not appear as a reactant is inferred to be a parameter. This means that if you want to model a reaction activated by a species (e.g. `kp*A, 0 --> P`), but that species does not occur as a reactant, it will be interpreted as a parameter. This can be handled by [manually declaring the system species](@ref dsl_advanced_options_declaring_species_and_parameters). A full example of how to do this for this example can be found [here](@ref ref). + +Above we used a simple example where the rate was the product of a species and a parameter. However, any valid Julia expression of parameters, species, and values can be used. E.g the following is a valid model: +```@example dsl_basics +rn_14 = @reaction_network begin + 2.0 + X^2, 0 --> X + Y + k1 + k2^k3, X --> ∅ + pi * X/(sqrt(2) + Y), Y → ∅ +end +``` + +### [Using functions in rates](@id dsl_description_nonconstant_rates_functions) +It is possible for the rate to contain Julia functions. These can either be functions from Julia's standard library: +```@example dsl_basics +rn_16 = @reaction_network begin + d, A --> 0 + kp*sqrt(A), 0 --> P +end +``` +or ones defined by the user: +```@example dsl_basics +custom_function(p1, p2, X) = (p1 + X) / (p2 + X) +rn_17 = @reaction_network begin + d, A --> 0 + custom_function(k1,k2,E), 0 --> P +end +``` + +### [Pre-defined functions](@id dsl_description_nonconstant_rates_available_functions) +Two functions frequently used within systems biology are the [*Michaelis-Menten*](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) and [*Hill*](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)) functions. These are pre-defined in Catalyst and can be called using `mm(X,v,K)` and `hill(X,v,K,n)`. E.g. a self-activation loop where `X` activates its own production through a Hill function can be created using: +```@example dsl_basics +rn_18 = @reaction_network begin + hill(X,v,K,n), 0 --> P + d, X --> 0 +end +``` + +Catalyst comes with the following predefined functions: +- The Michaelis-Menten function: $mm(X,v,K) = v * X/(X + K)$. +- The repressive Michaelis-Menten function: $mmr(X,v,K) = v * K/(X + K)$. +- The Hill function: $hill(X,v,K,n) = v * (X^n)/(X^n + K^n)$. +- The repressive Hill function: $hillr(X,v,K,n) = v * (K^n)/(X^n + K^n)$. +- The activating/repressive Hill function: $hillar(X,Y,v,K,n) = v * (X^n)/(X^n + Y^n + K^n)$. + +### [Time-dependant rates](@id dsl_description_nonconstant_rates_time) +Previously we have assumed that the rates are independent of the [time variable, $t$](@ref ref). However, time-dependent reactions are also possible. Here, simply use `t` to represent the time variable. E.g., to create a production/degradation model where the production rate decays as time progresses, we can use: +```@example dsl_basics +rn_14 = @reaction_network begin + kp/(1 + t), 0 --> P + d, P --> 0 +end +``` + +Like previously, `t` can be part of any valid expression. E.g. to create a reaction with a cyclic rate (e.g. to represent a [circadian system](https://en.wikipedia.org/wiki/Circadian_rhythm)) we can use: +```@example dsl_basics +rn_15 = @reaction_network begin + A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P + d, P --> 0 +end +``` + +!!! warn + Jump simulations cannot be performed for models with time-dependent rates without additional considerations, which are discussed [here](@ref ref). + +## [Non-standard stoichiometries](@id dsl_description_stoichiometries) + +### [Non-integer stoichiometries](@id dsl_description_stoichiometries_decimal) +Previously all stoichiometric constants have been integer numbers, however, decimal numbers are also permitted. Here we create a birth-death model where each production reaction produces 1.5 units of `X`: +```@example dsl_basics +rn_16 = @reaction_network begin + p, 0 --> 1.5X + d, X --> 0 +end +``` +It is also possible to have non-integer stoichiometric coefficients for substrates. However, in this case the [`combinatoric_ratelaw = false`](@ref ref) option must be used. We note that non-integer stoichiometric coefficients do not make sense in most fields, however, this feature is available for use for models where it does make sense. + +### [Parametric stoichiometries](@id dsl_description_stoichiometries_parameters) +It is possible for stoichiometric coefficients to be parameters. E.g. here we create a generic polymerisation system where `n` copies of `X` bind to form `Xn`: +```@example dsl_basics +rn_17 = @reaction_network begin + (kB,kD), n*X <--> Xn +end +``` +Now we can designate the value of `n` through a parameter when we e.g. create an `ODEProblem`: +```@example dsl_basics +u0 = [:X => 5.0, :Xn => 1.0] +ps = [:kB => 1.0, :kD => 0.1, :n => 4] +oprob = ODEProblem(rn_17, u0, (0.0, 1.0), ps) +nothing # hide +``` + +## [Using special symbols](@id dsl_description_symbols) +Julia permits any Unicode characters to be used in variable names, thus Catalyst can use these as well. Below we describe some cases where this can be useful. No functionality is, however, tied to this. + +### [Using ∅ in degradation/production reactions](@id dsl_description_symbols_empty_set) +Previously, we described how `0` could be used to [create degradation or production reactions](@ref dsl_description_reactions_degradation_and_production). Catalyst permits the user to instead use the `∅` symbol. E.g. the production/degradation system can alternatively be written as: +```@example dsl_basics +rn4 = @reaction_network begin + p, ∅ --> X + d, X --> ∅ +end +``` + +### [Using special arrow symbols](@id dsl_description_symbols_arrows) +Catalyst uses `-->`, `<-->`, and `<--` to denote forward, bi-directional, and backwards reactions, respectively. Several unicode representations of these arrows are available. Here, +- `>`, `→`, `↣`, `↦`, `⇾`, `⟶`, `⟼`, `⥟`, `⥟`, `⇀`, and `⇁` can be used to represent forward reactions. +- `↔`, `⟷`, `⇄`, `⇆`, `⇌`, `⇋`, , and `⇔` can be used to represent bi-directional reactions. +- `<`, `←`, `↢`, `↤`, `⇽`, `⟵`, `⟻`, `⥚`, `⥞`, `↼`, , and `↽` can be used to represent backwards reactions. + +E.g. the production/degradation system can alternatively be written as: +```@example dsl_basics +rn4 = @reaction_network begin + p, ∅ → X + d, X → ∅ +end +``` + +### [Using special symbols to denote species or parameters](@id dsl_description_symbols_special) +A range of possible characters are available which can be incorporated into species and parameter names. This includes, but is not limited to: +- Greek letters (e.g `α`, `σ`, `τ`, and `Ω`). +- Superscript and subscript characters (to create e.g. `k₁`, `k₂`, `Xₐ`, and `Xᴾ`). +- Non-latin, non-greek, letters (e.g. `ä`, `Д`, `س`, and `א`). +- Other symbols (e.g. `£`, `ℂ`, `▲`, and `♠`). + +An example of how this can be used to create a neat-looking model can be found in [Schwall et al. (2021)](https://www.embopress.org/doi/full/10.15252/msb.20209832) where it was used to model a sigma factor V circuit in the bacteria *Bacillus subtilis*: +```@example dsl_basics +σᵛ_model = @reaction_network begin + v₀ + hill(σᵛ,v,K,n), ∅ → σᵛ + A + kdeg, (σᵛ, A, Aσᵛ) → ∅ + (kB,kD), A + σᵛ ↔ Aσᵛ + L, Aσᵛ → σᵛ +end +nothing # hide +``` + +This functionality can also be used to create less serious models: +rn_13 = @reaction_network begin + 🍦, 😢 --> 😃 +end + +It should be noted that the following symbols are *not permitted* to be used as species or parameter names: +- `pi` and `π` (used in Julia to denote [`3.1415926535897...`](https://en.wikipedia.org/wiki/Pi)). +- `ℯ` (used in Julia to denote [Euler's constant](https://en.wikipedia.org/wiki/Euler%27s_constant)). +- `t` (used to denote the [time variable](@ref dsl_description_nonconstant_rates_time)). +- `∅` ([used for production/degradation reactions](@ref dsl_description_symbols_empty_set)). +- `im` (used in Julia to represent [complex numbers](https://docs.julialang.org/en/v1/manual/complex-and-rational-numbers/#Complex-Numbers)). +- `nothing` (used in Julia to denote [nothing](https://docs.julialang.org/en/v1/base/constants/#Core.nothing)). +- `Γ` (used by Catalyst to represent [conserved quantities](@ref ref)). + diff --git a/docs/src/model_creation/dsl_description.md b/docs/src/model_creation/dsl_description.md deleted file mode 100644 index e32bba7673..0000000000 --- a/docs/src/model_creation/dsl_description.md +++ /dev/null @@ -1,767 +0,0 @@ -# [The Reaction DSL](@id dsl_description) -This tutorial describes the syntax for building chemical reaction -network models using Catalyst's domain-specific language (DSL). Examples showing -how to both construct and solve ODE, SDE, and jump models are provided in [Basic -Chemical Reaction Network Examples](@ref basic_CRN_library). To learn more about the symbolic -[`ReactionSystem`](@ref)s generated by the DSL, and how to use them directly, see -the tutorial on [Programmatic Construction of Symbolic Reaction Systems](@ref programmatic_CRN_construction). - -We first load the `Catalyst` package, which is required for the code in this tutorial to run -```@example tut2 -using Catalyst -``` - -## [Basic syntax](@id basic_examples) - -The `@reaction_network` macro allows the (symbolic) specification of reaction -networks with a simple format. Its input is a set of chemical reactions, and -from them it generates a symbolic [`ReactionSystem`](@ref) reaction network -object. The `ReactionSystem` can be used as input to ModelingToolkit -`ODEProblem`, `NonlinearProblem`, `SteadyStateProblem`, `SDEProblem`, -`JumpProblem`, and more. `ReactionSystem`s can also be incrementally extended as -needed, allowing for programmatic construction of networks and network -composition. - -The basic syntax is: - -```@example tut2 -rn = @reaction_network begin - 2.0, X + Y --> XY - 1.0, XY --> Z1 + Z2 -end -``` -where each line of the [`@reaction_network`](@ref) macro corresponds to a -chemical reaction. Each reaction consists of a reaction rate (the expression on -the left-hand side of `,`), a set of substrates (the expression in-between `,` -and `-->`), and a set of products (the expression on the right-hand side of -`-->`). The substrates and the products may contain one or more reactants, -separated by `+`. The naming convention for these is the same as for normal -variables in Julia. - -The chemical reaction model is generated by the `@reaction_network` macro and -stored in the `rn` variable (a normal Julia variable, which does not need to be -called `rn`). It corresponds to a [`ReactionSystem`](@ref), a symbolic -representation of the chemical network. The generated `ReactionSystem` can be -converted to a symbolic differential equation model via -```@example tut2 -osys = convert(ODESystem, rn) -osys = complete(osys) -``` - -We can then convert the symbolic ODE model into a compiled, optimized -representation for use in the SciML ODE solvers by constructing an `ODEProblem`. -Creating an `ODEProblem` also requires our specifying the initial conditions for -the model. We do this by creating a mapping from each symbolic variable -representing a chemical species to its initial value -```@example tut2 -# define the symbolic variables -t = default_t() -@species X(t) Y(t) Z(t) XY(t) Z1(t) Z2(t) - -# create the mapping -u0 = [X => 1.0, Y => 1.0, XY => 1.0, Z1 => 1.0, Z2 => 1.0] -``` -Alternatively, we can create a mapping using Julia `Symbol`s for each variable, -and then convert them to a mapping involving symbolic variables like -```@example tut2 -u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]) -``` -Given the mapping, we can then create an `ODEProblem` from our symbolic `ODESystem` -```@example tut2 -tspan = (0.0, 1.0) # the time interval to solve on -oprob = ODEProblem(osys, u0, tspan, []) -``` - -Catalyst provides a shortcut to avoid having to explicitly `convert` to an -`ODESystem` and/or use `symmap_to_varmap`, allowing direct construction -of the `ODEProblem` like -```@example tut2 -u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0] -oprob = ODEProblem(rn, u0, tspan, []) -``` - -For more detailed examples, see the example implementations in the [Library of Basic Chemical Reaction Network Models](@ref basic_CRN_library) section. - -## Defining parameters and species -Numeric parameter values do not need to be set when the model is created, i.e. -Catalyst supports symbolic parameters too: -```@example tut2 -rn = @reaction_network begin - k1, X --> Y - k2, Y --> X -end -``` -All symbols that do not appear as a substrate or product in a reaction are -designated by Catalyst as a parameter (i.e. all symbols appearing only within -rate expressions and/or as [stoichiometric coefficients](@ref parametric_stoichiometry)). In this example `X` and `Y` -appear as a substrates and products, but neither `k1` nor `k2`. Hence `k1` and `k2` are -designated as parameters. Later in this tutorial, we will describe how to manually specify what should be -considered a species or parameter. - -## Production, Destruction, and Stoichiometry -Sometimes reactants are produced/destroyed from/to nothing. This can be -designated using either `0` or `∅`: -```@example tut2 -rn = @reaction_network begin - 2.0, 0 --> X - 1.0, X --> 0 -end -``` -If several molecules of the same reactant are involved in a reaction, the -stoichiometry of a reactant in a reaction can be set using a number. Here, two -molecules of species `X` form the dimer `X2`: -```@example tut2 -rn = @reaction_network begin - 1.0, 2X --> Y -end -``` -this corresponds to the differential equation: -```@example tut2 -convert(ODESystem, rn) -``` -Other numbers than 2 can be used, and parenthesis can be used to reuse the same -stoichiometry for several reactants: -```@example tut2 -rn = @reaction_network begin - 1.0, X + 2(Y + Z) --> W -end -``` -Note, one can explicitly multiply by integer coefficients too, i.e. -```@example tut2 -rn = @reaction_network begin - 1.0, X + 2*(Y + Z) --> W -end -``` - -## Arrow variants -A variety of Unicode arrows are accepted by the DSL in addition to `-->`. All of -these work: `>`, `→` `↣`, `↦`, `⇾`, `⟶`, `⟼`, `⥟`, `⥟`, `⇀`, `⇁`. Backwards -arrows can also be used to write the reaction in the opposite direction. For -example, these reactions are equivalent: -```@example tut2 -rn = @reaction_network begin - 1.0, X + Y --> XY - 1.0, X + Y → XY - 1.0, XY ← X + Y - 1.0, XY <-- X + Y -end -``` - -## Bi-directional arrows for reversible reactions -Bi-directional arrows, including bidirectional Unicode arrows like ↔, can be -used to designate a reversible reaction. For example, these two models are -equivalent: -```@example tut2 -rn = @reaction_network begin - 2.0, X + Y --> XY - 2.0, X + Y <-- XY -end -``` -```@example tut2 -rn2 = @reaction_network begin - (2.0,2.0), X + Y <--> XY -end -``` - -If the reaction rates in the backward and forward directions are different, they -can be designated in the following way: -```@example tut2 -rn = @reaction_network begin - (2.0,1.0), X + Y <--> XY -end -``` -which is identical to -```@example tut2 -rn = @reaction_network begin - 2.0, X + Y --> XY - 1.0, X + Y <-- XY -end -``` - -## Combining several reactions in one line -Several similar reactions can be combined in one line by providing a tuple of -reaction rates and/or substrates and/or products. If several tuples are provided, -they must all be of identical length. These pairs of reaction networks are all -identical. - -Pair 1: -```@example tut2 -rn1 = @reaction_network begin - 1.0, S --> (P1,P2) -end -``` -```@example tut2 -rn2 = @reaction_network begin - 1.0, S --> P1 - 1.0, S --> P2 -end -``` -Pair 2: -```@example tut2 -rn1 = @reaction_network begin - (1.0,2.0), (S1,S2) --> P -end -``` -```@example tut2 -rn2 = @reaction_network begin - 1.0, S1 --> P - 2.0, S2 --> P -end -``` -Pair 3: -```@example tut2 -rn1 = @reaction_network begin - (1.0,2.0,3.0), (S1,S2,S3) --> (P1,P2,P3) -end -``` -```@example tut2 -rn2 = @reaction_network begin - 1.0, S1 --> P1 - 2.0, S2 --> P2 - 3.0, S3 --> P3 -end -``` -This can also be combined with bi-directional arrows, in which case separate -tuples can be provided for the backward and forward reaction rates. -These reaction networks are identical -```@example tut2 -rn1 = @reaction_network begin - (1.0,(1.0,2.0)), S <--> (P1,P2) -end -``` -```@example tut2 -rn2 = @reaction_network begin - 1.0, S --> P1 - 1.0, S --> P2 - 1.0, P1 --> S - 2.0, P2 --> S -end -``` - -## Variable reaction rates -Reaction rates do not need to be a single parameter or a number, but can also be -expressions depending on time or the current amounts of system species (when, for -example, one species can activate the production of another). For instance, this -is a valid notation: -```@example tut2 -rn = @reaction_network begin - 1.0, X --> ∅ - k*X, Y --> ∅ -end -``` -corresponding to the ODE model -```@example tut2 -convert(ODESystem,rn) -``` - -With respect to the corresponding mass action ODE model, this is actually -equivalent to the reaction system -```@example tut2 -rn = @reaction_network begin - 1.0, X --> ∅ - k, X + Y --> X -end -``` -```@example tut2 -convert(ODESystem,rn) -``` -!!! note - While the ODE models corresponding to the preceding two reaction systems are - identical, in the latter example the `Reaction` stored in `rn` will be classified as - [`ismassaction`](@ref) while in the former it will not, which can impact optimizations - used in generating `JumpSystem`s. For this reason, it is recommended to use the - latter representation when possible. - -Most expressions and functions are valid reaction rates, e.g.: -```@example tut2 -using SpecialFunctions -rn = @reaction_network begin - 2.0*X^2, 0 --> X + Y - t*gamma(Y), X --> ∅ - pi*X/Y, Y --> ∅ -end -``` -where here `t` always denotes Catalyst's time variable. Please note that many -user-defined functions can be called directly, but others will require -registration with Symbolics.jl ([see the faq](@ref user_functions)). - -## [Explicit specification of network species and parameters](@id dsl_description_explicit_species) -Recall that the `@reaction_network` macro automatically designates symbols used -in the macro as either parameters or species, with symbols that appear as a -substrate or product being species, and all other symbols becoming parameters -(i.e. those that only appear within a rate expression and/or as [stoichiometric coefficients](@ref parametric_stoichiometry)). Sometimes, one might want to manually override this default -behavior for a given symbol. E.g one might want something to be considered as a -species, even if it only appears within a rate expression. In the following -network -```@example tut2 -rn = @reaction_network begin - k*X, Y --> 0 -end -``` -`X` (as well as `k`) will be considered a parameter. - -By using the `@species` and `@parameters` options within the `@reaction_network` -macro, one can manually declare that specified symbols should be considered a -species or parameter. E.g in: -```@example tut2 -rn = @reaction_network begin - @species X(t) Y(t) - k*X, Y --> 0 -end -``` -`X` and `Y` are set as species. Please note that when declaring species using -the `@species` option, their dependant variable (almost always `t`) also needs -to be designated. Similarly in -```@example tut2 -rn = @reaction_network begin - @parameters k - k*X, Y --> 0 -end -``` -both `X` and `k` will be considered as parameters. It is also possible to use -both options simultaneously, allowing users to fully specify which symbols are -species and/or parameters: -```@example tut2 -rn = @reaction_network begin - @species X(t) Y(t) - @parameters k - k*X, Y --> 0 -end -``` -Here, `X` and `Y` are designated as species and `k` as a parameter. - -The lists provided to the `@species` and `@parameters` options do not need to be extensive. Any symbol that appears in neither list will use the default option as determined by the macro. E.g. in the previous example, where we only want to change the default designation of `X` (making it a species rather than a parameter), we can simply write: -```@example tut2 -rn = @reaction_network begin - @species X(t) - k*X, Y --> 0 -end -``` - -Finally, note that the `@species` and `@parameters` options can also be used in -`begin ... end` block form, allowing more formatted lists of species/parameters: -```@example tut2 -rn = @reaction_network begin - @parameters begin - d1 - d2 - end - @species begin - X1(t) - X2(t) - end - d2, X2 --> 0 - d1, X1 --> 0 -end -``` -This can be especially useful when declaring default values for clarity of model -specification (see the next section). - -## [Setting default values for initial conditions and parameters](@id dsl_description_defaults) -When using the `@species` and ` @parameters` macros to declare species and/or -parameters, one can also provide default initial conditions for each species and -values for each parameter: -```@example tut2 -rn = @reaction_network begin - @species X(t)=1.0 - @parameters p=1.0 d=0.1 - p, 0 --> X - d, X --> ∅ -end -``` -This system can now be simulated without providing initial condition or -parameter vectors to the DifferentialEquations.jl solvers: -```@example tut2 -using DifferentialEquations, Plots -u0 = [] -tspan = (0.0, 10.0) -p = [] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob) -plot(sol) -``` - -When providing default values, it is possible to do so for only a subset of the -species or parameters, in which case the rest can be specified when constructing -the problem type to solve: -```@example tut2 -rn = @reaction_network begin - @species X(t) - @parameters p=1.0 d - p, 0 --> X - d, X --> 0 -end - -u0 = [:X => 1.0] -tspan = (0.0, 10.0) -p = [:d => .1] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob) -plot(sol) -``` - -Finally, default values can be overridden by passing mapping vectors to the -DifferentialEquations.jl problem being constructed. Only those initial conditions -or parameters for which we want to change their value from the default will need to be passed -```@example tut2 -u0 = [:X => 1.0] -tspan = (0.0, 10.0) -p = [:p => 2.0, :d => .1] # we change p to 2.0 -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob) -plot(sol) -``` - -## Constant/fixed species -It is possible to fix the amount of a species in a reaction. Without fixing -a species, a reaction could look like -```@example tut2 -rn = @reaction_network begin - k, X + Y --> 0 -end -``` - -```@example tut2 -ode_sys = convert(ODESystem, rn) -``` - -```@example tut2 -equations(ode_sys) -``` - -Fixing a species could either be achieved by modifying the reaction specification -and specifying constant species explicitly as species as described -[above](@ref dsl_description_explicit_species), i.e., -```@example tut2 -rn = @reaction_network begin - @species X(t) - k * X, Y --> 0 -end -``` - -```@example tut2 -ode_sys = convert(ODESystem, rn) -``` - -```@example tut2 -equations(ode_sys) -``` - -The species can of course also just be used as parameter - using the same modification -of the reaction, i.e., -```@example tut2 -rn = @reaction_network begin - k * X, Y --> 0 -end -``` - -```@example tut2 -ode_sys = convert(ODESystem, rn) -``` - -```@example tut2 -equations(ode_sys) -``` - -The same result can also be achieved by declaring a species as fixed/constant -without having to change the reaction itself, i.e., -```@example tut2 -rn = @reaction_network begin - @parameters X [isconstantspecies = true] - k, X + Y --> 0 -end -``` - -```@example tut2 -ode_sys = convert(ODESystem, rn) -``` - -```@example tut2 -equations(ode_sys) -``` - -## [Setting initial conditions that depend on parameters](@id dsl_description_parametric_initial_conditions) -It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system: -```@example tut2 -rn = @reaction_network begin - @parameters X0 - @species X(t)=X0 - p, 0 --> X - d, X --> ∅ -end -``` -We can now simulate the network without providing any initial conditions: -```@example tut2 -u0 = [] -tspan = (0.0, 10.0) -p = [:p => 2.0, :d => .1, :X0 => 1.0] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob) -plot(sol) -``` - -## Naming the generated `ReactionSystem` -ModelingToolkit uses system names to allow for compositional and hierarchical -models. To specify a name for the generated `ReactionSystem` via the -[`@reaction_network`](@ref) macro, just place the name before `begin`: -```@example tut2 -rn = @reaction_network production_degradation begin - p, ∅ --> X - d, X --> ∅ -end -ModelingToolkit.nameof(rn) == :production_degradation -``` - -## Pre-defined functions -Hill functions and a Michaelis-Menten function are pre-defined and can be used -as rate laws. Below, the pair of reactions within `rn1` are equivalent, as are -the pair of reactions within `rn2`: -```@example tut2 -rn1 = @reaction_network begin - hill(X,v,K,n), ∅ --> X - v*X^n/(X^n+K^n), ∅ --> X -end -``` -```@example tut2 -rn2 = @reaction_network begin - mm(X,v,K), ∅ --> X - v*X/(X+K), ∅ --> X -end -``` -Repressor Hill (`hillr`) and Michaelis-Menten (`mmr`) functions are also -provided: -```@example tut2 -rn1 = @reaction_network begin - hillr(X,v,K,n), ∅ --> X - v*K^n/(X^n+K^n), ∅ --> X -end -``` -```@example tut2 -rn2 = @reaction_network begin - mmr(X,v,K), ∅ --> X - v*K/(X+K), ∅ --> X -end -``` - -Please see the API [Rate Laws](@ref api_rate_laws) section for more details. - -## Including non-species variables -Non-species unknown variables can be specified in the DSL using the `@variables` -macro. These are declared similarly to species. For example, -```@example tut2 -rn_with_volume = @reaction_network begin - @variables V(t) - k*V, 0 --> A -end -``` -creates a network with one species -```@example tut2 -species(rn_with_volume) -``` -and one non-species -```@example tut2 -nonspecies(rn_with_volume) -``` -giving two unknown variables, always internally ordered by species and then -nonspecies: -```@example tut2 -unknowns(rn_with_volume) -``` - -`rn_with_volume` could then be extended with constraint equations for how `V(t)` -evolves in time, see the [associated tutorial](@ref constraint_equations). - -## Specifying alternative time variables and/or extra independent variables -While the DSL defaults to allowing `t` as the time variable, one can use the -`@ivs` macro to specify an alternative independent variable. For example, to -make `s` the default time variable one can say -```@example tut2 -rn_with_s = @reaction_network begin - @ivs s - @variables V(s) - @species B(s) - k, A + V*B --> C -end -show(stdout, MIME"text/plain"(), rn_with_s) # hide -``` -where we see all unknowns are now functions of `s`. - -Similarly, if one wants unknowns to be functions of more than one independent -variable, for example to encode a spatial problem, one can list more than one -variable, i.e. `@ivs t x y`. Here the first listed independent variable is -always chosen to represent time. For example, -```@example tut2 -rn_with_many_ivs = @reaction_network begin - @ivs s x - @variables V1(s) V2(s,x) - @species A(s) B(s,x) - k, V1*A --> V2*B + C -end -show(stdout, MIME"text/plain"(), rn_with_many_ivs) # hide -``` -Here again `s` will be the time variable, and any inferred species, `C` in this -case, are made functions of both variables, i.e. `C(s, x)`. - -## [Interpolation of Julia variables](@id dsl_description_interpolation_of_variables) -The DSL allows Julia variables to be interpolated for the network name, within -rate constant expressions, or for species/stoichiometry within reactions. Using -the lower-level symbolic interface we can then define symbolic variables and -parameters outside of the macro, which can then be used within expressions in -the DSL (see the [Programmatic Construction of Symbolic Reaction Systems](@ref programmatic_CRN_construction) -tutorial for details on the lower-level symbolic interface). For example, -```@example tut2 -t = default_t() -@parameters k α -@species A(t) -spec = A -par = α -rate = k*A -name = :network -rn = @reaction_network $name begin - $rate*B, 2*$spec + $par*B --> $spec + C - end -``` -As the parameters `k` and `α` were pre-defined and appeared via interpolation, -we did not need to declare them within the `@reaction_network` macro, -i.e. they are automatically detected as parameters: -```@example tut2 -parameters(rn) -``` -as are the species coming from interpolated variables -```@example tut2 -species(rn) -``` - -!!! note - When using interpolation, expressions like `2$spec` won't work; the - multiplication symbol must be explicitly included like `2*$spec`. - -## Including observables -Sometimes, one might want to include observable variables. These are variables that can be computed directly from the other system variables (rather than having their values implicitly given through some differential equation). These can be introduced through the `@observables` option. - -Let us consider a simple example where two species ($X$ and $Y$) are produced and degraded at constant rates. They can also bind, forming a complex ($XY$). If we want to access the total amount of $X$ in the system we can create an observable that denotes this quantity ($Xtot = X + XY$). Here, we create observables for the total amount of $X$ and $Y$: -```@example obs1 -using Catalyst # hide -rn = @reaction_network begin - @observables begin - Xtot ~ X + XY - Ytot ~ Y + XY - end - (pX,dX), 0 <--> X - (pY,dY), 0 <--> Y - (kB,kD), X + Y <--> XY -end -``` -The `@observables` option is followed by one line for each observable formula (enclosed by a `begin ... end` block). The left-hand sides indicate the observables' names, and the right-hand sides how their values are computed. The two sides are separated by a `~`. - -If we now simulate our model: -```@example obs1 -using DifferentialEquations # hide -u0 = [:X => 0.0, :Y => 0.0, :XY => 0.0] -tspan = (0.0, 10.0) -ps = [:pX => 1.0, :dX => 0.2, :pY => 1.0, :dY => 0.5, :kB => 1.0, :kD => 0.2] -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(oprob) -nothing # hide -``` -we can index the solution using our observables (just like for [other variables](@ref simulation_structure_interfacing_solutions)). E.g. we can receive a vector with all $Xtot$ values using -```@example obs1 -sol[:Xtot] -``` -similarly, we can plot the values of $Xtot$ and $Ytot$ using -```@example obs1 -using Plots -plot(sol; idxs = [rn.Xtot, rn.Ytot], label = ["Total X" "Total Y"]) -``` - -If we only wish to provide a single observable, the `begin ... end` block is note required. E.g., to record only the total amount of $X$ we can use: -```@example obs1 -using Catalyst # hide -rn = @reaction_network begin - @observables Xtot ~ X + XY - (pX,dX), 0 <--> X - (pY,dY), 0 <--> Y - (kB,kD), X + Y <--> XY -end -``` - -Finally, some general rules for creating observables: -- Observables can depend on any species, parameters, or variables, but not on other observables. -- All observables components appearing on the right side of the `~` must be declared somewhere (i.e., they cannot only appear as a part of the observables formula). -- Only a single `@observables` option block can be used in each `@reaction_network` call. -- The left-hand side of the observables expression must be a single symbol, indicating the observable's name. -- Metadata can, however, be provided, e.g through `@observables (Xtot, [description="Total amount of X"]) ~ X + XY`. -- The right-hand side of the observables expression can be any valid algebraic expression. -- Observables are (by default, but this can be changed) considered `variables` (and not `species`). This can be changed by e.g. pre-declaring them using the `@species` option: -```@example obs2 -using Catalyst # hide -rn = @reaction_network begin - @species Xtot(t) - @observables Xtot ~ X1 + X2 - (k1,k2), X1 <--> X2 -end -nothing # hide -``` - -## Incorporating (differential) equations into reaction network models -Some models cannot be purely described as reaction networks. E.g. consider the growth of a cell, where the rate of change in the cell's volume depends on some growth factor. Here, the cell's volume would be described by a normal ODE. Such equations can be incorporated into a model using the `@equations` option. Here, we create a model where a growth factor ($G$) is produced and degraded at a linear rates, and the rate of change in cell volume ($V$) is linear in the amount of growth factor: -```@example eqs1 -using Catalyst #hide -rn = @reaction_network begin - @equations begin - D(V) ~ G - end - (p,d), 0 <--> G -end -``` -Here, `D(V)` indicates the (time) derivative with respect to `D`. The differential equation left and right hand sides are separated by a `~`. The left-hand side should contain differential only, the right hand side can contain any algebraic expression. - -We can check the differential equation corresponding to this reaction network using latexify: -```@example eqs1 -using Latexify -latexify(rn; form=:ode) -``` -We can also simulate it using the normal syntax -```@example eqs1 -using DifferentialEquations, Plots # hide -u0 = [:G => 0.0, :V => 0.1] -ps = [:p => 1.0, :d => 0.5] -oprob = ODEProblem(rn, u0, (0.0, 1.0), ps) -sol = solve(oprob) -plot(sol) -``` -Here, growth is indefinite. To improve the model, [a callback](@ref advanced_simulations_callbacks) can be used to half the volume (cell division) once some threshold is reached. - -When creating differential equations this way, the subject of the differential is automatically inferred to be a variable, however, any component on the right-hand side must be declare somewhere in the macro. E.g. to add a scaling parameter ($k$), we must declare that $k$ is a parameter using the `@parameters` option: -```@example eqs1 -rn = @reaction_network begin - @parameters k - @equations begin - D(V) ~ k*G - end - (p,d), 0 <--> G -end -nothing #hide -``` -If the differential does not appear isolated on the lhs, its subject variable must also be explicitly declared (as it is not inferred for these cases). - -It is possible to add several equations to the model. In this case, each have a separate line. E.g. to keep track of a supply of nutrition ($N$) in the growth media, we can use: -```@example eqs1 -rn = @reaction_network begin - @equations begin - D(V) ~ G - D(N) ~ -G - end - (p,d), 0 <--> G -end -nothing #hide -``` - -When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using: -```@example eqs1 -rn = @reaction_network begin - @equations D(V) ~ G - (p,d), 0 <--> G -end -nothing # hide -``` diff --git a/docs/src/model_creation/examples/basic_CRN_examples.md b/docs/src/model_creation/examples/basic_CRN_examples.md deleted file mode 100644 index 762761cbbc..0000000000 --- a/docs/src/model_creation/examples/basic_CRN_examples.md +++ /dev/null @@ -1,274 +0,0 @@ -# [Library of Basic Chemical Reaction Network Models](@id basic_CRN_library) -Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the `@reaction_network` DSL, and basic simulations are performed. - -## Birth-death process -The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates: -```@example crn_library_birth_death -using Catalyst -bd_process = @reaction_network begin - (p,d), ∅ <--> X -end -``` -Next we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations). -```@example crn_library_birth_death -u0 = [:X => 1] -tspan = (0.0, 10.0) -ps = [:p => 1, :d => 0.2] -``` -We can now simulate our model using all three interpretations. First we perform a reaction rate equation-based ODE simulation: -```@example crn_library_birth_death -using OrdinaryDiffEq -oprob = ODEProblem(bd_process, u0, tspan, ps) -osol = solve(oprob, Tsit5()) -nothing # hide -``` -Next, a chemical Langevin equation-based SDE simulation: -```@example crn_library_birth_death -using StochasticDiffEq -sprob = SDEProblem(bd_process, u0, tspan, ps) -ssol = solve(sprob, ImplicitEM()) -nothing # hide -``` -Next, a stochastic chemical kinetics-based jump simulation: -```@example crn_library_birth_death -using JumpProcesses -dprob = DiscreteProblem(bd_process, u0, tspan, ps) -jprob = JumpProblem(bd_process, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -nothing # hide -``` -Finally, we plot the results: -```@example crn_library_birth_death -using Plots -oplt = plot(osol; title = "Reaction rate equation (ODE)") -splt = plot(ssol; title = "Chemical Langevin equation (SDE)") -jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)") -plot(oplt, splt, jplt; size=(800,700), layout = (3,1)) -``` - -## Michaelis-Menten enzyme kinetics -[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions it can be simplified to a singe function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model: -```@example crn_library_michaelis_menten -using Catalyst -mm_system = @reaction_network begin - kB, S + E --> SE - kD, SE --> S + E - kP, SE --> P + E -end -``` -Next, we perform ODE, SDE, and jump simulations of the model: -```@example crn_library_michaelis_menten -u0 = [:S => 301, :E => 100, :SE => 0, :P => 0] -tspan = (0., 100.) -ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1] - -using OrdinaryDiffEq -oprob = ODEProblem(mm_system, u0, tspan, ps) -osol = solve(oprob, Tsit5()) - -using StochasticDiffEq -sprob = SDEProblem(mm_system, u0, tspan, ps) -ssol = solve(sprob, ImplicitEM()) - -using JumpProcesses -dprob = DiscreteProblem(mm_system, u0, tspan, ps) -jprob = JumpProblem(mm_system, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) - -using Plots -oplt = plot(osol; title = "Reaction rate equation (ODE)") -splt = plot(ssol; title = "Chemical Langevin equation (SDE)") -jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)") -plot(oplt, splt, jplt; size=(800,700), layout = (3,1)) -``` - -## SIR infection model -The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery. -```@example crn_library_sir -using Catalyst -sir_model = @reaction_network begin - α, S + I --> 2I - β, I --> R -end -``` -First we perform a deterministic ODE simulation: -```@example crn_library_sir -using OrdinaryDiffEq, Plots -u0 = [:S => 99, :I => 1, :R => 0] -tspan = (0.0, 500.0) -ps = [:α => 0.001, :β => 0.01] - -# Solve ODEs. -oprob = ODEProblem(sir_model, u0, tspan, ps) -osol = solve(oprob, Tsit5()) -plot(osol; title = "Reaction rate equation (ODE)") -``` -Next we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of a outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak. -```@example crn_library_sir -using JumpProcesses -dprob = DiscreteProblem(sir_model, u0, tspan, ps) -jprob = JumpProblem(sir_model, dprob, Direct()) - -jsol1 = solve(jprob, SSAStepper()) -jsol2 = solve(jprob, SSAStepper()) -jsol3 = solve(jprob, SSAStepper()) -jsol1 = solve(jprob, SSAStepper(); seed=1) # hide -jsol2 = solve(jprob, SSAStepper(); seed=2) # hide -jsol3 = solve(jprob, SSAStepper(); seed=3) # hide - -jplt1 = plot(jsol1; title = "Outbreak") -jplt2 = plot(jsol2; title = "Outbreak") -jplt3 = plot(jsol3; title = "No outbreak") -plot(jplt1, jplt2, jplt3; size=(800,700), layout = (3,1)) -``` - -## The Wilhelm model -The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability. -```@example crn_library_wilhelm -wilhelm_model = @reaction_network begin - k1, Y --> 2X - k2, 2X --> X + Y - k3, X + Y --> Y - k4, X --> 0 -end -``` -We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states. -```@example crn_library_wilhelm -using OrdinaryDiffEq, Plots -u0_1 = [:X => 1.5, :Y => 0.5] -u0_2 = [:X => 2.5, :Y => 0.5] -tspan = (0., 10.) -ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5] - -oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps) -oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps) -osol1 = solve(oprob1, Tsit5()) -osol2 = solve(oprob2, Tsit5()) -oplt1 = plot(osol1; idxs = :X, label = "X(0) = 1.5") -oplt2 = plot!(osol2; idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,700)) -``` - -## Simple self-activation loop -The simplest self-activation loop consist of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability. -```@example crn_library_self_activation -using Catalyst -sa_loop = @reaction_network begin - v₀ + hill(X,v,K,n), ∅ --> X - d, X --> ∅ -end -``` -A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor. - -We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states. -```@example crn_library_self_activation -using OrdinaryDiffEq, Plots -u0 = [:X => 4] -tspan = (0.0, 1000.0) -ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1] - -oprob = ODEProblem(sa_loop, u0, tspan, ps) -osol = solve(oprob, Tsit5()) - -dprob = DiscreteProblem(sa_loop, u0, tspan, ps) -jprob = JumpProblem(sa_loop, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(); seed = 2091) # hide - -plot(osol; label = "Reaction rate equation (ODE)") -plot!(jsol; label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,600)) -``` - -## The Brusselator -The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator"). -```@example crn_library_brusselator -using Catalyst -brusselator = @reaction_network begin - A, ∅ --> X - 1, 2X + Y --> 3X - B, X --> Y - 1, X --> ∅ -end -``` -It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this results is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contain a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in he second case the system is able to generate oscillations. -```@example crn_library_brusselator -using OrdinaryDiffEq, Plots -u0 = [:X => 1.0, :Y => 1.0] -tspan = (0., 50.) -ps1 = [:A => 1.0, :B => 1.0] -ps2 = [:A => 1.0, :B => 1.8] - -oprob1 = ODEProblem(brusselator, u0, tspan, ps1) -oprob2 = ODEProblem(brusselator, u0, tspan, ps2) -osol1 = solve(oprob1, Rodas5P()) -osol2 = solve(oprob2, Rodas5P()) -oplt1 = plot(osol1; title = "No Oscillation") -oplt2 = plot(osol2; title = "Oscillation") - -plot(oplt1, oplt2; layout = (1,2), size(800,700)) -``` - -## The repressilator -The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that is able to generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) which production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)). -```@example crn_library_brusselator -using Catalyst -repressilator = @reaction_network begin - hillr(Z,v,K,n), ∅ --> X - hillr(X,v,K,n), ∅ --> Y - hillr(Y,v,K,n), ∅ --> Z - d, (X, Y, Z) --> ∅ -end -``` -Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model is able to sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*. -```@example crn_library_brusselator -using OrdinaryDiffEq, StochasticDiffEq, Plots -u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0] -tspan = (0., 200.) -ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1] -ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1] - -oprob1 = ODEProblem(repressilator, u0, tspan, ps1) -oprob2 = ODEProblem(repressilator, u0, tspan, ps2) -osol1 = solve(oprob1, Tsit5()) -osol2 = solve(oprob2, Tsit5()) -oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)") -oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)") - -sprob1 = SDEProblem(repressilator, u0, tspan, ps1) -sprob2 = SDEProblem(repressilator, u0, tspan, ps2) -ssol1 = solve(sprob1, ImplicitEM()) -ssol2 = solve(sprob2, ImplicitEM()) -ssol1 = solve(sprob1, ImplicitEM(); seed = 1) # hide -ssol2 = solve(sprob2, ImplicitEM(); seed = 100) # hide -splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)") -splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)") - -plot(oplt1, oplt2, splt1, splt2; layout = (2,2), size = (800,600)) -``` - -## The Willamowski–Rössler model -The Willamowski–Rössler model was introduced in [*Willamowski & Rössler (1979)*](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) as an example of a simple CRN model which exhibits [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory). This means that small changes in initial conditions can produce relatively large changes in the system's trajectory. -```@example crn_library_chaos -using Catalyst -wr_model = @reaction_network begin - k1, 2X --> 3X - k2, X --> 2X - k3, Z + 2X --> 2Z - k4, Y + X --> 2Y - k5, Y --> ∅ - k6, 2Z --> ∅ - k7, Z --> ∅ -end -``` -Here we first simulate the model for a single initial conditions, showing in both time-state space and phase space how how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/). -```@example crn_library_chaos -using OrdinaryDiffEq, Plots -u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5] -tspan = (0.0, 50.0) -p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7] -oprob = ODEProblem(wr_model, u0, tspan, p) -sol = solve(oprob, Rodas5P()) - -plt1 = plot(sol; title = "Time-state space") -plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space") -plot(plt1, plt2; layout = (1,2), size = (800,400)) -``` \ No newline at end of file diff --git a/docs/src/model_creation/model_visualisation.md b/docs/src/model_creation/model_visualisation.md new file mode 100644 index 0000000000..463c96ce5c --- /dev/null +++ b/docs/src/model_creation/model_visualisation.md @@ -0,0 +1,71 @@ +# [Model Visualisation](@id visualisation) +Catalyst-created `ReactionSystem` models can be visualised either as LaTeX code (of either the model reactions or its equations) or as a network graph. This section describes both functionalities. + +## [Displaying models using LaTeX](@id visualisation_latex) +Once a model has been created, the [Latexify.jl](https://github.com/korsbo/Latexify.jl) package can be used to generate LaTeX code of the model. This can either be used for easy model inspection (e.g. to check which equations are being simulated), or to generate code which can be directly pasted into a LaTeX document. + +Let us consider a simple [Brusselator model](@ref ref): +```@example visualisation_latex +using Catalyst +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end +``` +To display its reaction (using LaTeX formatting) we run `latexify` with our model as input: +```@example visualisation_latex +using Latexify +latexify(brusselator) +``` +Here, we note that the output of `latexify(brusselator)` is identical to how a model is displayed by default. Indeed, the reason is that Catalyst internally uses Latexify's `latexify` function to display its models. It is also possible to display the ODE equations a model would generate by adding the `form = :ode` argument: +```@example visualisation_latex +latexify(brusselator; form = :ode) +``` +!!! note + Internally, `latexify(brusselator; form = :ode)` calls `latexify(convert(ODESystem, brusselator))`. Hence, if you have already [generated the `ODESystem` corresponding to your model](@ref ref), it can be used directly as input to `latexify`. + +!!! note + It should be possible to also generate SDEs through the `form = :sde` input. This feature is, however, currently broken. + +If you wish to copy the output to your [clipboard](https://en.wikipedia.org/wiki/Clipboard_(computing)) (e.g. so that you can paste it into a LaTeX document), run `copy_to_clipboard(true)` before you run `latexify`. A more throughout description of Latexify's features can be found in [its documentation](https://korsbo.github.io/Latexify.jl/stable/). + +!!! note + For a model to be nicely displayed you have to use an IDE that actually supports this (such as a [notebook](https://jupyter.org/)). Other environments (such as [the Julia REPL]([@ref ref](https://docs.julialang.org/en/v1/stdlib/REPL/))) will simply return the full LaTeX code which would generate the desired expression. + +## [Displaying model networks](@id visualisation_graphs) +A network graph showing a Catalyst model's species and reactions can be displayed using the `Graph` function. This first requires [Graphviz](https://graphviz.org/) to be installed and command line accessible. Here, we first declare a [Brusselator model](@ref ref) and then displays its network topology: +```@example visualisation_graphs +using Catalyst +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end +Graph(brusselator) +``` +The network graph represents species as blue nodes and reactions as orange dots. Black arrows from species to reactions indicate substrates, and are labelled with their respective stoichiometries. Similarly, black arrows from reactions to species indicate products (also labelled with their respective stoichiometries). If there are any reactions where a species affect the rate, but does not participate as a reactant, this is displayed with a dashed red arrow. This can be seen in the following [repressilator model](@ref ref): +```@example visualisation_graphs +repressilator = @reaction_network begin + hillr(Z,v,K,n), ∅ --> X + hillr(X,v,K,n), ∅ --> Y + hillr(Y,v,K,n), ∅ --> Z + d, (X, Y, Z) --> ∅ +end +Graph(repressilator) +``` + +A generated graph can be saved using the `savegraph` function: +```@example visualisation_graphs +repressilator_graph = Graph(repressilator) +savegraph(repressilator_graph, "repressilator_graph.png") +rm("repressilator_graph.png") # hide +``` + +Finally, a [network's reaction complexes](@ref ref) (and the reactions in between these) can be displayed using the `complexgraph(brusselator)` function: +```@example visualisation_graphs +complexgraph(brusselator) +``` +Here, reaction complexes are displayed as blue nodes, and reactions in between these as black arrows. \ No newline at end of file diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md new file mode 100644 index 0000000000..78c0dfe8a6 --- /dev/null +++ b/docs/src/model_simulation/ode_simulation_performance.md @@ -0,0 +1,338 @@ +# [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. + +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. + +Generally, this short checklist provides a quick guide for dealing with ODE performance: +1. If performance is not critical, use [the default solver choice](@ref ode_simulation_performance_solvers) and do not worry further about the issue. +2. If improved performance would be useful, read about solver selection (both in [this tutorial](@ref ode_simulation_performance_solvers) and [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)) and then try a few different solvers to find one with good performance. +3. If you have a large ODE (approximately 100 variables or more), try the [various options for efficient Jacobian computation](@ref ode_simulation_performance_jacobian) (noting that some are non-trivial to use, and should only be investigated if truly required). +4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref investigating) (with preference for the former, which is easier to use). + +## [Regarding stiff and non-stiff problems and solvers](@id ode_simulation_performance_stiffness) +Generally, ODE problems can be categorised into [*stiff ODEs* and *non-stiff ODEs*](https://en.wikipedia.org/wiki/Stiff_equation). This categorisation is important due to stiff ODEs requiring specialised solvers. A common cause of failure to simulate an ODE is the use of a non-stiff solver for a stiff problem. There is no exact way to determine whether a given ODE is stiff or not, however, systems with several different time scales (e.g. a CRN with both slow and fast reactions) typically generate stiff ODEs. + +Here we simulate the (stiff) [Brusselator](@ref ref) model using the `Tsit5` solver (which is designed for non-stiff ODEs): +```@example ode_simulation_performance_1 +using Catalyst, OrdinaryDiffEq, Plots + +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end + +u0 = [:X => 1.0, :Y => 0.0] +tspan = (0.0, 20.0) +ps = [:A => 10.0, :B => 40.0] +oprob = ODEProblem(brusselator, u0, tspan, ps) + +sol1 = solve(oprob, Tsit5()) +plot(sol1) +``` +We note that we get a warning, indicating that an instability was detected (the typical indication of a non-stiff solver being used for a stiff ODE). Furthermore, the resulting plot ends at $t ≈ 10$, meaning that the simulation was not completed (as the simulation's endpoint is $t = 20$). Indeed, we can confirm this by checking the *return code* of the solution object: +```@example ode_simulation_performance_1 +sol1.retcode +``` +Next, we instead try the `Rodas5P` solver (which is designed for stiff problems): +```@example ode_simulation_performance_1 +sol2 = solve(oprob, Rodas5P()) +plot(sol2) +``` +This time the simulation was successfully completed, which can be confirmed by checking the return code: +```@example ode_simulation_performance_1 +sol2.retcode +``` + +Generally, ODE solvers can be divided into [*explicit* and *implicit* solvers](https://en.wikipedia.org/wiki/Explicit_and_implicit_methods). Roughly, explicit solvers are better for non-stiff problems, with implicit solvers being required for stiff problems. While we could use implicit solvers for all problems (to guarantee successful simulations irrespective of stiffness), these are typically less performant (as compared to the explicit solvers) on equations that are non-stiff. An important property of implicit solvers is that they require the *computation of a [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant)* as part of their routine. This means that the various options for efficient Jacobian computation [described later in this tutorial](@ref ode_simulation_performance_jacobian) are only relevant to implicit solvers. + +Finally, we should note that stiffness is not tied to the model equations only. If we change the parameter values of our previous Brusselator model to `[:A => 1.0, :B => 4.0]`, the non-stiff `Tsit5` solver can successfully simulate it. + + +## [ODE solver selection](@id ode_simulation_performance_solvers) +OrdinaryDiffEq implements an unusually large number of ODE solvers, with the performance of the simulation heavily depending on which one is chosen. These are provided as the second argument to the `solve` command, e.g. here we use the `Tsit5` solver to simulate a simple [birth-death process](@ref ref): +```@example ode_simulation_performance_2 +using Catalyst, OrdinaryDiffEq + +bd_model = @reaction_network begin + (p,d), 0 <--> X +end +u0 = [:X => 0.1] +tspan = (0.0, 1.0) +ps = [:p => 1.0, :d => 0.2] + +oprob = ODEProblem(bd_model, u0, tspan, ps) +solve(oprob, Tsit5()) +nothing # hide +``` +If no solver argument is provided to `solve`, one is automatically selected: +```@example ode_simulation_performance_2 +solve(oprob) +nothing # hide +``` +While the default choice is typically enough for most single simulations, if performance is important, it can be worthwhile exploring the available solvers to find one that is especially suited for the given problem. A complete list of possible ODE solvers, with advice on optimal selection, can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/). This section will give some general advice. + +The most important part of solver selection is to select one appropriate for [the problem's stiffness](@ref ode_simulation_performance_stiffness). Generally, the `Tsit5` solver is good for non-stiff problems, and `Rodas5P` for stiff problems. For large stiff problems (with many species), `FBDF` can be a good choice. We can illustrate the impact of these choices by simulating our birth-death process using the `Tsit5`, `Vern7` (an explicit solver yielding [low error in the solution](@ref ode_simulation_performance_error)), `Rodas5P`, and `FBDF` solvers (benchmarking their respective performance using [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)): +```@example ode_simulation_performance_2 +using BenchmarkTools +@btime solve(oprob, Tsit5()) +@btime solve(oprob, Vern7()) +@btime solve(oprob, Rodas5P()) +@btime solve(oprob, FBDF()) +``` +Here, we note that the fastest solver is several times faster than the slowest one (`FBDF`, which is a poor choice for this ODE), + +### [Simulation error, tolerance, and solver selection](@id ode_simulation_performance_error) +Numerical ODE simulations [approximate an ODEs' continuous solutions as discrete vectors](https://en.wikipedia.org/wiki/Discrete_time_and_continuous_time). This introduces errors in the computed solution. The magnitude of these errors can be controlled by setting solver *tolerances*. By reducing the tolerance, solution errors will be reduced, however, this will also increase simulation run times. The (absolute and relative) tolerance of a solver can be tuned through the `abstol` and `reltol` arguments. Here we see how run time increases with larger tolerances: +```@example ode_simulation_performance_2 +@btime solve(oprob, Tsit5(); abstol=1e-6, reltol=1e-6) +@btime solve(oprob, Tsit5(); abstol=1e-12, reltol=1e-12) +``` +It should be noted, however, that the result of the second simulation is a lot more accurate. Thus, ODE solver performance cannot be determined solely from run time, but is a composite of run +time and error. Benchmarks comparing solver performance (by plotting the run time as a function of the error) for various CRN models can be found in the [SciMLBenchmarks repository](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Bio/BCR/). + +Generally, whether solution error is a consideration depends on the application. If you want to compute the trajectory of an expensive space probe as it is sent from Earth, to slingshot Jupiter, and then reach Pluto a few years later, ensuring a minimal error will be essential. However, if you want to simulate a simple CRN to determine whether it oscillates for a given parameter set, a small error will not constitute a problem. An important aspect with regard to error is that it affects the selection of the optimal solver. E.g. if tolerance is low (generating larger errors) the `Rosenbrock23` method performs well for small, stiff, problems (again, more details can be found in [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)). + + +## [Jacobian computation options for implicit solvers](@id ode_simulation_performance_jacobian) +As [previously mentioned](@ref ode_simulation_performance_stiffness), implicit ODE solvers require the computation of the system's [*Jacobian*](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant). The reason is (roughly) that, in each time step, these solvers need to solve a non-linear equation to find the simulation's value at the next timestep (unlike explicit solvers, which compute the value at the next time step directly). Typically this is done using [the Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method), which requires the Jacobian. Especially for large systems, this can be computationally expensive (and a potential strain on available memory), in which case one might consider various Jacobian-computation options (as described below). A throughout tutorial on simulating a large, stiff, ODE can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#stiff). + +### [Building the Jacobian symbolically](@id ode_simulation_performance_symbolic_jacobian) +By default, OrdinaryDiffEq computes the Jacobian using [*automatic differentiation*](https://en.wikipedia.org/wiki/Automatic_differentiation) (however, using [*finite differences*](https://en.wikipedia.org/wiki/Finite_difference) is [also possible](https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/)). Since Catalyst builds its ODEs symbolically, it is able to *compute an analytic Jacobian symbolically*. Typically, this is only advantageous when you are also [using a sparse Jacobian](@ref ode_simulation_performance_sparse_jacobian). + +To use this option, simply set `jac = true` when constructing an `ODEProblem`: +```@example ode_simulation_performance_3 +using Catalyst, OrdinaryDiffEq + +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end +u0 = [:X => 1.0, :Y => 0.0] +tspan = (0.0, 20.0) +ps = [:A => 10.0, :B => 40.0] + +oprob = ODEProblem(brusselator, u0, tspan, ps; jac = true) +nothing # hide +``` + +### [Using a sparse Jacobian](@id ode_simulation_performance_sparse_jacobian) +For a system with $n$ variables, the Jacobian will be an $n\times n$ matrix. This means that, as $n$ becomes large, the Jacobian can become *very* large, potentially causing a significant strain on memory. In these cases, most Jacobian entries are typically $0$. This means that a [*sparse*](https://en.wikipedia.org/wiki/Sparse_matrix) Jacobian (rather than a *dense* one, which is the default) can be advantageous. To designate sparse Jacobian usage, simply provide the `sparse = true` option when constructing an `ODEProblem`: +```@example ode_simulation_performance_3 +oprob = ODEProblem(brusselator, u0, tspan, p; sparse = true) +nothing # hide +``` + +### [Linear solver selection](@id ode_simulation_performance_symbolic_jacobian_linear_solver) +When implicit solvers use e.g. the Newton-Raphson method to (at each simulation time step) solve a (typically non-linear) equation, they actually solve a linearised version of this equation. For this, they use a linear solver, the choice of which can impact performance. To specify one, we use the `linsolve` option (given to the solver function, *not* the `solve` command). E.g. to use the `KLUFactorization` linear solver (which requires loading the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) package) we run +```@example ode_simulation_performance_3 +using LinearSolve +solve(oprob, Rodas5P(linsolve = KLUFactorization())) +nothing # hide +``` +A full list of potential linear solvers can be found [here](https://docs.sciml.ai/LinearSolve/dev/solvers/solvers/#Full-List-of-Methods). Typically, the default choice performs well. + +A unique approach to the linear solvers is to use a *matrix-free Newton-Krylov method*. These do not actually compute the Jacobian, but rather *the effect of multiplying it with a vector*. They are typically advantageous for large systems (with large Jacobians), and can be designated using the `KrylovJL_GMRES` linear solver: +```@example ode_simulation_performance_3 +solve(oprob, Rodas5P(linsolve = KrylovJL_GMRES())) +nothing # hide +``` +Since these methods do not depend on a Jacobian, certain Jacobian options (such as [computing it symbolically](@ref ode_simulation_performance_symbolic_jacobian)) are irrelevant to them. + +### [Designating preconditioners for Jacobian-free linear solvers](@ref ode_simulation_performance_preconditioners) +When an implicit method solves a linear equation through an (iterative) matrix-free Newton-Krylov method, the rate of convergence depends on the numerical properties of the matrix defining the linear system. To speed up convergence, a [*preconditioner*](https://en.wikipedia.org/wiki/Preconditioner) can be applied to both sides of the linear equation, attempting to create an equation that converges faster. Preconditioners are only relevant when using matrix-free Newton-Krylov methods. + +In practice, preconditioners are implemented as functions with a specific set of arguments. How to implement these is non-trivial, and we recommend reading OrdinaryDiffEq's documentation pages [here](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#Preconditioners:-precs-Specification) and [here](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Adding-a-Preconditioner). In this example, we will define an [Incomplete LU](https://en.wikipedia.org/wiki/Incomplete_LU_factorization) preconditioner (which requires the [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) package): +```@example ode_simulation_performance_3 +using IncompleteLU +function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata) + if newW === nothing || newW + Pl = ilu(convert(AbstractMatrix, W), τ = 50.0) + else + Pl = Plprev + end + Pl, nothing +end +nothing # hide +``` +Next, `incompletelu` can be supplied to our solver using the `precs` argument: +```@example ode_simulation_performance_3 +solve(oprob, Rodas5P(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true)) +nothing # hide +``` +Finally, we note that when using preconditioners with a matrix-free method (like `KrylovJL_GMRES`, which is also the only case when these are relevant), the `concrete_jac = true` argument is required. + +Generally, the use of preconditioners is only recommended for advanced users who are familiar with the concepts. However, for large systems, if performance is essential, they can be worth looking into. + + +## [Parallelisation on CPUs and GPUs](@id ode_simulation_performance_parallelisation) +Whenever an ODE is simulated a large number of times (e.g. when investigating its behaviour for different parameter values), the best way to improve performance is to [parallelise the simulation over multiple processing units](https://en.wikipedia.org/wiki/Parallel_computing). Indeed, an advantage of the Julia programming language is that it was designed after the advent of parallel computing, making it well-suited for this task. Roughly, parallelisation can be divided into parallelisation on [CPUs](https://en.wikipedia.org/wiki/Central_processing_unit) and on [GPUs](https://en.wikipedia.org/wiki/General-purpose_computing_on_graphics_processing_units). CPU parallelisation is most straightforward, while GPU parallelisation requires specialised ODE solvers (which Catalyst have access to). + +Both CPU and GPU parallelisation require first building an `EnsembleProblem` (which defines the simulations you wish to perform) and then supplying this with the correct parallelisation options. These have [previously been introduced in Catalyst's documentation](@ref ref) (but in the context of convenient bundling of similar simulations, rather than to improve performance), with a more throughout description being found in [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#ensemble). Finally, a general documentation of parallel computing in Julia is available [here](https://docs.julialang.org/en/v1/manual/parallel-computing/). + +### [CPU parallelisation](@id ode_simulation_performance_parallelisation_CPU) +For this example (and the one for GPUs), we will consider a modified [Michaelis-Menten enzyme kinetics model](@ref ref), which describes an enzyme ($E$) that converts a substrate ($S$) to a product ($P$): +```@example ode_simulation_performance_4 +using Catalyst +mm_model = @reaction_network begin + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E + d, S --> ∅ +end +``` +The model can be simulated, showing how $P$ is produced from $S$: +```@example ode_simulation_performance_3 +using OrdinaryDiffEq, Plots +u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] +tspan = (0.0, 50.0) +p = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :d => 0.1] +oprob = ODEProblem(mm_model, u0, tspan, p) +sol = solve(oprob, Tsit5()) +plot(sol) +``` +Due to the degradation of $S$, if the production rate is not high enough, the total amount of $P$ produced is reduced. For this tutorial, we will investigate this effect for a range of values of $kP$. This will require a large number of simulations (for various $kP$ values), which we will parallelise on CPUs (this section) and GPUs ([next section](@ref ode_simulation_performance_parallelisation_GPU)). + +To parallelise our simulations, we first need to create an `EnsembleProblem`. These describe which simulations we wish to perform. The input to this is: +- The `ODEProblem` corresponds to the model simulation (`SDEProblem` and `JumpProblem`s can also be supplied, enabling the parallelisation of these problem types). +- A function, `prob_func`, describing how to modify the problem for each simulation. If we wish to simulate the same, unmodified problem, in each simulation (primarily relevant for stochastic simulations), this argument is not required. + +Here, `prob_func` takes 3 arguments: +- `prob`: The problem that it modifies at the start of each individual run (which will be the same as `EnsembleProblem`'s first argument). +- `i`: The index of the specific simulation (in the array of all simulations that are performed). +- `repeat`: The repeat of a specific simulation in the array. We will not use this option in this example, however, it is discussed in more detail [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem). + +and output the `ODEProblem` simulated in the i'th simulation. + +Let us assume that we wish to simulate our model 100 times, for $kP = 0.01, 0.02, ..., 0.99, 1.0$. We define our `prob_func` using [`remake`](@ref ref): +```@example ode_simulation_performance_4 +function prob_func(prob, i, repeat) + return remake(prob; p = [:kP => 0.01*i]) +end +nothing # hide +``` +Next, we can create our `EnsembleProblem`: +```@example ode_simulation_performance_4 +eprob = EnsembleProblem(oprob; prob_func) +nothing # hide +``` +We can now solve our `ODEProblem` using the same syntax we would normally use to solve the original `ODEProblem`, with the exception that an additional argument, `trajectories`, is required (which denotes how many simulations should be performed). +```@example ode_simulation_performance_4 +esol = solve(eprob, Tsit5(); trajectories = 100) +nothing # hide +``` +To access the i'th solution we use `esol.u[i]`. To e.g. plot the 47'th solution we use: +```@example ode_simulation_performance_4 +plot(esol.u[47]) +``` +To extract the amount of $P$ produced in each simulation, and plot this against the corresponding $kP$ value, we can use: +```@example ode_simulation_performance_4 +plot(0.01:0.01:1.0, map(sol -> sol[:P][end], esol.u), xguide = "kP", yguide = "P produced", label="") +``` + +Above, we have simply used `EnsembleProblem` as a convenient interface to run a large number of similar simulations. However, these problems have the advantage that they allow the passing of an *ensemble algorithm* to the `solve` command, which describes a strategy for parallelising the simulations. By default, `EnsembleThreads` is used. This parallelises the simulations using [multithreading](https://en.wikipedia.org/wiki/Multithreading_(computer_architecture)) (parallelisation within a single process), which is typically advantageous for small problems on shared memory devices. An alternative is `EnsembleDistributed` which instead parallelises the simulations using [multiprocessing](https://en.wikipedia.org/wiki/Multiprocessing) (parallelisation across multiple processes). To do this, we simply supply this additional solver to the solve command: +```@example ode_simulation_performance_4 +esol = solve(eprob, Tsit5(), EnsembleDistributed(); trajectories=100) +nothing # hide +``` +To utilise multiple processes, you must first give Julia access to these. You can check how many processes are available using the `nprocs` (which requires the [Distributed.jl](https://github.com/JuliaLang/Distributed.jl) package): +```@example ode_simulation_performance_4 +using Distributed +nprocs() +``` +Next, more processes can be added using `addprocs`. E.g. here we add an additional 4 processes: +```@example ode_simulation_performance_4 +addprocs(4) +nothing # hide +``` +Powerful personal computers and HPC clusters typically have a large number of available additional processes that can be added to improve performance. + +While `EnsembleThreads` and `EnsembleDistributed` cover the main cases, additional ensemble algorithms exist. A more throughout description of these can be found [here](https://docs.sciml.ai/DiffEqDocs/dev/features/ensemble/#EnsembleAlgorithms). + +Finally, it should be noted that OrdinaryDiffEq, if additional processes are available, automatically parallelises the [linear solve part of implicit simulations](@ref ode_simulation_performance_symbolic_jacobian_linear_solver). It is thus possible to see performance improvements from adding additional processes when running single simulations, not only multiple parallelised ones (this effect is primarily noticeable for large systems with many species). + +### [GPU parallelisation](@id ode_simulation_performance_parallelisation_GPU) +GPUs are different from CPUs in that they are much more restricted in what computations they can carry out. However, unlike CPUs, they are typically available in far larger numbers. Their original purpose is for rendering graphics (which typically involves solving a large number of very simple computations, something CPUs with their few, but powerful, cores are unsuited for). Recently, they have also started to be applied to other problems, such as ODE simulations. Generally, GPU parallelisation is only worthwhile when you have a very large number of parallel simulations to run (and access to good GPU resources, either locally or on a cluster). + +Generally, we can parallelise `EnsembleProblem`s across several GPUs in a very similar manner to how we parallelised them across several CPUs, but by using a different ensemble algorithm (such as `EnsembleGPUArray`). However, there are some additional requirements: +- GPU parallelisation requires the [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl) package. +- Depending on which GPU hardware is used, a specific back-end package has to be installed and imported (e.g. CUDA for NVIDIA's GPUs or Metal for Apple's). +- For some cases, we must use a special ODE solver supporting simulations on GPUs. + +Furthermore (while not required) to receive good performance, we should also make the following adaptations: +- By default, Julia's decimal numbers are implemented as `Float64`s, however, using `Float32`s is advantageous on GPUs. Ideally, all initial conditions and parameter values should be specified using these. +- We should designate all our vectors (i.e. initial conditions and parameter values) as [static vectors](https://github.com/JuliaArrays/StaticArrays.jl). + +We will assume that we are using the CUDA GPU hardware, so we will first load the [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl) backend package, as well as DiffEqGPU: +```@example ode_simulation_performance_5 +using CUDA, DiffEqGPU +``` +Which backend package you should use depends on your available hardware, with the alternatives being listed [here](https://docs.sciml.ai/DiffEqGPU/stable/manual/backends/). + +Next, we declare our model and `ODEProblem`. However, we make all values `Float64` (by appending `f0` to them) and all vectors static (by adding `@SVector` before their declaration, something which requires the [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) package). +```@example ode_simulation_performance_5 +using Catalyst, OrdinaryDiffEq, StaticArrays + +mm_model = @reaction_network begin + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E + d, S --> ∅ +end + +using OrdinaryDiffEq, Plots +u0 = @SVector [:S => 1.0f0, :E => 1.0f0, :SE => 0.0f0, :P => 0.0f0] +tspan = (0.0f0, 50.0f0) +p = @SVector [:kB => 1.0f0, :kD => 0.1f0, :kP => 0.5f0, :d => 0.1f0] +oprob = ODEProblem(mm_model, u0, tspan, p) +nothing # hide +``` +When we declare our `prob_func` and `EnsembleProblem` we need to ensure that the updated `ODEProblem` uses `Float32`: +```@example ode_simulation_performance_5 +function prob_func(prob, i, repeat) + return remake(prob; p = [:kP => 0.0001f0*i]) +end +eprob = EnsembleProblem(oprob; prob_func = prob_func) +nothing # hide +``` +Here have we increased the number of simulations to 10,000, since this is a more appropriate number for GPU parallelisation (as compared to the 100 simulations we performed in our CPU example). + +We can now simulate our model using a GPU-based ensemble algorithm. Currently, two such algorithms are available, `EnsembleGPUArray` and `EnsembleGPUKernel`. Their differences are that +- Only `EnsembleGPUKernel` requires arrays to be static arrays (although it is still advantageous for `EnsembleGPUArray`). +- While `EnsembleGPUArray` can use standard ODE solvers, `EnsembleGPUKernel` requires specialised versions (such as `GPUTsit5`). A list of available such solvers can be found [here](https://docs.sciml.ai/DiffEqGPU/dev/manual/ensemblegpukernel/#specialsolvers). + +Generally, it is recommended to use `EnsembleGPUArray` for large models (that have at least $100$ variables), and `EnsembleGPUKernel` for smaller ones. Here we simulate our model using both approaches (noting that `EnsembleGPUKernel` requires `GPUTsit5`): +```@example ode_simulation_performance_5 +esol1 = solve(eprob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()); trajectories = 10000) +esol2 = solve(eprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()); trajectories = 10000) +nothing # hide +``` +Note that we have to provide the `CUDA.CUDABackend()` argument to our ensemble algorithms (to designate our GPU backend, in this case, CUDA). + +Just like OrdinaryDiffEq is able to utilise parallel CPU processes to speed up the linear solve part of ODE simulations, GPUs can also be used. More details on this can be found [here](https://docs.sciml.ai/DiffEqGPU/stable/tutorials/within_method_gpu/). This is only recommended when ODEs are very large (at least 1,000 species), which is typically not the case for CRNs. + +For more information on differential equation simulations on GPUs in Julia, please read [DiffEqGPU's documentation](https://docs.sciml.ai/DiffEqGPU/stable/). Furthermore, if performance is critical, [this tutorial](https://docs.sciml.ai/DiffEqGPU/stable/tutorials/lower_level_api/) provides information on how to redesign your simulation code to make it more suitable for GPU simulations. + + +--- +## Citation +If you use GPU simulations in your research, please cite the following paper to support the authors of the DiffEqGPU package: +``` +@article{utkarsh2024automated, + title={Automated translation and accelerated solving of differential equations on multiple GPU platforms}, + author={Utkarsh, Utkarsh and Churavy, Valentin and Ma, Yingbo and Besard, Tim and Srisuma, Prakitr and Gymnich, Tim and Gerlach, Adam R and Edelman, Alan and Barbastathis, George and Braatz, Richard D and others}, + journal={Computer Methods in Applied Mechanics and Engineering}, + volume={419}, + pages={116591}, + year={2024}, + publisher={Elsevier} +} +``` + +--- +## References +[^1]: [E. Hairer, G. Wanner, *Solving Ordinary Differential Equations II*, Springer (1996).](https://link.springer.com/book/10.1007/978-3-642-05221-7) \ No newline at end of file diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md new file mode 100644 index 0000000000..113b0fa7b4 --- /dev/null +++ b/docs/src/model_simulation/simulation_introduction.md @@ -0,0 +1,335 @@ +# [Model Simulation Introduction](@id simulation_intro) +Catalyst's core functionality is the creation of *chemical reaction network* (CRN) models that can be simulated using ODE, SDE, and jump simulations. How such simulations are carried out has already been described in [Catalyst's introduction](@ref introduction_to_catalyst). This page provides a deeper introduction, giving some additional background and introducing various simulation-related options. + +Here we will focus on the basics, with other sections of the simulation documentation describing various specialised features, or giving advice on performance. Anyone who plans on using Catalyst's simulation functionality extensively is recommended to also read the documentation on [solution plotting](@ref ref), and on how to [interact with simulation problems, integrators, and solutions](@ref ref). Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for [ODEs](@ref ref), [SDEs](@ref ref), or [jump simulations](@ref ref). + +### [Background to CRN simulations](@id simulation_intro_theory) +This section provides some brief theory on CRN simulations. For details on how to carry out these simulations in actual code, please skip to the following sections. + +CRNs are defined by a set of *species* (with the amounts of these determining the system's state during simulations) and a set of *reaction events* (rules for how the state of the system changes). In real systems, the species amounts are *discrete copy-numbers*, describing the exact numbers of each species type present in the system (in systems biology this can e.g. be the number of a specific molecule present in a cell). Given rates for these reaction events, *stochastic chemical kinetics* provides a formula for simulating the system that recreates its real reaction process. During stochastic chemical kinetics simulations, the system's state is defined by discrete copy-numbers (denoting the number of each species present in the system). Next, at the occurrence of individual *reaction events*, the system's state is updated according to the occurred reaction. The result is a stochastic process. The most well-known approach for simulating stochastic chemical kinetics is [Gillespie's algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm). + +In practice, these jump simulations are computationally expensive. In many cases, copy-numbers are so large that they can be approximated as *continuous concentrations*, and the time-development of the system as a *deterministic process*. This creates an ordinary differential equation (ODE), and is the chemical reaction network form most people are most familiar with. The rule for how ODEs are generated from CRNs is called the [*reaction rate equation*](https://en.wikipedia.org/wiki/Rate_equation) (RRE). + +Here, the RRE enables fast, approximate, and deterministic simulations of CRNs, while stochastic chemical kinetics enables exact, stochastic, simulations of the true process. An intermediary approach is to use the [*chemical Langevin equation*](https://pubs.aip.org/aip/jcp/article/113/1/297/184125/The-chemical-Langevin-equation) (CLE) to formulate a stochastic differential equation (SDE). This approximates the system's state as continuous concentrations, but *does not* assume that its time development is deterministic. Generally, the CLE is used when copy-numbers are large enough that the continuous approximation holds, but not so large that the system's behaviour is deterministic. Generally, the advantage of SDE simulations (compared to jump ones) is that they are faster. Also, since the system state is continuous, interpretation of e.g. stability and steady state results from the deterministic (also continuous) domain is easier for SDEs (however one *should be careful* when making such interpretations). + +These three different approaches are summed up in the following table: +```@raw html + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
InterpretationReaction rate equationChemical Langevin equationStochastic chemical kinetics
Simulation formODE simulationsSDE simulationsJump simulations
Example simulation methods[Euler](https://en.wikipedia.org/wiki/Euler_method#:~:text=The%20Euler%20method%20is%20a,proportional%20to%20the%20step%20size.), [Runge-Kutta](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods)[Euler-Maruyama](https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method), [Milstein](https://en.wikipedia.org/wiki/Milstein_method)[Gillespie](https://en.wikipedia.org/wiki/Gillespie_algorithm), [Sorting direct](https://pubmed.ncbi.nlm.nih.gov/16321569/)
Species unitsConcentrationConcentrationCopy-numbers
Deterministic/StochasticDeterministicStochasticStochastic
ApplicabilityLarge species amountsNon-small species amountsAny species amounts
SpeedTypically fastTypically intermediateTypically slow
Simulation package[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)[StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl)[JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl)
+``` + +## [Performing (ODE) simulations](@id simulation_intro_ODEs) +The following section gives a (more throughout than [previous]) introduction of how to simulate Catalyst models. This is exemplified using ODE simulations (some ODE-specific options will also be discussed). Later on, we will describe options specific to [SDE](@ref ref) and [jump](@ref ref) simulations. All ODE simulations are performed using the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package, which full documentation can be found [here](https://docs.sciml.ai/OrdinaryDiffEq/stable/). A dedicated section giving advice on how to optimise ODE simulation performance can be found [here](@ref ref) + +To perform any simulation, we must first define our model, as well as the simulation's initial conditions, time span, and parameter values. Here we will use a simple [two-state model](@ref ref): +```@example simulation_intro_ode +using Catalyst +two_state_model = @reaction_network begin + (k1,k2), X1 <--> X2 +end +u0 = [:X1 => 100.0, :X2 => 200.0] +tspan = (0.0, 5.0) +ps = [:k1 => 2.0, :k2 => 5.0] +nothing # hide +``` +To simulate the model we first bundle these up into an `ODEProblem`: +```@example simulation_intro_ode +oprob = ODEProblem(two_state_model, u0, tspan, ps) +nothing # hide +``` +Next, we can simulate the model (requires loading the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package). Simulations are performed using the `solve` function. +```@example simulation_intro_ode +using OrdinaryDiffEq +sol = solve(oprob) +nothing # hide +``` +Finally, the result can be plotted using the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package's `plot` function: +```@example simulation_intro_ode +using Plots +plot(sol) +``` +More information on how to interact with solution structures is provided [here](@ref simulation_structure_interfacing) and on how to plot them [here](@ref ref). + +Some additional considerations: +- If a model without parameters has been declared, only the first three arguments must be provided to `ODEProblem`. +- While the first value of `tspan` will almost always be `0.0`, other starting times (both negative and positive) are possible. +- A discussion of various ways to represent species and parameters when designating their values in the `u0` and `ps` vectors can be found [here](@ref ref). + + +### [Designating solvers and solver options](@id simulation_intro_solver_options) +While good defaults are generally selected, OrdinaryDiffEq enables the user to customise simulations through a long range of options that can be provided to the `solve` function. This includes specifying a [solver algorithm](https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations), which can be provided as a second argument to `solve` (if none is provided, a suitable choice is automatically made). E.g. here we specify that the `Rodas5P` method should be used: +```@example simulation_intro_ode +sol = solve(oprob, Rodas5P()) +nothing # hide +``` +A full list of available solvers is provided [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/), and a discussion on optimal solver choices [here](@ref ref). + +Additional options can be provided as keyword arguments. E.g. the `adaptive` arguments determine whether adaptive time-stepping is used (for algorithms that permit this). This defaults to `true`, but can be disabled using +```@example simulation_intro_ode +sol = solve(oprob; adaptive = false) +nothing # hide +``` + +Here follows a list of solver options which might be of interest to the user. +- `adaptive`: Toggles adaptive time stepping for valid methods. Default to `true`. +- `dt`: For non-adaptive simulations, sets the step size (also sets the initial step size for adaptive methods). +- `saveat`: Determines the time points at which the simulation is saved. E.g. for `saveat = 2.0` the simulation is saved every second time unit. If not given, the solution is saved after each time step. +- `save_idxs`: Provides a vector of species whose values should be saved during the simulation. E.g. for `save_idxs = [:X1]`, only the value of species $X1$ is saved. +- `maxiters`: The maximum number of time steps of the simulation. If this number is reached, the simulation is terminated. +- `seed`: Sets a seed for stochastic simulations. Stochastic simulations with the same seed generate identical results. + +A full list of solver options can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). + +### [Alternative problem input forms](@id simulation_intro_ODEs_input_forms) +Throughout Catalyst's documentation, we typically provide initial condition and parameter values as vectors. However, these can also be provided as tuples: +```@example simulation_intro_ode +u0 = (:X1 => 100.0, :X2 => 200.0) +tspan = (0.0, 5.0) +ps = (:k1 => 2.0, :k2 => 5.0) +oprob = ODEProblem(two_state_model, u0, tspan, ps) +nothing # hide +``` +or dictionaries: +```@example simulation_intro_ode +u0 = Dict([:X1 => 100.0, :X2 => 200.0]) +tspan = (0.0, 5.0) +ps = Dict([:k1 => 2.0, :k2 => 5.0]) +oprob = ODEProblem(two_state_model, u0, tspan, ps) +nothing # hide +``` +The forms used for `u0` and `ps` does not need to be the same (but can e.g. be a vector and a tuple). + +!!! note + It [is possible](@ref ref) to designate specific types for parameters. When this is done, the tuple form for providing parameter values should be preferred. + +Throughout Catalyst's documentation, we typically provide the time span as a tuple. However, if the first time point is `0.0` (which is typically the case), this can be omitted. Here, we supply only the simulation endpoint to our `ODEProblem`: +```@example simulation_intro_ode +tend = 5.0 +oprob = ODEProblem(two_state_model, u0, tend, ps) +nothing # hide +``` + +## [Performing SDE simulations](@id simulation_intro_SDEs) +Catalyst uses the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package to perform SDE simulations. This section provides a brief introduction, with [StochasticDiffEq's documentation](https://docs.sciml.ai/StochasticDiffEq/stable/) providing a more extensive description. A dedicated section giving advice on how to optimise SDE simulation performance can be found [here](@ref ref). By default, Catalyst generates SDEs from CRN models using the chemical Langevin equation. + +SDE simulations are performed in a similar manner to ODE simulations. The only exception is that an `SDEProblem` is created (rather than an `ODEProblem`). Furthermore, the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package (rather than the OrdinaryDiffEq package) is required for performing simulations. Here we simulate the two-state model for the same parameter set as previously used: +```@example simulation_intro_sde +using Catalyst, StochasticDiffEq +two_state_model = @reaction_network begin + (k1,k2), X1 <--> X2 +end +u0 = [:X1 => 100.0, :X2 => 200.0] +tspan = (0.0, 1.0) +ps = [:k1 => 2.0, :k2 => 5.0] + +sprob = SDEProblem(two_state_model, u0, tspan, ps) +sol = solve(sprob, STrapezoid()) +sol = solve(sprob, STrapezoid(); seed = 123) # hide +plot(sol) +``` +we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations. + +!!! note + Unlike for ODE and jump simulations, there are no good heuristics for automatically selecting suitable SDE solvers. Hence, for SDE simulations a solver must be provided. `STrapezoid` will work for a large number of cases. When this is not the case, however, please check the list of [available SDE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) for a suitable alternative (making sure to select one compatible with non-diagonal noise and the [Ito interpretation]https://en.wikipedia.org/wiki/It%C3%B4_calculus). + +### [Common SDE simulation pitfalls](@id simulation_intro_SDEs_pitfalls) +Next, let us reduce species amounts (using [`remake`](@ref ref)), thereby also increasing the relative amount of noise, we encounter a problem when the model is simulated: +```@example simulation_intro_sde +sprob = remake(sprob; u0 = [:X1 => 0.33, :X2 => 0.66]) +sol = solve(sprob, STrapezoid()) +sol = solve(sprob, STrapezoid(); seed = 1234567) # hide +plot(sol) +``` +Here, we receive a warning that the simulation was aborted. In the plot, we also see that it is incomplete. In this case we also note that species concentrations are very low (and sometimes, due to the relatively high amount of noise, even negative). This, combined with the early termination, suggests that we are simulating our model for too low species concentration for the assumptions of the CLE to hold. Instead, [jump simulations](@ref simulation_intro_jumps) should be used. + +Next, let us consider a simulation for another parameter set: +```@example simulation_intro_sde +sprob = remake(sprob; u0 = [:X1 => 100.0, :X2 => 200.0], p = [:k1 => 200.0, :k2 => 500.0]) +sol = solve(sprob, STrapezoid()) +sol = solve(sprob, STrapezoid(); seed = 12345) # hide +plot(sol) +``` +Again, the simulation is aborted. This time, however, species concentrations are relatively large, so the CLE might still hold. What has happened this time is that the accuracy of the simulations has not reached its desired threshold. This can be deal with [by reducing simulation tolerances](@ref ref): +```@example simulation_intro_sde +sol = solve(sprob, STrapezoid(), abstol = 1e-1, reltol = 1e-1) +sol = solve(sprob, STrapezoid(); seed = 12345, abstol = 1e-1, reltol = 1e-1) # hide +plot(sol) +``` + +### [Scaling the noise in the chemical Langevin equation](@id simulation_intro_SDEs_noise_saling) +When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to scale the magnitude of the noise. This can be done by introducing a *noise scaling term*, with each noise term generated by the CLE being multiplied with this term. A noise scaling term can be set using the `@default_noise_scaling` option: +```@example simulation_intro_sde +two_state_model = @reaction_network begin + @default_noise_scaling 0.1 + (k1,k2), X1 <--> X2 +end +``` +Here, we set the noise scaling term to `0.1`, reducing the noise with a factor $10$ (noise scaling terms $>1.0$ increase the noise, while terms $<1.0$ reduce the noise). If we re-simulate the model using the low-concentration settings used previously, we see that the noise has been reduced (in fact by so much that the model can now be simulated without issues): +```@example simulation_intro_sde +u0 = [:X1 => 100.0, :X2 => 200.0] +tspan = (0.0, 1.0) +ps = [:k1 => 200.0, :k2 => 500.0] +sprob = SDEProblem(two_state_model, u0, tspan, ps) +sol = solve(sprob, STrapezoid()) +sol = solve(sprob, STrapezoid(); seed = 123) # hide +plot(sol) +``` + +The `@default_noise_scaling` option can take any expression. This can be used to e.g. designate a *noise scaling parameter*: +```@example simulation_intro_sde +two_state_model = @reaction_network begin + @parameters η + @default_noise_scaling η + (k1,k2), X1 <--> X2 +end +``` +Now we can tune the noise through $η$'s value. E.g. here we remove the noise entirely by setting $η = 0.0$ (thereby recreating an ODE simulation's behaviour): +```@example simulation_intro_sde +u0 = [:X1 => 0.33, :X2 => 0.66, :η => 0.0] +tspan = (0.0, 1.0) +ps = [:k1 => 2.0, :k2 => 5.0] +sprob = SDEProblem(two_state_model, u0, tspan, ps) +sol = solve(sprob, STrapezoid()) +sol = solve(sprob, STrapezoid(); seed = 123) # hide +plot(sol) +``` + +!!! note + Above, Catalyst is unable to infer that $η$ is a parameter from the `@default_noise_scaling η` option only. Hence, `@parameters η` is used to explicitly declare $η$ to be a parameter (as discussed in more detail [here](@ref ref)). + +It is possible to designate specific noise scaling terms for individual reactions through the `noise_scaling` [reaction metadata](@ref ref). Here, CLE noise terms associated with a specific reaction are multiplied by that reaction's noise scaling term. Here we use this to turn off the noise in the $X1 \to X2$ reaction: +```@example simulation_intro_sde +two_state_model = @reaction_network begin + k1, X1 <--> X2, [noise_scaling = 0.0] + k2, X2 --> X1 +end +nothing # hide +``` +If the `@default_noise_scaling` option is used, that term is only applied to reactions *without* `noise_scaling` metadata. + +While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref ref), the [`remake_reactionsystem`](@ref) function can be used to achieve a similar effect. + +## [Performing jump simulations using stochastic chemical kinetics](@id simulation_intro_jumps) + +Catalyst uses the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) package to perform jump simulations. This section provides a brief introduction, with [JumpProcesses's documentation](https://docs.sciml.ai/JumpProcesses/stable/) providing a more extensive description. A dedicated section giving advice on how to optimise jump simulation performance can be found [here](@ref ref). + +Jump simulations are performed using so-called `JumpProblem`s. Unlike ODEs and SDEs (for which the corresponding problem types can be created directly), jump simulations require first creating an intermediary `DiscreteProblem`. In this example, we first declare our two-state model and its initial conditions, time span, and parameter values. +```@example simulation_intro_jumps +using Catalyst +two_state_model = @reaction_network begin + (k1,k2), X1 <--> X2 +end +u0 = [:X1 => 5, :X2 => 10] +tspan = (0.0, 5.0) +ps = [:k1 => 2.0, :k2 => 5.0] +nothing # hide +``` +!!! note + Since jump simulations typically simulate the integer copy-numbers of each species present in the system, we designate our initial conditions for jump simulations as integers. Decimal-numbered initial conditions (and thus jump simulations) are, however, also possible. While ODE and SDE simulations accept integer initial conditions, these will be converted to decimal numbers. + +Next, we bundle these into a `DiscreteProblem` (similarly to how `ODEProblem`s and `SDEProblem`s are created): +```@example simulation_intro_jumps +dprob = DiscreteProblem(two_state_model, u0, tspan, ps) +nothing # hide +``` +This is then used as input to a `JumpProblem`. The `JumpProblem` also requires the CRN model as input. +```@example simulation_intro_jumps +jprob = JumpProblem(two_state_model, dprob, Direct()) +nothing # hide +``` +The `JumpProblem` can now be simulated using `solve` (just like any other problem type). +```@example simulation_intro_jumps +using JumpProcesses +sol = solve(jprob, SSAStepper()) +nothing # hide +``` +If we plot the solution we can see how the system's state does not change continuously, but instead in discrete jumps (due to the occurrence of the individual reactions of the system). +```@example simulation_intro_jumps +plot(sol) +``` + +### [Designating aggregators and simulation methods for jump simulations](@id simulation_intro_jumps_solver_designation) +Jump simulations (just like ODEs and SDEs) are performed using solver methods. Unlike ODEs and SDEs, jump simulations are carried out by two different types of methods acting in tandem. First, an *aggregator* method is used to (after each reaction) determine the time to, and type of, the next reaction. Next, a simulation method is used to actually carry out the simulation. + +Several different aggregators are available (a full list is provided [here](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation)). To designate a specific one, provide it as the third argument to the `JumpProblem`. E.g. to designate that Gillespie's direct method (`Direct`) should be used, use: +```@example simulation_intro_jumps +jprob = JumpProblem(brusselator, dprob, Direct()) +nothing # hide +``` +Especially for large systems, the choice of aggregator is relevant to simulation performance. A guide for aggregator selection is provided [here](@ref ref). + +Next, a simulation method can be provided (like for ODEs and SDEs) as the second argument to `solve`. Primarily two alternatives are available, `SSAStepper` and `FunctionMap` (other alternatives are only relevant when jump simulations are combined with ODEs/SDEs, which is described in more detail in JumpProcesses's documentation). Generally, `FunctionMap` is only used when a [continuous callback](@ref ref) is used (and `SSAStepper` otherwise). E.g. we can designate that the `FunctionMap` method should be used through: +```@example simulation_intro_jumps +sol = solve(jprob, FunctionMap()) +nothing # hide +``` + +### [Jump simulations where some rate depends on time](@id simulation_intro_jumps_variableratejumps) +For some models, the rate of some reactions depend on time. E.g. consider the following [circadian model](https://en.wikipedia.org/wiki/Circadian_rhythm), where the production rate of some protein ($P$) depends on a sinusoid function: +```@example simulation_intro_jumps +circadian_model = @reaction_network begin + A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P + d, P --> 0 +end +``` +This type of model will generate so called [*variable rate jumps*](@ref ref). Simulation of such model is non-trivial (and Catalyst currently lacks a good interface for this). A detailed description of how to carry out jump simulations for models with time-dependant rates can be found [here](https://docs.sciml.ai/JumpProcesses/stable/tutorials/simple_poisson_process/#VariableRateJumps-for-processes-that-are-not-constant-between-jumps). \ No newline at end of file diff --git a/docs/src/model_simulation/simulation_plotting.md b/docs/src/model_simulation/simulation_plotting.md new file mode 100644 index 0000000000..25cbec54ee --- /dev/null +++ b/docs/src/model_simulation/simulation_plotting.md @@ -0,0 +1,87 @@ +# [Simulation_plotting](@id simulation_plotting) +Catalyst uses the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package for performing all plots. This section provides a brief summary of some useful plotting options, while [Plots.jl's documentation](https://docs.juliaplots.org/stable/) provides a more throughout description of how to tune your plots. + +!!! note + [Makie.jl](https://github.com/MakieOrg/Makie.jl) is a popular alternative to the Plots.jl package. While it is not used within Catalyst's documentation, it is worth considering (especially for users interested in interactivity, or increased control over their plots). + +## [Common plotting options](@id simulation_plotting_options) +Let us consider the oscillating [Brusselator](@ref ref) model. We have previously shown how model simulation solutions can be plotted using the `plot` function. Here we plot an ODE simulation from the [Brusselator](@ref ref) model: +```@example simulation_plotting +using Catalyst, OrdinaryDiffEq, Plots + +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end +u0 = [:X => 1.0, :Y => 0.0] +tspan = (0.0, 50.0) +ps = [:A => 1.0, :B => 4.0] + +oprob = ODEProblem(brusselator, u0, tspan, ps) +sol = solve(oprob) +plot(sol) +``` + +Various plotting options can be provided as optional arguments to the `plot` command. Common options include: +- `lw`: Determine plot line widths. +- `la`: Determine plot line's transparency (at `la = 0.0` lines are fully transparent, i.e. not visible). +- `linestyle`: Determines plot line style. +- `color`: Determines the line colours. +- `legend`: Determines the position of the legend/labels. +- `label`: Determines label texts. +- `xguide`, `yguide`: Determines x and y axis labels. +- `title`: Determines plot title. +- `legendfontsize`, `guidefontsize`, `titlefontsize`: Determines the font size of the labels, x and y guides, and title, respectively. + +Here, we re-plot our simulations, utilising some of these options (`legend = :none` is used to disable the legends). +```@example simulation_plotting +plot(sol; lw = 4, linestyle = :dash, color = :green, xguide = "Time", yguide = "Concentration", guidefontsize = 14) +``` +Note that, by default, Catalyst uses `xguide = "t"`. Here, however, we modify this to `xguide = "Time"`. We also note that the `color = :green` change both lines' colours to green. To set different colours for each line, we provide these as *a vector without `,` in-between elements* (in Julia interpreted as a matrix with its first dimension equal to `1`): +```@example simulation_plotting +plot(sol; lw = 4, color = [:green :purple]) +``` +A full list of available colours can be found [here](https://juliagraphics.github.io/Colors.jl/stable/namedcolors/). A full list of possible plotting options can be found [here](https://docs.juliaplots.org/stable/attributes/) (look at the list of various plot attributes, e.g. "Series Attributes"). if there is some option(s) you intend to use multiple times, you can call the `default` function using these, in which case they will be used for all subsequent plots. E.g. here: +```@example simulation_plotting +default(framestyle = :box, grid = false) +``` +we designate a box-style frame, and remove the faint background grid, for all subsequent plots in this tutorial. + +A useful option unique to Catalyst (and other DifferentialEquations.jl-based) plots is `idxs`. Its input is a vector, listing all the species (or quantities) that should be plotted. I.e. +```@example simulation_plotting +plot(sol; idxs = [:X]) +``` +can be used to plot `X` only. When only a single argument is given, the vector form is unnecessary (e.g. `idxs = :X` could have been used instead). If [symbolic species representation is used](@ref ref), this can be used to designate any algebraic expression(s) that should be plotted. E.g. here we plot the total concentration of $X + Y$ throughout the simulation: +```@example simulation_plotting +plot(sol; idxs = brusselator.X + brusselator.Y) +``` + +## [Multi-plot plots](@id simulation_plotting_options) +It is possible to save plots in variables. These can then be used as input to the `plot` command. Here, the plot command can be used to create plots containing multiple plots (by providing multiple inputs). E.g. here we plot the concentration of `X` and `Y` in separate subplots: +```@example simulation_plotting +plt_X = plot(sol; idxs = [:X]) +plt_Y = plot(sol; idxs = [:Y]) +plot(plt_X, plt_Y) +``` + +When working with subplots, the [`layout`](https://docs.juliaplots.org/latest/layouts/) and [`size`](https://docs.juliaplots.org/latest/generated/attributes_plot/) options are typically useful. Here we use `layout` to put the first plot above the second one, and `size` to reshape the plot size: +```@example simulation_plotting +plot(plt_X, plt_Y; layout = (2,1), size = (700,500)) +``` + +## [Saving plots](@id simulation_plotting_options) +Once a plot has been saved to a variable, the `savefig` function can be used to save it to a file. Here we save our Brusselator plot simulation (the first argument) to a file called "saved_plot.png" (the second argument): +```@example simulation_plotting +plt = plot(sol) +savefig(plt, "saved_plot.png") +rm("saved_plot.png") # hide +``` +The plot file type is automatically determined from the extension (if none is given, a .png file is created). + +## [Phase-space plots](@id simulation_plotting_options) +By default, simulations are plotted as species concentrations over time. However, [phase space](https://en.wikipedia.org/wiki/Phase_space#:~:text=In%20dynamical%20systems%20theory%20and,point%20in%20the%20phase%20space.) plots are also possible. This is done by designating the axis arguments using the `idxs` option, but providing them as a tuple. E.g. here we plot our simulation in $X-Y$ space: +```@example simulation_plotting +plot(sol; idxs = (:X, :Y)) +``` \ No newline at end of file diff --git a/docs/src/steady_state_functionality/homotopy_continuation.md b/docs/src/steady_state_functionality/homotopy_continuation.md index b5fbe4f850..19a15d39e6 100644 --- a/docs/src/steady_state_functionality/homotopy_continuation.md +++ b/docs/src/steady_state_functionality/homotopy_continuation.md @@ -4,18 +4,18 @@ The steady states of a dynamical system ${dx \over dt} = f(x)$ can be found by solving $0 = f(x)$. This is typically a hard problem, and generally, there is no 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 with +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. This is implemented in Julia through the +*homotopy continuation* algorithm [^1]. This is implemented in Julia through the [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/) package. 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. -## [Basic example](@id homotopy_continuation_basic_example) -For this tutorial, we will use a model from Wilhelm (2009)[^1] (which + +For this tutorial, we will use the [Wilhelm model](@ref ref) (which demonstrates bistability in a small chemical reaction network). We declare the model and the parameter set for which we want to find the steady states: -```@example hc1 +```@example hc_basics using Catalyst import HomotopyContinuation @@ -26,12 +26,12 @@ wilhelm_2009_model = @reaction_network begin k4, X --> 0 end ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5] -nothing # hide +nothing # hide ``` Here, we only run `import HomotopyContinuation` as we do not require any of its functions, and just need it to be present in the current scope for the extension to be activated. Now we can find the steady states using: -```@example hc1 +```@example hc_basics hc_steady_states(wilhelm_2009_model, ps) ``` The order of the species in the output vectors are the same as in `species(wilhelm_2009_model)`. @@ -41,25 +41,25 @@ It should be noted that the steady state multivariate polynomials produced by re ## [Systems with conservation laws](@id homotopy_continuation_conservation_laws) Some systems are under-determined, and have an infinite number of possible steady states. These are typically systems containing a conservation law, e.g. -```@example hc3 -using Catalyst -import HomotopyContinuation - +```@example hc_claws +using Catalyst # hide two_state_model = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -Catalyst allows the conservation laws of such systems to be computed using the `conservationlaws` function. By using these to reduce the dimensionality of the system, as well specifying the initial amount of each species, HomotopyContinuation can again be used to find steady states. To find the steady states using the Catalyst interface to HomotopyContinuation, an initial condition must be provided (which is used to compute the system's conserved quantities, in this case `X1+X2`): -```@example hc3 +Catalyst allows the conservation laws of such systems to be computed [using the `conservationlaws` function](@ref ref). By using these to reduce the dimensionality of the system, as well as specifying the initial amount of each species, homotopy continuation can again be used to find steady states. Here we do this by designating such an initial condition (which is used to compute the system's conserved quantities, in this case $X1 + X2$): +```@example hc_claws +import HomotopyContinuation # hide ps = [:k1 => 2.0, :k2 => 1.0] u0 = [:X1 => 1.0, :X2 => 1.0] hc_steady_states(two_state_model, ps; u0) ``` ## Final notes -- `hc_steady_states` supports any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients). +- `hc_steady_states` support any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients). - When providing initial conditions to compute conservation laws, values are only required for those species that are part of conserved quantities. If this set of species is unknown, it is recommended to provide initial conditions for all species. -- Additional arguments provided to `hc_steady_states` are automatically passed to HomotopyContinuation's `solve` command. Use e.g. `show_progress=false` to disable the progress bar. +- Additional arguments provided to `hc_steady_states` are automatically passed to HomotopyContinuation's `solve` command. Use e.g. `show_progress = false` to disable the progress bar. + --- ## [Citation](@id homotopy_continuation_citation) @@ -75,9 +75,9 @@ If you use this functionality in your research, please cite the following paper } ``` + --- ## References -[^1]: [Thomas Wilhelm, *The smallest chemical reaction system with bistability*, BMC Systems Biology (2009).](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) +[^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) -[^3]: [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) -[^4]: [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) \ No newline at end of file +[^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) \ No newline at end of file diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index acf20d56c9..601fd51c93 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/docs/src/steady_state_functionality/nonlinear_solve.md @@ -1,16 +1,18 @@ -# [Finding Steady States using NonlinearSolve.jl](@id nonlinear_solve) +# [Finding Steady States using NonlinearSolve.jl and SteadyStateDiffEq.jl](@id steady_state_solving) -We have previously described how `ReactionSystem` steady states can be found through [homotopy continuation](@ref homotopy_continuation). Catalyst also supports the creation of `NonlinearProblems` corresponding to `ODEProblems` with all left-hand side derivatives set to 0. These can be solved using the variety of nonlinear equation-solving algorithms implemented by [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl), with the solutions corresponding to system steady states. Generally, using this approach is advantageous when: -- Only a single steady state solution is sought. -- The nonlinear system produced by the model does not correspond to a multivariate, rational, polynomial (homotopy continuation cannot be applied to these systems). Examples include models with non-integer hill coefficients or stoichiometric constants. +Catalyst `ReactionSystem` models can be converted to ODEs (through [the reaction rate equation](@ref ref)). 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 forward ODE simulation from an initial condition until a steady state has been reached. -However, if all (or multiple) steady states are sought, using homotopy continuation is better. +While these approaches only find a single steady state, they offer two advantages as compared to homotopy continuation: +- They are typically much faster. +- They can find steady states for models that do not produce multivariate, rational, polynomial systems (which is a requirement for homotopy continuation to work). Examples include models with non-integer hill coefficients. -This tutorial describes how to create `NonlinearProblem`s from Catalyst's `ReactionSystemn`s, and how to solve them using NonlinearSolve. More extensive descriptions of available solvers and options can be found in [NonlinearSolve's documentation](https://docs.sciml.ai/NonlinearSolve/stable/). +In practice, model steady states are found through [nonlinear system solving](@ref steady_state_solving_nonlinear) by creating a `NonlinearProblem`, and through [forward ODE simulation](@ref ref) by creating a `SteadyStateProblem`. These are then solved through solvers implemented in the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl), package (with the latter approach also requiring the [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl) package). This tutorial describes how to find steady states through these two approaches. More extensive descriptions of available solvers and options can be found in [NonlinearSolve's documentation](https://docs.sciml.ai/NonlinearSolve/stable/). -## Basic example -Let us consider a simple dimerisation network, where a protein ($P$) can exist in a monomer and a dimer form. The protein is produced at a constant rate from its mRNA, which is also produced at a constant rate. -```@example nonlinear_solve1 +## [Steady state finding through nonlinear solving](@ref steady_state_solving_nonlinear) +Let us consider a simple dimerisation system, where a protein ($P$) can exist in a monomer and a dimer form. The protein is produced at a constant rate from its mRNA, which is also produced at a constant rate. +```@example steady_state_solving_nonlinear using Catalyst dimer_production = @reaction_network begin pₘ, 0 --> mRNA @@ -36,68 +38,112 @@ To find its steady states we need to solve: \end{aligned} ``` -To solve this problem, we must first designate our parameter values, and also make an initial guess of the solution. Generally, for problems with a single solution (like this one), most arbitrary guesses will work fine (the exception typically being [systems with conservation laws](@ref nonlinear_solve_conservation_laws)). Using these, we can create the `NonlinearProblem` that we wish to solve. -```@example nonlinear_solve1 +To solve this problem, we must first designate our parameter values, and also make an initial guess of the solution. Generally, for problems with a single solution (like this one), most arbitrary guesses will work fine (the exception typically being [systems with conservation laws](@ref steady_state_solving_nonlinear_conservation_laws)). Using these, we can create the `NonlinearProblem` that we wish to solve. +```@example steady_state_solving_nonlinear p = [:pₘ => 0.5, :pₚ => 2.0, :k₁ => 5.0, :k₂ => 1.0, :d => 1.0] u_guess = [:mRNA => 1.0, :P => 1.0, :P₂ => 1.0] -nl_prob = NonlinearProblem(dimer_production, u_guess, p) +nlprob = NonlinearProblem(dimer_production, u_guess, p) ``` Finally, we can solve it using the `solve` command, returning the steady state solution: -```@example nonlinear_solve1 +```@example steady_state_solving_nonlinear using NonlinearSolve -sol = solve(nl_prob) +sol = solve(nlprob) ``` -NonlinearSolve provides [a wide range of potential solvers](https://docs.sciml.ai/NonlinearSolve/stable/solvers/NonlinearSystemSolvers/). If we wish to designate one, it can be supplied as a second argument to `solve`. Here, we use the Newton Trust Region method, and then check that the solution is equal to the previous one. -```@example nonlinear_solve1 -sol_ntr = solve(nl_prob, TrustRegion()) +Typically, a good default method is automatically selected for any problem. However, NonlinearSolve does provide [a wide range of potential solvers](https://docs.sciml.ai/NonlinearSolve/stable/solvers/NonlinearSystemSolvers/). If we wish to designate one, it can be supplied as a second argument to `solve`. Here, we use the [Newton Trust Region method](https://en.wikipedia.org/wiki/Trust_region), and then check that the solution is equal to the previous one. +```@example steady_state_solving_nonlinear +sol_ntr = solve(nlprob, TrustRegion()) sol ≈ sol_ntr ``` - -## [Finding steady states through ODE simulations](@id nonlinear_solve_ode_simulation_based) -The `NonlinearProblem`s generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. NonlinearSolve supports this through the `DynamicSS` method. To use it, an appropriate ODE solver (and any options you wish it to use) must also be supplied (with [a large number being available](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)). Here, we will use the `Tsit5` ODE solver to find the steady states of our dimerisation system. -```@example nonlinear_solve1 -using OrdinaryDiffEq, SteadyStateDiffEq -solve(nl_prob, DynamicSS(Tsit5())) -``` -Here, we had to import [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) (to use `Tsit5`) and [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl) (to use `DynamicSS`). - -Like previously, the found steady state is unaffected by the initial `u` guess. However, when [finding the steady states of systems with conservation laws](@ref nonlinear_solve_conservation_laws) it is important to select a `u` guess that corresponds to the initial condition for the ODE model for which a steady state is sought. - - -## [Systems with conservation laws](@id nonlinear_solve_conservation_laws) -As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](@ref homotopy_continuation_conservation_laws). E.g. consider the following two-state system: -```@example nonlinear_solve2 + +### [Systems with conservation laws](@id steady_state_solving_nonlinear_conservation_laws) +As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](homotopy_continuation_conservation_laws). E.g. consider the following [two-state system](@ref ref): +```@example steady_state_solving_claws using Catalyst, NonlinearSolve # hide two_state_model = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C=X1+X2$) must be provided. Since these can be computed from system initial conditions (`u0`, i.e. those provided when performing ODE simulations), designating an `u0` is often the best way. There are two ways to do this. First, one can perform [ODE simulation-based steady state finding](@ref nonlinear_solve_ode_simulation_based), using the initial condition as the initial `u` guess. Alternatively, any conserved quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by Catalyst's [conservation law finding and elimination feature](@ref network_analysis_deficiency). +It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C = X1 + X2$) must be provided. Since these can be computed from system initial conditions (`u0`, i.e. those provided when performing ODE simulations), designating an `u0` is often the best way. There are two ways to do this. First, one can perform [forward ODE simulation-based steady state finding](@ref steady_state_solving_simulation), using the initial condition as the initial `u` guess. Alternatively, any conserved quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by Catalyst's [conservation law finding and elimination feature](@ref ref). To eliminate conservation laws we simply provide the `remove_conserved = true` argument to `NonlinearProblem`: -```@example nonlinear_solve2 +```@example steady_state_solving_claws p = [:k1 => 2.0, :k2 => 3.0] u_guess = [:X1 => 3.0, :X2 => 1.0] nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true) nothing # hide ``` -here it is important that the quantities used in `u_guess` correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2= 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the calculated steady-state as well. We can now find the steady states using `solve` like before: -```@example nonlinear_solve2 +here it is important that the quantities used in `u_guess` correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2 = 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the computed steady state as well. We can now find the steady states using `solve` like before: +```@example steady_state_solving_claws sol = solve(nl_prob) ``` We note that the output only provides a single value. The reason is that the actual system solved only contains a single equation (the other being eliminated with the conserved quantity). To find the values of $X1$ and $X2$ we can [directly query the solution object for these species' values, using the species themselves as inputs](@ref simulation_structure_interfacing_solutions): -```@example nonlinear_solve2 -@unpack X1, X2 = two_state_model -sol[X1] +```@example steady_state_solving_claws +sol[[:X1, :X2]] +``` + +## [Finding steady states through ODE simulations](@id steady_state_solving_simulation) +The `NonlinearProblem`s generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. Here we do so for the dimerisation system considered in the previous section. First, we declare our model, initial condition, and parameter values. +```@example steady_state_solving_simulation +dimer_production = @reaction_network begin + pₘ, 0 --> mRNA + pₚ, mRNA --> mRNA + P + (k₁, k₂), 2P <--> P₂ + d, (mRNA, P, P₂) --> 0 +end +p = [:pₘ => 0.5, :pₚ => 2.0, :k₁ => 5.0, :k₂ => 1.0, :d => 1.0] +u0 = [:mRNA => 0.1, :P => 0.0, :P₂ => 0.0] +nothing # hide +``` +Next, we provide these as an input to a `SteadyStateProblem` +```@example steady_state_solving_simulation +ssprob = SteadyStateProblem(dimer_production, u0, p) +``` +Finally, we can find the steady states using the `solver` command (which requires loading the SteadyStateDiffEq package). +```@example steady_state_solving_simulation +using SteadyStateDiffEq +solve(ssprob) +``` +Note that, unlike for nonlinear system solving, `u0` is not just an initial guess of the solution, but the initial conditions from which the steady state simulation is carried out. This means that, for a system with multiple steady states, we can determine the steady states associated with specific initial conditions (which is not possible when the nonlinear solving approach is used). This also permits us to easily [handle the presence of conservation laws](@ref steady_state_solving_nonlinear_conservation_laws). The forward ODE simulation approach (unlike homotopy continuation and nonlinear solving) cannot find unstable steady states. + +The forward ODE solving approach uses the ODE solvers implemented by the [OrdinaryDiffEq.jl](@ref ref) package. If this package is loaded, it is possible to designate a specific solver to use. Any available ODE solver can be used, however, it has to be encapsulated by the `DynamicSS()` function. E.g. here we designate the `Rodas5P` solver: +```@example steady_state_solving_simulation +using OrdinaryDiffEqDiffEq +solve(ssprob, DynamicSS(Rodas5P())) +nothing # hide ``` -```@example nonlinear_solve2 -sol[X2] +Generally, `SteadyStateProblem`s can be solved using the [same options that are available for ODE simulations](@ref ref). E.g. here we designate a specific `dt` step size: +```@example steady_state_solving_simulation +solve(ssprob, DynamicSS(Rodas5P()); dt = 0.01) +nothing # hide ``` +It is possible to use solve `SteadyStateProblem`s using a nonlinear solver, and `NonlinearProblem`s using forward ODE simulation solvers: +```@example steady_state_solving_simulation +using NonlinearSolve +solve(ssprob, TrustRegion()) +nothing # hide +``` +```@example steady_state_solving_simulation +nlprob = NonlinearProblem(dimer_production, u0, p) +solve(nlprob, DynamicSS(Rodas5P())) +nothing # hide +``` +However, especially when the forward ODE simulation approach is used, it is recommended to use the problem type which corresponds to the intended solver. + + --- ## [Citations](@id nonlinear_solve_citation) -If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following papers to support the authors of the NonlinearSolve.jl package: +If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following paper to support the authors of the NonlinearSolve.jl package: ``` -A NonlinearSolve. jl-related publication is in preparation, once it is available, its details will be added here. +@article{pal2024nonlinearsolve, + title={NonlinearSolve. jl: High-Performance and Robust Solvers for Systems of Nonlinear Equations in Julia}, + author={Pal, Avik and Holtorf, Flemming and Larsson, Axel and Loman, Torkel and Schaefer, Frank and Qu, Qingyu and Edelman, Alan and Rackauckas, Chris and others}, + journal={arXiv preprint arXiv:2403.16341}, + year={2024} +} ``` + +--- +## References +[^1]: [J. Nocedal, S. J. Wright, *Numerical Optimization*, Springer (2006).](https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf) diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl index 5a5c4b509f..ab15dd68ae 100644 --- a/test/reactionsystem_core/events.jl +++ b/test/reactionsystem_core/events.jl @@ -243,10 +243,6 @@ let sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5()) sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5()) @test sol_dsl == sol_prog - - sol_dsl = solve(SDEProblem(rn_dsl, u0, tspan, ps), ImplicitEM(); seed = 1234) - sol_prog = solve(SDEProblem(rn_prog, u0, tspan, ps), ImplicitEM(); seed = 1234) - @test sol_dsl == sol_prog end # Checks that misformatted events yields errors in the DSL. @@ -313,6 +309,65 @@ end ### Additional Correctness Tests ### +# Tests that events are properly triggered for SDEs. +# Tests for continuous events, and all three types of discrete events. +let + # Creates model with all types of events. The `e` parameters track whether events are triggered. + rn = @reaction_network begin + @parameters e1=0 e2=0 e3=0 e4=0 + @continuous_events begin + [X ~ 1000.0] => [e1 ~ 1] + end + @discrete_events begin + [1.0] => [e2 ~ 1] + 1.0 => [e3 ~ 1] + (Y > 1000.0) & (e4==0) => [e4 ~ 1] + end + (p,d), 0 <--> X + (p,d), 0 <--> Y + end + + # Simulates the model for conditions where it *definitely* will cross `X = 1000.0` + u0 = [:X => 999.9, :Y => 999.9] + ps = [:p => 10.0, :d => 0.001] + sprob = SDEProblem(rn, u0, (0.0, 2.0), ps) + sol = solve(sprob, ImplicitEM(); seed) + + # Checks that all `e` parameters have been updated properly. + @test sol.ps[:e1] == 1 + @test sol.ps[:e2] == 1 + @test sol.ps[:e3] == 1 + @test sol.ps[:e4] == 1 +end + +# Tests that events are properly triggered for Jump simulations. +# Tests for all three types of discrete events. +let + # Creates model with all types of events. The `e` parameters track whether events are triggered. + rn = @reaction_network begin + @parameters e1=0 e2=0 e3=0 + @discrete_events begin + [1.0] => [e1 ~ 1] + # 1.0 => [e2 ~ 1] + (X > 1000.0) & (e3==0) => [e3 ~ 1] + end + (p,d), 0 <--> X + end + + # Simulates the model for conditions where it *definitely* will cross `X = 1000.0` + u0 = [:X => 999] + ps = [:p => 10.0, :d => 0.001] + dprob = DiscreteProblem(rn, u0, (0.0, 2.0), ps) + jprob = JumpProblem(rn, dprob, Direct(); rng) + sol = solve(jprob, SSAStepper(); seed) + + # Checks that all `e` parameters have been updated properly. + # Note that periodic discrete events are currently broken for jump processes. + @test sol.ps[:e1] == 1 + @test_broken sol.ps[:e2] == 1 + @test sol.ps[:e3] == 1 +end + # Compares simulations using MTK type events with those generated through callbacks. # Jump simulations must be handles differently (since these only accepts discrete callbacks). # Checks for all types of discrete callbacks, and for continuous callbacks. diff --git a/test/reactionsystem_core/higher_order_reactions.jl b/test/reactionsystem_core/higher_order_reactions.jl index f575cd5f55..7922905792 100644 --- a/test/reactionsystem_core/higher_order_reactions.jl +++ b/test/reactionsystem_core/higher_order_reactions.jl @@ -103,7 +103,7 @@ let dprob_alt2 = DiscreteProblem(u0_alt2, (0.0, 100.0), ps_alt2) jprob_alt2 = JumpProblem(dprob_alt2, Direct(), higher_order_network_alt2...; rng = StableRNG(1234)) - # Simualtes the models. + # Simulates the models. sol_base = solve(jprob_base, SSAStepper(); seed, saveat = 1.0) sol_alt1 = solve(jprob_alt1, SSAStepper(); seed, saveat = 1.0) sol_alt2 = solve(jprob_alt2, SSAStepper(); seed, saveat = 1.0) diff --git a/test/reactionsystem_core/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl index 745a9c50af..797a78f239 100644 --- a/test/reactionsystem_core/symbolic_stoichiometry.jl +++ b/test/reactionsystem_core/symbolic_stoichiometry.jl @@ -192,64 +192,68 @@ end # Tests symbolic stoichiometries in simulations. # Tests for decimal numbered symbolic stoichiometries. let + # Declares models. The references models have the `n` parameters so they can use the same + # parameter vectors as the non-reference ones. rs_int = @reaction_network begin - @parameters n1::Int64 n2::Int64 - p, 0 -->X - k, n1*X --> n2*Y - d, Y --> 0 + @parameters n::Int64 + (k1, k2), n*X1 <--> X2 end rs_dec = @reaction_network begin - @parameters n1::Float64 n2::Float64 - p, 0 -->X - k, n1*X --> n2*Y - d, Y --> 0 + @parameters n::Float64 + (k1, k2), n*X1 <--> X2 end rs_ref_int = @reaction_network begin - @parameters n1::Int64 n2::Int64 - p, 0 -->X - k, 3*X --> 4*Y - d, Y --> 0 + @parameters n::Int64 + (k1, k2), 3*X1 <--> X2 end rs_ref_dec = @reaction_network begin - @parameters n1::Float64 n2::Float64 - p, 0 -->X - k, 0.5*X --> 2.5*Y - d, Y --> 0 + @parameters n::Float64 + (k1, k2), 2.5*X1 <--> X2 end - u0 = [:X => 8, :Y => 8] - ps_int = (:p => 2.0, :k => 0.01, :n1 => 3, :n2 => 4, :d => 0.2) - ps_dec = [:p => 2.0, :k => 0.01, :n1 => 0.5, :n2 => 2.5, :d => 0.2] - tspan = (0.0, 10.0) + # Set simulation settings. Initial conditions are design to start, more or less, at + # steady state concentrations. + # Values are selected so that stochastic tests should always pass within the bounds (independent + # of seed). + u0_int = [:X1 => 150, :X2 => 600] + u0_dec = [:X1 => 100, :X2 => 600] + tspan_det = (0.0, 1.0) + tspan_stoch = (0.0, 10000.0) + ps_int = (:k1 => 0.00001, :k2 => 0.01, :n => 3) + ps_dec = (:k1 => 0.00001, :k2 => 0.01, :n => 2.5) # Test ODE simulations with integer coefficients. - oprob_int = ODEProblem(rs_int, u0, tspan, ps_int) - oprob_int_ref = ODEProblem(rs_ref_int, u0, tspan, ps_int) - @test solve(oprob_int, Tsit5()) ≈ solve(oprob_int_ref, Tsit5()) + oprob_int = ODEProblem(rs_int, u0_int, tspan_det, ps_int) + oprob_int_ref = ODEProblem(rs_ref_int, u0_int, tspan_det, ps_int) + @test solve(oprob_int, Tsit5())[:X1][end] ≈ solve(oprob_int_ref, Tsit5())[:X1][end] # Test ODE simulations with decimal coefficients. - oprob_dec = ODEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - oprob_dec_ref = ODEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - @test solve(oprob_dec, Tsit5()) ≈ solve(oprob_dec_ref, Tsit5()) + oprob_dec = ODEProblem(rs_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false) + oprob_dec_ref = ODEProblem(rs_ref_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false) + @test solve(oprob_dec, Tsit5())[:X1][end] ≈ solve(oprob_dec_ref, Tsit5())[:X1][end] # Test SDE simulations with integer coefficients. - sprob_int = SDEProblem(rs_int, u0, tspan, ps_int) - sprob_int_ref = SDEProblem(rs_ref_int, u0, tspan, ps_int) - @test solve(sprob_int, ImplicitEM(); seed) ≈ solve(sprob_int_ref, ImplicitEM(); seed) + sprob_int = SDEProblem(rs_int, u0_int, tspan_stoch, ps_int) + sprob_int_ref = SDEProblem(rs_ref_int, u0_dec, tspan_stoch, ps_int) + ssol_int = solve(sprob_int, ImplicitEM(); seed) + ssol_int_ref = solve(sprob_int_ref, ImplicitEM(); seed) + @test mean(ssol_int[:X1]) ≈ mean(ssol_int_ref[:X1]) atol = 2*1e0 # Test SDE simulations with decimal coefficients. - sprob_dec = SDEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - sprob_dec_ref = SDEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - @test solve(sprob_dec, ImplicitEM(); seed) ≈ solve(sprob_dec_ref, ImplicitEM(); seed) - - # Tests jump simulations with integer coefficients. - dprob_int = DiscreteProblem(rs_int, u0, (0.0, 100000.0), ps_int) - dprob_int_ref = DiscreteProblem(rs_ref_int, u0, (0.0, 100000.0), ps_int) - jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng) - jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng) - sol_int = solve(jprob_int, SSAStepper(); seed) - sol_int_ref = solve(jprob_int_ref, SSAStepper(); seed) - @test mean(sol_int[:Y]) ≈ mean(sol_int_ref[:Y]) atol = 1e-2 rtol = 1e-2 + sprob_dec = SDEProblem(rs_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false) + sprob_dec_ref = SDEProblem(rs_ref_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false) + ssol_dec = solve(sprob_dec, ImplicitEM(); seed) + ssol_dec_ref = solve(sprob_dec_ref, ImplicitEM(); seed) + @test mean(ssol_dec[:X1]) ≈ mean(ssol_dec_ref[:X1]) atol = 2*1e0 + + # Test Jump simulations with integer coefficients. + dprob_int = DiscreteProblem(rs_int, u0_int, tspan_stoch, ps_int) + dprob_int_ref = DiscreteProblem(rs_ref_int, u0_int, tspan_stoch, ps_int) + jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng, save_positions = (false, false)) + jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng, save_positions = (false, false)) + jsol_int = solve(jprob_int, SSAStepper(); seed, saveat = 1.0) + jsol_int_ref = solve(jprob_int_ref, SSAStepper(); seed, saveat = 1.0) + @test mean(jsol_int[:X1]) ≈ mean(jsol_int_ref[:X1]) atol = 1e-2 rtol = 1e-2 end # Check that jump simulations (implemented with and without symbolic stoichiometries) yield simulations diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl index 4b1c935435..5684db4132 100644 --- a/test/simulation_and_solving/simulate_SDEs.jl +++ b/test/simulation_and_solving/simulate_SDEs.jl @@ -128,17 +128,6 @@ let duW = zeros(size(rn_manual.nrp)) rn_manual.g(duW, u0_2, ps_2, 0.0) @test duW ≈ g_eval(rn_catalyst, u0_1, ps_1, 0.0) - - # Compares simulation with identical seed. - # Cannot test the third network (requires very fiddly parameters for CLE to be defined). - if nameof(rn_catalyst) != :rnc9 - sprob_1 = SDEProblem(rn_catalyst, u0_1, (0.0, 1.0), ps_1) - sprob_2 = SDEProblem(rn_manual.f, rn_manual.g, u0_2, (0.0, 1.0), ps_2; noise_rate_prototype = rn_manual.nrp) - seed = seed = rand(rng, 1:100) - sol1 = solve(sprob_1, ImplicitEM(); seed) - sol2 = solve(sprob_2, ImplicitEM(); seed) - @test sol1[u0_sym] ≈ sol2.u - end end end end diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index 6425758b89..9c8c02db8a 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -6,6 +6,7 @@ using Catalyst, JumpProcesses, Statistics, Test # Sets stable rng number. using StableRNGs rng = StableRNG(12345) +seed = rand(rng, 1:100) # Fetch test functions and networks. include("../test_functions.jl") @@ -16,12 +17,17 @@ include("../test_networks.jl") # Compares jump simulations generated through hCatalyst, and manually created systems. let - # Manually declares jumps to compare Catalyst-generated jump simulations to. + # Declares vectors to store information about each network in. Initial conditions are approximately + # at the steady state. catalyst_networks = [] manual_networks = [] u0_syms = [] ps_syms = [] + u0s = [] + ps = [] + sps = [] + # Manual declaration of the first network (and its simulation settings). rate_1_1(u, p, t) = p[1] rate_1_2(u, p, t) = p[2] * u[1] rate_1_3(u, p, t) = p[3] * u[2] @@ -47,12 +53,16 @@ let jump_1_7 = ConstantRateJump(rate_1_7, affect_1_7!) jump_1_8 = ConstantRateJump(rate_1_8, affect_1_8!) jumps_1 = (jump_1_1, jump_1_2, jump_1_3, jump_1_4, jump_1_5, jump_1_6, jump_1_7, - jump_1_8) + jump_1_8) push!(catalyst_networks, reaction_networks_standard[5]) push!(manual_networks, jumps_1) push!(u0_syms, [:X1, :X2, :X3, :X4]) push!(ps_syms, [:p, :k1, :k2, :k3, :k4, :k5, :k6, :d]) + push!(u0s, [:X1 => 40, :X2 => 30, :X3 => 20, :X4 => 10]) + push!(ps, [:p => 10.0, :k1 => 1.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.0, :k5 => 1.0, :k6 => 1.0, :d => 1.0]) + push!(sps, :X4) + # Manual declaration of the second network (and its simulation settings). rate_2_1(u, p, t) = p[1] / 10 + p[1] * (u[1]^p[3]) / (u[1]^p[3] + p[2]^p[3]) rate_2_2(u, p, t) = p[4] * u[1] * u[2] rate_2_3(u, p, t) = p[5] * u[3] @@ -83,7 +93,11 @@ let push!(manual_networks, jumps_2) push!(u0_syms, [:X1, :X2, :X3]) push!(ps_syms, [:v, :K, :n, :k1, :k2, :k3, :d]) + push!(u0s, [:X1 => 10, :X2 => 10, :X3 => 10]) + push!(ps, [:v => 20.0, :K => 2.0, :n => 2, :k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :d => 1.0]) + push!(sps, :X3) + # Manual declaration of the third network (and its simulation settings). rate_3_1(u, p, t) = p[1] * binomial(u[1], 1) rate_3_2(u, p, t) = p[2] * binomial(u[2], 2) rate_3_3(u, p, t) = p[3] * binomial(u[2], 2) @@ -107,33 +121,29 @@ let push!(manual_networks, jumps_3) push!(u0_syms, [:X1, :X2, :X3, :X4]) push!(ps_syms, [:k1, :k2, :k3, :k4, :k5, :k6]) + push!(u0s, [:X1 => 15, :X2 => 5, :X3 => 5, :X4 => 5]) + push!(ps, [:k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :k4 => 0.1, :k5 => 0.1, :k6 => 0.1]) + push!(sps, :X4) # Loops through all cases, checks that identical simulations are generated with/without Catalyst. - for (rn_catalyst, rn_manual, u0_sym, ps_sym) in zip(catalyst_networks, manual_networks, u0_syms, ps_syms) - for factor in [5, 50] - seed = rand(rng, 1:100) - u0_1 = rnd_u0_Int64(rn_catalyst, rng; n = factor) - ps_1 = rnd_ps(rn_catalyst, rng; factor = factor/100.0) - dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 100.0), ps_1) - jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng) - sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0) - - u0_2 = map_to_vec(u0_1, u0_sym) - ps_2 = map_to_vec(ps_1, ps_sym) - dprob_2 = DiscreteProblem(u0_2, (0.0, 100.0), ps_2) - jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng) - sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0) + for (rn_catalyst, rn_manual, u0_sym, ps_sym, u0_1, ps_1, sp) in + zip(catalyst_networks, manual_networks, u0_syms, ps_syms, u0s, ps, sps) - if nameof(rn_catalyst) == :rnh7 - # Have spent a few hours figuring this one out. For certain seeds it actually works, - # but others not. This feels weird, and I didn't get any longer. I tried using non-random - # parameters/initial conditions, and removing the non-hill function reactions. Problem - # still persists. - @test_broken sol1[u0_sym] == sol2.u - else - @test sol1[u0_sym] == sol2.u - end - end + # Simulates the Catalyst-created model. + dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 10000.0), ps_1) + jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng) + sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0) + + # Simulates the manually written model + u0_2 = map_to_vec(u0_1, u0_sym) + ps_2 = map_to_vec(ps_1, ps_sym) + dprob_2 = DiscreteProblem(u0_2, (0.0, 10000.0), ps_2) + jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng) + sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0) + + # Checks that the means are similar (the test have been check that it holds across a large + # number of simulates, even without seed). + @test mean(sol1[sp]) ≈ mean(sol2[findfirst(u0_sym .== sp),:]) rtol = 1e-1 end end