Structural Identifiability Extension

[TorkelE]

Creates an extension for Catalyst and the StructuralIdentifiability package (https://github.com/SciML/StructuralIdentifiability.jl).

Implements:

• dispatch for calling assess_local_identifiability on ReactionSystems
• dispatch for calling assess_identifiability on ReactionSystems
• dispatch for calling find_identifiable_functions on ReactionSystems
• The make_si_ode function, which takes a Catalyst ReactionSystem and creates a ODE model of the form used in StructuralIdentifiability.jl, that can be used for some of its more detailed features.

Example usage:

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])

here, we assess global identifiability for the model, given us being able to measure M.

[TorkelE]

THIS PR is now ready. There are two outstanding issues:

• Handling parameters in exponents (like occurs in e.g. Hill functions). This will be complicated, and require some major work. It will have to be a future update.
• • Handling conservation laws to speed up computation. Hopefully some additional functionality will come of in StructuralIdentfiability to simplify this. What to do about this depend on that.

[sumiya11]

Hi Torkel! Thanks again for working on this !

I think the tutorial looks very good.

I wanted to add a couple more simple tests but it looks like I cannot push to this PR. So I will just paste it here, feel free to add these to tests:

# Test some examples explicitly
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,
    )
    # Will probably be fixed in the 0.5 release of SI.jl
    @test_broken find_identifiable_functions(rs, measured_quantities = [:x3])
end

[sumiya11]

A couple of papers with examples that can be used for testing:

• https://www.pnas.org/doi/10.1073/pnas.1215508110. There are also some examples in the supporting information
• https://link.springer.com/article/10.1007/s10910-007-9307-x

[isaacsas]

@TorkelE does this work with Symbolics instead of symbols?

I'm concerned that across these PRs you are adding you are making the interface rely on Symbols all over the place. This is likely to be a major break from ModelingToolkit, and could lead us to have to change all this code down the line.

@ChrisRackauckas can you comment on this? Both the PETab intergration and now the structural identifiability extension are written to use symbols instead of symbolics as the workflow. Is this going to cause issues / inconsistencies with MTK?

[isaacsas]

Actually, I think the BifurcationKit extension also uses symbols instead of symbolics right?

[TorkelE]

@sumiya11 Thanks a lot for the test cases and input, I will inccoperate this.

@isaacsas Both this, the BifurcationKit PR, and the rest, can use both symbols and symbolic, and the tests covers both cases.

[isaacsas]

Great!

[sumiya11]

Handling conservation laws to speed up computation. Hopefully some additional functionality will come of in StructuralIdentfiability to simplify this. What to do about this depend on that.

Which functionality in SI would you need?

Note also, it is tricky to use conservation laws together with find_identifiable_functions (since if you eliminate x1 by x1 = C - x2, you would need to substitute x1 back after the analysis is done and probably simplify the expressions once again).

So, one feasible thing would be to do this in StructuralIdentfiability instead of Catalyst. There is some ongoing work on this

[TorkelE]

Solving this: SciML/StructuralIdentifiability.jl#236 should make it relatively easy to enable conservation laws. There is also some discussion in the Issue on alternative ways to do it, but it will be less straightforward, and I'm still not really sure how.

Yes, handling it all in StructuralIdentfiability would be preferable, since it could cover a wide range of cases. The advantage right now in Catalyst is that we have algorithms implemented that finds these laws on reaction networks, which we then can apply, while SI would need to implement more general methods (but which then would benefit all systems). Long term that would be preferable, but it might also take more time to get everything into place.

[TorkelE]

A question for everyone involved:
StructuralIdentifiabiltiy.jl currently return an output on the form of a dictionary taking components to identifiability. E.h.

:X => :globally
:p => :globally
:d => :globally

In Catalyst, the various components have a natural order (initial conditions followed by parameters, each with their own separate order). Currently, when identifability is printed, this can appear quite jumbled, e.g.:

:d => :globally
:X => :globally
:p => :globally

(mixing parameters and initial conditions etc.)

I think it would be nicer to impose order, so that when the user prints the entire identifiability dictionary, they are ordered initial conditions followed by parameters.

This would require the making the output an OrderedDict. This while the output of SI functions are simply Dicts. It seems undesirable that e.g. assess_identifiaility has a different output type when applied to a Catalyst ReactionSystem (as opposed to a SI type ODE). How bad do we think this is? Previously I decided not to make the output an OrderedDict, but there are advantages. Do anyone have any input?

[sumiya11]

I think it would be nicer to impose order, so that when the user prints the entire identifiability dictionary, they are ordered initial conditions followed by parameters.

I agree, some order can be helpful, especially if there are a lot of quantities. I opened an issue:
SciML/StructuralIdentifiability.jl#244

[pogudingleb]

A question for everyone involved: StructuralIdentifiabiltiy.jl currently return an output on the form of a dictionary taking components to identifiability. E.h.

:X => :globally
:p => :globally
:d => :globally

In Catalyst, the various components have a natural order (initial conditions followed by parameters, each with their own separate order). Currently, when identifability is printed, this can appear quite jumbled, e.g.:

:d => :globally
:X => :globally
:p => :globally

(mixing parameters and initial conditions etc.)

I think it would be nicer to impose order, so that when the user prints the entire identifiability dictionary, they are ordered initial conditions followed by parameters.

This would require the making the output an OrderedDict. This while the output of SI functions are simply Dicts. It seems undesirable that e.g. assess_identifiaility has a different output type when applied to a Catalyst ReactionSystem (as opposed to a SI type ODE). How bad do we think this is? Previously I decided not to make the output an OrderedDict, but there are advantages. Do anyone have any input?

This is fixed in the current dev version of StructuralIdentifiability: if funcs_to_check are not provided, states are listed first, otherwise, the results are listed in the same order as in funcs_to_check.

[TorkelE]

Thanks! This will make things much neater

[TorkelE]

now rebased on master