diff --git a/docs/src/PredatorPreyTutorial.md b/docs/src/PredatorPreyTutorial.md index 92ca5d8..98d9d0e 100644 --- a/docs/src/PredatorPreyTutorial.md +++ b/docs/src/PredatorPreyTutorial.md @@ -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 @@ -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} @@ -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 α β δ @@ -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"] @@ -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. --- diff --git a/examples/PredatorPreyTutorial.jl b/examples/PredatorPreyTutorial.jl index b4a70b2..6030722 100644 --- a/examples/PredatorPreyTutorial.jl +++ b/examples/PredatorPreyTutorial.jl @@ -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 @@ -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} @@ -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 α β δ @@ -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"] @@ -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.