@@ -88,10 +88,11 @@ plot(sol)
88
88
## Adding events
89
89
Our current model is unrealistic in assuming the cell will grow exponentially
90
90
forever. Let's modify it such that the cell divides in half every time its
91
- volume reaches a size of ` 2 ` . Note, we will only keep track of one cell, and
92
- hence follow a specific lineage of the system. To do this we can create a
93
- continuous event using the ModelingToolkit symbolic event interface and attach
94
- it to our system. Please see the associated [ ModelingToolkit
91
+ volume reaches a size of ` 2 ` . We also assume we lose half of the protein upon
92
+ division. Note, we will only keep track of one cell, and hence follow a specific
93
+ lineage of the system. To do this we can create a continuous event using the
94
+ ModelingToolkit symbolic event interface and attach it to our system. Please see
95
+ the associated [ ModelingToolkit
95
96
tutorial] ( https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/ ) for more
96
97
details on the types of events that can be represented symbolically. A
97
98
lower-level approach for creating events via the DifferentialEquations.jl
@@ -104,20 +105,20 @@ using Catalyst, DifferentialEquations, Plots
104
105
105
106
@parameters λ = 1.0
106
107
@variables t V(t) = 1.0
108
+ @species P(t) = 0.0
107
109
D = Differential(t)
108
110
eq = D(V) ~ λ * V
109
- rx1 = @reaction $V, 0 --> P
110
- rx2 = @reaction 1.0, P --> 0
111
+ rx1 = @reaction $V, 0 --> $ P
112
+ rx2 = @reaction 1.0, $ P --> 0
111
113
```
112
114
Before creating our ` ReactionSystem ` we make the event.
113
115
``` @example ceq3
114
- # every 1.0 time unit we half the volume of the cell
115
- continuous_events = [V ~ 2.0] => [V ~ V/2]
116
+ # every 1.0 time unit we half the volume of the cell and the number of proteins
117
+ continuous_events = [V ~ 2.0] => [V ~ V/2, P ~ P/2 ]
116
118
```
117
119
We can now create and simulate our model
118
120
``` @example ceq3
119
121
@named rs = ReactionSystem([rx1, rx2, eq], t; continuous_events)
120
- setdefaults!(rs, [:P => 0.0])
121
122
122
123
oprob = ODEProblem(rs, [], (0.0, 10.0))
123
124
sol = solve(oprob, Tsit5())
0 commit comments