Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Tree: ad43822291
Fetching contributors…

Cannot retrieve contributors at this time

172 lines (126 sloc) 4.754 kB
%% Plot_model_pdf.m
%
% Alistair Boettiger Date Begun: 04/06/09
% Levine Lab Functional Since: 10/09/09
% Functional computational code Last Modified: 11/07/09
%% Notes:
% compute moment generating function of first passage times for arbitrary
% pinch point decomposed (series) Markov Chains using Matlab's symbolic
% Algebra routines
% Compute laplace transform of complete probability distribution
%% Requires functions
% SeriesDecomp_pdist.m
%
clear all;
%
% Specific Model Constants
kab = sym('kab','real');
kba = sym('kba','real');
k12 = sym('k12','real');
k21 = sym('k21','real');
k23 = sym('k23','real');
k24 = sym('k24','real');
k32 = sym('k32','real');
k35 = sym('k35','real');
k53 = sym('k53','real');
k54 = sym('k54','real');
k42 = sym('k42','real');
k45 = sym('k45','real');
k56 = sym('k56','real');
k65 = sym('k65','real');
k67 = sym('k67','real');
k78 = sym('k78','real');
% Submatrices for initiation regulated system
GI{1} = zeros(3);
GI{2}=[[-kab,kab,0];
[kba, -kba-k12, k12];
[0, k21, -k21]];
GI{3}=[[-k23-k24, k23, k24, 0];
[k32, -k32-k35, 0, k35];
[k42, 0, -k42-k45, k45];
[0, k53, k54, -k53-k54]];
GI{4} =[[-k56, k56, 0, 0];
[k65, -k65-k67, k67, 0];
[0, 0, -k78, k78];
[0, 0, 0, 0]];
[vI] = SeriesDecomp_pdist(GI);
% Submatrices for elongation regulated system
p = kba/(kab+kba);
GE{1} = zeros(3);
GE{2}=[[-k12, k12];
[k21, -k21]];
GE{3}=[[-k23-k24, k23, k24, 0];
[k32, -k32-k35, 0, k35];
[k42, 0, -k42-k45, k45];
[0, k53, k54, -k53-k54]];
% 5 6 7A 7B 8
GE{4} = [[-k56, k56, 0, 0 ,0];
[k65,-k65-p*k67-(1-p)*k67, p*k67, (1-p)*k67, 0];
[0,0,-kab,kab,0 ];
[0,0,0,-k78,k78 ];
[0,0,0,0,0]];
% either the enhancer has already enabled the promoter to fire, and all
% gates are open, or the enhancer opens the gate after the chain
% reaches the paused state.
[vE] = SeriesDecomp_pdist(GE);
% % 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
% vars = [k12,k21,k23,k24,k32,k35,k53,k54,k42,k45,k56,k65,k67,k78,kab,kba];
save Model2b_pdist
clear all;
load Model2b_pdist;
% Specific Model Constants
kab = sym('kab','real');
kba = sym('kba','real');
k12 = sym('k12','real');
k21 = sym('k21','real');
k23 = sym('k23','real');
k24 = sym('k24','real');
k32 = sym('k32','real');
k35 = sym('k35','real');
k53 = sym('k53','real');
k54 = sym('k54','real');
k42 = sym('k42','real');
k45 = sym('k45','real');
k56 = sym('k56','real');
k65 = sym('k65','real');
k67 = sym('k67','real');
k78 = sym('k78','real');
lambda = sym('lambda','real');
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
pars = [ k12, k21, k23, k24, k32, k35, k53, k54, k42, k45, k56, k65, k67, k78, kab, kba];
vals = 10*[6*.0216, 6*.145, 6*.0216, 6*.0216, 6*.145, 6*.0216 6*.145, 6*.145, 6*.0216, 6*.0216, 3*.00159, .001, 3*.00159, 3*.00159,1E-4, 5E-5, ];
% vals = rand(16,1);
% pars = [ kab, kba, k12, k21, k23, k24, k32, k35, k42, k45, k56, k65 ,k67, k78];
% vals = 3*[1E-3, 3E-1, 6*.0216, 6*.145, 6*.0216, 6*.0216, 6*.145, 6*.0216 6*.145, 6*.0216, 3*.00159, .001, 3*.00159, 3*.00159];
fI = subs(vI,pars,vals);
fE = subs(vE,pars,vals);
fIs = char(fI(1));
fIs = strrep(fIs,'*','.*');
fIs = strrep(fIs,'/','./');
fIs = strrep(fIs,'^','.^');
lambda = linspace(1E-8,4E2,1000);
fIv = eval(fIs);
fIl = fIv(fIv>0);
% Here we actually numerically invert the Laplace space solutions of the
% model.
lam = linspace(10,3E3,1000);
FI = invlap2(fI(1), lam');
FE = invlap2(fE(1), lam');
% Check efficiency of inversion by confirming pdf integrates to 1.
dt = (max(lam)-min(lam)) / length(lam) ;
norm_E = sum(FE*dt) ; % integrate p(t)*dt
norm_I = sum(FI*dt) ;
mean_E = lam*FE*dt/60 ;% integrate t*p(t)*dt
mean_I = lam*FI*dt/60;
std_E = sqrt(lam.^2*FE*dt - (lam*FE*dt)^2)/60;
std_I = sqrt(lam.^2*FI*dt - (lam*FI*dt)^2)/60;
save compfull_pdist2_data;
% Plot results
figure(1); clf; plot(lam/60,FI,'b','LineWidth',2); hold on;
plot(lam/60,FE,'r','LineWidth',2); xlim([0,lam(end)/60]); %ylim([0,1.2E-3]);
title(['\mu_{ER} = ',num2str(mean_E,2), ' \sigma_{ER} = ', num2str(std_E,2), ...
' \mu_{IR} = ',num2str(mean_I,2), ' \sigma_{IR} = ', num2str(std_I,2),...
' minutes'],'FontSize',15);
xlabel('time (minutes)','FontSize',15);
set(gca,'FontSize',15); set(gcf,'color','w');
legend('IR Model','ER Model')
Jump to Line
Something went wrong with that request. Please try again.