the discovery of grid cells in the medial entorhinal cortex ( mec ) ( hafting et  al . , 2005 ) has been a major advance in cortical physiology ( burgess 2014 ) .
the assessment of single - unit activity in rats running in boxes has led to the discovery of a plethora of  functional  cell types in the mec : conjunctive ( head - directional ) grid cells ( sargolini et  al . , 2006 ) , border cells ( solstad et  al . , 2008 ) , boundary vector cells ( koenig et  al . , 2011 ) , speed cells ( kropff et  al . , 2015 ) , and cue cells ( kinkhabwala et  al .
grid and border cells also exist in areas neighboring the entorhinal cortex , such as the subiculum and pre- and parasubiculum ( lever et  al . , 2009 , boccara et  al . ,
computational models propose many different mechanisms to explain how grid cell discharges come about ( giocomo et  al .
, 2011 , zilli , 2012 ) . a better knowledge of the anatomy and spatio - temporal firing patterns of defined cell types is needed to constrain models and help prune the forest of different models .
two aspects of the temporal firing patterns were highlighted in recent work : burstiness and theta cycle skipping .
burstiness has been shown to be associated with grid cell firing ( newman and hasselmo , 2014 , latuske et  al . , 2015 ) and might serve important functions in parahippocampal microcircuits ( welday et  al . , 2011 , sheffield and dombeck , 2015 ) .
burstiness has also been linked to differences in extracellular spike shape ( newman and hasselmo , 2014 , latuske et  al . , 2015 ) .
theta cycle skipping might be related to the computation of head - directional information and grid firing ( brandon et  al . , 2013 ) .
previous investigations of burstiness and theta cycle skipping have analyzed mixed extracellular recordings from both the superficial medial entorhinal cortex and the parasubiculum ( brandon et  al .
it has thus remained unclear whether burstiness and theta cycle skipping map onto anatomical categories or whether bursty and non - bursty neurons are simply intermingled ( latuske et  al . , 2015 ) .
stellate cells ( stel ) in layer 2 ( l2 ) of the medial entorhinal cortex show a tendency to fire bursts of action potentials upon membrane depolarization in  vitro ( alonso and klink , 1993 , pastoll et  al . , 2012 ,
such findings led to the hypothesis that stellate cells might display bursty firing patterns in  vivo ( newman and hasselmo , 2014 , latuske et  al . , 2015 ) .
entorhinal grid cells phase - precess ; i.e. , they shift spike timing in a systematic way relative to the field potential during firing field transversals ( hafting et  al . , 2008 , jeewajee et  al . , 2013 ,
, it has been found that mec l2 cells phase - precess more strongly than mec layer 3 ( l3 ) cells ( hafting et  al . , 2008 ,
this difference between mec layers 2 and 3 has not been seen at the single run level ; however , it may arise because mec l3 cells are less correlated between runs ( reifenstein et  al .
recently , a single run analysis of phase precession revealed differences between pyramidal and stellate neurons in mec l2 ( reifenstein et  al . , 2016 ) .
parasubicular neurons provide specific input to mec l2 pyramidal neurons ( pyr ) ( tang et  al . , 2016 ) , but it is unknown whether parasubicular neurons phase - precess . here we analyze juxtacellular recordings from the medial entorhinal cortex ( ray et  al . , 2014 , tang et  al .
, 2014a , tang et  al . , 2015 ) and the parasubiculum ( tang et  al . , 2016 ) .
juxtacellular data offer two advantages ( pinault , 1996 , herfst et  al . , 2012 ) .
second , juxtacellular recording of the local field potential ( lfp ) and spikes has a very high temporal resolution and signal - to - noise ratio , which is crucial for investigating temporal patterns such as burstiness .
does burstiness differ between parasubicular neurons , mec l2 pyramids , mec l2 stellates , and mec l3 neurons ?
l2 of the mec contains patches of calbindin - positive pyramidal neurons arranged in a hexagonal grid ( ray et  al . , 2014 ;
figure  1a , top ) that are surrounded by calbindin - negative stellate cells ( figure  1a , top , black background ) .
the parasubiculum ( pas ) is a thin elongated structure that wraps around the mec mediodorsally and has high wolframin expression ( wfs1-positive cells ; figure  1a , top ;
tang et  al . , 2016 ) . axons from the parasubiculum specifically target the patches of mec l2 pyramidal cells ( burgalossi et  al .
mec l3 neurons are not arranged in a hexagonal grid but are visible as a homogenous band of purkinje cell protein 4 ( pcp4)-positive cells below layer 2 ( l3 ; red band in figure  1a , bottom ) .
figure  1b , left , shows reconstructions of example cells of the four neuron types : a parasubicular neuron ( blue ) , a mec l2 pyramidal neuron ( green ) , a mec l2 stellate cell ( black ) , and an mec l3 pyramidal neuron ( red ) , all recorded in freely moving rats .
all reconstructions are from tangential sections ( i.e. , a  top view  of the morphology ) .
in addition to the morphology , we also show juxtacellular recording traces from the reconstructed example cells ( figure  1b , right ) .
two signals are visible in the recordings : the spikes of the identified cells and the prominent theta rhythm in the lfp . to determine whether a neuron was discharging in a bursty pattern , we analyzed the interspike interval ( isi ) histogram using a similar approach as latuske et  al .
isis below 60  ms were binned in 2-ms bins ( normalized to area  = 1 to generate a probability distribution ) , which revealed that our dataset contained both non - bursty and bursty cells ( figure  2a ) .
we performed a principal component analysis on a matrix of the isi probability distributions of all neurons and found that the first three principal components ( pcs ; figure  2b , bottom ) explained 69% of the variance in the data . in agreement with latuske
( 2015 ) , we found that , when the first two principal components were plotted against each other , the neurons formed a c - shaped structure , indicative of a bimodal distribution ( figure  2b , top ) .
we assigned the neurons to two clusters using a k - means clustering algorithm on the first three principal components ( figure  2c , top ) .
the two clusters were well separated with little overlap ( figure  2d ) . to assess the separation quality of the two clusters , we calculated the projection of the neurons onto fisher s linear discriminant .
we can interpret the linear discriminant as a measure of  burstiness  because it is places the cells along an axis from non - bursty to bursty based on the shape of the isi histogram .
we plotted all cells sorted according to burstiness , and , in agreement with latuske et  al .
( 2015 ) , we found that bursty neurons were distinguished by a tendency to fire bursts at 125250  hz ( 4- to 8-ms bins ; figure  2d ) . to investigate differences in burstiness among cell types
, we plotted the median isi histogram of all recorded cells , resolved by cell type .
the median isi histograms of parasubicular as well as mec l2 pyramidal neurons indicated very bursty cells ( figure  3a , top ) .
the median isi histograms of mec l2 stellate and mec l3 neurons were flat with no obvious burstiness ( figure  3a , bottom ) . to assess whether this difference was statistically significant , we performed two tests : one based on categorical classifications of cells as 
non - bursty  and  bursty  with a guard zone ( experimental procedures ; latuske et  al . ,
when we compared the proportions of non - bursty , guard - zoned , and bursty cells among neuron types , we found no significant difference between parasubicular neurons and mec l2 pyramids , which both contained predominantly bursty cells ( pas versus pyr , bursty / guard / non - bursty : 11/11/0 versus 15/15/1 , p
we also found no difference between mec l2 stellate cells and mec l3 cells ( stel versus l3 , bursty / guard / non - bursty : 9/25/34 versus 3/5/24 , p > 0.05 ,  test ; figure  3b ) , which were both predominantly non - bursty . both parasubicular neurons and mec l2 pyramids contained significantly different proportions of bursty and non - bursty cells in comparison with both mec l2 stellates and mec l3 neurons ( all p  < 0.001 ,  tests ; figure  3b ) .
the width and placement of the guard zone is estimated from the bimodal fit , and thus the guard zone depends on the relative abundance of bursty and non - bursty cells , which is evidently not the same among neuron types ; i.e. , the guard zone might be either too wide or too narrow .
the guard zone also discards information telling us whether a neuron is near the guard zone or closer to the extremes .
these problems may inflate our estimated differences in burstiness among cell types . to make sure that no spurious results were imposed by the guard zone , we directly compared the burstiness of the neuron types and included all cells . in agreement with the estimations based on comparisons of the proportions , we found that the burstiness of parasubicular neurons and mec l2 pyramids was significantly higher than the burstiness in both mec l2 stellates and mec l3 neurons ( all p  < 0.001 , mann - whitney u tests ; figure  3c ) .
again , we did not find a significant difference between parasubicular neurons and mec l2 pyramids ( p > 0.05 , mann - whitney u test ; figure  3c ) , but we did find that mec l3 neurons had a significantly lower burstiness than mec l2 stellates ( p  = 0.0036 , mann - whitney u test ; figure  3c ) .
thus , parasubicular neurons and mec l2 pyramids are generally bursty , whereas mec l2 stellates and mec
l3 neurons are generally non - bursty ( figures 3a and 3b ) . furthermore , within the non - bursty neuron types , mec l3 neurons are more strictly non - bursty than mec l2 stellates ( figure  3c )
. it should be noted , however , that even though there are large and highly significant differences in burstiness among cell types , the distributions of burstiness among cell types are overlapping .
for example , a minority of l2 stellate cells and l3 neurons assume firing patterns that are otherwise classically parasubicular / pyramid - like .
our dataset includes mec l2 neurons that were classified as putatively pyramidal or stellate based on theta strength and preferred theta phase ( tang et  al .
we therefore also checked whether there was any correlation between burstiness and theta strength because such a correlation might introduce 
we fitted three generalized linear models to investigate whether burstiness might be related to theta strength ( model 1 , burstinessstrength ; figure  s3a , left ) , putative cell type ( model  2 , burstinesstype ; figure  s3a , middle ) , or both ( model  3 , burstinesstype  + strength ; figure  s3a , right ) .
both comparisons of the akaike information criterion ( aic ; akaike , 1974 ) and likelihood ratio tests of nested models indicated that model 2 is superior to the other models ( figure  s3b ) ; i.e. , the burstiness  depends only on putative cell type ( model 2 , ptype  = 0.0000076 ; figure  s3c , middle ) and not on theta strength ( model  2 versus model 3 , p  = 0.54 , likelihood ratio test ; figure  s3b ) .
second , we plotted the burstiness among cell types twice : once where we include the classified mec l2 cells ( figure  s3d , left ) and once where we only include identified mec l2 cells ( figure  s3d , right ) .
the pattern of burstiness among cell types remained the same when we only included the identified cells ( figure  s3d ) .
we thus conclude that cell type - specific differences in burstiness are not an artifact of our classification approach .
in tetrode recordings of parasubicular and mec l2/3 neurons , differences in spike shape have been linked to burstiness ( latuske et  al . , 2015 ) and theta phase preference of grid cells ( newman and hasselmo , 2014 ) .
we therefore investigated whether there was a difference in spike shape among our four anatomical categories of neurons .
first we removed a subset of cells for which the signal - to - noise ratio of spike waveforms was insufficient to reliably assess the spike shape .
second , we removed spikes that happened within 100  ms of the previous spike to disregard potential effects of spike shape adaptation during bursts ( experimental procedures ) . in figure  4a , we plot the remaining spike shapes ( normalized for display ; experimental procedures ) for all four neuron types .
we did not find any differences among neuron types in spike amplitude , peak - to - trough ratio , or spike half - width ( all p > 0.05 , kruskal - wallis test ) .
this was expected for two reasons : overall spike amplitude depends strongly on the particular recording pipette and relation to the  soma ( gold et  al . , 2009 ) , and narrow spikes and a small peak - to - trough ratio are indicative of interneurons ( mountcastle et  al .
we noticed , however , a large variability in the repolarization phase of the cell type : parasubicular neurons and mec l2 pyramids contained many cells that quickly reached the trough and repolarized , whereas mec l2 stellates and mec l3 neurons reached the trough more slowly ( figure  4a ) .
this tendency was also evident in the mean spike shape of the four neuron types ( figure  4b ) .
when we compared the peak - to - trough time of the cell types , we found significant differences ( p  = 0.0014 , kruskal - wallis test ) .
parasubicular neurons and mec l2 pyramids had significantly shorter peak - to - trough times than both mec l2 stellates and mec l3 neurons ( all p  < 0.05 , mann - whitney u test ; figure  4c ) . because parasubicular neurons and mec l2 pyramids have faster peak - to - trough times and are also the most bursty cell types , we wondered whether , as has been suggested ( latuske et  al . , 2015 ) , the burstiness of the cell predicts the spike shape .
alternatively , the spike shape might simply be different among neuron types , or it might depend on neuron types as well as burstiness . to figure this out , we decided to employ a generalized linear regression approach . because peak - to - trough time can not assume negative values , we modeled peak - to - trough time as a gamma - distributed variable ( experimental procedures ) .
we selected the appropriate model using the following approach : we first modeled peak - to - trough time as a function of only burstiness ( glm1 , peak - to - trough 1  + burstiness ) and found a significant dependence ( anova , pburstiness  = 0.0087 ; figure  s4a , dashed gray line ) .
( 2015 ) . also , when we modeled peak - to - trough time as a function of only neuron type ( glm2 , peak - to - trough 1  + type ) , we also found a significant dependence on neuron type ( anova , ptype  = 0.0015 ; figure  s4a , solid lines ) . however , when we modeled peak - to - trough time as a function of both burstiness and neuron type ( glm3 , peak - to - trough 1  + burstiness  + type ) , we found that the dependency on type but not the dependency on burstiness remained significant ( anova , pburstiness  = 0.22 , ptype  = 0.017 ; figure  s4c ) .
we also fitted a model where we allowed for interactions between burstiness and type ( glm4 , peak - to - trough 1  + burstiness  + type  + burstinesstype ) , where all effects became non - significant ( anova , all p > 0.05 ; figure  s4c ) . to determine which model
best explains the data , we calculated the aic of all models and found that , despite the four fitted parameters , glm2 had the lowest aic , indicating that the peak - to - trough time depends on neuron type , but not on burstiness ( figure  s4b ) . similarly , when comparing nested models , we found that glm3 better explains the data than glm1 ( p  = 0.0023 , likelihood ratio test ; figure  s4b ) ; i.e. , including neuron type as a predictor makes the model better
. we did not find that glm3 explains the data better than glm2 ( p  = 0.32 , likelihood ratio test ; figure  s4b ) ; i.e. , it is unnecessary to include burstiness as a predictor in addition to neuron type .
we thus infer that the differences in spike shape primarily reflect the anatomical type and not the burstiness of the neuron . to determine whether a neuron was theta cycle - skipping , we used a maximum likelihood estimation ( mle ) of a parametric model of the isi histogram ( climer et  al . , 2015 ; experimental procedures ) .
our dataset contained neurons that showed no theta modulation and also neurons that had strong rhythmic components ( figure  5a ) .
for every cell , we fitted three models to the isi distribution : a  flat  model with no rhythmic components ( figure  5a , left ) , a  rhythmic , non - skipping  model with a theta - rhythmic modulation of the isi histogram ( figure  5a , middle ) , and a  rhythmic , cycle - skipping  model with a theta - rhythmic modulation of the isi histogram and a second parameter introducing theta cycle skipping ( i.e. , a higher amplitude of every other peak in the isi histogram ; figure  5a , right ) .
the three fitted  models were compared using the appropriate  statistic ( calculated from the maximum log likelihood of the models ) to generate two p values : prhythmic ( comparing the flat and the rhythmic , non - skipping models ) and pskipping ( comparing the rhythmic , non - skipping and the rhythmic , cycle - skipping models ) .
the cells were classified using a two - level classification ( figure  5b ) : first , we determined whether a cell was  rhythmic  ( prhythmic  < 0.05 ) or  non - rhythmic  ( prhythmic > 0.05 ) .
then we classified the rhythmic cells as either rhythmic , cycle - skipping ( pskipping  < 0.05 ) or rhythmic , non - skipping ( pskipping > 0.05 ) . using the mle approach , we found that parasubicular neurons and mec l2 pyramids were overwhelmingly rhythmic ( 93% ; pas , 20/22 ; pyr , 29/31 ; figure  5c , left ) .
mec l2 stellates and mec l3 neurons were rarely rhythmic ( 26% ; stel , 16/68 ; l3 , 9/32 ) , both significantly less rhythmic than both parasubicular neurons and mec l2 pyramids ( all p  < 0.001 ,  tests ; figure  5c , left ) .
this is in agreement with previous observations in which evaluated spike train rhythmicity of cell types using a 
theta index  was used ( ray et  al . , 2014 , tang et  al .
we found that the generally rhythmic cell types were also significantly more likely to also be theta cycle - skipping than the generally non - rhythmic cell types ( p  = 0.018 , mixed - effects logistic regression ; figure  5c , right ; experimental procedures ) : approximately 49% of the rhythmic parasubicular neurons and rhythmic mec l2 pyramids were also theta cycle - skipping ( pas , 9/20 ; pyr , 15/29 ; figure  5c , right ) . of the mec l2 stellates and mec l3 neurons ,
which were classified as rhythmic using the mle approach , only 20% were also theta cycle - skipping ( stel , 4/16 ; l3 , 1/9 ; figure  5c , right ) .
our dataset includes mec l2 neurons that were classified as  putatively pyramidal or stellate based on theta strength and preferred theta phase ( tang et  al .
obviously , we expect a correlation between the theta rhythmicity ( which is calculated from the isi distribution ) and the theta strength ( locking to the lfp theta rhythm ) . however , the mle approach of climer et  al .
( 2015 ) returns a p value of the rhythmicity per cell and is sensitive to very low amounts of rhythmicity , which could potentially have been present in , e.g. , putative stellates with a low locking strength and locking to the peak of the lfp theta rhythm ( figure  s1 ; climer et  al . , 2015 , tang et  al . , 2014a ) .
more importantly , our classification procedure considers simply strength of locking to the local lfp , and there is no way of distinguishing a simply theta - rhythmic cell from a rhythmic and cycle - skipping cell based on theta strength because they might show equally strong locking . to be sure that the cell type differences were not an artifact of including the classified cells , we plotted the burstiness among cell types twice : once where we included the classified mec l2 cells ( figure  s5a , left ) and once where we only included identified mec l2 cells ( figure  s5a , right ) .
the proportions among cell types remained the same when restricting the analysis to identified cells only ( figure  s5a ) . to compare the magnitude of phase precession among cell types at the single - run level , we first selected single runs of high firing based on the firing rate ( figures 6a , top , and 6b ; experimental procedures ) . from these single runs , we determined the slope and range of phase precession by a circular - linear fit of time and theta phase angle of the spikes in each run ( figure  6a , bottom ; experimental procedures )
. figure  6c shows example single runs from example cells of the four neuron types .
the example mec l2 stellate and l3 neurons have steep phase precession slopes and cover larger ranges of theta phase angles during a single run .
in contrast , the example parasubicular neuron and mec l2 pyramid only weakly phase - precess . across the population
, we found the same result : first , identified and putative mec l2 stellate and l3 neurons had approximately 3-fold steeper phase precession slopes than parasubicular neurons and identified and putative mec l2 pyramids ( figure  6d ; median slopes : pas / pyr / stel / l3  = 16.7/25.9/76.7/64.8 degrees / s ; p(pas versus stel )  = 0.0001 ; p(pas versus l3 )  = 0.003 ; p(pyr versus stel )  = 0.0001 ; p(pyr versus l3 )  = 0.01 ; mann - whitney u tests ) .
second , identified and putative mec l2 stellate and l3 neurons covered a much larger range of theta phase angles per run than parasubicular neurons and identified and putative mec l2 pyramids ( approximately 2-fold ; figure  6e ;  median ranges : pas / pyr / stel / l3  = 63.2/48.7/127.5/114.2 degrees ; p(pas versus stel )  = 0.00008 ; p(pas versus l3 )  = 0.0007 ; p(pyr versus stel )  = 0.0000002 ; p(pyr versus l3 )  = 0.00005 ; mann - whitney u tests ) .
we did not find any differences in the circular - linear correlation coefficient among the cell types ( p  = 0.38 , kruskal - wallis test ) .
we used advanced statistical techniques to tease apart how  differences in burstiness , spike shape , theta modulation  ( rhythmicity , locking , skipping ) , and phase precession map onto regular spiking layer 3 medial entorhinal neurons , layer 2 medial entorhinal pyramidal neurons , layer 2 medial entorhinal stellate neurons , and parasubicular regular spiking cells .
we found significant differences in spike shape , burstiness , theta modulation ( rhythmicity , locking , cycle skipping , phase precession ) , and theta phase precession between the four groups of cells investigated .
thus , our data suggest that cell type is a major determinant of discharge patterns in the rat parasubiculum and superficial medial entorhinal cortex .
although our data emphasize the significance of cell types , the discharge patterns we observed do not directly match what is expected based on the analysis of intrinsic properties of these neurons in  vitro .
in  vitro recordings of parasubicular neurons have suggested an intrinsic disposition for theta rhythmicity ( glasgow and chapman , 2008 ) .
it is known that in  vitro measurements of l2 mec cell properties are very sensitive to recording conditions ( alonso and klink , 1993 , pastoll et  al . , 2012 )
. however , mec l2 stellates often display some intrinsic burstiness in  vitro ( alonso and klink , 1993 , pastoll et  al . , 2012 , alessi et  al . , 2016 , fuchs et  al . , 2016 ) , but they are generally not very bursty in  vivo ( ray et  al . , 2014 ; figure  2 ) .
thus , it is probably incorrect to assume that bursty cells recorded extracellularly in the superficial mec and the parasubiculum are mec l2 stellates ( newman and hasselmo , 2014 , latuske et  al . , 2015 ) because we show that bursty cells are more likely to be mec l2 pyramids or parasubicular neurons .
although phase precession is arguably the most intensely studied example of temporal coding in the brain , its underlying mechanism is still a matter of debate .
parasubicular neurons , which show only weak phase precession , project to pyramidal cells in mec l2 ( tang et  al . , 2016 ) .
also , these mec l2 pyramidal cells express only a low degree of phase precession .
conversely , stellate cells in mec l2 and pyramidal cells in mec l3 phase precess with steep slopes .
the latter finding is somewhat surprising because it challenges the long - held belief that cells in mec l3 do  not phase - precess ( hafting et  al . , 2008 ,
in contrast to that , we analyzed phase precession in single runs ( schmidt et  al . , 2009 ) .
we argue that the single - run approach is more appropriate because the animal needs to process information online and does not have the opportunity to pool over trials .
our finding of substantial mec l3 phase precession is in line with a previous single - run account ( reifenstein et  al . , 2014 ) .
mec l2 stellate cells project to the dentate gyrus , whereas mec l2 pyramidal cells send output to ca1 ( varga et  al .
, 2010 , kitamura et  al . , 2014 , ray et  al . , 2014 ) .
because mec l2 pyramidal cells show only weak phase precession , it seems unlikely that they substantially contribute to ca1 phase precession .
therefore , ca1 either generates phase precession de novo or inherits phase - precessing inputs via the strongly precessing stellate cells in mec l3 ( jaramillo et  al . , 2014 ) .
it is presently unknown how the functional categories ( grid cells , border cells , speed cells , cue cells , etc . )
for example , it is unknown whether mec l2 grid cells are predominantly pyramidal cells ( tang et  al . , 2014a ) or stellate cells ( domnisoru et  al . , 2013 ) or whether they show no preference for either cell type ( sun et  al . , 2015 )
similarly , some authors have reported that about a third to half of mec l3 neurons are grid cells ( sargolini et  al . , 2006 , boccara et  al . , 2010 ) , whereas others have estimated that if l3 grid cells exists , then they must be rare ( 1% ; tang et  al .
parasubicular neurons and mec l2 pyramids are more bursty , have narrower spikes , and are more likely to be theta - rhythmic , theta - locked , and theta cycle - skipping than mec l2 stellates and mec l3 neurons .
these differences remain even when we statistically control for interactions between spike shape , burstiness , and rhythmicity .
some studies have tried to elucidate the grid cell generation mechanism by characterizing the firing properties of the entorhinal network . from these studies we know that grid cells are bursty whereas border cells are not ( newman and hasselmo , 2014 , latuske et  al . ,
it has also been shown that theta cycle skipping is somehow necessary for maintaining grid cell firing ( brandon et  al . , 2013 ) .
in agreement with tang et  al . , 2014a , tang et  al . , 2016
, we conclude that , based on burstiness and theta cycle skipping , parasubicular neurons and mec l2 pyramids are likely to play a key role in generating grid cell activity in the parasubiculum and superficial medial entorhinal cortex . in line with the differences in temporal discharge patterns , we observed that parasubicular and mec l2 pyramidal cells had shorter spike durations than mec l2 stellates and mec l3 neurons .
several previous studies have noticed significant differences between mec l2 pyramidal and mec l2 stellate cells , most notably , that stellate cells have larger depolarizing afterpotentials ( alonso and klink , 1993 , alessi et  al .
in in  vivo recordings , it was generally observed that stellate cells had a shorter spike duration than pyramidal cells ( alonso and klink , 1993 ) .
interestingly , however , it was also found that the spike duration of both pyramidal and stellate cells varied depending on the depolarizing current pulse ( alonso and klink , 1993 ) .
thus , the juxtacellularly observed differences in spike shape are probably not primarily a reflection of differences in intrinsic cell properties .
however , the distributions of spike durations are largely overlapping ( figure  4c ) , probably precluding a classification of extracellularly recorded mec l2 regular spiking neurons into pyramidal and stellate cells based purely on spike shape .
we observed a strong similarity between spike shapes and firing patterns of parasubicular neurons and mec l2 pyramids .
these two neuron groups were different in spike shapes and firing patterns from layer 3 cells and layer 2 stellate cells , which were similar to each other , however .
parasubicular axons specifically target patches of mec l2 pyramidal cells ( tang et  al . , 2016 ) , which might be a pathway for head - directional information from the medial septum to reach the grid cell system ( winter et  al . , 2015 , unal et  al . , 2015 , tang et  al . ,
l3 cells and layer 2 stellate cells provide a massive direct ( l3 ) and indirect input to the hippocampus , whereas projections from both layer 2 pyramids and the parasubiculum are minor or absent ( varga et  al . , 2010 ,
thus , analysis of spike shapes and firing patterns , direct connectivity , and projection targets supports the distinction of layer 3 cells and layer 2 stellate cells on one hand and parasubiculum and layer 2 pyramids on the other hand as two distinct processing systems .
they both express wolframin ( ray and brecht , 2016 ) , and , in the early development stages , they also express calbindin ( ray and brecht , 2016 ) .
likewise , layer 3 neurons and layer 2 stellate cells also have an anatomical likeliness in their protein expression profile , with both expressing reelin in adult rats ( ray and brecht , 2016 ) .
this might allude to the electrophysiological and functional characteristics of these two groups being perhaps somewhat genetically determined , with the protein expression profiles of these respective cell groups shaping their inputs and outputs .
our results will constrain future modeling of network activity in the hippocampus and para - hippocampal cortices . because different anatomical cell types have different projection patterns , burstiness , and theta
rhythmicity / skipping might be passed on differentially to hippocampal subfields like the dentate gyrus , which receives massive mec l2 stellate input ( varga et  al . , 2010 ) , and ca1 , which receives some mec l2 pyramidal input ( kitamura et  al . , 2014 ) .
some grid cell models suggest that grid cells are generated by network mechanisms where a large number of similar ( stellate ) cells self - organize to generate symmetrical firing patterns either via continuous attractors or via oscillatory interference ( for reviews , see giocomo et  al . , 2011 ,
others have suggested mechanisms based on anatomical microcircuits ( brecht et  al . , 2013 ) .
our results do not resolve this question , but we add to the picture that the  network mechanism distributes firing patterns differentially according to cell type .
we conclude that the anatomical identity of the neuron is a strong determinant of the firing pattern .
analysis of burstiness , theta cycle skipping , and phase precession jointly suggest similarities between layer 3 cells and layer 2 stellate cells on one hand and layer 2 pyramidal cells and parasubicular cells on the other hand .
all experimental procedures were performed according to the german guidelines on animal welfare under the supervision of local ethics committees . in this paper , we analyzed a dataset of juxtacellular recordings from the superficial medial entorhinal cortex and the parasubiculum that we have published previously ( ray et  al .
2014a , tang et  al . , 2015 , tang et  al . ,
detailed descriptions of recording procedures ( pinault , 1996 , lee et  al . , 2006 , herfst et  al . , 2012 , tang et  al . , 2014b ) , quality control ( joshua et  al . , 2007 ) , tissue preparation , immunohistochemistry , and image acquisition ( naumann et  al . , 2016 , ray and brecht , 2016 )
in addition to labeled cells , we included a number of unlabeled , regularly spiking cells from mec l2 in our analysis .
these cells were assigned as either putatively calbindin - positive ( pcb+ ) pyramidal cells or putatively calbindin - negative ( pcb ) stellate cells based on their theta strength and preferred theta phase angle using the classification approach of tang et  al .
( 2014a ) ; i.e. , based on the theta strength and preferred theta phase angle of spiking activity . as in tang et  al .
( 2014a ) , we used a 0.1 guard zone and found that the cells were well separated with no cells in the guard zone ( figure  s1 ) . in the manuscript
, we refer to the pooled groups of identified and putative calbindin - positive pyramidal cells simply as 
 when we show example cells of the  four cell types ( figures 1a and 1b and 6a6c ) , we show only identified cb+/ cells . in figures
s3 and s5 , we show analysis of a dataset where we only included identified cells . to determine whether a neuron was discharging in a bursty pattern , we analyzed the isi histogram using a similar approach as latuske et  al .
isis below 60  ms were binned in 2-ms bins and normalized to area  = 1 to generate a probability distribution ( figure  2a ) . a principle component analysis ( pca )
was done on a matrix of the isi probability distribution of all neurons (  pca  in matlab , mathworks ) . for plotting ,
the density of cells in this space was estimated with a 2d gaussian kernel density estimator (  kde2d  ; botev et  al . , 2010 ) .
the neurons were assigned to two clusters using a k - means clustering algorithm on the first three principal components (  kmeans  ; matlab ; figure  2c , top ) . to assess the separation quality of the two clusters ,
we calculated the projection of the neurons onto fisher s linear discriminant ( the burstiness , using  lda  from scikit - learn in python ) and found that the two clusters ( non - bursty and bursty ) were well separated with little overlap ( figure  2c , top ) . to check whether the distribution of burstiness was bimodal , thus reflecting two distinct classes of isi histograms , we fitted probability density functions for gaussian mixture models with between one and three underlying gaussians and compared the models using the akaike information criterion ( akaike 1974 ; aic from  gmdistribution.fit  in matlab ) . a bimodal distribution best explained the data ( aicunimodal  = 622.7 , aicbimodal  = 609.6 , aictrimodal  = 614.7 ) . based on the mean and variance of the two gaussian distributions underlying the observed distribution of burstiness ( figure  2c , bottom , dashed red lines ) , we estimated that excluding cells where 0.4  < burstiness  < 1.5 would yield > 95% correct labeling of non - bursty and bursty neurons in the non - bursty and bursty categories and used this as a guard zone ( latuske et  al . , 2015 ) . during recording , the juxtacellular traces were digitized at 20 khz . to analyze the spike shapes , we first zero - phased high pass - filtered the raw signal at 100  hz with a finite impulse response filter of order 2 (  fir1  in matlab ) .
the spike times were detected by thresholding the filtered signal and saving each threshold crossing  2.5  ms .
spike sorting based on the first principal components was performed on these 5-ms snippets to remove any threshold crossings because of artifacts in the signal ( tang et  al . , 2014a ) . to align the spike shapes optimally after spike sorting , the 5-ms snippets were oversampled at five times their original sampling rate using a spline interpolation (  interp1  in matlab ) and
were then aligned to the peak sample . to ensure that we were only analyzing shapes free of distortions because of drift of the pipette and that the spikes were well above the noise floor , we only analyzed spikes for which the spike amplitude was in the top 60 - 90 percentile and where the z score of the spike amplitude was > 17 .
the noise floor was defined as the mean of the first and last 0.5  ms of each 5-ms spike snippet .
we also  removed any spikes where there was another spike in the preceding 100  ms . in the four cell groups ,
there were only a few cells where the spikes did not have sufficient quality to analyze the spike shape , and we could analyze 19/22 parasubicular cells , 24/31 mec l2 pyramidal cells , 58/68 mec l2 stellate cells , and 27/32 mec l3 cells .
we calculated the mean spike shape of every cell and determined the spike features from these traces . for plotting the comparison between cells and for illustrating the differences in peak - to - trough time ( figures 4a and 4b )
, we normalized the spike shape by subtracting the noise floor , dividing the mean spike by the peak - to - trough height , and setting the peak height to 1 . to determine whether a neuron was rhythmic and theta cycle - skipping , we used an mle of a parametric model of the isi histogram (  mle_rhythmicity  ; climer et  al . , 2015 )
for every cell , we fitted three models to the isi distribution : a flat model with no rhythmic components , a rhythmic , non - skipping model with a rhythmic modulation of the isi histogram , and a rhythmic , cycle - skipping model with a rhythmic modulation of the isi histogram and a second parameter introducing theta cycle skipping ( i.e. , a higher amplitude of every other peak in the isi histogram ) .
when fitting the models , we searched for a rhythmic component with a theta frequency between 5 and 13  hz and for cycle skippings > 0.01 .
the three fitted models were compared using the appropriate  statistic ( calculated from the maximum log likelihood of the models ) to generate two p  values : prhythmic ( comparing the flat and the rhythmic , non - skipping models ) and pskipping ( comparing the rhythmic , non - skipping and the rhythmic , cycle - skipping models ) .
the cells were classified using a two - level classification ( figure  5b ) .
first we determined whether a cell was rhythmic ( prhythmic  < 0.05 ) or  non - rhythmic ( prhythmic > 0.05 ) .
then we classified the rhythmic cells as either  rhythmic , cycle - skipping ( pskipping  < 0.05 ) or rhythmic , non - skipping ( pskipping > 0.05 ) . to statistically assess whether theta cycle skipping cells were rarer among rhythmic cells in the generally non - rhythmic cell types ( mec l2 stellates and mec l3 neurons ) than in the generally rhythmic cell types ( parasubicular neurons and mec l2 pyramids ) , we fitted a mixed - effects logistic regression .
we constructed a vector , isgenrhytm ( which takes the value 1 for parasubicular neurons and mec l2 pyramids and the value 0 for mec l2 stellates and mec l3 neurons ) .
we also constructed a vector type that simply dummy - coded the four neuron types from 1 , 2 , 3 , and 4 .
we dummy - coded when the neuron was theta cycle - skipping in the vector isskipping .
we then modeled the probability of being rhythmic as a function of being generally rhythmic while controlling for the different number of cells in the four categories of neurons : 
isskippingisgenrhytm  + ( 1|type )  using  fitglme  in matlab ( aarts et  al . ,
in addition to the mle approach , we also calculated the theta strength and preferred theta phase of every cell .
the local field potential was bandpass - filtered in the theta range ( 412  hz ) , and a hilbert transform was used to determine the instantaneous phase of the theta wave for every spike .
the theta locking strength and the preferred phase angle were calculated as the modulus and argument of the rayleigh average vector of the theta phase at all spike times .
statistical modeling ( generalized linear models ) was done in matlab using the  glmefit  function .
we modeled burstiness as a function of theta strength as a normally distributed variable ( figures s3a  s3c ) .
we modeled the peak - to - trough time as a gamma - distributed variable with a reciprocal link function in matlab because it can only assume positive values ( figures s4a  s4c ) . to compare models , we either calculated and compared the aic ( akaike , 1974 ) or , in the case of nested models , calculated the p value from likelihood ratio tests . in the manuscript , we describe all statistical models using standard wilkinson notation ( wilkinson and rogers , 1973 ) . to identify coherent periods of elevated firing (  single runs  ) , we follow a previously applied strategy based on the temporal structure of the recorded spike trains ( reifenstein et  al . , 2016 ) .
briefly , we convolved the spike train with a gaussian kernel to estimate the instantaneous firing rate .
we then used a firing rate threshold to locate periods of elevated firing ( figure  6a , top ) . for each of the single runs ,
we quantified phase precession by calculating the slope , phase range , and circular - linear correlation coefficient of the circular - linear regression line ( figure  6a , bottom ; kempter et  al . , 2012 ,
c.l.e . analyzed burstiness , spike shape , theta rhythmicity , and cycle skipping . e.t.r .
, r.k . , and m.b . provided expertise and feedback on the analysis and supervised the project .
 summarythe medial entorhinal cortex ( mec ) and the adjacent parasubiculum are known for their elaborate spatial discharges ( grid cells , border cells , etc . ) and the precessing of spikes relative to the local field potential . 
 we know little , however , about how spatio - temporal firing patterns map onto cell types . 
 we find that cell  type is a major determinant of spatio - temporal discharge properties . 
 parasubicular neurons and mec layer 2 ( l2 ) pyramids have shorter spikes , discharge spikes in bursts , and are theta - modulated ( rhythmic , locking , skipping ) , but spikes phase - precess only weakly . 
 mec l2 stellates and layer 3 ( l3 ) neurons have longer spikes , do not discharge in bursts , and are weakly theta - modulated ( non - rhythmic , weakly locking , rarely skipping ) , but spikes steeply phase - precess . 
 the similarities between mec l3 neurons and mec l2 stellates on one hand and parasubicular neurons and mec l2 pyramids on the other hand suggest two distinct streams of temporal coding in the parahippocampal cortex . 
Introduction
Results
Discussion
Experimental Procedures
Author Contributions