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
|
%%%%%%%%%%%
%
% In this example, I will show a simple implementation of the Gillespie
% algorithm. The example we will work on is a birth/death, or creation
% destruction, process. Two things can happen:
% 1 - the creation of a molecule/particle
% 2 - the destruction of a molecule/particle
% We will only consider one species of particles, and they will be
% independent of each other. We will monitor the overall count of
% particles, which we will call nn. (I am using two n, because otherwise
% when you want to find the variable in the code using the 'find' function,
% you will have a lot of trouble.)
%
% The general idea is that we first define the two 'reaction channels'.
% This means we define the events that will take place in the two cases of
% creation or destruction, which I will call 'reactions'. Each of these
% reactions has an associated 'reaction channel'. When the reaction channel
% is fired, a specific change happens in the system. The Gillespie
% simulation consists of stochastically determined series of individual
% firing events, where always exactly one reaction channel is executed
% exactly once. After each of these reactions, the 'propensities' of each
% channel to fire are calculated, to tell us how likely it is that for each
% reaction channel to be fired in the next reaction, and how long we will
% have to wait till then.
%
%%%%%%%%%%%
%%%% Definition of model parameters
% In the model we are using, two events can take place - the creation and
% the destruction of a molecule. The first one takes place at a constant
% rate:
kk = 1; % Constant rate of molecule creation
% The second one, destruction, happens at a rate that is proportional to
% the number of molecules that currently exist. This rate will later be
% multiplied with the number nn.
gamma = 0.1; % Destruction rate
%%%% Definition of general simulation parameters
% As a first thing, we have to set how many individual reactions should
% maximally be simulated
R_max = 1000; % maximal simulation steps
%%%% Preallocation of containers for the simulation time course
% To make the code more efficient, we 'pre-allocate' the variabels we will
% use to store the molecule counts and time points of our Gillespie
% simulation. This will make the code run faster - when MatLab has to
% change the size (number of elements) of a variable, what it will do is
% to destroy the old variable, make a new one in the right size and work
% with it. This is slow, and when it is done for every reaction step, we
% will make a very slow program. The way around it is to reserve enough
% space in a variable beforehand, and then save all the molecule counts and
% reaction times into there.
nn_cont = zeros(1,R_max+1); % Container for nn time course
tt_cont = zeros(1,R_max+1); % Container for time points of reactions
%%%% Pre-drawing of random numbers for simulation
% Very much like the pre-allocation of variables, we draw all the random
% numbers we will need beforehand. MatLab is fast at making large sets of
% numbers at once, while calling for these activities is very slow. So,
% what we want to do is call for the activity once, and have all the random
% numbers ready to work with later.
rand_nums = rand(2,R_max); % TWO random numbers for each reaction step
%%%% Definition of state changes for each reaction channel
% The next thing to define is the change in molecule number that occurs
% when a specific reaction channel is fired. We have two reaction channels,
% so we need two state changes, one for each reaction channel.
state_changes = [+1,-1];
% The first one is for the 'creation' reaction channel, where one molecule
% is added. The second one is for the 'destruction' channel, where one
% molecule is removed.
%%%% The actual simluation
% To simulate the time course, we now go step-by-steo from reaction to
% reaction, until we have reached the maximal number of reactions. As a
% first thing, we initialize our simulation at the time point 0
tt = 0; % current time
nn = 0; % current molecule count
% We also need a counter that tells us how many reactions have occured
% until now.
R_counter = 0; % reaction counter
% Now we can run a while loop, which is repeated until we reach the maximum
% number of reactions
while R_counter < R_max
R_counter = R_counter + 1; % Tell the counter one more reaction happened
% As a first thing, we need to calculate the propensities for both of
% the reaction channels to fire. Let us first allocate them, and then
% assign them, for clarity's sake
aa = zeros(1,2); % Pre-allocated propensity vector, no values yet
% Now, we assign the propensity for a molecule to be created. This
% happens at a constant rate, so that the propensity is also constant.
aa(1) = kk;
% Next, we assign the propensity for a molecule to be destroyed. This
% happens at a rate that is proportional to the number of currently
% existent molecules - two molecules are double as likely to lead to
% the destruction of a molecule than one molecule is.
aa(2) = gamma.*nn;
% Now we draw the waiting time till the next reaction channel is
% fired. For this we use one of the pre-drawn random numbers
deltat_t = -log(rand_nums(1,R_counter))./sum(aa);
% We also need to find out which of the reaction channels is fired at
% that time. For this we have to choose from the possible reaction
% channels, according to their relative propensity. The idea here is to
% construct an interval that ranges from 0 to 1, [0,1], which contains
% different segments that correspond to the different reaction
% channels. Then we draw a random number between 0 and 1, and dependent
% on which segments it falls into, we choose the reaction channel that
% is fired. It would look a little bit like this:
%
% rand_num=0.5 x
% 0[-----------------------,---]1 <- interval with segments
% channel 1 2
%
% Here the random number for choosing the reaction channel is roughly
% 0.5, which means we choose the weighted segment that covers the
% point 0.5 in the interval. Here, that would refer to reaction channel
% 1, so the next channel that would be fired would be channel 1.
%
% You might realize that channels that have a length of 0 could mess up
% this procedure, so we have to remove them from the following
% computation. It looks a little complicated, but it will avoid a lot
% of trouble in simulations that are more complicated than this simple
% example.
% So let us first get only the info from the valid reaction channels,
% that is the ones with >0 propensity.
valid_inds = aa > 0; % Find the indices that refer to valid reactions
valid_aa = aa(valid_inds); % Use only valid reactions
valid_changes = state_changes(valid_inds); % Use only valid changes
% Next, we can construct the intervals described above
selection_intvl = cumsum(valid_aa); % Cumulative sum
% Then, we have to normalize the interval to a [0,1]. We do this by
% division by the last element.
selection_intvl = selection_intvl./selection_intvl(end);
% What is left to do, is to pick the segment as described above, and
% execute the according reaction channel's change. The values in
% selection_intvl refer to the upper ends of the segments of the
% interval. So, what we really need to find is the segment, which has
% the smallest upper end that is greater than the random number.
% The following operation is simply an implementation of what I
% described above. It will probably not be immediately clear. Go to the
% the MatLab command line, and type 'help find', that should explain to
% you how this works in detail.
selected_ind = find(selection_intvl>rand_nums(2,R_counter),1,'first');
selected_change = valid_changes(selected_ind);
% To finish the loop, we apply the time and state change we
% stochastically determined above, and store the current state of the
% simulation into the pre-allocated storage variables.
tt = tt + deltat_t;
nn = nn + selected_change;
tt_cont(R_counter+1) = tt;
nn_cont(R_counter+1) = nn;
end
% Let us see what the outcome of this is
plot(tt_cont,nn_cont,'k-')
xlabel('t')
ylabel('n')
Code is from:
|
Saturday, December 12, 2015
Gillespie, MATLAB sample code
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment