%% Condensate looping model
% This version has simulations both from start on and start off
% the version in Fig3 is the start-off case. The longer the simulation
% runs, the closer these two converge. The hysteresis due to starting
% condition is explored in more detail in figure 5.
% size of the simulation
tSamples = 140;
nCells = 1000; % we'll track one particular promoter in N isolated stochastic instances, such as in N cells
nR = 21;
ti = 50; %start period
% model parameters
loopRates = linspace(0,.2,nR); %
addPol = .01; % can dedimensionalize on this one
losePol = .25;
tSteps = 140;
clusterMax = 6; % max cluster size
% if total PolII is conserved, and bound to diverse sites throughout the
% cell, there is an amount P bound at any given time, distributed across n
% spatially segregated promoters with an average of P/n. If the number of
% total promoters is large, our promoter will rarely have more the average
% available to it.
% A more complete model would have this as a drawn from a random
% distribution of available cluster sizes, which serves only to add spread
% to the on values (as would stochastic transcription models).
%% start ON
startOn = 1;
promoterPolLog1 = zeros(tSamples,nCells,nR);
for r=1:nR
promoterPol = startOn*clusterMax*ones(1,nCells);
e= loopRates(r); % .01; % looks promising for bistable in e
tt=0;
for t=1:tSteps
% promoter gains a Pol molecule
stoch = rand(nCells,clusterMax+1) < addPol + e;
% its actually faster in Matlab to generate more random numbers
% than we need (here for Nmax encounters), which saves us having to
% loop over encouters (even though many loops could exit before
% Nmax and need fewer total random number draws).
for c=1:nCells % would be better to do this without a loop
stoch(c,promoterPol(c)+2:end) = 0;
end
promoterPol(any(stoch,2)) = promoterPol(any(stoch,2))+1;
promoterPol(promoterPol>clusterMax) = clusterMax;
% promter loses a pol molecule
stoch = rand(1,nCells) < losePol;
promoterPol(stoch) = promoterPol(stoch) - 1;
promoterPol(promoterPol<0) = 0;
if rem(t,round(tSteps/tSamples))==0
tt=tt+1;
promoterPolLog1(tt,:,r) = promoterPol;
end
end
figure(1); clf;
imagesc(promoterPolLog1(:,:,r));
colorbar; colormap('default');
end
timeAveStON = squeeze(mean(promoterPolLog1(ti:end,:,:),1));
figure(1); clf; imagesc(timeAveStON);
popAveStON = mean(timeAveStON,1);
figure(2); clf; plot(loopRates,popAveStON,'r');
%% start off
startOn = 0;
promoterPolLog0 = zeros(tSamples,nCells,nR);
for r=1:nR
promoterPol = startOn*clusterMax*ones(1,nCells);
e= loopRates(r); % .01; % looks promising for bistable in e
tt=0;
for t=1:tSteps
stoch = rand(nCells,clusterMax+1) < addPol + e; % figure(1); clf; imagesc(stoch);
for c=1:nCells % would be better to do this without a loop
stoch(c,promoterPol(c)+2:end) = 0;
end
promoterPol(any(stoch,2)) = promoterPol(any(stoch,2))+1;
promoterPol(promoterPol>clusterMax) = clusterMax;
% promter loses pol
stoch = rand(1,nCells) < losePol;
promoterPol(stoch) = promoterPol(stoch) - 1;
promoterPol(promoterPol<0) = 0;
if rem(t,round(tSteps/tSamples))==0
tt=tt+1;
promoterPolLog0(tt,:,r) = promoterPol;
end
end
figure(1); clf;
imagesc(promoterPolLog0(:,:,r));
colorbar; colormap('default');
end
timeAveStOFF = squeeze(mean(promoterPolLog0(ti:end,:,:),1));
figure(1); clf; imagesc(timeAveStOFF);
popAveStOFF = mean(timeAveStOFF,1);
figure(2); hold on; plot(loopRates,popAveStOFF,'b');
%%
f3 = figure(3); clf;
plot(popAveStON,'r'); hold on;
plot(popAveStOFF,'b'); hold on;
violin(timeAveStON,'bandwidth',.02,'plotMean',false,'faceColor',[1 0 0],'alpha',.5); hold on;
violin(timeAveStOFF,'bandwidth',.02,'plotMean',false,'faceColor',[0 0 1],'alpha',.5); hold on;
set(gca,'XTick',1:2:nR,'XTickLabel',num2str(loopRates(1:2:nR)',2));
ylim(round([0,1.05*clusterMax]));
set(gcf,'color','w');
%%
savePath = SetFigureSavePath('C:\Shared\Documents\Jordan Looping Model\Revision1\Images\');
f3 = figure(3); clf;
plot(popAveStON,'r'); hold on;
plot(popAveStOFF,'b'); hold on;
violin(timeAveStON,'bandwidth',.02,'plotMean',false,'faceColor',[1 0 0],'alpha',1); hold on;
violin(timeAveStOFF,'bandwidth',.02,'plotMean',false,'faceColor',[0 0 1],'alpha',1); hold on;
set(gca,'XTick',1:2:nR,'XTickLabel',num2str(loopRates(1:2:nR)',2));
ylim(round([0,1.05*clusterMax]));
xlabel('loop rate');
ylabel('Promoter associated PolII');
set(gcf,'color','w');