Stuck at "Computing IO-equations"

[LinneaFranssen]

Hi,

I am new to both structural identifiability analysis and Julia. So I hope the answer to my question is not completely obvious. I am able to obtain the local SI results but when commuting the global ones, the output in the terminal always stalls at "Computing IO-equations". I have left it running on an HPC over night.

Here's the code:

#import Pkg; 
#using Pkg
#Pkg.add("StructuralIdentifiability")
#Pkg.add("OrdinaryDiffEq")

using StructuralIdentifiability, OrdinaryDiffEq

odeMyModel = @ODEmodel(
    Complex'(t) = 1/C*(((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) - (ke_Complex*Complex(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)),
    Complex_2'(t) = 1/C*(((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C) - (ke_Complex_2*Complex_2(t)*C)),
        Drug'(t) = 1/C*(-(ke_Drug*Drug(t)*C) - ((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (45/100*ka*Drug_SC(t))),
        free_receptor'(t) = 1/C*(-((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (66/2500*C) - (kdeg_free_receptor*free_receptor(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)),
    Drug_SC'(t) = -(45/100*ka*Drug_SC(t)) - ((1-45/100)*ka*Drug_SC(t)) + u_SC(t),
        y1(t) = Drug(t),                                   # output free Drug
        y2(t) = free_receptor(t) + Complex(t) + 2*Complex_2(t)      # output total receptor
)

local_id_all = assess_local_identifiability(odeMyModel, 0.99)

assess_identifiability(odeMyModel,0.99) 

The parameters, state variables, inputs and outputs seem to be correctly recognized:

[ Info: Summary of the model:
[ Info: State variables: Drug, free_receptor, Complex, Complex_2, Drug_SC
[ Info: Parameters: ke_Complex_2, kdeg_free_receptor, ka, kon, ke_Drug, ke_Complex, kon_2, koff, C
[ Info: Inputs: u_SC
[ Info: Outputs: y1, y2
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.091845932 seconds
[ Info: Assessing global identifiability
[ Info: Computing IO-equations

But it doesn't move past there.

Also, I was looking to include initial conditions but have failed in the attempt in multiple ways. I think it might not be possible yet. If that is the case, can you recommend an alternative for SI - maybe even additionally allowing bounds for parameters?

Thanks a lot in advance for your help!

[pogudingleb]

@LinneaFranssen , thanks for the example! Global identifiability is a much more challenging computation, so it may be stuck. Your system does not look too hard, I will check what causes the problem.

For the initial conditions, this tool can only include them into the local identifiability analysis. I suggest you to try this webapp (the same algorithm has also implementations in maple and julia)

[LinneaFranssen]

Dear @pogudingleb,

thank you very much for the quick response and for looking into this!

I had tried that app actually but it kept disconnecting quickly. I will try again. :)

[LinneaFranssen]

Just an update: it still happens with this code

odeMyModel = @ODEmodel(
    Complex'(t) = 1/C*(((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) - (ke_Complex*Complex(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)),
    Complex_2'(t) = 1/C*(((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C) - (ke_Complex_2*Complex_2(t)*C)),
    Drug'(t) = 1/C*(-(ke_Drug*Drug(t)*C) - ((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (Fa*ka*Drug_SC(t))),
    free_receptor'(t) = 1/C*(-((2*kon*free_receptor(t)*Drug(t)-koff*Complex(t))*C) + (66/2500*C) - (kdeg_free_receptor*free_receptor(t)*C) - ((kon_2*free_receptor(t)*Complex(t)-2*koff*Complex_2(t))*C)) + u_free_receptor(t),
    Drug_SC'(t) = -(Fa*ka*Drug_SC(t)) - ((1-Fa)*ka*Drug_SC(t)) + u_drug_SC(t),
    y1(t) = Drug(t),                                   # output free Drug
    y2(t) = free_receptor(t) + Complex(t) + 2*Complex_2(t)      # output total receptor
)

 SI_output = assess_identifiability(odeMyModel,0.99)

and Julia crashes out silently without error messages - even from the REPL.

[pogudingleb]

@LinneaFranssen Just two clarify, is the second system the same as the first one?

[LinneaFranssen]

No, I slightly modified it since, e.g. the input. I'm interested in the second one but both cause the same problem. Thanks!

On Thu, 15 Dec 2022, 10:31 pogudingleb, @.***> wrote:

@LinneaFranssen https://github.com/LinneaFranssen Just two clarify, is the second system the same as the first one?

— Reply to this email directly, view it on GitHub https://github.com/SciML/StructuralIdentifiability.jl/issues/132#issuecomment-1352781970, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHADLVC2YTQVVYTZD3F7HJ3WNLQORANCNFSM6AAAAAAS3F5BJ4 . You are receiving this because you were mentioned.Message ID: @.***>

[pogudingleb]

Very interesting case! The model is in some way special, and I did not optimise the code for models of this sort thinking that this does not happen in practice. I have an idea how to speed up the computation, I will try it over the weekend. I also did the computation a bit differently (by changing manually some things in the code) and I think that all the parameters in your original model are globally identifiable.

[pogudingleb]

Sorry, taking longer to fix this. Hope to do this by the end of the year)

[pogudingleb]

@LinneaFranssen Sorry, this has not been fixed last year as hoped, but I did not forget and hope to get to this asap

[pogudingleb]

@LinneaFranssen It took a bit longer than I anticipated, but it works now! The version on GitHub (can be installed by Pkg.add("https://github.com/SciML/StructuralIdentifiability.jl"); it will be in the release 0.4.9) computes both models in a couple of minutes on my laptop. Let me know if you encounter any further issues.

The model turned out to be very interesting algorithmically, I wonder if it is published somewhere?

[filipepfarias]

Same thing here. But for me I think it might be the dynamical system itself. I will let the system[1] here in case of some insight.

asm1_ode = @ODEmodel(SI'(t) = Q(t) * (QIN_SI(t) - SI(t)) / V - QEC * SI(t) / V,
       SS'(t) = Q(t) * (QIN_SS(t) - SS(t)) / V + QEC * (SEC - SS(t)) / V - μH * SS(t) / YH / (KS + SS(t)) * (SO(t) / (KOH + SO(t)) + νg * KOH * SNO(t) / (KOH + SO(t)) / (KNO + SNO(t))) * XBH(t) + kh * XS(t) / (KX * XBH(t) + XS(t)) * (SO(t) / (KOH + SO(t)) + νh * KOH * SNO(t) / (KOH + SO(t)) / (KNO + SNO(t))) * XBH(t),
       XI'(t) = Q(t) / V * (QIN_XI(t) - XI(t)) - QEC * XI(t) / V,
       XS'(t) = Q(t) / V * (QIN_XS(t) - XS(t)) - QEC * XS(t) / V - kh * XS(t) / (KX * XBH(t) + XS(t)) * (SO(t) / (KOH + SO(t)) + νh * KOH * SNO(t) / (KOH + SO(t)) / (KNO + SNO(t))) * XBH(t) + (1 - fP) * bH * XBH(t) + (1 - fP) * bA * XBA(t),
       XBH'(t) = Q(t) * (QIN_XBH(t) - XBH(t)) / V - QEC * XBH(t) / V + μH * SS(t) / (KS + SS(t)) * (SO(t) / (KOH + SO(t)) + νg * KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t))) * XBH(t) - bH * XBH(t),
       XBA'(t) = Q(t) * (QIN_XBA(t) - XBA(t)) / V - QEC * XBA(t) / V + μA * SNH(t) / (KNH + SNH(t)) * SO(t) / (KOA + SO(t)) * XBA(t) - bA * XBA(t),
       XP'(t) = Q(t) * (QIN_XP(t) - XP(t)) / V - QEC * XP(t) / V + fP * (bH * XBH(t) + bA * XBA(t)),
       SO'(t) = Q(t) * (QIN_SO(t) - SO(t)) / V - QEC * SO(t) / V + KLa4 * (Ssat - SO(t)) - (1 - YH) / YH * μH * SS(t) / (KS + SS(t)) * SO(t) / (KOH + SO(t)) * XBH(t) - (4.57 - YA) / YA * μA * SNH(t) / (KNH + SNH(t)) * SO(t) / (KOA + SO(t)) * XBA(t),
       SNO'(t) = Q(t) * (QIN_SNO(t) - SNO(t)) / V - QEC * SNO(t) / V - QEC * SNO(t) / V - (1 - YH) / (2.86 * YH) * μH * SS(t) / (KS + SS(t)) * KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t)) * νg * XBH(t) + μA / YA * SNH(t) / (KNH + SNH(t)) * SO(t) / (KOA + SO(t)) * XBA(t),
       SNH'(t) = Q(t) * (QIN_SNH(t) - SNH(t)) / V - QEC * SNH(t) / V - μH * SS(t) / (KS + SS(t)) * (SO(t) / (KOH + SO(t)) + νg * KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t))) * iXB * XBH(t) - μA * (iXB + 1 / YA) * SNH(t) / (KNH + SNH(t)) * SO(t) / (KOA + SO(t)) * XBA(t) + ka * SND(t) * XBH(t),
       SND'(t) = Q(t) * (QIN_SND(t) - SND(t)) / V - QEC * SND(t) / V - ka * SND(t) * XBH(t) + kh * XND(t) * XBH(t) / (KX * XBH(t) + XS(t)) * (SO(t) / (KOH + SO(t)) + νh * KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t))),
       XND'(t) = Q(t) * (QIN_XND(t) - XND(t)) / V - QEC * XND(t) / V - kh * XND(t) / (KX * XBH(t) + XS(t)) * (SO(t) / (KOH + SO(t)) + νh + KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t))) * XBH(t) + bH * (iXB - fP * iXP) * XBH(t) + bA * (iXB - fP * iXP) * XBA(t),
       SALK'(t) = Q(t) * (QIN_SALK(t) - SALK(t)) / V - QEC * SALK(t) / V - iXB / 14 * μH * SS(t) / (KS + SS(t)) * SO(t) / (KOH + SO(t)) * XBH(t) + 1 / 14 * ka * SND(t) * XBH(t) + ((1 - YH) / (14 * 2.86 * YH) - iXB / 14) * μH * SS(t) / (KS + SS(t)) * KOH / (KOH + SO(t)) * SNO(t) / (KNO + SNO(t)) * νg * XBH(t) - (iXB / 14 + 1 / (7 * YA)) * μA * SNH(t) / (KNH + SNH(t)) * SO(t) / (KOA + SO(t)) * XBA(t),
       y1(t) = SO(t),
       y2(t) = SNO(t),
       y3(t) = SNH(t),
       y4(t) = SALK(t)
       )

[1] Henze, Mogens, et al. Activated sludge models ASM1, ASM2, ASM2d and ASM3. IWA publishing, 2006.

[pogudingleb]

@filipepfarias Indeed, in this case this could be more like the size of the model, but I think we may manage to do something with it still. Just to confirm: QIN_SND, QIN_SS, QIN_XS, QIN_SNO, QIN_SNH, QIN_XND, QIN_XBH, QIN_XBA, QIN_XP, Q, QIN_XI, QIN_SI, QIN_SALK, QIN_SO all are external inputs, am I reading this correctly?

[filipepfarias]

@pogudingleb That’s correct.