Merge 8937143 into 4024e67

SciML · May 13, 2024 · 2f0f67d · 2f0f67d
2 parents 4024e67 + 8937143
commit 2f0f67d
Show file tree

Hide file tree

Showing 84 changed files with 8,061 additions and 2,444 deletions.
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -3,9 +3,11 @@ on:
   pull_request:
     branches:
       - master
+      - Catalyst_version_14
   push:
     branches:
       - master
+      - Catalyst_version_14
 jobs:
   test:
     runs-on: ubuntu-latest
@@ -14,7 +16,7 @@ jobs:
         group:
           - Core
         version:
-          - '1'
+          - '1.10.2'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

diff --git a/HISTORY.md b/HISTORY.md
@@ -1,6 +1,32 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
+
+## Catalyst 14.0
+- To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems.
+- Added CatalystStructuralIdentifiabilityExtension, which permits StructuralIdentifiability.jl function to be applied directly to Catalyst systems. E.g. use
+```julia
+using Catalyst, StructuralIdentifiability
+goodwind_oscillator = @reaction_network begin
+    (mmr(P,pₘ,1), dₘ), 0 <--> M
+    (pₑ*M,dₑ), 0 <--> E
+    (pₚ*E,dₚ), 0 <--> P
+end
+assess_identifiability(goodwind_oscillator; measured_quantities=[:M])
+```
+to assess (global) structural identifiability for all parameters and variables of the `goodwind_oscillator` model (under the presumption that we can measure `M` only).
+- Automatically handles conservation laws for structural identifiability problems (eliminates these internally to speed up computations).
+- Adds a tutorial to illustrate the use of the extension.
+- Enable adding metadata to individual reactions, e.g:
+```julia
+rn = @reaction_network begin
+    @parameters η
+    k, 2X --> X2, [noise_scaling=η]
+end
+get_noise_scaling(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.
 - Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
 - LatticeReactionSystem structure represents a spatial reaction network:
   ```julia
