From 6879390f5760946885fa2a0d3ab6687a3b298cc7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 27 Sep 2024 03:46:45 -0400 Subject: [PATCH 1/4] fix BK and HC tests --- test/extensions/bifurcation_kit.jl | 21 ++++++++++++--------- test/extensions/homotopy_continuation.jl | 4 ++-- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl index 7900ed436b..13d090ded4 100644 --- a/test/extensions/bifurcation_kit.jl +++ b/test/extensions/bifurcation_kit.jl @@ -56,32 +56,35 @@ end # Checks that the same bifurcation problem is created as for BifurcationKit. # Checks with Symbolics as bifurcation and plot vars. # Tries setting `jac=false`. +# Note: Only one parameter used, as tests technically depended on internal parameter ordering +# (Potentially the test should also be removed as it tests internal parameter stuff, however, test +# was written a while ago when I paid less attention to this kind of stuff.) let # Creates BifurcationProblem via Catalyst. bistable_switch = @reaction_network begin - 0.1 + hill(X,v,K,n), 0 --> X - d, X --> 0 + 0.1 + hill(X,5.0,K,3), 0 --> X + 1.0, X --> 0 end - @unpack X, v, K, n, d = bistable_switch + @unpack X, K = bistable_switch u0_guess = [X => 1.0] - p_start = [v => 5.0, K => 2.5, n => 3, d => 1.0] + p_start = [K => 2.5] bprob = BifurcationProblem(bistable_switch, u0_guess, p_start, K; jac=false, plot_var=X) # Creates BifurcationProblem via BifurcationKit. function bistable_switch_BK(u, p) X, = u - v, K, n, d = p - return [0.1 + v*(X^n)/(X^n + K^n) - d*X] + K, = p + return [0.1 + 5.0*(X^3)/(X^3 + K^3) - 1.0*X] end - bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [5.0, 2.5, 3, 1.0], (@lens _[1]); record_from_solution = (x, p) -> x[1]) + bprob_BK = BifurcationProblem(bistable_switch_BK, [1.0], [2.5], (@lens _[1]); record_from_solution = (x, p) -> x[1]) # Check the same function have been generated. bprob.u0 == bprob_BK.u0 bprob.params == bprob_BK.params for repeat = 1:20 u0 = rand(rng, 1) - p = rand(rng, 4) - @test bprob_BK.VF.F(u0, p) == bprob.VF.F(u0, p) + p = rand(rng, 1) + @test bprob_BK.VF.F(u0, p) ≈ bprob.VF.F(u0, p) end end diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 3a59f0f3a8..1cb31a1bfc 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -96,12 +96,12 @@ let v*(0.1/v + hill(X,1,K,n)), 0 --> X d, X --> 0 end - ps = [:v => 5.0, :K => 2.5, :n => 3, :d => 1.0] + ps = Dict([:v => 5.0, :K => 2.5, :n => 3, :d => 1.0]) sss = hc_steady_states(rs, ps; filter_negative = false, show_progress = false, seed = 0x000004d1) @test length(sss) == 4 for ss in sss - @test f_eval(rs,sss[1], last.(ps), 0.0)[1] ≈ 0.0 atol = 1e-12 + @test ps[:v]*(0.1/ps[:v] + ss[1]^ps[:n]/(ss[1]^ps[:n] + ps[:K]^ps[:n])) - ps[:d]*ss[1]≈ 0.0 atol = 1e-12 end ps = [:v => 5.0, :K => 2.5, :n => 2.7, :d => 1.0] From f237e8d0b0bd451b0b98c4500def2be7789163a2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 27 Sep 2024 03:47:20 -0400 Subject: [PATCH 2/4] up --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6ad68bf0b6..fd4e69c486 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,7 +60,7 @@ using SafeTestsets, Test # 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 + # @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end # Tests stability computation (uses HomotopyContinuation extension). @time @safetestset "Steady State Stability Computations" begin include("miscellaneous_tests/stability_computation.jl") end From 5e2812465d17caf6f3bcda45ed4459ca2aa9e2c4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 29 Sep 2024 12:52:12 -0400 Subject: [PATCH 3/4] disable SI docs --- .../structural_identifiability.md | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/docs/src/inverse_problems/structural_identifiability.md b/docs/src/inverse_problems/structural_identifiability.md index 649f68ab40..0a8e7a93a7 100644 --- a/docs/src/inverse_problems/structural_identifiability.md +++ b/docs/src/inverse_problems/structural_identifiability.md @@ -1,4 +1,7 @@ # [Structural Identifiability Analysis](@id structural_identifiability) +!!! note + Due to recent changes in upstream packages, this Catalyst feature is currently broken on Catalyst v14. + 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. @@ -31,7 +34,7 @@ Global identifiability can be assessed using the `assess_identifiability` functi - Unidentifiable. To it, we provide our `ReactionSystem` model and a list of quantities that we are able to measure. Here, we consider a Goodwind oscillator (a simple 3-component model, where the three species $M$, $E$, and $P$ are produced and degraded, which may exhibit oscillations)[^2]. Let us say that we are able to measure the concentration of $M$, we then designate this using the `measured_quantities` argument. We can now assess identifiability in the following way: -```@example si1 +```julia using Catalyst, Logging, StructuralIdentifiability gwo = @reaction_network begin (pₘ/(1+P), dₘ), 0 <--> M @@ -43,14 +46,14 @@ assess_identifiability(gwo; measured_quantities = [:M], loglevel = Logging.Error From the output, we find that `E(t)`, `pₑ`, and `pₚ` (the trajectory of $E$, and the production rates of $E$ and $P$, respectively) are non-identifiable. Next, `dₑ` and `dₚ` (the degradation rates of $E$ and $P$, respectively) are locally identifiable. Finally, `P(t)`, `M(t)`, `pₘ`, and `dₘ` (the trajectories of `P` and `M`, and the production and degradation rate of `M`, respectively) are all globally identifiable. We note that we also imported the Logging.jl package, and provided the `loglevel = Logging.Error` input argument. StructuralIdentifiability functions generally provide a large number of output messages. Hence, we will use this argument (which requires the Logging package) throughout this tutorial to decrease the amount of printed text. Next, we also assess identifiability in the case where we can measure all three species concentrations: -```@example si1 +```julia assess_identifiability(gwo; measured_quantities = [:M, :P, :E], loglevel = Logging.Error) ``` in which case all species trajectories and parameters become identifiable. ### [Indicating known parameters](@id structural_identifiability_gi_known_ps) 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 +```julia assess_identifiability(gwo; measured_quantities = [:M], known_p = [:pₑ], loglevel = Logging.Error) ``` Not only does this turn the previously non-identifiable `pₑ` (globally) identifiable (which is obvious, given that its value is now known), but this additional information improve identifiability for several other network components. @@ -59,14 +62,14 @@ To, in a similar manner, indicate that certain initial conditions are known is a ### [Providing non-trivial measured quantities](@id structural_identifiability_gi_nontrivial_mq) 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 si1 +```julia rs = @reaction_network begin (kA,kD), Eᵢ <--> Eₐ (Eₐ, d), 0 <-->P end ``` If we can measure the total amount of $E$ ($=Eᵢ+Eₐ$), as well as the amount of $P$, we can use the following to assess identifiability: -```@example si1 +```julia @unpack Eᵢ, Eₐ = rs assess_identifiability(rs; measured_quantities = [Eᵢ + Eₐ, :P], loglevel = Logging.Error) nothing # hide @@ -74,7 +77,7 @@ nothing # hide ### [Assessing identifiability for specified quantities only](@id structural_identifiability_gi_ftc) 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 +```julia @unpack pₘ, pₑ, pₚ, M, E, P = gwo assess_identifiability(gwo; measured_quantities = [:M], funcs_to_check = [pₘ, pₑ, pₚ, M + E + P], loglevel = Logging.Error) nothing # hide @@ -82,7 +85,7 @@ nothing # hide ### [Probability of correctness](@id structural_identifiability_gi_probs) 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 si1 +```julia assess_identifiability(gwo; measured_quantities=[:M], prob_threshold = 0.999, loglevel = Logging.Error) nothing # hide ``` @@ -90,14 +93,14 @@ giving a minimum bound of $99.9\%$ chance of correctness. In practise, the bound ## [Local identifiability analysis](@id structural_identifiability_lit) 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 +```julia assess_local_identifiability(gwo; measured_quantities = [:M], loglevel = Logging.Error) ``` We note that the results are consistent with those produced by `assess_identifiability` (with globally or locally identifiable quantities here all being assessed as at least locally identifiable). ## [Finding identifiable functions](@id structural_identifiability_identifiabile_funcs) Finally, StructuralIdentifiability provides the `find_identifiable_functions` function. Rather than determining the identifiability of each parameter and unknown of the model, it finds a set of identifiable functions, such as any other identifiable expression of the model can be generated by these. Let us again consider the Goodwind oscillator, using the `find_identifiable_functions` function we find that identifiability can be reduced to five globally identifiable expressions: -```@example si1 +```julia find_identifiable_functions(gwo; measured_quantities = [:M], loglevel = Logging.Error) ``` Again, these results are consistent with those produced by `assess_identifiability`. There, `pₑ` and `pₚ` where found to be globally identifiable. Here, they correspond directly to identifiable expressions. The remaining four parameters (`pₘ`, `dₘ`, `dₑ`, and `dₚ`) occur as part of more complicated composite expressions. @@ -106,19 +109,19 @@ Again, these results are consistent with those produced by `assess_identifiabili ## [Creating StructuralIdentifiability compatible ODE models from Catalyst `ReactionSystem`s](@id structural_identifiability_si_odes) 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 +```julia si_ode = make_si_ode(gwo; measured_quantities = [:M]) nothing # hide ``` and then used as input to various StructuralIdentifiability functions. In the following example we use StructuralIdentifiability's `print_for_DAISY` function, printing the model as an expression that can be used by the [DAISY](https://daisy.dei.unipd.it/) software for identifiability analysis[^3]. -```@example si1 +```julia print_for_DAISY(si_ode) nothing # hide ``` ## [Notes on systems with conservation laws](@id structural_identifiability_conslaws) Several reaction network models, such as -```@example si3 +```julia using Catalyst, Logging, StructuralIdentifiability # hide rs = @reaction_network begin (k1,k2), X1 <--> X2 From 4ed564f7b6fd3538c228303ddb09d4b1508597cc Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 30 Sep 2024 11:33:40 -0400 Subject: [PATCH 4/4] add note in history file --- HISTORY.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HISTORY.md b/HISTORY.md index b0fa05c1e6..22d6aea5ca 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,7 +1,7 @@ # Breaking updates and feature summaries across releases ## Catalyst unreleased (master branch) - +- Due to changes in upstream packages, the structural identifiability extension is currently broken. ## Catalyst 14.4.1 - Support for user-defined functions on the RHS when providing coupled equations