using AlgebraicDynamics
using Catlab
using LabelledArrays
using OrdinaryDiffEq
using Plots, Plots.PlotMeasures
There will be 16 parameters in to the total ecosystem. With four distinct types:
params = LVector(αr=0.3, βrf=0.015, γrf=0.015, δf=0.7,
βrh=0.01, γrh=0.01, δh=0.5,
γfishh=0.001, βfishh=0.003,
αfish=0.35, βfishF=0.015, γfishF=0.015,
δF=0.7, βFs=0.017, γFs=0.017, δs=0.35)
16-element LabelledArrays.LArray{Float64, 1, Vector{Float64}, (:αr, :βrf, :γrf, :δf, :βrh, :γrh, :δh, :γfishh, :βfishh, :αfish, :βfishF, :γfishF, :δF, :βFs, :γFs, :δs)}: :αr => 0.3 :βrf => 0.015 :γrf => 0.015 :δf => 0.7 :βrh => 0.01 :γrh => 0.01 ⋮ :γfishF => 0.015 :δF => 0.7 :βFs => 0.017 :γFs => 0.017 :δs => 0.35
A standard Lotka Volterra predator-prey model is the composition of three primitive resource sharers:
However, there are not two independent rabbit populations -- one that grows and one that gets eaten by foxes. Likewise, there are not two independent fox populations -- one that declines and one that feasts on rabbits. To capture these interactions between the trio of resource sharers, we compose them by identifying the exposed rabbit populations and identifying the exposed fox populations. The syntax for this undirected composition is defined by an undirected wiring diagram.
# Define the primitive systems
dotr(u,p,t) = p.αr*u
dotrf(u,p,t) = [-p.βrf*u[1]*u[2], p.γrf*u[1]*u[2]]
dotf(u,p,t) = -p.δf*u
rabbit_growth = ContinuousResourceSharer{Float64}(1, dotr)
rabbitfox_predation = ContinuousResourceSharer{Float64}(2, dotrf)
fox_decline = ContinuousResourceSharer{Float64}(1, dotf)
# Define the composition pattern
rabbitfox_pattern = @relation (rabbits, foxes) begin
rabbit_growth(rabbits)
rabbitfox_predation(rabbits,foxes)
fox_decline(foxes)
end
# Compose
submodels = Dict(
:rabbit_growth => rabbit_growth,
:rabbitfox_predation => rabbitfox_predation,
:fox_decline => fox_decline
)
rabbitfox_system = oapply(rabbitfox_pattern, submodels)
Previously, when we derived the Lotka-Volterra model via undirected composition, we by-hand defined the undirected wiring diagram that implements the composition pattern. In contrast, here we implement the same composition pattern as before but this time using the @relation
macro. This strategy simplifies the definition and explicitly names the boxes and variables. We visualize the composition pattern below.
to_graphviz(rabbitfox_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))
We can now construct an ODEProblem
from the resource sharer rabbitfox_system
and plot the solution.
u0 = [10.0, 100.0]
tspan = (0.0, 100.0)
prob = ODEProblem(rabbitfox_system, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, rabbitfox_system,
lw=2,
title = "Lotka-Volterra Predator-Prey Model",
xlabel = "time", ylabel = "population size")
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1
Suppose we now have a three species ecosystem containing rabbits, foxes, and hawks. Foxes and hawks both prey upon rabbits but do not interact with each other. This ecosystem consists of five primitive systems which share variables.
This means the desired composition pattern has five boxes and many ports and wires to keep track of. Instead of implementing this composition pattern by hand, we construct it as a pushout.
# Define the composition pattern for rabbit growth
rabbit_pattern = @relation (rabbits,) -> rabbit_growth(rabbits)
# Define the composition pattern for the rabbit/hawk Lotka Volterra model
rabbithawk_pattern = @relation (rabbits, hawks) begin
rabbit_growth(rabbits)
rabbithawk_predation(rabbits,hawks)
hawk_decline(hawks)
end
# Define transformations between the composition patterns
rabbitfox_transform = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbitfox_pattern)
rabbithawk_transform = ACSetTransformation((Box=[1], Junction=[1], Port=[1], OuterPort=[1]), rabbit_pattern, rabbithawk_pattern)
# Take the pushout to define the composition pattern for the rabbit, fox, hawk system
rabbitfoxhawk_pattern = ob(pushout(rabbitfox_transform, rabbithawk_transform))
# Visualize the compsition pattern
to_graphviz(rabbitfoxhawk_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".9"))
# Define the additional primitive systems
dotrh(u, p, t) = [-p.βrh*u[1]*u[2], p.γrh*u[1]*u[2]]
doth(u, p, t) = -p.δh*u
rabbithawk_predation = ContinuousResourceSharer{Float64}(2, dotrh)
hawk_decline = ContinuousResourceSharer{Float64}(1, doth)
# Compose
merge!(submodels, Dict(
:rabbithawk_predation => rabbithawk_predation,
:hawk_decline => hawk_decline
))
land_system = oapply(rabbitfoxhawk_pattern, submodels)
# Solve and plot
u0 = [10.0, 100.0, 50.0]
tspan = (0.0, 100.0)
prob = ODEProblem(land_system, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, land_system,
lw=2,
title = "Land Ecosystem",
xlabel = "time", ylabel = "population size")
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1
Unfortunately, the hawks are going extinct in this model. We'll have to give hawks something else to eat!
Consider a ocean ecosystem containing three species —- little fish, big fish, and sharks -— with two predation interactions —- sharks eat big fish and big fish eat little fish.
This ecosystem can be modeled as the composition of 3 machines:
# Define the primitive systems
dotfish(f, x, p, t) = [p.αfish*f[1] - p.βfishF*x[1]*f[1]]
dotFISH(F, x, p, t) = [p.γfishF*x[1]*F[1] - p.δF*F[1] - p.βFs*x[2]*F[1]]
dotsharks(s, x, p, t) = [p.γFs*s[1]*x[1]-p.δs*s[1]]
fish = ContinuousMachine{Float64}(1,1,1, dotfish, (f,p,t)->f)
FISH = ContinuousMachine{Float64}(2,1,1, dotFISH, (F,p,t)->F)
sharks = ContinuousMachine{Float64}(1,1,1, dotsharks, (s,p,t)->s)
We compose these machines by (1) sending the output of the big fish machine as the input to both the little fish and shark machines and (2) sending the output of the little fish and shark machines as the inputs to the big fish machine. The syntax for this directed composition is given by a directed wiring diagram.
# Define the composition pattern
ocean_pattern = WiringDiagram([], [:fish, :Fish, :shark])
fish_box = add_box!(ocean_pattern, Box(:fish, [:pop], [:pop]))
Fish_box = add_box!(ocean_pattern, Box(:Fish, [:pop, :pop], [:pop]))
shark_box = add_box!(ocean_pattern, Box(:shark, [:pop], [:pop]))
output_box = output_id(ocean_pattern)
add_wires!(ocean_pattern, Pair[
(fish_box, 1) => (Fish_box, 1),
(shark_box, 1) => (Fish_box, 2),
(Fish_box, 1) => (fish_box, 1),
(Fish_box, 1) => (shark_box, 1),
(fish_box, 1) => (output_box, 1),
(Fish_box, 1) => (output_box, 2),
(shark_box, 1) => (output_box, 3)
])
# Visualize the composition pattern
to_graphviz(ocean_pattern, orientation=TopToBottom)
# Compose
ocean_system = oapply(ocean_pattern, [fish, FISH, sharks])
# Solve and plot
u0 = [100.0, 10, 2.0]
tspan = (0.0, 100.0)
prob = ODEProblem(ocean_system, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ocean_system,
lw=2, title = "Ocean Ecosystem",
xlabel = "time", ylabel = "population size"
)
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1
We will introduce a final predation interaction -- hawks eat little fish -- which will combine the land and ocean ecosystems.
The composition will be as resource shareres so the first thing we will do is use the dynamics of the machine ocean_system
to define the dynamics of a resource sharer. We will also define a resource sharer that models hawk/little fish predation.
# Define the additional primitive systems
ocean_system_rs = ContinuousResourceSharer{Float64}(3, (u,p,t) -> eval_dynamics(ocean_system, u, [], p))
dothf(u,p,t) = [p.γfishh*u[1]*u[2], -p.βfishh*u[1]*u[2]]
fishhawk_predation = ContinuousResourceSharer{Float64}(2, dothf)
# Define the composition pattern
eco_pattern = @relation (rabbits, foxes, hawks, littlefish, BigFish, sharks) begin
turf(rabbits,foxes,hawks)
air(hawks, littlefish)
surf(littlefish, BigFish, sharks)
end
# Visualize the composition pattern
to_graphviz(eco_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))
# Compose
ecosystem = oapply(eco_pattern, Dict(
:turf => land_system,
:air => fishhawk_predation,
:surf => ocean_system_rs)
)
We can now plot the evolution of the total ecosystem.
# Solve and plot
u0 = [100.0, 50.0, 20.0, 100, 10, 2.0]
tspan = (0.0, 100.0)
prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
lw=2, title = "Total Ecosystem",
xlabel = "time", ylabel = "population size"
)
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1
Let's zoom in on a narrower time-window.
tspan = (0.0, 30.0)
prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
lw=2, title = "Total Ecosystem",
xlabel = "time", ylabel = "population size"
)
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1
As a sanity check we can define the rates for the hawk/little fish predation to be 0. This decouples the land and ocean ecosystems. As expected, the plot shows the original evolution of the land ecosystem overlayed with the original evolution of the ocean ecosystem. This shows that they two ecosystems now evolve independently.
tspan = (0.0, 100.0)
params.βfishh = 0; params.γfishh = 0
prob = ODEProblem(ecosystem, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, ecosystem,
lw=2, title = "Decoupled Ecosystem",
xlabel = "time", ylabel = "population size"
)
┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead. │ caller = ip:0x0 └ @ Core :-1