Lotka-Volterra Model

using AlgebraicPetri

using OrdinaryDiffEq
using Plots

using Catlab
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Programs.RelationalPrograms

display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

Step 1: Define the building block Petri nets needed to construct the model

birth_petri = Open(PetriNet(1, 1=>(1,1)));
Graph(birth_petri)
G S_1 S_1 T_1 T_1 S_1->T_1 1 T_1->S_1 2
predation_petri = Open(PetriNet(2, (1,2)=>(2,2)));
Graph(predation_petri)
G S_1 S_1 T_1 T_1 S_1->T_1 1 S_2 S_2 S_2->T_1 1 T_1->S_2 2
death_petri = Open(PetriNet(1, 1=>()));
Graph(death_petri)
G S_1 S_1 T_1 T_1 S_1->T_1 1

Step 2: Generate models using a relational syntax

lotka_volterra = @relation (wolves, rabbits) begin
  birth(rabbits)
  predation(rabbits, wolves)
  death(wolves)
end
display_uwd(lotka_volterra)
G n1 birth n7 rabbits n1--n7 n2 predation n6 wolves n2--n6 n2--n7 n3 death n3--n6 n4--n6 n5--n7
lv_dict = Dict(:birth=>birth_petri, :predation=>predation_petri, :death=>death_petri);
lotka_petri = apex(oapply(lotka_volterra, lv_dict))
Graph(lotka_petri)
G S_1 S_1 T_1 T_1 S_1->T_1 1 T_2 T_2 S_1->T_2 1 S_2 S_2 S_2->T_2 1 T_3 T_3 S_2->T_3 1 T_1->S_1 2 T_2->S_2 2

Generate appropriate vector fields, define parameters, and visualize solution

u0 = [100, 10];
p = [.3, .015, .7];
prob = ODEProblem(vectorfield(lotka_petri),u0,(0.0,100.0),p);
sol = solve(prob,Tsit5(),abstol=1e-8);
plot(sol)

Step 3: Extend your model to handle more complex phenomena

such as a small food chain between little fish, big fish, and sharks

dual_lv = @relation (fish, Fish, Shark) begin
  birth(fish)
  predation(fish, Fish)
  death(Fish)
  predation(Fish, Shark)
  death(Shark)
end
display_uwd(dual_lv)
G n1 birth n9 fish n1--n9 n2 predation n2--n9 n10 Fish n2--n10 n3 death n3--n10 n4 predation n4--n10 n11 Shark n4--n11 n5 death n5--n11 n6--n9 n7--n10 n8--n11
dual_lv_petri = apex(oapply(dual_lv, lv_dict))
Graph(dual_lv_petri)
G S_1 S_1 T_1 T_1 S_1->T_1 1 T_2 T_2 S_1->T_2 1 S_2 S_2 S_2->T_2 1 T_3 T_3 S_2->T_3 1 T_4 T_4 S_2->T_4 1 S_3 S_3 S_3->T_4 1 T_5 T_5 S_3->T_5 1 T_1->S_1 2 T_2->S_2 2 T_4->S_3 2

Generate a new solver, provide parameters, and analyze results

u0 = [100, 10, 2];
p = [.3, .015, .7, .017, .35];
prob = ODEProblem(vectorfield(dual_lv_petri),u0,(0.0,100.0),p);
sol = solve(prob,Tsit5(),abstol=1e-6);
plot(sol)