Symsen montecarlo.m

From 2006.igem.org

clear all;
close all;

anz = 10;

%Parameters to be tested: KSa, KSr, nSa, nSr, SA
%Ranges:
%KSa: 10-200 nM, da KSr 5 mal stärker
%KSr: 10-1000 nM
%nSa: 1-5
%nSr: 1-5
%SA: 10-1000 nM

%For the zincfingers:
%K: 10-1000 nM
%kd: 0.1-10 nM min^-1
%k: 100-1000 min^-1

%Set the random numbers:
K_r = unifrnd(200,1000,1,anz);
kd_r = unifrnd(0.1,10,1,anz);
k_r = unifrnd(100,1000,1,anz);
KS_r = unifrnd(100,1000,1,anz);
n_r = unifrnd(1,3,1,anz);
nS_r = unifrnd(1,3,1,anz);
SA_r = unifrnd(100,1000,1,anz);

%Order them:
n_r = unique(n_r);
nS_r = unique(nS_r);
K_r = unique(K_r);
KS_r = unique(KS_r);
kd_r = unique(kd_r);
k_r = unique(k_r);
SA_r = unique(SA_r);

   
%nominal values:

n  = 2;  %1 keine Kooperativität der Zinkfinger -> alle n immer = 1
nS  = 2; % Kooperativität von S (als Aktivator)
K  = 100; %500 Affinität (höheres K = höhere Affinität) 0.001-0.1 nM (für n=1)
KS  = 100; %100 Prm
k  = 700; %1000 100-1000 min^-1
kd  = 1.5; %1.5 Degradationsrate der Zinkfinger 0.1-10 nM min^-1
period = 30;
SA = 300; %Amplitude von S

% stepsize
tstep = 5;
% duration
tfinal = 120;

y0 = [2,0,0,0];
opt = odeset('RelTol',1e-3,'AbsTol',1e-6);


%Calculate the objective function for every set of parameters:
%The output function is the concentration of Zincfinger 4 (ZF4 -> y(4))

para1 = K_r
para2 = KS_r


para = [n,nS,nS,K,KS,KS,k,kd,period,SA];
[t,y] = ode15s('ode_counter',[0:tstep:tfinal],y0,opt,para);
f_obs = y(:,4);
f_obj = ones(length(para1),length(para2));
for i = 1:length(para1),
    for j = 1:length(para2),
        % Set here the right parameters!!
        para = [n,nS,nS,para1(i),para2(j),para2(j),k,kd,period,SA];
        [t,y] = ode15s('ode_counter',[0:tstep:tfinal],y0,opt,para);
        f_per = y(:,4);
        f_obj(i,j) = sum((f_obs - f_per).^2);
    end
end

figure(1);
surf(para2,para1,f_obj);
xlabel('para2');
ylabel('para1');

criterion = sum(sum(f_obj))/anz^2;

figure(7);
%hold on;
surf(para2,para1,f_obj);
xlabel('para2');
ylabel('para1');
hold on;
crit = ones(length(para1),length(para2))*criterion;
surf(para2,para1,crit);


f_para2 = 0.5*sign(f_obj - criterion)+0.5 %1, falls über criterion, sonst 0

para2_ueber = sum(f_para2);
para2_unter = anz-para2_ueber;

f_para1 = 0.5*sign(f_obj' - criterion)+0.5;

para1_ueber = sum(f_para1);
para1_unter = anz-para1_ueber;

figure(2);
hold on;

title('para2 Anzahl');
plot(para2,para2_ueber/anz);
plot(para2,para2_unter/anz,'--');
legend('above criterion','below criterion');

figure(3);
hold on;

title('para1 Anzahl');
plot(para1, para1_ueber/anz);
plot(para1, para1_unter/anz,'--');
legend('above criterion','below criterion');

para2_ueber_cum = zeros(1,anz);
para2_ueber_cum(1) = para2_ueber(1);
for i = 2:anz,
    para2_ueber_cum(i) = para2_ueber_cum(i-1) + para2_ueber(i);
end

para2_unter_cum = zeros(1,anz);
para2_unter_cum(1) = para2_unter(1);
for i = 2:anz,
    para2_unter_cum(i) = para2_unter_cum(i-1) + para2_unter(i);
end

para1_ueber_cum = zeros(1,anz);
para1_ueber_cum(1) = para1_ueber(1);
for i = 2:anz,
    para1_ueber_cum(i) = para1_ueber_cum(i-1) + para1_ueber(i);
end

para1_unter_cum = zeros(1,anz);
para1_unter_cum(1) = para1_unter(1);
for i = 2:anz,
    para1_unter_cum(i) = para1_unter_cum(i-1) + para1_unter(i);
end

figure(4);
hold on;

title('para2 kumuliert');
plot(para2, para2_ueber_cum);
plot(para2, para2_unter_cum,'--');
legend('above criterion','below criterion');

figure(5);
hold on;

title('para1 kumuliert');
plot(para1, para1_ueber_cum);
plot(para1, para1_unter_cum,'--');
legend('above criterion','below criterion');

para1_corr = crosscorr(para1_ueber_cum, para1_unter_cum,1)
para2_corr = crosscorr(para2_ueber_cum, para2_unter_cum,1)