diff --git a/HISTORY.md b/HISTORY.md index 65a8fd4374..e0479f8cb8 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,6 +1,21 @@ # Breaking updates and feature summaries across releases ## Catalyst unreleased (master branch) + +## Catalyst 14.0 +- 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. - 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 @@ -35,7 +50,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 index 5f1a189aff..6c0f521f2d 100644 --- a/Project.toml +++ b/Project.toml @@ -26,10 +26,12 @@ 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" @@ -48,8 +50,9 @@ 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.14" Unitful = "1.12.4" julia = "1.9" @@ -69,8 +72,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", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"] diff --git a/README.md b/README.md index 587a71b7e6..4359a4f6ec 100644 --- 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 index 2ea3b3fda9..1621bc3f26 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,6 +7,7 @@ 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" @@ -22,6 +23,7 @@ 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] @@ -34,17 +36,18 @@ Documenter = "0.27" HomotopyContinuation = "2.6" Latexify = "0.15, 0.16" ModelingToolkit = "8.47" -NonlinearSolve = "1.6.0, 2" +NonlinearSolve = "3.4.0" Optim = "1" Optimization = "3.19" 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 index 9b21e12b10..8bfd65c9cb 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -16,6 +16,7 @@ pages = Any["Home" => "index.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/petab_ode_param_fitting.md", + "inverse_problems/structural_identifiability.md"], "FAQs" => "faqs.md", "API" => "api/catalyst_api.md"] \ No newline at end of file diff --git a/docs/src/catalyst_applications/bifurcation_diagrams.md b/docs/src/catalyst_applications/bifurcation_diagrams.md index bb73ae368a..55eab485a2 100644 --- a/docs/src/catalyst_applications/bifurcation_diagrams.md +++ b/docs/src/catalyst_applications/bifurcation_diagrams.md @@ -1,6 +1,6 @@ # [Bifurcation Diagrams](@id bifurcation_diagrams) -Bifurcation diagrams describe how, for a dynamical system, the quantity and type of its steady states change as a parameter is varied[^1]. When using Catalyst-generated models, these can be computed with the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. Catalyst provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `ReactionSystem`s. If you use this feature in your research, please [cite the BifurcationKit.jl](@ref bifurcation_kit_citation) and [Catalyst.jl](@ref catalyst_citation) publications. +Bifurcation diagrams describe how, for a dynamical system, the quantity and type of its steady states change as a parameter is varied[^1]. When using Catalyst-generated models, these can be computed with the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. Catalyst provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `ReactionSystem`s. This tutorial briefly introduces how to use Catalyst with BifurcationKit through basic examples, with BifurcationKit.jl providing [a more extensive documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/). Especially for more complicated systems, where careful tuning of algorithm options might be required, reading the BifurcationKit documentation is recommended. Finally, BifurcationKit provides many additional features not described here, including [computation of periodic orbits](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/), [tracking of bifurcation points along secondary parameters](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/branchswitching/), and [bifurcation computations for PDEs](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials/#PDEs:-bifurcations-of-equilibria). diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 9b961ffb23..1fb751bc15 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -9,7 +9,7 @@ integer Hill exponents). The roots of these can reliably be found through a *homotopy continuation* algorithm. This is implemented in Julia through the [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/) package. -Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, `hc_steady_states`, that can be used to find the steady states of any `ReactionSystem` structure. If you use this in your research, please [cite the HomotopyContinuation.jl](@ref homotopy_continuation_citation) and [Catalyst.jl](@ref catalyst_citation) publications. +Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, `hc_steady_states`, that can be used to find the steady states of any `ReactionSystem` structure. ## Basic example For this tutorial, we will use a model from Wilhelm (2009)[^1] (which diff --git a/docs/src/catalyst_applications/nonlinear_solve.md b/docs/src/catalyst_applications/nonlinear_solve.md index 4692105d56..0d4a903714 100644 --- a/docs/src/catalyst_applications/nonlinear_solve.md +++ b/docs/src/catalyst_applications/nonlinear_solve.md @@ -6,7 +6,7 @@ We have previously described how `ReactionSystem` steady states can be found thr However, if all (or multiple) steady states are sought, using homotopy continuation is better. -This tutorial describes how to create `NonlinearProblem`s from Catalyst's `ReactionSystemn`s, and how to solve them using NonlinearSolve. More extensive descriptions of available solvers and options can be found in [NonlinearSolve's documentation](https://docs.sciml.ai/NonlinearSolve/stable/). If you use this in your research, please [cite the NonlinearSolve.jl](@ref nonlinear_solve_citation) and [Catalyst.jl](@ref catalyst_citation) publications. +This tutorial describes how to create `NonlinearProblem`s from Catalyst's `ReactionSystemn`s, and how to solve them using NonlinearSolve. More extensive descriptions of available solvers and options can be found in [NonlinearSolve's documentation](https://docs.sciml.ai/NonlinearSolve/stable/). ## Basic example Let us consider a simple dimerisation network, where a protein ($P$) can exist in a monomer and a dimer form. The protein is produced at a constant rate from its mRNA, which is also produced at a constant rate. 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 48c3593ec2..e1e8af412c 100644 --- a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md +++ b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md @@ -1,181 +1,176 @@ # [Introduction to Catalyst and Julia for New Julia users](@id catalyst_for_new_julia_users) The Catalyst tool for the modelling of chemical reaction networks is based in the Julia programming language. While experience in Julia programming is advantageous for using Catalyst, it is not necessary for accessing some of its basic features. This tutorial serves as an introduction to Catalyst for those unfamiliar with Julia, also introducing some basic Julia concepts. Anyone who plans on using Catalyst extensively is recommended to familiarise oneself more thoroughly with the Julia programming language. A collection of resources for learning Julia can be found [here](https://julialang.org/learning/), and a full documentation is available [here](https://docs.julialang.org/en/v1/). -Julia can be downloaded [here](https://julialang.org/downloads/). +Julia can be downloaded [here](https://julialang.org/downloads/). Generally, it is recommended to use the [*juliaup*](https://github.com/JuliaLang/juliaup) tool to install and update Julia. Furthermore, *Visual Studio Code* is a good IDE with [extensive Julia support](https://code.visualstudio.com/docs/languages/julia), and a good default choice. *Users who are already familiar with Julia can skip to the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial.* ## Basic Julia usage On the surface, Julia has many similarities to languages like MATLAB, Python, and R. -*Values* can be assigned to *variables* through the use of a `=` sign. Values (possibly stored in variables) can be used for most basic computations. +*Values* can be assigned to *variables* through `=` sign. Values (possibly stored in variables) can be used for most basic computations. ```@example ex1 length = 2.0 width = 4.0 -area = length*width +area = length * width ``` -*Functions* take one or more inputs (enclosed by `()`) and return some output. E.g. the `min` function returns the minimum of two values +*Functions* take one or more inputs (enclosed by `()`) and return some output. E.g. the `min` function returns the minimum of two values. ```@example ex1 min(1.0, 3.0) ``` -A line of Julia code is not required to end with `;`, however, if it does, the output of that line is not displayed. -```julia -min(1.0, 3.0); -``` - Each Julia variable has a specific *type*, designating what type of value it is. While not directly required to use Catalyst, this is useful to be aware of. To learn the type of a specific variable, use the `typeof` function. More information about types can be [found here](https://docs.julialang.org/en/v1/manual/types/). ```@example ex1 typeof(1.0) ``` - Here, `Float64` denotes decimal-valued numbers. Integer-valued numbers instead are of the `Int64` type. ```@example ex1 typeof(1) ``` +There exists a large number of Julia types (with even more being defined by various packages). Additional examples include `String`s (defined by enclosing text within `" "`): +```@example ex1 +"Hello world!" +``` +and `Symbol`s (defined by pre-appending an expression with `:`): +```@example ex1 +:Julia +``` -Finally, we note that the first time some code is run in Julia, it has to be *compiled*. However, this is only required once per Julia session. Hence, the second time the same code is run, it runs much faster. E.g. try running this line of code first one time, and then one additional time. You will note that the second run is much faster. +Finally, we note that the first time some code is run in Julia, it has to be [*compiled*](https://en.wikipedia.org/wiki/Just-in-time_compilation). However, this is only required once per Julia session. Hence, the second time the same code is run, it runs much faster. E.g. try running this line of code first one time, and then one additional time. You will note that the second run is much faster. ```@example ex1 rand(100, 100)^3.5 nothing # hide ``` (This code creates a random 100x100 matrix, and takes it to the power of 3.5) -This is useful to know when you e.g. declare, simulate, or plot, a Catalyst model. The first time you run a command there might be a slight delay. However, subsequent runs will execute much quicker. This holds even if you do minor adjustments before the second run (such as changing simulation initial conditions). +This is useful to know when you e.g. declare, simulate, or plot, a Catalyst model. The first time you run a command there might be a slight delay. However, subsequent runs will be much quicker. This holds even if you make minor adjustments before the second run (such as changing simulation initial conditions). -## Installing and activating packages -Except for some base Julia packages (such as Pkg, the package manager) that are available by default, Julia packages must be installed locally before they can be used. Most packages are registered with Julia, and can be added through the `Pkg.add("DesiredPackage")` command (where `DesiredPackage` is the name of the package you wish to install). We can thus install Catalyst: +## [Installing and activating packages](@id catalyst_for_new_julia_users_packages_intro) +Due to its native package manager (Pkg), and a registry of almost all packages of relevancy, package management in Julia is unusually easy. Here, we will briefly describe how to install and activate Catalyst (and two additional packages relevant to this tutorial). + +To import a Julia package into a session, you can use the `using PackageName` command (where `PackageName` is the name of the package you wish to import). However, before you can do so, it must first be installed on your computer. This is done through the `Pkg.add("PackageName")` command: ```julia using Pkg Pkg.add("Catalyst") ``` - -Here, the command `using Pkg` is required to activate the Pkg package manager. - -Next, we also wish to add 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 `DifferentialEquations` and `Plots` packages (for numeric simulation of models, and plotting, respectively). ```julia Pkg.add("DifferentialEquations") Pkg.add("Plots") ``` -Once a package has been installed through the `Pkg.add` command, this command does not have to be repeated in further Julia sessions on the same machine. - -Installing a Julia package is, however, not enough to use it. Before a package's features are used in a Julia session, it has to be loaded through the `using DesiredPackage` command (where `DesiredPackage` is the name of the package you wish to activate). This command has to be repeated whenever a Julia session is restarted. - -We thus activate our three desired packages: - -```@example ex1 +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 Plots ``` +Here, if we restart Julia, these commands *do need to be rerun. -For a more detailed introduction to Julia packages, please read [the Pkg documentation](https://docs.julialang.org/en/v1/stdlib/Pkg/). +A more comprehensive (but still short) introduction to package management in Julia is provided at [the end of this documentation page](catalyst_for_new_julia_users_packages). It contains some useful information and is hence highly recommended reading. For a more detailed introduction to Julia package management, please read [the Pkg documentation](https://docs.julialang.org/en/v1/stdlib/Pkg/). ## Simulating a basic Catalyst model -Now that we have some basic familiarity with Julia, and have installed and activated the required packages, we will create and simulate a basic chemical reaction network model through Catalyst. +Now that we have some basic familiarity with Julia, and have installed and imported the required packages, we will create and simulate a basic chemical reaction network model using Catalyst. -Catalyst models are created through the `@reaction_network` *macro*. For more information on macros, please read [the Julia documentation on macros](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros). This documentation is, however, rather advanced (and not required to use Catalyst). We instead recommend that you simply familiarise yourself with the Catalyst syntax, without studying in detail how macros work and what they are. +Catalyst models are created through the `@reaction_network` *macro*. For more information on macros, please read [the Julia documentation on macros](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros). This documentation is, however, rather advanced (and not required to use Catalyst). We instead recommend that you simply familiarise yourself with the Catalyst syntax, without studying in detail how macros work and what they are. -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 itself 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. +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`. - -```@example ex1 +Here, we create a simple *birth-death* model, where a single species ($X$) is created at rate $b$, and degraded at rate $d$. The model is stored in the variable `rn`. +```@example ex2 rn = @reaction_network begin b, 0 --> X d, X --> 0 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/). Next, we wish to simulate our model. To do this, we need to provide some additional information to the simulator. This is -* The initial condition. That is, the concentration or numbers of each species at the start of the simulation. +* The initial condition. That is, the concentration (or copy numbers) of each species at the start of the simulation. * The timespan. That is, the timeframe over which we wish to run the simulation. * The parameter values. That is, the values of the model's parameters for this simulation. -The initial condition is given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specific within brackets, `[]`, and separated by `,`. Since we only have one species, the vector holds a single element. In this element, we set the value of *X* using the `:X => 1.0` syntax. Here, we first denote the name of the species (with a `:` pre-appended), next follows a `=>` and then the value of *X*. Since we wish to simulate the *concentration* of X over time, we will let the initial condition be decimal valued. -```@example ex1 +The initial condition is given as a *Vector*. This is a type which collects several different values. To declare a vector, the values are specific within brackets, `[]`, and separated by `,`. Since we only have one species, the vector holds a single element. In this element, we set the value of $X$ using the `:X => 1.0` syntax. Here, we first denote the name of the species (with a `:` pre-appended, i.e. creating a `Symbol`), next follows a `=>` and then the value of $X$. Since we wish to simulate the *concentration* of X over time, we will let the initial condition be decimal valued. +```@example ex2 u0 = [:X => 1.0] ``` -The timespan sets the time point at which we start the simulation (typically `0.0` is used) and the final time point of the simulation. These are combined into a two-valued *Tuple*. Tuples are similar to vectors, but are enclosed by `()` and not `[]`. Again, we will let both time points be decimal valued. -```@example ex1 +The timespan sets the time point at which we start the simulation (typically `0.0` is used) and the final time point of the simulation. These are combined into a two-valued *tuple*. Tuples are similar to vectors, but are enclosed by `()` and not `[]`. Again, we will let both time points be decimal valued. +```@example ex2 tspan = (0.0, 10.0) ``` -Finally, the parameter values are, like the initial conditions, given as a vector. Since we have two parameters (*b* and *d*), the parameter vector has two values. We use a similar notation for setting the parameter values as the initial condition (first the parameter, then an arrow, then the value). -```@example ex1 +Finally, the parameter values are, like the initial conditions, given in a vector. Since we have two parameters ($b$ and $d$), the parameter vector has two values. We use a similar notation for setting the parameter values as the initial condition (first the parameter, then an arrow, then the value). +```@example ex2 params = [:b => 1.0, :d => 0.2] ``` -Please read here for more information on [Vectors](https://docs.julialang.org/en/v1/manual/arrays/) and [Tuples](https://docs.julialang.org/en/v1/manual/types/#Tuple-Types). +Please read here for more information on [vectors](https://docs.julialang.org/en/v1/manual/arrays/) and [tuples](https://docs.julialang.org/en/v1/manual/types/#Tuple-Types). Next, before we can simulate our model, we bundle all the required information together in a so-called `ODEProblem`. Note that the order in which the input (the model, the initial condition, the timespan, and the parameter values) is provided to the ODEProblem matters. E.g. the parameter values cannot be provided as the first argument, but have to be the fourth argument. Here, we save our `ODEProblem` in the `oprob` variable. - - -```@example ex1 +```@example ex2 oprob = ODEProblem(rn, u0, tspan, params) ``` We can now simulate our model. We do this by providing the `ODEProblem` to the `solve` function. We save the output to the `sol` variable. -```@example ex1 +```@example ex2 sol = solve(oprob) ``` Finally, we can plot the solution through the `plot` function. -```@example ex1 +```@example ex2 plot(sol) ``` -Here, the plot shows the time evolution of the concentration of the species *X* from its initial condition. +Here, the plot shows the time evolution of the concentration of the species $X$ from its initial condition. For more information about the numerical simulation package, please see the [DifferentialEquations documentation](https://docs.sciml.ai/DiffEqDocs/stable/). For more information about the plotting package, please see the [Plots documentation](https://docs.juliaplots.org/stable/). ## Additional modelling example -To make this introduction more comprehensive, we here provide another example, using a more complicated model. In addition, instead of simulating our model as concentrations evolve over time, we will simulate the individual reaction events through the [Gillespie algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm). This is a way to add *noise* to our model. +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. -This time, we will declare the 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 dynamics. 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](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. + It also has 2 reaction events: * Infection, where a susceptible individual meets an infected individual and also becomes infected. -* Recovery, where an infected individual recovers. +* Recovery, where an infected individual recovers from the infection. + Each reaction is also associated with a specific rate (corresponding to a parameter). * *b*, the infection rate. * *k*, the recovery rate. We declare the model using the `@reaction_network` macro, and store it in the `sir_model` variable. -```@example ex1 +```@example ex2 sir_model = @reaction_network begin b, S + I --> 2I k, I --> R end ``` - Note that the first reaction contains two different substrates (separated by a `+` sign). While there is only a single product (*I*), two copies of *I* are produced. The *2* in front of the product *I* denotes this. -Next, we declare our initial condition, time span, and parameter values. Since we want to simulate the individual reaction events, that discretely change the state of our model, we want our initial conditions to be integer-valued. We will start with a mostly susceptible population, but where a single individual has been infected through some means. -```@example ex1 -u0 = [:S => 50, :I => 1, :R => 0.0] +Next, we declare our initial condition, time span, and parameter values. Since we want to simulate the individual reaction events that discretely change the state of our model, we want our initial conditions to be integer-valued. We will start with a mostly susceptible population, but where a single individual has been infected through some means. +```@example ex2 +u0 = [:S => 50, :I => 1, :R => 0] tspan = (0.0, 10.0) params = [:b => 0.2, :k => 1.0] +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 denotes 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 ex1 +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 dprob = DiscreteProblem(sir_model, u0, tspan, params) jprob = JumpProblem(sir_model, dprob, Direct()) ``` -Again, the order in which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. The last argument to the `JumpProblem` (`Direct()`) denotes which simulation method we wish to use. For now, we recommend the user simply use the `Direct()` option, and then consider alternative ones (see the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/)) when they are more familiar with modelling in Catalyst and Julia. +Again, the order with which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. The last argument to the `JumpProblem` (`Direct()`) denotes which simulation method we wish to use. For now, we recommend that users simply use the `Direct()` option, and then consider alternative ones (see the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/)) when they are more familiar with modelling in Catalyst and Julia. Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. Here, the `solve` function also has a second argument (`SSAStepper()`). This is a time stepping algorithm that calls the `Direct` solver to advance a simulation. Again, we recommend at this stage you simply use this option, and then explore exactly what this means at a later stage. -```@example ex1 +```@example ex2 sol = solve(jprob, SSAStepper()) sol = solve(jprob, SSAStepper(); seed=1234) # hide plot(sol) @@ -183,6 +178,69 @@ plot(sol) **Exercise:** Try simulating the model several times. Note that the epidemic doesn't always take off, but sometimes dies out without spreading through the population. Try changing the infection rate (*b*), determining how this value affects the probability that the epidemic goes through the population. +## [Package management in Julia](@id catalyst_for_new_julia_users_packages) +We have previously introduced how to install and activate Julia packages. While this is enough to get started with Catalyst, for long-term users, there are some additional considerations for a smooth experience. These are described here. + +### [Setting up a new Julia environment](@id catalyst_for_new_julia_users_packages_environments) +Whenever you run Julia, it will run in a specific *environment*. You can specify any folder on your computer as a Julia environment. Some modes of running Julia will automatically use the environment corresponding to the folder you start Julia in. Others (or if you start Julia in a folder without an environment), will use your *default* environment. In these cases you can, during your session, switch to another environment. While it is possible to not worry about environments (and always use the default one), this can lead to long-term problems as more packages are installed. + +To activate your current folder as an environment, run the following commands: +```julia +using Pkg +Pkg.activate(".") +``` +This will: +1. If your current folder (which can be displayed using the `pwd()` command) is not designated as a possible Julia environment, designate it as such. +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). + +### [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): +```@example ex3 +using Pkg +``` +Latexify is a registered package, so it can be installed directly using: + ```julia +Pkg.add("Latexify") +``` +Finally, to import Latexify into our current Julia session we use: +```@example ex3 +using Latexify +``` +Here, `using Latexify` must be rerun whenever you restart a Julia session. However, you only need to run `Pkg.add("Latexify")` once to install it on your computer (but possibly additional times to add it to new environments, see the next section). + +### [Why environments are important](@id catalyst_for_new_julia_users_packages_environment_importance) +We have previously described how to set up new Julia environments, how to install Julia packages, and how to import them into a current session. Let us say that you were to restart Julia in a new folder and activate this as a separate environment. If you then try to import Latexify through `using Latexify` you will receive an error claiming that Latexify was not found. The reason is that the `Pkg.add("Latexify")` command actually carries out two separate tasks: +1. If Latexify is not already installed on your computer, install it. +2. Add Latexify as an available package to your current environment. + +Here, while Catalyst has previously been installed on your computer, it has not been added to the new environment you created. To do so, simply run +```julia +using Pkg +Pkg.add("Latexify") +``` +after which Catalyst can be imported through `using Latexify`. You can get a list of all packages available in your current environment using: +```julia +Pkg.status() +``` + +So, why is this required, and why cannot we simply import any package installed on our computer? The reason is that most packages depend on other packages, and these dependencies may be restricted to only specific versions of these packages. This creates complicated dependency graphs that restrict what versions of what packages are compatible with each other. When you use `Pkg.add("PackageName")`, only a specific version of that package is actually added (the latest possible version as permitted by the dependency graph). Here, Julia environments both define what packages are available *and* their respective versions (these versions are also displayed by the `Pkg.status()` command). By doing this, Julia can guarantee that the packages (and their versions) specified in an environment are compatible with each other. + +The reason why all this is important is that it is *highly recommended* to, for each project, define a separate environment. To these, only add the required packages. General-purpose environments with a large number of packages often produce package-incompatibility issues. While these might not prevent you from installing all desired package, they often mean that you are unable to use the latest version of some packages. + +!!! 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. + +Some additional useful Pkg commands are: +- `Pk.rm("PackageName")` removes a package from the current environment. +- `Pkg.update(PackageName")`: updates the designated package. + +!!! note + A useful feature of Julia's environment system is that enables the exact definition of what packages and versions were used to execute a script. This supports e.g. reproducibility in academic research. Here, by providing the corresponding Project.toml and Manifest.toml files, you can enable someone to reproduce the exact program just to perform some set of analyses. + + --- ## Feedback If you are a new Julia user who has used this tutorial, and there was something you struggled with or would have liked to have explained better, please [raise an issue](https://github.com/SciML/Catalyst.jl/issues). That way, we can continue improving this tutorial. diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/src/inverse_problems/petab_ode_param_fitting.md index 47f9056572..fce182e78e 100644 --- a/docs/src/inverse_problems/petab_ode_param_fitting.md +++ b/docs/src/inverse_problems/petab_ode_param_fitting.md @@ -1,5 +1,5 @@ # [Parameter Fitting for ODEs using PEtab.jl](@id petab_parameter_fitting) -The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data ([please cite the corresponding papers if you use the PEtab approach in your research](@ref petab_citations))[^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data and is a good default choice when fitting reaction rate equation ODE models. This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). +The [PEtab.jl package](https://github.com/sebapersson/PEtab.jl) implements the [PEtab format](https://petab.readthedocs.io/en/latest/) for fitting the parameters of deterministic CRN models to data [^1]. PEtab.jl both implements methods for creating cost functions (determining how well parameter sets fit to data), and for minimizing these cost functions. The PEtab approach covers most cases of fitting deterministic (ODE) models to data and is a good default choice when fitting reaction rate equation ODE models. This page describes how to combine PEtab.jl and Catalyst for parameter fitting, with the PEtab.jl package providing [a more extensive documentation](https://sebapersson.github.io/PEtab.jl/stable/) (this tutorial is partially an adaptation of this documentation). ## Introductory example @@ -127,7 +127,7 @@ fitted_sol = solve(oprob_fitted, Tsit5()) plot!(fitted_sol; idxs=:P, label="Fitted solution", linestyle=:dash, lw=6, color=:lightblue) ``` -Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, that even if PEtab.jl has found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to the concept of *identifiability*, which is very important for parameter fitting. +Here we use the `get_ps` function to retrieve a full parameter set using the optimal parameters. Alternatively, the `ODEProblem` or fitted simulation can be retrieved directly using the `get_odeproblem` or `get_odesol` [functions](https://sebapersson.github.io/PEtab.jl/dev/API_choosen/#PEtab.get_odeproblem), respectively (and the initial condition using the `get_u0` function). The calibration result can also be found in `res.xmin`, however, note that PEtab automatically ([unless a linear scale is selected](@ref petab_parameters_scales)) converts parameters to logarithmic scale, so typically `10 .^res.xmin` are the values of interest. If you investigate the result from this example you might note, that even if PEtab.jl has found the global optimum (which fits the data well), this does not actually correspond to the true parameter set. This phenomenon is related to the [concept of *identifiability*](@ref structural_identifiability), which is very important for parameter fitting. ### Final notes PEtab.jl also supports [multistart optimisation](@ref petab_multistart_optimisation), [automatic pre-equilibration before simulations](https://sebapersson.github.io/PEtab.jl/stable/Brannmark/), and [events](@ref petab_events). Various [plot recipes](@ref petab_plotting) exist for investigating the optimisation process. Please read the [PETab.jl documentation](https://sebapersson.github.io/PEtab.jl/stable/) for a more complete description of the package's features. Below follows additional details of various options and features (generally, PEtab is able to find good default values for most options that are not specified). diff --git a/docs/src/inverse_problems/structural_identifiability.md b/docs/src/inverse_problems/structural_identifiability.md new file mode 100644 index 0000000000..23aeb13e82 --- /dev/null +++ b/docs/src/inverse_problems/structural_identifiability.md @@ -0,0 +1,171 @@ +# [Structural Identifiability Analysis](@id structural_identifiability) +During parameter fitting, parameter values are inferred from data. Parameter identifiability refers to whether inferring parameter values for a given model is mathematically feasible. Ideally, parameter fitting should always be accompanied with an identifiability analysis of the problem. + +Identifiability can be divided into *structural* and *practical* identifiability[^1]. Structural identifiability considers only the mathematical model, and which parameters are and are not inherently identifiable due to model structure. Practical identifiability also considers the available data, and determines what system quantities can be inferred from it. In the idealised case of an infinite amount of non-noisy data, practical identifiability converges to structural identifiability. Generally, structural identifiability is assessed before parameters are fitted, while practical identifiability is assessed afterwards. + +Structural identifiability (which is what this tutorial considers) can be illustrated by the following differential equation: +${dx \over dt} = p1*p2*x(t)$ +where, however much data is collected on $x$, it is impossible to determine the distinct values of $p1$ and $p2$. Hence, these parameters are non-identifiable (however, their product, $p1*p2$, *is* identifiable). + +Catalyst contains a special extension for carrying out structural identifiability analysis of generated reaction rate equation ODE models using the [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) package. This enables StructuralIdentifiability's `assess_identifiability`, `assess_local_identifiability`, and `find_identifiable_functions` functions to be called directly on Catalyst `ReactionSystem`s. It also implements specialised routines to make these more efficient when applied to reaction network models (e.g. by improving runtimes). How to use these functions is described in the following tutorial, with [StructuralIdentifiability providing a more extensive documentation](https://docs.sciml.ai/StructuralIdentifiability/stable/). + +Structural identifiability can be divided into *local* and *global* identifiability. If a model quantity is locally identifiable, it means that its true value can be determined down to a finite-number of possible options. This also means that there is some limited region around the quantity's true value where this true value is the only possible value (and hence, within this region, the quantity is fully identifiable). Globally identifiable quantities' values, on the other hand, can be uniquely determined. Again, while identifiability can be confirmed structurally for a quantity, it does not necessarily mean that it is practically identifiable for some given data. + +Generally, there are three types of quantities for which identifiability can be assessed. +- Parameters (e.g. $p1$ and $p2$). +- Full variable trajectories (e.g. $x(t)$). +- Variable initial conditions (e.g. $x(0)$). + +StructuralIdentifiability currently assesses identifiability for the first two only (however, if $x(t)$ is identifiable, then $x(0)$ will be as well). + +!!! note + Currently, the StructuralIdentifiability.jl extension only considers structural identifiability for the ODE generated by the reaction rate equation. It is possible that for the SDE model (generated by the chemical Langevin equation) and the jump model (generated by stochastic chemical kinetics) the identifiability of model quantities is different. + +## Global identifiability analysis + +### Basic example + +Global identifiability can be assessed using the `assess_identifiability` function. For each model quantity (parameters and variables), it will assess whether they are: +- Globally identifiable. +- Locally identifiable. +- 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 +using Catalyst, Logging, StructuralIdentifiability +goodwind_oscillator = @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) +``` +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) +``` +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) +``` +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. + +To, in a similar manner, indicate that certain initial conditions are known is a work in progress. Hopefully this feature should be an available in the near future. + +### 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 + (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) +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) +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) +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) +``` +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 state 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) +``` +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. + +`find_identifiable_functions` tries to simplify its output functions to create nice expressions. The degree to which it does this can be adjusted using the `simplify` keywords. Using the `:weak`, `:standard` (default), and `:strong` arguments, increased simplification can be forced (at the expense of longer runtime). + +## 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]) +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 +print_for_DAISY(si_ode) +nothing # hide +``` + +## Notes on systems with conservation laws +Several reaction network models, such as +```example si2 +using Catalyst, Logging, StructuralIdentifiability # hide +rs = @reaction_network begin + (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 +conservedequations(rs) +``` +E.g. if you run: +```example si2 +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 +find_identifiable_functions(rs; measured_quantities = [:X1, :X2], remove_conserved = false) +``` + +## Systems with exponent parameters +Structural identifiability cannot currently be applied to systems with parameters (or species) in exponents. E.g. this +```julia +rn = @reaction_network begin + (hill(X,v,K,n),d), 0 <--> X +end +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). + +--- +## [Citation](@id structural_identifiability_citation) +If you use this functionality in your research, please cite the following paper to support the authors of the StructuralIdentifiability package: +``` +@article{structidjl, + author = {Dong, R. and Goodbrake, C. and Harrington, H. and Pogudin G.}, + title = {Differential Elimination for Dynamical Models via Projections with Applications to Structural Identifiability}, + journal = {SIAM Journal on Applied Algebra and Geometry}, + url = {https://doi.org/10.1137/22M1469067}, + year = {2023} + volume = {7}, + number = {1}, + pages = {194-235} +} +``` + +--- +## References +[^1]: [Guillaume H.A. Joseph et al., *Introductory overview of identifiability analysis: A guide to evaluating whether you have the right type of data for your modeling purpose*, Environmental Modelling & Software (2019).](https://www.sciencedirect.com/science/article/pii/S1364815218307278) +[^2]: [Goodwin B.C., *Oscillatory Behavior in Enzymatic Control Processes*, Advances in Enzyme Regulation (1965).](https://www.sciencedirect.com/science/article/pii/0065257165900671?via%3Dihub) +[^3]: [Bellu G., et al., *DAISY: A new software tool to test global identifiability of biological and physiological systems*, Computer Methods and Programs in Biomedicine (2007).](https://www.sciencedirect.com/science/article/abs/pii/S0169260707001605) \ No newline at end of file diff --git a/ext/CatalystStructuralIdentifiabilityExtension.jl b/ext/CatalystStructuralIdentifiabilityExtension.jl new file mode 100644 index 0000000000..026fbe6122 --- /dev/null +++ b/ext/CatalystStructuralIdentifiabilityExtension.jl @@ -0,0 +1,10 @@ +module CatalystStructuralIdentifiabilityExtension + +# Fetch packages. +using Catalyst +import StructuralIdentifiability as SI + +# Creates and exports hc_steady_states function. +include("CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl") + +end \ No newline at end of file diff --git a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl new file mode 100644 index 0000000000..7d99ecc59a --- /dev/null +++ b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl @@ -0,0 +1,218 @@ +### Structural Identifiability ODE Creation ### + +# For a reaction system, list of measured quantities and known parameters, generate a StructuralIdentifiability compatible ODE. +""" +make_si_ode(rs::ReactionSystem; measured_quantities=observed(rs), known_p = [], ignore_no_measured_warn=false) + +Creates a reaction rate equation ODE system of the form used within the StructuralIdentifiability.jl package. The output system is compatible with all StructuralIdentifiability functions. + +Arguments: +- `rs::ReactionSystem`; The reaction system we wish to convert to an ODE. +- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.x`, or the symbol `:x`). +- `known_p = []`: List of parameters for which their values are assumed to be known. +- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. +- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). + +Example: +```julia +using Catalyst, StructuralIdentifiability +rs = @reaction_network begin + (p,d), 0 <--> X +end +make_si_ode(rs; measured_quantities = [:X], known_p = [:p]) +``` + +Notes: +- This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it. +- `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X]) +""" +function Catalyst.make_si_ode(rs::ReactionSystem; measured_quantities = [], known_p = [], + ignore_no_measured_warn = false, remove_conserved = true) + # Creates a MTK ODESystem, and a list of measured quantities (there are equations). + # Gives these to SI to create an SI ode model of its preferred form. + osys, conseqs, _ = make_osys(rs; remove_conserved) + measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn) + return SI.mtk_to_si(osys, measured_quantities)[1] +end + +### Structural Identifiability Wrappers ### + +""" +assess_local_identifiability(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], remove_conserved = true, ignore_no_measured_warn=false, kwargs...) + +Applies StructuralIdentifiability.jl's `assess_local_identifiability` function to the reaction rate equation ODE model for the given Catalyst `ReactionSystem`. Automatically converts the `ReactionSystem` to an appropriate `ODESystem` and computes structural identifiability, + +Arguments: +- `rs::ReactionSystem`; The reaction system for which we wish to compute structural identifiability of the associated reaction rate equation ODE model. +- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). +- `known_p = []`: List of parameters which values are known. +- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. +- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). + +Example: +```julia +using Catalyst, StructuralIdentifiability +rs = @reaction_network begin + (p,d), 0 <--> X +end +assess_local_identifiability(rs; measured_quantities = [:X], known_p = [:p]) +``` + +Notes: +- This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it. +- `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X]) +""" +function SI.assess_local_identifiability(rs::ReactionSystem, args...; measured_quantities = [], + known_p = [], funcs_to_check = Vector(), remove_conserved = true, + ignore_no_measured_warn=false, prob_threshold = 0.99, kwargs...) + # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form. + osys, conseqs, vars = make_osys(rs; remove_conserved) + measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn) + funcs_to_check = make_ftc(funcs_to_check, conseqs, vars) + + # Computes identifiability and converts it to a easy to read form. + out = SI.assess_local_identifiability(osys, args...; measured_quantities, funcs_to_check, p=prob_threshold, kwargs...) + return make_output(out, funcs_to_check, reverse.(conseqs)) +end + +""" +assess_identifiability(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], remove_conserved = true, ignore_no_measured_warn=false, kwargs...) + +Applies StructuralIdentifiability.jl's `assess_identifiability` function to a Catalyst `ReactionSystem`. Internally it is converted ot a `ODESystem`, for which structural identifiability is computed. + +Arguments: +- `rs::ReactionSystem`; The reaction system we wish to compute structural identifiability for. +- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). +- `known_p = []`: List of parameters which values are known. +- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. +- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). + +Example: +```julia +using Catalyst, StructuralIdentifiability +rs = @reaction_network begin + (p,d), 0 <--> X +end +assess_identifiability(rs; measured_quantities = [:X], known_p = [:p]) +``` + +Notes: +- This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it. +- `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X]) +""" +function SI.assess_identifiability(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], + funcs_to_check = Vector(), remove_conserved = true, ignore_no_measured_warn=false, + prob_threshold = 0.99, kwargs...) + # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form. + osys, conseqs, vars = make_osys(rs; remove_conserved) + measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn) + funcs_to_check = make_ftc(funcs_to_check, conseqs, vars) + + # Computes identifiability and converts it to a easy to read form. + out = SI.assess_identifiability(osys, args...; measured_quantities, funcs_to_check, p=prob_threshold, kwargs...) + return make_output(out, funcs_to_check, reverse.(conseqs)) +end + +""" +find_identifiable_functions(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], remove_conserved = true, ignore_no_measured_warn=false, kwargs...) + +Applies StructuralIdentifiability.jl's `find_identifiable_functions` function to a Catalyst `ReactionSystem`. Internally it is converted ot a `ODESystem`, for which structurally identifiable functions are computed. + +Arguments: +- `rs::ReactionSystem`; The reaction system we wish to compute structural identifiability for. +- `measured_quantities=[]`: The quantities of the system we can measure. May either be equations (e.g. `x1 + x2`), or single species (e.g. the symbolic `x`, `rs.s`, or the symbol `:x`). +- `known_p = []`: List of parameters which values are known. +- `ignore_no_measured_warn = false`: If set to `true`, no warning is provided when the `measured_quantities` vector is empty. +- `remove_conserved = true`: Whether to eliminate conservation laws when computing the ode (this can reduce runtime of identifiability analysis significantly). + +Example: +```julia +using Catalyst, StructuralIdentifiability +rs = @reaction_network begin + (p,d), 0 <--> X +end +find_identifiable_functions(rs; measured_quantities = [:X], known_p = [:p]) +``` + +Notes: +- This function is part of the StructuralIdentifiability.jl extension. StructuralIdentifiability.jl must be imported to access it. +- `measured_quantities` and `known_p` input may also be symbolic (e.g. measured_quantities = [rs.X]) +""" +function SI.find_identifiable_functions(rs::ReactionSystem, args...; measured_quantities = [], + known_p = [], remove_conserved = true, ignore_no_measured_warn=false, + prob_threshold = 0.99, kwargs...) + # Creates a ODESystem, and list of measured quantities, of SI's preferred form. + osys, conseqs = make_osys(rs; remove_conserved) + measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn) + + # Computes identifiable functions and converts it to a easy to read form. + out = SI.find_identifiable_functions(osys, args...; measured_quantities, p=prob_threshold, kwargs...) + return vector_subs(out, reverse.(conseqs)) +end + +### Helper Functions ### + +# From a reaction system, creates the corresponding MTK-style ODESystem for SI application +# Also compute the, later needed, conservation law equations and list of system symbols (states and parameters). +function make_osys(rs::ReactionSystem; remove_conserved=true) + # Creates the ODESystem corresponding to the ReactionSystem (expanding functions and flattening it). + # Creates a list of the systems all symbols (states and parameters). + rs = Catalyst.expand_registered_functions(flatten(rs)) + osys = convert(ODESystem, rs; remove_conserved) + vars = [states(rs); parameters(rs)] + + # Computes equations for system conservation laws. + # If there are no conserved equations, the `conseqs` variable must still have the `Vector{Pair{Any, Any}}` type. + if !remove_conserved + conseqs = Vector{Pair{Any, Any}}[] + else + conseqs = [ceq.lhs => ceq.rhs for ceq in conservedequations(rs)] + isempty(conseqs) && (conseqs = Vector{Pair{Any, Any}}[]) + end + + return osys, conseqs, vars +end + +# Creates a list of measured quantities of a form that SI can read. +# Each measured quantity must have a form like: +# `obs_var ~ X` # (Here, `obs_var` is a variable, and X is whatever we can measure). +function make_measured_quantities(rs::ReactionSystem, measured_quantities::Vector{T}, known_p::Vector{S}, + conseqs; ignore_no_measured_warn=false) where {T,S} + # Warning if the user didn't give any measured quantities. + if ignore_no_measured_warn || isempty(measured_quantities) + @warn "No measured quantity provided to the `measured_quantities` argument, any further identifiability analysis will likely fail. You can disable this warning by setting `ignore_no_measured_warn=true`." + end + + # Appends the known parameters to the measured_quantities vector. Converts any Symbols to symbolics. + mqiterator = Iterators.flatten((measured_quantities, known_p)) + mqs = [(q isa Symbol) ? Catalyst._symbol_to_var(rs, q) : q for q in mqiterator] + mqs = vector_subs(mqs, conseqs) + + # Creates one internal observation variable for each measured quantity (`___internal_observables`). + # Creates a vector of equations, setting each measured quantity equal to one observation variable. + @variables t (___internal_observables(Catalyst.get_iv(rs)))[1:length(mqs)] + return Equation[(q isa Equation) ? q : (___internal_observables[i] ~ q) for (i,q) in enumerate(mqs)] +end + +# Creates the functions that we wish to check for identifiability. +# If no `funcs_to_check` are given, defaults to checking identifiability for all states and parameters. +# Also, for conserved equations, replaces these in (creating a system without conserved quantities). +# E.g. for `X1 <--> X2`, replaces `X2` with `Γ[1] - X2`. +# Removing conserved quantities makes SI's algorithms much more performant. +function make_ftc(funcs_to_check, conseqs, vars) + isempty(funcs_to_check) && (funcs_to_check = vars) + return vector_subs(funcs_to_check, conseqs) +end + +# Processes the outputs to a better form. +# Replaces conservation law equations back in the output (so that e.g. Γ are not displayed). +# Sorts the output according to their input order (defaults to the `[states; parameters]` order). +function make_output(out, funcs_to_check, conseqs) + funcs_to_check = vector_subs(funcs_to_check, conseqs) + out = Dict(zip(vector_subs(keys(out), conseqs), values(out))) + sortdict = Dict(ftc => i for (i,ftc) in enumerate(funcs_to_check)) + return sort(out; by = x -> sortdict[x]) +end + +# For a vector of expressions and a conservation law, substitutes the law into every equation. +vector_subs(eqs, subs) = [substitute(eq, subs) for eq in eqs] \ No newline at end of file diff --git a/src/Catalyst.jl b/src/Catalyst.jl index fa2b3e51c2..233b2274bf 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -112,6 +112,10 @@ export balance_reaction function hc_steady_states end export hc_steady_states +# StructuralIdentifiability +function make_si_ode end +export make_si_ode + ### Spatial Reaction Networks ### # spatial reactions diff --git a/test/extensions/structural_identifiability.jl b/test/extensions/structural_identifiability.jl new file mode 100644 index 0000000000..683203e134 --- /dev/null +++ b/test/extensions/structural_identifiability.jl @@ -0,0 +1,322 @@ +### Fetch Packages ### + +using Catalyst, Test +using StructuralIdentifiability + + +### Helper Function ### + +# Converts the output dicts from StructuralIdentifiability functions from "weird symbol => stuff" to "symbol => stuff" (the output have some strange meta data which prevents equality checks, this enables this). +# Structural identifiability also provides variables like x (rather than x(t)). This is a bug, but we have to convert to make it work (now just remove any (t) to make them all equal). +function sym_dict(dict_in) + dict_out = Dict{Symbol,Any}() + for key in keys(dict_in) + sym_key = Symbol(key) + sym_key = Symbol(replace(String(sym_key), "(t)" => "")) + dict_out[sym_key] = dict_in[key] + end + return dict_out +end + + +### Run Tests ### + +# Tests for Goodwin model (model with both global, local, and non identifiable components). +# Tests for system using Catalyst function (in this case, Michaelis-Menten function) +let + # Identifiability analysis for Catalyst model. + goodwind_oscillator_catalyst = @reaction_network begin + (mmr(P,pₘ,1), dₘ), 0 <--> M + (pₑ*M,dₑ), 0 <--> E + (pₚ*E,dₚ), 0 <--> P + end + gi_1 = assess_identifiability(goodwind_oscillator_catalyst; measured_quantities=[:M]) + li_1 = assess_local_identifiability(goodwind_oscillator_catalyst; measured_quantities=[:M]) + ifs_1 = find_identifiable_functions(goodwind_oscillator_catalyst; measured_quantities=[:M]) + + # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model. + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M]) + gi_2 = assess_identifiability(si_catalyst_ode) + li_2 = assess_local_identifiability(si_catalyst_ode) + ifs_2 = find_identifiable_functions(si_catalyst_ode) + + # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.). + goodwind_oscillator_si = @ODEmodel( + M'(t) = pₘ / (1 + P(t)) - dₘ*M(t), + E'(t) = -dₑ*E(t) + pₑ*M(t), + P'(t) = -dₚ*P(t) + pₚ*E(t), + y1(t) = M(t) + ) + gi_3 = assess_identifiability(goodwind_oscillator_si) + li_3 = assess_local_identifiability(goodwind_oscillator_si) + ifs_3 = find_identifiable_functions(goodwind_oscillator_si) + + # Check outputs. + @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3) + @test sym_dict(li_1) == sym_dict(li_2) == sym_dict(li_3) + @test length(ifs_1) == length(ifs_2) == length(ifs_3) + + # Checks output to manually checked correct answers. + @test isequal(collect(keys(gi_1)), [states(goodwind_oscillator_catalyst); parameters(goodwind_oscillator_catalyst)]) + @test isequal(collect(values(gi_1)), [:globally, :nonidentifiable, :globally, :globally, :globally, :nonidentifiable, :locally, :nonidentifiable, :locally]) + @test isequal(collect(keys(li_1)), [states(goodwind_oscillator_catalyst); parameters(goodwind_oscillator_catalyst)]) + @test isequal(collect(values(li_1)), [1, 0, 1, 1, 1, 0, 1, 0, 1]) +end + +# Tests on a made-up reaction network with mix of identifiable and non-identifiable components. +# Tests for symbolics input. +# Tests using known_p argument. +let + # Identifiability analysis for Catalyst model. + rs_catalyst = @reaction_network begin + (p1, d), 0 <--> X1 + k1, X1 --> X2 + (k2f,k2b), X2 <--> X3 + k3, X3 --> X4 + d, X4 --> 0 + end + @unpack X2, X3 = rs_catalyst + gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f]) + li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f]) + ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f]) + + # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model. + rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f]) + gi_2 = assess_identifiability(rs_ode) + li_2 = assess_local_identifiability(rs_ode) + ifs_2 = find_identifiable_functions(rs_ode) + + # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.). + rs_si = @ODEmodel( + X1'(t) = p1 - d*X1(t) - k1*X1(t), + X2'(t) = k1*X1(t) + k2b*X3(t) - k2f*X2(t), + X3'(t) = -k2b*X3(t) + k2f*X2(t) - k3*X3(t), + X4'(t) = d*X4(t) + k3*X3(t), + y1(t) = X2, + y2(t) = X3, + y3(t) = k2f + ) + gi_3 = assess_identifiability(rs_si) + li_3 = assess_local_identifiability(rs_si) + ifs_3 = find_identifiable_functions(rs_si) + + # Check outputs. + @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3) + @test sym_dict(li_1) == sym_dict(li_2) == sym_dict(li_3) + @test length(ifs_1) == length(ifs_2) == length(ifs_3) + + # Checks output to manually checked correct answers. + @test isequal(collect(keys(gi_1)),[states(rs_catalyst); parameters(rs_catalyst)]) + @test isequal(collect(values(gi_1)),[:nonidentifiable, :globally, :globally, :nonidentifiable, :nonidentifiable, :nonidentifiable, :nonidentifiable, :globally, :globally, :globally]) + @test isequal(collect(keys(li_1)),[states(rs_catalyst); parameters(rs_catalyst)]) + @test isequal(collect(values(li_1)),[0, 1, 1, 0, 0, 0, 0, 1, 1, 1]) +end + +# Tests on a made-up reaction network with mix of identifiable and non-identifiable components. +# Tests for system with conserved quantity. +# Tests for symbolics known_p +# Tests using an equation for measured quantity. +let + # Identifiability analysis for Catalyst model. + rs_catalyst = @reaction_network begin + p, 0 --> X1 + k1, X1 --> X2 + k2, X2 --> X3 + k3, X3 --> X4 + k3, X3 --> X5 + d, (X4,X5) --> 0 + (kA*X3, kD), Yi <--> Ya + end + @unpack X1, X2, X3, X4, k1, k2, Yi, Ya, k1, kD = rs_catalyst + gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD]) + li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD]) + ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD]) + + # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model. + rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD], remove_conserved=false) + gi_2 = assess_identifiability(rs_ode) + li_2 = assess_local_identifiability(rs_ode) + ifs_2 = find_identifiable_functions(rs_ode) + + # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.). + rs_si = @ODEmodel( + X1'(t) = p - k1*X1(t), + X2'(t) = k1*X1(t) - k2*X2(t), + X3'(t) = k2*X2(t) - 2k3*X3(t), + X4'(t) = -d*X4(t) + k3*X3(t), + X5'(t) = -d*X5(t) + k3*X3(t), + Yi'(t) = kD*Ya(t) - kA*Yi(t)*X3(t), + Ya'(t) = -kD*Ya(t) + kA*Yi(t)*X3(t), + y1(t) = X1 + Yi, + y2(t) = Ya, + y3(t) = k1, + y4(t) = kD + ) + gi_3 = assess_identifiability(rs_si) + li_3 = assess_local_identifiability(rs_si) + ifs_3 = find_identifiable_functions(rs_si) + + # Check outputs. + @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3) + @test sym_dict(li_1) == sym_dict(li_2) == sym_dict(li_3) + @test length(ifs_1[2:end]) == length(ifs_2) == length(ifs_3) # In the first case, the conservation law parameter is also identifiable. + + # Checks output to manually checked correct answers. + correct_gi = Pair.([states(rs_catalyst); parameters(rs_catalyst)], [:globally, :locally, :locally, :nonidentifiable, :nonidentifiable, :globally, :globally, :globally, :globally, :locally, :locally, :nonidentifiable, :locally, :globally]) + correct_li = Pair.([states(rs_catalyst); parameters(rs_catalyst)], [1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1]) + @test issetequal(gi_1, correct_gi) + @test issetequal(li_1, correct_li) +end + +# Tests that various inputs types work. +let + goodwind_oscillator_catalyst = @reaction_network begin + (mmr(P,pₘ,1), dₘ), 0 <--> M + (pₑ*M,dₑ), 0 <--> E + (pₚ*E,dₚ), 0 <--> P + end + @unpack M, E, P, pₑ, pₚ, pₘ = goodwind_oscillator_catalyst + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p=[:pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M], known_p=[:pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M, :E], known_p=[:pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M], known_p=[:pₑ, :pₚ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M, :E], known_p=[:pₑ, :pₚ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p=[pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M], known_p=[pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M, E], known_p=[pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M], known_p=[pₑ, pₚ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M, E], known_p=[pₑ, pₚ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M + pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M + E, pₑ*M], known_p=[:pₑ]) + si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[pₑ, pₚ], known_p=[pₑ]) + + # Tests using model.component style (have to make system complete first). + gw_osc_complt = complete(goodwind_oscillator_catalyst) + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M]) isa ODE + @test make_si_ode(gw_osc_complt; known_p=[gw_osc_complt.pₑ]) isa ODE + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p=[gw_osc_complt.pₑ]) isa ODE + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M, gw_osc_complt.E], known_p=[gw_osc_complt.pₑ]) isa ODE + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p=[gw_osc_complt.pₑ, gw_osc_complt.pₚ]) isa ODE + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p = [:pₚ]) isa ODE + @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M*gw_osc_complt.E]) isa ODE +end + +# Check that `prob_threshold` alternative kwarg works. +let + rs = @reaction_network begin + p, X --> 0 + end + @unpack X = rs + + assess_identifiability(rs; measured_quantities=[X], prob_threshold=0.9) + assess_identifiability(rs; measured_quantities=[X], prob_threshold=0.999) +end + +# Tests for hierarchical model with conservation laws at both top and internal levels. +let + # Identifiability analysis for Catalyst model. + rs1 = @reaction_network rs1 begin + (k1, k2), X1 <--> X2 + end + rs2 = @reaction_network rs2 begin + (k3, k4), X3 <--> X4 + end + @named rs_catalyst = compose(rs1, [rs2]) + @unpack X1, X2, k1, k2 = rs1 + gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1]) + li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1]) + ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1]) + + # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model. + rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1]) + gi_2 = assess_identifiability(rs_ode) + li_2 = assess_local_identifiability(rs_ode) + ifs_2 = find_identifiable_functions(rs_ode) + + # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.). + rs_si = @ODEmodel( + X1'(t) = -k1*X1(t) + k2*X2(t), + X2'(t) = k1*X1(t) - k2*X2(t), + rs2₊X3'(t) = -rs2₊k3*rs2₊X3(t) + rs2₊k4*rs2₊X4(t), + rs2₊X4'(t) = rs2₊k3*rs2₊X3(t) - rs2₊k4*rs2₊X4(t), + y1(t) = X1, + y2(t) = X2, + y3(t) = rs2₊X3, + y4(t) = k1 + ) + gi_3 = assess_identifiability(rs_si) + li_3 = assess_local_identifiability(rs_si) + ifs_3 = find_identifiable_functions(rs_si) + + # Check outputs. + @test sym_dict(gi_1) == sym_dict(gi_3) + @test sym_dict(li_1) == sym_dict(li_3) + @test length(ifs_1)-2 == length(ifs_2)-2 == length(ifs_3) # In the first case, the conservation law parameter is also identifiable. + + # Checks output for the SI converted version of the catalyst model. + # For nested systems with conservation laws, conserved quantities like Γ[1], cannot be replaced back. + # Hence, here you display identifiability for `Γ[1]` instead of X2. + gi_1_no_cq = filter(x -> !occursin("X2",String(x[1])) && !occursin("X4",String(x[1])), sym_dict(gi_1)) + gi_2_no_cq = filter(x -> !occursin("Γ",String(x[1])), sym_dict(gi_2)) + li_1_no_cq = filter(x -> !occursin("X2",String(x[1])) && !occursin("X4",String(x[1])), sym_dict(li_1)) + li_2_no_cq = filter(x -> !occursin("Γ",String(x[1])), sym_dict(li_2)) + @test gi_1_no_cq == gi_2_no_cq + @test li_1_no_cq == li_2_no_cq +end + +# Tests directly on reaction systems with known identifiability structures. +# Test provided by Alexander Demin. +let + rs = @reaction_network begin + k1, x1 --> x2 + end + # Measure the source + id_report = assess_identifiability(rs, measured_quantities = [:x1]) + @test sym_dict(id_report) == Dict( + :x1 => :globally, + :x2 => :nonidentifiable, + :k1 => :globally + ) + # Measure the target instead + id_report = assess_identifiability(rs, measured_quantities = [:x2]) + @test sym_dict(id_report) == Dict( + :x1 => :globally, + :x2 => :globally, + :k1 => :globally + ) + + # Example from + # Identifiability of chemical reaction networks + # DOI: 10.1007/s10910-007-9307-x + # The rate constants a, b, c are not identifiable even if all of the species + # are observed. + rs = @reaction_network begin + a, A0 --> 2A1 + b, A0 --> 2A2 + c, A0 --> A1 + A2 + end + id_report = assess_identifiability(rs, measured_quantities = [:A0, :A1, :A2]) + @test sym_dict(id_report) == Dict( + :A0 => :globally, + :A1 => :globally, + :A2 => :globally, + :a => :nonidentifiable, + :b => :nonidentifiable, + :c => :nonidentifiable + ) + + # Test with no parameters + rs = @reaction_network begin + 1, x1 --> x2 + 1, x2 --> x3 + end + id_report = assess_identifiability(rs, measured_quantities = [:x3]) + @test sym_dict(id_report) == Dict( + :x1 => :globally, + :x2 => :globally, + :x3 => :globally, + ) + @test length(find_identifiable_functions(rs, measured_quantities = [:x3])) == 1 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index e855dfd3d8..cd2b6550e3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,5 +54,6 @@ using SafeTestsets ### Tests extensions. ### @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end + @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end end # @time