
    % SIMULATION OF LINEAGE GENERATION THROUGH LAYER-BY-LAYER CELL DIVISIONS 
    % Llorca et al, 2019
    % written by Miguel Maravall, University of Sussex
    % m.maravall@sussex.ac.uk for info

    clear all
    close all

    % PRINCIPAL PARAMETERS 
    
    nreps = 100; 
    nclones_master = 100; 
    nlayers = 4; 
    capacity = 12; % determined by the maximum observed clonal size
                   % (each clone can have as many neurons as its capacity permits)
    
    % VALUES FROM EXPERIMENTAL DATA

    clone_sizes_Emx = [6     8     7     7     4     8     9     5     8     6     3     4    11     4     8     4    12     8      8     6     6     4     6     3    11     8     4    10     6     8     8     3     9     5     9     7     8     4     8     3     9     6     7    11     7     3    12     3    12     3     3     4     7     5     3     7     8     6    12     5     4     4     5    11     5     9     8    10     9     4     3     9     7    11     9    10     4     6    12     8     8    11     4     7     8     7     9     5     9     6     3     7    12     8    12     6     9     9    10     7     4     8     4]; 
     
    r_real = -.18; p_real = 0.066; 

    % distribution of Fucci and Emx counts (order: L6 L5 L4 L2/3)
    f_Fref = [0.364607267 0.157029189 0.234043911 0.244319632]; % Fucci
    f_Eref = [0.2458333333 0.1486111111 0.3069444444 0.2986111111]; % Emx

    % SIMULATION PARAMETERS    
    
    % each clone gets a chance to populate each of a number of slots per layer
    % number of slots per layer is given by nslots
    nslots = [7 6 10 9]; % in order: L6 L5 L4 L2/3
    % (determined based on maximum count observed experimentally)
    % Noppstot = sum(nslots); Noppstot = Noppstot(length(Noppstot)); 

    % probability of division, per layer
    % below is the vector of values used in final simulation data for 2 rules
    probval = 1.35*[0.3 0.2 0.33 0.38]; % in order: L6 L5 L4 L2/3
    % determined starting from a uniform probability (= 0.25 per layer) and empirically tuning
    % (optimal for 3rd probmaster rule with parameters +.4, .7, .8)

    % modulation of probability - this parameter introduces variants of the
    % main probability of division rule
    % probmaster = .15*(2*rand(1)-1) % to introduce anticorrelation between layers

    % INITIALIZATION OF VALUES
    
    probdiv = zeros(nlayers,nclones_master); 
    nclones_reps = zeros(1,nreps); 
    clone_sizes_total = []; 
    clone_sizes_repno = []; 
    clones_all_total = [];
    clone_sizes_max_reps = zeros(1,nreps); 
    df_reps = zeros(nreps,4); 
    fs_modelo_Fucci_reps = zeros(1,nreps); 
    fs_modelo_Emx_reps = zeros(1,nreps);
    r_reps = zeros(1,nreps);
    p_reps = zeros(1,nreps);
    frac_UPPER_reps = zeros(1,nreps);
    mean_neurons_UPPER_reps = zeros(1,nreps);
    frac_LOWER_reps = zeros(1,nreps);
    mean_neurons_LOWER_reps = zeros(1,nreps);
    frac_TRANSLAMINAR_reps = zeros(1,nreps);
    frac_ALL_LAYER_reps = zeros(1,nreps);
    frac_L5_EXCL_reps = zeros(1,nreps);
    frac_L6_EXCL_reps = zeros(1,nreps);
    frac_L4_EXCL_reps = zeros(1,nreps);
    frac_L2_EXCL_reps = zeros(1,nreps);
    frac_L2_L5_EXCL_reps = zeros(1,nreps);
    frac_L4_L6_EXCL_reps = zeros(1,nreps);
    frac_L4_L5_EXCL_reps = zeros(1,nreps);
    frac_L2_L6_EXCL_reps = zeros(1,nreps);
    tot_clones_cells_layer = zeros(nlayers,nreps*nclones_master);
     
    % % if all clones are equal:
    % probdiv(1:4,1) = probval; 
    % probdiv(1,2:nclones) = probdiv(1,1);
    % probdiv(2,2:nclones) = probdiv(2,1);
    % probdiv(3,2:nclones) = probdiv(3,1);
    % probdiv(4,2:nclones) = probdiv(4,1);

    % SIMULATION LOOP STARTS HERE
    
    for irep = 1:nreps; % each repeat of a population of progenitors of size nclones 
        
        % initialize
        nclones = nclones_master; 
        clones.cells_layer = zeros(nlayers,nclones);
        clones.ncells = ones(1,nclones);

        % if probability varies with clone but its distribution per layer is uniform across clones:
        for iclone = 1:nclones;
            
                probdiv(:,iclone) = probval;
            
        % introduce variants of the division rule, modulating probability 

        %%%%% first variant: each clone has +/- modulation of div probability across layers, introducing an effective anticorrelation

    %            probmaster = .25*(2*rand(1)-1); 
                % probmaster = .1*(2*rand(1)-1);
    %             probmaster = 0; 
    %             probdiv(1,iclone) = probval(1)+probmaster;
    %             probdiv(2,iclone) = probval(2)+probmaster;
    %             probdiv(3,iclone) = probval(3)-probmaster;
    %             probdiv(4,iclone) = probval(4)-probmaster;
                
        %%%%% second variant: some clones with modulation in superficial layers, some in deep layers; again, effective anticorrelation

