Modeling of signaling pathways is an important part of cancer research, as it is essential to understand how proteins interact with each other to provide or impair a specific cell function. There are two main approaches to tackle that problem mathematically: 1) deterministic approach, in which we assume that the number of molecules is large and stochastic effects are negligible, 2) stochastic approach (typically using Gillespie algorithm), in which we model exact number of molecules and take into account inherent system noise. Today I’ll focus on the stochastic approach. You can find online plenty of packages to perform the stochastic simulations using Gillespie algorithm, but I haven’t seen any that is harnessing the vectorization in MATLAB. Today I will show that using vectorization we can increase the speed of the simulations about 100 fold. Such an increase in computational speed is very important, as more and more frequently agent based models incorporate signaling pathways on the cell level.
Let us start with basic information about the Gillespie algorithm. Using the algorithm we describe the change in the number of molecules from K different species (x1(t),…,xK(t)) that at each moment can undergo one of L different biochemical reactions. Each reaction has a specific propensity function, e.g. biding of two molecules from different species can have propensity r*x1(t)*x2(t) where r is biding rate. The propensities determine the time to next reaction (exponentially distributed) and which reaction occurs next (discrete distribution with propensity based probabilites). After a reaction occurs the system is updated according to stoichiometric matrix, i.e. matrix describing how many reactants are removed from the system and what products are created.
As a working example we will consider a very basic gene expression pathway – we will describe the number of mRNAs and corresponding proteins. Four considered reactions will be: transcription, translation, mRNA degradation, and protein degradation. Ok, let us start coding by defining considered pathway in the main script.
%DEFINING REACTION RATES
par.kr = 1/10; %transcription rate
par.kp = 10; %translation rate
par.gr = 1/150; %mRNA degradation rate
par.gp = 1/30; %protein degradation rate
%DEFINING PROPENSITY FUNCTIONS AS A FUNCTION HANDLE
prop = @(x,par)([par.kr,... %transcription, one mRNA molecule created
par.kp*x(1),... %translation, one protein created
par.gr*x(1),... %mRNA degradation, one mRNA molecule removed
par.gp*x(2)]); %protein degradation, one protein molecule removed
%DEFINING INITIAL CONDITION, order [mRNA, Protein]
init = [0;1];
%DEFINING STOICHIOMETRIC MATRIX
%column corresponds to the reaction, row corresponds to the molecule
%order as in prop and init variables
stoch = [1 0 -1 0;... 0 1 0 -1];
%DEFINING TIME MESH FOR THE OUTPUT TRAJECTORY
tmesh = linspace(0,1000,100) ;
Now in a few lies we can code the whole Gillespie algorithm in a separate function file (without using any vectorization techniques first).
function trajectory = SSA(tmesh, par,prop,stoch, init )
%tmesh - time mesh on which solution should be returned
%per - parameters of the pathway
%prop - definition of propensity functions
%stoch - stochiometric matrix
%init - initial condition for the pathway
t = 0; %current time
state = init(:); %variable with current system state
trajectory = zeros(length(init),length(tmesh)); %preparing output trajectory
trajectory(:,1) = init(:); %setting initial value as the first element in trajectory
cindx = 2; %current trajectory index
N = length(tmesh); %number of time points
while t<tmesh(end)
Q = feval(prop,state,par); %calculating propensities of the reactions
Qs = sum(Q); %total propensity
dt = -log(rand())/Qs; %generating time to the next reaction
R = sum(rand >= cumsum([0 Q])/Qs); %selecting reaction
state = state + stoch(:,R); %updating state
t = t + dt; %updating time
%writing the output
while cindx<=N && t>tmesh(cindx)
trajectory(:,cindx) = state;
cindx = cindx+1;
end
end
end
and add the following lines in the main script to simulate and plot trajectory.
%simulating
traj = SSA(tmesh, par,prop,stoch, init );
%plotting
figure(1)
plot(tmesh,traj(1,:))
xlabel('Time, min')
ylabel('#mRNAs')
figure(2)
plot(tmesh,traj(2,:),'r')
xlabel('Time, min')
ylabel('#Proteins')
Below are the resulting plots after one stochastic realization of the whole pathway.
So now we can simulate our simple pathway, but there is one problem: it takes about 3.5 seconds to generate trajectory presented in the above plot. If we want to simulate that pathway for each agent present in the system and we reach the level of 2000 agents, then simulating the pathway for each of the agents using simple for loop will take about 2 hours… What can we do to decrease computational time? We could of course use multiple CPUs – using 100 of them will allow to do the computations in about a minute. But who has access to 100 CPUs? Only few of us. What about using the vectorization option in MATLAB and modifying the SSA function to perform those 2000 simulations simultaneously? That is definitely worth trying. We need only to add additional argument to SSAv function – number of trajectories to generate, and modify few lines to do vectorized calculations:
function trajectory = SSAv(tmesh, par,prop,stoch, init, nSim )
%tmesh - time mesh on which solution should be returned
%par - parameters of the pathway
%prop - definition of propensity functions
%stoch - stochiometric matrix
%init - initial condition for the pathway
%nSim - number of simulations to perform
tmesh = tmesh(:); %reshaping mesh to be vertical vector
t = zeros(nSim,1); %current time for each simulation
state = repmat(init(:)',nSim,1); %current state variable for each simulation
trajectory = zeros(nSim,length(init),length(tmesh)); %preparing output trajectory
trajectory(:,:,1) = state;%setting initial value as the first element in trajectory
cindx = 2*ones(nSim,1);%current trajectory indices
N = length(tmesh); %number of time points
aux = 1:nSim; %
while ~isempty(t)
Q = feval(prop,state,par); %calculating propensities of the reactions
Qs = sum(Q,2); %total propensities
dt = -log(rand(size(Qs,1),1))./Qs; %generating time to the next reaction
P = bsxfun(@rdivide, cumsum([zeros(size(Qs,1),1) Q],2),Qs); %probabilities for each reaction
R = sum(bsxfun(@ge,rand(size(Qs,1),1),P),2); %selecting reaction
state = state + stoch(:,R)'; %updating state
t = t + dt; %updating time
%writing the output
update = t > tmesh(cindx);
while any(update)
%updating state
iupdt = find(update);
for i = 1:length(iupdt)
trajectory(aux(iupdt(i)),:,cindx(iupdt(i))) = state(iupdt(i),:);
end
cindx = cindx+update;
%removing finished simulations from the variables
indx = cindx > N;
if any(indx)
cindx(indx) = [];
t(indx) = [];
aux(indx) = [];
state(indx,:) = [];
if isempty(cindx)
break;
end
end
update = t > tmesh(cindx);
end
end
end
We need also to make sure that in the main script that the propensity functions are defined in the vectorized manner (everything else remains the same).
%DEFINING PROPENSITY FUNCTIONS AS A FUNCTION HANDLE IN VECTORIZED MANNER
prop = @(x,par)([par.kr*ones(size(x,1),1),...
par.kp*x(:,1),...
par.gr*x(:,1),...
par.gp*x(:,2)]);
Below are the plots of trajectories generated using the vectorized function for nSim = 100.
Now the important question is how much time did it take to simulate those 100 trajectories. We know that using the simple non-vectorized code and a for loop it would take about 360 seconds. The vectorized function did it in 16 seconds! That is 22 fold decrease in computational time. Plot below shows how computational time changes with the number of trajectories to simulate.
Vectorized code generated 2000 trajectories in 62 seconds! So we don’t need those 100 CPUs to obtain reasonable time of computations.
I hope that this post will convince some people that MATLAB is not bad after all.
COMPUTE CANCER 