diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 7e9af3e1dd..fc1871b19a 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -14,7 +14,7 @@ jobs: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: - version: '1' + version: '1.10.2' - name: Install dependencies run: julia --project=docs/ -e 'ENV["JULIA_PKG_SERVER"] = ""; using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy diff --git a/HISTORY.md b/HISTORY.md index d2c610a262..270d64b177 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -24,7 +24,7 @@ rn = @reaction_network begin @parameters η k, 2X --> X2, [noise_scaling=η] end -get_noise_scaling(rn) +getnoisescaling(rn) ``` - `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling). - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions. diff --git a/docs/Project.toml b/docs/Project.toml index 28d446f257..ebb086fda6 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -29,45 +29,45 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" -Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] -BenchmarkTools = "1" +BenchmarkTools = "1.5" BifurcationKit = "0.3" CairoMakie = "0.12" Catalyst = "13" -DataFrames = "1" -DiffEqParamEstim = "2.1" -DifferentialEquations = "7.7" +DataFrames = "1.6" +DiffEqParamEstim = "2.2" Distributions = "0.25" Documenter = "1.4.1" -DynamicalSystems = "3" -GlobalSensitivity = "2.4.0" -HomotopyContinuation = "2.6" +DynamicalSystems = "3.3" +GlobalSensitivity = "2.6" +HomotopyContinuation = "2.9" IncompleteLU = "0.2" -JumpProcesses = "9" -Latexify = "0.15, 0.16" -LinearSolve = "2" -ModelingToolkit = "9.5" -NonlinearSolve = "3.4.0" -Optim = "1" -Optimization = "3.19" -OptimizationBBO = "0.1.5, 0.2" -OptimizationNLopt = "0.1.8" -OptimizationOptimJL = "0.1.14" -OptimizationOptimisers = "0.1.1" -OrdinaryDiffEq = "6" -Plots = "1.36" -SciMLBase = "2.13" -SciMLSensitivity = "7.19" -Setfield = "1.1" -SpecialFunctions = "2.1" -SteadyStateDiffEq = "2.0.1" -StochasticDiffEq = "6" -StructuralIdentifiability = "0.5.1" -Symbolics = "5.14" +JumpProcesses = "9.11" +Latexify = "0.16" +LinearSolve = "2.30" +ModelingToolkit = "9.15" +NonlinearSolve = "3.12" +Optim = "1.9" +Optimization = "3.25" +OptimizationBBO = "0.2.1" +OptimizationNLopt = "0.2.1" +OptimizationOptimJL = "0.3.1" +OptimizationOptimisers = "0.2.1" +OrdinaryDiffEq = "6.80.1" +Plots = "1.40" +QuasiMonteCarlo = "0.3" +SciMLBase = "2.39" +SciMLSensitivity = "7.60" +SpecialFunctions = "2.4" +StaticArrays = "1.9" +SteadyStateDiffEq = "2.2" +StochasticDiffEq = "6.65" +StructuralIdentifiability = "0.5.7" +Symbolics = "5.28" diff --git a/docs/pages.jl b/docs/pages.jl index 941e906c3a..bd716f4327 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,38 +1,36 @@ pages = Any[ - "Home" => "home.md", + "Home" => "index.md", "Introduction to Catalyst" => Any[ "introduction_to_catalyst/catalyst_for_new_julia_users.md", - "introduction_to_catalyst/introduction_to_catalyst.md" + # "introduction_to_catalyst/introduction_to_catalyst.md" # Advanced introduction. ], "Model Creation and Properties" => Any[ "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", + #"model_creation/programmatic_CRN_construction.md", + #"model_creation/compositional_modeling.md", + #"model_creation/constraint_equations.md", # Events. - "model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions. + #"model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions. # Loading and writing models to files. - # Model visualisation. - "model_creation/network_analysis.md", + #"model_creation/model_visualisation.md", + #"model_creation/network_analysis.md", "model_creation/chemistry_related_functionality.md", "Model creation examples" => Any[ #"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_creation/examples/hodgkin_huxley_equation.md", + #"model_creation/examples/smoluchowski_coagulation_equation.md" ] ], "Model simulation" => Any[ "model_simulation/simulation_introduction.md", - # Simulation introduction. "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. + "model_simulation/ode_simulation_performance.md", # SDE Performance considerations/advice. # Jump Performance considerations/advice. # Finite state projection @@ -51,7 +49,7 @@ pages = Any[ # ODE parameter fitting using Turing. # SDE/Jump fitting. "inverse_problems/behaviour_optimisation.md", - "inverse_problems/structural_identifiability.md", + #"inverse_problems/structural_identifiability.md", # Broken on Julia v1.10.3, requires v1.10.2 or 1.10.4. # Practical identifiability. "inverse_problems/global_sensitivity_analysis.md", "Inverse problem examples" => Any[ @@ -67,6 +65,6 @@ pages = Any[ # # Contributor's guide. # # Repository structure. # ], - "FAQs" => "faqs.md", - "API" => "api.md" -] \ No newline at end of file + #"FAQs" => "faqs.md", + #"API" => "api.md" +] diff --git a/docs/src/api.md b/docs/src/api.md index 49b74adaa1..9e7086e6f8 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -152,16 +152,13 @@ accessor functions. ```@docs species nonspecies -reactionsystemparams reactions nonreactions numspecies numparams numreactions -numreactionsystemparams speciesmap paramsmap -reactionsystemparamsmap isspecies isautonomous Catalyst.isconstant diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml index e7454091d2..ebb086fda6 100644 --- a/docs/src/assets/Project.toml +++ b/docs/src/assets/Project.toml @@ -1,18 +1,26 @@ [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" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" 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" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b" OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1" @@ -21,37 +29,45 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" -Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] +BenchmarkTools = "1.5" BifurcationKit = "0.3" +CairoMakie = "0.12" Catalyst = "13" -DataFrames = "1" -DiffEqParamEstim = "2.1" -DifferentialEquations = "7.7" +DataFrames = "1.6" +DiffEqParamEstim = "2.2" Distributions = "0.25" Documenter = "1.4.1" -HomotopyContinuation = "2.6" -Latexify = "0.15, 0.16" -ModelingToolkit = "9.5" -NonlinearSolve = "3.4.0" -Optim = "1" -Optimization = "3.19" -OptimizationNLopt = "0.1.8" -OptimizationOptimJL = "0.1.14" -OptimizationOptimisers = "0.1.1" -OrdinaryDiffEq = "6" -Plots = "1.36" -SciMLBase = "2.13" -SciMLSensitivity = "7.19" -Setfield = "1.1" -SpecialFunctions = "2.1" -SteadyStateDiffEq = "2.0.1" -StochasticDiffEq = "6" -StructuralIdentifiability = "0.5.1" -Symbolics = "5.14" +DynamicalSystems = "3.3" +GlobalSensitivity = "2.6" +HomotopyContinuation = "2.9" +IncompleteLU = "0.2" +JumpProcesses = "9.11" +Latexify = "0.16" +LinearSolve = "2.30" +ModelingToolkit = "9.15" +NonlinearSolve = "3.12" +Optim = "1.9" +Optimization = "3.25" +OptimizationBBO = "0.2.1" +OptimizationNLopt = "0.2.1" +OptimizationOptimJL = "0.3.1" +OptimizationOptimisers = "0.2.1" +OrdinaryDiffEq = "6.80.1" +Plots = "1.40" +QuasiMonteCarlo = "0.3" +SciMLBase = "2.39" +SciMLSensitivity = "7.60" +SpecialFunctions = "2.4" +StaticArrays = "1.9" +SteadyStateDiffEq = "2.2" +StochasticDiffEq = "6.65" +StructuralIdentifiability = "0.5.7" +Symbolics = "5.28" diff --git a/docs/src/assets/long_ploting_times/model_creation/mm_kinetics.svg b/docs/src/assets/long_ploting_times/model_creation/mm_kinetics.svg new file mode 100644 index 0000000000..824a5fd376 --- /dev/null +++ b/docs/src/assets/long_ploting_times/model_creation/mm_kinetics.svg @@ -0,0 +1,128 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/long_ploting_times/model_creation/sir_outbreaks.svg b/docs/src/assets/long_ploting_times/model_creation/sir_outbreaks.svg new file mode 100644 index 0000000000..3e213ebbdd --- /dev/null +++ b/docs/src/assets/long_ploting_times/model_creation/sir_outbreaks.svg @@ -0,0 +1,128 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg b/docs/src/assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg new file mode 100644 index 0000000000..4f9f01fedf --- /dev/null +++ b/docs/src/assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/home.md b/docs/src/index.md similarity index 100% rename from docs/src/home.md rename to docs/src/index.md 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 82c297b33d..bfeee426cf 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 @@ -55,15 +55,15 @@ To import a Julia package into a session, you can use the `using PackageName` co using Pkg Pkg.add("Catalyst") ``` -Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we also wish to install the `DifferentialEquations` and `Plots` packages (for numeric simulation of models, and plotting, respectively). +Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we also wish to install the `OrdinaryDiffEq` and `Plots` packages (for numeric simulation of models, and plotting, respectively). ```julia -Pkg.add("DifferentialEquations") +Pkg.add("OrdinaryDiffEq") Pkg.add("Plots") ``` Once a package has been installed through the `Pkg.add` command, this command does not have to be repeated if we restart our Julia session. We can now import all three packages into our current session with: ```@example ex2 using Catalyst -using DifferentialEquations +using OrdinaryDiffEq using Plots ``` Here, if we restart Julia, these `using` commands *must be rerun*. @@ -77,11 +77,11 @@ Catalyst models are created through the `@reaction_network` *macro*. For more in The `@reaction_network` command is followed by the `begin` keyword, which is followed by one line for each *reaction* of the model. Each reaction consists of a *reaction rate*, followed by the reaction itself. The reaction contains a set of *substrates* and a set of *products* (what is consumed and produced by the reaction, respectively). These are separated by a `-->` arrow. Finally, the model ends with the `end` keyword. -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`. +Here, we create a simple [*birth-death* model](@ref basic_CRN_library_bd), 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/). @@ -130,12 +130,16 @@ 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. However, we do need to install, and then import, the JumpProcesses package (just to perform Gillespie, and other jump, simulations) +```julia +Pkg.add("JumpProcesses") +using JumpProcesses +``` -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. -* $I$, the amount of *infected* individuals. -* $R$, the amount of *recovered* (or *removed*) individuals. +This time, we will declare a so-called [SIR model for an infectious disease](@ref basic_CRN_library_sir). 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. +* *I*, the amount of *infected* individuals. +* *R*, the amount of *recovered* (or *removed*) individuals. It also has 2 reaction events: * Infection, where a susceptible individual meets an infected individual and also becomes infected. @@ -148,8 +152,8 @@ 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. @@ -164,12 +168,14 @@ nothing # hide Previously we have bundled this information into an `ODEProblem` (denoting a deterministic *ordinary differential equation*). Now we wish to simulate our model as a jump process (where each reaction event corresponds to a single jump in the state of the system). We do this by first creating a `DiscreteProblem`, and then using this as an input to a `JumpProblem`. ```@example ex2 +using JumpProcesses # hide dprob = DiscreteProblem(sir_model, u0, tspan, params) jprob = JumpProblem(sir_model, dprob, Direct()) +nothing # hide ``` 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. -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. For jump simulations, the `solve` function also requires 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 +200,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): @@ -231,7 +237,7 @@ So, why is this required, and why cannot we simply import any package installed 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. @@ -239,7 +245,7 @@ Some additional useful Pkg commands are: - `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 used 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. --- diff --git a/docs/src/inverse_problems/behaviour_optimisation.md b/docs/src/inverse_problems/behaviour_optimisation.md index 41c1e9eecf..a555037119 100644 --- a/docs/src/inverse_problems/behaviour_optimisation.md +++ b/docs/src/inverse_problems/behaviour_optimisation.md @@ -1,14 +1,14 @@ # [Optimization for non-data fitting purposes](@id behaviour_optimisation) In previous tutorials we have described how to use [PEtab.jl](@ref petab_parameter_fitting) and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. -## Maximising the pulse amplitude of an incoherent feed forward loop. +## [Maximising the pulse amplitude of an incoherent feed forward loop](@id behaviour_optimisation_IFFL_example) Incoherent feedforward loops (network motifs where a single component both activates and deactivates a downstream component) are able to generate pulses in response to step inputs[^2]. In this tutorial we will consider such an incoherent feedforward loop, attempting to generate a system with as prominent a response pulse as possible. Our model consists of 3 species: $X$ (the input node), $Y$ (an intermediary), and $Z$ (the output node). In it, $X$ activates the production of both $Y$ and $Z$, with $Y$ also deactivating $Z$. When $X$ is activated, there will be a brief time window where $Y$ is still inactive, and $Z$ is activated. However, as $Y$ becomes active, it will turn $Z$ off. This creates a pulse of $Z$ activity. To trigger the system, we create [an event](@ref ref), which increases the production rate of $X$ ($pX$) by a factor of $10$ at time $t = 10$. ```@example behaviour_optimization using Catalyst incoherent_feed_forward = @reaction_network begin - @discrete_event [10.0] ~ [p ~ 10*p] + @discrete_events [10.0] => [pX ~ 10*pX] pX, 0 --> X pY*X, 0 --> Y pZ*X/Y, 0 --> Z @@ -34,19 +34,20 @@ function pulse_amplitude(p, _) ps = Dict([:pX => p[1], :pY => p[2], :pZ => p[2]]) u0_new = [:X => ps[:pX], :Y => ps[:pX]*ps[:pY], :Z => ps[:pZ]/ps[:pY]^2] oprob_local = remake(oprob; u0= u0_new, p = ps) - sol = solve(oprob_local, verbose = false, maxiters = 10000) + sol = solve(oprob_local, Tsit5(); verbose = false, maxiters = 10000) SciMLBase.successful_retcode(sol) || return Inf return -(maximum(sol[:Z]) - sol[:Z][1]) end nothing # here ``` -This cost function takes two arguments (a parameter value `p`, and an additional one which we will ignore here but discuss later). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the cost function provided, input parameter set. While we could create a new `ODEProblem` within the cost function, cost functions are often called a large number of times during the optimisation process (making performance important). Here, using [`remake` on a previously created `ODEProblem`](@ref ref) is more performant than creating a new one. Just like [when using Optimization.jl to fit parameters to data](@ref optimization_parameter_fitting), we use the `verbose = false` option to prevent unnecessary simulation printouts, and a reduced `maxiters` value to reduce time spent simulating (for the model) unsuitable parameter sets. We also use `SciMLBase.successful_retcode(sol)` to check whether the simulation return code indicates a successful simulation (and if it did not, returns a large cost function value). Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our cost function return the negative pulse amplitude. +This cost function takes two arguments (a parameter value `p`, and an additional one which we will ignore here but discuss later). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the cost function provided, input parameter set. While we could create a new `ODEProblem` within the cost function, cost functions are often called a large number of times during the optimisation process (making performance important). Here, using [`remake` on a previously created `ODEProblem`](@ref simulation_structure_interfacing_problems_remake) is more performant than creating a new one. Just like [when using Optimization.jl to fit parameters to data](@ref optimization_parameter_fitting), we use the `verbose = false` option to prevent unnecessary simulation printouts, and a reduced `maxiters` value to reduce time spent simulating (for the model) unsuitable parameter sets. We also use `SciMLBase.successful_retcode(sol)` to check whether the simulation return code indicates a successful simulation (and if it did not, returns a large cost function value). Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our cost function return the negative pulse amplitude. Just like for [parameter fitting](@ref optimization_parameter_fitting), we create a `OptimizationProblem` using our cost function, and some initial guess of the parameter value. We also set upper and lower bounds for each parameter using the `lb` and `ub` optional arguments (in this case limiting each parameter's value to the interval $(0.1,10.0)$). ```@example behaviour_optimization using Optimization initial_guess = [1.0, 1.0, 1.0] opt_prob = OptimizationProblem(pulse_amplitude, initial_guess; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1]) +nothing # hide ``` !!! note As described in a [previous section on Optimization.jl](@ref optimization_parameter_fitting), `OptimizationProblem`s do not support setting parameter values using maps. We must instead set `initial_guess` values using a vector. Next, in the first line of our cost function, we reshape the parameter values to the common form used across Catalyst (e.g. `[:pX => p[1], :pY => p[2], :pZ => p[2]]`, however, here we use a dictionary to easier compute the steady state initial condition). We also note that the order used in this array corresponds to the order we give each parameter's bounds in `lb` and `ub`, and the order in which their values occur in the output solution. @@ -55,13 +56,14 @@ As [previously described](@ref optimization_parameter_fitting), Optimization.jl ```@example behaviour_optimization using OptimizationBBO opt_sol = solve(opt_prob, BBO_adaptive_de_rand_1_bin_radiuslimited()) +nothing # hide ``` Finally, we plot a simulation using the found parameter set (stored in `opt_sol.u`): ```@example behaviour_optimization ps_res = Dict([:pX => opt_sol.u[1], :pY => opt_sol.u[2], :pZ => opt_sol.u[2]]) u0_res = [:X => ps_res[:pX], :Y => ps_res[:pX]*ps_res[:pY], :Z => ps_res[:pZ]/ps_res[:pY]^2] oprob_res = remake(oprob; u0 = u0_res, p = ps_res) -sol_res = solve(oprob_res) +sol_res = solve(oprob_res, Tsit5()) plot(sol_res; idxs=:Z) ``` For this model, it turns out that $Z$'s maximum pulse amplitude is equal to twice its steady state concentration. Hence, the maximisation of its pulse amplitude is equivalent to maximising its steady state concentration. @@ -71,16 +73,17 @@ For this model, it turns out that $Z$'s maximum pulse amplitude is equal to twic There are several modifications to our problem where it would actually have parameters. E.g. our model might have had additional parameters (e.g. a degradation rate) which we would like to keep fixed throughout the optimisation process. If we then would like to run the optimisation process for several different values of these fixed parameters, we could have made them parameters to our `OptimizationProblem` (and their values provided as a third argument, after `initial_guess`). -## Utilising automatic differentiation +## [Utilising automatic differentiation](@id behaviour_optimisation_AD) Optimisation methods can be divided into differentiation-free and differentiation-based optimisation methods. E.g. consider finding the minimum of the function $f(x) = x^2$, given some initial guess of $x$. Here, we can simply compute the differential and descend along it until we find $x=0$ (admittedly, for this simple problem the minimum can be computed directly). This principle forms the basis of optimisation methods such as [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent), which utilises information of a function's differential to minimise it. When attempting to find a global minimum, to avoid getting stuck in local minimums, these methods are often augmented by additional routines. While the differentiation of most algebraic functions is trivial, it turns out that even complicated functions (such as the one we used above) can be differentiated computationally through the use of [*automatic differentiation* (AD)](https://en.wikipedia.org/wiki/Automatic_differentiation). Through packages such as [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl), [ReverseDiff.jl](https://github.com/JuliaDiff/ReverseDiff.jl), and [Zygote.jl](https://github.com/FluxML/Zygote.jl), Julia supports AD for most code. Specifically for code including simulation of differential equations, differentiation is supported by [SciMLSensitivity.jl](https://github.com/SciML/SciMLSensitivity.jl). Generally, AD can be used without specific knowledge from the user, however, it requires an additional step in the construction of our `OptimizationProblem`. Here, we create a [specialised `OptimizationFunction` from our cost function](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#optfunction). To it, we will also provide our choice of AD method. There are [several alternatives](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Automatic-Differentiation-Construction-Choice-Recommendations), and in our case we will use `AutoForwardDiff()` (a good choice for small optimisation problems). We can then create a new `OptimizationProblem` using our updated cost function: ```@example behaviour_optimization opt_func = OptimizationFunction(pulse_amplitude, AutoForwardDiff()) opt_prob = OptimizationProblem(opt_func, initial_guess; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1]) +nothing # hide ``` Finally, we can find the optimum using some differentiation-based optimisation methods. Here we will use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s `BFGS` method: -```@example behaviour_optimization +```julia using OptimizationOptimJL opt_sol = solve(opt_prob, OptimizationOptimJL.BFGS()) ``` diff --git a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md index 47f4623514..f5eaa4c65f 100644 --- a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md +++ b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md @@ -2,7 +2,7 @@ In this example we will use [Optimization.jl](https://github.com/SciML/Optimization.jl) to fit the parameters of an oscillatory system (the Brusselator) to data. Here, special consideration is taken to avoid reaching a local minimum. Instead of fitting the entire time series directly, we will start with fitting parameter values for the first period, and then use those as an initial guess for fitting the next (and then these to find the next one, and so on). Using this procedure is advantageous for oscillatory systems, and enables us to reach the global optimum. First, we fetch the required packages. -```@example pe1 +```@example pe_osc_example using Catalyst using OrdinaryDiffEq using Optimization @@ -11,7 +11,7 @@ using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic diff ``` Next, we declare our model, the Brusselator oscillator. -```@example pe1 +```@example pe_osc_example brusselator = @reaction_network begin A, ∅ --> X 1, 2X + Y --> 3X @@ -24,7 +24,7 @@ nothing # hide We simulate our model, and from the simulation generate sampled data points (to which we add noise). We will use this data to fit the parameters of our model. -```@example pe1 +```@example pe_osc_example u0 = [:X => 1.0, :Y => 1.0] tspan = (0.0, 30.0) @@ -37,7 +37,7 @@ nothing # hide ``` We can plot the real solution, as well as the noisy samples. -```@example pe1 +```@example pe_osc_example using Plots default(; lw = 3, framestyle = :box, size = (800, 400)) @@ -49,7 +49,7 @@ Next, we create a function to fit the parameters using the `ADAM` optimizer. For a given initial estimate of the parameter values, `pinit`, this function will fit parameter values, `p`, to our data samples. We use `tend` to indicate the time interval over which we fit the model. -```@example pe1 +```@example pe_osc_example function optimise_p(pinit, tend) function loss(p, _) newtimes = filter(<=(tend), sample_times) @@ -71,12 +71,12 @@ nothing # hide ``` Next, we will fit a parameter set to the data on the interval `(0, 10)`. -```@example pe1 +```@example pe_osc_example p_estimate = optimise_p([5.0, 5.0], 10.0) ``` We can compare this to the real solution, as well as the sample data -```@example pe1 +```@example pe_osc_example newprob = remake(prob; tspan = (0., 10.), p = p_estimate) sol_estimate = solve(newprob, Rosenbrock23()) plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2) @@ -88,7 +88,7 @@ plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash, Next, we use this parameter estimate as the input to the next iteration of our fitting process, this time on the interval `(0, 20)`. -```@example pe1 +```@example pe_osc_example p_estimate = optimise_p(p_estimate, 20.) newprob = remake(prob; tspan = (0., 20.), p = p_estimate) sol_estimate = solve(newprob, Rosenbrock23()) @@ -101,20 +101,20 @@ plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash, Finally, we use this estimate as the input to fit a parameter set on the full time interval of the sampled data. -```@example pe1 +```@example pe_osc_example p_estimate = optimise_p(p_estimate, 30.0) newprob = remake(prob; tspan = (0., 30.0), p = p_estimate) sol_estimate = solve(newprob, Rosenbrock23()) plot(sol_real; color = [:blue :red], label = ["X real" "Y real"], linealpha = 0.2) scatter!(sample_times, sample_vals'; color = [:blue :red], - label = ["Samples of X" "Samples of Y"], alpha = 0.4) + label = ["Samples of X" "Samples of Y"], alpha = 0.4) plot!(sol_estimate; color = [:darkblue :darkred], linestyle = :dash, label = ["X estimated" "Y estimated"], xlimit = tspan) ``` The final parameter estimate is then -```@example pe1 +```@example pe_osc_example p_estimate ``` which is close to the actual parameter set of `[1.0, 2.0]`. @@ -125,7 +125,7 @@ then extend the interval, is to avoid getting stuck in a local minimum. Here specifically, we chose our initial interval to be smaller than a full cycle of the oscillation. If we had chosen to fit a parameter set on the full interval immediately we would have obtained poor fit and an inaccurate estimate for the parameters. -```@example pe1 +```@example pe_osc_example p_estimate = optimise_p([5.0,5.0], 30.0) newprob = remake(prob; tspan = (0.0,30.0), p = p_estimate) diff --git a/docs/src/inverse_problems/global_sensitivity_analysis.md b/docs/src/inverse_problems/global_sensitivity_analysis.md index d611997a9f..ee20b6a802 100644 --- a/docs/src/inverse_problems/global_sensitivity_analysis.md +++ b/docs/src/inverse_problems/global_sensitivity_analysis.md @@ -11,7 +11,7 @@ A related concept to global sensitivity is *local sensitivity*. This, rather tha While local sensitivities are primarily used as a subroutine of other methodologies (such as optimisation schemes), it also has direct uses. E.g., in the context of fitting parameters to data, local sensitivity analysis can be used to, at the parameter set of the optimal fit, [determine the cost function's sensitivity to the system parameters](@ref ref). ## [Basic example](@id global_sensitivity_analysis_basic_example) -We will consider a simple [SEIR model of an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). This is an expansion of the classic [SIR model](@ref ref) with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection to others. +We will consider a simple [SEIR model of an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). This is an expansion of the classic [SIR model](@ref basic_CRN_library_sir) with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection to others. ```@example gsa_1 using Catalyst seir_model = @reaction_network begin @@ -53,7 +53,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0 !!! note We should make a couple of notes about the example above: - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_logarithmic_scale), this is advantageous in the context of inverse problems such as this one. - - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref ref) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set. + - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set. - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`. - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`. diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md index c579a68aa7..6dd50c3e4a 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) +# [Parameter Fitting for ODEs using Optimization.jl and DiffEqParamEstim.jl](@id optimization_parameter_fitting) Fitting parameters to data involves solving an optimisation problem (that is, finding the parameter set that optimally fits your model to your data, typically by minimising a cost function). The SciML ecosystem's primary package for solving optimisation problems is [Optimization.jl](https://github.com/SciML/Optimization.jl). It provides access to a variety of solvers via a single common interface by wrapping a large number of optimisation libraries that have been implemented in Julia. -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/). +This tutorial demonstrates both how to create parameter fitting cost functions using the [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl) package, and how to use Optimization.jl to minimise these. Optimization.jl can also be used in other contexts, such as [finding parameter sets that maximise the magnitude of some system behaviour](@ref behaviour_optimisation). More details on how to use these packages can be found in their [respective](https://docs.sciml.ai/Optimization/stable/) [documentations](https://docs.sciml.ai/DiffEqParamEstim/stable/). ## [Basic example](@id optimization_parameter_fitting_basics) -Let us consider a [Michaelis-Menten enzyme kinetics model](@ref ref), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$): +Let us consider a [Michaelis-Menten enzyme kinetics model](@ref basic_CRN_library_mm), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$): ```@example diffeq_param_estim_1 using Catalyst rn = @reaction_network begin @@ -30,8 +30,8 @@ data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] # 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) +plot(true_sol; idxs = :P, label = "True solution", lw = 8) +plot!(data_ts, data_vals; label = "Measurements", seriestype=:scatter, ms = 6, color = :blue) ``` Next, we will use DiffEqParamEstim to build a loss function to measure how well our model's solutions fit the data. @@ -40,12 +40,12 @@ using DiffEqParamEstim, Optimization 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](@ref ref) we wish to simulate our model with. +- `Tsit5()`: The [numeric solver](@ref simulation_intro_solver_options) 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. @@ -63,27 +63,27 @@ 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 [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationNLopt 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 NLopt.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()) +using OptimizationNLopt +optsol = solve(optprob, NLopt.LN_NELDERMEAD()) nothing # hide ``` We can now simulate our model for the corresponding parameter set, checking that it fits our data. ```@example diffeq_param_estim_1 -oprob_fitted = remake(oprob; p=optsol.u) +oprob_fitted = remake(oprob; p = optsol.u) fitted_sol = solve(oprob_fitted, Tsit5()) -plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=6, color=:lightblue) +plot!(fitted_sol; idxs = :P, label = "Fitted solution", linestyle = :dash, lw = 6, color = :lightblue) ``` !!! note Here, a good exercise is to check the resulting parameter set and note that, while it creates a good fit to the data, it does not actually correspond to the original parameter set. [Identifiability](@ref structural_identifiability) is a concept that studies how to deal with this problem. -Say that we instead would like to use the [Broyden–Fletcher–Goldfarb–Shannon](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) algorithm, as implemented by the [NLopt.jl](https://github.com/JuliaOpt/NLopt.jl) package. In this case we would run: +Say that we instead would like to use the [Broyden–Fletcher–Goldfarb–Shannon](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) algorithm, as implemented by the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) package. In this case we would run: ```@example diffeq_param_estim_1 -using OptimizationNLopt -sol = solve(optprob, NLopt.LD_LBFGS()) +using OptimizationOptimJL +sol = solve(optprob, Optim.LBFGS()) nothing # hide ``` diff --git a/docs/old_files/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md similarity index 100% rename from docs/old_files/petab_ode_param_fitting.md rename to docs/src/inverse_problems/petab_ode_param_fitting.md diff --git a/docs/src/inverse_problems/structural_identifiability.md b/docs/src/inverse_problems/structural_identifiability.md index 95e5de615a..8b118b367c 100644 --- a/docs/src/inverse_problems/structural_identifiability.md +++ b/docs/src/inverse_problems/structural_identifiability.md @@ -31,27 +31,27 @@ Global identifiability can be assessed using the `assess_identifiability` functi - Unidentifiable. To it, we provide our `ReactionSystem` model and a list of quantities that we are able to measure. Here, we consider a Goodwind oscillator (a simple 3-component model, where the three species $M$, $E$, and $P$ are produced and degraded, which may exhibit oscillations)[^2]. Let us say that we are able to measure the concentration of $M$, we then designate this using the `measured_quantities` argument. We can now assess identifiability in the following way: -```example si1 +```@example si1 using Catalyst, Logging, StructuralIdentifiability -goodwind_oscillator = @reaction_network begin +gwo = @reaction_network begin (pₘ/(1+P), dₘ), 0 <--> M (pₑ*M,dₑ), 0 <--> E (pₚ*E,dₚ), 0 <--> P end -assess_identifiability(goodwind_oscillator; measured_quantities=[:M], loglevel=Logging.Error) +assess_identifiability(gwo; measured_quantities = [:M], loglevel = Logging.Error) ``` -From the output, we find that `E(t)`, `pₑ`, and `pₚ` (the trajectory of $E$, and the production rates of $E$ and $P$, respectively) are non-identifiable. Next, `dₑ` and `dₚ` (the degradation rates of $E$ and $P$, respectively) are locally identifiable. Finally, `P(t)`, `M(t)`, `pₘ`, and `dₘ` (the trajectories of `P` and `M`, and the production and degradation rate of `M`, respectively) are all globally identifiable. We note that we also imported the Logging.jl package, and provided the `loglevel=Logging.Error` input argument. StructuralIdentifiability functions generally provide a large number of output messages. Hence, we will use this argument (which requires the Logging package) throughout this tutorial to decrease the amount of printed text. +From the output, we find that `E(t)`, `pₑ`, and `pₚ` (the trajectory of $E$, and the production rates of $E$ and $P$, respectively) are non-identifiable. Next, `dₑ` and `dₚ` (the degradation rates of $E$ and $P$, respectively) are locally identifiable. Finally, `P(t)`, `M(t)`, `pₘ`, and `dₘ` (the trajectories of `P` and `M`, and the production and degradation rate of `M`, respectively) are all globally identifiable. We note that we also imported the Logging.jl package, and provided the `loglevel = Logging.Error` input argument. StructuralIdentifiability functions generally provide a large number of output messages. Hence, we will use this argument (which requires the Logging package) throughout this tutorial to decrease the amount of printed text. Next, we also assess identifiability in the case where we can measure all three species concentrations: -```example si1 -assess_identifiability(goodwind_oscillator; measured_quantities=[:M, :P, :E], loglevel=Logging.Error) +```@example si1 +assess_identifiability(gwo; measured_quantities = [:M, :P, :E], loglevel = Logging.Error) ``` in which case all species trajectories and parameters become identifiable. ### Indicating known parameters In the previous case we assumed that all parameters are unknown, however, this is not necessarily true. If there are parameters with known values, we can supply these using the `known_p` argument. Providing this additional information might also make other, previously unidentifiable, parameters identifiable. Let us consider the previous example, where we measure the concentration of $M$ only, but now assume we also know the production rate of $E$ ($pₑ$): -```example si1 -assess_identifiability(gwo; measured_quantities=[:M], known_p=[:pₑ], loglevel=Logging.Error) +```@example si1 +assess_identifiability(gwo; measured_quantities = [:M], known_p = [:pₑ], loglevel = Logging.Error) ``` Not only does this turn the previously non-identifiable `pₑ` (globally) identifiable (which is obvious, given that its value is now known), but this additional information improve identifiability for several other network components. @@ -59,47 +59,46 @@ To, in a similar manner, indicate that certain initial conditions are known is a ### Providing non-trivial measured quantities Sometimes, ones may not have measurements of species, but rather some combinations of species (or possibly parameters). To account for this, `measured_quantities` accepts any algebraic expression (and not just single species). To form such expressions, species and parameters have to first be `@unpack`'ed from the model. Say that we have a model where an enzyme ($E$) is converted between an active and inactive form, which in turns activates the production of a product, $P$: -```example si2 -using Catalyst, StructuralIdentifiability # hide -enzyme_activation = @reaction_network begin +```@example si1 +rs = @reaction_network begin (kA,kD), Eᵢ <--> Eₐ (Eₐ, d), 0 <-->P end ``` If we can measure the total amount of $E$ ($=Eᵢ+Eₐ$), as well as the amount of $P$, we can use the following to assess identifiability: -```example si2 -@unpack Eᵢ, Eₐ = enzyme_activation -assess_identifiability(enzyme_activation; measured_quantities=[Eᵢ+Eₐ, :P], loglevel=Logging.Error) +```@example si1 +@unpack Eᵢ, Eₐ = rs +assess_identifiability(rs; measured_quantities = [Eᵢ + Eₐ, :P], loglevel = Logging.Error) nothing # hide ``` ### Assessing identifiability for specified quantities only By default, StructuralIdentifiability assesses identifiability for all parameters and variables. It is, however, possible to designate precisely which quantities you want to check using the `funcs_to_check` option. This both includes selecting a smaller subset of parameters and variables to check, or defining customised expressions. Let us consider the Goodwind from previously, and say that we would like to check whether the production parameters ($pₘ$, $pₑ$, and $pₚ$) and the total amount of the three species ($P + M + E$) are identifiable quantities. Here, we would first unpack these (allowing us to form algebraic expressions) and then use the following code: -```example si1 -@unpack pₘ, pₑ, pₚ, M, E, P = goodwind_oscillator -assess_identifiability(goodwind_oscillator; measured_quantities=[:M], funcs_to_check=[pₘ, pₑ, pₚ, M + E + P], loglevel=Logging.Error) +```@example si1 +@unpack pₘ, pₑ, pₚ, M, E, P = gwo +assess_identifiability(gwo; measured_quantities = [:M], funcs_to_check = [pₘ, pₑ, pₚ, M + E + P], loglevel = Logging.Error) nothing # hide ``` ### Probability of correctness The identifiability methods used can, in theory, produce erroneous results. However, it is possible to adjust the lower bound for the probability of correctness using the argument `prob_threshold` (by default set to `0.99`, that is, at least a $99\%$ chance of correctness). We can e.g. increase the bound through: -```example si2 -assess_identifiability(goodwind_oscillator; measured_quantities=[:M], prob_threshold=0.999, loglevel=Logging.Error) +```@example si1 +assess_identifiability(gwo; measured_quantities=[:M], prob_threshold = 0.999, loglevel = Logging.Error) nothing # hide ``` giving a minimum bound of $99.9\%$ chance of correctness. In practise, the bounds used by StructuralIdentifiability are very conservative, which means that while the minimum guaranteed probability of correctness in the default case is $99\%$, in practise it is much higher. While increasing the value of `prob_threshold` increases the certainty of correctness, it will also increase the time required to assess identifiability. ## Local identifiability analysis Local identifiability can be assessed through the `assess_local_identifiability` function. While this is already determined by `assess_identifiability`, assessing local identifiability only has the advantage that it is easier to compute. Hence, there might be models where global identifiability analysis fails (or takes a prohibitively long time), where instead `assess_local_identifiability` can be used. This function takes the same inputs as `assess_identifiability` and returns, for each quantity, `true` if it is locally identifiable (or `false` if it is not). Here, for the Goodwind oscillator, we assesses it for local identifiability only: -```example si1 -assess_local_identifiability(goodwind_oscillator; measured_quantities=[:M], loglevel=Logging.Error) +```@example si1 +assess_local_identifiability(gwo; measured_quantities = [:M], loglevel = Logging.Error) ``` We note that the results are consistent with those produced by `assess_identifiability` (with globally or locally identifiable quantities here all being assessed as at least locally identifiable). ## Finding identifiable functions Finally, StructuralIdentifiability provides the `find_identifiable_functions` function. Rather than determining the identifiability of each parameter and unknown of the model, it finds a set of identifiable functions, such as any other identifiable expression of the model can be generated by these. Let us again consider the Goodwind oscillator, using the `find_identifiable_functions` function we find that identifiability can be reduced to five globally identifiable expressions: -```example si1 -find_identifiable_functions(goodwind_oscillator; measured_quantities=[:M], loglevel=Logging.Error) +```@example si1 +find_identifiable_functions(gwo; measured_quantities = [:M], loglevel = Logging.Error) ``` Again, these results are consistent with those produced by `assess_identifiability`. There, `pₑ` and `pₚ` where found to be globally identifiable. Here, they correspond directly to identifiable expressions. The remaining four parameters (`pₘ`, `dₘ`, `dₑ`, and `dₚ`) occur as part of more complicated composite expressions. @@ -107,34 +106,34 @@ Again, these results are consistent with those produced by `assess_identifiabili ## Creating StructuralIdentifiability compatible ODE models from Catalyst `ReactionSystem`s While the functionality described above covers the vast majority of analysis that user might want to perform, the StructuralIdentifiability package supports several additional features. While these does not have inherent Catalyst support, we do provide the `make_si_ode` function to simplify their use. Similar to the previous functions, it takes a `ReactionSystem`, lists of measured quantities, and known parameter values. The output is a [ODE of the standard form supported by StructuralIdentifiability](https://docs.sciml.ai/StructuralIdentifiability/stable/tutorials/creating_ode/#Defining-the-model-using-@ODEmodel-macro). It can be created using the following syntax: -```example si1 -si_ode = make_si_ode(goodwind_oscillator; measured_quantities=[:M]) +```@example si1 +si_ode = make_si_ode(gwo; measured_quantities = [:M]) nothing # hide ``` and then used as input to various StructuralIdentifiability functions. In the following example we use StructuralIdentifiability's `print_for_DAISY` function, printing the model as an expression that can be used by the [DAISY](https://daisy.dei.unipd.it/) software for identifiability analysis[^3]. -```example si1 +```@example si1 print_for_DAISY(si_ode) nothing # hide ``` ## Notes on systems with conservation laws Several reaction network models, such as -```example si2 +```@example si3 using Catalyst, Logging, StructuralIdentifiability # hide rs = @reaction_network begin - (k1,k2), X1 <--> X2 + (k1,k2), X1 <--> X2 end ``` contain conservation laws (in this case $Γ = X1 + X2$, where $Γ = X1(0) + X2(0)$ is a constant). Because the presence of such conservation laws makes structural identifiability analysis prohibitively computationally expensive (for all but the simplest of cases), these are automatically eliminated by Catalyst (removing one ODE from the resulting ODE system for each conservation law). For the `assess_identifiability` and `assess_local_identifiability` functions, this will be unnoticed by the user. However, for the `find_identifiable_functions` and `make_si_ode` functions, this may result in one, or several, parameters of the form `Γ[i]` (where `i` is an integer) appearing in the produced expressions. These correspond to the conservation law constants and can be found through -```example si2 +```@example si3 conservedequations(rs) ``` E.g. if you run: -```example si2 +```@example si3 find_identifiable_functions(rs; measured_quantities = [:X1, :X2]) ``` we see that `Γ[1]` (`= X1(0) + X2(0)`) is detected as an identifiable expression. If we want to disable this feature for any function, we can use the `remove_conserved = false` option: -```example si2 +```@example si3 find_identifiable_functions(rs; measured_quantities = [:X1, :X2], remove_conserved = false) ``` @@ -144,7 +143,7 @@ Structural identifiability cannot currently be applied to systems with parameter rn = @reaction_network begin (hill(X,v,K,n),d), 0 <--> X end -assess_identifiability(rn; measured_quantities=[:X]) +assess_identifiability(rn; measured_quantities = [:X]) ``` is currently not possible. Hopefully this will be a supported feature in the future. For now, such expressions will have to be rewritten to not include such exponents. For some cases, e.g. `10^k` this is trivial. However, it is also possible generally (but more involved and often includes introducing additional variables). diff --git a/docs/src/model_creation/chemistry_related_functionality.md b/docs/src/model_creation/chemistry_related_functionality.md index ff6dfd9373..f87402d91c 100644 --- a/docs/src/model_creation/chemistry_related_functionality.md +++ b/docs/src/model_creation/chemistry_related_functionality.md @@ -1,9 +1,10 @@ # [Chemistry-related functionality](@id chemistry_functionality) -While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems. +While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common in chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems. - The `@compound` option, which enables the user to designate that a specific species is composed of certain subspecies. - The `balance_reaction` function, which enables the user to balance a reaction so the same number of components occur on both sides. + ## Modelling with compound species ### Creating compound species programmatically @@ -11,13 +12,14 @@ We will first show how to create compound species through [programmatic model co ```@example chem1 using Catalyst t = default_t() -@species C(t) O(t) +@species C(t) O(t) +nothing # hide ``` Next, we create the `CO2` compound species: ```@example chem1 @compound CO2 ~ C + 2O ``` -Here, the compound is the first argument to the macro, followed by its component (with the left-hand and right-hand sides separated by a `~` sign). While non-compound species (such as `C` and `O`) have their independent variable (in this case `t`) designated, independent variables are generally not designated for compounds (these are instead directly inferred from their components). Components with non-unit stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is also a species, and can be used wherever e.g. `C` can be used: +Here, the compound is the first argument to the macro, followed by its component (with the left-hand and right-hand sides separated by a `~` sign). While non-compound species (such as `C` and `O`) have their independent variable (in this case `t`) designated, independent variables are generally not designated for compounds (these are instead directly inferred from their components). Components with non-unitary stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is also a species, and can be used wherever e.g. `C` can be used: ```@example chem1 isspecies(CO2) ``` @@ -32,7 +34,7 @@ Alternatively, we can retrieve the components and their stoichiometric coefficie ```@example chem1 component_coefficients(CO2) ``` -Finally, it is possible to check whether a species is a compound or not using the `iscompound` function: +Finally, it is possible to check whether a species is a compound using the `iscompound` function: ```@example chem1 iscompound(CO2) ``` @@ -68,7 +70,7 @@ end ``` When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` or `@compound` options, or by appearing in a reaction*. E.g. the following is not valid ```julia -rn = @reaction_network begin`` +rn = @reaction_network begin @compounds begin C2O ~ C + 2O H2O ~ 2H + O @@ -77,25 +79,28 @@ rn = @reaction_network begin`` (k1,k2), H2O+ CO2 <--> H2CO3 end ``` -as the components `C`, `H`, and `O` are not declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`. +as the components `C`, `H`, and `O` are not declared as species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`. ### Designating metadata and default values for compounds Just like for normal species, it is possible to designate metadata and default values for compounds. Metadata is provided after the compound name, but separated from it by a `,`: ```@example chem1 @compound (CO2, [unit="mol"]) ~ C + 2O +nothing # hide ``` Default values are designated using `=`, and provided directly after the compound name.: ```@example chem1 @compound (CO2 = 2.0) ~ C + 2O +nothing # hide ``` If both default values and meta data are provided, the metadata is provided after the default value: ```@example chem1 @compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O +nothing # hide ``` -In all of these cases, the side to the left of the `~` must be enclosed within `()`. +In all of these cases, the left-hand side must be enclosed within `()`. ### Compounds with multiple independent variables -While we generally do not need to specify independent variables for compound, if the components (together) have more than one independent variable, this have to be done: +While we generally do not need to specify independent variables for compound, if the components (together) have more than one independent variable, this *must be done*: ```@example chem1 t = default_t() @variables s @@ -113,7 +118,7 @@ Here, the reaction rate (`k`) is not involved in the reaction balancing. We use ```@example chem1 balance_reaction(rx) ``` -which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible depending on the given reaction. +which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that, in some cases, the reaction balancing problem does not have a single obvious solution. Typically, a single solution is the obvious candidate (in which case this is the vector's only element). However, when this is not the case, the vector instead contain several reactions (from which a balanced reaction cab be generated). Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction: ```@example chem2 diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md index 8a1026bbae..1369efe04e 100644 --- a/docs/src/model_creation/dsl_advanced.md +++ b/docs/src/model_creation/dsl_advanced.md @@ -9,10 +9,10 @@ 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: +Previously, 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 dsl_description_stoichiometries_parameters)). 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ₐ + k*E, Pᵢ --> Pₐ end ``` `E` (as well as `k`) will be considered a parameter, something we can confirm directly: @@ -22,19 +22,19 @@ 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ₐ + @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)). + 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 dsl_advanced_options_ivs)). 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ₐ + @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. @@ -44,37 +44,37 @@ While designating something which would default to a parameter as a species is s 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ₐ + @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 + (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 + @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. + 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 can be found [here](@ref optimization_parameter_fitting_basics)). 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. @@ -85,7 +85,7 @@ Generally, there are four main reasons for specifying species/parameters using t 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 +!!! 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) @@ -93,9 +93,9 @@ When declaring species/parameters using the `@species` and `@parameters` options ```@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 + @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: @@ -120,12 +120,12 @@ It is also possible to declare a model with default values for only some initial ```@example dsl_advanced_defaults using Catalyst # hide rn = @reaction_network begin - @species X(t)=1.0 - (p,d), 0 <--> X + @species X(t)=1.0 + (p,d), 0 <--> X end tspan = (0.0, 10.0) -p = [:p => 1.0, :D => 0.2] +p = [:p => 1.0, :d => 0.2] oprob = ODEProblem(rn, u0, tspan, p) sol = solve(oprob) plot(sol) @@ -136,9 +136,9 @@ API for checking the default values of species and parameters can be found [here 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 + @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: @@ -166,41 +166,41 @@ Whenever a species/parameter is declared using the `@species`/`@parameters` opti ```@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 + @species Xᵢ(t) [description="X's inactive form"] Xₐ(t) [description=" X's active form"] + @parameters kA [description="Activation rate"] kD [description="Deactivation rate"] + (ka,kD), Xᵢ <--> Xₐ 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), +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](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 + @parameters kA [description="Activation rate", bounds=(0.01,10.0)] + (ka,kD), Xᵢ <--> Xₐ 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 + @parameters kA=1.0 [description="Activation rate", bounds=(0.01,10.0)] + (ka,kD), Xᵢ <--> Xₐ 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 + @parameters begin + kA, [description="Activation rate", bounds=(0.01,10.0)] + kD = 1.0, [description="Deactivation rate"] + end + (kA,kD), Xᵢ <--> Xₐ 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): +Each metadata has its own getter functions. E.g. we can get the description of the parameter `kA` using `ModelingToolkit.getdescription` (here we use [system indexing](@ref ref) to access the parameter): ```@example dsl_advanced_metadata -getdescription(two_state_system.kA) +ModelingToolkit.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). @@ -211,8 +211,8 @@ Catalyst enables the designation of parameters as `constantspecies`. These param ```@example dsl_advanced_constant_species using Catalyst # hide rn = @reaction_network begin - @parameters X [isconstantspecies=true] - k, X --> Xᴾ + @parameters X [isconstantspecies=true] + k, X --> Xᴾ end ``` We can confirm that $Xᴾ$ is the only species of the system: @@ -222,14 +222,14 @@ 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ᴾ + 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`) +Sometimes it is desired to designate that a parameter should have a specific [type](https://docs.julialang.org/en/v1/manual/types/). 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 @@ -238,7 +238,7 @@ polymerisation_network = @reaction_network begin 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.: +Generally, when simulating models with mixed parameter types, it is recommended to [declare parameter values as tuples, rather than vectors](@ref simulation_intro_ODEs_input_forms), e.g.: ```@example dsl_advanced_parameter_types ps = (:kB => 0.2, :kD => 1.0, :n => 2) nothing # hide @@ -247,14 +247,14 @@ 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."] + @parameters n::Int64 = 2 [description="Parameter n, an integer with defaults 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: +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 basic_CRN_library_two_states). 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 @@ -279,14 +279,14 @@ Each reaction network model has a name. It can be accessed using the `nameof` fu ```@example dsl_advanced_names using Catalyst # hide rn = @reaction_network begin - (p,d), 0 <--> X + (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 + (p,d), 0 <--> X end nameof(rn) ``` @@ -294,10 +294,10 @@ 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 + (p,d), 0 <--> X end rn2 = @reaction_network begin - (p,d), 0 <--> X + (p,d), 0 <--> X end rn1 == rn2 ``` @@ -308,15 +308,16 @@ 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 + (p,d), 0 <--> X end rn2 = @reaction_network my_network begin - (p,d), 0 <--> X + (p,d), 0 <--> X end rn1 == rn2 ``` +If you wish to check for identity, and wish that models that have different names but are otherwise identical, should be considered equal, you can use the [`isequivalent`](@ref) function. -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. +Setting model names is primarily useful for [hierarchical modelling](@ref compositional_modeling), 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. @@ -325,11 +326,11 @@ Let us consider a model where two species (`X` and `Y`) can bind to form a compl ```@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 + @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): @@ -346,30 +347,32 @@ 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] +nothing # hide ``` 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]) +plot(sol; idxs = :Xtot) +plot!(ylimit = (minimum(sol[:Xtot])*0.95, maximum(sol[:Xtot])*1.05)) # hide ``` 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: +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 dsl_description_stoichiometries_parameters)) `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 + @observables Xtot ~ X + n*Xn + (kB,kD), n*X <--> Xn end nothing # hide ``` -!!! +!!! note 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 + @observables (Xtot, [description="The total amount of X in the system."]) ~ X + n*Xn + (kB,kD), n*X <--> Xn end nothing # hide ``` @@ -377,9 +380,9 @@ 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 + @species Xtot(t) + @observables Xtot ~ X + n*Xn + (kB,kD), n*X <--> Xn end nothing # hide ``` @@ -396,12 +399,12 @@ As [described elsewhere](@ref ref), Catalyst's `ReactionSystem` models depend on ```@example dsl_advanced_ivs using Catalyst # hide rn = @reaction_network begin - @ivs τ - (ka,kD), Xi <--> Xa + @ivs τ + (ka,kD), Xᵢ <--> Xₐ end nothing # hide ``` -We can confirm that `Xi` and `Xa` depend on `τ` (and not `t`): +We can confirm that `Xᵢ` and `Xₐ` depend on `τ` (and not `t`): ```@example dsl_advanced_ivs species(rn) ``` @@ -409,19 +412,19 @@ 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 + @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 + @ivs t x + @species Xᵢ(t,x) Xₐ(t,x) + (ka,kD), Xᵢ <--> Xₐ end species(rn) ``` @@ -435,8 +438,8 @@ It is possible to supply reactions with *metadata*, containing some additional i ```@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"] + p, 0 --> X, [description="Production reaction"] + d, X --> 0, [description="Degradation reaction"] end nothing # hide ``` @@ -444,7 +447,7 @@ 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"]) + (p,d), 0 <--> X, ([description="Production reaction"], [description="Degradation reaction"]) end nothing # hide ``` @@ -452,8 +455,8 @@ 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"] + p, 0 --> X, [description="Production reaction", misc=:value] + d, X --> 0, [description="Degradation reaction"] end nothing # hide ``` diff --git a/docs/src/model_creation/dsl_basics.md b/docs/src/model_creation/dsl_basics.md index 19057d2d65..8f8029585f 100644 --- a/docs/src/model_creation/dsl_basics.md +++ b/docs/src/model_creation/dsl_basics.md @@ -1,7 +1,7 @@ # [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. +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 simulation_intro) 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 @@ -9,21 +9,22 @@ 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 +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 basic_CRN_library_mm) can be written as ```@example dsl_basics_intro rn = @reaction_network begin - (kB,kD), S + E <--> SE - kP, SE --> P + E + (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). +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 simulation_intro). ## [Basic syntax](@id dsl_description_basic_syntax) The basic syntax of the DSL is ```@example dsl_basics +using Catalyst # hide rn = @reaction_network begin - 2.0, X --> Y - 1.0, Y --> X + 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 @@ -36,23 +37,22 @@ Each reaction line declares, in order, the rate, the substrate(s), and the produ 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: +Typically, the rates are not constants, but rather parameters (which values can be set e.g. at [the beginning of each simulation](@ref simulation_intro_ODEs)). 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 + 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 + 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. +Generally, anything that is a [permitted Julia variable name](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) @@ -60,14 +60,14 @@ Generally, anything that is a [permitted Julia variable name](@id https://docs.j 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 + 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 + k, X + Y + Z --> A + B + C + D end ``` @@ -75,8 +75,8 @@ end 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 + p, 0 --> X + d, X --> 0 end ``` @@ -84,19 +84,18 @@ end 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 + 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). +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 dsl_description_stoichiometries). 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 + 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) @@ -105,14 +104,14 @@ nothing # hide As is the case for the following two-state model: ```@example dsl_basics rn7 = @reaction_network begin - k1, X1 --> X2 - k2, X2 --> X1 + 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 + (k1,k2), X1 <--> X2 end ``` Here, the first rate (`k1`) denotes the *forward rate* and the second rate (`k2`) the *backwards rate*. @@ -120,13 +119,13 @@ Here, the first rate (`k1`) denotes the *forward rate* and the second rate (`k2` 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 + 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 + k, Y --> X end ``` Generally, using forward reactions is clearer than backwards ones, with the latter typically never being used. @@ -135,49 +134,49 @@ Generally, using forward reactions is clearer than backwards ones, with the latt 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 + 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 + 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 + dX, X --> 0 + dY, Y --> 0 end ``` This can be represented using: ```@example dsl_basics rn9 = @reaction_network begin - (dX,dY), (X,Y) --> 0 + (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) + 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 + 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) + (kf, (kb_1, kb_2)), S <--> (P1,P2) end ``` @@ -199,7 +198,7 @@ rn_13 = @reaction_network begin 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): +Here, `P`'s production rate will be reduced as `A` decays. We can [print the ODE this model produces with `Latexify`](@ref visualisation_latex): ```@example dsl_basics using Latexify latexify(rn_13; form=:ode) @@ -221,14 +220,14 @@ Here, while these models will generate identical ODE, SDE, and jump simulations, 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). + 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). 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 → ∅ + 2.0 + X^2, 0 --> X + Y + k1 + k2^k3, X --> ∅ + pi * X/(sqrt(2) + Y), Y → ∅ end ``` @@ -266,7 +265,7 @@ Catalyst comes with the following predefined functions: - 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: +Previously we have assumed that the rates are independent of the time variable, $t$. 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 @@ -319,8 +318,8 @@ Julia permits any Unicode characters to be used in variable names, thus Catalyst 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 --> ∅ + p, ∅ --> X + d, X --> ∅ end ``` @@ -333,8 +332,8 @@ Catalyst uses `-->`, `<-->`, and `<--` to denote forward, bi-directional, and ba E.g. the production/degradation system can alternatively be written as: ```@example dsl_basics rn4 = @reaction_network begin - p, ∅ → X - d, X → ∅ + p, ∅ → X + d, X → ∅ end ``` @@ -348,18 +347,19 @@ A range of possible characters are available which can be incorporated into spec 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σᵛ → σᵛ + 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: +```@example dsl_basics 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)). @@ -369,4 +369,3 @@ It should be noted that the following symbols are *not permitted* to be used as - `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/examples/basic_CRN_library.md b/docs/src/model_creation/examples/basic_CRN_library.md index 4955aa522d..4955bae3a6 100644 --- a/docs/src/model_creation/examples/basic_CRN_library.md +++ b/docs/src/model_creation/examples/basic_CRN_library.md @@ -86,9 +86,9 @@ Catalyst has special methods for working with conserved quantities, which are de ```@example crn_library_michaelis_menten using Catalyst mm_system = @reaction_network begin - kB, S + E --> SE - kD, SE --> S + E - kP, SE --> P + E + kB, S + E --> SE + kD, SE --> S + E + kP, SE --> P + E end ``` Next, we perform ODE, SDE, and jump simulations of the model: @@ -110,7 +110,7 @@ using JumpProcesses dprob = DiscreteProblem(mm_system, u0, tspan, ps) jprob = JumpProblem(mm_system, dprob, Direct()) jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(); seed = 12) # hide +jsol = solve(jprob, SSAStepper(), seed = 12) # hide using Plots oplt = plot(osol; title = "Reaction rate equation (ODE)") @@ -118,8 +118,10 @@ splt = plot(ssol; title = "Chemical Langevin equation (SDE)") jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)") plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1)) plot!(bottom_margin = 3Plots.Measures.mm) # hide +nothing # hide ``` -Note that, due to the large amounts of the species involved, teh stochastic trajectories are very similar to the deterministic one. +![MM Kinetics](../../assets/long_ploting_times/model_creation/mm_kinetics.svg) +Note that, due to the large amounts of the species involved, the stochastic trajectories are very similar to the deterministic one. ## [SIR infection model](@id basic_CRN_library_sir) 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. @@ -146,20 +148,22 @@ Next, we perform 3 different Jump simulations. Note that for the stochastic mode ```@example crn_library_sir using JumpProcesses dprob = DiscreteProblem(sir_model, u0, tspan, ps) -jprob = JumpProblem(sir_model, dprob, Direct()) +jprob = JumpProblem(sir_model, dprob, Direct(); save_positions = (false, false)) -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 +jsol1 = solve(jprob, SSAStepper(); saveat = 1.0) +jsol2 = solve(jprob, SSAStepper(); saveat = 1.0) +jsol3 = solve(jprob, SSAStepper(); saveat = 1.0) +jsol1 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 1) # hide +jsol2 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 2) # hide +jsol3 = solve(jprob, SSAStepper(); saveat = 1.0, seed = 3) # hide jplt1 = plot(jsol1; title = "Outbreak") jplt2 = plot(jsol2; title = "Outbreak") jplt3 = plot(jsol3; title = "No outbreak") plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1)) +nothing # hide ``` +![SIR Outbreak](../../assets/long_ploting_times/model_creation/sir_outbreaks.svg) ## [Chemical cross-coupling](@id basic_CRN_library_cc) In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the [Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm). @@ -241,9 +245,9 @@ oprob = ODEProblem(sa_loop, u0, tspan, ps) osol = solve(oprob) dprob = DiscreteProblem(sa_loop, u0, tspan, ps) -jprob = JumpProblem(sa_loop, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(); seed = 12) # hide +jprob = JumpProblem(sa_loop, dprob, Direct(); save_positions = (false,false)) +jsol = solve(jprob, SSAStepper(); saveat = 10.0) +jsol = solve(jprob, SSAStepper(); saveat = 10.0, seed = 12) # hide plot(osol; lw = 3, label = "Reaction rate equation (ODE)") plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350)) @@ -279,7 +283,7 @@ oplt2 = plot(osol2; title = "Oscillation") plot(oplt1, oplt2; lw = 3, size = (800,600), layout = (2,1)) ``` -## [The Repressilator](@id basic_CRN_library_) +## [The Repressilator](@id basic_CRN_library_repressilator) The Repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can 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$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)). ```@example crn_library_brusselator using Catalyst diff --git a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md index ddd7993ed0..ae2be87b48 100644 --- a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md +++ b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md @@ -1,5 +1,5 @@ # [Programmatic, generative, modelling of a linear pathway](@id programmatic_generative_linear_pathway) -This example will show how to use programmatic, generative, modelling to model a system implicitly. I.e. rather than listing all system reactions explicitly, the reactions are implicitly generated from a simple set of rules. This example is specifically designed to show how [programmatic modelling](@ref ref) enables *generative workflows* (demonstrating one of its advantages as compared to [DSL-based modelling](@ref ref)). In our example, we will model linear pathways, so we will first introduce these. Next, we will model them first using the DSL, and then using a generative programmatic workflow. +This example will show how to use programmatic, generative, modelling to model a system implicitly. I.e. rather than listing all system reactions explicitly, the reactions are implicitly generated from a simple set of rules. This example is specifically designed to show how [programmatic modelling](@ref programmatic_CRN_construction) enables *generative workflows* (demonstrating one of its advantages as compared to [DSL-based modelling](@ref dsl_description)). In our example, we will model linear pathways, so we will first introduce these. Next, we will model them first using the DSL, and then using a generative programmatic workflow. ## [Linear pathways](@id programmatic_generative_linear_pathway_intro) Linear pathways consists of a series of species ($X_0$, $X_1$, $X_2$, ..., $X_n$) where each activates the subsequent one. These are often modelled through the following reaction system: @@ -75,11 +75,13 @@ plot!(sol_n10; idxs = :X10, label = "n = 10") ``` ## [Modelling linear pathways using programmatic, generative, modelling](@id programmatic_generative_linear_pathway_generative) -Above, we investigated the impact of linear pathways' lengths on their behaviours. Since the models were implemented using the DSL, we had to implement a new model for each pathway (in each case writing out all reactions). Here, we will instead show how [programmatic modelling](@ref ref) can be used to generate pathways of arbitrary lengths. +Above, we investigated the impact of linear pathways' lengths on their behaviours. Since the models were implemented using the DSL, we had to implement a new model for each pathway (in each case writing out all reactions). Here, we will instead show how [programmatic modelling](@ref programmatic_CRN_construction) can be used to generate pathways of arbitrary lengths. First, we create a function, `generate_lp`, which creates a linear pathway model of length `n`. It utilises [*vector variables*](@ref ref) to create an arbitrary number of species, and also creates an [observable](@ref ref) for the final species of the chain. ```@example programmatic_generative_linear_pathway_generative using Catalyst # hide +t = default_t() +@parameters τ function generate_lp(n) # Creates a vector `X` with n+1 species. @species X(t)[1:n+1] @@ -115,6 +117,7 @@ nothing # hide ``` We can now simulate linear pathways of arbitrary lengths using a simple syntax. We use this to recreate our previous result from the DSL: ```@example programmatic_generative_linear_pathway_generative +using OrdinaryDiffEq, Plots # hide sol_n3 = solve(generate_oprob(3)) sol_n10 = solve(generate_oprob(10)) plot(sol_n3; idxs = :Xend, label = "n = 3") diff --git a/docs/src/model_creation/model_visualisation.md b/docs/src/model_creation/model_visualisation.md index 463c96ce5c..6efb71b4bb 100644 --- a/docs/src/model_creation/model_visualisation.md +++ b/docs/src/model_creation/model_visualisation.md @@ -4,7 +4,7 @@ Catalyst-created `ReactionSystem` models can be visualised either as LaTeX code ## [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): +Let us consider a simple [Brusselator model](@ref basic_CRN_library_brusselator): ```@example visualisation_latex using Catalyst brusselator = @reaction_network begin @@ -32,10 +32,10 @@ latexify(brusselator; form = :ode) 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. + 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](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: +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 basic_CRN_library_brusselator) and then displays its network topology: ```@example visualisation_graphs using Catalyst brusselator = @reaction_network begin @@ -46,7 +46,7 @@ brusselator = @reaction_network begin 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): +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 basic_CRN_library_repressilator): ```@example visualisation_graphs repressilator = @reaction_network begin hillr(Z,v,K,n), ∅ --> X @@ -64,7 +64,7 @@ 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: +Finally, a [network's reaction complexes](@ref network_analysis_reaction_complexes) (and the reactions in between these) can be displayed using the `complexgraph(brusselator)` function: ```@example visualisation_graphs complexgraph(brusselator) ``` diff --git a/docs/src/model_creation/network_analysis.md b/docs/src/model_creation/network_analysis.md index 0f212e2447..792ed5477d 100644 --- a/docs/src/model_creation/network_analysis.md +++ b/docs/src/model_creation/network_analysis.md @@ -101,7 +101,7 @@ Let's check these are equal symbolically isequal(odes, odes2) ``` -## Reaction complex representation +## [Reaction complex representation](@id network_analysis_reaction_complexes) We now introduce a further decomposition of the RRE ODEs, which has been used to facilitate analysis of a variety of reaction network properties. Consider a simple reaction system like diff --git a/docs/src/model_simulation/ensemble_simulations.md b/docs/src/model_simulation/ensemble_simulations.md index 47689b7050..79a8dce1ec 100644 --- a/docs/src/model_simulation/ensemble_simulations.md +++ b/docs/src/model_simulation/ensemble_simulations.md @@ -3,10 +3,10 @@ In many contexts, a single model is re-simulated under similar conditions. Examp - Performing Monte Carlo simulations of a stochastic model to gain insight in its behaviour. - Scanning a model's behaviour for different parameter values and/or initial conditions. -While this can be handled using `for` loops, it is typically better to first create an `EnsembleProblem`, and then perform an ensemble simulation. Advantages include a more concise interface and the option for [automatic simulation parallelisation](@ref ref). Here we provide a short tutorial on how to perform parallel ensemble simulations, with a more extensive documentation being available [here](@ref https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/). +While this can be handled using `for` loops, it is typically better to first create an `EnsembleProblem`, and then perform an ensemble simulation. Advantages include a more concise interface and the option for [automatic simulation parallelisation](@ref ode_simulation_performance_parallelisation). Here we provide a short tutorial on how to perform parallel ensemble simulations, with a more extensive documentation being available [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/). ## [Monte Carlo simulations using unmodified conditions](@id ensemble_simulations_monte_carlo) -We will first consider Monte Carlo simulations where the simulation conditions are identical in-between simulations. First, we declare a [simple self-activation loop](@ref ref) model +We will first consider Monte Carlo simulations where the simulation conditions are identical in-between simulations. First, we declare a [simple self-activation loop](@ref basic_CRN_library_self_activation) model ```@example ensemble using Catalyst sa_model = @reaction_network begin @@ -16,23 +16,28 @@ end u0 = [:X => 10.0] tspan = (0.0, 1000.0) ps = [:v0 => 0.1, :v => 2.5, :K => 40.0, :n => 4.0, :deg => 0.01] +nothing # hide ``` We wish to simulate it as an SDE. Rather than performing a single simulation, however, we want to perform multiple ones. Here, we first create a normal `SDEProblem`, and use it as the single input to a `EnsembleProblem` (`EnsembleProblem` are created similarly for ODE and jump simulations, but the `ODEProblem` or `JumpProblem` is used instead). ```@example ensemble +using StochasticDiffEq sprob = SDEProblem(sa_model, u0, tspan, ps) eprob = EnsembleProblem(sprob) nothing # hide ``` -Next, the `EnsembleProblem` can be used as input to the `solve` command. Here, we use exactly the same inputs that we use for single simulations, however, we add a `trajectories` argument to denote how many simulations we wish to carry out. Here we perform 100 simulations: +Next, the `EnsembleProblem` can be used as input to the `solve` command. Here, we use exactly the same inputs that we use for single simulations, however, we add a `trajectories` argument to denote how many simulations we wish to carry out. Here we perform 10 simulations: ```@example ensemble -sols = solve(eprob, STrapezoid(); trajectories = 100) +sols = solve(eprob, STrapezoid(); trajectories = 10) nothing # hide ``` Finally, we can use our ensemble simulation solution as input to `plot` (just like normal simulations): ```@example ensemble -plot(sols; la = 0.5) +using Plots +plot(sols) ``` -Here, each simulation is displayed as an individual trajectory. We also use the [`la` plotting option](@ref ref) to reduce the transparency of each individual line, improving the plot visual. +Here, each simulation is displayed as an individual trajectory. +!!! note + While not used here, the [`la` plotting option](@ref simulation_plotting_options) (which modifies line transparency) can help improve the plot visual when a large number of (overlapping) lines are plotted. Various convenience functions are available for analysing and plotting ensemble simulations (a full list can be found [here]). Here, we use these to first create an `EnsembleSummary` (retrieving each simulation's value at time points `0.0, 1.0, 2.0, ... 1000.0`). Next, we use this as an input to the `plot` command, which automatically plots the mean $X$ activity across the ensemble, while also displaying the 5% and 95% quantiles as the shaded area: ```@example ensemble @@ -45,7 +50,8 @@ Previously, we assumed that each simulation used the same initial conditions and Here, we first create an `ODEProblem` of our previous self-activation loop: ```@example ensemble -oprob = ODEProblem(sa_model, u0, tspan, p) +using OrdinaryDiffEq +oprob = ODEProblem(sa_model, u0, tspan, ps) nothing # hide ``` Next, we wish to simulate the model for a range of initial conditions of $X$`. To do this we create a problem function, which takes the following arguments: @@ -53,7 +59,7 @@ Next, we wish to simulate the model for a range of initial conditions of $X$`. T - `i`: The number of this specific Monte Carlo iteration in the interval `1:trajectories`. - `repeat`: The iteration of the repeat of the simulation. Typically `1`, but potentially higher if [the simulation re-running option](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem) is used. -Here we will use the following problem function (utilising [remake](@ref ref)), which will provide a uniform range of initial concentrations of $X$: +Here we will use the following problem function (utilising [remake](@ref simulation_structure_interfacing_problems_remake)), which will provide a uniform range of initial concentrations of $X$: ```@example ensemble function prob_func(prob, i, repeat) remake(prob; u0 = [:X => i * 5.0]) diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md index 78c0dfe8a6..83e5fc7061 100644 --- a/docs/src/model_simulation/ode_simulation_performance.md +++ b/docs/src/model_simulation/ode_simulation_performance.md @@ -13,7 +13,7 @@ Generally, this short checklist provides a quick guide for dealing with ODE perf ## [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): +Here we simulate the (stiff) [Brusselator](@ref basic_CRN_library_brusselator) model using the `Tsit5` solver (which is designed for non-stiff ODEs): ```@example ode_simulation_performance_1 using Catalyst, OrdinaryDiffEq, Plots @@ -26,13 +26,16 @@ end u0 = [:X => 1.0, :Y => 0.0] tspan = (0.0, 20.0) -ps = [:A => 10.0, :B => 40.0] +ps = [:A => 400.0, :B => 2000.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: +nothing # hide +``` +![Incomplete Brusselator Simulation](../assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg) + +We get a warning, indicating that the simulation was terminated. Furthermore, the resulting plot ends at $t ≈ 12$, 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 ``` @@ -52,7 +55,7 @@ Finally, we should note that stiffness is not tied to the model equations only. ## [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): +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 basic_CRN_library_bd): ```@example ode_simulation_performance_2 using Catalyst, OrdinaryDiffEq @@ -75,18 +78,18 @@ 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 +```julia 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), +If you perform the above benchmarks on your machine, and check the results, you will 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 +```julia @btime solve(oprob, Tsit5(); abstol=1e-6, reltol=1e-6) @btime solve(oprob, Tsit5(); abstol=1e-12, reltol=1e-12) ``` @@ -123,7 +126,7 @@ 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) +oprob = ODEProblem(brusselator, u0, tspan, ps; sparse = true) nothing # hide ``` @@ -143,7 +146,7 @@ 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) +### [Designating preconditioners for Jacobian-free linear solvers](@id 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): @@ -172,10 +175,10 @@ Generally, the use of preconditioners is only recommended for advanced users who ## [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/). +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. `EnsembleProblem`s have [previously been introduced in Catalyst's documentation](@ref ensemble_simulations) (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$): +For this example (and the one for GPUs), we will consider a modified [Michaelis-Menten enzyme kinetics model](@ref basic_CRN_library_mm), 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 @@ -186,12 +189,12 @@ mm_model = @reaction_network begin end ``` The model can be simulated, showing how $P$ is produced from $S$: -```@example ode_simulation_performance_3 +```@example ode_simulation_performance_4 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) +ps = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :d => 0.1] +oprob = ODEProblem(mm_model, u0, tspan, ps) sol = solve(oprob, Tsit5()) plot(sol) ``` @@ -208,7 +211,7 @@ Here, `prob_func` takes 3 arguments: 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): +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 simulation_structure_interfacing_problems_remake): ```@example ode_simulation_performance_4 function prob_func(prob, i, repeat) return remake(prob; p = [:kP => 0.01*i]) @@ -232,22 +235,21 @@ 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="") +plot!(left_margin = 3Plots.Measures.mm) # hide ``` 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 +```julia +esol = solve(eprob, Tsit5(), EnsembleDistributed(); trajectories = 100) ``` 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 +```julia 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 +```julia 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. @@ -268,7 +270,7 @@ Furthermore (while not required) to receive good performance, we should also mak - 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 +```julia 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/). @@ -283,11 +285,12 @@ mm_model = @reaction_network begin kP, SE --> P + E d, S --> ∅ end +@unpack S, E, SE, P, kB, kD, kP, d = mm_model using OrdinaryDiffEq, Plots -u0 = @SVector [:S => 1.0f0, :E => 1.0f0, :SE => 0.0f0, :P => 0.0f0] +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] +p = @SVector [kB => 1.0f0, kD => 0.1f0, kP => 0.5f0, d => 0.1f0] oprob = ODEProblem(mm_model, u0, tspan, p) nothing # hide ``` @@ -300,16 +303,17 @@ 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). +!!! note + Currently, declaration of static vectors requires [symbolic, rather than symbol, form](@ref ref) for species and parameters. Hence, we here first [`@unpack` these](@ref ref) before constructing `u0` and `ps` using `@SVector`. 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 +```julia 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). diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md index 113b0fa7b4..81894ac74e 100644 --- a/docs/src/model_simulation/simulation_introduction.md +++ b/docs/src/model_simulation/simulation_introduction.md @@ -1,7 +1,7 @@ # [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). +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 simulation_plotting), and on how to [interact with simulation problems, integrators, and solutions](@ref simulation_structure_interfacing). Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for [ODEs](@ref ode_simulation_performance), [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. @@ -40,9 +40,9 @@ These three different approaches are summed up in the following table: 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/) + Euler, Runge-Kutta + Euler-Maruyama, Milstein + Gillespie, Sorting direct Species units @@ -70,22 +70,22 @@ These three different approaches are summed up in the following table: 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) + OrdinaryDiffEq.jl + StochasticDiffEq.jl + 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) +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 things specific to [SDE](@ref simulation_intro_SDEs) and [jump](@ref simulation_intro_jumps) 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 ode_simulation_performance) -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): +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 basic_CRN_library_two_states): ```@example simulation_intro_ode using Catalyst two_state_model = @reaction_network begin - (k1,k2), X1 <--> X2 + (k1,k2), X1 <--> X2 end u0 = [:X1 => 100.0, :X2 => 200.0] tspan = (0.0, 5.0) @@ -108,7 +108,7 @@ Finally, the result can be plotted using the [Plots.jl](https://github.com/Julia 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). +More information on how to interact with solution structures is provided [here](@ref simulation_structure_interfacing) and on how to plot them [here](@ref simulation_plotting). Some additional considerations: - If a model without parameters has been declared, only the first three arguments must be provided to `ODEProblem`. @@ -122,11 +122,11 @@ While good defaults are generally selected, OrdinaryDiffEq enables the user to c 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). +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 ode_simulation_performance_solvers). -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 +Additional options can be provided as keyword arguments. E.g. the `maxiters` arguments determines the maximum number of simulation time steps (before the simulation is terminated). This defaults to `1e5`, but can be modified through: ```@example simulation_intro_ode -sol = solve(oprob; adaptive = false) +sol = solve(oprob; maxiters = 1e4) nothing # hide ``` @@ -160,7 +160,7 @@ 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. + It is possible to [designate specific types for parameters](@ref dsl_advanced_options_parameter_types). 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 @@ -174,9 +174,9 @@ Catalyst uses the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffE 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 +using Catalyst, StochasticDiffEq, Plots two_state_model = @reaction_network begin - (k1,k2), X1 <--> X2 + (k1,k2), X1 <--> X2 end u0 = [:X1 => 100.0, :X2 => 200.0] tspan = (0.0, 1.0) @@ -193,23 +193,30 @@ we can see that while this simulation (unlike the ODE ones) exhibits some fluctu 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: +Next, let us reduce species amounts (using [`remake`](@ref simulation_structure_interfacing_problems_remake)), 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 +nothing # hide +``` +Here, we receive a warning that the simulation was terminated. next, if we plot the solution: +```@example simulation_intro_sde 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. +we note that the simulation didn't reach the designated final time point ($t = 1.0$). 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()) +nothing # hide +``` +```@example simulation_intro_sde 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): +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 ode_simulation_performance_error): ```@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 @@ -220,9 +227,10 @@ plot(sol) 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 + @default_noise_scaling 0.1 + (k1,k2), X1 <--> X2 end +nothing # hide ``` 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 @@ -238,10 +246,11 @@ 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 + @parameters η + @default_noise_scaling η + (k1,k2), X1 <--> X2 end +nothing # hide ``` 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 @@ -255,19 +264,19 @@ 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)). + 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 dsl_advanced_options_declaring_species_and_parameters)). -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: +It is possible to designate specific noise scaling terms for individual reactions through the `noise_scaling` [reaction metadata](@ref dsl_advanced_options_reaction_metadata). 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 + 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. +While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the [`remake_reactionsystem`](@ref) function can be used to achieve a similar effect. ## [Performing jump simulations using stochastic chemical kinetics](@id simulation_intro_jumps) @@ -275,9 +284,9 @@ Catalyst uses the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) 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 +using Catalyst, JumpProcesses, Plots two_state_model = @reaction_network begin - (k1,k2), X1 <--> X2 + (k1,k2), X1 <--> X2 end u0 = [:X1 => 5, :X2 => 10] tspan = (0.0, 5.0) @@ -292,37 +301,33 @@ Next, we bundle these into a `DiscreteProblem` (similarly to how `ODEProblem`s a 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. +This is then used as input to a `JumpProblem`. The `JumpProblem` also requires the CRN model and [an aggregator](@ref simulation_intro_jumps_solver_designation) 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 +using Plots 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: +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 the sorting direct method (`SortingDirect`) should be used, use: ```@example simulation_intro_jumps -jprob = JumpProblem(brusselator, dprob, Direct()) +jprob = JumpProblem(two_state_model, dprob, SortingDirect()) 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 -``` +Next, a simulation method can be provided (like for ODEs and SDEs) as the second argument to `solve`. Currently, the only relevant solver is `SSAStepper()` (which is the one used throughout Catalyst's documentation). Other choices are primarily relevant to combined ODE/SDE + jump simulations, or inexact simulations. These situations are described in more detail [here](https://docs.sciml.ai/JumpProcesses/stable/jump_solve/). ### [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: diff --git a/docs/src/model_simulation/simulation_plotting.md b/docs/src/model_simulation/simulation_plotting.md index 25cbec54ee..b1e471cd04 100644 --- a/docs/src/model_simulation/simulation_plotting.md +++ b/docs/src/model_simulation/simulation_plotting.md @@ -1,11 +1,11 @@ -# [Simulation_plotting](@id simulation_plotting) +# [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: +Let us consider the oscillating [Brusselator](@ref basic_CRN_library_brusselator) 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: ```@example simulation_plotting using Catalyst, OrdinaryDiffEq, Plots diff --git a/docs/src/model_simulation/simulation_structure_interfacing.md b/docs/src/model_simulation/simulation_structure_interfacing.md index a993841262..aa70d51ed4 100644 --- a/docs/src/model_simulation/simulation_structure_interfacing.md +++ b/docs/src/model_simulation/simulation_structure_interfacing.md @@ -5,7 +5,7 @@ Generally, when we have a structure `simulation_struct` and want to interface wi ## [Interfacing problem objects](@id simulation_structure_interfacing_problems) -We begin by demonstrating how we can interface with problem objects. First, we create an `ODEProblem` representation of a [chemical cross-coupling model](@ref ref) (where a catalyst, $C$, couples two substrates, $S₁$ and $S₂$, to form a product, $P$). +We begin by demonstrating how we can interface with problem objects. First, we create an `ODEProblem` representation of a [chemical cross-coupling model](@ref basic_CRN_library_cc) (where a catalyst, $C$, couples two substrates, $S₁$ and $S₂$, to form a product, $P$). ```@example structure_indexing using Catalyst cc_system = @reaction_network begin @@ -21,7 +21,7 @@ oprob = ODEProblem(cc_system, u0, tspan, ps) nothing # hide ``` -We can find a specie's (or [variable's](@ref ref)) initial condition value by simply indexing with the species of interest as input. Here we check the initial condition value of $C$: +We can find a species's (or [variable's](@ref ref)) initial condition value by simply indexing with the species of interest as input. Here we check the initial condition value of $C$: ```@example structure_indexing oprob[:C] ``` @@ -53,7 +53,7 @@ oprob[[:S₁, :S₂]] Generally, when updating problems, it is often better to use the [`remake` function](@ref simulation_structure_interfacing_problems_remake) (especially when several values are updated). !!! warn - Indexing *should not* be used not modify `JumpProblem`s. Here, [remake](@ref ref) should be used exclusively. + Indexing *should not* be used not modify `JumpProblem`s. Here, [remake](@ref simulation_structure_interfacing_problems_remake) should be used exclusively. A problem's time span can be accessed through the `tspan` field: ```@example structure_indexing @@ -162,7 +162,7 @@ get_S(oprob) ``` ## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation) -As [previously described](@ref ref), when e.g. [programmatic modelling is used](@ref ref), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref ref) programmatically, and use its symbolic variables to check, and update, an initial condition: +As [previously described](@ref ref), when e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref rbasic_CRN_library_two_statesef) programmatically, and use its symbolic variables to check, and update, an initial condition: ```@example structure_indexing_symbolic_variables using Catalyst t = default_t() @@ -195,7 +195,7 @@ Just like symbolic variables can be used to directly interface with a structure, oprob[two_state_model.X1 + two_state_model.X2] ``` This can be used to form symbolic expressions using model quantities when a model has been created using the DSL (as an alternative to [@unpack] -(@ref ref)). Alternatively, [creating an observable](@ref ref), and then interface using its `Symbol` representation, is also possible. +(@ref ref)). Alternatively, [creating an observable](@ref dsl_advanced_options_observables), and then interface using its `Symbol` representation, is also possible. !!! warn With interfacing with a simulating structure using symbolic variables stored in a `ReactionSystem` model, ensure that the [model is complete](@ref ref). \ No newline at end of file diff --git a/docs/src/steady_state_functionality/bifurcation_diagrams.md b/docs/src/steady_state_functionality/bifurcation_diagrams.md index d869141e03..cebe51dc1c 100644 --- a/docs/src/steady_state_functionality/bifurcation_diagrams.md +++ b/docs/src/steady_state_functionality/bifurcation_diagrams.md @@ -37,21 +37,21 @@ nothing # hide BifurcationKit computes bifurcation diagrams using the `bifurcationdiagram` function. From an initial point in the diagram, it tracks the solution (using a continuation algorithm) until the entire diagram is computed (BifurcationKit's continuation can be used for other purposes, however, this tutorial focuses on bifurcation diagram computation). The continuation settings are provided in a `ContinuationPar` structure. In this example, we will only specify three settings, `p_min` and `p_max` (which sets the minimum and maximum values over which the bifurcation parameter is varied) and `max_steps` (the maximum number of continuation steps to take as the bifurcation diagram is tracked). We wish to compute a bifurcation diagram over the interval $(2.0,20.0)$, and will use the following settings: ```@example ex1 p_span = (2.0, 20.0) -opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps=1000) +opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000) nothing # hide ``` Finally, we compute our bifurcation diagram using: ```@example ex1 -bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true) +bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) nothing # hide ``` -Where `PALC()` designates that we wish to use the pseudo arclength continuation method to track our solution. The third argument (`2`) designates the maximum number of recursions when branches of branches are computed (branches appear as continuation encounters certain bifurcation points). For diagrams with highly branched structures (rare for many common small chemical reaction networks) this input is important. Finally, `bothside=true` designates that we wish to perform continuation on both sides of the initial point (which is typically the case). +Where `PALC()` designates that we wish to use the pseudo arclength continuation method to track our solution. The third argument (`2`) designates the maximum number of recursions when branches of branches are computed (branches appear as continuation encounters certain bifurcation points). For diagrams with highly branched structures (rare for many common small chemical reaction networks) this input is important. Finally, `bothside = true` designates that we wish to perform continuation on both sides of the initial point (which is typically the case). We can plot our bifurcation diagram using the Plots.jl package: ```@example ex1 using Plots -plot(bif_dia; xguide="k1", yguide="X") +plot(bif_dia; xguide = "k1", yguide = "X") ``` Here, the steady state concentration of $X$ is shown as a function of $k1$'s value. Stable steady states are shown with thick lines, unstable ones with thin lines. The two [fold bifurcation points](https://en.wikipedia.org/wiki/Saddle-node_bifurcation) are marked with "bp". @@ -69,7 +69,7 @@ opt_newton = NewtonPar(tol = 1e-9, max_iterations = 1000) opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], dsmin = 0.001, dsmax = 0.01, max_steps = 1000, newton_options = opt_newton) -bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true) +bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) nothing # hide ``` (however, in this case these additional settings have no significant effect on the result) @@ -79,8 +79,8 @@ Let's consider the previous case, but instead compute the bifurcation diagram ov ```@example ex1 p_span = (2.0, 15.0) opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000) -bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true) -plot(bif_dia; xguide="k1", yguide="X") +bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) +plot(bif_dia; xguide = "k1", yguide = "X") ``` Here, in the bistable region, we only see a single branch. The reason is that the continuation algorithm starts at our initial guess (here made at $k1 = 4.0$ for $(X,Y) = (5.0,2.0)$) and tracks the diagram from there. However, with the upper bound set at $k1=15.0$ the bifurcation diagram has a disjoint branch structure, preventing the full diagram from being computed by continuation alone. In this case it could be solved by increasing the bound from $k1=15.0$, however, this is not possible in all cases. In these cases, *deflation* can be used. This is described in the [BifurcationKit documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2/#Snaking-computed-with-deflation). @@ -99,12 +99,12 @@ end u_guess = [:K => 1.0, :X => 1.0, :Xp => 1.0] p_start = [:p => 1.0, :d => 0.5, :kP => 2.0, :kD => 5.0] u0 = [:X => 1.0, :Xp => 0.0] -bprob = BifurcationProblem(kinase_model, u_guess, p_start, :d; plot_var=:Xp, u0=u0) +bprob = BifurcationProblem(kinase_model, u_guess, p_start, :d; plot_var = :Xp, u0) p_span = (0.1, 10.0) opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000) -bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true) -plot(bif_dia; xguide="d", yguide="Xp") +bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) +plot(bif_dia; xguide = "d", yguide = "Xp") ``` This bifurcation diagram does not contain any interesting features (such as bifurcation points), and only shows how the steady state concentration of $Xp$ is reduced as $d$ increases. diff --git a/docs/src/steady_state_functionality/dynamical_systems.md b/docs/src/steady_state_functionality/dynamical_systems.md index fcd56129e0..45847a12f8 100644 --- a/docs/src/steady_state_functionality/dynamical_systems.md +++ b/docs/src/steady_state_functionality/dynamical_systems.md @@ -54,7 +54,7 @@ More information on how to compute basins of attractions for ODEs using Dynamica While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviour if its attractor(s) have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not. -First, let us consider the [Willamowski–Rössler model](@ref ref), which is known to exhibit chaotic behaviour. +First, let us consider the [Willamowski–Rössler model](@ref basic_CRN_library_wr), which is known to exhibit chaotic behaviour. ```@example dynamical_systems_lyapunov using Catalyst wr_model = @reaction_network begin @@ -89,7 +89,7 @@ Here, the `autodiff = false` argument is required when Lyapunov spectrums are co ```@example dynamical_systems_lyapunov lyapunovspectrum(ds, 100) ``` -Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one chaotic attractor, to which is approached from the initial condition $(1.5,1.5,1.5)). +Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one chaotic attractor, to which is approached from the initial condition $(1.5,1.5,1.5)$). Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle: ```@example dynamical_systems_lyapunov diff --git a/docs/src/steady_state_functionality/homotopy_continuation.md b/docs/src/steady_state_functionality/homotopy_continuation.md index 19a15d39e6..ee81f88121 100644 --- a/docs/src/steady_state_functionality/homotopy_continuation.md +++ b/docs/src/steady_state_functionality/homotopy_continuation.md @@ -12,7 +12,7 @@ integer Hill exponents). The roots of these can reliably be found through a 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. -For this tutorial, we will use the [Wilhelm model](@ref ref) (which +For this tutorial, we will use the [Wilhelm model](@ref basic_CRN_library_wilhelm) (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 hc_basics diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index 601fd51c93..3a97392fdf 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/docs/src/steady_state_functionality/nonlinear_solve.md @@ -8,9 +8,9 @@ While these approaches only find a single steady state, they offer two advantage - 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. -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/). +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 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/). -## [Steady state finding through nonlinear solving](@ref steady_state_solving_nonlinear) +## [Steady state finding through nonlinear solving](@id 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 @@ -43,6 +43,7 @@ To solve this problem, we must first designate our parameter values, and also ma 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] nlprob = NonlinearProblem(dimer_production, u_guess, p) +nothing # hide ``` Finally, we can solve it using the `solve` command, returning the steady state solution: ```@example steady_state_solving_nonlinear @@ -57,7 +58,7 @@ sol ≈ sol_ntr ``` ### [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): +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](@ref basic_CRN_library_two_states): ```@example steady_state_solving_claws using Catalyst, NonlinearSolve # hide two_state_model = @reaction_network begin @@ -85,6 +86,7 @@ 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 +using Catalyst # hide dimer_production = @reaction_network begin pₘ, 0 --> mRNA pₚ, mRNA --> mRNA + P @@ -98,21 +100,18 @@ nothing # hide Next, we provide these as an input to a `SteadyStateProblem` ```@example steady_state_solving_simulation ssprob = SteadyStateProblem(dimer_production, u0, p) +nothing # hide ``` -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. +Finally, we can find the steady states using the `solver` command. Since `SteadyStateProblem`s are solved through forward ODE simulation, we must load the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package, and [select an ODE solver](@ref simulation_intro_solver_options). 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: -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: +(which requires loading the SteadyStateDiffEq package). ```@example steady_state_solving_simulation -using OrdinaryDiffEqDiffEq +using SteadyStateDiffEq, OrdinaryDiffEq solve(ssprob, DynamicSS(Rodas5P())) -nothing # hide ``` -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: +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. + +Generally, `SteadyStateProblem`s can be solved using the [same options that are available for ODE simulations](@ref simulation_intro_solver_options). E.g. here we designate a specific `dt` step size: ```@example steady_state_solving_simulation solve(ssprob, DynamicSS(Rodas5P()); dt = 0.01) nothing # hide @@ -144,6 +143,7 @@ If you use this functionality in your research, [in addition to Catalyst](@ref c } ``` + --- ## References -[^1]: [J. Nocedal, S. J. Wright, *Numerical Optimization*, Springer (2006).](https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf) +[^1]: [J. Nocedal, S. J. Wright, *Numerical Optimization*, Springer (2006).](https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf) \ No newline at end of file diff --git a/docs/src/steady_state_functionality/steady_state_stability_computation.md b/docs/src/steady_state_functionality/steady_state_stability_computation.md index 0cc60ebf94..8e43dcfe28 100644 --- a/docs/src/steady_state_functionality/steady_state_stability_computation.md +++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md @@ -19,7 +19,7 @@ steady_state = [:X => 4.0] steady_state_stability(steady_state, rn, ps) ``` -Next, let us consider the following [self-activation loop](@ref ref): +Next, let us consider the following [self-activation loop](@ref basic_CRN_library_self_activation): ```@example stability_1 sa_loop = @reaction_network begin (hill(X,v,K,n),d), 0 <--> X @@ -51,11 +51,11 @@ ss_jac = steady_state_jac(sa_loop) ps_1 = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0] steady_states_1 = hc_steady_states(sa_loop, ps) -stability_1 = [steady_state_stability(state, sa_loop, ps_1; ss_jac = ss_jac) for state in steady_states_1] +stabs_1 = [steady_state_stability(st, sa_loop, ps_1; ss_jac) for st in steady_states_1] ps_2 = [:v => 4.0, :K => 1.5, :n => 2, :d => 1.0] steady_states_2 = hc_steady_states(sa_loop, ps) -stability_2 = [steady_state_stability(state, sa_loop, ps_2; ss_jac = ss_jac) for state in steady_states_2] +stabs_2 = [steady_state_stability(st, sa_loop, ps_2; ss_jac) for st in steady_states_2] nothing # hide ``` diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 2aab60fb54..42433bcd2f 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -93,7 +93,6 @@ end include("reaction.jl") export isspecies export Reaction -export get_noise_scaling, has_noise_scaling # The `ReactionSystem` structure and its functions. include("reactionsystem.jl") diff --git a/src/reaction.jl b/src/reaction.jl index 95eb7c0e71..23538214aa 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -377,8 +377,8 @@ function ModelingToolkit.get_variables!(set, rx::Reaction) for stoichs in (rx.substoich, rx.prodstoich), stoich in stoichs (stoich isa BasicSymbolic) && get_variables!(set, stoich) end - if has_noise_scaling(rx) - get_variables!(set, get_noise_scaling(rx)) + if hasnoisescaling(rx) + get_variables!(set, getnoisescaling(rx)) end return (set isa AbstractVector) ? unique!(set) : set end @@ -440,7 +440,7 @@ end ### Reaction Metadata Implementation ### -# These are currently considered internal, but can be used by public accessor functions like get_noise_scaling. +# These are currently considered internal, but can be used by public accessor functions like getnoisescaling. """ getmetadata_dict(reaction::Reaction) @@ -507,7 +507,7 @@ end # Noise scaling. """ -has_noise_scaling(reaction::Reaction) +hasnoisescaling(reaction::Reaction) Returns `true` if the input reaction has the `noise_scaing` metadata field assigned, else `false`. @@ -517,15 +517,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [noise_scaling=0.0] -has_noise_scaling(reaction) +hasnoisescaling(reaction) ``` """ -function has_noise_scaling(reaction::Reaction) +function hasnoisescaling(reaction::Reaction) return hasmetadata(reaction, :noise_scaling) end """ -get_noise_scaling(reaction::Reaction) +getnoisescaling(reaction::Reaction) Returns `noise_scaing` metadata field for the input reaction. @@ -535,11 +535,11 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [noise_scaling=0.0] -get_noise_scaling(reaction) +getnoisescaling(reaction) ``` """ -function get_noise_scaling(reaction::Reaction) - if has_noise_scaling(reaction) +function getnoisescaling(reaction::Reaction) + if hasnoisescaling(reaction) return getmetadata(reaction, :noise_scaling) else error("Attempts to access noise_scaling metadata field for a reaction which does not have a value assigned for this metadata.") @@ -548,7 +548,7 @@ end # Description. """ -has_description(reaction::Reaction) +hasdescription(reaction::Reaction) Returns `true` if the input reaction has the `description` metadata field assigned, else `false`. @@ -558,15 +558,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [description="A reaction"] -has_description(reaction) +hasdescription(reaction) ``` """ -function has_description(reaction::Reaction) +function hasdescription(reaction::Reaction) return hasmetadata(reaction, :description) end """ -get_description(reaction::Reaction) +getdescription(reaction::Reaction) Returns `description` metadata field for the input reaction. @@ -576,11 +576,11 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [description="A reaction"] -get_description(reaction) +getdescription(reaction) ``` """ -function get_description(reaction::Reaction) - if has_description(reaction) +function getdescription(reaction::Reaction) + if hasdescription(reaction) return getmetadata(reaction, :description) else error("Attempts to access `description` metadata field for a reaction which does not have a value assigned for this metadata.") @@ -589,7 +589,7 @@ end # Misc. """ -has_misc(reaction::Reaction) +hasmisc(reaction::Reaction) Returns `true` if the input reaction has the `misc` metadata field assigned, else `false`. @@ -599,15 +599,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [misc="A reaction"] -misc(reaction) +hasmisc(reaction) ``` """ -function has_misc(reaction::Reaction) +function hasmisc(reaction::Reaction) return hasmetadata(reaction, :misc) end """ -get_misc(reaction::Reaction) +getmisc(reaction::Reaction) Returns `misc` metadata field for the input reaction. @@ -617,7 +617,7 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [misc="A reaction"] -get_misc(reaction) +getmisc(reaction) ``` Notes: @@ -628,8 +628,8 @@ creating a `ReactionSystem` programmatically (in which case any symbolic variabl `misc` metadata field should also be explicitly provided to the `ReactionSystem` constructor). """ -function get_misc(reaction::Reaction) - if has_misc(reaction) +function getmisc(reaction::Reaction) + if hasmisc(reaction) return getmetadata(reaction, :misc) else error("Attempts to access `misc` metadata field for a reaction which does not have a value assigned for this metadata.") diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 5c04b2e642..3be6b6c06c 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -497,7 +497,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; spa end # Extract all quantities encountered in relevant `Reaction` metadata. - has_noise_scaling(rx) && findvars!(ps, us, get_noise_scaling(rx), ivs, vars) + hasnoisescaling(rx) && findvars!(ps, us, getnoisescaling(rx), ivs, vars) end # Extracts any species, variables, and parameters that occur in (non-reaction) equations. diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index 38a634ff77..2a9688cd73 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -121,7 +121,7 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true, for (j, rx) in enumerate(get_rxs(rs)) rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws))) - has_noise_scaling(rx) && (rlsqrt *= get_noise_scaling(rx)) + hasnoisescaling(rx) && (rlsqrt *= getnoisescaling(rx)) remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap)) for (spec, stoich) in rx.netstoich diff --git a/test/dsl/dsl_basic_model_construction.jl b/test/dsl/dsl_basic_model_construction.jl index 72dd01f1ab..6c3d5466c3 100644 --- a/test/dsl/dsl_basic_model_construction.jl +++ b/test/dsl/dsl_basic_model_construction.jl @@ -437,3 +437,5 @@ let @test_throws LoadError @eval @reaction k, 0 --> im @test_throws LoadError @eval @reaction k, 0 --> nothing end + + diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index b6c8ada2e3..93ce90b541 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -584,7 +584,7 @@ let @observables (X, [description="my_description"]) ~ X1 + X2 k, 0 --> X1 + X2 end - @test getdescription(observed(rn)[1].lhs) == "my_description" + @test ModelingToolkit.getdescription(observed(rn)[1].lhs) == "my_description" end # Declares observables implicitly/explicitly. @@ -629,7 +629,7 @@ let (k1, k2), X1 <--> X2 end @test isequal(observed(rn1)[1].lhs, X) - @test getdescription(rn1.X) == "An observable" + @test ModelingToolkit.getdescription(rn1.X) == "An observable" @test isspecies(rn1.X) @test length(unknowns(rn1)) == 2 diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl index ab15dd68ae..ab8488eaf7 100644 --- a/test/reactionsystem_core/events.jl +++ b/test/reactionsystem_core/events.jl @@ -97,9 +97,9 @@ let @test Symbolics.unwrap(rs_ce_de.α) isa Symbolics.BasicSymbolic{Int64} @test Symbolics.unwrap(rs_de.α) isa Symbolics.BasicSymbolic{Int64} @test Symbolics.unwrap(rs_ce_de.α) isa Symbolics.BasicSymbolic{Int64} - @test getdescription(rs_ce_de.A) == "A species" - @test getdescription(rs_de.A) == "A species" - @test getdescription(rs_ce_de.A) == "A species" + @test ModelingToolkit.getdescription(rs_ce_de.A) == "A species" + @test ModelingToolkit.getdescription(rs_de.A) == "A species" + @test ModelingToolkit.getdescription(rs_ce_de.A) == "A species" # Tests that species/variables/parameters can be accessed correctly one a MTK problem have been created. u0 = [X => 1] diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl index eff2ab20b3..4760641c4b 100644 --- a/test/reactionsystem_core/reaction.jl +++ b/test/reactionsystem_core/reaction.jl @@ -116,10 +116,10 @@ let r1 = Reaction(k, [X], [X2], [2], [1]) r2 = Reaction(k, [X], [X2], [2], [1]; metadata=[:noise_scaling => η]) - @test !Catalyst.has_noise_scaling(r1) - @test Catalyst.has_noise_scaling(r2) - @test_throws Exception Catalyst.get_noise_scaling(r1) - @test isequal(Catalyst.get_noise_scaling(r2), η) + @test !Catalyst.hasnoisescaling(r1) + @test Catalyst.hasnoisescaling(r2) + @test_throws Exception Catalyst.getnoisescaling(r1) + @test isequal(Catalyst.getnoisescaling(r2), η) end # Tests the description metadata. @@ -131,10 +131,10 @@ let r1 = Reaction(k, [X], [X2], [2], [1]) r2 = Reaction(k, [X], [X2], [2], [1]; metadata=[:description => "A reaction"]) - @test !Catalyst.has_description(r1) - @test Catalyst.has_description(r2) - @test_throws Exception Catalyst.get_description(r1) - @test isequal(Catalyst.get_description(r2), "A reaction") + @test !Catalyst.hasdescription(r1) + @test Catalyst.hasdescription(r2) + @test_throws Exception Catalyst.getdescription(r1) + @test isequal(Catalyst.getdescription(r2), "A reaction") end # Tests the misc metadata. @@ -146,8 +146,8 @@ let r1 = Reaction(k, [X], [X2], [2], [1]) r2 = Reaction(k, [X], [X2], [2], [1]; metadata=[:misc => ('M', :M)]) - @test !Catalyst.has_misc(r1) - @test Catalyst.has_misc(r2) - @test_throws Exception Catalyst.get_misc(r1) - @test isequal(Catalyst.get_misc(r2), ('M', :M)) + @test !Catalyst.hasmisc(r1) + @test Catalyst.hasmisc(r2) + @test_throws Exception Catalyst.getmisc(r1) + @test isequal(Catalyst.getmisc(r2), ('M', :M)) end \ No newline at end of file diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl index 5684db4132..3e71947495 100644 --- a/test/simulation_and_solving/simulate_SDEs.jl +++ b/test/simulation_and_solving/simulate_SDEs.jl @@ -2,6 +2,7 @@ # Fetch packages. using Catalyst, Statistics, StochasticDiffEq, Test +using Catalyst: getnoisescaling # Sets stable rng number. using StableRNGs @@ -243,8 +244,8 @@ let u0 = [:X1 => 500.0, :X2 => 500.0] p = [:p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, :η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] - @test getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1" - @test getdescription(parameters(noise_scaling_network)[5]) == "Parameter η2" + @test ModelingToolkit.getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1" + @test ModelingToolkit.getdescription(parameters(noise_scaling_network)[5]) == "Parameter η2" sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p) @test sprob.ps[:η1] == sprob.ps[:η2] == sprob.ps[:η3] == sprob.ps[:η4] == 0.0 @@ -325,8 +326,8 @@ let @test issetequal([X, H], species(rs)) @test issetequal([X, H, h], unknowns(rs)) @test issetequal([p, d, η], parameters(rs)) - @test isequal(get_noise_scaling(reactions(rs)[1]), η*H + 1) - @test isequal(get_noise_scaling(reactions(rs)[2]), h) + @test isequal(getnoisescaling(reactions(rs)[1]), η*H + 1) + @test isequal(getnoisescaling(reactions(rs)[2]), h) end # Tests the `remake_noise_scaling` function. @@ -369,9 +370,9 @@ let # Checks that systems have the correct noise scaling terms. rn = set_default_noise_scaling(rn, 0.5) - rn1_noise_scaling = [get_noise_scaling(rx) for rx in Catalyst.get_rxs(rn)] - rn2_noise_scaling = [get_noise_scaling(rx) for rx in Catalyst.get_rxs(Catalyst.get_systems(rn)[1])] - rn_noise_scaling = [get_noise_scaling(rx) for rx in reactions(rn)] + rn1_noise_scaling = [getnoisescaling(rx) for rx in Catalyst.get_rxs(rn)] + rn2_noise_scaling = [getnoisescaling(rx) for rx in Catalyst.get_rxs(Catalyst.get_systems(rn)[1])] + rn_noise_scaling = [getnoisescaling(rx) for rx in reactions(rn)] @test issetequal(rn1_noise_scaling, [2.0, 0.5]) @test issetequal(rn2_noise_scaling, [5.0, 0.5]) @test issetequal(rn_noise_scaling, [2.0, 0.5, 5.0, 0.5])