Skip to content

Commit

Permalink
Merge pull request #48 from ChrisRackauckas/patch-1
Browse files Browse the repository at this point in the history
Check Catalyst 13
  • Loading branch information
augustinas1 authored May 6, 2023
2 parents 0e89f6f + c3ca652 commit e35caf7
Show file tree
Hide file tree
Showing 35 changed files with 333 additions and 419 deletions.
7 changes: 7 additions & 0 deletions .github/dependabot.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
version: 2
updates:
- package-ecosystem: "github-actions"
directory: "/" # Location of package manifests
schedule:
interval: "weekly"
4 changes: 2 additions & 2 deletions .github/workflows/Documentation.yml
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ jobs:
build:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v2
- uses: actions/checkout@v3
- uses: julia-actions/setup-julia@latest
with:
version: '1'
Expand All @@ -20,5 +20,5 @@ jobs:
- name: Build and deploy
env:
GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
#DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
run: julia --project=docs/ docs/make.jl
6 changes: 3 additions & 3 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,12 @@ jobs:
arch:
- x64
steps:
- uses: actions/checkout@v2
- uses: actions/checkout@v3
- uses: julia-actions/setup-julia@v1
with:
version: ${{ matrix.version }}
arch: ${{ matrix.arch }}
- uses: actions/cache@v1
- uses: actions/cache@v3
env:
cache-name: cache-artifacts
with:
Expand All @@ -34,6 +34,6 @@ jobs:
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-runtest@v1
- uses: julia-actions/julia-processcoverage@v1
- uses: codecov/codecov-action@v1
- uses: codecov/codecov-action@v3
with:
file: lcov.info
14 changes: 7 additions & 7 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "MomentClosure"
uuid = "01a1b25a-ecf0-48c5-ae58-55bfd5393600"
authors = ["Augustinas Sukys"]
version = "0.2.0"
version = "0.3.0"

[deps]
Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
Expand All @@ -19,17 +19,17 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"

[compat]
Catalyst = "11, 12"
Catalyst = "13"
Combinatorics = "1"
Cumulants = "1"
DataStructures = "0.18"
Distributions = "0.23, 0.24, 0.25"
Distributions = "0.25"
DocStringExtensions = "0.8, 0.9"
Latexify = "0.15.14"
ModelingToolkit = "6.5, 7, 8"
Latexify = "0.15.14, 0.16"
ModelingToolkit = "8"
SciMLBase = "1.8.3"
SymbolicUtils = "0.19"
Symbolics = "3, 4"
SymbolicUtils = "1"
Symbolics = "5"
TupleTools = "1.2.0"
julia = "1.6"

Expand Down
169 changes: 41 additions & 128 deletions docs/src/assets/derivative_matching_cumulant.svg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
19 changes: 11 additions & 8 deletions docs/src/tutorials/LMA_example.md
Original file line number Diff line number Diff line change
Expand Up @@ -20,20 +20,22 @@ using Catalyst
# how the nonlinear reactions are to be transformed using LMA

rn_nonlinear = @reaction_network begin
@parameters σ_b σ_u ρ_b ρ_u
σ_b, g + p 0
σ_u*(1-g), 0 g + p
ρ_u, g g + p
ρ_b*(1-g), 0 p
1, p 0
end σ_b σ_u ρ_b ρ_u
end

rn_linear = @reaction_network begin
@parameters σ_b_LMA σ_u ρ_b ρ_u
σ_b_LMA, g 0 # typing ̄σ_b is not allowed it seems
σ_u*(1-g), 0 g
ρ_u, g g+p
(ρ_b*(1-g)), 0 p
1, p 0
end σ_b_LMA σ_u ρ_b ρ_u
end
```
We can now apply the LMA to find the effective parameter $\bar{σ}_b$ and generate the corresponding moment equations of the linear GRN using MomentClosure's [`linear_mapping_approximation`](@ref):
```julia
Expand All @@ -45,7 +47,7 @@ using MomentClosure
@variables g(t)
binary_vars = [speciesmap(rn_nonlinear)[g]]

LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false)
LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false)
display(effective_params)
```
```julia
Expand Down Expand Up @@ -79,7 +81,7 @@ u₀map = deterministic_IC(u₀, LMA_eqs)
oprob_LMA = ODEProblem(LMA_eqs, u₀map, tspan, p)
sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)