%                 probmaster = .18*rand(1);
%                 if (rand(1) > 0.5)
%                     probdiv(1,iclone) = probval(1)-probmaster;
%                     probdiv(2,iclone) = probval(2)-probmaster;
%                     probdiv(3,iclone) = probval(3);
%                     probdiv(4,iclone) = probval(4);
%                 else
%                     probdiv(1,iclone) = probval(1);
%                     probdiv(2,iclone) = probval(2);
%                     probdiv(3,iclone) = probval(3)-probmaster;
%                     probdiv(4,iclone) = probval(4)-probmaster;
%                 end;

        %%%%% second variant, another version: sometimes + modulation, sometimes -
    %             
    %             probmaster = .3*(2*rand(1)-1);
    % 
    %             probdiv(1,iclone) = probval(1)-probmaster;
    %             probdiv(2,iclone) = probval(2)-probmaster;
    %             probdiv(3,iclone) = probval(3);
    %             probdiv(4,iclone) = probval(4);
    %             
    %         %%%%%

            %    probdiv(:,iclone) = probval*(0.5+1*rand(1));

        %%%%% third variant: one category of clones simply has the same probability vector for all, another category has its probability vector modulated 
        
%                             if (rand(1) > 0.35)
%             %                     if (rand(1) > 0.85)
%             %                         fact = 1; 
%             %                     else
%             %                         fact = .3+.3*rand(1); % (2*rand(1)-1);
%             %                         if (fact < 0.05) fact = 0.05; end;
%             %                     end;  
%                                 fact = .55+.45*(2*rand(1)-1);
%                                 probdiv(1,iclone) = fact*probval(1);
%                                 probdiv(2,iclone) = fact*probval(2);
%                                 probdiv(3,iclone) = fact*probval(3);
%                                 probdiv(4,iclone) = fact*probval(4);
%                             else
%                                 probmaster = .6*rand(1);
%                                 if (rand(1) > 0.85)
%                                     probdiv(1,iclone) = probval(1)-probmaster;
%                                     probdiv(2,iclone) = probval(2)-probmaster;
%                                     probdiv(3,iclone) = probval(3);
%                                     probdiv(4,iclone) = probval(4);
%                                 else
%                                     probdiv(1,iclone) = probval(1);
%                                     probdiv(2,iclone) = probval(2);
%                                     probdiv(3,iclone) = probval(3)-probmaster;
%                                     probdiv(4,iclone) = probval(4)-probmaster;
%                                 end;
%                             end;
% 

        %%%%% fourth variant: one homogeneous population and another to which a modulation rule (both +/-) applies

            %    probmaster = .14*(1.8*rand(1)-1);
    %         probmaster = .2*(2*rand(1)-1);
    %         if (rand(1) > 0.35)
    %             probdiv(1,iclone) = probval(1);
    %             probdiv(2,iclone) = probval(2);
    %             probdiv(3,iclone) = probval(3);
    %             probdiv(4,iclone) = probval(4);
    %         else
    %             probdiv(1,iclone) = probval(1)+probmaster;
    %             probdiv(2,iclone) = probval(2)+probmaster;
    %             probdiv(3,iclone) = probval(3)-probmaster;
    %             probdiv(4,iclone) = probval(4)-probmaster;
    %         end;
    
        %%%%% fourth variant, second version: one homogeneous population and another for which modulation follows a gradient across layers
        
