function dose_response_normal_dist
% This program estimate the probability of infection at certain doses:
% dose-response
tic
clc; clear all;
global r a  beta K delta sigma

delta=20;


u=10.^linspace(1,12,12);

sn=4;% number of repeatation
yy=nan(length(u),sn)
for s1=1:sn
count=zeros(length(u),1);
for j=1:1:length(u)
    IP=u(j);%10^(j);
      if(j<=6)k1=10^4;end
      if(j>6)k1=10^3;end
      if(j<=4)k1=10^5;end
      if(j<=3)k1=10^6;end
      if(j<=2)k1=10^8;end
      n0=0;n1=0;n2=0;n3=0;n4=0;n5=0;n6=0;
      mLstr=zeros(k1,1);
  %    k1=1000;
 for i=1:k1
        
% n0 is unrealistice; n2 to n5 are positive responses; n1 and n6 are negative
% response. sum of these should equal to the total number of samples k1
% n0: unrealistic parameters; n4: threshold in negative quadrant!;n1: L_str does not exist

r=random('norm',0.863,.43);% uniform distribution
ab=[ -3.351849051	4.904051711	0.529674595	0.750450028];%% estimated data collecting from the excel sheet

meana0=ab(1);meanb=ab(2);SDa0=ab(3);SDb=ab(4);
a0=random('norm',meana0,SDa0);%-3.69 1
beta=random('norm',meanb,SDb); % uniform distribution 6.74,.75

sigma=random('norm',0.4832,0.1405);% Hellmig2006
K0=random('norm',8,1);% uniform distribution
if(r<0||0<a0||sigma<=0||K0<=0||beta<=0)
   %    [r,a0,sigma,K0,beta,n0];
    n0=n0+1;% testing if all cases sum to k1
    continue;
end
a=10^a0;K=10^K0;
D=(r*K*a)^2+2*r^2*K*a+r*r-4*r*a*beta*K; % L_str exists if D>0. 
%%-----folloing is based on vector field analysis- Figure 2 of manuscript
%% (1) case D<0.  If D<0, only Lo exists and stable
if(r<(4*a*K*beta)/(1+a*K)^2)% negative response, Case I of Section 3.1; 
    % this is equivalent to D<0
    n1=n1+1;% testing if all cases sum to k1
    continue;
end
%%% (2) next case D>0
if(r>beta) % positive response, Case II of Section 3.1
    count(j)=count(j)+1; 
     n2=n2+1;% testing if all cases sum to k1; also postive response
    continue;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



L_str=(r*K*a-r-sqrt(D))/(2*r*a);
if(L_str<0)% only two equilibria (L0 and L+) exist; L+ is table
      count(j)=count(j)+1;
     n3=n3+1;% testing if all cases sum to k1
     continue;
end
m=(-(delta+sigma)+r*L_str/K-a*beta*L_str/(1+a*L_str)^2)/sigma;% slope of the eigenvector threshold
     if(m>0) % not sure? only small cases, should not affect the result
          count(j)=count(j)+1;
         n4=n4+1;% testing if all cases sum to k1
         continue;
     end
     mLstr(i)=L_str;
 if (IP>-m*L_str) % if L_str is -ve L+ is always stable=> positive count
     count(j)=count(j)+1;
     n5=n5+1;% positive response out of bistable scenerio
     continue;
 end % if 
 n6=n6+1;% negative response out of bistable scenerio; testing if all cases sum to k1
%----
%-------
end %%i loop
%[n0 n1 n2 n3 n4 n5 n6, count(j)]
count(j)=count(j)/(k1-n0);% k1-p

end % j loop
%hold on;
yy(:,s1)=count;
end % s1 loop

format long
y=mean(yy,2)% mean of columns
format short
figure(1)
hold on;

    %%-- data writing in a file
    filename = 'data1.xlsx';
    sheet = 2;
    xlswrite(filename,[yy'],sheet,'A1')
    %------end data write ----------
 plot(log10(u), y,'Bo--','LineWidth',2, 'MarkerSize',6)% response data 
  xlabel('Inoculum (Log dose)','fontsize',16);ylabel('Probability of infection','fontsize',16)
toc

end

