ruebefrei

mann, hab ich nen ohrwurm! seit ich 'say hey' in dieses script eingebaut habe... gern mal ausfuehren (matlab + cftool required)...

function main
%% Uses ODE45 to solve
%% EXTENSIONS QSS NETWORK
%% and returns oscillations
clear all
%% INITIAL CONDITION VECTOR Zo and INITIAL CONDITION VALUES
Zo = [];
Zo(1, 1) = 0;
Zo(1, 2) = 0;
Zo(1, 3) = 0;
Zo(1, 4) = 1; %% SIC1 S
Zo(1, 5) = 0;
%% Constant Production rates v
v = [];
v(1) = 0.05;
v(2) = 0;
v(3) = 0;
v = v';
%% timespan
tspan = [0,15000];
%% REACTION RATE CONSTANTS
k_selfinh = 1;%.5;
k_act = .5;%.75;
k_inh = 1;%2;

v_scale = 1;
k_degs = 0;


part = [1 1 1 1 1 1 1 1 1 1 0.06 0.002];

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%_oscillations_finder_%%%%%%%%%%%%%%%%%%
par = part(1,:);
par1 = par(1);
par2 = par(2);
par3 = par(3);
par4 = par(4);
par5 = par(5);
par6 = par(6);
par7 = par(7);
par8 = par(8);
par9 = par(9);
par10 = par(10);
q = par(11);

k2counter = 0;
kk = [];
f = []; A = []; Asic1 = []; meanmax = []; meanmin = [];
meanmaxsic1 = []; meanminsic1 = [];
for k2index = .1:.1:10

fAcounter = 0;
k2 = k2index*par(12);
k2counter = k2counter + 1;
kk_indices = [];
%%%----------------parameters-------------------%
K_selfinh = diag([0,par9,par10]);

K_act = [0, 0, 0;...
par1, par4, 0;...
par2, par3, par5 ];


K_inh = [0, par6, par7;...
0, 0, par8;...
0, 0, 0 ];

for index1 = -4:0.1:-2.8
for index2 = -2:0.1:-1
disp([k2counter])
% for index1 = -3.5:0.02:-2.5;...-3;...
% for index2 = -2:0.02:-1;...-1.14;...

k_prod = 10^index1;
k_degr = 10^index2;


try
%%%----------------solver--------------------------%
options = odeset('RelTol',1e-3); %%%, 'OutputFcn',@odephas3);
[t,Z] = ode113( @input_ODEs,tspan,Zo,options,...
v_scale*v,q,k2,k_degs,k_selfinh*K_selfinh,k_act*K_act,k_inh*K_inh,k_prod,k_degr);

for indexZ = 1:size(Z,1)
Z(indexZ,6) = sum(Z(indexZ,2:3))/3;
if round(Z(indexZ,6)) == 0, Z(indexZ,6) = 0; end
end

peaks = diff(Z(:,6));
start = 0;
for indexpeaks = 1:size(peaks,1)-7
if sign(peaks(indexpeaks)) == 1 ...
&& sign(peaks(indexpeaks + 1)) == -1 ...
&& sign(peaks(indexpeaks + 2)) == -1 ...
&& sign(peaks(indexpeaks + 3)) == -1 ...
&& sign(peaks(indexpeaks + 4)) == -1 ...
&& sign(peaks(indexpeaks + 5)) == -1 ...
&& sign(peaks(indexpeaks + 6)) == -1 ...
&& sign(peaks(indexpeaks + 7)) == -1
start = start +1;
end
end

