Generation of granule cell dendritic morphology (Schneider et al. 2014)

 Download zip file 
Help downloading and running models
Accession:167638
The following code was used to generate a complete population of 1.2 million granule cell dendritic morphologies within a realistic three-dimensional context. These generated dendritic morphologies match the known biological variability and context-dependence of morphological features.
Reference:
1 . Schneider CJ, Cuntz H, Soltesz I (2014) Linking macroscopic with microscopic neuroanatomy using synthetic neuronal populations. PLoS Comput Biol 10:e1003921 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Dendrite;
Brain Region(s)/Organism:
Cell Type(s): Dentate gyrus granule GLU cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Detailed Neuronal Models;
Implementer(s):
Search NeuronDB for information about:  Dentate gyrus granule GLU cell;
% Analyze Stats Based on Position in Granule Cell Layer

function dentate_9_p_analysis(input,directory)

global trees
trees       = {};
subset_size = 1000;

% Load Somata
load('Outputs/Somata.mat')
load('Outputs/Somata_pos.mat')


% Load trees and angles associated with input
load(sprintf('%s/Angles/%s.mat',directory,input))
load_tree(sprintf('%s/Trees_Tapered/%s.mtr',directory,input))
trees       = trees{1};
input2      = str2double(input);
start_soma  = 1+(input2-1)*subset_size;
if input2*subset_size < length(Somata)
    end_soma  = input2*subset_size;
else
    end_soma  = length(Somata);
end
range           = start_soma:end_soma;    

% Preallocate variables
Stats = struct;

% Get stats from TREES toolbox
s = stats_tree(trees,[],[],'-x');

% Rotate trees and get spanning information
trees_rotated   = rotate_created(trees,angles{1},angles{2}); 
[spanning,~]    = gscale_tree2(trees_rotated,' ');

% Calculate dendritic segments for each tree and append to matrix
% Allocate variables
Stats.type          = uint8(Somata_pos(range,1:2));
Stats.n_dend        = zeros(length(trees),1,'uint8');
Stats.n_branches    = zeros(length(trees),1,'uint8');
Stats.max_BO        = zeros(length(trees),1,'uint8');
Stats.BO            = cell(length(trees),1);
Stats.tip_dist      = zeros(length(trees),1);
Stats.tip_path      = zeros(length(trees),1);
Stats.sholl         = cell(length(trees),1);
Stats.asym          = zeros(length(trees),1);
Stats.inter_plen    = zeros(length(trees),1);
Stats.term_plen     = zeros(length(trees),1);
Stats.percents      = zeros(length(trees),16,'single');
Stats.length        = single(s.gstats.len);
Stats.spread_trans  = single(spanning.xlims{1});
Stats.height        = single(spanning.ylims{1});
Stats.spread_long   = single(spanning.zlims{1});
BO_all              = cell(length(trees),1);
B                   = cell(4,1);
T                   = cell(4,1);
X                   = cell(4,1);
Y                   = cell(4,1);
Z                   = cell(4,1);
S1                  = cell(4,1);
all_B               = cell(4,1);
all_T               = cell(4,1);
B_pts               = cell(length(trees),1);
T_pts               = cell(length(trees),1);
B{1}                = cell(length(trees),1);
B{2}                = cell(length(trees),1);
B{3}                = cell(length(trees),1);
B{4}                = cell(length(trees),1);
T{1}                = cell(length(trees),1);
T{2}                = cell(length(trees),1);
T{3}                = cell(length(trees),1);
T{4}                = cell(length(trees),1);
eucdist             = cell(length(trees),1);
pathdist            = cell(length(trees),1);
contraction         = cell(length(trees),1);
branch_indices      = cell(length(trees),1);

for sublayer = 1:4
    % Define surface for GCL and ML thirds
    if sublayer == 1
        [x_i,y_i,z_i]       = layer_eq_GCL(-1.95);
        [x_m,y_m,z_m]       = layer_eq_GCL(-1.0);
        [x_o,y_o,z_o]       = layer_eq_GCL(0);
        X{sublayer}         = [x_o;x_m;x_i];
        Y{sublayer}         = [y_o;y_m;y_i];
        Z{sublayer}         = [z_o;z_m;z_i];
        [~,S1{sublayer}]    = alphavol([X{sublayer}(:),Y{sublayer}(:),Z{sublayer}(:)],120);    
    elseif sublayer > 1
        [x_i,y_i,z_i]       = layer_eq_ML(sublayer-2);
        [x_o,y_o,z_o]       = layer_eq_ML(sublayer-1);
        X{sublayer}         = [x_o;x_i];
        Y{sublayer}         = [y_o;y_i];
        Z{sublayer}         = [z_o;z_i];
        [~,S1{sublayer}]    = alphavol([X{sublayer}(:),Y{sublayer}(:),Z{sublayer}(:)],150);
    end
end

