CBE 504 HW2

birthdeath.m

function [nlist, tlist] = birthdeath(l,K,numsteps)
    nlist =[0]; % list of numbers of transcripts at corresponding place in t
    tlist =[0]; % list of times. drawn from bernoulli distribution.
    evlist=[0]; % list of events. degradation is -1, transcription is 1,neither is 0.
    t  = 0;
    n  = 0;
    ev = 1;
    for j=[1:numsteps]
        t = t - log(1-rand(1))/(K*n+l);
        tlist = [tlist t];
        p_tr = l/(l+K*n); % bernoulli probability of a transcription event.
        p_degr = 1-p_tr;  % bernoulli probability of a degradation event.
        transcription = binornd(1,p_tr);
        degradation = binornd(1,p_degr);
        ev = transcription - degradation;
        evlist = [evlist ev];
        n = n + ev;
        nlist = [nlist n];
    end

hw2.m

%%  CBE 504 - HW2 - Tom Bertalan

%%
clear;

%%  #1C: two number-of-mRNA-transcripts trajectories
% n molecules at time t
% 
% # Chose T according to:
% 
% |P(T<=t) = 1-exp( -(K*n + l)*tau ) -> t=t+T|
% 
% using a uniformly-distributed variable |u=rand(1)|:
%    
% # Choose event according to Bernoulli distribution:
% 
% |p(transcription) = 1 - p(degradation) = l / (l + K*n) -> n=n-1|
% 
% |p(degradation) = 1 - p(transcription) = K*n / (l + k*n) -> n=n+1|
% 

nlist =[0]; % list of numbers of transcripts at corresponding place in t
tlist =[0]; % list of times. drawn from bernoulli distribution.
evlist=[0]; % list of events. degradation is -1, transcription is 1,neither is 0.

K = 0.1; % [min^-1]
l = 2;   % [min^-1]
numtrials=2;
figure(1);
hold on;
trialcolors = ['r', 'b'];
for i=[1:numtrials]
    [nlist, tlist] = birthdeath(l,K,256);
    plot(tlist,nlist,trialcolors(i));
end
xlim([0 5/K]);
title('#1C');
xlabel('time [min]');
ylabel('number of transcripts');



%%  #1D: Stochastic empirical time-dependent distribution function for number-of-mRNA-transcripts

K = 0.1; % [min^-1]
l = 2;   % [min^-1]
numtrials = 200; % 200
nlists = []; % each row is a trial. nlists(i,j) is the number of
             % transcripts existing at time tlists(i,j) in trial i.
tlists = []; % each row is a trial.

events = 256; % 256
h = waitbar(0,sprintf('#1D: running %d birth/death simulations', numtrials));
for i=[1:numtrials]
    [nlist, tlist] = birthdeath(l,K,events);
    tlists = [tlists; tlist];
    nlists = [nlists; nlist];
    waitbar(i/numtrials,h);
end
close(h);
numsteps = 64; % 64 This shouldn't be too large relative to 'events', or the graph gets choppy.
nlistsize = size(nlists);
rows = nlistsize(1);
cols = nlistsize(2);
maxn = max(max(nlists));
pt = zeros([numsteps+1 maxn+1]);
ymax = 5/K;
dt = ymax/numsteps;
tvec = [0:dt:ymax];
nvec = [0:maxn];
for n=[0:maxn]
    nvec(:,n+1) = n;
end
row=1;
h = waitbar(0, sprintf('#1D: compiling %d time-dependent distributions', numsteps));
for t=tvec
    for trial=[1:numtrials]
        for event=[1:events]
            if tlists(trial,event) > t
                if tlists(trial,event) < t+dt
                    n = nlists(trial,event);
                    pt(row,n+1) = pt(row,n+1) + 1;
                end
            end
        end
    end
    waitbar(row/numsteps,h);
    row = row+1;
end
close(h);
figure(2);
% surf(nvec,tvec,pt,'LineStyle','none') % 'LineStyle' stuff is only relevant for surf()
image(nvec,tvec,pt)%,'LineStyle','none') % 'LineStyle' stuff is only relevant for surf()
colormap('Bone') % likewise, 'colormap()' is only relevant for image()
title('#1D (brightness corresponds to count of cases, i.e., frequency)');
ylim([0 ymax]);
xlabel('n');
ylabel('t [sec]');
zlabel('count');

