Fig 7.23 - MMSB

Hasty synthetic oscillator model

using OrdinaryDiffEq
using ComponentArrays: ComponentArray
using SimpleUnPack
using CairoMakie

_dx723(u, p, t) = (1 + u.x^2 + p.alpha * p.sigma * u.x^4) / ((1 + u.x^2 + p.sigma * u.x^4) * (1 + u.y^4)) - p.gammax * u.x
_dy723(u, p, t) = p.ax * ((1 + u.x^2 + p.alpha * p.sigma * u.x^4) / ((1 + u.x^2 + p.sigma * u.x^4) * (1 + u.y^4))) - p.gammay * u.y
function model723!(D, u, p, t)
    @unpack alpha, sigma, gammax, gammay, ax = p
    @unpack x, y = u
    D.x = _dx723(u, p, t)
    D.y = _dy723(u, p, t)
    nothing
end

model723! (generic function with 1 method)

ps723 = ComponentArray(
    alpha=11.0,
    sigma=2.0,
    gammax=0.2,
    gammay=0.012,
    ax=0.2,
)

ics723 = ComponentArray(
    x=0.3963,
    y=2.3346,
)

tend = 300.0
prob723 = ODEProblem(model723!, ics723, (0.0, tend), ps723)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(x = 1, y = 2)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 300.0)
u0: ComponentVector{Float64}(x = 0.3963, y = 2.3346)

@time sol723 = solve(prob723, KenCarp47())

  2.667936 seconds (9.36 M allocations: 630.825 MiB, 14.97% gc time, 101.51% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 35-element Vector{Float64}:
   0.0
   0.1137759772306821
   1.251535749537503
   4.728509657194026
   9.398926062359777
  15.765295184700044
  23.0402525522803
  33.517438731491964
  42.87100163244791
  50.3085152864516
   ⋮
 221.81036514721126
 232.71417808246605
 243.61799101772084
 256.71738925440116
 263.7058433577334
 271.4496994209871
 281.0258141849401
 290.6019289488931
 300.0
u: 35-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(x = 1, y = 2)}}}}:
 ComponentVector{Float64}(x = 0.3963, y = 2.3346)
 ComponentVector{Float64}(x = 0.3925291595045122, y = 2.332454815145531)
 ComponentVector{Float64}(x = 0.3578014572116769, y = 2.310864368649203)
 ComponentVector{Float64}(x = 0.28461117805517727, y = 2.245305376069942)
 ComponentVector{Float64}(x = 0.24725396497717678, y = 2.163162524723254)
 ComponentVector{Float64}(x = 0.25715235889067345, y = 2.066916334924095)
 ComponentVector{Float64}(x = 0.30817460110557143, y = 1.9820949412830156)
 ComponentVector{Float64}(x = 0.45276512615690084, y = 1.9216232934247819)
 ComponentVector{Float64}(x = 1.0424482991509572, y = 2.066792564347545)
 ComponentVector{Float64}(x = 1.4116823267075949, y = 2.3362733187517306)
 ⋮
 ComponentVector{Float64}(x = 0.5209317094223064, y = 2.390623550743751)
 ComponentVector{Float64}(x = 0.25582141005310094, y = 2.1909599961214967)
 ComponentVector{Float64}(x = 0.2740927394141302, y = 2.0297851920325023)
 ComponentVector{Float64}(x = 0.4177610762875748, y = 1.92470096674037)
 ComponentVector{Float64}(x = 0.6492668919534699, y = 1.951596560876399)
 ComponentVector{Float64}(x = 1.3853830447461177, y = 2.2246676036774264)
 ComponentVector{Float64}(x = 1.2219113559423684, y = 2.444887775577022)
 ComponentVector{Float64}(x = 0.7972780487074158, y = 2.463653518199335)
 ComponentVector{Float64}(x = 0.38925190396096115, y = 2.3302719613174285)

fig = Figure()
ax = Axis(fig[1, 1],
    xlabel = "Time",
    ylabel = "Concentration",
    title = "Fig 7.23 A"
)
lines!(ax, 0..tend, t-> sol723(t).x, label = "x")
lines!(ax, 0..tend, t-> sol723(t).y, label = "y")
axislegend(ax, position = :rc)
fig

xx = range(0, 1.5, 51)
yy = range(0, 3, 51)
tt = 0:1:tend
∂X723 = [_dx723((; x=x, y=y), ps723, nothing) for x in xx, y in yy]
∂Y723 = [_dy723((; x=x, y=y), ps723, nothing) for x in xx, y in yy]
aa = [sol723(t).x for t in tt]
bb = [sol723(t).y for t in tt];

Vector field

fig = Figure(size=(600, 600))
ax = Axis(fig[1, 1],
    xlabel = "X",
    ylabel = "Y",
    title = "Fig 7.23 B\nVector field",
    aspect = 1,
)
contour!(ax, xx, yy, ∂X723, levels=[0], color=:black, linestyle=:solid, linewidth=2, label="X nullcline")
contour!(ax, xx, yy, ∂Y723, levels=[0], color=:black, linestyle=:dash, linewidth=2, label="Y nullcline")
lines!(ax, aa, bb, color=:tomato, label="Trajectory")
axislegend(ax, position = :rb)
limits!(ax, 0.0, 1.5, 1.5, 3.0)
fig

This notebook was generated using Literate.jl.