Example 3: Spiny stellates in barrel cortex, but partially generated
====================================================================

This is a more realistic model of thalamocortical projections
to L4 spiny stellate cells. Thalamocortical projection axon denisties are 
again modelled as quasi randomly distributed spherical Gaussians.
Spiny Stellate cell densities on the other hand are now derived from 
real data, by taking six reconstructed spiny stellate cells and fitting 
a 14-component Gaussian mixture model to each of them.
Two versions of this model are used as examples. This version version
@@link@@ uses all realistic values, but only generates a part of the
barrel, The other limits the number of postsynaptic neurons and the number
of voxels to be able to generate a whole neuron in a whole barrel.

References:
[1] Dimensions of a projection column and architecture of VPM and POm axons in rat vibrissal cortex. Wimmer VC, Bruno RM, de Kock CP, Kuner T, Sakmann B. Cereb Cortex. 2010 Oct;20(10):2265-76
[2] Number and laminar distribution of neurons in a thalamocortical projection column of rat vibrissal cortex. Meyer HS, Wimmer VC, Oberlaender M, de Kock CP, Sakmann B, Helmstaedter M. Cereb Cortex. 2010 Oct;20(10):2277-86
[3] Cell type-specific thalamic innervation in a column of rat vibrissal cortex. Meyer HS, Wimmer VC, Hemberger M, Bruno RM, de Kock CP, Frick A, Sakmann B, Helmstaedter M. Cereb Cortex. 2010 Oct;20(10):2287-303
[4] Functional diversity of layer IV spiny neurons in rat somatosensory cortex: quantitative morphology of electrophysiologically characterized and biocytin labeled cells. Staiger JF, Flagmeyer I, Schubert D, Zilles K, Koetter R, Luhmann HJ.Cereb Cortex. 2004 Jun;14(6):690-701

Define lpar
-----------

.. code-block:: matlab

  clear lpar
  lpar.nDim = 3;
  lpar.nSigRange = 2 ;
  
  lpar.vox.voxelSize = 1; %voxel size is set to three micron
  lpar.vox.voxPerDim = [370 370 263] ; %L4 of one barrel is about this size [2] (in number of voxels)
  
  %in this script, we will only generate a part of the barrel, therefore we
  %define the area (in world coordinates) which to generate
  area = [100 150;100 150;100 150] ;
  
  %definition of a neural population
  %Seed for this population
  postPop.seed = 27552 ;
  postPop.name = 'Spiny Stellate Dendrites' ;
  postPop.nr = 1 ;
  postPop.isPre = 0 ;
  
  postPop.Nn = 2374 ; %should be about 2374 [2]&[4]
  postPop.soma.edge = 0.3 ;
  postPop.soma.distr = 'qrand' ; 
  postPop.soma.seed = 2421 ;
  
  %Set the typeSpec from fitted neurons. This also requires changing Yinfo
  %after prepareForYgen. (see below). The mat file can be generated from
  %neuron reconstructions using the trees toolbox to load and resample them,
  %and MATLABS gmdistribution.fit to generate the fit. 
  [postPop.typeSpec, wVecC]= loadTypeSpecFromMatFile('GMfit_BlueBrain_MB') ;
  %make sure to use all neurons, for more information see 'help randSample'
  postPop.typeSpec.neuronNumberDistr = 'alli' ;
  
  %Trees balancing factor
  postPop.mst.bf = [.55 0.025] ;
  postPop.mst.bfDistr = 'normal' ;
  %Trees max branch threshold
  postPop.mst.thr = [100] ;
  
  lpar.neurPop{1} = postPop ;
  % define second population
  prePop.isPre = 1 ;
  prePop.Nn = 50 ; %@@Actual number not available
  prePop.seed = 18713 ;
  prePop.name = 'VPM projection neurons' ; %[1] ventral postmedial axons project to the barrels
  prePop.soma.edge = 0.1 ;
  prePop.soma.distr = 'qrand' ; %quasi random, like random, but more uniform
  prePop.soma.seed = 294312;
  
  prePop.typeSpec.generator = @(x) genGauss(x); %genDen1(x);, genGauss(x);
  prePop.typeSpec.sig = [20 40] ;
  
  %Trees balancing factor
  prePop.mst.bf = [.55 0.025] ;
  prePop.mst.bfDistr = 'normal' ;
  %Trees max branch threshold
  prePop.mst.thr = [100] ;
  
  lpar.neurPop{2} = prePop ;
  
  [lpar,Yinfo] = prepareForYgen(lpar) ;
  
  % Because of the design of the program, the weight vector from fits should be set after prepareForYgen
  newWvec = arrayfun(@(x) wVecC{1}{x},Yinfo{1}.sampled.typeSpec.neuronNumber,'uni',0);
  Yinfo{1}.sampled.typeSpec.wVec = cell2mat(newWvec);

Generate everything (Y, X & connectivity matrix)
------------------------------------------------

.. code-block:: matlab

  [lpar,Yinfo] = preparePartialY(lpar,Yinfo) ;
  [YprePart,YpostPart,neurInd] = generatePartialY(lpar,Yinfo,area) ;
  
  nSynPreMean = 5000 ;
  nSynPreSig = 300 ;
  nSynPostMean = 300 ;
  nSynPostSig = 70 ;
  removeOverlap = true ;
  TreeJitterAmplitude = 0.3 ;
  
  disp('Generating presynaptic elements') ;
  [XprePart] = genRand1(YprePart,nSynPreMean,nSynPreSig,removeOverlap) ;
  disp('Generating postsynaptic elements') ;
  [XpostPart] = genRand1(YpostPart,nSynPostMean,nSynPostSig,removeOverlap) ;

Growing trees in partially generated worlds is currently not possible

.. code-block:: matlab

  % disp('Generating presynaptic trees') ;
  % [ TreesPre,synBOpre ] = growTrees(lpar,XprePart,Yinfo,true,TreeJitterAmplitude) ;
  % disp('Generating postsynaptic trees') ;
  % [ TreesPost,synBOpost ] = growTrees(lpar,XpostPart,Yinfo,false,TreeJitterAmplitude) ;
  % disp('Generating Connectivity matrix') ;
  
  [conMat,synNeurMap] = generateConnectivity(XprePart,XpostPart) ;
  % Calculate measures and visualize some things
  
  disp('Calculating measures') ;
  figure(1)
  Ymsrs = Ymeasures(YprePart,YpostPart,true)

.. image:: _images/part_1.png
  :width: 60%

Selected Y measures. These graphs quickly show wether the generated neurons
have a significant portion of their density in the generated part of the
world. Because only a fraction of the whole space is generated, only a
fraction of the densities are within this area, yet the fraction of voxels
occupied is significantly higher.

.. code-block:: matlab

  XmsrsPre = Xmeasures(lpar,XprePart) 
  XmsrsPost = Xmeasures(lpar,XpostPart) 
  
  figure(2)
  somaColor = [1 0 0]; 
  %Somasizes depend on number of synaptic elements
  showSoma(lpar,Yinfo,false,1,somaColor)
  %Other visualizations for partially generated densities are under construction
  hold on ;
  axis equal
  title('Postsynaptic Soma');

.. image:: _images/part_2.png
  :width: 60%

All Soma