function AK_SensitivityAnalysis_20151120
%% data for the model
Filename = '20151120_20151119_RedParam_TCell_1500_1';
load(Filename);
p = pmin;
pinit = pmin;

x0(1) = 0.54; 
x0(2) = 0.07;x0(3) = 10;
x0(4)=0.05; x0(5)=0.1; x0(6)=p(32); x0(7)=0.1;
CaEx=1500; %uM
tdata=1:1200;



ratio=[0.8 0.9 1/1.05 1/1.04 1/1.03 1/1.02 1/1.01 1 1.01 1.02 1.03 1.04 1.05 1.1 1.25] ;

tspan = [0:0.1:1*max(tdata)];
q=zeros(length(ratio),length(p));
peaktime=zeros(length(ratio),length(p));
tau1=zeros(length(ratio),length(p));

ff = fittype('exp2');

%General model Exp2:
%f2(x) = a*exp(b*x) + c*exp(d*x), Tau1 is b, which will be used below

for k=1:length(p)
    k
    for kk=1:length(ratio)
        kk
        p(k)=ratio(kk)*p(k);
        [tsol1,xsol1] = ode23s(@(t,x)AK_f1_20150402_TCell_pub_tofit(t,x,p),tspan,x0);
        qq=max(xsol1(:,4))
        q(kk,k)=qq;
        peakT=(find(xsol1(:,4)==qq))/10;
        peaktime(kk,k)=peakT;
%         hold on, plot(tsol1,xsol1(:,4));
        decay = xsol1(10*peakT:end,4);
        size(decay)
        size(tsol1(10*peakT:end))
        [c2,gof2] = fit(tsol1(10*peakT:end),decay,ff,'Robust', 'LAR',...
            'Lower',[0,-1,0,-1],...
            'Upper',[10,0,10,0]); %,...
        tau1(kk,k)=c2.b;
        p(k)=pinit(k);
    end
end      

figure, imagesc(peaktime)
figure, imagesc(q)

save([Filename,'_Peak_Young'],'peakT')
save([Filename,'_Max_Young'],'qq')
save([Filename,'_Tau1_Young'],'tau1')
save(['PeakTimeTau_' Filename])