Lotka-Volterra Model

In [1]:
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

In [2]:
birth_petri = Open(PetriNet(1, 1=>(1,1)));
Graph(birth_petri)
Out[2]:
G s1 1 t1 1 s1->t1 1 t1->s1 2
In [3]:
predation_petri = Open(PetriNet(2, (1,2)=>(2,2)));
Graph(predation_petri)
Out[3]:
G s1 1 t1 1 s1->t1 1 s2 2 s2->t1 1 t1->s2 2
In [4]:
death_petri = Open(PetriNet(1, 1=>()));
Graph(death_petri)
Out[4]:
G s1 1 t1 1 s1->t1 1

Step 2: Generate models using a relational syntax

In [5]:
lotka_volterra = @relation (wolves, rabbits) begin
  birth(rabbits)
  predation(rabbits, wolves)
  death(wolves)
end
display_uwd(lotka_volterra)
Out[5]:
G n1 birth n7 rabbits n1--n7 n2 predation n6 wolves n2--n6 n2--n7 n3 death n3--n6 n4--n6 n5--n7
In [6]:
lv_dict = Dict(:birth=>birth_petri, :predation=>predation_petri, :death=>death_petri);
lotka_petri = apex(oapply(lotka_volterra, lv_dict))
Graph(lotka_petri)
Out[6]:
G s1 1 t1 1 s1->t1 1 t2 2 s1->t2 1 s2 2 s2->t2 1 t3 3 s2->t3 1 t1->s1 2 t2->s2 2

Generate appropriate vector fields, define parameters, and visualize solution

In [7]:
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)
Out[7]:

Step 3: Extend your model to handle more complex phenomena

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

In [8]:
dual_lv = @relation (fish, Fish, Shark) begin
  birth(fish)
  predation(fish, Fish)
  death(Fish)
  predation(Fish, Shark)
  death(Shark)
end
display_uwd(dual_lv)
Out[8]:
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
In [9]:
dual_lv_petri = apex(oapply(dual_lv, lv_dict))
Graph(dual_lv_petri)
Out[9]:
G s1 1 t1 1 s1->t1 1 t2 2 s1->t2 1 s2 2 s2->t2 1 t3 3 s2->t3 1 t4 4 s2->t4 1 s3 3 s3->t4 1 t5 5 s3->t5 1 t1->s1 2 t2->s2 2 t4->s3 2

Generate a new solver, provide parameters, and analyze results

In [10]:
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)
Out[10]: