diff --git a/docs/pages.jl b/docs/pages.jl index 69c81e3216..ea66c7e088 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -19,7 +19,7 @@ pages = Any[ "model_creation/examples/basic_CRN_library.md", "model_creation/examples/programmatic_generative_linear_pathway.md", "model_creation/examples/hodgkin_huxley_equation.md", - #"model_creation/examples/smoluchowski_coagulation_equation.md" + "model_creation/examples/smoluchowski_coagulation_equation.md" ] ], "Model simulation" => Any[ @@ -32,7 +32,7 @@ pages = Any[ "Steady state analysis" => Any[ "steady_state_functionality/homotopy_continuation.md", "steady_state_functionality/nonlinear_solve.md", - "steady_state_functionality/steady_state_stability_computation.md", + "steady_state_functionality/steady_state_stability_computation.md", "steady_state_functionality/bifurcation_diagrams.md", "steady_state_functionality/dynamical_systems.md" ], diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml index ebb086fda6..83ecae3cc6 100644 --- a/docs/src/assets/Project.toml +++ b/docs/src/assets/Project.toml @@ -5,7 +5,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c" -DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" @@ -38,7 +37,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] BenchmarkTools = "1.5" -BifurcationKit = "0.3" +BifurcationKit = "0.3.4" CairoMakie = "0.12" Catalyst = "13" DataFrames = "1.6" diff --git a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md index 8e9dd34ecb..23fa2c4f63 100644 --- a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md +++ b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md @@ -4,72 +4,88 @@ This tutorial shows how to programmatically construct a [`ReactionSystem`](@ref) The Smoluchowski coagulation equation describes a system of reactions in which monomers may collide to form dimers, monomers and dimers may collide to form trimers, and so on. This models a variety of chemical/physical processes, including polymerization and flocculation. We begin by importing some necessary packages. -```julia +```@example smcoag1 using ModelingToolkit, Catalyst, LinearAlgebra -using DiffEqBase, JumpProcesses +using JumpProcesses using Plots, SpecialFunctions ``` Suppose the maximum cluster size is `N`. We assume an initial concentration of monomers, `Nₒ`, and let `uₒ` denote the initial number of monomers in the system. We have `nr` total reactions, and label by `V` the bulk volume of the system (which plays an important role in the calculation of rate laws since we have bimolecular reactions). Our basic parameters are then -```julia -## Parameter -N = 10 # maximum cluster size -Vₒ = (4π/3)*(10e-06*100)^3 # volume of a monomers in cm³ -Nₒ = 1e-06/Vₒ # initial conc. = (No. of init. monomers) / bulk volume -uₒ = 10000 # No. of monomers initially -V = uₒ/Nₒ # Bulk volume of system in cm³ - -integ(x) = Int(floor(x)) -n = integ(N/2) -nr = N%2 == 0 ? (n*(n + 1) - n) : (n*(n + 1)) # No. of forward reactions +```@example smcoag1 +# maximum cluster size +N = 10 + +# volume of a monomers in cm³ +Vₒ = (4π / 3) * (10e-06 * 100)^3 + +# initial conc. = (No. of init. monomers) / bulk volume +Nₒ = 1e-06 / Vₒ + +# No. of monomers initially +uₒ = 10000 + +# Bulk volume of system in cm³ +V = uₒ / Nₒ +n = floor(Int, N / 2) + +# No. of forward reactions +nr = ((N % 2) == 0) ? (n*(n + 1) - n) : (n*(n + 1)) +nothing #hide ``` The [Smoluchowski coagulation equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation) Wikipedia page illustrates the set of possible reactions that can occur. We can easily enumerate the `pair`s of multimer reactants that can combine when allowing a maximal cluster size of `N` monomers. We initialize the volumes of the reactant multimers as `volᵢ` and `volⱼ` -```julia +```@example smcoag1 # possible pairs of reactant multimers pair = [] for i = 2:N - push!(pair, [1:integ(i/2) i .- (1:integ(i/2))]) + halfi = floor(Int, i/2) + push!(pair, [(1:halfi) (i .- (1:halfi))]) end pair = vcat(pair...) -vᵢ = @view pair[:,1] # Reactant 1 indices -vⱼ = @view pair[:,2] # Reactant 2 indices -volᵢ = Vₒ*vᵢ # cm⁻³ -volⱼ = Vₒ*vⱼ # cm⁻³ +vᵢ = @view pair[:, 1] # Reactant 1 indices +vⱼ = @view pair[:, 2] # Reactant 2 indices +volᵢ = Vₒ * vᵢ # cm⁻³ +volⱼ = Vₒ * vⱼ # cm⁻³ sum_vᵢvⱼ = @. vᵢ + vⱼ # Product index +nothing #hide ``` We next specify the rates (i.e. kernel) at which reactants collide to form products. For simplicity, we allow a user-selected additive kernel or constant kernel. The constants(`B` and `C`) are adopted from Scott's paper [2](https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469_1968_025_0054_asocdc_2_0_co_2.xml) -```julia +```@example smcoag1 # set i to 1 for additive kernel, 2 for constant i = 1 -if i==1 - B = 1.53e03 # s⁻¹ - kv = @. B*(volᵢ + volⱼ)/V # dividing by volume as its a bi-molecular reaction chain +if i == 1 + B = 1.53e03 # s⁻¹ + + # dividing by volume as it is a bimolecular reaction chain + kv = @. B * (volᵢ + volⱼ) / V elseif i==2 - C = 1.84e-04 # cm³ s⁻¹ - kv = fill(C/V, nr) + C = 1.84e-04 # cm³ s⁻¹ + kv = fill(C / V, nr) end ``` -We'll store the reaction rates in `pars` as `Pair`s, and set the initial condition that only monomers are present at ``t=0`` in `u₀map`. -```julia -# unknown variables are X, pars stores rate parameters for each rx +We'll set the parameters and the initial condition that only monomers are present at ``t=0`` in `u₀map`. +```@example smcoag1 +# k is a vector of the parameters, with values given by the vector kv +@parameters k[1:nr] = kv + +# create the vector of species X_1,...,X_N t = default_t() -@species k[1:nr] (X(t))[1:N] -pars = Pair.(collect(k), kv) +@species (X(t))[1:N] # time-span if i == 1 - tspan = (0. ,2000.) + tspan = (0.0, 2000.0) elseif i == 2 - tspan = (0. ,350.) + tspan = (0.0, 350.0) end # initial condition of monomers u₀ = zeros(Int64, N) u₀[1] = uₒ -u₀map = Pair.(collect(X), u₀) # map variable to its initial value +u₀map = Pair.(collect(X), u₀) # map species to its initial value +nothing #hide ``` Here we generate the reactions programmatically. We systematically create Catalyst `Reaction`s for each possible reaction shown in the figure on [Wikipedia](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation). When `vᵢ[n] == vⱼ[n]`, we set the stoichiometric coefficient of the reactant multimer to two. -```julia +```@example smcoag1 # vector to store the Reactions in rx = [] for n = 1:nr @@ -82,19 +98,20 @@ for n = 1:nr end end @named rs = ReactionSystem(rx, t, collect(X), collect(k)) +rs = complete(rs) ``` We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/jump_types/) documentation -```julia +```@example smcoag1 # solving the system -jumpsys = convert(JumpSystem, rs) -dprob = DiscreteProblem(jumpsys, u₀map, tspan, pars) -jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false,false)) -jsol = solve(jprob, SSAStepper(), saveat = tspan[2]/30) +jumpsys = complete(convert(JumpSystem, rs)) +dprob = DiscreteProblem(rs, u₀map, tspan) +jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions = (false, false)) +jsol = solve(jprob, SSAStepper(), saveat = tspan[2] / 30) ``` Lets check the results for the first three polymers/cluster sizes. We compare to the analytical solution for this system: -```julia +```@example smcoag1 # Results for first three polymers...i.e. monomers, dimers and trimers -v_res = [1;2;3] +v_res = [1; 2; 3] # comparison with analytical solution # we only plot the stochastic solution at a small number of points @@ -114,15 +131,15 @@ elseif i == 2 end # plotting normalised concentration vs analytical solution -default(lw=2, xlabel="Time (sec)") -scatter(ϕ, jsol(t)[1,:]/uₒ, label="X1 (monomers)", markercolor=:blue) +default(lw = 2, xlabel = "Time (sec)") +scatter(ϕ, jsol(t)[1,:] / uₒ, label = "X1 (monomers)", markercolor = :blue) plot!(ϕ, sol[1,:]/Nₒ, line = (:dot,4,:blue), label="Analytical sol--X1") -scatter!(ϕ, jsol(t)[2,:]/uₒ, label="X2 (dimers)", markercolor=:orange) -plot!(ϕ, sol[2,:]/Nₒ, line = (:dot,4,:orange), label="Analytical sol--X2") +scatter!(ϕ, jsol(t)[2,:] / uₒ, label = "X2 (dimers)", markercolor = :orange) +plot!(ϕ, sol[2,:] / Nₒ, line = (:dot, 4, :orange), label = "Analytical sol--X2") -scatter!(ϕ, jsol(t)[3,:]/uₒ, label="X3 (trimers)", markercolor=:purple) -plot!(ϕ, sol[3,:]/Nₒ, line = (:dot,4,:purple), label="Analytical sol--X3", +scatter!(ϕ, jsol(t)[3,:] / uₒ, label = "X3 (trimers)", markercolor = :purple) +plot!(ϕ, sol[3,:] / Nₒ, line = (:dot, 4, :purple), label = "Analytical sol--X3", ylabel = "Normalized Concentration") ``` For the **additive kernel** we find