Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear Interpolation as an algebraic function in a reaction network #647

Closed
danilodurs opened this issue May 11, 2023 · 2 comments
Closed

Comments

@danilodurs
Copy link

Hello ,

I am currently using Catalyst to build up a reaction network from which I generate a system of ODEs .
The model I am implementing implies interpolated data as an algebraic function, as described in Bergman's minimal model (Bergman et al 1979) :

Iintpol(t) = linear_interpolation(tticks,datamean2)(t)
minimalmodel = @reaction_network minimalmodel begin
    B_0,     ∅   ⇒   G_p            # hepatic glucose production
    k_gelim, G_p --> ∅              # glucose utilization/clearance
    (k_iut + k_imet)*I_r, G_p --> ∅ # insulin-dependent glucose regulation
    k_gmet, G_p --> ∅               # glucose hepatic metabolization
    k_ipan*Iintpol(t), ∅ ⇒ I_r        # blood insulin --> remote insulin            
    k_ielim, I_r --> ∅              # insulin clearance
end B_0 k_gelim k_iut k_imet k_gmet k_ipan k_ielim

However, unlike other algebraic functions I have previously used, some implying the variable t, in different models, linear interpolation does not seem to work inside an ODE:

TypeError: non-boolean (Num) used in boolean context

I also have kind of the same problem when I directly implement the system of ODEs:

function fminimalmodel(u, p, t)
    B_0,k_gelim,k_iut,k_imet,k_gmet, k_ipan ,k_ielim = p
    G_p , I_r = u
    dGp = B_0 - k_gelim * G_p - k_gmet * G_p - (k_iut + k_imet) * G_p * I_r
    dIr = k_ipan * Iintpol(t) - k_ielim* I_r
    [dGp, dIr]
end
prm = [datamean1[1],0.03,0.04,0.005,0.09,0.9,0.005]
oprob = ODEProblem(fminimalmodel, [datamean1[1],0.0], tspan, prm)
sol = solve(oprob)

output:

BoundsError: attempt to access 9-element extrapolate(interpolate((::Vector{Union{Missing, Int64}},), ::Vector{Float64}, Gridded(Linear())), Throw()) with element type Float64 at index [305.0204580512005]

Is there any workaround in Catalyst or DifferentialEquations to solve this issue?

@isaacsas
Copy link
Member

Can you post a complete minimum working example that loads the appropriate package and builds the interpolant so we can run the code ourselves? Your DifferentialEquations error seems to be related to the interpolant -- are you sure you are building / using it correctly? Have you checked you can call Iintpol(t) and plot it yourself from how you've defined it?

Regarding Catalyst, you probably need to register your function using @register_symbolic. See https://symbolics.juliasymbolics.org/stable/manual/functions/

@danilodurs
Copy link
Author

For the DifferentialEquations object, it was just a misuse indeed (tspan[end] > tticks[end]): issue solved.
The register_symbolic macro did the trick, the model works as expected, I've never used it before, good to know.
Thank you

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants