Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update the Lotka-Volterra example #4

Merged
merged 1 commit into from
Feb 14, 2024
Merged

Update the Lotka-Volterra example #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 51 additions & 16 deletions docs/src/PredatorPreyTutorial.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,10 @@ EditURL = "../../examples/PredatorPreyTutorial.jl"

# Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model

This example demonstrates how to use `LinearNoiseApproximation.jl` package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic
trajectories obtained using the Gillespie algorithm.
This example demonstrates how to use `LinearNoiseApproximation.jl` package
to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA)
and compare it to the stochastic trajectories obtained using the Gillespie
algorithm.

````@example PredatorPreyTutorial
using LinearNoiseApproximation
Expand All @@ -16,8 +18,15 @@ using Plots
````

## Define the Lotka-Volterra model
The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population $U$ grows at a rate represented by the parameter $\alpha$ in the absence of predators, but decreases as they are consumed by the predator population $V$ at a rate determined by the interaction strength parameter, $\beta$. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, $\delta$.
The Lotka-Volterra model describes the dynamics of predator-prey
interactions in an ecosystem. It assumes that the prey population $U$ grows
at a rate represented by the parameter $\alpha$ in the absence of predators,
but decreases as they are consumed by the predator population $V$ at a rate
determined by the interaction strength parameter, $\beta$. The predator
population decreases in size if they cannot find enough prey to consume,
which is represented by the mortality rate parameter, $\delta$.
The corresponding ODE system reads

```math
\begin{equation}
\begin{aligned}
Expand All @@ -27,6 +36,8 @@ The corresponding ODE system reads
\end{equation}
```

This system of ODEs can be converted to following reaction network.

````@example PredatorPreyTutorial
rn = @reaction_network begin
@parameters α β δ
Expand All @@ -37,33 +48,57 @@ rn = @reaction_network begin
end
````

## Convert the model to a JumpSystem for the Gillespie algorithm
Using `Catalyst.jl` we can convert `ReactionSystem` to a `JumpSystem` which can be used to simulate the stochastic trajectories using the Gillespie algorithm.
## Solve the model using the Linear Noise Approximation (LNA)
First of all, we define the initial conditions, parameters, and time span

````@example PredatorPreyTutorial
jumpsys = convert(JumpSystem, rn)

#define the initial conditions, parameters, and time span
u0 = [120.0, 140.0]
ps = [0.8, 0.005, 0.4]
duration = 30.0
tspan = (0.0, duration)
nothing # hide
````

#solve the model using the Gillespie algorithm
dprob = DiscreteProblem(jumpsys, u0, tspan, ps)
jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))
jsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)
Now we can solve the model, using `DifferentialEquations.jl` and the Runge-
Kutta variant `Vern7` as the solver. We want to save at every 0.5th
timestep.

# solve the model using the Linear Noise Approximation (LNA)
````@example PredatorPreyTutorial
lna = LNASystem(rn)
alg = Vern7()
prob = ODEProblem(lna, u0, tspan, ps)
sol = solve(prob, alg, saveat=0.5, seed=2)
nothing # hide
````

Finally, we save the indices of the mean and variance states for later.

#find the indices of the mean and variance states
````@example PredatorPreyTutorial
mean_idxs, var_idxs = find_states_cov_number([1, 2], lna)
nothing # hide
````

## Convert the model to a JumpSystem for the Gillespie algorithm
Using `Catalyst.jl` we can convert `ReactionSystem` to a `JumpSystem` which
can be used to simulate the stochastic trajectories using the Gillespie
algorithm. This will be our reference to which we compare the solution from
the LNA.

````@example PredatorPreyTutorial
jumpsys = convert(JumpSystem, rn)
dprob = DiscreteProblem(jumpsys, u0, tspan, ps)
jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))
jsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)
nothing # hide
````

## Plotting the results
In this final step, we visualize and compare the results from the LNA and
the stochastic trajectories. The LNA calculated mean solution is plotted as
a solid line (LNA) and the standard deviations around it as a shaded area.
The stochastic jumps are plotted as dots (SSA).

````@example PredatorPreyTutorial
#plot the results
plt = Plots.plot(; size=(800, 350))
color_palette=["#1f78b4" "#ff7f00"]
Expand Down Expand Up @@ -103,8 +138,8 @@ end
plt
````