plot(sol_LMA, vars=(0, [2]), label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg")
plot(sol_LMA, idxs=[2], label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg")
```
![LMA feedback loop mean protein number](../assets/LMA_feedback_loop_mean_protein_number.svg)

Expand Down Expand Up @@ -141,11 +143,12 @@ As pointed out by the authors in [1], a more efficient way is to expand the gene
However, TaylorSeries only supports elementary function operations at the time and hence evaluating the Kummer's function $M(\cdot,\cdot,\cdot)$ requires some more work (these specialised numerics are readily available in more established scientific computing frameworks such as Mathematica but there's no fun in that). We can extend the TaylorSeries framework by constructing a function `t_pFq` that implements a recurrence relation between the Taylor coefficients for the generalized hypergeometric function `pFq` as defined in [HypergeometricFunctions.jl](https://github.com/JuliaMath/HypergeometricFunctions.jl). This can be done as follows (note that our construction is valid only for a single-variable Taylor series [`Taylor1`](https://juliadiff.org/TaylorSeries.jl/stable/api/#TaylorSeries.Taylor1)):
```julia
using TaylorSeries, HypergeometricFunctions
using HypergeometricFunctions: pFqweniger

# please let me know if a simpler and more efficient way to do this exists!
function t_pFq::AbstractVector, β::AbstractVector, a::Taylor1)
order = a.order
aux = pFq(α, β, constant_term(a))
aux = pFqweniger(α, β, constant_term(a))
c = Taylor1(aux, order)

iszero(order) && return c
Expand All @@ -170,13 +173,13 @@ oprob_LMA = remake(oprob_LMA, tspan=tspan)
sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)

# rebuild the symbolic expression for the effective parameter as a function of raw moments
μ_sym = LMA_eqs.odes.states
p_sub = Pair.(LMA_eqs.odes.ps, p)
μ_sym = [LMA_eqs.odes.states...]
p_sub = Dict.(Pair.(LMA_eqs.odes.ps, p)...)
avg_σ_b_sym = collect(values(effective_params))[1]
fn = build_function(substitute(avg_σ_b_sym, p_sub), μ_sym)
avg_σ_b = eval(fn)
# evaluate the time-averaged value of the effective parameter
σ_b_avg = sum(avg_σ_b.(sol_LMA[:])) * dt / T
@time σ_b_avg = sum(avg_σ_b.(sol_LMA.u)) * dt / T
```
We proceed with the very last steps of the LMA to obtain the probability distribution:
```julia
Expand Down
31 changes: 16 additions & 15 deletions docs/src/tutorials/P53_system_example.md
Original file line number Diff line number Diff line change
Expand Up @@ -40,21 +40,22 @@ with parameters

We begin by loading all the packages we will need
```julia
using Catalyst, MomentClosure, OrdinaryDiffEq, DiffEqJump,
using Catalyst, MomentClosure, OrdinaryDiffEq, JumpProcesses,
DiffEqBase.EnsembleAnalysis, Plots, Plots.PlotMeasures
```
and then build the model using Catalyst and set its parameters as follows:
```julia
# → for mass-actions rate
# ⇒ for non mass-actions rate
rn = @reaction_network begin
@parameters k₁ k₂ k₃ k₄ k₅ k₆ k₇
(k₁), 0 x
(k₂), x 0
(k₃*x*y/(x+k₇)), x 0
(k₄*x), 0 y₀
(k₅), y₀ y
(k₆), y 0
end k₁ k₂ k₃ k₄ k₅ k₆ k₇
end

# parameters [k₁, k₂, k₃, k₄, k₅, k₆, k₇]
p = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01]
Expand Down Expand Up @@ -106,7 +107,7 @@ h2 = histogram(data[2], normalize=true, xlabel="y₀", ylabel="P(y₀)")
h3 = histogram(data[3], normalize=true, xlabel="y", ylabel="P(y)")
using Plots.PlotMeasures
plot(h1, h2, h3, legend=false, layout=(1,3), size = (1050, 250),
left_margin = 5mm, bottom_margin = 7mm, guidefontsize=10)
left_margin = 5PlotMeasures.mm, bottom_margin = 7PlotMeasures.mm, guidefontsize=10)
```
![P53-Mdm2 distribution](../assets/p53-Mdm2_distribution.svg)

