function [dF0,dalpha,dA0,dA1,dA2,pF0,palpha,pA0,pA1,pA2]=shuffGPP(x1,x2,y1,y2,b01,b02,N)
inoF0=0;
inoalpha=0;
inoA0=0;
inoA1=0;
inoA2=0;

modelfun='y~(b1*b2*x1/(b1+b2*x1))*(b3+b4*x2+b5*x3)';

f1 = fitnlm(x1,y1,modelfun,b01);%nlinfit
fb1=f1.Coefficients.Estimate;
fb1=fb1';
f2 = fitnlm(x2,y2,modelfun,b02);
fb2=f2.Coefficients.Estimate;
fb2=fb2';
F01=fb1(1,1);
F02=fb2(1,1);
alpha1=fb1(1,2);
alpha2=fb2(1,2);
A01=fb1(1,3);
A02=fb2(1,3);
A11=fb1(1,4);
A12=fb2(1,4);
A21=fb1(1,5);
A22=fb2(1,5);

dF0=F01-F02;
dalpha=alpha1-alpha2;
dA0=A01-A02;
dA1=A11-A12;
dA2=A21-A22;

count=1;
Nd1=size(x1,1);
Nd2=size(x2,1);

 while count<=N
  
  arn1=[];arn2=[];y1r=[];y2r=[];
  ip=randperm(Nd1+Nd2);
  xx=[x1;x2];
  yy=[y1;y2];
  arn1=xx(ip(1:Nd1),:);
  arn2=xx(ip(Nd1+1:Nd1+Nd2),:);
  y1r=yy(ip(1:Nd1));
  y2r=yy(ip(Nd1+1:Nd1+Nd2));

    f1rn = fitnlm(arn1,y1r,modelfun,b01);
    fb1rn=f1rn.Coefficients.Estimate;
    fb1rn=fb1rn';
    f2rn = fitnlm(arn2,y2r,modelfun,b02);
    fb2rn=f2rn.Coefficients.Estimate;
    fb2rn=fb2rn';

    F01rn=fb1rn(1,1);
    F02rn=fb2rn(1,1);
    alpha1rn=fb1rn(1,2);
    alpha2rn=fb2rn(1,2);
    A01rn=fb1rn(1,3);
    A02rn=fb2rn(1,3);
    A11rn=fb1rn(1,4);
    A12rn=fb2rn(1,4);
    A21rn=fb1rn(1,5);
    A22rn=fb2rn(1,5);

  if(max(abs(fb1rn))>100 || max(abs(fb2rn))>100 || min(abs(fb1rn))<1.e-5 || min(abs(fb2rn))<1.e-5)
      continue
  end


  dF0rn=F01rn-F02rn;
  dalpharn=alpha1rn-alpha2rn;
  dA0rn=A01rn-A02rn;
  dA1rn=A11rn-A12rn;
  dA2rn=A21rn-A22rn;
   
  if(abs(dF0rn)>=abs(dF0))
      inoF0=inoF0+1;
  end 
  if(abs(dalpharn)>=abs(dalpha))
      inoalpha=inoalpha+1;
  end 
  if(abs(dA0rn)>=abs(dA0))
      inoA0=inoA0+1;
  end
  if(abs(dA1rn)>=abs(dA1))
      inoA1=inoA1+1;
  end
  if(abs(dA2rn)>=abs(dA2))
      inoA2=inoA2+1;
  end
  
  count=count+1;
 
 end

pF0=inoF0/N;
palpha=inoalpha/N;
pA0=inoA0/N;
pA1=inoA1/N;
pA2=inoA2/N;