% Find statistics for individual cells
for i = 1:length(trees)
    % Define topological points and remove soma branch points
    Soma        = Somata(range(i),:);
    iB          = B_tree(trees{i});
    iT          = T_tree(trees{i}); 
    B_pts{i}    = [trees{i}.X(iB) trees{i}.Y(iB) trees{i}.Z(iB)];
    T_pts{i}    = [trees{i}.X(iT) trees{i}.Y(iT) trees{i}.Z(iT)];

    % Find asymmetry
    asym                = asym_tree(trees{i});
    Stats.asym(i)       = mean(asym(iB));
    
    % Remove soma as a branch point
    if ismember(Soma,B_pts{i},'rows') > 0
        iB(1,1)     = 0;
        B_pts{i}(1,:)  = [];
    end
    iBT     = iB | iT;
    
    Stats.n_dend(i,1)     = size(find(PL_tree(trees{i})==1),1);
    Stats.n_branches(i,1) = size(B_pts{i},1)*2 + Stats.n_dend(i,1);
               
    % Branch order
    BO_all{i}           = BO_tree(trees{i});
    Stats.BO{i}         = uint8(BO_all{i}(iBT));
    Stats.max_BO(i,1)   = max(Stats.BO{i})+1;
    
    % Find euclidean and path distance to tips
    eucl                = eucl_tree(trees{i});
    Plen                = Pvec_tree(trees{i});
    Stats.tip_dist(i)   = mean(eucl(iT));
    Stats.tip_path(i)   = mean(Plen(iT));
    
    % Find sholl distribution
    diam                = 50;
    Stats.sholl{i}      = uint8(sholl_tree(trees{i},diam));
    
    % Break down tree into branches
    [sect,~]        = dissect_tree(trees{i});
    len             = len_tree(trees{i});
    sect            = sect(sect(:,1) ~= sect(:,2),:);
    n_branches      = length(sect);
    % Get indices for each branch
    ipar            = ipar_tree(trees{i});    
    
    % Allocate variables
    eucdist{i}           = zeros(n_branches,1);
    pathdist{i}          = zeros(n_branches,1);
    contraction{i}       = zeros(n_branches,1,'single');
    branch_indices{i}    = zeros(n_branches,1);
    
    for j = 1:n_branches
        % Determine all nodes for each branch from ipar
        start               = sect(j,1);
        stop                = sect(j,2);
        
        parent_line         = ipar(stop,:);
        stop_index          = find(parent_line == start);
        if isempty(stop_index) == 1
            stop_index = find(parent_line > 0,1,'last')+1;
        end
        branch_indices{j}   = transpose(parent_line(1:stop_index-1));
        
        % Get euclidean distance from branch beginning and end
        startxyz                = [trees{i}.X(start) trees{i}.Y(start) trees{i}.Z(start)];    
        stopxyz                 = [trees{i}.X(stop) trees{i}.Y(stop) trees{i}.Z(stop)];    
        eucdist{i}(j,1)         = pdist([startxyz;stopxyz],'euclidean');
        pathdist{i}(j,1)        = sum(len(branch_indices{j}));
        contraction{i}(j,1)     = eucdist{i}(j,1)/pathdist{i}(j,1);
    end    
    
    % Get intermediate and terminal branch lengths
    sect                = dissect_tree(trees{i});
    blen                = diff (Plen(sect), [], 2);
    [~,term_i]          = intersect(sect(:,2),find(iT));
    [~,inter_i]         = setdiff(sect(:,2),find(iT));
    inter_plen1         = blen(inter_i);
    Stats.inter_plen(i) = mean(inter_plen1(inter_plen1>0.01));
    Stats.term_plen(i)  = mean(blen(term_i));
end

all_B_pts = vertcat(B_pts{:});
all_T_pts = vertcat(T_pts{:});

% Determine which layer branch and termination points are in
for sublayer = 1:4
    IN_B            = inpolyhedron(S1{sublayer}.bnd,[X{sublayer},Y{sublayer},Z{sublayer}],all_B_pts,'FLIPNORMALS','True');
    IN_T            = inpolyhedron(S1{sublayer}.bnd,[X{sublayer},Y{sublayer},Z{sublayer}],all_T_pts,'FLIPNORMALS','True');
    all_B{sublayer} = all_B_pts(IN_B,:);
    all_T{sublayer} = all_T_pts(IN_T,:);
end

% Write out branch and termination laminar distribution percentages
for i = 1:length(trees)
    for sublayer = 1:4
       B{sublayer}{i} = intersect(B_pts{i},all_B{sublayer},'rows');
       T{sublayer}{i} = intersect(T_pts{i},all_T{sublayer},'rows');
    end
    B_Total = size(B{1}{i},1) + size(B{2}{i},1) + size(B{3}{i},1) + size(B{4}{i},1);
    T_Total = size(T{1}{i},1) + size(T{2}{i},1) + size(T{3}{i},1) + size(T{4}{i},1);

    % Write out percentages to matrix
    Stats.percents(i,:) = [100*size(B{1}{i},1)/B_Total 100*size(B{2}{i},1)/B_Total 100*size(B{3}{i},1)/B_Total 100*size(B{4}{i},1)/B_Total...
                     100*size(T{1}{i},1)/B_Total 100*size(T{2}{i},1)/T_Total 100*size(T{3}{i},1)/T_Total 100*size(T{4}{i},1)/T_Total...
                     100*(size(B{1}{i},1)+size(T{1}{i},1))/(B_Total+T_Total) 100*(size(B{2}{i},1)+size(T{2}{i},1))/(B_Total+T_Total) 100*(size(B{3}{i},1)+size(T{3}{i},1))/(B_Total+T_Total) 100*(size(B{4}{i},1)+size(T{4}{i},1))/(B_Total+T_Total)...
                     size(B{1}{i},1)+size(T{1}{i},1) size(B{2}{i},1)+size(T{2}{i},1) size(B{3}{i},1)+size(T{3}{i},1) size(B{4}{i},1)+size(T{4}{i},1)];
end

% Save outputs
save(sprintf('%s/Stats_Tapered/%s.mat',directory,input),'Stats','-v7.3')
save(sprintf('%s/Contraction/%s.mat',directory,input),'contraction','-v7.3')