Hello!
I am new to Julia and planning to use it (particularly BifurcationKit.jl) as a scriptable alternative to XPPAUT, which I have been using successfully for several years to calculate bifurcation diagrams for nonlinear ODE systems. However, I seem to be running into issues with some aspects of continuation and would be really grateful for any tips on how I might fix them.
For reference, here is the ODE system I am working with (it’s quite stiff):
using BifurcationKit, ModelingToolkit, Plots, DifferentialEquations
using ModelingToolkit: t_nounits as t, D_nounits as D
# Parameters
@parameters a b kwee eps kscyce_prime kscyce_double kdcyce_prime kdcyce_double
@parameters kscyca_prime kscyca_double kdcyca_prime kdcyca_double kdcyca
@parameters Rbtot JRb kprb kdprb alpha SK E2FT kdpe2f kpe2f Kdrbe2f
@parameters ksemi1_prime ksemi1_double kdemi1_prime kdemi1_double kdemi1 Kdc1e1
@parameters Cdh1tot kacdh1 kicdh1_prime kicdh1_double kicdh1
@parameters kscycb_prime kscycb kdcycb_prime kdcycb_double kdcycb
@parameters kpyme_prime kpyme kdpyme Jyme k25_prime k25 kwee_prime
@parameters kicdc20 kacdc20 SAC
@parameters kspolo_prime kspolo kdpolo_prime kdpolo kapolo_prime kapolo kipolo Jpolo
@parameters ENSAtot B55tot kass kdiss kcatB55
@parameters kGwENSA kppx_prime kCdkGwl_prime kCdkGwl kB55Gwl Gwtot
# Dynamical variables
@variables CycE(t) CycA(t) E2FPt(t) Rb(t) Emi1(t) CycB(t) Cdk1(t)
@variables Cdh1(t) Cdc20(t) PoloT(t) Polo(t) pENSAt(t) pGwl(t) PP2AB55(t)
# Algebraic / auxiliary variables
@variables Rbt(t) BB2(t) Comp2(t) E2F(t) BB1(t) Comp1(t) Cdh1t(t) YMEP(t) V25(t) Vwee(t) Vdcycb(t) Complex(t)
# Goldbeter-Koshland function
GB(arg1,arg2,arg3,arg4) = arg2 - arg1 + arg2*arg3 + arg1*arg4
GK(arg1,arg2,arg3,arg4) = 2*arg1*arg4 / (GB(arg1,arg2,arg3,arg4) + sqrt(GB(arg1,arg2,arg3,arg4)^2 - 4*(arg2-arg1)*arg1*arg4))
# Differential equations
eqs = [
D(CycE) ~ kscyce_prime + kscyce_double*E2F - (kdcyce_prime + kdcyce_double*CycA)*CycE,
D(CycA) ~ kscyca_prime + kscyca_double*E2F - (kdcyca_prime + kdcyca_double*Cdc20 + kdcyca*Cdh1)*CycA,
D(E2FPt) ~ kpe2f*(CycA+b*eps*Cdk1)*(E2FT - E2FPt) - kdpe2f*E2FPt,
D(Rb) ~ kdprb*(Rbt - Rb)/(JRb + Rbt - Rb) - kprb*(CycE + CycA + a*eps*Cdk1)*Rb/(JRb + Rb),
D(Emi1) ~ ksemi1_prime + ksemi1_double*E2F - (kdemi1_prime + kdemi1_double*Cdh1 + kdemi1*Polo)*Emi1,
D(CycB) ~ kscycb_prime + kscycb*CycA - Vdcycb*CycB,
D(Cdk1) ~ kscycb_prime + kscycb*CycA + V25*(CycB - Cdk1) - Vwee*Cdk1 - Vdcycb*Cdk1,
D(Cdh1) ~ kacdh1*(Cdh1t - Cdh1) - (kicdh1_prime*CycE + kicdh1_double*CycA + kicdh1*eps*Cdk1)*Cdh1,
D(Cdc20) ~ kacdc20*eps*Cdk1*(1 - Cdc20) - kicdc20*PP2AB55*Cdc20,
D(PoloT) ~ kspolo_prime + kspolo*CycA - (kdpolo_prime + kdpolo*Cdh1)*PoloT,
D(Polo) ~ (kapolo_prime*CycA + kapolo*eps*Cdk1)*(PoloT - Polo)/(Jpolo + PoloT - Polo) - kipolo*Polo/(Jpolo + PoloT - Polo),
D(pENSAt) ~ kGwENSA*pGwl*(ENSAtot - pENSAt) - kcatB55*Complex,
D(pGwl) ~ (kCdkGwl_prime*CycA + kCdkGwl*eps*Cdk1)*(Gwtot - pGwl) - (kppx_prime + kB55Gwl*PP2AB55)*pGwl,
D(PP2AB55) ~ (kdiss + kcatB55)*Complex - kass*(pENSAt - Complex)*(B55tot - Complex)
]
# Algebraic equations
algs = [
Rbt ~ Rbtot/(1 + alpha*SK),
BB2 ~ Rb + E2FT + Kdrbe2f,
Comp2 ~ (BB2 - sqrt(BB2^2 - 4*Rb*E2FT))/2,
E2F ~ (E2FT - E2FPt)*(E2FT - Comp2)/E2FT,
BB1 ~ Emi1 + Cdh1tot + Kdc1e1,
Comp1 ~ (BB1 - sqrt(BB1^2 - 4*Emi1*Cdh1tot))/2,
Cdh1t ~ Cdh1tot - Comp1,
YMEP ~ GK(kpyme_prime*CycA + kpyme*eps*Cdk1, kdpyme*PP2AB55, Jyme, Jyme),
V25 ~ k25_prime + k25*YMEP,
Vwee ~ kwee_prime + kwee*(1 - YMEP),
Vdcycb ~ kdcycb_prime + SAC*kdcycb_double*Cdc20 + kdcycb*Cdh1,
Complex ~ B55tot - PP2AB55
]
# Combine all equations
all_eqs = [eqs; algs]
# Initial conditions (example values, adjust as needed)
u0 = [
CycE => 0.0,
CycA => 0.0,
E2FPt => 0.0,
Rb => 0.0,
Emi1 => 0.0,
CycB => 0.0,
Cdk1 => 0.0,
Cdh1 => 0.0,
Cdc20 => 0.0,
PoloT => 0.0,
Polo => 0.0,
pENSAt => 0.0,
pGwl => 0.0,
PP2AB55 => 1.0
]
# Parameter values (example, use the values from XPPAUT)
pvals = [
a=>1, b=>1, kwee=>1, eps=>1,
kscyce_prime=>0, kscyce_double=>0.1, kdcyce_prime=>0.04, kdcyce_double=>0.1,
kscyca_prime=>0, kscyca_double=>0.03, kdcyca_prime=>0.003, kdcyca_double=>0.05, kdcyca=>0.25,
Rbtot=>1.75, JRb=>0.1, kprb=>1, kdprb=>0.7, alpha=>1, SK=>1,
E2FT=>1, kdpe2f=>0.02, kpe2f=>0.1, Kdrbe2f=>0.001,
ksemi1_prime=>0, ksemi1_double=>0.1, kdemi1_prime=>0.012, kdemi1_double=>0.2, kdemi1=>0.5, Kdc1e1=>0.0175,
Cdh1tot=>1, kacdh1=>1, kicdh1_prime=>1, kicdh1_double=>2, kicdh1=>400,
kscycb_prime=>0, kscycb=>0.02, kdcycb_prime=>0.004, kdcycb_double=>0.05, kdcycb=>0.25,
kpyme_prime=>0, kpyme=>2, kdpyme=>0.4, Jyme=>0.1, k25_prime=>0.03, k25=>1, kwee_prime=>0.01,
kicdc20=>1, kacdc20=>0.2, SAC=>1,
kspolo_prime=>0.01, kspolo=>0, kdpolo_prime=>0.01, kdpolo=>1,
kapolo_prime=>0, kapolo=>1, kipolo=>0.5, Jpolo=>0.01,
ENSAtot=>4, B55tot=>1, kass=>500, kdiss=>0.3, kcatB55=>1,
kGwENSA=>1.5, kppx_prime=>0.4, kCdkGwl_prime=>0, kCdkGwl=>2, kB55Gwl=>4, Gwtot=>1
]
# Build the system
@named cellcycle_system = ODESystem(all_eqs, t)
# Simplify / complete for ODEProblem
sys_simplified = structural_simplify(cellcycle_system)
And here is my attempt at running a bifurcation problem. First, I set the parameter I want to do continuation over to a value for which the system doesn’t oscillate, I run a simulation to steady state and use the vector of values at the last time point as guess for the bifurcation problem:
pvals2 = [
a=>1, b=>1, kwee=>1, eps=>1,
kscyce_prime=>0, kscyce_double=>0.1, kdcyce_prime=>0.04, kdcyce_double=>0.1,
kscyca_prime=>0, kscyca_double=>0.03, kdcyca_prime=>0.003, kdcyca_double=>0.05, kdcyca=>0.25,
Rbtot=>1.75, JRb=>0.1, kprb=>1, kdprb=>0.7, alpha=>1, SK=>1,
E2FT=>1, kdpe2f=>0.02, kpe2f=>0.1, Kdrbe2f=>0.001,
ksemi1_prime=>0, ksemi1_double=>0.1, kdemi1_prime=>0.012, kdemi1_double=>0.2, kdemi1=>0.5, Kdc1e1=>0.0175,
Cdh1tot=>1, kacdh1=>1, kicdh1_prime=>1, kicdh1_double=>2, kicdh1=>400,
kscycb_prime=>0, kscycb=>0.02, kdcycb_prime=>0.004, kdcycb_double=>0.05, kdcycb=>0.25,
kpyme_prime=>0, kpyme=>2, kdpyme=>0.4, Jyme=>0.1, k25_prime=>0.03, k25=>1, kwee_prime=>0.01,
kicdc20=>1, kacdc20=>0.2, SAC=>1,
kspolo_prime=>0.01, kspolo=>0, kdpolo_prime=>0.01, kdpolo=>1,
kapolo_prime=>0, kapolo=>1, kipolo=>0.5, Jpolo=>0.01,
ENSAtot=>4, B55tot=>5, kass=>500, kdiss=>0.3, kcatB55=>1,
kGwENSA=>1.5, kppx_prime=>0.4, kCdkGwl_prime=>0, kCdkGwl=>2, kB55Gwl=>4, Gwtot=>1
]
prob2 = ODEProblem(sys_simplified, u0, (0.0, 4000.0), pvals2)
sol2 = solve(prob2, Rosenbrock23())
plot(sol2, vars=[CycB, Cdk1, Cdc20, PP2AB55],
xlabel="Time", ylabel="Concentration",
lw=3)
bprob = BifurcationProblem(sys_simplified, sol2.u[end], pvals2, B55tot, plot_var = Cdk1)
br = continuation(bprob, PALC(),
ContinuationPar(p_min = -1., p_max = 10., ds = -1e-3, dsmin = 1e-9, dsmax = 1e-2, max_steps = 100000),
normC = norminf)
After running this continuation problem, 3 Hopf (there should be 4) and 2 saddle node points are found. Subsequently, I tried to run continuations from the Hopf points with various sets of settings, and I just can’t seem to get the right output. Here is one of the ‘better’ attempts:
newton_opts_po = NewtonPar(tol = 1e-5, max_iterations = 8)
br_po = continuation(br, 5,
ContinuationPar(p_min = -1., p_max = 5., ds = -0.0001, dsmin = 1e-4,
dsmax = 0.1, max_steps = 10000,
newton_options = newton_opts_po),
PeriodicOrbitOCollProblem(50, 3; meshadapt = true);
usedeflation = true, plot = true)
The problem is that, when plotted, the limit cycle continuation starts very far away from the Hopf bifurcation. For reference, this is an AUTO plot that shows the expected output for continuation from the same Hopf point.