%             probmaster = .2*(2*rand(1)-1);
%             if (rand(1) > 0.4)
%                 probdiv(1,iclone) = probval(1);
%                 probdiv(2,iclone) = probval(2);
%                 probdiv(3,iclone) = probval(3);
%                 probdiv(4,iclone) = probval(4);
%             else
%                 for ilayer = 1:nlayers;
%                     probdiv(ilayer,iclone) = probval(ilayer)-probmaster*ilayer+probmaster*4;
%                 end; 
%             end;

        %%%%% fifth variant: modulation set uniformly for all layers

%                 probmaster = .4*rand(1);
%                 fact = 1; 
%                 if (rand(1) > 0.7)
%                         probdiv(1,iclone) = fact*probval(1);
%                         probdiv(2,iclone) = fact*probval(2);
%                         probdiv(3,iclone) = fact*probval(3);
%                         probdiv(4,iclone) = fact*probval(4);
%                 else
%                         probdiv(1,iclone) = probval(1)-probmaster;
%                         probdiv(2,iclone) = probval(2)-probmaster;
%                         probdiv(3,iclone) = probval(3)-probmaster;
%                         probdiv(4,iclone) = probval(4)-probmaster;
%                 end;

            %     %%%%%
            
               
            % finally: 
            % restricted populations where modulation of the probability of
            % division occurs in specific ways
            
            if (rand(1) > 0.9)
                capacity = 4;
                nslots = [7 6 4 8];
                for islots = 1:4; nopps(islots) = 1 + round(nslots(islots)*rand(1,1)); end;
                
                probdiv(1,iclone) =  0.05;
                %if probdiv(1,iclone) < 0; probdiv(1,iclone) = 0; end;
                probdiv(2,iclone) =  0.01;
                %if probdiv(2,iclone) < 0; probdiv(2,iclone) = 0; end;
                probdiv(3,iclone) = 0.6;
                probdiv(4,iclone) = 0.7;
                
            else
                if (rand(1) > 0.93)
                    capacity = 4;
                    nslots = [5 6 10 9];
                    for islots = 1:4; nopps(islots) = 1 + round(nslots(islots)*rand(1,1)); end;
                    probdiv(1,iclone) =  0.4;
                    probdiv(2,iclone) =  0.6;
                    probdiv(3,iclone) = 0.001;
                    probdiv(4,iclone) = 0.001;
                    
                else
                    capacity = 12;
                    nslots = [7 6 10 9];
                    for islots = 1:4; nopps(islots) = 1 + round(nslots(islots)*rand(1,1)); end;
                    
                end;
            end;

        % cells per clone:
        % all clones can have as many cells as "capacity" allows
        clones.ncells(iclone) = capacity; 
        % OR
        % each clone is limited by its own capacity (with max capacity set as
        % above)
        % clones.ncells = round(5*randn(1,nclones)+14);
        % clones.ncells = zeros(1,nclones);
        % for iclones = 1:nclones;
        %     nsearch = find(clonesize_cdf > rand(1));
        %     clones.ncells(iclones) = clonesize_vals(nsearch(1));
        % end;

        % GENERATE CLONES

            actualclonesize = 0; 
            
            while actualclonesize < 3; % sets minimal clone size
                
                clonesize = 0;
                
                % here define stochasticity of transitions between layers, according to a capacity (nopps, number of division opportunities per layer) which is either
                % fixed or stochastic
                % 
                % nopps = nslots % means fixed
                % for islots = 1:3; nopps(islots) = round(nslots(islots)+0.8*randn(1,1)); end; nopps(4) = Noppstot - (nopps(1)+nopps(2)+nopps(3)); nopps
                %for islots = 1:4; nopps(islots) = 1 + round(nslots(islots)*rand(1,1)); end; 
                %nopps = nslots;
                
                death = 0; 
                
                for ilayer = 1:nlayers;
                    cellgenerator = 0;
                    
                    if death > 0; break; end; % this refers to "premature" clone death