figure(1); subplot(2,1,1); plot(t,Z(:,1:3),t,Z(:,4),'c'); title([num2str([k_prod, k_degr])]); legend('clb5,6tot','clb3,4tot','clb1,2tot',3); set(gca,'XTick', 0:150:tspan(2));
subplot(2,1,2); plot(t,Z(:,4),'c'); legend('s'); set(gca,'XTick', 0:150:tspan(2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if start > 4 ...&& sum(minima(:,2))/size(minima,1) < 3


disp('say hey!')


for index3 = index1-0.1:0.06:index1+0.1;...
for index4 = index2-0.1:0.06:index2+0.1;...

k_prod = 10^index3;
k_degr = 10^index4;


options = odeset('RelTol',1e-3); %%%, 'OutputFcn',@odephas3);
[t,Z] = ode113( @input_ODEs,tspan,Zo,options,...
v_scale*v,q,k2,k_degs,k_selfinh*K_selfinh,k_act*K_act,k_inh*K_inh,k_prod,k_degr);


for indexZ = 1:size(Z,1)
Z(indexZ,6) = sum(Z(indexZ,2:3))/3;
if round(Z(indexZ,6)) == 0, Z(indexZ,6) = 0; end
end

peaks = diff(Z(:,6));
maxcounter = 0;
maxima = [];
for indexpeaks = 1:size(peaks,1)-7
if sign(peaks(indexpeaks)) == 1 ...
&& sign(peaks(indexpeaks + 1)) == -1 ...
&& sign(peaks(indexpeaks + 2)) == -1 ...
&& sign(peaks(indexpeaks + 3)) == -1 ...
&& sign(peaks(indexpeaks + 4)) == -1 ...
&& sign(peaks(indexpeaks + 5)) == -1 ...
&& sign(peaks(indexpeaks + 6)) == -1 ...
&& sign(peaks(indexpeaks + 7)) == -1
maxcounter = maxcounter +1;
maxima(maxcounter,1:2) = [maxcounter, Z(indexpeaks+1,6)];
end
end

if ~isempty(maxima)
exclude = 3;
fit = expfit(maxima(exclude:size(maxima,1),1),maxima(exclude:size(maxima,1),2));
[stretch] = coeffvalues(fit);

if abs(stretch(2)) < 0.001 && maxcounter > 5
kk_indices(size(kk_indices,1)+1,1:2) = [index3, index4];
fAcounter = fAcounter + 1;
[f A Asic1 meanmax meanmin meanmaxsic1 meanminsic1] = fA(t, Z, peaks, k2counter, fAcounter, ...
f, A, Asic1, meanmax, meanmin, meanmaxsic1, meanminsic1);
end
end
clear t Z peaks maxima minima
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

catch
end

end
end

kk(1, size(kk,2)+1:size(kk,2)+2) = [q, k2];
kk(2:size(kk_indices,1)+1, size(kk,2)-1:size(kk,2)) = kk_indices;
clear kk_indices
save fAaddkirr kk f A Asic1 meanmax meanmin meanmaxsic1 meanminsic1
end


function [f A Asic1 meanmax meanmin meanmaxsic1 meanminsic1] = fA(t, Z, peaks, k2counter, fAcounter, ...
f, A, Asic1, meanmax, meanmin, meanmaxsic1, meanminsic1)

peaksic1 = diff(Z(:,4));

tmax_max = [];
tmin_min = [];
tmax_maxsic1 = [];
tmin_minsic1 = [];
for indexpeaks = 1:size(peaks,1)-7
if sign(peaks(indexpeaks)) == 1 ...
&& sign(peaks(indexpeaks + 1)) <= 0 ...
&& sign(peaks(indexpeaks + 2)) <= 0 ...
&& sign(peaks(indexpeaks + 3)) <= 0 ...
tmax_max(size(tmax_max,1)+1,1:2) = [t(indexpeaks+1,1),Z(indexpeaks+1,6)];
end
if sign(peaksic1(indexpeaks)) == 1 ...
&& sign(peaksic1(indexpeaks + 1)) <= 0 ...
&& sign(peaksic1(indexpeaks + 2)) <= 0 ...
&& sign(peaksic1(indexpeaks + 3)) <= 0 ...
tmax_maxsic1(size(tmax_maxsic1,1)+1,1:2) = [t(indexpeaks+1,1),Z(indexpeaks+1,4)];
end
if sign(peaks(indexpeaks)) == -1 ...
&& sign(peaks(indexpeaks + 1)) >= 0 ...
&& sign(peaks(indexpeaks + 2)) >= 0 ...
&& sign(peaks(indexpeaks + 3)) >= 0 ...
tmin_min(size(tmin_min,1)+1,1:2) = [t(indexpeaks+1,1),Z(indexpeaks+1,6)];
end
if sign(peaksic1(indexpeaks)) == -1 ...
&& sign(peaksic1(indexpeaks + 1)) >= 0 ...
&& sign(peaksic1(indexpeaks + 2)) >= 0 ...
&& sign(peaksic1(indexpeaks + 3)) >= 0 ...
tmin_minsic1(size(tmin_minsic1,1)+1,1:2) = [t(indexpeaks+1,1),Z(indexpeaks+1,4)];
end
end


f(fAcounter, k2counter) = (sum(diff(tmax_max(3:size(tmax_max,1),1)))/(size(tmax_max,1)-3) ...
+ sum(diff(tmin_min(3:size(tmin_min,1),1)))/(size(tmin_min,1)-3))/2;

meanmax(fAcounter, k2counter) = sum(tmax_max(3:size(tmax_max,1),2))/(size(tmax_max,1)-2);
meanmin(fAcounter, k2counter) = sum(tmin_min(3:size(tmin_min,1),2))/(size(tmin_min,1)-2);
meanmaxsic1(fAcounter, k2counter) = sum(tmax_maxsic1(3:size(tmax_maxsic1,1),2))/(size(tmax_maxsic1,1)-2);
meanminsic1(fAcounter, k2counter) = sum(tmin_minsic1(3:size(tmin_minsic1,1),2))/(size(tmin_minsic1,1)-2);

A(fAcounter, k2counter) = meanmax(fAcounter, k2counter) - meanmin(fAcounter, k2counter);
Asic1(fAcounter, k2counter) = meanmaxsic1(fAcounter, k2counter) - meanminsic1(fAcounter, k2counter);
return


function output = expfit(t,data)
ok_ = isfinite(t) & isfinite(data);
if ~all( ok_ )
warning( 'GenerateMFile:IgnoringNansAndInfs', ...
'Ignoring NaNs and Infs in data' );
end
st_ = [0 0];
ft_ = fittype('exp1');
cf_ = fit(t(ok_),data(ok_),ft_,'Startpoint',st_);
% ft_ = fittype('poly1');
% cf_ = fit(t(ok_),data(ok_),ft_);
output = cf_;
return



%% a function which returns a rate change vector
function [dz_dt]= input_ODEs(t,z,v,q,k2,k_degs,K_selfinh,K_act,K_inh,k_prod,k_degr)


if any(z(1:5) > 200), return, end

z = z(1:5); %% variable vector z
eta = 1/(1+z(4)/q); %% eta calculation

for m = 1:size(v,1) %% QSS ODE equations loop
dz_dt(m) = v(m) ...
+ K_act(m, 1:3)*eta*z(1:3) ...
- K_selfinh(m, 1:3)*eta*z(1:3) ...
- eta*z(m)*K_inh(m, 1:3)*eta*z(1:3);
end

dz_dt(4) = - k_degs*z(4) ... %% Sic1 degradation equation
- k2*(1-eta)*sum(z(1:3)) ...
+ z(5);

dz_dt(5) = k_prod - k_degr*z(5)*eta*sum(z(1:3));

dz_dt = dz_dt'; %% transpose dz_dt so it is a column vector
return