1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
|
<html>
<head>
<title>Birth-Death</title>
</head>
<body>
<h1>Birth-Death</h1>
<p>
Consider the linear birth-death process presented in Section 1.3 of
<a href="http://www.staff.ncl.ac.uk/d.j.wilkinson/smfsb/">Stochastic
Modelling for Systems Biology</a>. The text introduces some differences
between continuous, deterministic modelling and discrete, stochastic modelling.
Let <em>X(t)</em> be the population of bacteria which reproduce at a rate
of λ and die at a rate of μ. The continuous model of this process
is the differential equation
</p>
<p align="center">
<em>
X'(t) = (λ - μ) X(t),
X(0) = x<sub>0</sub>
</em>
</p>
<p>
which has the solution
</p>
<p align="center">
<em>X(t) = x<sub>0</sub></em> e<sup>(λ - μ)<em>t</em></sup>.
</p>
<p>
We can numerically solve the continuous, deterministic model in Cain.
Open the file <tt>examples/cain/BirthDeath.xml</tt> and select the
<tt>Birth Death 0 1</tt> model, which has parameter values λ = 0
and μ = 1. In the method list, clone the <tt>Direct</tt> method
by selecting it and then clicking the clone button <img src="editcopy.png">.
Click twice (don't double-click) on the result and change its name to
"Deterministic."
In the method editor, select the <tt>Time
Homogeneous</tt> and then the
<tt>Time Series, Deterministic</tt> categories with the default
method (<tt>ODE, Integrate Reactions</tt>) and options
(<tt>Runge-Kutta, Cash-Karp</tt>).
</p>
<p>
Check that the species <em>X</em>
has been selected in the recorder panel. In the launcher panel set the
number of trajectories to 1 and
click the launch button <img src="launch.png">
to generate the solution. You can visualize the solution by clicking the
plot button <img src="plot.png"> in the output list to bring up
the plotting configuraton window (shown below).
</p>
<p align="center">
<img src="ExamplesBirthDeathDeterministicSpeciesPlotWindow.png">
</p>
<p>
Add axes labels as shown above and then click the plot button.
The population as a function of time is plotted below.
</p>
<p align="center">
<img src="BirthDeath-0-1-ODE-Populations.jpg">
</p>
<p>
The ordinary differential equation (ODE) integration method in Cain
actually solves a different formulation of the model than the population-based
formulation above.
for the method.)
As the method name "<tt>ODE, Integrate Reactions</tt>"
suggests, it integrates the reaction counts. Instead of combining the
effect of the two reactions to formulate the problem in terms of the
species population <em>X</em>, the birth reaction <em>B</em> and the
death reaction <em>D</em> are each represented.
The birth-death process can be modelled with the following
set of differential equations:
</p>
<p align="center">
<em>
B'(t) = λ X(t), B(0) = 0<br>
D'(t) = μ X(t), D(0) = 0<br>
X'(t) = B'(t) - D'(t), X(0) = x<sub>0</sub>
</em>
</p>
<p>
It is perhaps instructive to translate the equations into plain English.
<ul>
<li> The birth rate is &lambda. We start counting the number of births
at the initial time <em>t = 0</em>.
<li> The death rate is &mu. We start counting the number of deaths
at the initial time <em>t = 0</em>.
<li> The rate of change is the population <em>X</em> is the difference
between the birth rate and the death rate. At the initial time
<em>t = 0</em>, the population is <em>x<sub>0</sub></em>.
</ul>
For λ ≠ μ the system of equations has the solution
</p>
<p align="center">
<em>B(t) = λ x<sub>0</sub></em>
(1 - e<sup>(λ - μ)<em>t</em></sup>) / (μ - λ)<br>
<em>D(t) = μ x<sub>0</sub></em>
(1 - e<sup>(λ - μ)<em>t</em></sup>) / (μ - λ)<br>
<em>X(t) = x<sub>0</sub></em> e<sup>(λ - μ)<em>t</em></sup>.
</p>
<p>
While the population solutions are the same, the reaction-based
formulation of the model carries more information than the
population-based formulation. The population depends only
on the difference λ - μ. However, the reaction counts depend
on the two parameters separately. For λ = 0 and μ = 1
no birth reactions occur.
</p>
<p>
While the reaction counts are calculated, by default they are not
recorded in the simulation output. Delete the solution
that we generated previously by clicking the delete button
<img src="cancel.png"> in the simulation output panel.
Then select the reactions tab in the recorder panel and select both
the birth and death reactions. Now run the simulation again.
In the plotting window select the <tt>Cumulative reactions</tt>
radio button, move the legend to the upper left corner,
and change the Y axes label to "Reaction Counts."
</p>
<p align="center">
<img src="ExamplesBirthDeathDeterministicReactionsPlotWindow.png">
</p>
<p>
Plotting will result in a figure like the one below.
</p>
<p align="center">
<img src="BirthDeath-0-1-ODE-Reactions.jpg">
</p>
<p>
For λ = 10 and μ = 11, the population is the same, but the
reaction counts differ. Below is a plot of the reaction counts for this case.
</p>
<p align="center">
<img src="BirthDeath-10-11-ODE-Reactions.jpg">
</p>
<p>
Now consider the discrete stochastic model, which has reaction
propensities instead of deterministic reaction rates.
This model is composed of the birth reaction
<em>X → 2 X</em> and the death reaction <em>X →</em> 0 which
have propensities λ<em>X</em> and μ<em>X</em>, respectively.
First we will generate a trajectory that records all of the reactions.
Select the <tt>Birth Death 0 1</tt> model. Then clone the
"Deterministic" method (with the clone button
<img src="editcopy.png">) and rename it "DirectAll."
For this method select the <tt>Time Homogeneous</tt> and
<tt>Time Series, All Reactions</tt> categories with the default method
and options. Then generate a trajectory with
the launch button. Below is a plot of the species
populations. We see that the population changes by discrete amounts.
</p>
<p align="center">
<img src="BirthDeath-0-1-1-DAR-Populations.jpg">
</p>
<p>
Below we use Cain to reproduce the results in the text that demonstrate how
increasing λ + μ while holding λ - μ = -1 increases the
volatility in the system. For each test, we generate an ensemble of five
trajectories and plot these populations along with the
deterministic solution. Note that hitting the "Plot" button
in the plotting window adds to the current figure while
"New plot" creates a new figure. Also note that you can
click on the grid elements and columns in the plotting window to
customize the appearance of the curves.
</p>
<p align="center">
<img src="BirthDeath-0-1-DAR-ODE-Populations.jpg"><br>
λ = 0, μ = 1
</p>
<p align="center">
<img src="BirthDeath-3-4-DAR-ODE-Populations.jpg"><br>
λ = 3, μ = 4
</p>
<p align="center">
<img src="BirthDeath-7-8-DAR-ODE-Populations.jpg"><br>
λ = 7, μ = 8
</p>
<p align="center">
<img src="BirthDeath-10-11-DAR-ODE-Populations.jpg"><br>
λ = 10, μ = 11
</p>
<p>
For a simple problem like this we can store and visualize all of the
reactions. However, for more complicated models (or longer running
times) generating a suite of trajectories may involve billions of reaction
events. Storing, and particularly plotting, that much data could be
time consuming or just impossible on your computer. Thus instead of
storing all of the reaction events, one typically
stores snapshots of the populations and reaction counts at set points
in time. Again clone the "Deterministic" method and rename
it to "Direct."
For this method select the <tt>Time Homogeneous</tt> and
<tt>Time Series, Uniform</tt> categories with the default method
and options. In the reactions tab of the recorder panel select
both reactions.
For each test, we generate an ensemble of ten
trajectories and plot the species populations and the cumulative
reaction counts. Note that because we are only sampling the
state, we don't see the same "noisiness" in the trajectories.
</p>
<p align="center">
<img src="BirthDeath-0-1-Populations.jpg"><br>
<img src="BirthDeath-0-1-Reactions.jpg"><br>
λ = 0, μ = 1
</p>
<p align="center">
<img src="BirthDeath-3-4-Populations.jpg"><br>
<img src="BirthDeath-3-4-Reactions.jpg"><br>
λ = 3, μ = 4
</p>
<p align="center">
<img src="BirthDeath-7-8-Populations.jpg"><br>
<img src="BirthDeath-7-8-Reactions.jpg"><br>
λ = 7, μ = 8
</p>
<p align="center">
<img src="BirthDeath-10-11-Populations.jpg"><br>
<img src="BirthDeath-10-11-Reactions.jpg"><br>
λ = 10, μ = 11
</p>
<p>
In the plotting window you may choose between a number of
visualization options. First you choose what to plot:
<ul>
<li> <tt>Species</tt> - The species populations.
<li> <tt>Cumulative Reactions</tt> - The cumulative reaction counts.
<li> <tt>Binned Reactions</tt> - The reaction counts binned
for each frame.
</ul>
Next you choose how to display the data: either plot the trajectories
or plot their mean (and optionally their standard deviation).
Thus there are six ways of generating a plot for a given simulation
output. Below we show each of these for the
birth-death model with λ = 3 and μ = 4.
</p>
<p align="center">
<img src="BirthDeath-3-4-PopulationStatistics.jpg"><br>
Population Statistics
</p>
<p align="center">
<img src="BirthDeath-3-4-PopulationTrajectories.jpg"><br>
Population Trajectories
</p>
<p align="center">
<img src="BirthDeath-3-4-BinnedReactionCountStatistics.jpg"><br>
Binned Reaction Count Statistics
</p>
<p align="center">
<img src="BirthDeath-3-4-BinnedReactionCountTrajectories.jpg"><br>
Binned Reaction Count Trajectories
</p>
<p align="center">
<img src="BirthDeath-3-4-CumulativeReactionCountStatistics.jpg"><br>
Cumulative Reaction Count Statistics
</p>
<p align="center">
<img src="BirthDeath-3-4-CumulativeReactionCountTrajectories.jpg"><br>
Cumulative Reaction Count Trajectories
</p>
<p>
After selecting the plotting method you can customize the appearance
by editing the table that lists the species or reactions. See the
<a href="VisualizationPlottingTimeSeries.htm">Plotting Time Series Data</a>
section of the
<a href="Visualization.htm">Visualization and Analysis</a> chapter
for details.
</p>
<p>
Consider how the value of λ affects the population of X at time
t = 2.5. From the plots above it appears that with increasing λ there
is greater variance in the population, and also a greater likelihood of
extinction (X = 0). However, it is not possible to quantify these
observations by looking at trajectory plots. Recording histograms of the state
is the right tool for this. Once more clone the
"Deterministic" method. This time rename the result
"Histogram." Select the <tt>Time Homogeneous</tt> and
<tt>Histograms, Transient Behavior</tt> categories with the default
solver method and options. Since we are only interested in the population at
t = 2.5, we set the recording time to that value and set the number of frames
to 1. We set the number of bins to 128 to obtain a high resolution
histogram and launch a simulation with
100,000 trajectories for each value of λ. In the plotting
window we switch to the histograms tab. When plotting histograms
you choose the species and the frame (time). You can also choose colors
and enter a title and axes labels if you like. The plot configuration
window is shown below.
</p>
<p align="center">
<img src="ExamplesBirthDeathPlotWindowHistograms.png">
</p>
<p>
The histograms for each value of λ are shown below. We see that for
λ = 0, the mode (most likely population) is 4, but for λ =
3, 7, or 10, the mode is 0. The likelihood of extinction increases with
increasing λ.
</p>
<p align="center">
<img src="BirthDeath2.5Lambda0.jpg"><br>
λ = 0, μ = 1
</p>
<p align="center">
<img src="BirthDeath2.5Lambda3.jpg"><br>
λ = 3, μ = 4
</p>
<p align="center">
<img src="BirthDeath2.5Lambda7.jpg"><br>
λ = 7, μ = 8
</p>
<p align="center">
<img src="BirthDeath2.5Lambda10.jpg"><br>
λ = 10, μ = 11
</p>
</body>
</html>
|
