DomainError while solving ODE

Specific Domains Modelling & Simulations
1

I have been trying to estimate the parameters of the model included in the code below. When I run the code below, it solves for random parameters it generates. However, sometimes I get this: ERROR: DomainError with -0.11794106990249296: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

From what I understand G_b/G̃ might yield a negative number during integration, the problem might be caused by that. On the other hand, optimization run is interrupted by this error while trying to estimate the parameters. How can I avoid this error while trying to estimate the parameters?

function simo(du,u,p,t)
    S, J, R, L, G, I, G_prod = u
    k_js, k_gj, k_rj, k_lr, k_gl,
    k_xg, k_xgi, η, I_b, G_b, f_gi, β, γ, k_xi, k_λ = p
    du[1] = dS = -k_js*S
    du[2] = dJ = k_js*S - k_gj*J - k_rj*J
    du[3] = dR = -k_lr*R + k_rj*J
    du[4] = dL = k_lr*R - k_gl*L
    G0_prod = (k_xg + k_xgi*I_b)*G_b
    G_prod = k_λ/((k_λ/G0_prod) + (G - G_b))
    du[5] = dG = -k_xg*G - k_xgi*I*G + G_prod + η*(k_gj*J + k_gl*L)
    G̃ = G + f_gi*(k_gj*J + k_gl*L)
    ins_rhs = (β.^γ + 1)/((β.^γ)*((G_b/G̃).^γ) + 1) - I/I_b
    du[6] = dI = k_xi*I_b*ins_rhs
end

time_points = [0.0,30.0,60.0,90.0,120.0]
tspan = (0.0,120.0)

u0 = [
    416.31; 
    0.; 
    0.; 
    0.;
    G0;
    I0;
    0.01;
]

p0 = rand(15)
p0[9] = I0
p0[10] = G0
prob = ODEProblem(simo, u0, tspan, p0)
sol = solve(prob,Tsit5(),saveat=time_points)
plot(sol)
2

You had a value go negative. Use a lower tolerance or isoutofdomain.

3

Thank you for your reply @ChrisRackauckas! I decreased both absolute and relative tolerance to 1e-2 but it did not seem to solve the problem. I will try isoutofdomain.

4

If it doesn’t solve it, check the model. Is the model correct? Likely the model has a true solution which goes negative.

5

Thanks again @ChrisRackauckas, I set the isoutofdomain assuming that the negative number is G̃ which is defined as G̃ = G + f_gi*(k_gj*J + k_gl*L) (for calculation of (G_b/G̃).^γ), however it does not seem to help. I tried the following:

sol = solve(prob,Tsit5(),saveat=time_points,isoutofdomain=(u,p,t)->u[5] + p[11]*(p[2]*u[2] + p[5]*u[4]) < 0)

Both β and γ are positive numbers thus I do not think they cause the problem (for the term β.^γ).

6

For that to work, you need to make it not fail in your rhs, i.e. abs(...)^gamma, and then isoutofdomain will cause a rejection and not accept a negative step. But if it errors before it rejects it will halt.

7

Alright, I will check my model and try to understand at what part the variables go negative. Using abs did not seem to help as well.

It would have been simpler if I had more ideas about parameter values (and get one of the solutions) but currently the best I can do is to search for the parameters by fitting to data. I will report back when I find out the issue.

8

@ChrisRackauckas, you mentioned earlier that there is a solution which goes negative and that seems to be the case. However, since this is a biological system, none of the variables should be allowed to go negative. To avoid it, I tried PositiveDomain as such:

sol = solve(prob,Tsit5(),saveat=time_points,cb=PositiveDomain(zeros(7), abstol=1.0e-8))

This did not help. Instead I added the following lines inside the model function (following the DiffEqCallbacks docs):

    u⁺ = float.(u .> 0.0)
    du = u⁺.*du

I expected this one to work but it did not, then I removed this part as well.

Eventually, I added u = max.(0.0,u) in the beginning of the model function and that seems to do the trick. However, I am not sure how correct this approach is.

9

That’s not the point. You need to allow it to go negative in order to reject the negative solution.

If you have this with PositiveDomain or isoutofdomain, then it’s correct. Those will force the solution to be positive, but to act they have to allow the solution to go slightly negative to trigger. If you error during that triggering, then it can’t work.

10

Hi @zamk @ChrisRackauckas and sorry for reviving this old thread.

When negative solutions are forbidden by the equations, I understand that isoutofdomain can only enforce that the approximation be positive if dudt does not fail given a negative u as input. Only then can it be triggered and correctly reject negative approximations.

However, I am unsure how to determine the correct way to make the RHS not failing. For example in the above, failure happens when raising u^g with a non-integer value of g, and possible fixes include:

# Nullify negative input before it gets even raised to the power g.
u = max(u, 0.0)
p = u^g

# Raise only the magnitude of u to g.
p = abs(u)^g

# Use the `abs` trick, but allow negative output to let `isoutofdomain` triggered.
p = sign(u) * abs(u)^g

# Refuse to do the calculation and produce dummy output instead.
if u >= 0
    p = u^g
else
    p = -1
end

What difference, if any, does it make to pick either of the above possible infallible versions of dudt? If any, then how to determine the most adequate?

11

Those should all be rather similar.