@@ -24,6 +50,13 @@
   ```
   X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
 
+- Update how compounds are created. E.g. use
+```julia
+@variables t C(t) O(t)
+@compound CO2 ~ C + 2O
+```
+to create a compound species `CO2` that consists of `C` and 2 `O`.
+- Added documentation for chemistry related functionality (compound creation and reaction balancing).
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst
@@ -35,7 +68,6 @@ wilhelm_2009_model = @reaction_network begin
     k5, 0 --> X
 end
 
-
 using BifurcationKit
 bif_par = :k1
 u_guess = [:X => 5.0, :Y => 2.0]

diff --git a/Project.toml b/Project.toml
@@ -6,6 +6,8 @@ version = "13.5.1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -26,41 +28,47 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [weakdeps]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 
 [extensions]
 CatalystBifurcationKitExtension = "BifurcationKit"
 CatalystHomotopyContinuationExtension = "HomotopyContinuation"
+CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
 BifurcationKit = "0.3"
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
+DynamicQuantities = "0.13.2"
 Graphs = "1.4"
-JumpProcesses = "9.3.2"
 HomotopyContinuation = "2.9"
+JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "8.73"
+ModelingToolkit = "9.11.0"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 Setfield = "1"
+StructuralIdentifiability = "0.5.1"
 SymbolicUtils = "1.0.3"
-Symbolics = "5.0.3"
+Symbolics = "5.27"
 Unitful = "1.12.4"
 julia = "1.9"
 
 [extras]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -69,8 +77,9 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
+test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
diff --git a/README.md b/README.md
@@ -31,10 +31,7 @@ etc).
 
 ## Breaking changes and new features
 
-**NOTE:** version 13 is a breaking release, with changes to simplify the DSL
-notation while also adding more features, changes to how chemical species are
-specified symbolically when directly building `ReactionSystem`s, and changes that
-simplify how to include ODE or algebraic constraint equation.
+**NOTE:** version 14 is a breaking release, prompted by the release of ModelingToolkit.jl version 9. This caused several breaking changes in how Catalyst models are represented and interfaced with.
 
 Breaking changes and new functionality are summarized in the
 [HISTORY.md](HISTORY.md) file.
@@ -48,6 +45,8 @@ the current master branch.
 
 Several Youtube video tutorials and overviews are also available, but note these use older 
 Catalyst versions with slightly different notation (for example, in building reaction networks):
+- From JuliaCon 2023: A short 15 minute overview of Catalyst as of version 13 is
+available in the talk [Catalyst.jl, Modeling Chemical Reaction Networks](https://www.youtube.com/watch?v=yreW94n98eM&ab_channel=TheJuliaProgrammingLanguage).
 - From JuliaCon 2022: A three hour tutorial workshop overviewing how to use
   Catalyst and its more advanced features as of version 12.1. [Workshop
   video](https://youtu.be/tVfxT09AtWQ), [Workshop Pluto.jl
@@ -60,6 +59,8 @@ Catalyst.jl](https://www.youtube.com/watch?v=5p1PJE5A5Jw).
   Modelling of Biochemical Reaction
   Networks](https://www.youtube.com/watch?v=s1e72k5XD6s)
 
+Finally, an overview of the package and its features (as of version 13) can also be found in its corresponding research paper, [Catalyst: Fast and flexible modeling of reaction networks](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530).
+
 ## Features
 
 - A DSL provides a simple and readable format for manually specifying chemical

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -2,18 +2,21 @@
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 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"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+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"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -22,29 +25,33 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 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]
 BifurcationKit = "0.3"
 Catalyst = "13"
 DataFrames = "1"
+DiffEqParamEstim = "2.1"
 DifferentialEquations = "7.7"
 Distributions = "0.25"
 Documenter = "0.27"
 HomotopyContinuation = "2.6"
 Latexify = "0.15, 0.16"
-ModelingToolkit = "8.47"
-NonlinearSolve = "1.6.0, 2"
+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"
-PEtab = "2"
 Plots = "1.36"
-SciMLBase = "~2.5"
+SciMLBase = "2.13"
 SciMLSensitivity = "7.19"
 Setfield = "1.1"
 SpecialFunctions = "2.1"
-SteadyStateDiffEq = "1"
+SteadyStateDiffEq = "2.0.1"
 StochasticDiffEq = "6"
-Symbolics = "5.0.3"
+StructuralIdentifiability = "0.5.1"
+Symbolics = "5.14"
diff --git a/docs/pages.jl b/docs/pages.jl
@@ -1,12 +1,13 @@
 pages = Any["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"],
             "Catalyst Functionality" => Any["catalyst_functionality/dsl_description.md",
                                             "catalyst_functionality/programmatic_CRN_construction.md",
                                             "catalyst_functionality/compositional_modeling.md",
                                             "catalyst_functionality/constraint_equations.md",
                                             "catalyst_functionality/parametric_stoichiometry.md",
                                             "catalyst_functionality/network_analysis.md",
+                                            "catalyst_functionality/chemistry_related_functionality.md",
                                             "Model creation examples" => Any["catalyst_functionality/example_networks/basic_CRN_examples.md",
                                                                              "catalyst_functionality/example_networks/hodgkin_huxley_equation.md",
                                                                              "catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md"]],
@@ -15,7 +16,9 @@ pages = Any["Home" => "index.md",
                                            "catalyst_applications/homotopy_continuation.md",
                                            "catalyst_applications/nonlinear_solve.md",
                                            "catalyst_applications/bifurcation_diagrams.md"],
-            "Inverse Problems" => Any["inverse_problems/parameter_estimation.md",
-                                      "inverse_problems/petab_ode_param_fitting.md"],
+            "Inverse Problems" => Any["inverse_problems/optimization_ode_param_fitting.md",
+                                        #"inverse_problems/petab_ode_param_fitting.md",
+                                        "inverse_problems/structural_identifiability.md",
+                                        "Inverse problem examples" => Any["inverse_problems/examples/ode_fitting_oscillation.md"]],
             "FAQs" => "faqs.md",
             "API" => "api/catalyst_api.md"]
diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
@@ -36,8 +36,8 @@ models, and stochastic chemical kinetics jump process models.
 
 ```@example ex1
 using Catalyst, DifferentialEquations, Plots
+t = default_t()
 @parameters β γ
-@variables t
 @species S(t) I(t) R(t)
 
 rxs = [Reaction(β, [S,I], [I], [1,1], [2])
@@ -95,11 +95,11 @@ retrieve info from just a base [`ReactionSystem`](@ref) `rn`, ignoring
 sub-systems of `rn`, one can use the ModelingToolkit accessors (these provide
 direct access to the corresponding internal fields of the `ReactionSystem`)
 
-* `ModelingToolkit.get_states(rn)` is a vector that collects all the species
+* `ModelingToolkit.get_unknowns(rn)` is a vector that collects all the species
   defined within `rn`, ordered by species and then non-species variables.
 * `Catalyst.get_species(rn)` is a vector of all the species variables in the system. The
   entries in `get_species(rn)` correspond to the first `length(get_species(rn))`
-  components in `get_states(rn)`.
+  components in `get_unknowns(rn)`.
 * `ModelingToolkit.get_ps(rn)` is a vector that collects all the parameters
   defined *within* reactions in `rn`.
 * `ModelingToolkit.get_eqs(rn)` is a vector that collects all the
@@ -120,10 +120,10 @@ To retrieve information from the full reaction network represented by a system
 `rn`, which corresponds to information within both `rn` and all sub-systems, one
 can call:
 
-* `ModelingToolkit.states(rn)` returns all species *and variables* across the
+* `ModelingToolkit.unknowns(rn)` returns all species *and variables* across the
   system, *all sub-systems*, and all constraint systems. Species are ordered
-  before non-species variables in `states(rn)`, with the first `numspecies(rn)`
-  entires in `states(rn)` being the same as `species(rn)`.
+  before non-species variables in `unknowns(rn)`, with the first `numspecies(rn)`
+  entires in `unknowns(rn)` being the same as `species(rn)`.
 * [`species(rn)`](@ref) is a vector collecting all the chemical species within
   the system and any sub-systems that are also `ReactionSystems`.
 * `ModelingToolkit.parameters(rn)` returns all parameters across the
@@ -181,22 +181,14 @@ reactionrates
 ```
 
 ## [Functions to extend or modify a network](@id api_network_extension_and_modification)
-`ReactionSystem`s can be programmatically extended using
-[`@add_reactions`](@ref), [`addspecies!`](@ref), [`addparam!`](@ref) and
-[`addreaction!`](@ref), or using [`ModelingToolkit.extend`](@ref) and
+`ReactionSystem`s can be programmatically extended using [`ModelingToolkit.extend`](@ref) and
 [`ModelingToolkit.compose`](@ref).
 
 ```@docs
-@add_reactions
-addspecies!
-addparam!
-addreaction!
 setdefaults!
 ModelingToolkit.extend
 ModelingToolkit.compose
 Catalyst.flatten
-merge!(network1::ReactionSystem, network2::ReactionSystem)
-reorder_states!
 ```
 
 ## Network analysis and representations
@@ -251,6 +243,8 @@ latexify(rn; form=:sde)
 (As of writing this, an upstream bug causes the SDE form to be erroneously
 displayed as the ODE form)
 
+Finally, another optional argument (`expand_functions=true`) automatically expands functions defined by Catalyst (such as `mm`). To disable this, set `expand_functions=false`.
+
 If [Graphviz](https://graphviz.org/) is installed and commandline accessible, it
 can be used to create and save network diagrams using [`Graph`](@ref) and
 [`savegraph`](@ref).
@@ -280,6 +274,17 @@ Base.convert
 ModelingToolkit.structural_simplify
 ```
 
+## Chemistry-related functionalities
+Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology).
+```@docs
+@compound
+@compounds
+iscompound
+components
+coefficients
+component_coefficients
+```
+
 ## Unit validation
 ```@docs
 validate(rx::Reaction; info::String = "")