Trouble Getting Newton solver to converge for ContinuationPar

New to Julia
,
1

I’m trying to use Julia to analyze the Luo-Rudy model of a cardiac action potential (similar to what Tran et al did at PMID 19659123). I have the model producing reasonable-looking trajectories and I would like to make a bifurcation diagram (I’m expecting to see a Hopf bifurcation as I vary the background current-- when I manually vary it, I do see the system transition from a normal action potential to one with early after-depolarizations). However, I’m having trouble getting the Newton solver to converge in ContinuationPar.

Here is the model:

using DynamicalSystems
using DifferentialEquations
using CairoMakie
using Revise, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit

# CONSTANTS
C_m = 1  #muF/cm2
Esi = 80 #mV per Tran-- in LR had been dependent on Ca_i
Ek = -77 #mV
EK1 = -84  
Ko = 5.4 #mM
ENa = 54.4 #per LR
G_K = 0.282*sqrt(Ko/5.4) #1 per Life/Death SN --  in LR  #mS/cm2
G_Si = 0.07
G_K1 = 0.6047*sqrt(Ko/5.4) #Life/Death added 0.4 factor 

# TIME DEPENDENT POTASSIUM CURRENT
function alpha_w(V)
    return 0.0005*exp(0.083*(V+50))/(1 + exp(0.057*(V+50)))
end

function beta_w(V)
    return 0.0013*exp(-0.06*(V+20))/(1 + exp(-0.04*(V+20)))
end

function xi(V)
    if V <= -100
        return 1
    else
        return 2.837*(exp(0.04*(V+77))-1)/((V+77)*exp(0.04*(V+35)))
    end
end

function I_K(V, w, xi)
    return G_K*w*xi*(V-Ek)
end

# TIME INDEPENDENT POTASSIUM CURRENT
function alpha_K1(V)
    return 1.02/(1 + exp(0.2385*(V-EK1-59.215)))
end

function beta_K1(V)
    return (0.49124*exp(0.08032*(V-EK1+5.476))+exp(0.06175*(V-EK1-594.31)))/(1 + exp(-0.5143*(V-EK1+4.753)))
end

function I_K1(V, K1x)
    return G_K1*K1x*(V-EK1)
end

# PLATEAU POTASSIUM CURRENT
function Kp(V)
    return 1/(1 + exp((7.488-V)/5.98))
end

function I_Kp(V)
    return 0.0183*Kp(V)*(V-EK1) #EK1 
end

# BACKGROUND CURRENT
function I_b(V, Ib_param) 
    return 0.03921*(V + 59.87)*Ib_param # + Ib_param
end

# CALCIUM CURRENT (SLOW INWARD CURRENT)
function alpha_d(V)
    return 0.095*exp(-0.01*(V-5))/(1 + exp(-0.072*(V-5)))
end

function beta_d(V)
    return 0.07*exp(-0.017*(V+44))/(1 + exp(0.05*(V+44)))
end

function alpha_f(V)
    return 0.012*exp(-0.008*(V+28))/(1 + exp(0.15*(V+28)))
end

function beta_f(V)
    return 0.0065*exp(-0.02*(V+30))/(1 + exp(-0.2*(V+30)))
end

function I_si(d, f, V, gsi_param)  # in Tran paper they have a GSI param instead of 0.09
    return 0.09*gsi_param*d*f*(V-Esi)  #in LR, ESI was dependent on Ca_i, Qu fixes it to a constant
end


# SODIUM CURRENT 
function alpha_h(V)
    if V >= -40
        return 0
    else
        return 0.135*exp((80+V)/-6.8)
    end
end

function alpha_j(V)
    if V >= -40
        return 0
    else
        return (-1.2714*10^5*exp(0.2444*V)-3.474*10^-5*exp(-0.04391*V))*(V+37.78)/(1+exp(0.311*(V+79.23)))
    end
end

