% % Likelihood_path.m % % Run Gillespie on the model % % A -> B -> C % % with rates k1 = theta, k2 = 1 so as to compute the % derivative of E C(tend) with respect to theta. % % usage: Likelihood_path(tend,initial,theta) % % tend = max time for the program to run % % initial = initial condition, ex: [30, 0, 0] % % theta = value of parameter at which we are computing the % derivative. % % example: Likelihood_path(3,[300, 0, 0],0.5) % function [x,M] = Likelihood_path(tend,initial,theta) %rate constants k = [theta 1]; %maximum number of steps to allow. maxi = 10^4; %reaction vectors zeta = [-1 0; 1 -1; 0 1]; %"M" stores the weights used to compute the derivatives. M = 0; %initialize. T represents time and x is the state of the system. T = 0; x = initial; %intensity functions lambda = zeros([2,1]); %main loop for i=1:maxi-1, %update intensity functions lambda(1) = k(1)*x(1); lambda(2) = k(2)*x(2); lambda0 = lambda(1) + lambda(2); %calculate when the next reaction will take place. r = rand; if lambda0 == 0 break else delta = log(1/r)/lambda0; end %update time T = T + delta; %find which reaction fired. r = rand; if r < lambda(1)/lambda0; loc = 1; else loc = 2; end %update x x = x + zeta(:,loc); %update M if loc == 1 M = M + (1/theta)*(1 - lambda(1)*delta); else M = M - (1/theta)*lambda(1)*delta; end if T + delta > tend break end end end %The program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%