%% lifetimes_simulator.m
% Written for Berro Lab at Yale University
% Last edited by Mike Lacy, March 2019

% This script runs a set of simulations of accumulation and disassembly
% dynamics, as a simple illustration of the logic of molecular assembly and
% turnover. Binding events are represented by a randomly generated list of
% arrival times, and unbinding events are represented by another list of
% departure times. Depending on the logical rules of a model you wish to
% illustrate, these lists can be calculated and sorted in different ways.
%
% The script then computes and plots some features that resemble measured
% data, such as the number of "molecules" at the site (how many objects
% have already arrive but have not yet departed), the distribution of 
% "residence times" (the difference between a departure time and arrival
% time for a given event), and the relation between the number of
% "molecules" consumed relative to the maximum number present (the peak #
% height).

%% Initialize parameters

molecules = 900;        %the number of binding/unbinding events
patchLifetime = 20;     %total timespan
nruns = 50;             %simulator iterations

% initialize arrays for: arrival times (*a), departure times (*d), current
% number of molecules (*c), and lifetimes (*lt) for each selected model:
% random assembly/random disassembly with Gaussian profile (randN*), random
% assembly with first-in/first-out disassembly (randF*), and turnover (toN*)

randNa = zeros(molecules,nruns);     %arrival times
randNd = zeros(molecules,nruns);     %departure times
randNc = zeros(patchLifetime,nruns); %current number of molecules at each time point
randNlt = zeros(molecules,nruns);    %lifetimes of molecules
randFa = zeros(molecules,nruns);
randFd = zeros(molecules,nruns);
randFc = zeros(patchLifetime,nruns);
randFlt = zeros(molecules,nruns);
toNa = zeros(molecules,nruns);
toNd = zeros(molecules,nruns);
toNc = zeros(patchLifetime,nruns);
toNlt = zeros(molecules,nruns);

for iter = 1:nruns
    [randNa(:,iter), randNd(:,iter), randNc(:,iter), randNlt(:,iter)] = ltSimulator(molecules, patchLifetime, 'Normal');
    [randFa(:,iter), randFd(:,iter), randFc(:,iter), randFlt(:,iter)] = ltSimulator(molecules, patchLifetime, 'FIFO');
    [toNa(:,iter), toNd(:,iter), toNc(:,iter), toNlt(:,iter)] = ltSimulator(molecules, patchLifetime, 'Turnover');
end
% see function definition at end of script.

