Skip to content

Commit 4440867

Browse files
dmitrisvetlovdmitrisvetlov
committed
Add initial draft of (in)direct Gillespie.
1 parent 5a56454 commit 4440867

File tree

2 files changed

+118
-39
lines changed

2 files changed

+118
-39
lines changed

‎src/+ssit/+ssa/runSingleSsa.m

+108-32
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,27 @@
1-
function [X_array] = runSingleSsa(x0,S,W,T_array,isTimeVarying,signalUpdateRate,parameters)
1+
function [X_array] = runSingleSsa(x0, S, W, T_array, isTimeVarying, ...
2+
signalUpdateRate, delayedReactions, parameters)
3+
24
% Start the simulation.
5+
6+
%% Step 1: Initialization
7+
% Input values for the initial state from the provided vector, set time t
8+
% to be the initial time of simulation, and set the reaction counter i = 1.
9+
310
t = min(T_array); % initial time of simulation.
411
x = x0;
512
iprint = 1; % The next time at which to record results.
613
Nsp = size(S,1); %Number of species.
714
Nt = length(T_array);
815
X_array = zeros(Nsp,Nt);
916
props = [];
17+
18+
hasDelayedReactions = ~isempty(delayedReactions);
19+
20+
% To avoid hints/warnings, we will initialize an array that will hold a
21+
% record of delayed reactions: the first column will contain the scheduled
22+
% time; the second, the reaction number.
23+
scheduledDelayedReactions = [NaN, NaN];
24+
1025
if isTimeVarying
1126
S = [zeros(Nsp,1),S];
1227
props(1) = signalUpdateRate;
@@ -17,49 +32,110 @@
1732

1833
isAFun = isa(W{1}, 'function_handle');
1934

20-
while t<max(T_array)
21-
22-
%% Choose time of reaction
35+
while t < max(T_array)
36+
%% Step 2: Compute propensities for all reactions
2337
if ~isAFun
24-
%% Command Line Code Approach
38+
% Command Line Code Approach
2539
WT = W{1}.hybridFactorVector(t,parameters);
2640
for i=length(W):-1:1
41+
% Evaluate the propensity functions at the current state:
2742
if ~W{i}.isTimeDependent||W{i}.isFactorizable
28-
props(i+jt) = WT(i)*W{i}.stateDependentFactor(x,parameters); % evaluate the propensity functions at the current state.
43+
props(i+jt) = WT(i)*W{i}.stateDependentFactor(x,parameters);
2944
else
30-
props(i+jt) = W{i}.jointDependentFactor(t,x,parameters); % evaluate the propensity functions at the current state.
45+
props(i+jt) = W{i}.jointDependentFactor(t,x,parameters);
3146
end
3247
end
3348
else
34-
%% GUI approach
49+
% GUI approach
3550
for i=length(W):-1:1
36-
props(i+jt) = W{i}(x,t); % evaluate the propensity functions at the current state.
51+
% Evaluate the propensity functions at the current state.
52+
props(i+jt) = W{i}(x,t);
3753
end
3854
end
39-
w0 = sum(props); % sum the propensity functions (inverse of ave. waiting time).
40-
tau = -1/w0*log(rand); % The time until the next reaction.
41-
42-
%% update time
43-
t = t+tau; % The time of the next reaction.
4455

45-
while iprint<=Nt&&t>T_array(iprint)
46-
X_array(:,iprint) = x;
47-
iprint=iprint+1;
48-
end
49-
50-
if t>max(T_array)
51-
break;
52-
end
53-
54-
%% Choose which reaction.
55-
r2 = w0*rand; % this is a uniform r.v. betweeen zero and w0;
56-
j = 1;
57-
while sum(props(1:j))<r2
58-
j=j+1;
56+
%% Step 3: Generate uniform random number(s) in [0, 1].
57+
u1 = rand;
58+
u2 = rand;
59+
60+
%% Step 4: Compute the time interval until the next reaction:
61+
% Note that the sum of the propensities is the inverse of the average
62+
% waiting time.
63+
64+
delta_t = -1 * log(u1) / sum(props);
65+
66+
if (t + delta_t) > max(T_array)
67+
break % Do not simulate beyond the time boundaries.
5968
end
60-
% At this point j is the chosen reaction.
69+
70+
%% Step 5: Are there upcoming delayed reactions in [t, t + delta_t]?
6171