Expand All @@ -120,15 +121,15 @@ closures = ["normal", "log-normal", "gamma"]
plts = [plot() for i in 1:length(closures)]

for q in 3:6
eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false)
eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)
for (closure, plt) in zip(closures, plts)
closed_eqs = moment_closure(eqs, closure)

u₀map = deterministic_IC(u₀, closed_eqs)
oprob = ODEProblem(closed_eqs, u₀map, tspan, p)

sol = solve(oprob, Tsit5(), saveat=0.1)
plt = plot!(plt, sol, vars=(0, 1), lw=3, label = "q = "*string(q))
plt = plot!(plt, sol, idxs=[1], lw=3, label = "q = "*string(q))
end
end

Expand All @@ -139,13 +140,13 @@ end
```
Normal closure:
```julia
plot(plts[1], size=(750, 450), leftmargin=2mm)
plot(plts[1], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
![P53-Mdm2 normal means 2nd order expansion](../assets/p53-Mdm2_normal_2nd_order.svg)

Log-normal closure:
```julia
plot(plts[2], size=(750, 450), leftmargin=2mm)
plot(plts[2], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
![P53-Mdm2 log-normal means 2nd order expansion](../assets/p53-Mdm2_log-normal_2nd_order.svg)

Expand All @@ -157,7 +158,7 @@ plot(plts[2], xlims=(0., 50.), lw=3)

Gamma closure:
```julia
plot(plts[3], size=(750, 450), leftmargin=2mm)
plot(plts[3], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
![P53-Mdm2 gamma means 2nd order expansion](../assets/p53-Mdm2_gamma_2nd_order.svg)

Expand All @@ -168,7 +169,7 @@ We can also plot the variance predictions:
# rerunning the same calculations as they are reasonably fast
plt = plot()
for q in [4,6]
eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false)
eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)
for closure in closures
closed_eqs = moment_closure(eqs, closure)

Expand All @@ -177,7 +178,7 @@ for q in [4,6]
sol = solve(oprob, Tsit5(), saveat=0.1)

# index of M₂₀₀ can be checked with `u₀map` or `closed_eqs.odes.states`
plt = plot!(plt, sol, vars=(0, 4), lw=3, label = closure*" q = "*string(q))
plt = plot!(plt, sol, idxs=[4], lw=3, label = closure*" q = "*string(q))
end
end

Expand All @@ -200,7 +201,7 @@ q_vals = [4, 6]

for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars)

eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaw=false)
eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaws=false)
for closure in closures

closed_eqs = moment_closure(eqs, closure)
Expand All @@ -209,8 +210,8 @@ for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars)
oprob = ODEProblem(closed_eqs, u₀map, tspan, p)

sol = solve(oprob, Tsit5(), saveat=0.1)
plt_m = plot!(plt_m, sol, vars=(0, 1), label = closure)
plt_v = plot!(plt_v, sol, vars=(0, 4), label = closure)
plt_m = plot!(plt_m, sol, idxs=[1], label = closure)
plt_v = plot!(plt_v, sol, idxs=[4], label = closure)

end

Expand Down Expand Up @@ -270,7 +271,7 @@ Finally, we can check whether better estimates can be obtained using even higher
plt = plot()
closures = ["zero", "normal", "log-normal", "gamma"]

eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaw=false)
eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaws=false)
# faster to store than recompute in case we want to try different solvers/params
oprobs = Dict()

Expand All @@ -281,7 +282,7 @@ for closure in closures
oprobs[closure] = ODEProblem(closed_eqs, u₀map, tspan, p)
sol = solve(oprobs[closure], Tsit5(), saveat=0.1)

plt = plot!(plt, sol, vars=(0, 1), label = closure)
plt = plot!(plt, sol, idxs=[1], label = closure)
end

plt = plot!(plt, xlabel = "Time [h]", ylabel = "Mean p53 molecule number")
Expand Down
Loading

0 comments on commit e35caf7

Please sign in to comment.