%% combined plots
figure
% First row: plot binding and unbinding event frequencies
subplot(3,3,1)
hold on;
histogram(randNa(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 1],'LineWidth',1);
histogram(randNd(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [1 0 0],'LineWidth',1);
ylabel('Frequency');
xlabel('Time, sec');
xlim([0 patchLifetime]);
ylim([0 0.06]);
title(sprintf('Random assembly\n and disassembly'));

subplot(3,3,2)
hold on;
histogram(randFa(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 1],'LineWidth',1);
histogram(randFd(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [1 0 0],'LineWidth',1);
xlabel('Time, sec');
xlim([0 patchLifetime]);
ylim([0 0.06]);
title(sprintf('Random assembly,\n first-in-first-out'));

subplot(3,3,3)
hold on;
histogram(toNa(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 1],'LineWidth',1);
histogram(toNd(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [1 0 0],'LineWidth',1);
xlabel('Time, sec');
xlim([0 patchLifetime]);
ylim([0 0.06]);
legend('Binding','Unbinding');
title(sprintf('Continuous\n turnover'));

% Second row: plot number of molecules present over time
subplot(3,3,4)
hold on
plot(1:length(randNc(:,1)), mean(randNc,2)./max(mean(randNc,2)),'-','LineWidth',1,'Color',[0 0 0]);
% plot(1:length(randNc(:,1)), randNc(:,2),'-s','LineWidth',1,'Color',[0 0.5 0]);
ylabel(sprintf('Number of molecules,\nnormalized'));
xlabel('Time, sec');

subplot(3,3,5)
hold on
plot(1:length(randFc(:,1)), mean(randFc,2)./max(mean(randFc,2)),'-','LineWidth',1,'Color',[0 0 0]);
% plot(1:length(randFc(:,1)), randFc(:,2),'-s','LineWidth',1,'Color',[0 0.5 0]);
xlabel('Time, sec');
% legend('average of 50 runs','single run','Location','north');

subplot(3,3,6)
hold on
plot(1:length(toNc(:,1)), mean(toNc,2)./max(mean(toNc,2)),'-','LineWidth',1,'Color',[0 0 0]);
% plot(1:length(toNc(:,1)), toNc(:,2),'-s','LineWidth',1,'Color',[0 0.5 0]);
xlabel('Time, sec');

% Third row: plot molecule residence time distribution
subplot(3,3,7)
hold on;
histogram(randNlt(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 0],'LineWidth',1);
ylabel('Frequency');
xlabel(sprintf('Residence time, sec\nMean = %.2g +/- %.2g s',mean(randNlt(:)), std(randNlt(:))));
ylim([0 0.05]);

subplot(3,3,8)
hold on;
histogram(randFlt(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 0],'LineWidth',1);
xlabel(sprintf('Residence time, sec\nMean = %.2g +/- %.2g s',mean(randFlt(:)), std(randFlt(:))));
ylim([0 0.5]);

subplot(3,3,9)
hold on;
histogram(toNlt(:),0:0.25:patchLifetime,'Normalization','probability','DisplayStyle','stairs', 'EdgeColor', [0 0 0],'LineWidth',1);
xlabel(sprintf('Residence time, sec\nMean = %.2g +/- %.2g s',mean(toNlt(:)), std(toNlt(:))));

f=gcf;
ax = findobj(f,'Type','Axes');
for i=1:numel(ax)
    ax(i).Box = 'on';
    ax(i).FontSize = 10;
end

% plot cumulative total molecules consumed in each model
figure
hold on
histogram(randNa(:),0:0.25:patchLifetime,'Normalization','cdf','DisplayStyle','stairs', 'EdgeColor', [0.75 0 0.25],'LineWidth',1);
histogram(randFa(:),0:0.25:patchLifetime,'Normalization','cdf','DisplayStyle','stairs', 'EdgeColor', [1 0 0.5],'LineWidth',1);
histogram(toNa(:),0:0.25:patchLifetime,'Normalization','cdf','DisplayStyle','stairs', 'EdgeColor', [0 0.75 0],'LineWidth',1);
xlabel('Time, sec')
ylabel('Total molecules consumed')
legend('random','fifo','continuous turnover','Location','southeast')
a=gca; a.FontSize = 12;
xlim([0 patchLifetime])
title(sprintf('Molecules consumed,\nrelative to total consumed'));

%get the CDF values and normalize them to the peak #molecules numbers
%instead of total consumed.
randNCDF = a.Children(3).Values;
randFCDF = a.Children(2).Values;
toCDF = a.Children(1).Values;
% multiply CDF by factor of total molecules over peak #
randNCDFnorm = randNCDF.*molecules./max(mean(randNc,2));
randFCDFnorm = randFCDF.*molecules./max(mean(randFc,2));
toCDFnorm = toCDF.*molecules./max(mean(toNc,2));

figure
hold on
plot(0:0.25:patchLifetime-0.25, randNCDFnorm,'LineWidth',1,'Color',[0.7 0 0.3]);
plot(0:0.25:patchLifetime-0.25, randFCDFnorm,'LineWidth',1,'Color',[1 0 0.5]);
plot(0:0.25:patchLifetime-0.25, toCDFnorm,'LineWidth',1,'Color',[0 0.7 0]);
ylabel(sprintf('Molecules consumed,\nnormalized to peak # present'));
xlabel('Time, sec');
title(sprintf('Molecules consumed,\nrelative to peak number present'));
a=gca;
a.FontSize =12;
a.Box = 'on';
legend('Random','First-in/First-out','Continuous turnover','Location','east');

 
%% simulations of simple assembly/disassembly models and their lifetime distributions
% random in/out (linear assembly/disassembly or normally-distributed)
% first-in-first-out
% continuous on/off with assembly off-switch (constant off-rate, on-rate is a switch)

function [arrives, departs, current, lifetimes] = ltSimulator(num, patchLT, model)
% Inputs
% num: number of molecules consumed, e.g. 900
% patchLT: overall lifetime of patch, sec, e.g. 20
% model: choose from 'Normal', 'FIFO' (first-in/first-out), or 'Turnover'

% parameters for single molecule residence time (for Turnover model only):
smrtShape = 2;  %gamma distribution shape parameter
smrtScale = 1;  %gamma distribution scale parameter

if(~contains({'Normal', 'FIFO', 'Turnover'},model))
    error('Invalid model input. Choose from ''Normal'', ''FIFO'', or ''Turnover''.');
end

if(strcmp(model, 'Uniform'))
    arrives = rand(num,1)*patchLT/2;    %here stdev is arbitrary choice, this result looks fine
    departs = rand(num,1)*patchLT/2 + patchLT/2;
else
    arrives = randn(num,1)*patchLT/10 + patchLT/4;  
    departs = randn(num,1)*patchLT/10 + 3*patchLT/4;
end

if(strcmp(model, 'FIFO'))
    arrives = sort(arrives);
    departs = sort(departs);
% arrival and departure are paired for a given event. The first arrival now
% has the same index as the first departure.
end

lifetimes = departs - arrives;

if(strcmp(model, 'Turnover'))
    % continuous turnover model. also requires different window for arrival times
    arrives = abs(randn(num,1)*patchLT/5 + patchLT/3);
    
    % lifetimes distribution is set by the disassembly rates
    lifetimes = gamrnd(smrtShape,smrtScale,num,1);
    % the departures are arrival + random selection from lifetimes distribution.
    departs = arrives + lifetimes;
end

current = zeros(patchLT,1);
for i=1:patchLT
    % valid for all model types: 
    % at every time point (1 sec), how many molecules have arrived AND not yet departed?
    current(i) = sum(arrives<i & departs>i);
end

end % end of function