function beta_h(V)
    if V >= -40
        return 1/(0.13*(1+exp((V+10.66)/-11.1)))
    else
        return 3.56*exp(0.079*V)+3.1*10^5*exp(0.35*V)
    end
end

function beta_j(V)
    if V >= -40
        return 0.3*exp(-2.535*10^-7*V)/(1+exp(-0.1*(V+32)))
    else
        return 0.1212*exp(-0.01052*V)/(1+exp(-0.1378*(V+40.14)))
    end
end

function alpha_m(V)
    return 0.32*(V+47.13)/(1-exp(-0.1*(V+47.13)))
end

function beta_m(V)
    return 0.08*exp(-V/11)
end

function I_Na(m, h, j, V)
    return 23*m^3*h*j*(V-ENa)
end
 
# STIMULUS 
function I_stim(t)
    if t%3000 < 2
        return -30
    else
        return 0
    end
end

# DIFFERENTIAL EQUATIONS
function d_dot(d, V, tau_d_param)
    return (alpha_d(V)*(1-d) - beta_d(V)*d)/tau_d_param
end

function f_dot(f, V, tau_f_param)
    return (alpha_f(V)*(1-f) - beta_f(V)*f)/tau_f_param
end

function w_dot(w, V, tau_w_param)
    return (alpha_w(V)*(1-w) - beta_w(V)*w)/tau_w_param
end

function K1x_dot(K1x, V)
    return alpha_K1(V)*(1-K1x) - beta_K1(V)*K1x
end

function h_dot(h, V)
    return alpha_h(V)*(1-h) - beta_h(V)*h
end

function j_dot(j, V)
    return alpha_j(V)*(1-j) - beta_j(V)*j
end

function m_dot(m, V)
    return alpha_m(V)*(1-m) - beta_m(V)*m
end

function V_dot(d, f, w, K1x, m, h, j, V, gsi_param, Ib_param, t)
    return -1/C_m*(I_si(d,f,V,gsi_param)+ I_K(V, w, xi(V)) + I_b(V, Ib_param) + I_K1(V, K1x) + I_Kp(V) + I_Na(m, h, j, V) + I_stim(t))
end

function LR_EOM!(du, u, p, t)
    @unpack gsi_param, tau_d_param, tau_f_param, tau_w_param, Ib_param = p
    d, f, w, K1x, m, h, j, V = u
    du[1] = d_dot(d, V, tau_d_param)
    du[2] = f_dot(f, V, tau_f_param)
    du[3] = w_dot(w, V, tau_w_param)
    du[4] = K1x_dot(K1x, V)
    du[5] = m_dot(m, V)
    du[6] = h_dot(h, V)
    du[7] = j_dot(j, V)
    du[8] = V_dot(d, f, w, K1x, m, h, j, V, gsi_param, Ib_param, t)
    du
end

The code above is working as I expect. Here’s what I tried to do for the bifurcation diagram (modeling it after Neural mass equation (Hopf aBS) · Bifurcation Analysis in Julia)

#OOP function for ContinuationPar
function LR_EOM(u, p, t = 21000)
    @unpack gsi_param, tau_d_param, tau_f_param, tau_w_param, Ib_param = p
    d, f, w, K1x, m, h, j, V = u
    ddt = d_dot(d, V, tau_d_param)
    dft = f_dot(f, V, tau_f_param)
    dwt = w_dot(w, V, tau_w_param)
    dK1t = K1x_dot(K1x, V)
    dmt = m_dot(m, V)
    dht = h_dot(h, V)
    djt = j_dot(j, V)
    dVt = V_dot(d, f, w, K1x, m, h, j, V, gsi_param, Ib_param, t)
    return SVector(ddt, dft, dwt, dK1t, dmt, dht, djt, dVt, gsi_param, tau_d_param, tau_f_param, tau_w_param, Ib_param)
end

# BIFURCATION DIAGRAM
u0 = [0.7, 0.4, 0.1, 0.1, 0.9, 0.1, 0.1, -15]  # [d, f, w, K1x, m, h, j, V]