%                     if clonesize >= clones.ncells(iclone); break; end; 
                      
                    for islots = 1:nopps(ilayer);
                        
                        if rand(1) < 0.015; 
                            death = 1; 
                            break; 
                        end; 
                        
                        if rand(1) < probdiv(ilayer,iclone);
                            cellgenerator = cellgenerator + 1;
                            clonesize = clonesize + 1;
                        end;
                        
%                          if clonesize >= clones.ncells(iclone);
%                             break;
%                          end;
                         
                    end;                    
                    clones.cells_layer(ilayer,iclone) = cellgenerator;
                end;
                
                actualclonesize = clonesize; 
                if actualclonesize > capacity; 
                    actualclonesize = 0;
                    clones.cells_layer(:,iclone) = zeros(nlayers,1); 
                end; 
                
                end; 

        end; % end of clone generation

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

        % set up parameters for analysis
        clone_sizes = cumsum(clones.cells_layer);
        clone_sizes = clone_sizes(end,:);
        clones_ind = find(clone_sizes > 2);
        nclones = length(clones_ind);
        clone_sizes = clone_sizes(clones_ind);
        clone_sizes_max_reps(irep) = max(clone_sizes);
        clone_sizes_total = [clone_sizes_total clone_sizes]; 
        clone_sizes_repno = [clone_sizes_repno clone_sizes*0+irep]; 
        clones.cells_layer = clones.cells_layer(:,clones_ind);
        
        tot_clones_cells_layer(:,(irep-1)*nclones+1:irep*nclones) = clones.cells_layer; 
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % CALCULATION OF MEASURES FOR COMPARISON TO EXPERIMENTAL DATA
        
        % ------ Fractional distribution of neurons per layer; scores for fits to expt data ------
        distrib = cumsum(clones.cells_layer,2);
        SUMA = cumsum(distrib(:,end)); df = distrib(:,end)'/SUMA(end); display(['MODEL = ' num2str(df) ' (L6 L5 L4 L2/3)']); 
        df_reps(irep,:) = df; 
       %  fs_modelo_Fucci_reps(irep) = (f_Fref-df)*(f_Fref-df)';
        fs_modelo_Emx_reps(irep) = (f_Eref-df)*(f_Eref-df)';

        % ------ Anticorrelation in cell counts per layer ------
        cell_numbers_56 = cumsum(clones.cells_layer(1:2,:),1); cell_numbers_56 = cell_numbers_56(2,:);
        cell_numbers_24 = cumsum(clones.cells_layer(3:4,:),1); cell_numbers_24 = cell_numbers_24(2,:);
        [r,p] = corr(cell_numbers_24',cell_numbers_56','type','Spearman');
        r_reps(irep) = r;
        p_reps(irep) = p;

        % ------ Distribution of clonal types ------

        % clones with no neurons in L5 & L6 (UPPER)
        indclones_UPPER = find(clones.cells_layer(1,:)+clones.cells_layer(2,:) < 1);
        l_UPPER = length(indclones_UPPER);
        frac_UPPER_reps(irep) = l_UPPER/nclones;
        nums_UPPER = clones.cells_layer(3,indclones_UPPER)+clones.cells_layer(4,indclones_UPPER);
        nums_UPPER = nums_UPPER(find(nums_UPPER > 2));
        mean_neurons_UPPER_reps(irep) = mean(nums_UPPER);

        % clones with no neurons in L2/3 & L4 (LOWER)
        indclones_LOWER = find(clones.cells_layer(3,:)+clones.cells_layer(4,:) < 1);
        l_LOWER = length(indclones_LOWER);
        frac_LOWER_reps(irep) = l_LOWER/nclones;
        nums_LOWER = clones.cells_layer(1,indclones_LOWER)+clones.cells_layer(2,indclones_LOWER);
        nums_LOWER = nums_LOWER(find(nums_LOWER > 2));
        mean_neurons_LOWER_reps(irep) = mean(nums_LOWER);

        % clones with cells in all layers (TRANSLAMINAR)
        frac_TRANSLAMINAR_reps(irep) = (nclones - length(find(clones.cells_layer(1,:)+clones.cells_layer(2,:) < 1)) - length(find(clones.cells_layer(3,:)+clones.cells_layer(4,:) < 1)))/nclones;

        frac_ALL_LAYER_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*clones.cells_layer(2,:).*clones.cells_layer(3,:).*clones.cells_layer(4,:) > 0));

        frac_L5_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*(clones.cells_layer(2,:)-1).*clones.cells_layer(3,:).*clones.cells_layer(4,:) < 0));
        frac_L6_EXCL_reps(irep) = 1/nclones*length(find( (clones.cells_layer(1,:)-1).*clones.cells_layer(2,:).*clones.cells_layer(3,:).*clones.cells_layer(4,:) < 0));
        frac_L4_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*clones.cells_layer(2,:).*(clones.cells_layer(3,:)-1).*clones.cells_layer(4,:) < 0));
        frac_L2_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*clones.cells_layer(2,:).*clones.cells_layer(3,:).*(clones.cells_layer(4,:)-1) < 0));

        frac_L2_L5_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*clones.cells_layer(3,:).*(clones.cells_layer(2,:)+clones.cells_layer(4,:)-1) < 0));
        frac_L4_L6_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(2,:).*clones.cells_layer(4,:).*(clones.cells_layer(1,:)+clones.cells_layer(3,:)-1) < 0));
        frac_L4_L5_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(1,:).*clones.cells_layer(4,:).*(clones.cells_layer(2,:)+clones.cells_layer(3,:)-1) < 0));
        frac_L2_L6_EXCL_reps(irep) = 1/nclones*length(find(clones.cells_layer(2,:).*clones.cells_layer(3,:).*(clones.cells_layer(1,:)+clones.cells_layer(4,:)-1) < 0));

        ind_all = find(clones.cells_layer(1,:).*clones.cells_layer(2,:).*clones.cells_layer(3,:).*clones.cells_layer(4,:) > 0);
        clones_all = clones.cells_layer(:,ind_all);
        clones_all_total = [clones_all_total clones_all];

    end; 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % DISPLAY OF RESULTS
    
    display(' ')
    display(' ')
    display('------ Distribution of clonal types ------')
    display(' ')
    display(['clonal size (MODEL) = ' num2str(mean(clone_sizes_total)) ' (mean), ' num2str(std(clone_sizes_total)) ' (std), ' num2str(max(clone_sizes_total)) ' (max)']) 
    display(['mean max clonal size (MODEL) = ' num2str(mean(clone_sizes_max_reps)) ' (mean max clonal size in a simulation)']) 
    display(['clonal size (expt from Emx) = ' num2str(mean(clone_sizes_Emx)) ' (mean), ' num2str(std(clone_sizes_Emx)) ' (std), ' num2str(max(clone_sizes_Emx)) ' (max)']) 
    figure; 
    hold on
    %possible_sizes = min(clone_sizes_total):1:max(clone_sizes_total); 
    possible_sizes = min(clone_sizes_total):1:12; 
    histogram(clone_sizes_total,possible_sizes,'Normalization','probability','DisplayStyle','stairs'); 
    histogram(clone_sizes_Emx,possible_sizes,'Normalization','probability','DisplayStyle','stairs');
    title('Distribution of clonal sizes, simul and expt')
    [h,p] = kstest2(clone_sizes_total,clone_sizes_Emx)

    frac_clone_size_distrib = zeros(nreps,length(possible_sizes)); 
    for irep = 1:nreps; 
        clone_sizes_irep = find(clone_sizes_repno == irep);
        x = hist(clone_sizes_total(clone_sizes_irep),possible_sizes);
        xc = cumsum(x);
        frac_clone_size_distrib(irep,:) = x/xc(end);
    end; 
    clear x xc; 

    xEmx = hist(clone_sizes_Emx,possible_sizes); 
    xcEmx = cumsum(xEmx); 
    frac_clone_size_distrib_Emx = xEmx/xcEmx(end); 
    x = cumsum((frac_clone_size_distrib_Emx-mean(frac_clone_size_distrib)).*(frac_clone_size_distrib_Emx-mean(frac_clone_size_distrib))./var(frac_clone_size_distrib)); 
    display(' ')
    display(['z-score distrib clonal sizes = ' num2str(sqrt(1/length(possible_sizes)*x(end))) ])
    display(' ')

    figure(10)
    subplot(1,2,1)
    y = [frac_clone_size_distrib_Emx ; mean(frac_clone_size_distrib)]; 
    bar(y')
    ax = gca; 
    ax.XTickLabel = 3:1:12;
    % legend('observed','simulated')
    hold on
    errorbar(1.15:1:10.15,mean(frac_clone_size_distrib),0*std(frac_clone_size_distrib),std(frac_clone_size_distrib),'.k')
    xlim([0 11])


    display(' ')
    display(' ')
    display('------ Fractional distribution of neurons per layer; scores for fits to expt data ------')
    display(' ')
    display(' (numbers for somatosensory cortex) ')
    display(' ')
    display('EXP (Emx) =            0.24583     0.14861     0.30694     0.29861    (L6 L5 L4 L2/3)')
    display('EXP (Fucci complete) = 0.36461     0.15703     0.23404     0.24432    (L6 L5 L4 L2/3)')
    display(' ')
    display(['MODEL (mean) =        ' num2str(mean(df_reps)) '    (L6 L5 L4 L2/3)'])
    display(['Score, fit to MODEL (model vs Emx) : ' num2str(mean(fs_modelo_Emx_reps)) ' (mean), ' num2str(min(fs_modelo_Emx_reps)) ' (min), '  num2str(max(fs_modelo_Emx_reps)) ' (max)'])
    display(' ')
    fs_Emx_vs_Fucci = (f_Eref-f_Fref)*(f_Eref-f_Fref)'; 
    display(['Score, fit to EXP (Emx vs Fucci complete) = ' num2str(fs_Emx_vs_Fucci)])
    figure
    hist(fs_modelo_Emx_reps); 
    hold on; 
    plot([fs_Emx_vs_Fucci fs_Emx_vs_Fucci],[0 nreps/2],'--','color',[.5 .5 .5]) 
    title('Histogram of scores (model vs Emx; --: Emx vs Fucci)') 

    x = cumsum((f_Eref-mean(df_reps)).*(f_Eref-mean(df_reps))./var(df_reps)); 
    display(' ')
    display(['z-score distrib by layer = ' num2str(sqrt(1/4*x(end))) ])
    display(' ')

    figure(10)
    subplot(1,2,2)
    y = [mean(df_reps) ; f_Eref]; 
    bar(y')
    ax = gca; 
    ax.XTickLabel = {'VI','V','IV','II-III'};
    % legend('observed','simulated')
    hold on
    errorbar([.85,1.85,2.85,3.85],mean(df_reps),0*std(df_reps),std(df_reps),'.k')
    xlim([0 5])
    view([90 -90])

    display(' ')
    display(' ')
    display('------ Anticorrelation between counts per layer ------')
    display(' ')
    display(' (see figure)')
    figure;
    subplot(1,2,1)
    hold on 
    title('Spearman (anti)correlation (r)')
    for irep = 1:nreps; 
        plot(r_reps(irep),irep,'k.'); 
    end; 
    plot([r_real r_real],[1 nreps],'--','color',[.5 .5 .5])
    subplot(1,2,2)
    hold on 
    title('p-value')
    for irep = 1:nreps; 
        plot(p_reps(irep),irep,'k.'); 
    end; 
    plot([0.05 0.05],[1 nreps],'--','color',[1 0 0])
    plot([p_real p_real],[1 nreps],'--','color',[.5 .5 .5])
    ax = gca; 
    ax.XScale = 'log';

    display(' ')
    display(['z-score (anti)correlation = ' num2str(abs(r_real-mean(r_reps))/std(r_reps))])
    display(' ')


    frac_real_1 = [0.116504854	0.106796117 0.7767]; 
    bardata_1 = [mean(frac_UPPER_reps) mean(frac_LOWER_reps) mean(frac_TRANSLAMINAR_reps)];
    barerror_1 = [std(frac_UPPER_reps) std(frac_LOWER_reps) std(frac_TRANSLAMINAR_reps)];
    figure; 
    bar([1 2 3],bardata_1,'r')
    hold on
    errorbar([1 2 3],bardata_1,[0 0 0],barerror_1,'k.')
    plot([1 2 3],frac_real_1,'k*')
    title('Upper-only, lower-only, and translaminar clones')
    ax = gca;
    ax.XTickLabel = {'UPPER','LOWER','TRANS'};

    frac_real_2 = [0.223300971 0.203883495 0.067961165 0.048543689 0.077669903 0.019417476 0.077669903 0.009708738 0.048543689]; 
    bardata_2 = [mean(frac_ALL_LAYER_reps) mean(frac_L5_EXCL_reps) mean(frac_L6_EXCL_reps) mean(frac_L4_EXCL_reps) mean(frac_L2_L5_EXCL_reps) mean(frac_L4_L6_EXCL_reps) mean(frac_L2_EXCL_reps) mean(frac_L4_L5_EXCL_reps) mean(frac_L2_L6_EXCL_reps)];
    barerror_2 = [std(frac_ALL_LAYER_reps) std(frac_L5_EXCL_reps) std(frac_L6_EXCL_reps) std(frac_L4_EXCL_reps) std(frac_L2_L5_EXCL_reps) std(frac_L4_L6_EXCL_reps) std(frac_L2_EXCL_reps) std(frac_L4_L5_EXCL_reps) std(frac_L2_L6_EXCL_reps)];
    figure; 
    bar([1 2 3 4 5 6 7 8 9],bardata_2,'b')
    hold on
    errorbar([1 2 3 4 5 6 7 8 9],bardata_2,[0 0 0 0 0 0 0 0 0],barerror_2,'k.')
    plot([1 2 3 4 5 6 7 8 9],frac_real_2,'k*')
    title('Layer configurations for translaminar clones')
    ax = gca;
    ax.XTickLabel = {'ALL','L2L4L6','L2L4L5','L2L5L6','L4L6','L2L5','L4L5L6','L2L6','L4L5'};


    % distribution by clonal type
    figure
    hold on
    boxplot([frac_UPPER_reps' frac_LOWER_reps' frac_TRANSLAMINAR_reps'],'colors',[.5 .5 .5],'symbol','k.');
    plot([1 2 3],frac_real_1,'k*')
    
    figure
    hold on
    boxplot(                    [frac_ALL_LAYER_reps'   frac_L5_EXCL_reps'  frac_L4_EXCL_reps'  frac_L6_EXCL_reps'  frac_L4_L6_EXCL_reps'   frac_L2_L5_EXCL_reps'   frac_L2_EXCL_reps'  frac_L4_L5_EXCL_reps'   frac_L2_L6_EXCL_reps'],'colors',[.5 .5 .5],'symbol','k.');    
    plot([1 2 3 4 5 6 7 8 9],   [0.223300971            0.203883495         0.048543689         0.067961165         0.019417476             0.077669903             0.077669903         0.009708738             0.048543689],'k*')

    
    nsize = nreps*nclones_master; 
    % distributions of neurons by layer
    sim_sizes_l23 = tot_clones_cells_layer(4,:);
    sim_sizes_l4 = tot_clones_cells_layer(3,:);
    sim_sizes_l5 = tot_clones_cells_layer(2,:);
    sim_sizes_l6 = tot_clones_cells_layer(1,:);
    uiopen('neuron counts per layer.fig',1)
    subplot(2,4,1)
    hold on
    plot(0:1:12,histc(sim_sizes_l23,-.5:1:11.5)*103/nsize,'g')
    subplot(2,4,2)
    hold on
    plot(0:1:12,histc(sim_sizes_l4,-.5:1:11.5)*103/nsize,'g')
    subplot(2,4,3)
    hold on
    plot(0:1:12,histc(sim_sizes_l5,-.5:1:11.5)*103/nsize,'g')
    subplot(2,4,4)
    hold on
    plot(0:1:12,histc(sim_sizes_l6,-.5:1:11.5)*103/nsize,'g')

    
    
    