Revision 8af1e7fcf1fa126908465558c4ef6d127000a9c8 authored by AlistairBoettiger on 27 April 2021, 19:02:44 UTC, committed by AlistairBoettiger on 27 April 2021, 19:02:44 UTC
1 parent 191524d
Fig5_PTM_TxVsLoopRateStartOn.m
% Stochastic model, Fig 4A
%
%
% Corresponding ODE model:
% da/dt = (Et*(k1*k2*r2*a.^2 + k1*r1*tau2*a + r0*tau1*tau2))./(k1*k2*a.^2 + k1*tau2*a + tau1*tau2) +loopRate;
% Model Parameters
k1 = .003; % Rate of conversion from state 0 to state 1.
k2 = .5; % Rate of conversion from state 1 to state 2.
tau1 = 20; % Rate of conversion from state 1 to state 0.
tau2 = .05; % Rate of conversion from state 2 to state 1.
hdac = 2; % Rate of tag removal.
r0 = .01; % Basal tagging rate (state 0)
r1 = .5; % Stimulated tagging rate (state 1)
r2 =5.0; % Maximally stimulated tagging rate (state 2)
loopRate = 9; % .65 ;% .95; % the linear shift to the sigmoid. % 20% change, 3.5 to 4.5
Et= 20; % Maximum tag-transferases associated with promoter
% initial conditions
ac_init = 40; % No tags (40 = start high)
Pr_init = 3; % all enzymes in the unstimulated state 0
% E-P loop rates explored
loopRates = 0:.2:20;
% number of simulations, length of simulation, time scale of sampling
tSteps = 3e5; % length of simulation
nRuns = 3e3; % number of replicates
t_scale = 1000; % time scale of sampling
numLoops = length(loopRates);
%% Run the simulation
saveFolder = SetFigureSavePath('U:\Data\Analysis\2020-04-23_JordanLoopSims\');
saveName = [saveFolder,'loopRatesVsExpStoch_startON.mat'];
% only need to repeat the simulation if a saved file is not found
if exist(saveName,'file') == false
acVals3 = cell(numLoops,1);
for r=1:numLoops
disp(r/numLoops);
loopRate = loopRates(r);
acR = [r0 r1 r2]; % acetylation rates for the different enhancer states
compTx = false;
nZymes =Et;
clear ac_values;
ac_values = zeros(nRuns,tSteps,'uint8');
parfor n=1:nRuns
acR = [r0 r1 r2];
ac = ac_init; % 0 ... inf
Pr = Pr_init*ones(1,nZymes); % 1, 2, 3
for t=1:tSteps % straight discreet time simulation
% enhancer mediated actyl
if rand*t_scale < loopRate
ac = ac + 1;
end
% lose actyl
if (rand*t_scale < hdac*ac) && ac > 1
ac = ac-1;
end
% handle promoter state
for e=1:nZymes
% gain actyl
if rand*t_scale < acR(Pr(e)) % rate depends on enzyme state
ac = ac + 1;
end
if Pr(e)==1
if rand*t_scale < ac*k1
Pr(e)=2;
end
elseif Pr(e)==2
if rand*t_scale < ac*k2
Pr(e) = Pr(e)+1;
end
if rand*t_scale < tau1
Pr(e) = Pr(e)-1;
end
elseif Pr(e) == 3
if rand*t_scale < tau2
Pr(e) = Pr(e) - 1;
end
end
end
ac_values(n,t) = ac;
end
end
acVals3{r} = ac_values;
acv = ac_values(:,end-500:end);
end
%% Save or load results
save(saveName,'acVals3','expR','-v7.3');
else
load(saveName,'acVals3');
end
%% plot stoch results
% reorganizing data into matrix for plotting
numLoops = length(acVals3);
[nRuns, tSteps] = size(acVals3{1});
highDist = zeros(numLoops,nRuns);
for r=1:numLoops
highDist(r,:) = acVals3{r}(:,end);
end
% plot results as a violin plot
f3 = figure(3); clf;
subsample = 1:5:numLoops-10;
violin(highDist(subsample,:)','bandwidth',.5,'plotMean',false,'faceColor',[1 .5 .5],'alpha',1); hold on;
highMed = median(highDist(subsample,:),2);
plot(highMed,'r'); hold on;
ylim([0,70]);
set(gca,'xTickLabel',loopRates(subsample));
set(gcf,'color','w');
xlabel('loop rate');
ylabel('transcription rate');
%%
% SaveFigure(f3,'name','violin_stochBistable_low','formats',{'png','eps'});
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...