Hi all! I’m trying to set up a pipeline for creating an ODE system based on a network of interactions. As a simple example, the corresponding system for the following network would be as such:
\frac{dA}{dt}=-A+\Phi(FB)
\\
\frac{dB}{dt}=-B+\Phi(A)
Where A and B are processes that depend on and amplify each other, and F is some parameter (that I have been trying to bifurcate with). \Phi(x) is a sigmoid function, in my case I use:
Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))
This is my code for generating my custom ODE equations of motion:
# create equations of motion
function create_viability_model(odenetworkmodel::ODENetworkModel)
#should be in place
return function model!(du, u, p, t=0)
Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))
for i in range(1, odenetworkmodel.num_nodes)
du[i] = -u[i]
input = 1.0
for j in range(1, odenetworkmodel.num_nodes)
if odenetworkmodel.adj_matrix[i,j] == 1
input *= u[j]
end
end
for k in range(1, length(odenetworkmodel.parameters))
if odenetworkmodel.par_adj_matrix[k,i] != 0
input *= p[k]
end
end
du[i] += Φ(input)
end
return du
end
end
ODENetworkModel is a struct I created that holds the information that the user inputs to make such a model. If odenetworkmodel.adj_matrix[i,j] == 1, then node i depends on node j, so we multiply the input by the current value of node j (I’ve been instructed that if A depends on C, D, E, then the product of the inputs is used for the sigmoid function, i.e. \frac{dA}{dt}=-A+\Phi(CDE). The parameter adjacency matrix has a row for each parameter, so for the kth parameter, if the ith column == 1, then that parameter is being fed into node i.
struct ODENetworkModel
num_nodes::Int64
adj_matrix::Array{Int64}
par_adj_matrix::Array{Int64}
parameters::Array{Float64}
p_index::Int64
p_start::Float64
p_end::Float64
u₀::Array{Float64}
end
I have been able to use interactive_trajectory from DynamicalSystems.jl to visualise the trajectories for when there are 2 nodes (such as the example above), and even find_attractors(ds) from Attractors.jl seems to be doing fine and finding attractors that make sense for the parameters and context, BUT I can’t seem to get BifurcationKit.continuation to work! I’ve been wrestling with some dimension errors, and am currently stuck on a conversion error:
function create_bifurcation_diagram(problem)
p = problem.parameters
u₀ = problem.u₀
model! = create_viability_model(problem)
ds = ContinuousDynamicalSystem(model!, u₀, p)
# initial_conditions = vec(collect.(Iterators.product(0:0.1:2, 0:0.1:2)))
attractors = find_attractors(ds)
equilibrium_guess = vec(attractors[1])
!line 76! bfprob = BifurcationProblem(model!, equilibrium_guess, vec(p), (@lens _[problem.p_index]), record_from_solution = recordFromSolution)
println(bfprob)
# continuation options
opts_br = ContinuationPar(p_min = problem.p_start, p_max = problem.p_end)
br = BifurcationKit.continuation(bfprob, PALC(), opts_br; bothside = true)
# scene = Plots.plot(diagram; code = (), title="Viability Bifurcation Diagram", ylim = (0.0,1.0), color=:blue, ylabel = "A")
# display(scene)
end
function main()
# don't load from file- does that fix it?
model = ODENetworkModel(2, [0 1; 1 0], [1 0; 0 0], [0.4, 0.9], 1, 0.1, 0.9, [0.4, 0.5])
create_bifurcation_diagram(model)
end
Mapping initial conditions to attractors: 100%|███████████████████████████████████████████████████████████████| Time: 0:00:00
ERROR: MethodError: Cannot `convert` an object of type Float64 to an object of type SVector{2, Float64}
Closest candidates are:
convert(::Type{SVector{N, T}}, ::CartesianIndex{N}) where {N, T}
@ StaticArrays C:\Users\smahn\.julia\packages\StaticArrays\yXGNL\src\SVector.jl:8
convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
@ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
@ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
...
Stacktrace:
[1] create_bifurcation_diagram(problem::ODENetworkModel)
@ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:76
[2] main()
@ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:93
[3] top-level scope
@ REPL[5]:1
Does anyone have any tips for when you’re defining an arbitrarily large system like this? When I explicitly write out the equations of motion, I can get BifurcationKit.continuation running smoothly, but not with this method.
Many thanks for your help!!