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/model_creation/examples/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md index 8e9dd34ecb..5ec69e5467 100644 --- a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md +++ b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md @@ -6,61 +6,73 @@ The Smoluchowski coagulation equation describes a system of reactions in which m We begin by importing some necessary packages. ```julia 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 +# 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)) ``` 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 # possible pairs of reactant multimers -pair = [] +pair = Int[] 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 ``` 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 # 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 t = default_t() -@species k[1:nr] (X(t))[1:N] +@parameters k[1:nr] +@species (X(t))[1:N] pars = Pair.(collect(k), kv) # 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 @@ -82,19 +94,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 # 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) +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 # 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 +127,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