% Take slices of these average emperical distribution functions at several
% times.
emperical_slices = zeros(3,maxn+1);
slice_times = [5 25 45];
slicesize = size(slice_times);
num_slices = slicesize(2);
for i=1:num_slices
    diffs = abs(tvec-slice_times(i));
    ind = find(diffs == min(diffs));
    emperical_slices(i,:) = pt(ind,:);
end
analytical_slices = zeros(3,maxn+1);
for i=1:num_slices
    t = slice_times(i);
    mu = l/K^(1-exp(-1*K*t))
    for n=0:maxn
        analytical_slices(i,n+1) = mu^n/factorial(n)*exp(-1*mu);
    end
    analytical_slices(i,:) = analytical_slices(i,:) * max(emperical_slices(i,:))/max(analytical_slices(i,:));
end
figure(4)
hold on;
legends=[];
scatter([0:maxn],emperical_slices(1,:),'dr');
l1 = sprintf('t=%d, emperical',slice_times(1));
scatter([0:maxn],emperical_slices(2,:),'+b');
l2 = sprintf('t=%d, emperical',slice_times(2));
scatter([0:maxn],emperical_slices(3,:),'og');
l3 = sprintf('t=%d, emperical',slice_times(3));
plot([0:maxn],analytical_slices(1,:),'r');
l4 = sprintf('t=%d, analytical',slice_times(1));
plot([0:maxn],analytical_slices(2,:),'b');
l5 = sprintf('t=%d, analytical',slice_times(2));
plot([0:maxn],analytical_slices(3,:),'g');
l6 = sprintf('t=%d, analytical',slice_times(3));
legend(l1,l2,l3,l4,l5,l6);


%%  #2: 1D random walk.
numtrials=1000;
% figure(4);
% hold on;
N = 100; % number of steps
finalpositions = [];
l = 1;
h = waitbar(0, '#2: 1D random walks');
modthis = numtrials/100;
for i=1:numtrials;
    if mod(i,modthis) == 0
        waitbar(i/numtrials,h);
    end
    positions = randwalk_1d(N);
%     plot(positions);
    finalpositions = [finalpositions; positions(end)];
end
close(h);

figure(5);
hold on;
[nout,xout] = hist(finalpositions);
bar(xout,nout,'histc');
minpos = min(finalpositions);
maxpos = max(finalpositions);
xlist = [minpos:maxpos];
p_true = [];
p_norm = [];
sigma = sqrt(N*l^2);
% sigma=10;
mu=0;
for x=xlist
    p_true = [p_true 1/sqrt(2*pi*N*l^2)*exp(-1*x^2/2/N/l^2)];
    p_norm = [p_norm 1/sigma/sqrt(2*pi) * exp( -1*(x-mu)^2/2/sigma^2 ) ];
end
p_true = p_true * max(nout) / max(p_true);
p_norm = p_norm * max(nout) / max(p_norm);
scatter(xlist,p_true,'b')
plot(xlist,p_norm,'g');
xlabel('final position');
ylabel(sprintf('count in %d trials', numtrials))
legend('simulated','exact','normal');


%%  #3: 3D random walk.
numtrials=256;
N = 512; % number of steps
finalpositions = [];
displacements = zeros(numtrials,N);
h = waitbar(0, '#3: 3D random walks');
modthis = numtrials/100;
for i=1:numtrials;
    if mod(i,modthis) == 0
        waitbar(i/numtrials,h);
    end
    positions = randwalk_3d(N);
    x = positions(1,:);
    y = positions(2,:);
    z = positions(3,:);
%     plot3(x,y,z);
    finalposition = [x(end) y(end) z(end)];
    finalpositions = [finalpositions; finalposition];
    displacement=zeros(1,N);
    for n=[1:N]
        displacement(n) = x(n)^2 + y(n)^2 + z(n)^2;
    end
    displacements(i,:) = displacement;
end
close(h);
ms_displacements = zeros(N,1);
for n=1:N;
    ms_displacements(n,1) = mean(displacements(:,n));
end
figure(6);
hold on;
plot(ms_displacements);
title('average displacement over time');
ylabel(sprintf('squared displacement over time, mean for %d trials', numtrials));
xlabel('time (number of steps)');

randwalk_1d.m

function positions = randwalk_1d(N)
    positions = [];
    position = 0;
    for i=0:N
        position = position + round(rand(1)*2)-1;
        positions = [positions position];
    end

randwalk_3d.m

function positions = randwalk_3d(N)
    positions = zeros(3,N+1);
    position = 0;
    for d=1:3
        currentresult = randwalk_1d(N);
        positions(d,:) = currentresult;
    end