62-
%% Update state
63-
x = x + S(:,j);
64-
end
72+
scheduledDelayedReactions = sortrows(scheduledDelayedReactions);
73+
delayedReactionTimes = scheduledDelayedReactions(:, 1);
74+
rows = find(...
75+
delayedReactionTimes >= t & delayedReactionTimes <= (t + delta_t));
76+
rowsize = size(rows);
77+
if rowsize(1, 1) > 0
78+
% There IS a delayed reaction in that time interval.
79+
% Advance t to the time of the next scheduled delayed reaction.
80+
81+
t = scheduledDelayedReactions(rows, 1);
82+
rxn = scheduledDelayedReactions(rows, 2);
83+
84+
% Add the current state to the record of states for all intervening
85+
% times:
86+
while (iprint <= Nt) & (t > T_array(iprint))
87+
X_array(:, iprint) = x;
88+
iprint = iprint + 1;
89+
end
90+
91+
% Update the state according to the delayed reaction:
92+
93+
x = x + S(:, rxn);
94+
else % rowsize(1, 1) > 0
95+
%% Step 6: Find the channel of the next reaction:
96+
97+
threshold = u2 * sum(props); % threshold is in [0, sum(props)]
98+
rxn = 1;
99+
while sum(props(1: rxn)) < threshold
100+
rxn = rxn + 1;
101+
end
102+
% At this point, rxn is the number of the chosen reaction.
103+
104+
%% Step 7: Determine whether the next reaction is delayed:
105+
106+
rxnIsDelayed = false;
107+
rxnDelayIndex = 0;
108+
if hasDelayedReactions
109+
rxnDelayIndex = find(delayedReactions(:, 1) == rxn, 1);
110+
rxnIsDelayed = ~isempty(rxnDelayIndex);
111+
end
112+
113+
if rxnIsDelayed
114+
% If the reaction is delayed, postpone the update to
115+
% t_d = t + tau, where tau is the delay corresponding to the
116+
% reaction.
117+
118+
nextRxn = [t + delayedReactions(rxnDelayIndex, 1), rxn];
119+
scheduledDelayedReactions = ...
120+
[scheduledDelayedReactions; nextRxn];
121+
else
122+
% If the reaction is not delayed, perform it; i.e.
123+
124+
% Add the current state to the record of states for all
125+
% intervening times:
126+
127+
while (iprint <= Nt) & (t > T_array(iprint))
128+
X_array(:, iprint) = x;
129+
iprint = iprint + 1;
130+
end
131+
132+
% Update the state according to the reaction:
133+
134+
x = x + S(:, rxn);
135+
136+
% Advance t to the time of the reaction:
65137

138+
t = t + delta_t;
139+
end % Step 7 ...
140+
end % rowsize(1, 1) == 0 [No delayed reactions were upcoming]
141+
end % while t < max(T_array) ...

‎src/CommandLine/SSIT.m

+10-7
Original file line numberDiff line numberDiff line change
@@ -16,9 +16,11 @@
1616
'constantJacobianTime',1.1); % Options for FSP solver.
1717
sensOptions = struct('solutionMethod','forward','useParallel',true);
1818
% Options for FSP-Sensitivity solver.
19-
ssaOptions = struct('Nexp',1,'nSimsPerExpt',100,'useTimeVar',false,...
20-
'signalUpdateRate',[],'useParalel',false,...
21-
'verbose',false); % Options for SSA solver
19+
ssaOptions = struct('Nexp', 1, 'nSimsPerExpt', 100, ...
20+
'isTimeVarying', false, ...
21+
'signalUpdateRate', [], 'useParallel', false, ...
22+
'isDelayed', false, ...
23+
'verbose', false); % Options for SSA solver
2224
pdoOptions = struct('unobservedSpecies',[],'PDO',[]);
2325
% Options for FIM analyses
2426
fittingOptions = struct('modelVarsToFit','all','pdoVarsToFit',[],...
@@ -34,7 +36,7 @@
3436
useHybrid = false
3537
hybridOptions = struct('upstreamODEs',[]);
3638
propensitiesGeneral = [];% Processed propensity functions for use in solvers
37-
39+
delayedReactions = [];
3840
end
3941

4042
properties (Dependent)
@@ -897,9 +899,9 @@ function summarizeModel(obj)
897899
case 'SSA'
898900
Solution.T_array = obj.tSpan;
899901
Nt = length(Solution.T_array);
900-
nSims = obj.ssaOptions.Nexp*obj.ssaOptions.nSimsPerExpt*Nt;
902+
nSims = obj.ssaOptions.Nexp * obj.ssaOptions.nSimsPerExpt * Nt;
901903
W = obj.propensitiesGeneral;
902-
if obj.ssaOptions.useParalel
904+
if obj.ssaOptions.useParallel
903905
trajs = zeros(length(obj.species),...
904906
length(obj.tSpan),nSims);% Creates an empty Trajectories matrix from the size of the time array and number of simulations
905907
parfor isim = 1:nSims
@@ -910,8 +912,9 @@ function summarizeModel(obj)
910912
obj.stoichiometry,...
911913
W,...
912914
obj.tSpan,...
913-
obj.ssaOptions.useTimeVar,...
915+
obj.ssaOptions.isTimeVarying,...
914916
obj.ssaOptions.signalUpdateRate,...
917+
obj.delayedReactions, ...
915918
[obj.parameters{:,2}]');
916919
end
917920
Solution.trajs = trajs;

0 commit comments

Comments
 (0)