param_LR = (gsi_param = 1.1, tau_d_param = 1.0, tau_f_param = 1.0, tau_w_param = 8.0, Ib_param = 1.0)

prob = BifurcationProblem(LR_EOM, u0, param_LR, (@lens _.Ib_param);
    record_from_solution = (x,p) -> (V = x[8]))

opts_br = ContinuationPar(p_min = 1.0, p_max = 1.2, max_steps = 10000,
    newton_options = NewtonPar(tol = 1e-6, max_iterations = 100))

br = BifurcationKit.continuation(prob, PALC(), opts_br; normC = norminf)

Based on what it says in 4.2.3 of the Datseris textbook, I’ve tried playing around with the dsmin, dsmax, and a parameters in ContinuationPar and the alpha parameter in NewtonPar, but it keeps failing to converge. It also doesn’t converge when I reset the guess solution to the resting potential of the cell where I know there must be a stable equilibrium point (using [0.1, 1, 0.1, 0.5, 0.1, 1., 1., -80] as the initial guess).

If it’s helpful, here is the stack trace for the u0 parameters shown above (8 is the lowest residual I’ve been able to achieve):
Newton failed to converge the initial guess on the branch.
101-element Vector{Float64}:

  • 37.91588404085118*
  • 363.124780002418*
  • 591.4408386814689*
  • 112.12446959498052*
  • 27.55160911874179*
  • 9.061111540232528*
  • 200.92761204072872*
  • 66.09356729757563*
  • 5.796350056800617e7*
  • 3.742416383692493e11*
  • 5.546936220325984e9*
  • 7.267771207898572e7*
  • 1168.6041717904889*
  • ⋮*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
  • 8.0*
    “”

Stacktrace:

  • [1] iterate(it::ContIterable{BifurcationKit.EquilibriumCont, BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#101#102"}, PALC{Secant, MatrixBLS{DefaultLS}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, Float64, DefaultLS, DefaultEig{typeof(real)}, typeof(norminf), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing}; _verbosity::Int64)*
  • @ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:275*
  • [2] iterate*
  • @ ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:250 [inlined]*
  • [3] continuation(it::ContIterable{BifurcationKit.EquilibriumCont, BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#101#102"}, PALC{Secant, MatrixBLS{DefaultLS}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, Float64, DefaultLS, DefaultEig{typeof(real)}, typeof(norminf), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing})*
  • @ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:471*
  • [4] continuation(prob::BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#101#102"}, alg::PALC{Secant, MatrixBLS{Nothing}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, contparams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}; linear_algo::Nothing, bothside::Bool, kwargs::Base.Pairs{Symbol, typeof(norminf), Tuple{Symbol}, NamedTuple{(:normC,), Tuple{typeof(norminf)}}})*
  • @ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:537*
  • [5] top-level scope*
  • @ In[376]:2*

I’m very new to Julia and don’t have much experience with continuation algorithms to begin with, so I apologize if this is a silly question, but I would greatly appreciate any advice about what I’m doing wrong or what else I could try. Thank you!

2

I’m afraid I can’t help with your bifurcation problem, but @datseris is active on this forum. When posting examples like this, it’s important to include all the necessary packages at the top of the code block (using ...) so readers don’t have to guess.

3

Fixed it, thank you for the pointer!

4

but that is not what you pass to BifurcationProblem.

I also recommand to remove all global variables. Put them in param_LR for example and change the vector field accordingly.

5

Sorry if I wasn’t clear, the code above has the initial guess placed where I think there is an unstable fixed point. Here’s the stack trace for what happens when I use the repolarized state (ie, the stable equilibrium) as the initial guess:

Newton failed to converge the initial guess on the branch.
20-element Vector{Float64}:
34.30982691295988
2431.890131667749
860.271324867491
309.44143297740834
113.04055037368715
41.59915715090505
15.557874821967618
8.0
14273.193282071685
498.73964466260077
466.1050862482896
1128.7473464513098
140.33678184913302
20.586298380570163
8.0
343.258288535251
171.86552491078655
29.32914957538595
1.35230523147794e141
NaN
“”

Stacktrace:
[1] iterate(it::ContIterable{BifurcationKit.EquilibriumCont, BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#103#104"}, PALC{Secant, MatrixBLS{DefaultLS}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, Float64, DefaultLS, DefaultEig{typeof(real)}, typeof(norminf), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing}; _verbosity::Int64)
@ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:275
[2] iterate
@ ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:250 [inlined]
[3] continuation(it::ContIterable{BifurcationKit.EquilibriumCont, BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#103#104"}, PALC{Secant, MatrixBLS{DefaultLS}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, Float64, DefaultLS, DefaultEig{typeof(real)}, typeof(norminf), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:471
[4] continuation(prob::BifurcationProblem{BifFunction{typeof(LR_EOM), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, NamedTuple{(:gsi_param, :tau_d_param, :tau_f_param, :tau_w_param, :Ib_param), NTuple{5, Float64}}, Setfield.PropertyLens{:Ib_param}, typeof(BifurcationKit.plot_default), var"#103#104"}, alg::PALC{Secant, MatrixBLS{Nothing}, Float64, BifurcationKit.DotTheta{BifurcationKit.var"#278#280", BifurcationKit.var"#279#281"}}, contparams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}; linear_algo::Nothing, bothside::Bool, kwargs::Base.Pairs{Symbol, typeof(norminf), Tuple{Symbol}, NamedTuple{(:normC,), Tuple{typeof(norminf)}}})
@ BifurcationKit ~/.julia/packages/BifurcationKit/5DxOO/src/Continuation.jl:537
[5] top-level scope
@ In[380]:2

6

Your first residual is 34.30, this is far from an initial guess. If you simulate and make it converge to a stationary point using DifferentialEquations, give the last value of the ODESolution to BifurcationProblem and that should converge.

Also you should use the same function:

LR_EOM(u, p) = LR_EOM!(similar(u),u,p, 21000) and pass it to BifurcationProblem.

7

This works:

# BIFURCATION DIAGRAM
u0 = [0.7, 0.4, 0.1, 0.1, 0.9, 0.1, 0.1, -15]  # [d, f, w, K1x, m, h, j, V]

param_LR = (gsi_param = 1.1, tau_d_param = 1.0, tau_f_param = 1.0, tau_w_param = 8.0, Ib_param = 1.0)


prob_ode = ODEProblem(LR_EOM!, u0, (0,10000), param_LR)
sol_ode = solve(prob_ode, Tsit5());
plot(sol_ode)

LR_EOM(u, p) = LR_EOM!(similar(u),u,p, 2000)

prob = BifurcationProblem(LR_EOM, sol_ode(2000), param_LR, (@lens _.Ib_param);
    record_from_solution = (x,p) -> (V = x[8]))

opts_br = ContinuationPar(p_min = 1.0, p_max = 1.2, max_steps = 10000,
    newton_options = NewtonPar(tol = 1e-6, max_iterations = 100))

sol = newton(prob, NewtonPar(verbose=true))

br = BifurcationKit.continuation(prob, PALC(), opts_br; normC = norminf, verbosity = 2)

plot(br)
8

did you give up?

9

Haha no, did not give up. Thank you very much for helping me figure out how to get solvable starting parameters for the function. I can reproduce what you did above to get the stable point where the cell repolarizes. I haven’t yet been able to get it to find the unstable fixed point that should appear under the right conditions during the plateau phase of the cardiac AP that causes the oscillations of the early after depolarizations (EADs). I’ve been finding that the solver converges from those initial conditions to bifurcation points up around 330 mV which makes no physical sense, so I’m still in the process of trying to sort that out.

10

Turns out I’d forgotten to set the slow variable in the system to a parameter, but it’s working beautifully now & I was able to replicate the dynamics from the paper as I’d hoped. Thank you so much for your help!