function organoid_growth(ageDivH,Nneigh,maxRad)

% -- to remember random number  
  seed=randi(5000);
  rng(seed);

  
% -- showing images and saving the data:
  ShowImages   =    1;   % 0=dont show images; 1=show images
  HowOftenShow =  672;   % choose an iteration step for drawing
  SaveData     =    1;   % 0=dont show images; 1=show images

  
% -- cell & microenvironment characteristics 
  cellLine = 3; StiffRep = 0.45; StiffAdh = 0.03;
  microEnv = 1;  
  drug     = 1;   
  concDrug = 0;   
  nu=10;          % viscosity of the media

  pathsave=['cell',num2str(cellLine),'-mE',num2str(microEnv)];



  
% -- domain boundary:
  xmin=-200; xmax=-xmin; ymin=xmin; ymax= xmax; zmin=xmin; zmax=xmax;  

% -- time progression:
  Ndays     = 14;          % simulation time [days]
  dt        = 30;          % interval for diffusion in minutes
  Niter     = floor(Ndays*24*60/dt); % number of iterations corresponing  
                                     % to 14 days
                                  
% -- initial cell and its coordinates:
  Ncell    = 1;           % 1 initial cell
  cell     = [0,0,0];     % cell coordinates [x,y,z];
  neigh    = 0;           % number of neighbors