save the plot
Plots.savefig(plt, "predator_prey_fig.png")
Upon visual inspection, we can see that most stochastic jumps are within the
standard deviation calculated with the LNA.

---

Expand Down
59 changes: 43 additions & 16 deletions examples/PredatorPreyTutorial.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
# # Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model
#
# This example demonstrates how to use `LinearNoiseApproximation.jl` package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic
# trajectories obtained using the Gillespie algorithm.
# This example demonstrates how to use `LinearNoiseApproximation.jl` package
# to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA)
# and compare it to the stochastic trajectories obtained using the Gillespie
# algorithm.

using LinearNoiseApproximation
using DifferentialEquations
Expand All @@ -11,8 +13,15 @@ using Plots


# ## Define the Lotka-Volterra model
# The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population $U$ grows at a rate represented by the parameter $\alpha$ in the absence of predators, but decreases as they are consumed by the predator population $V$ at a rate determined by the interaction strength parameter, $\beta$. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, $\delta$.
# The Lotka-Volterra model describes the dynamics of predator-prey
# interactions in an ecosystem. It assumes that the prey population $U$ grows
# at a rate represented by the parameter $\alpha$ in the absence of predators,
# but decreases as they are consumed by the predator population $V$ at a rate
# determined by the interaction strength parameter, $\beta$. The predator
# population decreases in size if they cannot find enough prey to consume,
# which is represented by the mortality rate parameter, $\delta$.
# The corresponding ODE system reads
#
# ```math
# \begin{equation}
# \begin{aligned}
Expand All @@ -21,6 +30,8 @@ using Plots
# \end{aligned}\quad t\in (0, t_{\text{end}}).
# \end{equation}
# ```
#
# This system of ODEs can be converted to following reaction network.

rn = @reaction_network begin
@parameters α β δ
Expand All @@ -30,33 +41,46 @@ rn = @reaction_network begin
δ, V --> 0
end

# ## Convert the model to a JumpSystem for the Gillespie algorithm
# Using `Catalyst.jl` we can convert `ReactionSystem` to a `JumpSystem` which can be used to simulate the stochastic trajectories using the Gillespie algorithm.

jumpsys = convert(JumpSystem, rn)
## Solve the model using the Linear Noise Approximation (LNA)
# First of all, we define the initial conditions, parameters, and time span

#define the initial conditions, parameters, and time span
u0 = [120.0, 140.0]
ps = [0.8, 0.005, 0.4]
u0 = [120.0, 140.0] # U0, V0
ps = [0.8, 0.005, 0.4] # α, β, δ
duration = 30.0
tspan = (0.0, duration)

#solve the model using the Gillespie algorithm
dprob = DiscreteProblem(jumpsys, u0, tspan, ps)
jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))
jsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)
# Now we can solve the model, using `DifferentialEquations.jl` and the Runge-
# Kutta variant `Vern7` as the solver. We want to save at every 0.5th
# timestep.

## solve the model using the Linear Noise Approximation (LNA)
lna = LNASystem(rn)
alg = Vern7()
prob = ODEProblem(lna, u0, tspan, ps)
sol = solve(prob, alg, saveat=0.5, seed=2)

#find the indices of the mean and variance states
# Finally, we save the indices of the mean and variance states for later.

mean_idxs, var_idxs = find_states_cov_number([1, 2], lna)

# ## Convert the model to a JumpSystem for the Gillespie algorithm
# Using `Catalyst.jl` we can convert `ReactionSystem` to a `JumpSystem` which
# can be used to simulate the stochastic trajectories using the Gillespie
# algorithm. This will be our reference to which we compare the solution from
# the LNA.

jumpsys = convert(JumpSystem, rn)
dprob = DiscreteProblem(jumpsys, u0, tspan, ps)
jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))
jsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)


# ## Plotting the results
# In this final step, we visualize and compare the results from the LNA and
# the stochastic trajectories. The LNA calculated mean solution is plotted as
# a solid line (LNA) and the standard deviations around it as a shaded area.
# The stochastic jumps are plotted as dots (SSA).

#plot the results
plt = Plots.plot(; size=(800, 350))
color_palette=["#1f78b4" "#ff7f00"]
label = ["U" "V"]
Expand Down Expand Up @@ -95,3 +119,6 @@ end
plt
# save the plot
# Plots.savefig(plt, "predator_prey_fig.png")

# Upon visual inspection, we can see that most stochastic jumps are within the
# standard deviation calculated with the LNA.
Loading