Skip to content

Commit c07dbd9

Browse files
committed
up
1 parent fe29787 commit c07dbd9

File tree

2 files changed

+29
-22
lines changed

2 files changed

+29
-22
lines changed

docs/Project.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
77
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
88
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
99
HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
10+
JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
1011
Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
1112
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
1213
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

docs/src/catalyst_functionality/example_networks/basic_CRN_library.md

Lines changed: 28 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -9,13 +9,14 @@ bd_process = @reaction_network begin
99
(p,d), ∅ <--> X
1010
end
1111
```
12-
Next we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
12+
Next, we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
1313
```@example crn_library_birth_death
1414
u0 = [:X => 1]
1515
tspan = (0.0, 10.0)
1616
ps = [:p => 1, :d => 0.2]
17+
nothing # hide
1718
```
18-
We can now simulate our model using all three interpretations. First we perform a reaction rate equation-based ODE simulation:
19+
We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:
1920
```@example crn_library_birth_death
2021
using OrdinaryDiffEq
2122
oprob = ODEProblem(bd_process, u0, tspan, ps)
@@ -43,11 +44,11 @@ using Plots
4344
oplt = plot(osol; title = "Reaction rate equation (ODE)")
4445
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
4546
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
46-
plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
47+
plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))
4748
```
4849

4950
## Michaelis-Menten enzyme kinetics
50-
[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions it can be simplified to a singe function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
51+
[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
5152
```@example crn_library_michaelis_menten
5253
using Catalyst
5354
mm_system = @reaction_network begin
@@ -79,7 +80,8 @@ using Plots
7980
oplt = plot(osol; title = "Reaction rate equation (ODE)")
8081
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
8182
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
82-
plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
83+
plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))
84+
plot!(bottom_margin = 3Plots.Measures.mm) # hide
8385
```
8486

8587
## SIR infection model
@@ -91,7 +93,7 @@ sir_model = @reaction_network begin
9193
β, I --> R
9294
end
9395
```
94-
First we perform a deterministic ODE simulation:
96+
First, we perform a deterministic ODE simulation:
9597
```@example crn_library_sir
9698
using OrdinaryDiffEq, Plots
9799
u0 = [:S => 99, :I => 1, :R => 0]
@@ -101,9 +103,9 @@ ps = [:α => 0.001, :β => 0.01]
101103
# Solve ODEs.
102104
oprob = ODEProblem(sir_model, u0, tspan, ps)
103105
osol = solve(oprob, Tsit5())
104-
plot(osol; title = "Reaction rate equation (ODE)")
106+
plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))
105107
```
106-
Next we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of a outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
108+
Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
107109
```@example crn_library_sir
108110
using JumpProcesses
109111
dprob = DiscreteProblem(sir_model, u0, tspan, ps)
@@ -119,12 +121,13 @@ jsol3 = solve(jprob, SSAStepper(); seed=3) # hide
119121
jplt1 = plot(jsol1; title = "Outbreak")
120122
jplt2 = plot(jsol2; title = "Outbreak")
121123
jplt3 = plot(jsol3; title = "No outbreak")
122-
plot(jplt1, jplt2, jplt3; size=(800,700), layout = (3,1))
124+
plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
123125
```
124126

125127
## The Wilhelm model
126128
The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
127129
```@example crn_library_wilhelm
130+
using Catalyst
128131
wilhelm_model = @reaction_network begin
129132
k1, Y --> 2X
130133
k2, 2X --> X + Y
@@ -144,12 +147,13 @@ oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
144147
oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
145148
osol1 = solve(oprob1, Tsit5())
146149
osol2 = solve(oprob2, Tsit5())
147-
oplt1 = plot(osol1; idxs = :X, label = "X(0) = 1.5")
148-
oplt2 = plot!(osol2; idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,700))
150+
plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
151+
plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))
152+
plot!(bottom_margin = 3Plots.Measures.mm) # hide
149153
```
150154

151155
## Simple self-activation loop
152-
The simplest self-activation loop consist of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
156+
The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
153157
```@example crn_library_self_activation
154158
using Catalyst
155159
sa_loop = @reaction_network begin
@@ -161,7 +165,7 @@ A simple example of such a loop is a transcription factor which activates its ow
161165

162166
We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
163167
```@example crn_library_self_activation
164-
using OrdinaryDiffEq, Plots
168+
using JumpProcesses, OrdinaryDiffEq, Plots
165169
u0 = [:X => 4]
166170
tspan = (0.0, 1000.0)
167171
ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
@@ -174,12 +178,13 @@ jprob = JumpProblem(sa_loop, dprob, Direct())
174178
jsol = solve(jprob, SSAStepper())
175179
jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
176180
177-
plot(osol; label = "Reaction rate equation (ODE)")
178-
plot!(jsol; label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,600))
181+
plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
182+
plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide="X", size = (800,350))
183+
plot!(bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
179184
```
180185

181186
## The Brusselator
182-
The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
187+
The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well-known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
183188
```@example crn_library_brusselator
184189
using Catalyst
185190
brusselator = @reaction_network begin
@@ -189,7 +194,7 @@ brusselator = @reaction_network begin
189194
1, X --> ∅
190195
end
191196
```
192-
It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this results is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contain a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in he second case the system is able to generate oscillations.
197+
It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
193198
```@example crn_library_brusselator
194199
using OrdinaryDiffEq, Plots
195200
u0 = [:X => 1.0, :Y => 1.0]
@@ -204,11 +209,11 @@ osol2 = solve(oprob2, Rodas5P())
204209
oplt1 = plot(osol1; title = "No Oscillation")
205210
oplt2 = plot(osol2; title = "Oscillation")
206211
207-
plot(oplt1, oplt2; layout = (1,2), size(800,700))
212+
plot(oplt1, oplt2; lw = 3, layout = (2,1), size = (800,600))
208213
```
209214

210215
## The repressilator
211-
The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that is able to generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) which production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
216+
The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
212217
```@example crn_library_brusselator
213218
using Catalyst
214219
repressilator = @reaction_network begin
@@ -218,7 +223,7 @@ repressilator = @reaction_network begin
218223
d, (X, Y, Z) --> ∅
219224
end
220225
```
221-
Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model is able to sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
226+
Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
222227
```@example crn_library_brusselator
223228
using OrdinaryDiffEq, StochasticDiffEq, Plots
224229
u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
@@ -242,7 +247,7 @@ ssol2 = solve(sprob2, ImplicitEM(); seed = 100) # hide
242247
splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
243248
splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
244249
245-
plot(oplt1, oplt2, splt1, splt2; layout = (2,2), size = (800,600))
250+
plot(oplt1, oplt2, splt1, splt2; lw = 2, layout = (2,2), size = (800,600))
246251
```
247252

248253
## The Willamowski–Rössler model
@@ -259,7 +264,7 @@ wr_model = @reaction_network begin
259264
k7, Z --> ∅
260265
end
261266
```
262-
Here we first simulate the model for a single initial conditions, showing in both time-state space and phase space how how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
267+
Here we first simulate the model for a single initial condition, showing in both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
263268
```@example crn_library_chaos
264269
using OrdinaryDiffEq, Plots
265270
u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
@@ -271,4 +276,5 @@ sol = solve(oprob, Rodas5P())
271276
plt1 = plot(sol; title = "Time-state space")
272277
plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
273278
plot(plt1, plt2; layout = (1,2), size = (800,400))
279+
plot!(bottom_margin = 3Plots.Measures.mm) # hide
274280
```

0 commit comments

Comments
 (0)