# Using ModelingToolkit.jl With StructuralIdentifiability.jl

In this tutorial, we will cover examples of solving identifiability problems for models defined with the syntax of ModelingToolkit.jl.

## Input System

Let us consider the following ODE model with two outputs:

$$$\begin{cases} \dot{S} = -b \, S \, (I + J + q \, A) \, N_{inv},\\ \dot{E} = b \, S \, (I + J + q \, A) \, N_{inv} - k \, E,\\ \dot{A} = k \, (1 - r) \, E - g_1 \, A,\\ \dot{I} = k \, r \, E - (\alpha + g_1) \, I,\\ \dot{J} = \alpha \, I - g_2 \, J,\\ \dot{C} = \alpha \, I,\\ y_1 = C,\\ y_2 = N_{inv} \end{cases}$$$

This is an infectious desease model defined in [1].

The main difference between the input formats in ModelingToolkit.jl and StructuralIdentifiability.jl is that the output (measured values/functions) must be specified separately in ModelingToolkit.jl. In this example, measured quantities are presented by $y_1$, $y_2$.

First, let us define the ODE. We will use @parameters and @variables macro to define parameters and time-depended functions in the ODE.

using StructuralIdentifiability, ModelingToolkit

@parameters b q N_inv k r alpha g1 g2
@variables t S(t) E(t) A(t) I(t) J(t) C(t) y1(t) y2(t)

The actual ODE will be defined using ODESystem structure from ModelingToolkit.jl:

D = Differential(t)

eqs = [
D(S) ~ -b * S * (I + J + q * A) * N_inv,
D(E) ~ b * S * (I + J + q * A) * N_inv - k * E,
D(A) ~ k * (1 - r) * E - g1 * A,
D(I) ~ k * r * E - (alpha + g1) * I,
D(J) ~ alpha * I - g2 * J,
D(C) ~ alpha * I,
]

ode = ODESystem(eqs, t, name = :SEIAJRCmodel)

Finally, let us define the array of measured quantities and call the assess_identifiability function. This is the main function that determines local/global identifiability properties of each parameter and state. We will use the probability of correctness $p=0.99$.

For ModelingToolkit.jl, both assess_identifiability and assess_local_identifiability functions accept keyword arguments:

• measured_quantities, also called "output functions" in identifiability literature; these are crucial for answering identifiability questions.
• p, probability of correctness. This value equals 0.99 by default.
• funcs_to_check, functions of parameters of which we wish to check identifiability.
measured_quantities = [y1 ~ C, y2 ~ N_inv]
@time global_id = assess_identifiability(ode, measured_quantities=measured_quantities)

Let us put all of the code above together:

julia> using StructuralIdentifiability, ModelingToolkitjulia> @parameters b q N_inv k r alpha g1 g28-element Vector{Num}:
b
q
N_inv
k
r
alpha
g1
g2julia> @variables t S(t) E(t) A(t) I(t) J(t) C(t) y1(t) y2(t)9-element Vector{Num}:
t
S(t)
E(t)
A(t)
I(t)
J(t)
C(t)
y1(t)
y2(t)julia> D = Differential(t)(::Differential) (generic function with 2 methods)julia> eqs = [
D(S) ~ -b * S * (I + J + q * A) * N_inv,
D(E) ~ b * S * (I + J + q * A) * N_inv - k * E,
D(A) ~ k * (1 - r) * E - g1 * A,
D(I) ~ k * r * E - (alpha + g1) * I,
D(J) ~ alpha * I - g2 * J,
D(C) ~ alpha * I,
]6-element Vector{Equation}:
Differential(t)(S(t)) ~ -N_inv*b*(q*A(t) + I(t) + J(t))*S(t)
Differential(t)(E(t)) ~ N_inv*b*(q*A(t) + I(t) + J(t))*S(t) - k*E(t)
Differential(t)(A(t)) ~ k*(1 - r)*E(t) - g1*A(t)
Differential(t)(I(t)) ~ k*r*E(t) - (alpha + g1)*I(t)
Differential(t)(J(t)) ~ alpha*I(t) - g2*J(t)
Differential(t)(C(t)) ~ alpha*I(t)julia> ode = ODESystem(eqs, t, name = :SEIAJRCmodel)Model SEIAJRCmodel with 6 equations
States (6):
S(t)
E(t)
A(t)
I(t)
⋮
Parameters (8):
N_inv
q
b
k
⋮julia> measured_quantities = [y1 ~ C, y2 ~ N_inv]2-element Vector{Equation}:
y1(t) ~ C(t)
y2(t) ~ N_invjulia> @time global_id = assess_identifiability(ode, measured_quantities=measured_quantities)[ Info: Preproccessing ModelingToolkit.ODESystem object
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.1248456 seconds
[ Info: Assessing global identifiability
[ Info: Computing IO-equations
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Selected tablesize 512
[ Info: 1: selected lucky prime 2147483647
[ Info: CRT modulo (1, 2147483647)
[ Info: Reconstructing modulo 2147483647
[ Info: Heuristic check failed.
[ Info: 2: selected lucky prime 2147483659
[ Info: CRT modulo (2147483647, 2147483659)
[ Info: Reconstructing modulo 4611686039902224373
[ Info: Heuristic check failed.
[ Info: 3: selected lucky prime 2147483693
[ Info: CRT modulo (4611686039902224373, 2147483693)
[ Info: Reconstructing modulo 9903520567925774155444649489
[ Info: 4: selected lucky prime 2147483713
[ Info: CRT modulo (9903520567925774155444649489, 2147483713)
[ Info: Reconstructing modulo 21267649120981110191733715050621272657
[ Info: 5: selected lucky prime 2147483743
[ Info: CRT modulo (21267649120981110191733715050621272657, 2147483743)
[ Info: Reconstructing modulo 45671930739135174346839766056203605080877915151
[ Info: 6: selected lucky prime 2147483777
[ Info: CRT modulo (45671930739135174346839766056203605080877915151, 2147483777)
[ Info: Reconstructing modulo 98079730326560405919904968824172512120100095704355005327
[ Info: Heuristic check failed.
[ Info: 7: selected lucky prime 2147483783
[ Info: CRT modulo (98079730326560405919904968824172512120100095704355005327, 2147483783)
[ Info: 8: selected lucky prime 2147483813
[ Info: CRT modulo (210624630317301765882893117451031048172285903861850336414611112041, 2147483813)
[ Info: Reconstructing modulo 452312984225514596069828623335196996110407213751397785678981819797976892333
[ Info: Heuristic check failed.
[ Info: 9: selected lucky prime 2147483857
[ Info: CRT modulo (452312984225514596069828623335196996110407213751397785678981819797976892333, 2147483857)
[ Info: 10: selected lucky prime 2147483867
[ Info: CRT modulo (971334831935788242577832613368869049061991281227475155931159242212638377544144568381, 2147483867)
[ Info: 11: selected lucky prime 2147483869
[ Info: CRT modulo (2085925881037261630864178029036094802896257759330662854505473835259586299281105540913875809327, 2147483869)
[ Info: 12: selected lucky prime 2147483887
[ Info: CRT modulo (4479492181457132340213454847299227207974448018628682717127999033421525005319580445599067818799052246163, 2147483887)
[ Info: Reconstructing modulo 9619637281621671881834996425177135896677125027724171971067756840824299427891388292524278203091200389506200075581
[ Info: 13: selected lucky prime 2147483893
[ Info: CRT modulo (9619637281621671881834996425177135896677125027724171971067756840824299427891388292524278203091200389506200075581, 2147483893)
[ Info: 14: selected lucky prime 2147483929
[ Info: CRT modulo (20658016118784845285971654106780479009986238218584837754630069827310747864405871311604659752589335646499890885945580116833, 2147483929)
[ Info: 15: selected lucky prime 2147483951
[ Info: CRT modulo (44362757620113410260375536343857928604867277085576528201140520294297636327782679774913158020198714437645340777821705269480809876857, 2147483951)
[ Info: 16: selected lucky prime 2147483993
[ Info: CRT modulo (95268310011296503334035195531452119103056298026345647894248167227793970831147880232398299267103673085675359550297968805662169313032693822007, 2147483993)
[ Info: 17: selected lucky prime 2147483999
[ Info: CRT modulo (204587170789420890086711714502418553819742917389414771138112095891274737061798978614956467676528769422992752228784566390572836365393516268430023633951, 2147483999)
[ Info: 18: selected lucky prime 2147484007
[ Info: CRT modulo (439347675670961559937551129419800091128618245387347072993342744994866141553145580730162196416906240199037398104286443541908350538727893524799664605101605650049, 2147484007)
[ Info: 19: selected lucky prime 2147484037
[ Info: CRT modulo (943492107016012944277662969173807884835850261977729359411465162345954336091178275160750699321298855245933309223847255653156617041667985469306137618419668743501066266343, 2147484037)
[ Info: 20: selected lucky prime 2147484041
[ Info: CRT modulo (2026134238852383501221651721966775511189722802919441848842357150919550408286738322428905735729076025747015490724796841241411843757904162649282884001681435793496387299450782866691, 2147484041)
[ Info: Reconstructing modulo 4351090942859175723685201076583819542509546642483289579016496306421962976690804699359187424571554264967420869711284739533142562778484656896803273488065100532499663316725644271175152978331
[ Info: Heuristic check passed!
[ Info: Randomized check passed!
[ Info: Success!
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
[ Info: Computing in lex order, result is in lex order
[ Info: Using fglm: false
┌ Info: Computed in 28.119437063 seconds
│   :ioeq_time = :ioeq_time
└   ioeq_time = 28.119437063
[ Info: Computing Wronskians
┌ Info: Computed in 23.338631096 seconds
│   :wrnsk_time = :wrnsk_time
└   wrnsk_time = 23.338631096
[ Info: Dimensions of the wronskians [1282, 2]
┌ Info: Ranks of the Wronskians computed in 0.486833654 seconds
│   :rank_time = :rank_time
└   rank_times = 0.486833654
[ Info: Assessing global identifiability using the coefficients of the io-equations
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Selected tablesize 65536
[ Info: 1: selected lucky prime 2147483647
[ Info: CRT modulo (1, 2147483647)
[ Info: Reconstructing modulo 2147483647
[ Info: Heuristic check failed.
[ Info: 2: selected lucky prime 2147483659
[ Info: CRT modulo (2147483647, 2147483659)
[ Info: Reconstructing modulo 4611686039902224373
[ Info: 3: selected lucky prime 2147483693
[ Info: CRT modulo (4611686039902224373, 2147483693)
[ Info: Reconstructing modulo 9903520567925774155444649489
[ Info: 4: selected lucky prime 2147483713
[ Info: CRT modulo (9903520567925774155444649489, 2147483713)
[ Info: Reconstructing modulo 21267649120981110191733715050621272657
[ Info: 5: selected lucky prime 2147483743
[ Info: CRT modulo (21267649120981110191733715050621272657, 2147483743)
[ Info: Reconstructing modulo 45671930739135174346839766056203605080877915151
[ Info: 6: selected lucky prime 2147483777
[ Info: CRT modulo (45671930739135174346839766056203605080877915151, 2147483777)
[ Info: Reconstructing modulo 98079730326560405919904968824172512120100095704355005327
[ Info: 7: selected lucky prime 2147483783
[ Info: CRT modulo (98079730326560405919904968824172512120100095704355005327, 2147483783)
[ Info: 8: selected lucky prime 2147483813
[ Info: CRT modulo (210624630317301765882893117451031048172285903861850336414611112041, 2147483813)
[ Info: Reconstructing modulo 452312984225514596069828623335196996110407213751397785678981819797976892333
[ Info: 9: selected lucky prime 2147483857
[ Info: CRT modulo (452312984225514596069828623335196996110407213751397785678981819797976892333, 2147483857)
[ Info: 10: selected lucky prime 2147483867
[ Info: CRT modulo (971334831935788242577832613368869049061991281227475155931159242212638377544144568381, 2147483867)
[ Info: 11: selected lucky prime 2147483869
[ Info: CRT modulo (2085925881037261630864178029036094802896257759330662854505473835259586299281105540913875809327, 2147483869)
[ Info: 12: selected lucky prime 2147483887
[ Info: CRT modulo (4479492181457132340213454847299227207974448018628682717127999033421525005319580445599067818799052246163, 2147483887)
[ Info: Reconstructing modulo 9619637281621671881834996425177135896677125027724171971067756840824299427891388292524278203091200389506200075581
[ Info: 13: selected lucky prime 2147483893
[ Info: CRT modulo (9619637281621671881834996425177135896677125027724171971067756840824299427891388292524278203091200389506200075581, 2147483893)
[ Info: 14: selected lucky prime 2147483929
[ Info: CRT modulo (20658016118784845285971654106780479009986238218584837754630069827310747864405871311604659752589335646499890885945580116833, 2147483929)
[ Info: 15: selected lucky prime 2147483951
[ Info: CRT modulo (44362757620113410260375536343857928604867277085576528201140520294297636327782679774913158020198714437645340777821705269480809876857, 2147483951)
[ Info: 16: selected lucky prime 2147483993
[ Info: CRT modulo (95268310011296503334035195531452119103056298026345647894248167227793970831147880232398299267103673085675359550297968805662169313032693822007, 2147483993)
[ Info: 17: selected lucky prime 2147483999
[ Info: CRT modulo (204587170789420890086711714502418553819742917389414771138112095891274737061798978614956467676528769422992752228784566390572836365393516268430023633951, 2147483999)
[ Info: 18: selected lucky prime 2147484007
[ Info: CRT modulo (439347675670961559937551129419800091128618245387347072993342744994866141553145580730162196416906240199037398104286443541908350538727893524799664605101605650049, 2147484007)
[ Info: 19: selected lucky prime 2147484037
[ Info: CRT modulo (943492107016012944277662969173807884835850261977729359411465162345954336091178275160750699321298855245933309223847255653156617041667985469306137618419668743501066266343, 2147484037)
[ Info: 20: selected lucky prime 2147484041
[ Info: CRT modulo (2026134238852383501221651721966775511189722802919441848842357150919550408286738322428905735729076025747015490724796841241411843757904162649282884001681435793496387299450782866691, 2147484041)
[ Info: Reconstructing modulo 4351090942859175723685201076583819542509546642483289579016496306421962976690804699359187424571554264967420869711284739533142562778484656896803273488065100532499663316725644271175152978331
[ Info: 21: selected lucky prime 2147484043
[ Info: CRT modulo (4351090942859175723685201076583819542509546642483289579016496306421962976690804699359187424571554264967420869711284739533142562778484656896803273488065100532499663316725644271175152978331, 2147484043)
[ Info: 22: selected lucky prime 2147484061
[ Info: CRT modulo (9343898369431904662746946467210173419530811589897090265086113451809603917180284036703267319713678895726130232570165995176834923210921544406214927525784754338733829875540776081243005879049747272233, 2147484061)
[ Info: 23: selected lucky prime 2147484083
[ Info: CRT modulo (20065872815958904888120648014754106555488283987698069974550693430198816018867824032273005555707316512243945695654654538766455874023612957733850266201892926459231494579209430390039396192988626093198595378213, 2147484083)
[ Info: 24: selected lucky prime 2147484109
[ Info: CRT modulo (43091142483774136629349987395329992986797046156565373080167829217964628925964079790151157741452032016686947994534732886864670443787561992885995196973877924041489183021153034486106085067373870735182458132669782483679, 2147484109)
[ Info: 25: selected lucky prime 2147484161
[ Info: CRT modulo (92537543722559748756723920930821460750228103429363664709316796298604937942589598854157665957720569341594443645972757723101574611555767151575045734151727729985007077233318752416041798970387581765624476055466171628187204357011, 2147484161)
[ Info: 26: selected lucky prime 2147484167
[ Info: CRT modulo (198722909442042038831204062448755463679998029251628252272172489182517530628099090862647336640933818324976266215333588647851055272475737526211496916941452070927287465831455722277730245602853439872770976603037361042839602499551311802771, 2147484167)
[ Info: 27: selected lucky prime 2147484197
[ Info: CRT modulo (426754301646960082538409909654781607107539322409074530724372195212435170223780362934629527141124338947740992348205894244551079771883517229186937122501082388305710841130604654152974879129149132148262166672150176948960655088360066700530949226757, 2147484197)
[ Info: 28: selected lucky prime 2147484203
[ Info: CRT modulo (916448118788617850341050926491841226749703574429575524125780251964903585942573283001041453585147106702345389916870079192426697169406218174456174629403727544281531636179551107848163973467832883344657663941422585008646682377020881885256144973870027059129, 2147484203)
[ Info: 29: selected lucky prime 2147484221
[ Info: CRT modulo (1968057857967624329811225027059743242829129461020148174055558476143990161229728990414584954122261127074442147896353978269095329406764668419616133132453884210660572173340329275735081455575882245626594137757087718653493368813010834049716429940063730884662074439187, 2147484221)
[ Info: 30: selected lucky prime 2147484223
[ Info: CRT modulo (4226373196000532377125305654291096638286926916707002766866273404862023795220088963021581437241564675234040404984768511763459111845824416091422652778980019352554618729080033982585445601998919590637157168804446206721264376054074265623815562820738837809202175775281506568327, 2147484223)
[ Info: 31: selected lucky prime 2147484233
[ Info: CRT modulo (9076069758921229979477279986642822280089513300782323555462669287745687592085723860745276544985924779899220602249320934319218350594500341984517472467666697591845718466479684281906101179717417083938913537578896001186111806188063480297455174219402031398617556694689828739153104004921, 2147484233)
[ Info: 32: selected lucky prime 2147484239
[ Info: CRT modulo (19490716704891452469894372353041911449113339642073826400460583315527204217747827575342369009581988671757571572319181041407349996785955660924859202237325835377667749775322061010211279469945852583243654357100932156093924402364018156741891136660433945116702140998830002042801560623577001910593, 2147484239)
[ Info: 33: selected lucky prime 2147484259
[ Info: CRT modulo (41856006930568408384915786622954848543404547406007443469411205010841035033347783794537322476999360650875929380469889963809890496854540438388963300098770768981050105221099917168425140721732992605933183228637929537419990457734223332313042807032076992038708863452543146827254154843734623405871354643727, 2147484259)
[ Info: 34: selected lucky prime 2147484271
[ Info: CRT modulo (89885116027990562929290264813398305314690343823420266887392910759003047085381905771304190207363016550820052906554710720743819510508814604119258006290803371636130570252605786785396841519781500822204900989262751691960580979914449411732785488494559848479195602948116801329854089719268708537401702276410285593293, 2147484271)
[ Info: 35: selected lucky prime 2147484277
[ Info: CRT modulo (193026872867119729627086528880197610721353218596377024563298404051865715256930036671879871626704306429998735768214074073752425819280598569201198376700319292542357934919731423885219369656810508372458592413554099144663984806388047037321359691559118743677215828063442059926693540397152357206553570847215982654214620594403, 2147484277)
[ Info: 36: selected lucky prime 2147484331
[ Info: CRT modulo (414522174520617529650659394088730785677072665099063669413726113960574716029616269041895432351125911386612286192117240443496672837513928878508269963521858821614477021746312489856330849033851452698231687041658909602065055819885220173384052133884777118022814153901777182193066150600349022674561434752601891973030625510000836701631, 2147484331)
[ Info: Reconstructing modulo 890179874635073581368718952623503385918832774248649792837340786055854554428375469710150823514513384910844367769659429567368595769194471260844912400581133395411213526959752328996066939452122483251040219329658252616980152615843868540827355202804673020381331096029089011812941775260735769324784393428091424493037942865855680613642294643861
[ Info: Heuristic check passed!
[ Info: Randomized check passed!
[ Info: Success!
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
[ Info: Computing in degrevlex order, result is in degrevlex order
[ Info: Using fglm: false
┌ Info: Computed in 38.932922401 seconds
│   :check_time = :check_time
└   check_time = 38.932922401
[ Info: Global identifiability assessed in 90.879020625 seconds
102.768808 seconds (206.16 M allocations: 19.365 GiB, 14.02% gc time, 12.02% compilation time)
Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Symbol} with 8 entries:
r     => :nonidentifiable
N_inv => :globally
g1    => :globally
alpha => :globally
q     => :nonidentifiable
g2    => :globally
k     => :globally
b     => :globally