Example 1: Long range, zeroth order pyramidal density model =========================================================== The Matlab script :mat:func:`Examples/LongRangeZerothOrder.m` examplifies the typical workflow of the synapse generation pipeline. Before running it, make sure that your Matlab path points to the synGenPipeline folders, by running :mat:func:`pathme.m`. The main steps of LongRangeZerothOrder are: 1. Setting the model parameters (lpar). 2. Generating axonal and dendritic densities in a discrete space. 3. Placing a synapses and assigning them to neurons. 4. Generate axonal and dendritic trees 5. Visualize results Implemented is a very crude model of Thalamocortical projections on L4 pyramidal neurons. Pyramidal dendrite densities are given by two Gaussians, rotation symmetric around the vertical axis: 1. A wide flat Gaussian to represent basal dendrites; 2. A tall and narrow Gaussian to represent apical dendrites. Together their shape can second best be described as a lava lamp. Thalamocortical axon arborizations are modelled as quasi-randomly placed spherical Gaussians. The pyramidal densities are placed upright in the neuropil, with the orientation of their apical dendrites perturbed by an angle normally distributed around zero with a standard deviation of 0.15 rad. Initializing ------------ .. code-block:: matlab clear variables % if not already done: add synGenPipeline functions to the Matlab path run(['..' filesep 'pathme']) Model settings (lpar) --------------------- .. code-block:: matlab %number of dimensions lpar.nDim = 3; %cutoff [number of sigma's] for generating target distribution lpar.nSigRange = 3 ; %voxel size and number lpar.vox.voxelSize = 1; lpar.vox.voxPerDim = 150 ; %Defining the first post-synaptic neural population (lpar) %Seed for this population postPop.seed = 27552 ; %Identification name postPop.name = 'Pyramidal cells' ; %Identification number postPop.nr = 1 ; %whether or not this is a presynaptic population postPop.isPre = 0 ; %Number of neurons in this population postPop.Nn = 20 ; %Where these soma exist, in normalized coordinates postPop.soma.area = [.2 .8;.2 .8;.1 .7]; %soma distribution function postPop.soma.distr = 'rand' ; %'qrand' 'unifRand' ,'unif' ||(debug: 'cen' 'diag') %Seed for the distribution of the soma postPop.soma.seed = 2421 ; %anonymous function to the blob postPop.typeSpec.generator = @(x) genPyr0(x); %See 'help blobs' %paramters for the generator postPop.typeSpec.tsz = 20 ; postPop.typeSpec.tsr = 0.4; postPop.typeSpec.bsxyr = 25; postPop.typeSpec.bsr = 0.15 ; postPop.typeSpec.dr = 1.55 ; postPop.typeSpec.f = 0.3 ; postPop.typeSpec.wVec = [0.5,1.5] ; %Rotation angles, perturb around the vertical axis is a random direction %Rotations are applied in the order x, y, then z. postPop.rot.X = 0 ; postPop.rot.Y = [0,0.15] ; postPop.rot.Z = [0,pi] ; %rotation distribution See 'help randSample' postPop.rot.YDistr = 'normRand' ; postPop.rot.ZDistr = 'unifRand' ; %Trees balancing factor. See 'help MST_tree' postPop.mst.bf = 0.55 ; %Trees max branch threshold postPop.mst.thr = 100 ; %Assign postPop to the neurPop lpar.neurPop{1} = postPop ; %Define the second population (lpar) prePop.isPre = true ; prePop.Nn = 30 ; prePop.seed = 18713 ; prePop.name = 'Thalamocortical axons' ; %Identifier %To specify the area where soma reside, you can also define the size of %the edge where there are no soma. This is, like area, in normalized %coordinates, so edge [0,0,5) ; prePop.soma.edge = 0.1 ; prePop.soma.distr = 'qrand' ; %'qrand' 'unifRand' ,'unif' ||(debug: 'cen' 'diag') prePop.soma.seed = 294312; prePop.typeSpec.generator = @(x) genGauss(x); %See 'help blobs' %parameter for this generator prePop.typeSpec.sig = [25 35] ; prePop.typeSpec.sigDistr = 'unifRand' ; %Trees balancing factor. See 'help MST_tree' prePop.mst.bf = [0.4 0.6 ] ; %Because there is no definition of bfDistr, unifRand will be used by default %Trees max branch threshold.See 'help MST_tree' prePop.mst.thr = 1000 ; %assign lpar.neurPop{2} = prePop ; Generating the target density (lpar -> Y) ----------------------------------------- .. code-block:: matlab %First step, fixing lpar (e.g. translating soma.edge to soma.area) and %sampling all stuff for Yinfo disp('Preparation') ; [lpar,Yinfo] = prepareForYgen(lpar) ; %Generates the target density disp('Generating target densities') ; [Ypre,Ypost] = generateY(lpar,Yinfo) ; Generates synaptic elements using the genRand1 Method. (Y->X) ------------------------------------------------------------- The synaptic element generator (Y->X) that is used here is genRand1. This method needs a specified number of synaptic elements, an aimed in-population standard deviation. The method is able to remove overlapping synaptic elements. The following variables are used as inputs to the function genRand1 .. code-block:: matlab nSynPreMean = 3000 ; nSynPreSig = sqrt(3000) ; nSynPostMean = 300; nSynPostSig = 40 ; removeOverlap = true ; disp('Generating presynaptic elements') ; [Xpre] = genRand1(Ypre,nSynPreMean,nSynPreSig,removeOverlap) ; disp('Generating postsynaptic elements') ; [Xpost] = genRand1(Ypost,nSynPostMean,nSynPostSig,removeOverlap) ; Growing the trees using the TREES toolbox' mimimum spanning tree (X->Tree) -------------------------------------------------------------------------- The TreeJitterAmplitude is used in growTrees. The generated synaptic elements are placed somewhere in the near the center of their voxel, normally distributed in 3 dimensions independently, with a standard deviation given by the TreeJitterAmplitude. Grow trees is a gateway function to the trees toolbox' MST_tree. It generates branches between the generated synaptic elements in X by creating a weighted minimum spanning tree. The output if this function is a cell array of trees, and list which associates a branch order with each synaptic element. .. code-block:: matlab disp('Generating presynaptic trees') ; TreeJitterAmplitude = 0.3 ; [ TreesPre,synBOpre ] = growTrees(lpar,Xpre,Yinfo,true,TreeJitterAmplitude*4) ; disp('Generating postsynaptic trees') ; [ TreesPost,synBOpost ] = growTrees(lpar,Xpost,Yinfo,false,TreeJitterAmplitude) ; Generating the connectivity matrix (X -> Conmat) ------------------------------------------------ A connectivity matrix is generated from two sets of synatic elements. Everywhere where a presynaptic and a postsynaptic element occur in a voxel, a synapse is formed. The function returns a connectivity matrix a cell matrix containing the voxel numbers of the synapses between two neurons .. code-block:: matlab disp('Generating Connectivity matrix') ; [conMat,synNeurMap] = generateConnectivity(Xpre,Xpost) ; Calculating measures, for more information see the respecitve help (e.g. type:"help Xmeasures") ----------------------------------------------------------------------------------------------- .. code-block:: matlab disp('Calculating and showing measures') ; Ymsrs = Ymeasures(Ypre,Ypost) XmsrsPre = Xmeasures(lpar,Xpre) XmsrsPost = Xmeasures(lpar,Xpost) ; nBins = 10 ; CMmsrs = conMatMeasures(conMat,nBins,false); Branch order distribution ------------------------- Associates branch orders to synapses by comparing the synaptic element branch order association in the synBO variable, to the the neuron-synapse association in the synNeurMap .. code-block:: matlab disp('Calculating branch order distributions') ; boSynListPre = branchOrderDistr(synBOpre,synNeurMap,true) ; boSynListPost = branchOrderDistr(synBOpost,synNeurMap,false) ; %Show the branch order histogram figure(1) allBO = cell2mat(boSynListPost) ; hist(allBO,1:max(allBO)) ; title('branch order histrogram for all postsynaptic neurons'); xlabel('Branch Order') ; ylabel('Number of synapses') ; .. image:: _images/lrzo_fig1.png :width: 60% Histogram with the postsynaptic branch orders of all synapses. Calculate projections --------------------- .. code-block:: matlab disp('Calculating Projections') ; %Target density projections YprojPre = densityProjection(lpar,Ypre,true) ; % YprojPost = densityProjection(lpar,Ypre,true) ; %synaptic element projections XprojPre = densityProjection(lpar,Xpre) ; % XprojPost = densityProjection(lpar,Xpost) ; % Calculating soma projections %Soma projection resolution (in voxels) SomaDensResolution = 5 ; %Soma density projections, SomaDensResolution is defined above somDensPre = somaDens3D(lpar,Yinfo(lpar.preIndex),SomaDensResolution) ; somDensPost = somaDens3D(lpar,Yinfo(lpar.postIndex),SomaDensResolution) ; %Special one, synaptic element projectiosn, yse genSynX to generate a %matrix of X's shape from the synNeurMap synProjPre = densityProjection(lpar,genSynX(lpar,synNeurMap,[],[],true)) ; Plot result figures ------------------- .. code-block:: matlab %Showing the connectivity matrix figure(2) showConMat(conMat,CMmsrs) ; .. image:: _images/lrzo_fig2.png :width: 60% Connectivity matrix and associated measures. Each square has a color based on the number of synapses between the associated two neurons. The Degree is the number of connected neurons, and the strength is the total number of synapses. In and out correspond to the postsynaptic and presynaptic population respecitively. The multiplicity is the number of synapses between any two neurons, consequently the last histogram is a histogram over all the colors of the squares. .. code-block:: matlab %Showing Projections figure(3) ; showProjections(lpar,YprojPre,XprojPre,synProjPre) ; .. image:: _images/lrzo_fig3.png :width: 60% Projections (not slices) of the target density, the genereated synaptic elements and the generated synapses of the presynaptic population. .. code-block:: matlab figure(4); showProjections(lpar,Ypost,Xpost) .. image:: _images/lrzo_fig4.png :width: 60% Projections (not slices) of the target density, the genereated synaptic elements of the postsynaptic population. .. code-block:: matlab %Showing soma density projections figure(5) subplot(1,3,1) ; imagesc(somDensPre.X); xlabel('Z (bins)'); ylabel('Y (bins)'); title('Soma') ; subplot(1,3,2) ; imagesc(somDensPre.Y); xlabel('X (bins)'); ylabel('Z (bins)'); title('Soma') ; subplot(1,3,3) ; imagesc(somDensPre.Z); xlabel('Y (bins)'); ylabel('X (bins)'); title('Soma') ; .. image:: _images/lrzo_fig5.png :width: 60% Projected soma locations .. code-block:: matlab % RecenterPlot figure(6) recenterPlot(lpar,Xpost,Yinfo{lpar.postIndex}) ; title('Synaptic elements of all postsynaptic neurons, relative to soma and unrotated') ; .. image:: _images/lrzo_fig6.png :width: 60% Recenterplot shows the synaptic elements of one population. Unlike showY, the synaptic elements are not plotted by their location in world coordinates, but relative to the soma, with the axes of the generator aligned. Visualize two neurons in full detail ------------------------------------ The following sections uses several plotting functions starting with show, namely, showY, showSoma,showTrees and showSynapses. It is left as an excersize for the reader to figure out what their function is. Another function that belongs to this list is showX. For More information see the respective help files. .. code-block:: matlab %for what neurons do we show everything whichNeuronExtensively = [2 8] ; figure(7);clf %soma plot color somaColor = [0 0 1] ; showY(lpar,Ypost(:,whichNeuronExtensively)) ; hold on ; %show alls soma showSoma(lpar,Yinfo,false,[],somaColor,'.') %Plotting functions clear the axes if hold is off %show some trees showTrees(TreesPost(whichNeuronExtensively),[],'-p') axis equal %Synapse colormap otherCc= [1 1 1; somaColor;distinguishable_colors(numel(whichNeuronExtensively))]; synapseCC = distinguishable_colors(lpar.postN,otherCc) ; %add the synapses showSynapses(lpar,synNeurMap,[],whichNeuronExtensively,true,synapseCC) ; .. image:: _images/lrzo_fig6.png :width: 60% Two neurons in their full glory; Soma, Target density, generated tree (where each corner is a generated synaptic element), and synapses. The synapse colors represent different presynaptic partners.