% -- define cell size (originally 8 for #4 cell line and 9.5 for others)
  rad(Ncell,1)=maxRad;

 
% -- other parameters associated with cell lifespan 
  ageDiv=ageDivH*60;
  age            = [ageDiv/2,ageDiv];  % age [current,mature]  
  cell_phase     = [0,0];              % phase of cell cycle ([age-based,numeric-value]) 
  cellhist       = 0;                  %

% -- variables needed to calculate area of the projection of the spheroid
  maxd    =   2*maxRad;  
  maxdiam = 1.5*maxRad;

  
  
%% ------ drawing figure
  if (ShowImages==1)
    ff=figure('Position',[150,150,1000,550]);
    set(ff,'PaperPositionMode','auto')
    set(ff,'InvertHardcopy','off');
    set(gca,'color','white')
	set(gcf,'color','white')
    clf
    colordef white
  end

  
 %% ------ experimetnal data & fitting
    time=[0,144,336,480,672];
    aver=[2*maxRad,53.5104430000000,149.4803700000000,...
          219.0245300000000,247.1622400000000];
    std =[0,17.824574999999999,37.986648000000002,...
          55.869179000000003,60.979132000000000];
    datafit=polyfit(time,aver,3);
    timefit=0:671;
    valfit =polyval(datafit,timefit);
    
    
%% ======= MAIN LOOP ======

  iter=0; stopcond=0;
  while ((iter<Niter)&&(stopcond==0))
    iter=iter+1;
       
    % --- CELL DIVISION ---
    age(1:Ncell,1) = age(1:Ncell,1)+dt;         %% -------- CELL AGE
    for icell=1:Ncell                           %% -------- DIVISION
	  if (cell_phase(icell,2)==4) && (age(icell,1)>=age(icell,2)) &&...
         (neigh(icell,1)<= Nneigh) 
        Ncell  = Ncell+1;      % new cell index
        thet   = pi*rand;      % random angle for location of daughter cells 
        phi    = 2*pi*rand;    % random angle for location of daughter cells
 
        newrad = rad(icell,1); % radius of a mother cell 
        cell(Ncell,1:3)=cell(icell,1:3)+0.5*newrad*...    % cell location
             [sin(thet)*cos(phi),sin(thet)*sin(phi),cos(thet)];
        cell(icell,1:3)=cell(icell,1:3)-0.5*newrad*...
             [sin(thet)*cos(phi),sin(thet)*sin(phi),cos(thet)];
        rad(icell,1) = 0.65*newrad;  % radius of two newborn cells 
        rad(Ncell,1) = 0.65*newrad;  
        age(icell,1) = 0;            % age of new daughter cells
        age(Ncell,1) = 0;
        agemem       = age(icell,2);                 % mother maturation age
        age(icell,2) = ageDiv+0.1*(0.5-rand)*ageDiv; % maturation age of new  
        age(Ncell,2) = ageDiv+0.1*(0.5-rand)*ageDiv; % daughter cells  +/- 5%
        neigh(icell,1) = 0;          % number of neighbors for daughters
        neigh(Ncell,1) = 0;                               
      end  % if
    end  % for

    
    % --- phases of cell cycle: G0, G1, S, G2, M; % of time ---
    for icell=1:Ncell
      cell_phase(icell,1)=100*age(icell,1)/age(icell,2);
      if (cell_phase(icell,1) >=0) && (cell_phase(icell,1) < 5)
        cell_phase(icell,2) = 0;
      elseif (cell_phase(icell,1) >=5) && (cell_phase(icell,1) < 45)
        cell_phase(icell,2) = 1;
      elseif (cell_phase(icell,1) >=45) && (cell_phase(icell,1) < 80)
        cell_phase(icell,2) = 2;
      elseif (cell_phase(icell,1) >=80) && (cell_phase(icell,1) < 95)
		cell_phase(icell,2) = 3;
      else
        cell_phase(icell,2) = 4;
      end
    end


    % --- cell-cell physical interactions ---
    [forceR,forceA,neigh] = DefineForce(rad,cell,StiffRep,StiffAdh,Ncell,maxRad);

    % --- cell relocation due to forces ---
    cell(1:Ncell,1:3)=cell(1:Ncell,1:3)+dt*(forceR(1:Ncell,1:3)+forceA(1:Ncell,1:3))/nu;    

    % --- cell growth and size increase ---
    for icell=1:Ncell
      if (cell_phase(icell,2)<4) && (rad(icell,1)<maxRad)
        alph=0.35*dt/(0.75*age(icell,2));  % cells grows in 75% of its maturation time
        rad(icell,1)=rad(icell,1)+alph*maxRad;    
      end 
    end
    
    
    % --- spheroid diameter calculation ---
    maxd_old=maxd;
    [maxar,maxd]=maxDiameter(cell,maxRad,maxd_old);
 
    % --- remember diameter, area and cell count for all iterations
    maxdiam(iter)  = maxd; 
    cellhist(iter) = Ncell;  


  
    %%% conditions for stoping the code 
    if (iter==144)         % day 3, data 20% outside std
      if (abs(maxd-aver(2))<1.2*std(2))
      else  stopcond=1; end
    end
    if (iter==336)         % day 7, data 20% outside std
      if (abs(maxd-aver(3))<1.2*std(3))
      else  stopcond=1;  end
    end
    if (iter==480)         % day 10, data 20% outside std
      if (abs(maxd-aver(4))<1.2*std(4))
      else  stopcond=1; end
    end
    if (iter==624)         % day 14, data 10% outside std
      if (abs(maxd-aver(5))<1.1*std(5))
      else  stopcond=1;  end
    end        
    %%% end of conditions for stopping the code

    

       
    % ===============
    % --- DRAWING ---
    if (ShowImages==1)&&((mod(iter,HowOftenShow)==0)||(iter==Niter))
      clf
      col=['g','m','r','b'];
      
      % --- subplot: SPHEROID
      subplot(2,2,[1,3])  
        set(gca,'color','white')
        DrawSpheroid(cell,Ncell,rad,xmin,xmax,ymin,ymax,zmin,zmax,iter,...
        Niter,cellLine,microEnv,drug,concDrug,ageDivH,Nneigh,cell_phase)

      % --- subplot: CELL COUNTER    
      subplot(2,2,2)
        set(gca,'color','white'); set(gcf,'color','white')
        plot(1:iter,cellhist(1:iter),'d','Color',col(cellLine))
        title(['Number of cells: ',num2str(Ncell)],'FontSize',14,'Color',col(cellLine))
        ax = gca;
        xtop=Niter+20;
        set(ax,'XTickLabel',{'0','2','4','6','8','10','12','14'});
        if (concDrug==0) axis([0,xtop,0,1.5*Ncell])  
        else axis([0,xtop,0,1.25*Ncell]) 
        end
        xlabel('Days','FontSize',14)
        ylabel('# Cells','FontSize',14)
        set(gca,'fontsize',14)
        grid
        hold on
        
      % --- subplot: DIAMETER 
      subplot(2,2,4)
        set(gca,'color','white'); set(gcf,'color','white')
     
        % --- experimental data ---     
        axis([0,xtop,0,aver(4)+2.5*std(4)])
        hold on
        plot(time,aver,'ko','Linewidth',4,'Markersize',10)
        plot(time,aver+std,'k-.')
        plot(time,aver-std,'k-.')
        errorbar(time,aver,std,'ko')
        plot(timefit,valfit,'k.','Linewidth',1)
        

        % --- computational data ---      
        plot(1:iter, maxdiam(1:iter),'*','color',col(cellLine))
        title(['Diameter = ',num2str(maxd,'%5.3f'),' \mum'],'FontSize',14,'color',col(cellLine))
        ax = gca;
        xtop=Niter+20;
        set(ax,'XTickLabel',{'0','2','4','6','8','10','12','14'});         
        set(gca,'fontsize',14)
        grid
        xlabel('Days','FontSize',14) 
        ylabel('Diameter   (\mum)','FontSize',14) 
        hold on

      % compute r2 & norm for the last oteration        
      if (iter==Niter)       
        r2_s2ef=rsquare(maxdiam,valfit);       % r2 simulated data to exp fit      
        [RgyrSph,RgyrNorm] = rgyration(cell,Ncell,maxdiam,Niter);
     
     
        % add title to the spheroid image
        subplot(2,2,[1,3])        
          if (cellLine==1) text='-AT1'; elseif (cellLine==2) text='-CA1a';
          elseif (cellLine==3) text='-CA1d'; else text='A'; end
  
          title({['CELL LINE: MCF10',text,'     iter=',num2str(iter),...
            ' [',num2str(Niter),']'];['cell: #',num2str(cellLine),...
            ', mE: #',num2str(microEnv),', drug: #',num2str(drug),...
            ', conc: ',num2str(concDrug),'ng/mL, divAge: ',...
            num2str(ageDivH),', neigh #',num2str(Nneigh),', rad: ',...
            num2str(maxRad)];['r^2-s2ef=',num2str(r2_s2ef),...
            ';   Rgyr=',num2str(RgyrSph),';   RgyrNorm=',num2str(RgyrNorm)] ;...
            ['f-data=',num2str(datafit(1)),'x^3+',num2str(datafit(2)),...
            'x^2+',num2str(datafit(3)),'x+',num2str(datafit(4))];...
            },'FontSize',10)  
      end    
      
      pause(0.1)
    end  % end of DRAWING
    
  end % --- end of "for iter loop" or "while loop"

  
  
 
  % =================================================
  % --- SAVING DATA ONLY IF SIMULATION SUCCESSFUL ---
      
  if (stopcond==0)

    if (SaveData==1)
      % --- saving parameters ---  
      pathdir=MakeDirectory(pathsave,Nneigh,ageDivH,maxRad);     
      param=[xmin,xmax,ymin,ymax,zmin,zmax,dt,StiffRep,StiffAdh,nu,...
             Niter,seed,ageDivH,Nneigh,maxRad,cellLine,microEnv,drug,...
             concDrug,r2_s2ef,RgyrSph,RgyrNorm]';
      save([pathsave,'/',pathdir,'/parameters.txt'],'param','-ascii');
 
      % --- saving figure ---
      if (ShowImages==1) 
        mkdir([pathsave,'/',pathdir,'/Figs']);
        filename=[pathsave,'/',pathdir,'/Figs/fig',num2str(iter)]; 
        saveas(ff,filename,'jpg');
        savefig(ff,filename);
      end  

      % --- saving last data set ---   
      filename=[pathsave,'/',pathdir,'/tumor_',num2str(iter),'.txt']; 
      save(filename,'cell','-ascii');           % tumor data
         filename=[pathsave,'/',pathdir,'/cell_phase_',num2str(iter),'.txt']; 
         save(filename,'cell_phase','-ascii');  % cell phases
      filename=[pathsave,'/',pathdir,'/cell_radius_',num2str(iter),'.txt']; 
      save(filename,'rad','-ascii');            % cell radii
        filename=[pathsave,'/',pathdir,'/diameter_',num2str(iter),'.txt']; 
        save(filename,'maxd','-ascii');         % diameter data

      % --- saving data history ---   
      filename=[pathsave,'/',pathdir,'/CellHistory.txt']; 
      save(filename,'cellhist','-ascii');      % cell history data
      filename=[pathsave,'/',pathdir,'/DiamHistory.txt']; 
      save(filename,'maxdiam','-ascii');       % diameter history data
      filename=[pathsave,'/',pathdir,'/ExpFitting.txt']; 
      save(filename,'valfit','-ascii');        % experimental data fitting

      % --- saving correlations
       maxval=max(max(maxdiam),max(valfit));

    fg2=figure;
      corrfit2=fitlm(valfit,maxdiam);
      plot(corrfit2)
      grid on
      axis([0,maxval+5,0,maxval+5])
      axis equal
      axis([0,maxval+5,0,maxval+5])
      ylabel('simulated data','fontsize',15)
      xlabel('fitting to experimental data','fontsize',15)
      title(['correlation between fitting to experimental data and ',...
            ' simulated data'],'fontsize',15)
      filename=[pathsave,'/',pathdir,'/Figs/fig_correlation_f2d']; 
      saveas(fg2,filename,'jpg');

    end  % if SaveData
  end  % if stopcond  
  
end   % --- end of the main program



% ====================================================================== 
% --- FUNCTIONS

function path=MakeDirectory(pathsave,Nneigh,ageDivH,maxRad)    
  path=['ageDiv',num2str(ageDivH),'-neigh',num2str(Nneigh),'-rad',...
        num2str(maxRad)];
  mkdir([pathsave,'/',path])
end % function MakeDirectory

% ----------------------------------------------------------------------

function [forceR,forceA,neigh]=DefineForce(rad,cell,StiffRep,StiffAdh,...
Ncell,maxRad)

   forceR=zeros(Ncell,3); forceA=zeros(Ncell,3); neigh=zeros(Ncell,1);     
   for ii=1:Ncell
     for jj=ii+1:Ncell
       dx=cell(jj,1)-cell(ii,1);
       dy=cell(jj,2)-cell(ii,2);    
       dz=cell(jj,3)-cell(ii,3);
       dxyz=sqrt(dx*dx+dy*dy+dz*dz);
       if (dxyz>0)&&(dxyz<(rad(jj,1)+rad(ii,1)))  
 	     forceR(ii,1:3)=forceR(ii,1:3)+(StiffRep*(dxyz-(rad(ii,1)+...
                        rad(jj,1)))/dxyz)*[dx,dy,dz];  
         forceR(jj,1:3)=forceR(jj,1:3)-(StiffRep*(dxyz-(rad(ii,1)+...
                        rad(jj,1)))/dxyz)*[dx,dy,dz];  
       end
	   if (dxyz<2.25*maxRad)&&(dxyz>2*maxRad)
         forceA(ii,1:3)=forceA(ii,1:3)+(StiffAdh*(dxyz-2*maxRad)/dxyz)*[dx,dy,dz];
         forceA(jj,1:3)=forceA(jj,1:3)-(StiffAdh*(dxyz-2*maxRad)/dxyz)*[dx,dy,dz];
       end
	   if (dxyz>0)&&(dxyz<4*maxRad)
		 neigh(ii,1)=neigh(ii,1)+1;
         neigh(jj,1)=neigh(jj,1)+1;
       end
     end % for jj
   end % for ii  
end % end of function DefineForce

% ----------------------------------------------------------------------

function [maxar,maxd]=maxDiameter(cell,maxRad,maxd_old)
   distanceMxy=pdist2(cell(:,1:2),cell(:,1:2),'euclidean');
   distanceMyz=pdist2(cell(:,2:3),cell(:,2:3),'euclidean');
   distanceMxz=pdist2(cell(:,1:2:3),cell(:,1:2:3),'euclidean');
 
   uniqueMxy=unique(distanceMxy);
   uniqueMyz=unique(distanceMyz);
   uniqueMxz=unique(distanceMxz);

   maxd1=uniqueMxy(end,end);     % diameter from xy plane
   maxd2=uniqueMyz(end,end);     % diameter from yz plane
   maxd3=uniqueMxz(end,end);     % diameter from xz plane
   
   maxd=mean([maxd1,maxd2])+2*maxRad;  
  
   if maxd<maxd_old maxd=maxd_old; end
   maxar = pi*(maxd1/2)*(maxd2/2);

end % end of function maxDiameter

% ----------------------------------------------------------------------

function DrawSpheroid(cell,Ncell,rad,xmin,xmax,ymin,ymax,zmin,zmax,iter,...
Niter,cellLine,microEnv,drug,concDrug,ageDivH,Nneigh,cell_phase)
   set(gca,'color','white')
   set(gcf,'color','white')
   colordef white
   set(gca,'fontsize',14)
   axis([xmin,xmax,ymin,ymax,zmin,zmax])
   axis equal
   axis([xmin,xmax,ymin,ymax,zmin,zmax])
   hold on
   grid  
   
   % define cell line name
   if (cellLine==1)     text='-AT1';  
       colphase=[1 1 1;0.8 1 0.4;0.5 0.8 0.4;0.3 0.7 0.4;0.2 0.5 0.4];
   elseif (cellLine==2) text='-CA1a';
       colphase=[1 1 1;1 0.9 1;1 0.7 1;1 0.3 1;1 0.1 1];
   elseif (cellLine==3) text='-CA1d'; 
       colphase=[1 1 1;1 1 0.6;1 0.8 0.2;0.9 0.5 0.3;1 0.1 0.1];
   else                 text='A';     
       colphase=[1,1,1;0.5,1,1;0.1,0.5,1;0.1,0.3,1;0.1 0.1 1]; 
   end
   
   
   title({['CELL LINE: MCF10',text,'     iter=',num2str(iter),... 
      ' [',num2str(Niter),']'];...
      ['cell line: #',num2str(cellLine),', micoEnv: #',num2str(microEnv),...
      ', drug: #',num2str(drug),', conc: ',num2str(concDrug),'ng/mL'];...
      ['div age: ',num2str(ageDivH),';    neigh #',num2str(Nneigh)]},...
      'FontSize',15)  

     
   [xsph,ysph,zsph]=sphere(35); 
   for ii=1:Ncell
     crad=1.3*rad(ii,1);
     surf(cell(ii,1)+crad*xsph,cell(ii,2)+crad*ysph,cell(ii,3)+crad*zsph,...
          'facecolor',colphase(cell_phase(ii,2)+1,1:3),'edgecolor','none')
     hold on
   end % for Ncell
   view(-80+iter/10,20)
   alpha(0.6)
   light('Position',[-1 -1 1]*300, 'Style','infinite')
 
  
   % draw cell nuclei
   th=0:0.1:2*pi+0.2;
   for kk=1:Ncell
     xx=cell(kk,1)+crad*cos(th)*1.5;
     yy=cell(kk,2)+crad*sin(th)*1.5;
     fill3(xx,yy,zmin*ones(size(th,2),1),[1,1,1]*0.75,'EdgeColor','none');
   end
       
end  % end of function DrawSpheroid

% ----------------------------------------------------------------------

function [RgyrSpheroid,RgyrNorm] = rgyration(cell,Ncell,maxdiam,Niter)
  distOfCell = 0;
  Rgyr      = 0;
  centerSph = [sum(cell(1:Ncell,1)),sum(cell(1:Ncell,2)),...
               sum(cell(1:Ncell,3))] / Ncell;

  for nn=1:Ncell
	distOfCell(nn)=sqrt(...
	     (cell(nn,1)-centerSph(1,1))^2 +...
		 (cell(nn,2)-centerSph(1,2))^2 +...
		 (cell(nn,3)-centerSph(1,3))^2 );
  end
  meand = mean(distOfCell);

  Rg=distOfCell;
  Rgsq=Rg.*Rg;
  RgyrSpheroid = sqrt((1/Ncell) * sum(Rgsq));
  RgyrSphere   = sqrt((3/5) * (maxdiam(Niter)/2)^2);
  RgyrNorm = RgyrSpheroid / RgyrSphere;
end

% ----------------------------------------------------------------------
  
function r2 = rsquare(y,f,varargin)
  if isempty(varargin); c = true; 
  elseif length(varargin)>1; error 'Too many input arguments';
  elseif ~islogical(varargin{1}); error 'C must be logical (TRUE||FALSE)'
  else c = varargin{1}; 
  end
  
% Compare inputs
  if ~all(size(y)==size(f)); error 'Y and F must be the same size'; end

% Check for NaN
  tmp = ~or(isnan(y),isnan(f));
  y = y(tmp);
  f = f(tmp);
  
  n=length(y);
  sx=sum(y);  sy=sum(f);  % elements of corr coef calculations
  sxy=sum(y.*f);  sxx=sum(y.*y);  syy=sum(f.*f);

  r=(sxy-sx*sy/n) / sqrt((sxx-(sx*sx)/n)*(syy-(sy*sy)/n) );
  r2=r*r;
end

% =======================================================================