{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `/opt/julia/registries/General`\n",
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PolyesterWeave ─────────────────── v0.1.7\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LayoutPointers ─────────────────── v0.1.8\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ConstructionBase ───────────────── v1.4.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m VectorizationBase ──────────────── v0.21.33\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m BandedMatrices ─────────────────── v0.16.13\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DifferentialEquations ──────────── v7.1.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DEDataArrays ───────────────────── v0.2.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StochasticDiffEq ───────────────── v6.49.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NLSolversBase ──────────────────── v7.8.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m libblastrampoline_jll ──────────── v3.1.0+2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FunctionWrappers ───────────────── v1.1.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LevyArea ───────────────────────── v1.0.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PoissonRandom ──────────────────── v0.4.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m CloseOpenIntervals ─────────────── v0.1.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DiffEqNoiseProcess ─────────────── v5.12.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m OrdinaryDiffEq ─────────────────── v6.11.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m CPUSummary ─────────────────────── v0.1.25\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Graphs ─────────────────────────── v1.5.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ThreadingUtilities ─────────────── v0.5.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FastClosures ───────────────────── v0.3.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArrayInterfaceStaticArraysCore ─── v0.1.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m VertexSafeGraphs ───────────────── v0.2.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m BoundaryValueDiffEq ────────────── v2.8.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m UnPack ─────────────────────────── v1.0.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SparseDiffTools ────────────────── v1.21.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m CommonSolve ────────────────────── v0.2.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m GPUArraysCore ──────────────────── v0.1.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SteadyStateDiffEq ──────────────── v1.9.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m HostCPUFeatures ────────────────── v0.1.8\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Polyester ──────────────────────── v0.6.12\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SLEEFPirates ───────────────────── v0.6.33\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Parameters ─────────────────────── v0.12.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NLsolve ────────────────────────── v4.5.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Optim ──────────────────────────── v1.7.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ResettableStacks ───────────────── v1.1.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m MuladdMacro ────────────────────── v0.2.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SimpleTraits ───────────────────── v0.9.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FiniteDiff ─────────────────────── v2.13.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m KLU ────────────────────────────── v0.3.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DiffEqBase ─────────────────────── v6.84.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ExponentialUtilities ───────────── v1.18.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m IterativeSolvers ───────────────── v0.9.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArrayInterfaceCore ─────────────── v0.1.16\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DiffEqCallbacks ────────────────── v2.23.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m RecursiveFactorization ─────────── v0.2.11\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LabelledArrays ─────────────────── v1.9.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ManualMemory ───────────────────── v0.1.8\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Setfield ───────────────────────── v0.8.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m KrylovKit ──────────────────────── v0.5.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m RandomNumbers ──────────────────── v1.5.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Inflate ────────────────────────── v0.1.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m BitTwiddlingConvenienceFunctions ─ v0.1.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m GenericSchur ───────────────────── v0.5.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DiffEqJump ─────────────────────── v8.3.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArnoldiMethod ──────────────────── v0.2.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DelayDiffEq ────────────────────── v5.37.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SciMLBase ──────────────────────── v1.45.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LoopVectorization ──────────────── v0.12.110\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m TreeViews ──────────────────────── v0.3.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StaticArraysCore ───────────────── v1.0.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m IfElse ─────────────────────────── v0.1.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArrayLayouts ───────────────────── v0.7.5\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArrayInterface ─────────────────── v3.1.32\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Sundials ───────────────────────── v4.9.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Krylov ─────────────────────────── v0.8.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LinearSolve ────────────────────── v1.20.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ZygoteRules ────────────────────── v0.2.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SIMDDualNumbers ────────────────── v0.1.1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Static ─────────────────────────── v0.3.3\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Random123 ──────────────────────── v1.6.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m TriangularSolve ────────────────── v0.1.12\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PreallocationTools ─────────────── v0.3.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StrideArraysCore ───────────────── v0.3.15\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m RecursiveArrayTools ────────────── v2.31.2\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NonlinearSolve ─────────────────── v0.3.22\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PositiveFactorizations ─────────── v0.2.4\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FastBroadcast ──────────────────── v0.1.17\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m CpuId ──────────────────────────── v0.3.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Sundials_jll ───────────────────── v5.2.0+1\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SIMDTypes ──────────────────────── v0.1.0\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m LineSearches ───────────────────── v7.1.1\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Project.toml`\n",
" \u001b[90m [0c46a032] \u001b[39m\u001b[92m+ DifferentialEquations v7.1.0\u001b[39m\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Manifest.toml`\n",
" \u001b[90m [ec485272] \u001b[39m\u001b[92m+ ArnoldiMethod v0.2.0\u001b[39m\n",
" \u001b[90m [4fba245c] \u001b[39m\u001b[92m+ ArrayInterface v3.1.32\u001b[39m\n",
" \u001b[90m [30b0a656] \u001b[39m\u001b[92m+ ArrayInterfaceCore v0.1.16\u001b[39m\n",
" \u001b[90m [dd5226c6] \u001b[39m\u001b[92m+ ArrayInterfaceStaticArraysCore v0.1.0\u001b[39m\n",
" \u001b[90m [4c555306] \u001b[39m\u001b[92m+ ArrayLayouts v0.7.5\u001b[39m\n",
" \u001b[90m [aae01518] \u001b[39m\u001b[92m+ BandedMatrices v0.16.13\u001b[39m\n",
" \u001b[90m [62783981] \u001b[39m\u001b[92m+ BitTwiddlingConvenienceFunctions v0.1.4\u001b[39m\n",
" \u001b[90m [764a87c0] \u001b[39m\u001b[92m+ BoundaryValueDiffEq v2.8.0\u001b[39m\n",
" \u001b[90m [2a0fbf3d] \u001b[39m\u001b[92m+ CPUSummary v0.1.25\u001b[39m\n",
" \u001b[90m [fb6a15b2] \u001b[39m\u001b[92m+ CloseOpenIntervals v0.1.3\u001b[39m\n",
" \u001b[90m [38540f10] \u001b[39m\u001b[92m+ CommonSolve v0.2.1\u001b[39m\n",
" \u001b[90m [187b0558] \u001b[39m\u001b[92m+ ConstructionBase v1.4.0\u001b[39m\n",
" \u001b[90m [adafc99b] \u001b[39m\u001b[92m+ CpuId v0.3.0\u001b[39m\n",
" \u001b[90m [754358af] \u001b[39m\u001b[92m+ DEDataArrays v0.2.4\u001b[39m\n",
" \u001b[90m [bcd4f6db] \u001b[39m\u001b[92m+ DelayDiffEq v5.37.0\u001b[39m\n",
" \u001b[90m [2b5f629d] \u001b[39m\u001b[92m+ DiffEqBase v6.84.0\u001b[39m\n",
" \u001b[90m [459566f4] \u001b[39m\u001b[92m+ DiffEqCallbacks v2.23.1\u001b[39m\n",
" \u001b[90m [c894b116] \u001b[39m\u001b[92m+ DiffEqJump v8.3.0\u001b[39m\n",
" \u001b[90m [77a26b50] \u001b[39m\u001b[92m+ DiffEqNoiseProcess v5.12.0\u001b[39m\n",
" \u001b[90m [0c46a032] \u001b[39m\u001b[92m+ DifferentialEquations v7.1.0\u001b[39m\n",
" \u001b[90m [d4d017d3] \u001b[39m\u001b[92m+ ExponentialUtilities v1.18.0\u001b[39m\n",
" \u001b[90m [7034ab61] \u001b[39m\u001b[92m+ FastBroadcast v0.1.17\u001b[39m\n",
" \u001b[90m [9aa1b823] \u001b[39m\u001b[92m+ FastClosures v0.3.2\u001b[39m\n",
" \u001b[90m [6a86dc24] \u001b[39m\u001b[92m+ FiniteDiff v2.13.1\u001b[39m\n",
" \u001b[90m [069b7b12] \u001b[39m\u001b[92m+ FunctionWrappers v1.1.2\u001b[39m\n",
" \u001b[90m [46192b85] \u001b[39m\u001b[92m+ GPUArraysCore v0.1.1\u001b[39m\n",
" \u001b[90m [c145ed77] \u001b[39m\u001b[92m+ GenericSchur v0.5.3\u001b[39m\n",
" \u001b[90m [86223c79] \u001b[39m\u001b[92m+ Graphs v1.5.0\u001b[39m\n",
" \u001b[90m [3e5b6fbb] \u001b[39m\u001b[92m+ HostCPUFeatures v0.1.8\u001b[39m\n",
" \u001b[90m [615f187c] \u001b[39m\u001b[92m+ IfElse v0.1.1\u001b[39m\n",
" \u001b[90m [d25df0c9] \u001b[39m\u001b[92m+ Inflate v0.1.2\u001b[39m\n",
" \u001b[90m [42fd0dbc] \u001b[39m\u001b[92m+ IterativeSolvers v0.9.2\u001b[39m\n",
" \u001b[90m [ef3ab10e] \u001b[39m\u001b[92m+ KLU v0.3.0\u001b[39m\n",
" \u001b[90m [ba0b0d4f] \u001b[39m\u001b[92m+ Krylov v0.8.3\u001b[39m\n",
" \u001b[90m [0b1a1467] \u001b[39m\u001b[92m+ KrylovKit v0.5.4\u001b[39m\n",
" \u001b[90m [2ee39098] \u001b[39m\u001b[92m+ LabelledArrays v1.9.0\u001b[39m\n",
" \u001b[90m [10f19ff3] \u001b[39m\u001b[92m+ LayoutPointers v0.1.8\u001b[39m\n",
" \u001b[90m [2d8b4e74] \u001b[39m\u001b[92m+ LevyArea v1.0.0\u001b[39m\n",
" \u001b[90m [d3d80556] \u001b[39m\u001b[92m+ LineSearches v7.1.1\u001b[39m\n",
" \u001b[90m [7ed4a6bd] \u001b[39m\u001b[92m+ LinearSolve v1.20.0\u001b[39m\n",
" \u001b[90m [bdcacae8] \u001b[39m\u001b[92m+ LoopVectorization v0.12.110\u001b[39m\n",
" \u001b[90m [d125e4d3] \u001b[39m\u001b[92m+ ManualMemory v0.1.8\u001b[39m\n",
" \u001b[90m [46d2c3a1] \u001b[39m\u001b[92m+ MuladdMacro v0.2.2\u001b[39m\n",
" \u001b[90m [d41bc354] \u001b[39m\u001b[92m+ NLSolversBase v7.8.2\u001b[39m\n",
" \u001b[90m [2774e3e8] \u001b[39m\u001b[92m+ NLsolve v4.5.1\u001b[39m\n",
" \u001b[90m [8913a72c] \u001b[39m\u001b[92m+ NonlinearSolve v0.3.22\u001b[39m\n",
" \u001b[90m [429524aa] \u001b[39m\u001b[92m+ Optim v1.7.1\u001b[39m\n",
" \u001b[90m [1dea7af3] \u001b[39m\u001b[92m+ OrdinaryDiffEq v6.11.2\u001b[39m\n",
" \u001b[90m [d96e819e] \u001b[39m\u001b[92m+ Parameters v0.12.3\u001b[39m\n",
" \u001b[90m [e409e4f3] \u001b[39m\u001b[92m+ PoissonRandom v0.4.1\u001b[39m\n",
" \u001b[90m [f517fe37] \u001b[39m\u001b[92m+ Polyester v0.6.12\u001b[39m\n",
" \u001b[90m [1d0040c9] \u001b[39m\u001b[92m+ PolyesterWeave v0.1.7\u001b[39m\n",
" \u001b[90m [85a6dd25] \u001b[39m\u001b[92m+ PositiveFactorizations v0.2.4\u001b[39m\n",
" \u001b[90m [d236fae5] \u001b[39m\u001b[92m+ PreallocationTools v0.3.2\u001b[39m\n",
" \u001b[90m [74087812] \u001b[39m\u001b[92m+ Random123 v1.6.0\u001b[39m\n",
" \u001b[90m [e6cf234a] \u001b[39m\u001b[92m+ RandomNumbers v1.5.3\u001b[39m\n",
" \u001b[90m [731186ca] \u001b[39m\u001b[92m+ RecursiveArrayTools v2.31.2\u001b[39m\n",
" \u001b[90m [f2c3362d] \u001b[39m\u001b[92m+ RecursiveFactorization v0.2.11\u001b[39m\n",
" \u001b[90m [ae5879a3] \u001b[39m\u001b[92m+ ResettableStacks v1.1.1\u001b[39m\n",
" \u001b[90m [3cdde19b] \u001b[39m\u001b[92m+ SIMDDualNumbers v0.1.1\u001b[39m\n",
" \u001b[90m [94e857df] \u001b[39m\u001b[92m+ SIMDTypes v0.1.0\u001b[39m\n",
" \u001b[90m [476501e8] \u001b[39m\u001b[92m+ SLEEFPirates v0.6.33\u001b[39m\n",
" \u001b[90m [0bca4576] \u001b[39m\u001b[92m+ SciMLBase v1.45.0\u001b[39m\n",
" \u001b[90m [efcf1570] \u001b[39m\u001b[92m+ Setfield v0.8.2\u001b[39m\n",
" \u001b[90m [699a6c99] \u001b[39m\u001b[92m+ SimpleTraits v0.9.4\u001b[39m\n",
" \u001b[90m [47a9eef4] \u001b[39m\u001b[92m+ SparseDiffTools v1.21.0\u001b[39m\n",
" \u001b[90m [aedffcd0] \u001b[39m\u001b[92m+ Static v0.3.3\u001b[39m\n",
" \u001b[90m [1e83bf80] \u001b[39m\u001b[92m+ StaticArraysCore v1.0.1\u001b[39m\n",
" \u001b[90m [9672c7b4] \u001b[39m\u001b[92m+ SteadyStateDiffEq v1.9.0\u001b[39m\n",
" \u001b[90m [789caeaf] \u001b[39m\u001b[92m+ StochasticDiffEq v6.49.1\u001b[39m\n",
" \u001b[90m [7792a7ef] \u001b[39m\u001b[92m+ StrideArraysCore v0.3.15\u001b[39m\n",
" \u001b[90m [c3572dad] \u001b[39m\u001b[92m+ Sundials v4.9.4\u001b[39m\n",
" \u001b[90m [8290d209] \u001b[39m\u001b[92m+ ThreadingUtilities v0.5.0\u001b[39m\n",
" \u001b[90m [a2a6695c] \u001b[39m\u001b[92m+ TreeViews v0.3.0\u001b[39m\n",
" \u001b[90m [d5829a12] \u001b[39m\u001b[92m+ TriangularSolve v0.1.12\u001b[39m\n",
" \u001b[90m [3a884ed6] \u001b[39m\u001b[92m+ UnPack v1.0.2\u001b[39m\n",
" \u001b[90m [3d5dd08c] \u001b[39m\u001b[92m+ VectorizationBase v0.21.33\u001b[39m\n",
" \u001b[90m [19fa3120] \u001b[39m\u001b[92m+ VertexSafeGraphs v0.2.0\u001b[39m\n",
" \u001b[90m [700de1a5] \u001b[39m\u001b[92m+ ZygoteRules v0.2.2\u001b[39m\n",
" \u001b[90m [fb77eaff] \u001b[39m\u001b[92m+ Sundials_jll v5.2.0+1\u001b[39m\n",
" \u001b[90m [8e850b90] \u001b[39m\u001b[92m+ libblastrampoline_jll v3.1.0+2\u001b[39m\n",
" \u001b[90m [4536629a] \u001b[39m\u001b[92m+ OpenBLAS_jll\u001b[39m\n",
" \u001b[90m [bea87d4a] \u001b[39m\u001b[92m+ SuiteSparse_jll\u001b[39m\n",
"\u001b[32m\u001b[1mPrecompiling\u001b[22m\u001b[39m project...\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mCommonSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSIMDTypes\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mIfElse\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mMuladdMacro\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStaticArraysCore\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mConstructionBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mUnPack\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mFastClosures\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mTreeViews\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPositiveFactorizations\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mArrayInterfaceCore\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mCpuId\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mManualMemory\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPoissonRandom\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mFunctionWrappers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mGPUArraysCore\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mInflate\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOpenBLAS_jll\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mParameters\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mZygoteRules\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mGenericSchur\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStatic\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mKrylovKit\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mRandomNumbers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mArrayInterfaceStaticArraysCore\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLevyArea\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mlibblastrampoline_jll\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSuiteSparse_jll\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSimpleTraits\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mBitTwiddlingConvenienceFunctions\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mResettableStacks\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mThreadingUtilities\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mIterativeSolvers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mCPUSummary\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mHostCPUFeatures\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSetfield\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mKLU\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mRandom123\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mArnoldiMethod\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPolyesterWeave\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mArrayInterface\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSundials_jll\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mExponentialUtilities\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mFiniteDiff\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mKrylov\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mCloseOpenIntervals\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLayoutPointers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mArrayLayouts\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mRecursiveArrayTools\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStrideArraysCore\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLabelledArrays\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNLSolversBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPolyester\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mFastBroadcast\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mBandedMatrices\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPreallocationTools\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mGraphs\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLineSearches\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mVertexSafeGraphs\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSciMLBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNLsolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSparseDiffTools\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDEDataArrays\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOptim\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mVectorizationBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSLEEFPirates\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSIMDDualNumbers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLoopVectorization\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mTriangularSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mRecursiveFactorization\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNonlinearSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLinearSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqCallbacks\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqNoiseProcess\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqJump\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mBoundaryValueDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSundials\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSteadyStateDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOrdinaryDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDelayDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStochasticDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39mDifferentialEquations\n",
" 83 dependencies successfully precompiled in 218 seconds (214 already precompiled)\n",
"\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
"\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m SpecialFunctions ─ v1.8.3\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Project.toml`\n",
" \u001b[90m [276daf66] \u001b[39m\u001b[92m+ SpecialFunctions v1.8.3\u001b[39m\n",
"\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m `/opt/julia/environments/v1.6/Manifest.toml`\n",
" \u001b[90m [276daf66] \u001b[39m\u001b[93m↑ SpecialFunctions v1.6.1 ⇒ v1.8.3\u001b[39m\n",
" \u001b[90m [05823500] \u001b[39m\u001b[92m+ OpenLibm_jll\u001b[39m\n",
"\u001b[32m\u001b[1mPrecompiling\u001b[22m\u001b[39m project...\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOpenLibm_jll\u001b[39m\n",
"\u001b[32m ✓ \u001b[39mSpecialFunctions\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffRules\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLevyArea\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStatsFuns\u001b[39m\n",
"\u001b[32m ✓ \u001b[39mForwardDiff\n",
"\u001b[32m ✓ \u001b[39mDistributions\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mPreallocationTools\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNLSolversBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSparseDiffTools\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSIMDDualNumbers\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mKernelDensity\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLineSearches\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNLsolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOptim\u001b[39m\n",
"\u001b[32m ✓ \u001b[39mJuMP\n",
"\u001b[32m ✓ \u001b[39mGadfly\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLoopVectorization\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mTriangularSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mRecursiveFactorization\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mNonlinearSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mLinearSolve\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqBase\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqCallbacks\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqJump\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDiffEqNoiseProcess\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mBoundaryValueDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSundials\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mSteadyStateDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mOrdinaryDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mDelayDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39m\u001b[90mStochasticDiffEq\u001b[39m\n",
"\u001b[32m ✓ \u001b[39mDifferentialEquations\n",
" 33 dependencies successfully precompiled in 202 seconds (265 already precompiled)\n"
]
}
],
"source": [
"using Pkg\n",
"Pkg.add(\"DifferentialEquations\");\n",
"Pkg.add(\"SpecialFunctions\");\n",
"#Pkg.add(\"StatsBase\");"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"using DifferentialEquations, Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Ordinary Differential Equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ODE 1\n",
"\n",
"$$\\frac{du(t)}{dt} = - \\kappa u(t) \\quad u(t_0) = 0.1, $$ \n",
"where $u: \\mathbb{R}^+\\to\\mathbb{R}^+$ and $\\kappa = 0.25$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define your differential equation as a function that updates the LHS"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ode1! (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ode1!(du, u, p, t)\n",
" κ = p\n",
" du .= -κ.*u\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set the initial condition, the time interval of interest and the parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"u₀ = [0.1]\n",
"\n",
"tspan = (0, 10)\n",
"\n",
"κ = 0.25\n",
"p = [κ];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create an ordinary differential equation problem with the variables defined above."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[36mODEProblem\u001b[0m with uType \u001b[36mVector{Float64}\u001b[0m and tType \u001b[36mInt64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n",
"timespan: (0, 10)\n",
"u0: 1-element Vector{Float64}:\n",
" 0.1"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prob = ODEProblem(ode1!, u₀, tspan, p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve the ODE problem"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"sol = solve(prob)\n",
"# OR sol = solve(####, #### (algorithm name()); #### (parameters))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot your numerical solution "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol, label = \"Numeric solution\", linewidth= 6.0) \n",
"plot!(t -> u₀[1]*exp(-κ*t), label = \"Analytical solution\", linewidth= 6.0, linestyle=:dash)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### ODE 2\n",
"\n",
"$$d\\begin{bmatrix}u_1(t) \\\\ u_2(t) \\end{bmatrix} = \\begin{bmatrix}-\\kappa_1u_1(t)+u_2(t) \\\\ -\\kappa_2u_1(t)u_2(t) \\end{bmatrix}dt, \\quad\\quad \\begin{bmatrix}u_1(0) \\\\ u_2(0) \\end{bmatrix} = \\begin{bmatrix}1 \\\\ 2 \\end{bmatrix}.$$\n",
"\n",
"Here, $u = (u_1, u_2)^\\top$ and $\\kappa_1 = 0.5$ and $\\kappa_2 = 0.25$."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ode2! (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ode2!(du, u, p, t)\n",
" κ₁, κ₂ = p\n",
" du .= [-κ₁*u[1]+u[2], \n",
" -κ₂*u[1]*u[2]]\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Set the initial condition, the time interval of interest and the parameters"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"u₀ = [1, 2]\n",
"\n",
"tspan = (0, 10)\n",
"\n",
"κ₁ = 0.5\n",
"κ₂ = 0.25\n",
"p = [κ₁, κ₂];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create an ordinary differential equation problem with the variables defined above."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[36mODEProblem\u001b[0m with uType \u001b[36mVector{Int64}\u001b[0m and tType \u001b[36mInt64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n",
"timespan: (0, 10)\n",
"u0: 2-element Vector{Int64}:\n",
" 1\n",
" 2"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prob = ODEProblem(ode2!, u₀, tspan, p)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve the ODE problem"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"sol = solve(prob)\n",
"# OR sol2 = solve(prob2, #### (algorithm name); #### (parameters))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol, label = [\"u1(t)\" \"u2(t)\"], linewidth= 6.0) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The object \"sol\" behaves like a function; for an arbitrary t, you can ask for sol(t). "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" 0.35486822476777363\n",
" 0.11976761491105567"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol(10.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2.Degradation and production (revisited)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) numerical solution to the chemical master equation "
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ME_degr_and_prod! (generic function with 1 method)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ME_degr_and_prod!(du, u, p, t)\n",
" \n",
" κ₁, κ₂, V = p\n",
"\n",
" u ./= sum(u) \n",
"\n",
" u_padded = [0; u; 0]\n",
" \n",
" n = 0:length(u)-1\n",
" \n",
" a = -κ₁.*n .- κ₂*V\n",
" b = (n.+1)*κ₁\n",
" c = κ₂*V\n",
"\n",
" du .= a.*u + b.*u_padded[3:end] + c*u_padded[1:end-2]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"n_upperbound=24\n",
"u₀ = [1; zeros(n_upperbound)]\n",
"\n",
"tspan = (0, 50)\n",
"\n",
"κ₁, κ₂, V = [0.1, 1000, 0.001]\n",
"p = [κ₁, κ₂, V];"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[36mODEProblem\u001b[0m with uType \u001b[36mVector{Float64}\u001b[0m and tType \u001b[36mInt64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n",
"timespan: (0, 50)\n",
"u0: 25-element Vector{Float64}:\n",
" 1.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prob = ODEProblem(ME_degr_and_prod!, u₀, tspan, p)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"retcode: Success\n",
"Interpolation: specialized 4th order \"free\" interpolation, specialized 2nd order \"free\" stiffness-aware interpolation\n",
"t: 57-element Vector{Float64}:\n",
" 0.0\n",
" 0.0009990005004983772\n",
" 0.010989005505482149\n",
" 0.04857417416890508\n",
" 0.10652138158604843\n",
" 0.1759574101735114\n",
" 0.2669280280363636\n",
" 0.37864878022829973\n",
" 0.5191537297722533\n",
" 0.6955260633200333\n",
" 0.9153390708720146\n",
" 1.183078908419348\n",
" 1.5090413790635058\n",
" ⋮\n",
" 24.138912733435\n",
" 24.73573216325989\n",
" 25.33255295277502\n",
" 26.526196697827384\n",
" 28.072969468323404\n",
" 30.526332052963134\n",
" 33.00377825555365\n",
" 36.320417212384775\n",
" 39.83393724905629\n",
" 44.236375992799225\n",
" 49.09303441904753\n",
" 50.0\n",
"u: 57-element Vector{Vector{Float64}}:\n",
" [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
" [0.9990015481829916, 0.000997953198016281, 4.98452973993036e-7, 1.659766328213407e-10, 4.14506153330629e-14, 8.281093359105346e-18, 1.4236437908660304e-21, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
" [0.9890771228087425, 0.01086300417409281, 5.965402341651524e-5, 2.1839277749333684e-7, 5.996514634314703e-10, 1.3167984317841413e-12, 2.4714878093675965e-15, 1.9029220878165408e-18, 8.633990916434816e-22, 2.736886433766084e-25 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
" [0.9526988832825186, 0.046164350954624296, 0.0011184789586409364, 1.806581200547898e-5, 2.1885531046300414e-7, 2.1193831312241796e-9, 1.741455337882848e-11, 1.0218925116220429e-13, 4.287762445173533e-16, 1.3570997255086707e-18 … 1.1438512391522709e-36, 2.7532872476228315e-40, 5.342038189812808e-44, 8.182929338874347e-48, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]\n",
" [0.8994641782529833, 0.09530367181754885, 0.005049000414446518, 0.00017832394900178184, 4.72369121427449e-6, 1.0007610456174494e-7, 1.7717595241517965e-9, 2.659919631212679e-11, 3.384561283491325e-13, 3.645800033232894e-15 … 9.276946171647651e-29, 2.2130718539590556e-31, 4.296862957673146e-34, 6.438038755481103e-37, 4.619036082306904e-40, 2.0312598638165287e-43, 6.329101782177271e-47, 1.51726219368882e-50, 2.933171860104335e-54, 4.437220109562198e-58]\n",
" [0.8399453913713946, 0.14650193074198922, 0.012776316447958631, 0.0007428078753593389, 3.239006666006225e-5, 1.1298040324779254e-6, 3.2855965854291456e-8, 8.185303952421183e-10, 1.776440400598095e-11, 3.391306301694056e-13 … 1.2710681431759963e-24, 1.004032534527242e-26, 6.807512174532853e-29, 3.875119813586576e-31, 1.753547154375641e-33, 6.226273359737543e-36, 1.7632330301170553e-38, 4.0612149065354745e-41, 7.611164244544021e-44, 1.0831320395856838e-46]\n",
" [0.7684367774735404, 0.20240393122180286, 0.02665629421340624, 0.0023403942750446325, 0.00015411421075179065, 8.118301952384433e-6, 3.564349365224978e-7, 1.3414137796310352e-8, 4.41228583639819e-10, 1.2855980451229141e-11 … 8.815645636897414e-22, 1.2670194771739073e-23, 1.6492643111001122e-25, 1.9383138445637152e-27, 2.0465372366545033e-29, 1.9280493932814874e-31, 1.6076260869122005e-33, 1.1744856497993753e-35, 7.408457154675495e-38, 3.939048010072965e-40]\n",
" [0.689650959051583, 0.25625337452845315, 0.04760799526139763, 0.005896556963279582, 0.000547748692965707, 4.070435409347011e-5, 2.520840252777415e-6, 1.3382692064740629e-7, 6.214989355387301e-9, 2.5625378081542633e-10 … 1.6215940346897408e-19, 3.511948347038737e-21, 7.011440177415021e-23, 1.2906687087062292e-24, 2.189784126348927e-26, 3.420632794122616e-28, 4.9108356992160734e-30, 6.463189865711893e-32, 7.77182165019118e-34, 8.499572348364683e-36]\n",
" [0.6029578975750265, 0.3050411607377206, 0.07716138508769596, 0.013012166443283488, 0.0016457501962976785, 0.00016651674035091704, 1.4040421075910448e-5, 1.0148221213855532e-6, 6.41793909686899e-8, 3.6061987606657374e-9 … 1.550453855416268e-17, 4.703795667947446e-19, 1.3255996137151687e-20, 3.474972922653493e-22, 8.481423397362861e-24, 1.9283592490967279e-25, 4.0847100532915764e-27, 8.058863460309328e-29, 1.4799475780013198e-30, 2.525689042396936e-32]\n",
" [0.5107426754063648, 0.343162558457868, 0.11528370182987882, 0.025819237427733642, 0.004336950687287632, 0.0005827833664879016, 6.526032424552504e-5, 6.264256019725602e-6, 5.26159163081341e-7, 3.927755016933663e-8 … 9.567454284287619e-16, 3.911930196851786e-17, 1.491745813723129e-18, 5.3156075711190263e-20, 1.7726839054164284e-21, 5.53913861811611e-23, 1.6231629863437256e-24, 4.4631710720715697e-26, 1.1519214914606803e-27, 2.787412262440226e-29]\n",
" [0.41698867565712666, 0.36473833595491356, 0.1595178178255132, 0.04650969554655074, 0.0101705502630515, 0.0017792177546794485, 0.00025937469011710077, 3.2411127130819535e-5, 3.544017374631234e-6, 3.4446021238902135e-7 … 4.162296003136637e-14, 2.2357250196192727e-15, 1.1227970413565487e-16, 5.284014020854873e-18, 2.334610964478046e-19, 9.698711998831352e-21, 3.7931340686938376e-22, 1.3979544768435237e-23, 4.858754175497383e-25, 1.5898338286808975e-26]\n",
" [0.32766120366030604, 0.36559599734774856, 0.20396198313047423, 0.075858136242124, 0.021160326863598162, 0.004722036348062286, 0.0008781071144324539, 0.00013996631814661143, 1.9522241336375655e-5, 2.42048797080862e-6 … 1.2762104926042221e-12, 8.802031919748255e-14, 5.686937720948686e-15, 3.4505600855883125e-16, 1.970248978729974e-17, 1.0605474148523195e-18, 5.3895001348684834e-20, 2.5888156206301745e-21, 1.1765576576999068e-22, 5.043443625944879e-24]\n",
" [0.2464248455259147, 0.34516628176066305, 0.24173815814232222, 0.11286651584384173, 0.039523240589176924, 0.011072100426409905, 0.0025847499247592777, 0.0005171987182266043, 9.05559801721761e-5, 1.4094451531274138e-5 … 2.93078791071675e-11, 2.548718107964117e-12, 2.079315399230564e-13, 1.595677117344258e-14, 1.1544703893940108e-15, 7.889893371516392e-17, 5.1017655902629246e-18, 3.125595300496123e-19, 1.8163445437750293e-20, 9.956901391932576e-22]\n",
" ⋮\n",
" [0.0001110777270874526, 0.0010113998334755525, 0.0046045671738486746, 0.01397537650137537, 0.031812624724360236, 0.05793292882134088, 0.08791646590966831, 0.11435842241578591, 0.13015897064035484, 0.13168231271482891 … 0.020822320526807454, 0.011848822453616077, 0.006344048024315218, 0.0032081723510380207, 0.0015338009833826436, 0.0006969494224590138, 0.0002976301620954064, 0.00012078513487631653, 4.316926828535015e-5, 1.2528133615487988e-5]\n",
" [0.00010546766710534948, 0.0009657851526510975, 0.004421928766589505, 0.013497449812558274, 0.0308996085015002, 0.056590590390652014, 0.08636827024387497, 0.11298412019498348, 0.1293268196314751, 0.13158523285726267 … 0.02152755039624733, 0.012319681911303554, 0.006633533149684619, 0.003373426532956713, 0.0016218455654494065, 0.0007409015819250732, 0.000318116795161608, 0.00012964389430803608, 4.656168076381886e-5, 1.352736156691957e-5]\n",
" [0.00010044217466731296, 0.000924670544810527, 0.004256258197029929, 0.013061037118631174, 0.030059976276295076, 0.05534642719122799, 0.08491986430394272, 0.11168174025284507, 0.1285177285945382, 0.13145927930570828 … 0.022203964436233525, 0.012774335188145344, 0.006914830556172608, 0.003534982807261165, 0.0017084194234339594, 0.0007843579381462943, 0.00033847477644934006, 0.0001384866564832959, 4.99598762061808e-5, 1.4530680673482391e-5]\n",
" [9.184299779455575e-5, 0.0008537432315464183, 0.003968043979196659, 0.012295116915292557, 0.028572511185418217, 0.053119322480199425, 0.0822950759823298, 0.10928156312639149, 0.12697743418260515, 0.13114514056596352 … 0.023471715299131318, 0.013633649092993883, 0.007451972295115505, 0.0038447033775700344, 0.0018773675071522406, 0.0008679168841873169, 0.00037970247467498, 0.00015502425300904546, 5.729642946814793e-5, 1.629572967794528e-5]\n",
" [8.300218892228221e-5, 0.000779982113485516, 0.0036647567069569924, 0.011479171999210561, 0.02696705461180441, 0.050680733011364484, 0.07937205038990049, 0.10654724751075643, 0.12514736896296308, 0.13066085576364803 … 0.02495031677279393, 0.014649273538069768, 0.008093112754142897, 0.004220370181627749, 0.0020821472537543434, 0.0009728421072503188, 0.00042933936087840414, 0.00017714360016716192, 6.572893824977244e-5, 1.8873338995091638e-5]\n",
" [7.274294293220214e-5, 0.0006932096665325266, 0.0033029157744570207, 0.01049134106727797, 0.024992963904163955, 0.04763068037669603, 0.07564271698042833, 0.1029659089611688, 0.12263679854872732, 0.12983405327936132 … 0.026944963145725816, 0.016040048783113974, 0.00898389828968203, 0.004748949875492894, 0.002374521786847042, 0.0011238288197750023, 0.0005021530938294264, 0.00020940092347639143, 7.84373957778738e-5, 2.2636826775144225e-5]\n",
" [6.557334595183026e-5, 0.0006316955331787115, 0.003042614286939104, 0.009769827870652899, 0.023527747507612346, 0.04532691164203148, 0.07276845865190253, 0.10013272752541305, 0.12056175798520863, 0.12902770567103564 … 0.028572741357280607, 0.01719276667938004, 0.009732894711395552, 0.0051995246236521695, 0.002626884781929832, 0.0012557936736619217, 0.0005664047583788163, 0.0002381951146600946, 8.982242513296768e-5, 2.604203193912821e-5]\n",
" [5.908097972297847e-5, 0.0005753209310242901, 0.0028010951252595818, 0.009091739061313718, 0.022131947937617618, 0.04309974177250355, 0.06994244827149917, 0.09728648227711462, 0.11840342673007132, 0.12808994400024723 … 0.030253423626318457, 0.01839898433655682, 0.010526447624071464, 0.0056824802607487405, 0.0029003690092184262, 0.001400239574849708, 0.000637395008304187, 0.00027024456233873585, 0.00010258104913936635, 2.9869033545987638e-5]\n",
" [5.463500926115185e-5, 0.0005362977277322803, 0.002632063780834687, 0.008611727482406828, 0.021131947356915853, 0.0414832477289767, 0.0678607351702997, 0.09515050401749464, 0.11673595713683114, 0.12730273393440114 … 0.031546146383514, 0.0193378343480595, 0.01115088615466354, 0.006066406715479799, 0.003119849590630032, 0.0015171842952738592, 0.0006953198690663203, 0.0002965727218072246, 0.00011311523060237527, 3.303976600748797e-5]\n",
" [5.113821406283989e-5, 0.0005053587915100225, 0.0024969220907984325, 0.008224613331134405, 0.020318086771559336, 0.04015460841936493, 0.06613047587106706, 0.0933502483113847, 0.11530041767314075, 0.12658608120497547 … 0.032656719678993154, 0.020151792296875495, 0.011696809204803888, 0.0064046803186089555, 0.0033146349932013765, 0.0016216612208947807, 0.0007473769418369354, 0.0003203526149657739, 0.00012266695883636174, 3.5921777585330676e-5]\n",
" [4.884963042407508e-5, 0.00048497979629751894, 0.002407325339105334, 0.007966227731199916, 0.019771002967849914, 0.03925464409406495, 0.06494833043876316, 0.09210713005978047, 0.1142932430994254, 0.1260630065535964 … 0.033435116024404465, 0.020726336760383778, 0.012084650696907163, 0.006646439722936009, 0.0034546165264241384, 0.001697123177418708, 0.0007851460370756075, 0.00033767128486040854, 0.00012964367781351675, 3.803074157375056e-5]\n",
" [4.854652544487309e-5, 0.00048226083812587207, 0.0023953391153449773, 0.007931545908588166, 0.019697313421565305, 0.03913296675802833, 0.06478782421122545, 0.09193746687092204, 0.11415472291756043, 0.12598973030245025 … 0.03354213827580439, 0.02080560602905177, 0.01213833009810171, 0.006679998249209769, 0.003474099657798998, 0.0017076520306958964, 0.0007904272546898002, 0.00034009737678657246, 0.00013062238927784753, 3.8326854796012915e-5]"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = solve(prob);\n",
"#sol = solve(prob, Tsit5(); dt=0.01)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1001×25 Matrix{Float64}:\n",
" 1.0 0.0 0.0 … 0.0 0.0\n",
" 0.951348 0.0474487 0.00118326 7.26426e-55 1.4127e-58\n",
" 0.905288 0.0900777 0.00448144 -5.35743e-53 -1.04593e-56\n",
" 0.861672 0.128286 0.00954967 -3.47055e-44 -6.34263e-47\n",
" 0.820359 0.162442 0.0160828 2.08843e-38 1.51351e-40\n",
" 0.781218 0.192883 0.0238115 … -2.41008e-38 -2.0081e-40\n",
" 0.744126 0.219923 0.0324984 8.55186e-35 1.18805e-36\n",
" 0.708967 0.243846 0.041935 -4.09686e-35 -9.33308e-37\n",
" 0.675633 0.264919 0.0519382 5.27187e-32 1.08597e-33\n",
" 0.64402 0.283385 0.0623482 3.71941e-32 6.06261e-34\n",
" 0.614034 0.299468 0.0730261 … 4.62317e-31 6.84043e-33\n",
" 0.585582 0.313374 0.0838506 3.57045e-29 1.0072e-30\n",
" 0.558582 0.325293 0.0947176 4.53676e-29 1.09576e-30\n",
" ⋮ ⋱ \n",
" 4.87265e-5 0.000483872 0.00240244 0.000130037 3.81498e-5\n",
" 4.87096e-5 0.00048372 0.00240177 … 0.000130091 3.81661e-5\n",
" 4.86928e-5 0.000483569 0.0024011 0.000130145 3.81824e-5\n",
" 4.8676e-5 0.000483418 0.00240044 0.000130199 3.81986e-5\n",
" 4.86594e-5 0.000483269 0.00239978 0.000130252 3.82148e-5\n",
" 4.86428e-5 0.00048312 0.00239912 0.000130305 3.82308e-5\n",
" 4.86262e-5 0.000482972 0.00239847 … 0.000130358 3.82468e-5\n",
" 4.86098e-5 0.000482825 0.00239782 0.00013041 3.82627e-5\n",
" 4.85934e-5 0.000482678 0.00239717 0.000130463 3.82785e-5\n",
" 4.85772e-5 0.000482533 0.00239653 0.000130515 3.82943e-5\n",
" 4.8561e-5 0.000482388 0.00239589 0.000130566 3.83099e-5\n",
" 4.85448e-5 0.000482244 0.00239526 … 0.000130618 3.83255e-5"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tvec = 0:0.05:50\n",
"xvec = sol(tvec).u'\n",
"xvec = permutedims(reshape(hcat(xvec...),(length(xvec[1]), length(xvec))))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"surface(tvec, 0:n_upperbound, xvec', \n",
" camera = (25, 30), \n",
" colour = :coolwarm, \n",
" xlabel = \"time (sec)\", \n",
" ylabel = \"n (# of molecules)\",\n",
" zlabel = \"p(t, n)\",\n",
" title = \"Degradation & Production\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) validation with the steady-state solution $\\phi(n)$ and Gillespie SSA"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\phi(n)$ (Exercise 1 - Problem 3c)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"compute_ϕ = function(n_upperbound, κ₁, κ₂)\n",
" \n",
" ϕ = zeros(n_upperbound+1)\n",
"\n",
" ϕ[1] = 1.0\n",
" ϕ[2] = κ₂*V/κ₁ * ϕ[1];\n",
"\n",
" for n in 1:n_upperbound-1\n",
" ϕ[n+2] = (κ₁*n*ϕ[n+1] + κ₂*V*ϕ[n+1] - κ₂*V*ϕ[n])/(κ₁*(n+1));\n",
" end\n",
" return ϕ./(sum(ϕ))\n",
"end\n",
"\n",
"ϕ = compute_ϕ(n_upperbound, κ₁, κ₂);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gillespie algorithm (Exercise 1 - Problem 3b)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#7 (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Compute the total propensity\n",
"using Distributions, StatsBase \n",
"\n",
"propensities = function(κ, x, ν, Nᵣ)\n",
" ν_eval = zeros(Nᵣ)\n",
" for r in 1:Nᵣ\n",
" ν_eval[r] = ν[r](κ[r], x) \n",
" end\n",
" return ν_eval\n",
"end\n",
"\n",
"gillespie_alg = function(SM, κ, ν, x₀, stoptime)\n",
" t = 0.0\n",
" x = x₀\n",
" \n",
" Nₛ, Nᵣ = size(SM)\n",
" \n",
" tt = [t]\n",
" xx = copy(x) \n",
" \n",
" k=1\n",
" \n",
" while (t <= stoptime)\n",
" \n",
" # step 1\n",
" ν_eval = propensities(κ, x, ν, Nᵣ) \n",
" α = sum(ν_eval) \n",
"\n",
" # step 2\n",
" τ = rand(Exponential(1/α)) \n",
" t += τ\n",
"\n",
" # step 3\n",
" index_j = StatsBase.sample(1:Nᵣ, Weights(ν_eval./α))\n",
" \n",
" x = x + SM[:, index_j]\n",
" append!(tt, t)\n",
" xx = cat(xx, x, dims=2)\n",
" end\n",
" \n",
" return tt, xx'\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#11 (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Number of reactions Nᵣ = 2\n",
"\n",
"# Compute the stoichiometric matrix \n",
"Products = [0; 1]\n",
"Reactants = [1; 0]\n",
"SM = (Products-Reactants)'\n",
"\n",
"# Compute the propensity function for each of the reactions\n",
"ν₁ = function(κ₁, x) return κ₁*x[1] end # degradation\n",
"ν₂ = function(κ₂, x) return κ₂*V end # production"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"N = 1e3 #number of realisations\n",
"n_upperbound = 24\n",
"p = zeros(n_upperbound+1) #vector of probabilities\n",
"\n",
"for n in 1:N \n",
" tt, xx = gillespie_alg(SM, [κ₁, κ₂], [ν₁, ν₂], [0], 50)\n",
" A = xx[end] \n",
" if (A <= n_upperbound) \n",
" p[xx[end]+1] += 1\n",
" end\n",
"end\n",
"\n",
"p /= N;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Histogram (Gillespie)"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"bar(collect(0:24), p, bar_width = 0.1,\n",
" colour = :lightgrey,\n",
" label = \"Gillespie SSA\",\n",
" xlabel = \"number of molecules\",\n",
" ylabel = \"stationary distribution\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Numerical solution $p(t, n)$ at $t=T$."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"plot(collect(0:n_upperbound), sol(tspan[2]), \n",
" label = \"numeric solution\",\n",
" title = \"p(t, n) @ t = \".*string(tspan[2]).*\"sec\",\n",
" linewidth = 2,\n",
" colour = :black);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Steady-state solution $\\phi(n)$"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!(collect(0:n_upperbound), ϕ,\n",
" linewidth = 3,\n",
" colour = :red,\n",
" label = \"steady-state solution phi(n)\", linestyle=:dash)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Dimerization and production"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) numerical solution to the chemical master equation "
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"ME_dim_and_prod! (generic function with 1 method)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ME_dim_and_prod!(du, u, p, t)\n",
" \n",
" κ₁, κ₂, V = p\n",
"\n",
" u ./= sum(u) \n",
"\n",
" u_padded = [0; u; 0; 0]\n",
" \n",
" n = 0:length(u)-1 \n",
" A = -κ₁ * (1/V) * (n.*(n.-1)/2 .- κ₂*V)\n",
" B = κ₁ * (1/V) * (n.+2) .* (n.+1) / 2 \n",
" C = κ₂ * (1/V)\n",
"\n",
" du .= A .* u_padded[2:end-2] + B .* u_padded[4:end] + C * u_padded[1:end-3]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"n_upperbound = 24\n",
"u₀ = [1; zeros(n_upperbound)]\n",
"\n",
"tspan = (0, 50)\n",
"κ₁, κ₂, V = [0.005, 1, 1]\n",
"p = [κ₁, κ₂, V];"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"prob = ODEProblem(ME_dim_and_prod!, u₀, tspan, p)\n",
"#sol = solve(prob)\n",
"sol = solve(prob, Tsit5(); dt=0.1, adaptive=false);"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5001×25 Matrix{Float64}:\n",
" 1.0 0.0 0.0 … 0.0 0.0 0.0\n",
" 1.00005 0.00995058 4.94922e-5 0.0 0.0 0.0\n",
" 1.0001 0.0198036 0.000196038 0.0 0.0 0.0\n",
" 1.00015 0.029561 0.000436807 0.0 0.0 0.0\n",
" 1.0002 0.0392246 0.00076906 0.0 0.0 0.0\n",
" 1.00024 0.0487962 0.00119014 … 0.0 0.0 0.0\n",
" 1.00029 0.0582776 0.00169749 0.0 0.0 0.0\n",
" 1.00034 0.0676704 0.00228863 0.0 0.0 0.0\n",
" 1.00039 0.0769762 0.00296117 0.0 0.0 0.0\n",
" 1.00043 0.0861968 0.00371282 0.0 0.0 0.0\n",
" 1.00048 0.0953337 0.00454136 … 0.0 0.0 0.0\n",
" 0.909157 0.0956781 0.00503361 0.0 0.0 0.0\n",
" 0.909202 0.10464 0.00602032 0.0 0.0 0.0\n",
" ⋮ ⋱ \n",
" 6.27812e-8 1.26663e-6 1.24572e-5 0.00806832 0.00355944 0.00149434\n",
" 6.33553e-8 1.27822e-6 1.25713e-5 … 0.0081428 0.00359231 0.00150814\n",
" 5.81508e-8 1.17322e-6 1.15387e-5 0.00747459 0.00329752 0.00138439\n",
" 5.87249e-8 1.18481e-6 1.16527e-5 0.00754907 0.00333039 0.00139819\n",
" 5.92989e-8 1.1964e-6 1.17667e-5 0.00762356 0.00336326 0.00141199\n",
" 5.98729e-8 1.20799e-6 1.18807e-5 0.00769805 0.00339613 0.00142579\n",
" 6.04469e-8 1.21957e-6 1.19947e-5 … 0.00777255 0.003429 0.0014396\n",
" 6.10208e-8 1.23116e-6 1.21087e-5 0.00784704 0.00346187 0.0014534\n",
" 6.15947e-8 1.24274e-6 1.22227e-5 0.00792154 0.00349474 0.0014672\n",
" 6.21686e-8 1.25433e-6 1.23367e-5 0.00799603 0.00352761 0.00148101\n",
" 6.27425e-8 1.26591e-6 1.24507e-5 0.00807053 0.00356049 0.00149481\n",
" 6.33163e-8 1.27749e-6 1.25647e-5 … 0.00814503 0.00359336 0.00150861"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tvec = 0:0.01:50\n",
"xvec = sol(tvec).u'\n",
"xvec = permutedims(reshape(hcat(xvec...),(length(xvec[1]), length(xvec))))"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"surface(tvec, 0:n_upperbound, xvec', \n",
" camera = (25, 30), \n",
" colour = :coolwarm, \n",
" xlabel = \"time (sec)\", \n",
" ylabel = \"n (# of molecules)\",\n",
" zlabel = \"p(t, n)\",\n",
" title = \"Dimerization & Production\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) steady-state solution to the chemical master equation "
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"25-element Vector{Float64}:\n",
" 5.353340645375288e-8\n",
" 1.0843089071676345e-6\n",
" 1.0706681290750576e-5\n",
" 6.87183667142588e-5\n",
" 0.0003225301930012232\n",
" 0.0011808492184844418\n",
" 0.003513169761026681\n",
" 0.008736776061023249\n",
" 0.018541487675851714\n",
" 0.03411648770224447\n",
" 0.05511342173087868\n",
" 0.07896715568289732\n",
" 0.10120440581236947\n",
" 0.11684538914055229\n",
" 0.12227407482978915\n",
" 0.11659178112952047\n",
" 0.1017723162279849\n",
" 0.08165939910695519\n",
" 0.06044957494804622\n",
" 0.04142160733532141\n",
" 0.02635168492246991\n",
" 0.01560714919902379\n",
" 0.008626994167528079\n",
" 0.0044608465893508755\n",
" 0.002162335675361753"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SpecialFunctions\n",
"compute_ϕ = function(n_upperbound, κ₁, κ₂)\n",
" \n",
" ϕ = zeros(n_upperbound+1)\n",
"\n",
" for n in 0:n_upperbound\n",
" ϕ[n+1] = (1/factorial(big(n))) * sqrt(κ₂*2*V^2/κ₁)^n * besseli(n-1, 2*sqrt(κ₂*2*V^2/κ₁))\n",
" end\n",
" return ϕ./(sum(ϕ))\n",
"end\n",
"\n",
"ϕ = compute_ϕ(n_upperbound, κ₁, κ₂)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) validation of (a) and (b) with the Gillespie algorithm"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#19 (generic function with 1 method)"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Number of reactions Nᵣ = 2\n",
"# Compute the stoichiometric matrix \n",
"Products = [0; 1]\n",
"Reactants = [2; 0]\n",
"SM = (Products-Reactants)'\n",
"\n",
"# Compute the propensity function for each of the reactions\n",
"ν₁ = function(κ₁, x) return κ₁*x[1]*(x[1]-1)/(2*V) end\n",
"ν₂ = function(κ₂, x) return κ₂*V end "
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"N = 1e4 #number of realisations\n",
"n_upperbound = 24\n",
"p = zeros(n_upperbound+1) #vector of probabilities\n",
"\n",
"for n in 1:N \n",
" tt, xx = gillespie_alg(SM, [κ₁, κ₂], [ν₁, ν₂], [0], 50.0)\n",
" A = xx[end] \n",
" if (A <= n_upperbound) \n",
" p[xx[end]+1] += 1\n",
" end\n",
"end\n",
"\n",
"p /= N;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Histogram (Gillespie)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"bar(collect(0:24), p, bar_width = 0.1,\n",
" colour = :lightgrey,\n",
" label = \"Gillespie SSA\",\n",
" xlabel = \"number of molecules\",\n",
" ylabel = \"stationary distribution\",\n",
" legend = :topleft);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Numerical solution $p(t, n)$ at $t=T$."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"plot!(collect(0:n_upperbound), sol(tspan[2]), \n",
" label = \"numeric solution\",\n",
" title = \"p(t, n) @ t = \".*string(tspan[2]).*\"sec\",\n",
" linewidth = 2,\n",
" colour = :black); "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Steady-state solution $\\phi(n)$"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!(collect(0:n_upperbound), ϕ, \n",
" linewidth = 3,\n",
" colour = :red,\n",
" label = \"steady-state solution phi(n)\", linestyle=:dash)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}