BACKGROUND
spines are protrusions that occur on the dendrites of most mammalian neurons.
they contain the post-synaptic apparatus and have a role in learning and memory storage.
spine distribution is a critically important question for multiple reasons.
changes in spine distributions and shape have been linked to neurological disorders such as fragile x syndrome
<cit> .
spine distributions determine the extent to which the neuropil will be electrically sampled, i.e.
dense distributions will sample the neural connectivity map more fully
<cit> .
furthermore, the nature of optimal sampling is unknown and likely depends on the surrounding anatomy and the total information content available to dendrites.
because pruning takes place during development in an activity dependent manner, spine distributions may reflect activity within neural circuits.
distributions of spine types are biologically important because the electrical properties of spines, such as the spine neck resistance, promote nonlinear dendritic processing and associated forms of plasticity and storage
<cit>  to enhance the computational capabilities of neurons.
the shapes and types of dendritic spines contribute to synaptic plasticity.
because neighboring spines on the same short segment of dendrite can express a full range of structural dimensions, individual spines might act as separate computational units
<cit> .
nevertheless, the dendrite acts in a coordinated manner and thus the spatio-temporal distributions of different spine types is likely to be significant.
little is known about this population level organization of dendritic spines.
our aim was to perform an exploratory analysis of neuronal data from different time periods during the growth of rat dissociated hippocampal neurons, a well-established model system
<cit> .
the observations here pertain only to the culture system and not necessarily to in vivo settings although the analytical tools used here could be adapted to in vivo analyses.
by quantifying populations of dendritic spines with automated tools at a global level, we were able to perform a much larger and more comprehensive analysis than most previous studies.
many studies only analyze a small region of interest on the largest dendrites, for example the 50– <dig> μm closest to the soma
<cit> , or  <dig> μm segments
<cit> , making it easier to measure manually the spine type counts and dimensions.
other works determine spine lengths and widths by manually drawing a line along the maximal length and measuring the length of that line
<cit> , and therefore are only able to analyze a few neurons and a few hundred spines at a time.
in this study we determined the ratios of spine types along the dendrites as a function of time in culture, clustering or repulsion of spines in space, and how best to model spine type distributions.
a model that fits the spatial distribution of spine types in healthy cultured neurons would be useful to assess neuronal cultures following gene manipulation such as rnai and to study features of induced pluripotent stem cells differentiated to neurons.
log-linear models  and multinomial logistic regressions  are two basic and essential statistical methods, and have an extensive history of being used in biological studies.
however, these tools have not been used thus far in the analysis of spine distributions.
we fit a log-linear model  to the contingency table of spine features to determine the dependence between spine types , distance from the cell body along the dendrite , the branch order of the dendritic segment on which it lies , and the day in vitro  on which it was imaged.
once we determined which of these attributes contributed to the overall dendritic spine model, we then asked whether these attributes can predict the occurrence of spines and of spine types.
to answer this question we used a multinomial logistic regression  model, which predicted the spine type, using the attributes that were found to be important through the llm and associated contingency tables.
finally, to understand how the dendritic spine density varied over the length of the neuron or whether the appearance of spines was completely spatially random i.e uniformly distributed over the neurites, we made use of spatial point processes.
spatial point processes have been used before in biological studies to model the locations of entire neurons
<cit> , locations of ants nests
<cit>  or xylem conduits
<cit> .
there have also been other more ad-hoc methods created to study the number of "clustered spines" on each dendritic segment in monkey brains, where a cluster is defined as a group of  <dig> or more spines
<cit>  however we believe our use of the linear network k-function
<cit>  is the first work to analyze the locations of dendritic spines and their clustering properties in such a principled manner.
our analysis indicated that the density of spines is generally completely spatially random  over the dendritic length probably due to the absence of instructive directional signals found an in vivo setting in which spine distributions are unlikely to be csr.
methods
cell imaging
dissociated hippocampal neurons from embryonic rat brains  were plated onto poly-l-lysine coated coverslips.
once neurons adhered to the coverslip, they were placed face-down on glial cells grown in vitro for  <dig> days.
these neurons were a primary neuronal culture system, and no cell line was used.
neurons were grown for specific time periods up to  <dig> days in a neuronal medium containing b <dig>  this co-culture of neurons and glia mimic the physiological conditions of neuronal growth and development in mammalian brain
<cit> .
work with the neuronal cultures was approved by the ucsb animal care committee.
to fill the neuronal processes including dendritic spines green fluorescent protein  was expressed from a plasmid containing the beta-actin promoter
<cit> .
of this plasmid,  <dig> μg was transfected into each coverslip containing about  <dig>  neurons .
transfection was performed as described in the manufacturer’s protocol  with minor changes.
the transfection mix and neurons were incubated for two hours to avoid toxicity caused by lipo <dig>  following transfection, coverslips were flipped back onto the glial dish, where they were originally cultured.
gfp-actin transfected into the neurons at div <dig>  and neurons were studied at three time points- div <dig>   <dig> and  <dig>  these time points survey the maturation period over which synapses and spines emerge
<cit> .
note that these were not the same neurons studied over time, but each time point represents a different population of neurons which were grown in culture up until the point of imaging.
in this way our analysis represents a study at the population level.
at each time point the number of images taken per plate depended on the transfection efficiency of that plate.
on average approximately 1% of cells were transfected.
the plating density was set so that neurons were relatively isolated in order to capture one neuron per image.
an olympus fluoview laser scanning confocal microscope was used.
image slices were  <dig> by  <dig> pixels at 154nm per pixel resolution.
there were 7– <dig> z-slices per stack depending on the depth of the neuron, taken at 200nm steps.
this means that the stacks were  <dig>  μm ×  <dig>  μm ×  <dig> – <dig>  μm.
the z dimension slices were used to capture each depth level at the optimal focus, however we cannot claim to have accurate volumetric information at this resolution.
a 40x oil objective lens with no optical zoom was used.
numerical aperture  was  <dig> , and illumination conditions were kept constant.
deconvolution of the raw data before processing was not necessary because the images were clear enough to manually annotate the neuron traces and manually edit all the spine detections and types as described in the following section.
we performed three biological replicates, the results of which are detailed below.
although there are other higher resolution, full volume methods, the analysis of this data is broadly applicable to imaged neurons in other systems
<cit> .
we attempted to capture the entire neuron in each image, however because of limits in available imaging techniques we found that this does not always happen.
in the cases where dendrites were truncated at the end of the image plane we assumed that the proportion of spines in the missing data was similar to what had already been observed, and therefore the resulting distributions did not change.
we verified this assumption visually by taking tiled mosaics of a few neurons imaged in their entirety from each div and checking that the branch orders, distances to soma and spine type counts were unchanged as compared to those of the same div.
there was an observed increase in the dendritic length truncated by the image plane as the div increased.
however in our particular analyses the methods used, such as the log-linear model and multinomial logistic regression, were focused on trends between spine characteristics such as distance to soma and type and these trends are innately unaffected by the truncation of dendrites given the above assumption.
in addition, spatial point process analyses such as the linear network k-function always include the specification of an observation window
<cit> , which in our case was the image plane.
we verified  that the overall spine density and the density of each spine type did not vary with distance from the soma so that we could assume spine density at the ends of the dendrites which were truncated was similar to the dendritic length which was observed.
we recognize that we cannot see the proximity of labeled cells to other neurons which haven’t taken up the gfp labeling.
these unlabeled neighboring neurons may cause some difference in spine distributions which we cannot quantify.
for this reason we have attempted to quantify our biological findings statistically over entire experiments and div time points instead of by individual neurons, although in certain cases showing results from individual randomly sampled neurons was necessary.
neuronal reconstruction
there exist many automated methods for studying neuronal growth and morphometry and therefore we present a brief review of available software for tracing dendrites and detecting and classifying spines.
in particular, neuronj
<cit>  is a widely used software; however it is only semi- automatic and one must click several points to trace each neurite.
the labeling is done manually and the statistics output only include lengths of neurites and not spine data.
hca-vision
<cit>  is a costly software with similar goals, however the parameters of the neurite tracing are set manually with a sliding bar and thus results require much hand-tuning.
in addition, it is also focused on tracing neurites as opposed to spine analysis.
for a full review of existing methods and softwares for neuron tracing and spine detection see
<cit> .
we found neuronstudio
<cit>  to be the most user-friendly, and for this reason we used it to annotate dendrites and spines for this analysis.
despite the abundance of automated softwares, neuronal reconstructions are still largely performed by hand
<cit>  and this is is especially essential for a study like this one, where the traversed distance of the dendrites and number of spines and their shapes were analyzed in such detail.
using automated reconstruction algorithms on raw data is prone to both false positive and false negative detections of spines, as well as misleading spine shape measurements.
in cases where neurites from neighboring neurons enter into an image , neuronstudio often incorrectly traces these neurites as belonging to the neuron of interest.
for this reason we manually traced each dendritic branch and soma of each neuron, ran neuronstudio’s automated spine detection/classification algorithm and then manually inspected and verified each spine’s location and type.
the verification and tracing were done by the primary author and an undergraduate biology student working in the kosik lab .
they were both familiar with dendrite and spine morphology and the resulting annotations from each were cross-checked by the other.
relevant spine attributes output from the neuronstudio software include branch order , type , distance to soma along dendrite , length  and width at widest point .
however since neuronstudio uses the length and width of the spines to determine the spine type, we chose to make use of spine type and discard the other  <dig> measurements.
neuronstudio uses centrifugal labeling for branch orders, meaning it starts at  <dig> at the cell body and moves outwards, incrementing at every y-shaped bifurcation regardless of the diameter of the daughter branches.
note that the entire image stack with z-dimension information was loaded into neuronstudio for the spine classification, and that the software has interpolation algorithms to estimate the spine type in 3d.
for spine detection the default cut-offs were used, i.e.
a required spine height between  <dig> – <dig> μm, a maximum spine width of  <dig> μm, a minimum stubby size of  <dig> voxels , a minimum non-stubby size of  <dig> voxels, and automatic z-smear compensation.
for spine classification, the default settings were also used, i.e.
a head-to-neck ratio threshold of  <dig>  μm, an aspect ratio  threshold of  <dig>  μm and a minimum mushroom head size of  <dig>  μm.
neuronstudio delineates spine types by these  <dig> thresholds.
it is generally known that mushroom spines have a large head and a narrow neck, thin spines have a small head and a narrow neck, and stubby spines display no obvious subdivision in head and neck.
if the head-to-neck ratio is above the threshold and the minimum mushroom head size is met, the spine is considered mushroom.
if both the head-to-neck and aspect ratios are lower than the respective thresholds then the spine is considered stubby.
the remaining cases result in thin spines.
for further information on neuronstudio reconstruction, detection, and spine classification algorithms please refer to
<cit> .
in addition to the spine information, a trace file is output which labels the cell body, branch points and end points of the dendrites.
the trace provides a skeletonization, or centerline, of the dendrite which we used to compute the linear network distances in the following analyses.
log-linear model as a tool for exploring important features and their dependencies
to find the most influential attributes with regard to prediction and spatio-temporal modeling of spines we fit a log-linear model to the feature data, which is a type of generalized linear model
<cit> .
the co-occurrence frequencies of the features in question are essentially a large multidimensional contingency table of counts.
the standard linear models assume that data is normally distributed around a certain mean, which means that the observations can take any real value, positive, negative, integer or fractional.
log-linear models, on the other hand, assume that data is intrinsically non-negative, typically counts that could be poisson distributed, and allow us to model the association and interaction patterns among categorical variables.
the attributes under consideration are bo, type, sd and div.
again, since the type of spine was quite directly dependent on the length and the head diameter of the spine, we left these latter variables out of the modeling.
in order to analyze the data using a log-linear model, the various features must be in a categorical form or discretized.
in an exploratory analysis such as this, one does not know what dependencies among features to expect; however we would like to note that these dependencies were not lost in the discretization process since trends in increasing and decreasing feature values would be preserved.
to ensure that there were a reasonable number of observations at the higher branch orders, we pooled bo values of  <dig> or higher into a single category called "higher-order branches".
we created a categorical variable to represent the continuous variable soma distance  where categories were determined using the  <dig> quartiles of the sd spine data pooled over all  <dig> experiments.
specifically, sd values of less than  <dig>  μm were classified into the first group, from this value to less than  <dig>  μm the second, from this to less than  <dig>  μm the third, and the rest  fell in the fourth group.
binning the observed data for the continuous variables is the best way to get a general feel for how these quantities relate to each other.
after this post-processing of the data we arrived at  <dig> categories of branch order,  <dig> categories of soma distance,  <dig> spine types , and  <dig> divs .
using the observed frequencies for the aforementioned attributes, we created a four-way contingency table and fit the model using the 'glm’ function in the r package 'stats’.
the table of the frequency of occurrences of the four attributes was modeled as poisson with each entry being a simple count of the co-occurrences of that bin.
we called this count f
ijkl
with each of the subscripts i,j,k,l corresponding to a different attribute.
the method uses the link function y
ijkl
= log, and treats the model as a regular linear model.
each entry y
ijkl
is modeled by a combination of coefficients: the intercept, plus main effects, plus every combination of interactions between these four attributes, as shown below.
yijkl=μ+αi+βj+γk+δl+ij+ik+il+…+error.
we estimated this full interaction model using the least-squares maximum-likelihood approach.
we also used a stepwise fit algorithm, which begins with a model that includes only the constant term, and at each step chooses whether or not to add one additional term.
the algorithm begins with the main effects then tries each possible 2-way interaction, aiming to minimize the akaike information criterion .
the aic is defined as
aic=2k-2ln)
where k is the number of parameters i.e the total number of coefficients being estimated, and
l=maxθ∏n=1neynθ′xne-eθ′xnyn!
is the maximized value of the likelihood function for the estimated poisson model.
in the above equations
x=x <dig> …,xn∈r <dig> are the input vectors,
θ
= θ1…θ
k
are the parameter values , and
y=y <dig> …,yn∈r is the output.
the aic is a commonly used goodness-of-fit measure for a model given the observed data.
adding or subtracting terms, whether they be main effects, pairwise interactions, or up to 3-way interactions between attributes, will change the aic value for the model.
a lower aic criterion indicates a better fit to the data and therefore a better model.
to compute the stepwise fit we used the r function 'step’.
for more information on the stepwise fit algorithm as well as the aic criterion we ask that the readers refer to the 'step’ function reference .
we ran both of these llm fitting procedures for all  <dig> experiments separately expecting to find general agreement between coefficients of the corresponding models created.
multinomial logistic regression to predict spine type from neighbor types
in order to predict spine type we first determined which attributes contributed most to spine type prediction.
given the complexity of the multidimensional llm and the various interactions and conditional frequencies that would impinge on this issue, we decided to determine these attributes by analyzing 2-way contingency tables for spine type vs. sd, bo, div, as well as the spine types of the  <dig> nearest neighbors.
this analysis helped us pick attributes that would be useful as the predictors in the multinomial logistic regression
<cit>  explained below.
when the response variable of a regression takes binary values "logistic regression" is used.
this is an approach which uses a linear combination of the predictor variables to predict the log-odds of a success .
since our response variable was spine type and it can take  <dig> values , we needed to use a "multinomial logistic regression"  which attempts to model the probability of any of multiple possible outcomes.
we did not use the attributes sd or bo as predictors variables since the results of both the llm analysis and 2-way contingency tables mentioned above told us that these quantities were not as relevant for spine type prediction.
therefore our model consisted of spine type as the output variable and the div, 1st, 2nd and 3rd nearest neighbor type along the dendrite as the predictor variables.
we tried using only  <dig> or  <dig> nearest neighbors, however the results proved inconclusive because the prediction probabilities for each of the  <dig> types were predominantly close to 1/ <dig>  if we used more than the  <dig> nearest neighbors we sometimes ended up spanning a segment of dendrite which we did not consider to be "local", so we decided that  <dig> nearest neighbors provide the most useful information in the case of this study.
the mlr analysis we performed in this paper does disregard the actual inter-spine distances, meaning that if the  <dig> nearest neighbors are very close or very far apart we still treat them the same.
we did this partially because adding the distance variables would complicate the model significantly, but also because we believe that over a large population of spines such as the one we have, these differences in distance will average out and we will still get a general picture of the trends between neighboring spine types.
to verify that this was true we computed a histogram showing the distribution of 3rd nearest neighbor distances for each spine, shown in figure
<dig>  although the maximum distance to any 3rd nearest neighbor is extremely high  we can see from the histogram as well as the fact that the median 3rd nearest neighbor distance was  <dig>  μm that this distance is clearly an outlier case and that the majority of 3rd nearest neighbor distances lie below  <dig> μm.
suppose the output variable categories are denoted by  <dig> , <dig> corresponding to mushroom, stubby or thin spines, with  <dig> being the reference category.
if y
i
denotes the observed outcome of the output variable , and x
i
is the corresponding vector of the  <dig> neighbor types and div for the ith observation, one regression is run for the logit probability of each category k, with β
k
representing the vector of regression coefficients in the kth regression .
this is done for all but the reference category, whose probability is then obtained by subtracting all other probabilities from one .
note that because the predictor variables were spine types, which were nominal as opposed to ordinal variables, the predictor variables x
i
must be represented with a "dummy coding".
this means each neighbor type was represented by  <dig> predictor variables, where  corresponded to mushroom type,  corresponded to stubby type and  corresponded to thin type.
this does not need to be done for the output variable y. with the addition of the div, which does not have to be dummy coded since it is an ordinal variable, this made each x
i
vector of length  <dig>
the regressions are then written as:
p=exp1+exp+exp
p=exp1+exp+exp
and
p=1-p-p=11+exp+exp
the parameters are estimated typically by using an iterative procedure such as "iteratively re-weighted least squares"  or, more commonly by a numerical approach  such as the broyden-fletcher-goldfarb-shanno  method.
in our case we create an mlr using the command multinom in the r package nnet
<cit>  which uses bfgs by calling the r function optim.
it can be seen that
logp)=β1xi
logp)=β2xi
so that the beta coefficients represent the change in the log odds of the dependent variable being in a particular category with respect to the reference category i.e.
the thin type, for a unit change of the corresponding independent variable.
to check if the models created from all three experiments were in agreement, we ran the mlr separately for each experiment.
to satisfy one of the major assumptions of this analysis, namely that the data must be a set of independent observations, we took  <dig> randomly sampled spines of each type from each experiment  to use for the parameter estimation.
we chose to select equal proportions of each spine type in order to remove any bias in the model towards the less frequent thin spines, and  <dig> was the largest number we could justify using since there were only  <dig> thin spines in experiment  <dig>  we verified that these randomly sampled spines did not lie within  <dig> μm of the image border so that we were fairly certain their nearest neighbors did not fall outside of the image plane.
note that due to the tortuosity of the dendritic structure this did not mean that our sample was necessarily biased towards spines which were proximal to the soma.
we did not verify explicitly that the sampled spines were not neighbors of each other, since we assumed that the variation captured by the random sampling was enough to ensure some level of independence.
the idea was to aim for an independent set of observations which represented the entire "population" of spines in that experiment.
to be clear we used all  <dig>  spines for the llm model and k-function analysis, only the mlr model required random sampling since we were using neighbor information which would have been redundant if we considered every spine.
to verify that the prediction of spine type provided by the mlr was better than what we would get purely by their relative abundance i.e.
without neighboring spine type information, we computed something similar to a "bayes factor"
<cit> .
bayes factor is a method of choosing between two models on the basis of the observed data.
in our case, the first prediction model was simply the prior global probability of finding a given spine type based on its frequency in the particular experiment under consideration.
the second model was the mlr prediction model using the neighbor type information.
we computed p/p and reasoned that values considerably larger than one indicated the neighboring spine type information was helpful in the prediction of the central spine type.
linear network k-function as a tool for testing spatial point patterns
originally proposed by ripley in 1981
<cit> , the purpose of the k-function is to estimate whether or not there is clustering or repulsion present in a given spatial point process.
the common null hypothesis is that the points within the observation window are distributed as a homogeneous poisson process, which is also termed "completely spatially random" or csr.
this means that the density of points does not vary depending on the spatial parameters i.e.
x and y in the 2d euclidean case, or the location along the dendritic network in our case.
in order to determine if this is a valid null hypothesis for our data, we created q-q plots
<cit>  for individual dendrites which compared the quantiles of the sd values of observed spines to the theoretical quantiles for the csr case.
if the two distributions  being compared were similar, the points in the q-q plot would approximately lie on the line y = x. in order to create the theoretical quantiles it is necessary to know the values of sd at any location on the given network, not just at the spine locations.
once we have this we can partition the network into epsilon small segments and assign each segment a value  <dig> if it contains a spine and  <dig> otherwise based on the csr assumptions.
we did this using code provided to us by adrian baddeley and gopal nair at the commonwealth scientific and industrial research organization , australia.
the k-function computes the expected number of points within a distance t of an arbitrary point p, therefore the empirical value in 2d euclidean space for the csr case will be proportional to the circular area, λπt <dig>  the proportionality constant λ represents the density of points in the homogeneous poisson case, and can be estimated by finding the total number of points n divided by the total area of the observation window a. ripley’s k-function, which is a function of t, is a very useful tool because it describes the 2
nd
order characteristics of the point process at several scales t. if we ignore the edge effects due to the observation window, the observed
k^ can be written as:
k^=|a|n2∑i∑j≠ii
where i stands for the indicator function, and d
ij
stands for the euclidean distance between two points p
i
and p
j
.
in the above equation, we see that the expectation is normalized by 1/λ since
λ=n|a|, so we infer that theoretically k = πt <dig> implies spatial independence of points, or a csr point process.
therefore, if k is the theoretical csr value of the function and
k^ is the observed function, then
k^>k implies clustering between points and
k^<k implies repulsion.
it is possible to extend this function to multi-type point patterns  or to higher dimensional data .
since our particular point process consists of spines which lie along the "linear network" of the dendritic tree we were primarily concerned with inter-spine distances along the dendrite as opposed to in euclidean space.
therefore we used a version of the k-function developed recently for linear networks by okabe and yamada
<cit> .
this modified version of the k-function takes into account the structure of the linear network on which the point process resides and imitates the euclidean space k-function described above.
the linear network k-function is calculated as follows:
k^=ℓtn2∑i=1n∑j≠ii
where ℓ
t
is the length of the total network l
t
.
the theoretical csr for this case is described as follows:
k=1ℓt∫p∈ltℓpdp
where p is a point belonging to the set of all points p = {p <dig> …,p
n
}, and ℓ
p
is the length of the subset of the network l
p
where the distance between p and any other point is ≤ t. note that here the distance d
ij
stands for the linear network distance along the dendrite.
accounting for variability in the length ℓ
p
means the formula takes into account the edge effects due to the observation window  inherently, but at the cost of added complexity.
the computation of the theoretical linear network k-function requires us to find
lpi, the subset of l
t
where the network distance between a specific point p
i
and any other point is ≤ t, and
ℓpi, the length of that subset, for every point p
i
.
a visualization of the quantities d
ij
, l
t
, ℓ
t
,
lpi, and
ℓpi is shown in figure
<dig>
lpi, and
ℓpi which were used to compute the linear network k-function.
here dij is the linear network distance shown by the gray line between points pi and pj.
lt  is the entirety of the single dendritic network and ℓt is the length of lt. similarly,
lpi  is subset of the network where the distance between a point pi and any other point is ≤ t and
ℓpi is the length of
lpi.
in this particular example there are  <dig> spines which fall within
lpi and would be counted in determining the empirical function value
k^, however point pj falls outside this radius and would therefore not be counted.
note that although many biological applications of point processes treat individual observations as replicate patterns coming from the same underlying distribution, we cannot do that using the above definition of the network linear k-function due to the change in linear network structure from dendrite to dendrite.
the term "dendrite" here refers to the entire dendritic tree resulting from a single root branch of a neuron.
other in-vivo studies
<cit>  focus on clustering of spines which lie on the same unbranched section of the dendrite, however we focus on the entire dendritic tree under the hypothesis that it follows rule-based distributions of spines due to anatomical constraints and integration of the a signal over the entire dendrite.
one can infer from figure
<dig> that since the geometry of the linear network changes from dendrite to dendrite, so do the total lengths of the networks ℓ
t
, the ranges of possible t-values and the amount of dendritic length that is present within a given distance of any point.
we did not simply normalize the lengths of the networks to a  scale because it is desirable for the t-axis to retain its real physical values in order to make conclusions about the scale  of clustering or repulsion among spines.
however, we did desire to compare the linear network k-functions of various dendrites in a meaningful way.
for this reason we used a corrected version of the network k-function that intrinsically compensates for the geometry of the network called ang’s correction
<cit> .
the observed k-function then becomes:
k^=ℓtn∑i=1n∑j≠iim
where m is the number of points of l lying at the exact distance t away from the point i measured by the shortest path.
that is, the contribution to the function from each pair of points  is weighted by the reciprocal of the number of points that are situated at the same distance from i as j is.
as a result, the theoretical csr case is simply k = t for all 0 ≤ t < t. this enables direct comparison of t-values across dendrites, as we will see in the results section.
simulations and q-values
to test the null hypothesis that the locations of spines on the dendrites were indeed csr, we created a summary statistic which encompasses the difference between the empirical
k^ and the theoretical k under csr.
the summary or "test statistic" we used, is the max absolute difference  over t, viz.
d=maxt|k-k^|.
one method for obtaining a distribution of d proposed by diggle
<cit>  is to bootstrap the residuals, or differences between the observed and theoretical values.
however a more heuristic and intuitive way is to simulate the csr case for each dendrite, compute the k-function for each of these simulations, and find the simulated distribution of our test statistic.
we then found the p-value of the observed difference d from this simulated distribution.
specifically, we carried out  <dig> csr simulations for each dendrite by placing uniform points on a line [ <dig> ℓ
t
], and mapping them to that specific dendrite’s linear network structure.
the number of points simulated per dendrite equaled the number of observed spines for that dendrite, thus preserving the overall density λ. this means the same number of spines that existed on each dendrite were randomly placed along the linear network specific to that dendrite.
we used these simulations to obtain  <dig> values of the summary statistic, say d. then the p-value for each dendrite was simply the proportion of simulated values that fell above the observed or experimental value of d, i.e.
the rank of this d within the  <dig> values of d, or nrank/.
this p-value approach is similar to the test which rejects the null hypothesis if the graph of the observed k-function lies outside the "point-wise simulation envelope" at any value of t. a simulation envelope is essentially a graphical measure of how far a function can deviate from the theoretical value without being considered significant at a given level.
as mentioned above in our case the envelope is calculated by first creating the  <dig> csr simulations of a point pattern on a given dendritic network with the same observed network intensity, then calculating the linear k-function for each of these  <dig> simulations.
to perform a two-sided significance test at the 10% level, the 5% and 95% percentiles are then calculated based off the  <dig> lowest and  <dig> highest linear k-function values per t-value, hence the term "point-wise".
plotting these values as a function of t gives one a visual idea of the spread that is produced by chance mechanisms alone.
if the observed k-function for a given t-value does not fall outside these percentiles, it is considered insignificant for that t-value at the 10% significance level.
we make use of the r package 'spatstat’
<cit>  for obtaining the point-wise simulation envelope.
because we have a multitude of hypothesis tests and p-values , to reach a conclusion about the general trend for each div and experiment, we used the concept of false discovery rate
<cit> .
the fdr is defined as
π0=#truenulltests#totaltests
controlling the overall fdr, or expected proportion of incorrectly rejected null hypotheses termed "false discoveries", is a statistical method commonly used in multiple hypothesis testing which increases the statistical power of each test.
what is more general and useful however, is a test-specific fdr measure.
this essentially allows us to look at all possible significance thresholds at once, as well as provide each test with a measure of significance that can be easily interpreted.
this is accomplished by calculating an analogue of the p-value for each test called a "q-value"
<cit> .
a p-value of  <dig>  implies that 5% of all tests will result in false positives, whereas a q-value of  <dig>  implies that 5% of significant tests will result in false positives.
since the latter is clearly a far smaller quantity, q-values generally indicate fewer significant tests than p-values for a given significance threshold and provide a far more accurate indication of the level of false positives in the case of multiple hypothesis testing.
for q-value estimation we used the qvalue package available from
<cit> .
RESULTS
data analyzed
we performed three biological replicate experiments resulting in a total of  <dig> neurons from the following time points: div  <dig>  div  <dig>  and div  <dig> .
this provided a rich and complete data set resulting in  <dig> dendritic branches and  <dig>  spines.
example images from each div along with zoomed in dendritic segments where spines and annotations are visible are shown in figure
<dig>  scale bars are shown in red in panels a-c and the yellow rectangular boxes in panels a-c show the region of interest which has been zoomed in on in panels d-f respectively.
panels d-f are all at the same resolution.
the number of spines per μm, or λ, for each dendrite in different experiments and time points is shown in figure
<dig>  we chose to include this in order to help the reader compare these neuronal culture results with other experimental paradigms with which they may be more familiar.
it is clear from the histograms that the distribution of spine density for div <dig> is skewed toward lower values as compared to the density for div <dig>  as expected.
the image data as well as spine and trace annotations are made publicly available through the bisque system
<cit>  and the url is given in the section titled "availability of supporting data".
we chose bisque over other databases like neuromorpho.org
<cit>  because it allows us to upload multiple layers of annotations as opposed to only the digital reconstruction files.
we calculated a 2-way contingency table over all experiments and spine types and obtained table
<dig>  from this table we note the high frequency of mushroom and stubby spines as compared to thin spines, and also the fact that the ratio of types does not remain the same per experiment even though they were indeed biological replicates.
in fact, a pearson’s chi-squared test on table
<dig> shows dependence between the spine type counts and experiment number, χ2 =  <dig> , p <  <dig> .
we believe that the large experimental variation between spine type proportions and counts in each experiment was a positive result because this meant that statistical agreement across all  <dig> experiments relating to spine type clustering and density estimation carries heavier weight than if the  <dig> experiments were more uniform in these quantities, or if we had pooled data from all  <dig> experiments together.
also, if all  <dig> experiments were unusually homogeneous there could be a possibility that it is a result of our specific culturing, imaging or spine extraction methods rather than a true representation of the underlying biological process.
the various biological systems to which these techniques will be applied will certainly have this type of variability.
spine type is independent of distance from soma
as described in the methods section, we calculated a stepwise-fit of the log-linear model starting with just a constant term, and at each step choosing to add the main effects  and possible 2-way interactions between main effects one-by-one if they decreased the corresponding aic value.
the captions above tables
<dig>
<dig> and
<dig> show the final models arrived at for each of the  <dig> experiments as well as their corresponding aic values.
the tables indicate the change in the aic value that would occur from adding or omitting each of the terms in the first column.
this gives us an idea of how important that term was to the model.
the rows of the table are ordered by their overall contribution to the model, i.e.
the term in the first column of the first row of each table had the lowest aic value and was therefore the most important to the overall model.
if the reader requires further information on the aic criterion or how to interpret this table we ask them to refer to chapter  <dig> of
<cit> .
despite the fact that they were included in the final stepwise fit model for experiments  <dig> and  <dig>  the aic values in tables
<dig>
<dig> and
<dig> show that in all  <dig> experiments the interaction between spine type and soma distance  as well as spine type and branch order  were the least important in modeling the overall frequency table of occurrences.
this implies that the correlation between these quantities was not very high, therefore we reason that it was not necessary to use either sd or bo to predict the spine type in the mlr created in the following section.
we also noticed that the term marking the interaction between bo and sd was the most important pairwise term in all stepwise fit models.
it is expected that bo and sd are correlated because both necessarily increase as we move away from the cell body.
indeed, running a 2-way chi-square test on the contingency table of the discretized versions of these variables showed us high dependence, χ2 =  <dig> ,p< <dig> .
we also saw a high level of dependence between div and sd  =  <dig> , p <  <dig> ) and between div and bo  =  <dig> , p <  <dig> ).
this was intuitive as well since we expect both bo and sd to generally increase with div.
it is possible that the type vs. sd relationship could have also been estimated using a sholl-type analysis  where we count the number of each type occurring within concentric circles from the soma and verify that it is constant, however this would not necessarily produce the same results.
the crucial difference between our approach and the sholl approach is that in our approach the "distance from soma measures" the actual distance along the centerline of the dendrite instead of the radial distance from the cell center.
this is especially important for dendrites with high tortuosity , since the radial distance in those cases will not correspond to the dendritic distance from the cell body.
many studies of cultured neurons use sholl analysis, however they use it in its original form for counting dendritic intersections and do not comment on the relation to spine density or type.
to our knowledge this is the first study to quantify the spine density vs. distance to the soma in dissociated neuronal cultures.
three-way and 4-way interactions are generally known to be weak  and difficult to interpret, however in the interest of exploring all possibilities we computed the maximum likelihood fit using all  <dig> attributes as well as a stepwise fit model which allows for 3-way interactions between attributes.
the table presented in additional file
<dig> results from the llm which models all possible interactions of all  <dig> attributes, i.e.
up to the fourth order.
the coefficients presented in the table are those which were significant at the  <dig> % level, and the corresponding p-values are shown in the last column.
the table contains the interactions which were more important to the model, and shows that of these interactions only one  between type and either bo or sd, was shown as being significant over all experiments.
this verifies once again that neither bo nor sd were highly correlated with the spine type.
in addition to this, the stepwise fit models in additional file
<dig> show that if we did allow 3rd order interactions, the strongest 3rd order correlation over all experiments was that of div, sd and bo, again affirming that all  <dig> of these quantities should intuitively increase together.
spines tend to cluster with other spines of the same type
in creating a regression model, we first ascertain that the predictor variables used are not only useful in predicting the output variable, but also that they do not provide redundant information as this can throw off the model fitting process.
using all spines in the dataset, we performed a chi-square test on the 2-way contingency tables of spine type versus binned sd and bo, div, and the types of the  <dig> nearest neighbors  as described in the log-linear model section above.
due to the aforementioned dependence between the type and experiment number we performed the test separately for each experiment and the results are shown in table
<dig>  from the table we can see that the div and the  <dig> nearest neighbors showed clear dependency with spine type in all experiments, whereas sd and bo showed independence at the 5% significance level in experiments  <dig> and  <dig> respectively.
since we expected sd and bo to have a similar relationship with type due to the high correlation mentioned above, and we had found this was not a very strong relationship, we chose to use only div, n <dig>  n <dig> and n <dig> as predictors for spine type in the mlr model.
n
n
n
the resulting beta coefficients for each of the predictor variables are shown in table
<dig>  here "n1-var1" refers to the beta coefficent of the first dummy variable for the type of the first nearest neighbor; "n1-var2" refers to the second dummy variable, and so on.
the "mushroom" row is omitted because it is the reference category and its probability is obtained as shown in eqn.
<dig>  we computed the prediction probabilities for each spine type given each combination of neighbor types for each experiment separately to determine the agreement between experiments.
a selected set of results are shown below in tables
<dig>
<dig> and
<dig>  the highest probability for each row is marked by an asterisk.
note that in these tables all divs in all experiments predicted the spine type to be mushroom when its  <dig> nearest neighbors were mushroom type, and stubby when the  <dig> nearest neighbors were stubby type.
thin types were the most probable when the three nearest neighbors were thin type in all but experiment  <dig> div <dig> and div <dig>  the probabilities for cases where all  <dig> of the nearest neighbors were not of the same type have been omitted for brevity and because they did not show any clear trends.
* denotes the highest probability per row.
*denotes the highest probability per row.
* denotes the highest probability per row.
the bayes factor results in table
<dig> show that the proportional gain in information for the spine type in question was always greater than one for the prediction of a particular type when the neighborhood types were all of that same type.
due to the low frequency of thin spines, their corresponding bayes factors were higher than that of other types, meaning that their prediction probabilities benefit more than other types from neighborhood type information.
dendritic spine densities are completely spatially random
we created q-q plots as described above based on the quantiles of spine counts vs. distance from the soma and found that upon visual inspection almost all dendrites follow the theoretical uniform distribution closely enough to assume that the density of the spines was homogeneous and therefore the csr case was a viable null hypothesis.
we selected  <dig>  example dendrites and their q-q plots are shown in figure
<dig>  we randomly selected  <dig> dendrite from each div and each biological replicate  to ensure the diversity of the set.
the y = x line is marked in red, and the observed q-q values are marked as black circles.
note that because this is a graphical method for comparing two probability distributions there was no p-value or significance level associated.
of all the  <dig> dendrites analyzed, only three of them  were considered non-csr at the 5% significance level.
figure
<dig> shows histograms of the p-values of all  <dig> dendrites separated into each div and experiment number.
the 5% significance level is shown by the red vertical line in each case.
we then computed the q-values for each dendrite and found that they are all equal to  <dig>  this is not surprising according to the explanation of the q-value above.
recall that q-values equal to  <dig> imply that 100% of the significant tests resulted in false positives, i.e.
there were no significant tests.
we therefore conclude that regardless of the maturity of the neuron, or the variation over biological replicate experiments, the locations of spines along all of the dendrites we analyzed were completely spatially random.
as mentioned above, the k-function is a function of the inter-point distance, t, that we consider around each observed point.
the range of t-values is determined by the total length of the network ℓ
t
, therefore because each dendrite has a different network length it also has a different range of t-values.
our chosen summary statistic throws away this information by computing the maximum absolute deviation  over all t in order to determine whether that value deviates significantly from the spatially random case.
however it may be of interest to determine whether clustering or repulsion between spines occured at specific inter-point distances t. ang’s correction normalizes the k-function such that the theoretical k = t for all t, so we can easily use this as a reference point.
figure
<dig> shows the k-function for the same  <dig> example dendrites used for the q-q plots of figure
<dig>  each graph shows the observed
k^ function , the theoretical k function  as well as the two-sided 5% and 95% point-wise simulation envelopes as a function of the radius t. following the description of the point-wise simulation envelope above we calculated these lower and upper envelopes at the 5% and 95% percentiles per t-value in the interest of checking if any t-value fell outside of this range.
since the black curves do not leave the gray shaded area for any value of t, the deviation from spatially random was insignificant at the 10% level for every t-value and is in agreement with our previous conclusion using the mad statistic.
this observation holds for almost all of the  <dig> dendrites we inspected visually, with no specific t-value evidencing either repulsion or clustering.
<dig>  we randomly selected  <dig> dendrite from each div and each biological replicate  to ensure the diversity of the set.
each graph shows the observed
k^ function , the theoretical k function  as well as the two-sided 5% and 95% point-wise simulation envelopes as a function of the radius t. we see here that the black curves do not leave the gray shaded area for any value of t, which means that the deviation from spatially random is insignificant at the 10% level for every t-value.
CONCLUSIONS
the models used in this work allow spatial prediction of spine types, which has not previously been studied.
the conclusions presented here relate to qualities of neurons in dissociated culture.
we acknowledge that some of these results will most likely not hold for in vivo settings due to neuronal interactions not modeled here, but maintain that the statistical methods used here will be useful and easily applicable.
specifically, we found here that spine type and density are not dependent on the distance from the cell body, and these observations are likely to change for in vitro slices or micro-injection of fixed brain tissue.
we also note that we chose not to deconvolve our data because of its high contrast.
we acknowledge that this choice may have precluded the image analysis software from detecting some stubby spines among the halo of the bright dendrites, but we do not feel this significantly impacted our results.
as a partial compensation for this effect we used neuronstudio’s in-built automatic z-smear compensation, and for more details on this we refer the reader to
<cit> .
although in this study the spine distributions seemed to be completely spatially random it is possible that we will find studies using different neuronal types and treatments where this is not true.
in these cases, where spine density may vary with distance from the cell body, it would be interesting to test for inhomogeneous patterns of points such as the hard core strauss process used in
<cit> .
we could also place an exponentially decaying function to model the interaction between spine types within a certain radius or experiment with other pairwise interaction functions such as those used by diggle, gates and stibbard
<cit>  or diggle and gratton in
<cit> .
we find it an interesting result that spines were not spatially clustered when type was disregarded, as shown by the linear network k-function analysis, however spine types do tend to group together as shown by the mlr analysis.
we would like to note that these results are not contradictory because they are in fact measuring different quantities.
the mlr results tells us that, regardless of their densities along the dendrite, if we have a spine which is of a given type, its  <dig> nearest neighbors are likely to be of the same type.
the k-function, on the other hand, tells us that regardless of type the spines’ locations along the dendritic network are spatially random.
these two results provide complementary information and together could aid us in future modeling tasks such as simulation of neuronal growth.
for example, we could first place spines uniformly along the dendritic network, and then decide the types of those spines based on the type of information given by the mlr model.
as future work we plan to analyze the network cross k-function
<cit>  of the dendritic network, which models the spine distribution as a multi-type point process and therefore provides information about repulsion and clustering of each spine type with each other spine type, modeling both density and type simultaneously.
generally previous studies such as
<cit>  have relied on physiology or biochemical markers to validate their neuronal properties.
the quantitative morphological features described here provide an additional phenotypic dimension for these analyses.
likewise these approaches can be applied to phenotypic analyses of neuronal cultures following over-expression or suppression of specific genes to capture their effect on a complex phenotype.
as mentioned in the introduction section, the only other study we are aware of which analyzes clustering of dendritic spines in monkey brains is
<cit> .
the authors of this work study the number of "clustered spines" on each dendritic segment, where a cluster is defined as a group of  <dig> or more spines.
the method used here defines clustering as a statistically significant positive deviation in the linear k-function from the theoretical value of the spatially random linear k-function.
we believe our method to be more principled and our results easier to interpret than those of
<cit>  due to the more formal statistical definition of clustering.
we chose to use dissociated hippocampal cultures because they are widely used and they allow us to perform an in-depth and automated analysis with larger spine populations than most previous studies.
these approaches will be important in assessing features of neurons derived from human induced pluripotent stem cells which have so far not been characterized by detailed morphological features.
our paper utilized a highly simplified neuronal culture system to develop the statistical and computational tools for more advanced in vivo studies needed to address the aforementioned bigger biological questions.
our overall hypothesis was that we can utilize imaging and statistical analyses to capture features of spine distributions that can be used for testing hypotheses in in-vivo settings.
indeed, we have been conservative about hypotheses and findings concerning spine type clustering because any conclusions we might reach on the specifics of spine distribution would be limited to the neuronal culture system we studied.
availability of supporting data
all the image stacks and neuronstudio annotation files supporting the results of this article are available in the bisque repository,
http://bisque.ece.ucsb.edu/client_service/view?resource=http://bisque.ece.ucsb.edu/data_service/dataset/ <dig>
abbreviations
div: day in vitro; llm: log-linear model; sd: distance from the soma to the spine along the centerline of the dendrite; bo: branch order of the dendritic segment on which a spine lies; n <dig>  n <dig>  n3: the spine types of the 1st, 2nd, and 3rd nearest neighbors along the dendrite, respectively; mlr: multinomial logistic regression; csr: completely spatially random; fdr: false discovery rate; mad: maximum absolute deviation.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
aj designed the study, wrote the software to perform the statistical analyses, interpreted the results and drafted the manuscript.
sb cultured and collected the data and drafted the section on cell imaging.
kk conceived of the study and helped to draft the manuscript.
bm coordinated the effort between departments and reviewed the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
table of 4-way llm coefficients.
this table shows the 4-way interaction llm coefficients which are significant at the  <dig> % level.
note that only one interaction between type and either branch order or soma distance  is significant in the entire table.
this further proves the result that these interactions are not very important to the overall model of frequencies.
click here for file
additional file 2
aic stepwise models for 3-way llm.
this table shows the results of the aic stepwise algorithm using an llm with up to 3-way interactions.
the models arrived at by this method are shown in the caption above the table.
from this table we can see that if we do allow 3rd order interactions, the strongest 3rd order correlation over all experiments is that of div, sd and bo, which makes sense because all three of these quantities should intuitively increase together.
click here for file
acknowledgements
we would like to acknowledge sruti aiyaswamy for her help with the neuron annotation, and the statlab at ucsb for checking over the statistical analyses included here.
this work was supported by the larry l. hillblom foundation, the tau consortium and an nsf award iii- <dig>

 spines are protrusions that occur on the dendrites of most mammalian neurons. they contain the post-synaptic apparatus and have a role in learning and memory storage. spine distribution is a critically important question for multiple reasons. changes in spine distributions and shape have been linked to neurological disorders such as fragile x syndrome. spine distributions determine the extent to which the neuropil will be electrically sampled, i.e. dense distributions will sample the neural connectivity map more fully. furthermore, the nature of optimal sampling is unknown and likely depends on the surrounding anatomy and the total information content available to dendrites. because pruning takes place during development in an activity dependent manner, spine distributions may reflect activity within neural circuits. distributions of spine types are biologically important because the electrical properties of spines, such as the spine neck resistance, promote nonlinear dendritic processing and associated forms of plasticity and storage to enhance the computational capabilities of neurons. the shapes and types of dendritic spines contribute to synaptic plasticity. because neighboring spines on the same short segment of dendrite can express a full range of structural dimensions, individual spines might act as separate computational units. nevertheless, the dendrite acts in a coordinated manner and thus the spatio-temporal distributions of different spine types is likely to be significant. little is known about this population level organization of dendritic spines. our aim was to perform an exploratory analysis of neuronal data from different time periods during the growth of rat dissociated hippocampal neurons, a well-established model system. the observations here pertain only to the culture system and not necessarily to in vivo settings although the analytical tools used here could be adapted to in vivo analyses. by quantifying populations of dendritic spines with automated tools at a global level, we were able to perform a much larger and more comprehensive analysis than most previous studies. many studies only analyze a small region of interest on the largest dendrites, for example the 50– μm closest to the soma, or μm segments, making it easier to measure manually the spine type counts and dimensions. other works determine spine lengths and widths by manually drawing a line along the maximal length and measuring the length of that line, and therefore are only able to analyze a few neurons and a few hundred spines at a time. in this study we determined the ratios of spine types along the dendrites as a function of time in culture, clustering or repulsion of spines in space, and how best to model spine type distributions. a model that fits the spatial distribution of spine types in healthy cultured neurons would be useful to assess neuronal cultures following gene manipulation such as rnai and to study features of induced pluripotent stem cells differentiated to neurons. log-linear models and multinomial logistic regressions are two basic and essential statistical methods, and have an extensive history of being used in biological studies. however, these tools have not been used thus far in the analysis of spine distributions. we fit a log-linear model to the contingency table of spine features to determine the dependence between spine types, distance from the cell body along the dendrite, the branch order of the dendritic segment on which it lies, and the day in vitro on which it was imaged. once we determined which of these attributes contributed to the overall dendritic spine model, we then asked whether these attributes can predict the occurrence of spines and of spine types. to answer this question we used a multinomial logistic regression model, which predicted the spine type, using the attributes that were found to be important through the llm and associated contingency tables. finally, to understand how the dendritic spine density varied over the length of the neuron or whether the appearance of spines was completely spatially random i.e uniformly distributed over the neurites, we made use of spatial point processes. spatial point processes have been used before in biological studies to model the locations of entire neurons, locations of ants nests or xylem conduits. there have also been other more ad-hoc methods created to study the number of "clustered spines" on each dendritic segment in monkey brains, where a cluster is defined as a group of or more spines however we believe our use of the linear network k-function is the first work to analyze the locations of dendritic spines and their clustering properties in such a principled manner. our analysis indicated that the density of spines is generally completely spatially random over the dendritic length probably due to the absence of instructive directional signals found an in vivo setting in which spine distributions are unlikely to be csr. methods cell imaging dissociated hippocampal neurons from embryonic rat brains were plated onto poly-l-lysine coated coverslips. once neurons adhered to the coverslip, they were placed face-down on glial cells grown in vitro for days. these neurons were a primary neuronal culture system, and no cell line was used. neurons were grown for specific time periods up to days in a neuronal medium containing b this co-culture of neurons and glia mimic the physiological conditions of neuronal growth and development in mammalian brain. work with the neuronal cultures was approved by the ucsb animal care committee. to fill the neuronal processes including dendritic spines green fluorescent protein was expressed from a plasmid containing the beta-actin promoter. of this plasmid, μg was transfected into each coverslip containing about neurons. transfection was performed as described in the manufacturer’s protocol with minor changes. the transfection mix and neurons were incubated for two hours to avoid toxicity caused by lipo following transfection, coverslips were flipped back onto the glial dish, where they were originally cultured. gfp-actin transfected into the neurons at div and neurons were studied at three time points- div and these time points survey the maturation period over which synapses and spines emerge. note that these were not the same neurons studied over time, but each time point represents a different population of neurons which were grown in culture up until the point of imaging. in this way our analysis represents a study at the population level. at each time point the number of images taken per plate depended on the transfection efficiency of that plate. on average approximately 1% of cells were transfected. the plating density was set so that neurons were relatively isolated in order to capture one neuron per image. an olympus fluoview laser scanning confocal microscope was used. image slices were by pixels at 154nm per pixel resolution. there were 7– z-slices per stack depending on the depth of the neuron, taken at 200nm steps. this means that the stacks were μm ×  μm ×  – μm. the z dimension slices were used to capture each depth level at the optimal focus, however we cannot claim to have accurate volumetric information at this resolution. a 40x oil objective lens with no optical zoom was used. numerical aperture was, and illumination conditions were kept constant. deconvolution of the raw data before processing was not necessary because the images were clear enough to manually annotate the neuron traces and manually edit all the spine detections and types as described in the following section. we performed three biological replicates, the results of which are detailed below. although there are other higher resolution, full volume methods, the analysis of this data is broadly applicable to imaged neurons in other systems. we attempted to capture the entire neuron in each image, however because of limits in available imaging techniques we found that this does not always happen. in the cases where dendrites were truncated at the end of the image plane we assumed that the proportion of spines in the missing data was similar to what had already been observed, and therefore the resulting distributions did not change. we verified this assumption visually by taking tiled mosaics of a few neurons imaged in their entirety from each div and checking that the branch orders, distances to soma and spine type counts were unchanged as compared to those of the same div. there was an observed increase in the dendritic length truncated by the image plane as the div increased. however in our particular analyses the methods used, such as the log-linear model and multinomial logistic regression, were focused on trends between spine characteristics such as distance to soma and type and these trends are innately unaffected by the truncation of dendrites given the above assumption. in addition, spatial point process analyses such as the linear network k-function always include the specification of an observation window, which in our case was the image plane. we verified that the overall spine density and the density of each spine type did not vary with distance from the soma so that we could assume spine density at the ends of the dendrites which were truncated was similar to the dendritic length which was observed. we recognize that we cannot see the proximity of labeled cells to other neurons which haven’t taken up the gfp labeling. these unlabeled neighboring neurons may cause some difference in spine distributions which we cannot quantify. for this reason we have attempted to quantify our biological findings statistically over entire experiments and div time points instead of by individual neurons, although in certain cases showing results from individual randomly sampled neurons was necessary. neuronal reconstruction there exist many automated methods for studying neuronal growth and morphometry and therefore we present a brief review of available software for tracing dendrites and detecting and classifying spines. in particular, neuronj is a widely used software; however it is only semi- automatic and one must click several points to trace each neurite. the labeling is done manually and the statistics output only include lengths of neurites and not spine data. hca-vision is a costly software with similar goals, however the parameters of the neurite tracing are set manually with a sliding bar and thus results require much hand-tuning. in addition, it is also focused on tracing neurites as opposed to spine analysis. for a full review of existing methods and softwares for neuron tracing and spine detection see. we found neuronstudio to be the most user-friendly, and for this reason we used it to annotate dendrites and spines for this analysis. despite the abundance of automated softwares, neuronal reconstructions are still largely performed by hand and this is is especially essential for a study like this one, where the traversed distance of the dendrites and number of spines and their shapes were analyzed in such detail. using automated reconstruction algorithms on raw data is prone to both false positive and false negative detections of spines, as well as misleading spine shape measurements. in cases where neurites from neighboring neurons enter into an image, neuronstudio often incorrectly traces these neurites as belonging to the neuron of interest. for this reason we manually traced each dendritic branch and soma of each neuron, ran neuronstudio’s automated spine detection/classification algorithm and then manually inspected and verified each spine’s location and type. the verification and tracing were done by the primary author and an undergraduate biology student working in the kosik lab. they were both familiar with dendrite and spine morphology and the resulting annotations from each were cross-checked by the other. relevant spine attributes output from the neuronstudio software include branch order, type, distance to soma along dendrite, length and width at widest point. however since neuronstudio uses the length and width of the spines to determine the spine type, we chose to make use of spine type and discard the other measurements. neuronstudio uses centrifugal labeling for branch orders, meaning it starts at at the cell body and moves outwards, incrementing at every y-shaped bifurcation regardless of the diameter of the daughter branches. note that the entire image stack with z-dimension information was loaded into neuronstudio for the spine classification, and that the software has interpolation algorithms to estimate the spine type in 3d. for spine detection the default cut-offs were used, i.e. a required spine height between – μm, a maximum spine width of μm, a minimum stubby size of voxels, a minimum non-stubby size of voxels, and automatic z-smear compensation. for spine classification, the default settings were also used, i.e. a head-to-neck ratio threshold of μm, an aspect ratio threshold of μm and a minimum mushroom head size of μm. neuronstudio delineates spine types by these thresholds. it is generally known that mushroom spines have a large head and a narrow neck, thin spines have a small head and a narrow neck, and stubby spines display no obvious subdivision in head and neck. if the head-to-neck ratio is above the threshold and the minimum mushroom head size is met, the spine is considered mushroom. if both the head-to-neck and aspect ratios are lower than the respective thresholds then the spine is considered stubby. the remaining cases result in thin spines. for further information on neuronstudio reconstruction, detection, and spine classification algorithms please refer to. in addition to the spine information, a trace file is output which labels the cell body, branch points and end points of the dendrites. the trace provides a skeletonization, or centerline, of the dendrite which we used to compute the linear network distances in the following analyses. log-linear model as a tool for exploring important features and their dependencies to find the most influential attributes with regard to prediction and spatio-temporal modeling of spines we fit a log-linear model to the feature data, which is a type of generalized linear model. the co-occurrence frequencies of the features in question are essentially a large multidimensional contingency table of counts. the standard linear models assume that data is normally distributed around a certain mean, which means that the observations can take any real value, positive, negative, integer or fractional. log-linear models, on the other hand, assume that data is intrinsically non-negative, typically counts that could be poisson distributed, and allow us to model the association and interaction patterns among categorical variables. the attributes under consideration are bo, type, sd and div. again, since the type of spine was quite directly dependent on the length and the head diameter of the spine, we left these latter variables out of the modeling. in order to analyze the data using a log-linear model, the various features must be in a categorical form or discretized. in an exploratory analysis such as this, one does not know what dependencies among features to expect; however we would like to note that these dependencies were not lost in the discretization process since trends in increasing and decreasing feature values would be preserved. to ensure that there were a reasonable number of observations at the higher branch orders, we pooled bo values of or higher into a single category called "higher-order branches". we created a categorical variable to represent the continuous variable soma distance where categories were determined using the quartiles of the sd spine data pooled over all experiments. specifically, sd values of less than μm were classified into the first group, from this value to less than μm the second, from this to less than μm the third, and the rest fell in the fourth group. binning the observed data for the continuous variables is the best way to get a general feel for how these quantities relate to each other. after this post-processing of the data we arrived at categories of branch order, categories of soma distance, spine types, and divs. using the observed frequencies for the aforementioned attributes, we created a four-way contingency table and fit the model using the'glm’ function in the r package'stats’. the table of the frequency of occurrences of the four attributes was modeled as poisson with each entry being a simple count of the co-occurrences of that bin. we called this count f ijkl with each of the subscripts i,j,k,l corresponding to a different attribute. the method uses the link function y ijkl  = log, and treats the model as a regular linear model. each entry y ijkl is modeled by a combination of coefficients: the intercept, plus main effects, plus every combination of interactions between these four attributes, as shown below. yijkl=μ+αi+βj+γk+δl+ij+ik+il+…+error. we estimated this full interaction model using the least-squares maximum-likelihood approach. we also used a stepwise fit algorithm, which begins with a model that includes only the constant term, and at each step chooses whether or not to add one additional term. the algorithm begins with the main effects then tries each possible 2-way interaction, aiming to minimize the akaike information criterion. the aic is defined as aic=2k-2ln) where k is the number of parameters i.e the total number of coefficients being estimated, and l=maxθ∏n=1neynθ′xne-eθ′xnyn! is the maximized value of the likelihood function for the estimated poisson model. in the above equations x=x …,xn∈r are the input vectors, θ  = θ1…θ k are the parameter values, and y=y …,yn∈r is the output. the aic is a commonly used goodness-of-fit measure for a model given the observed data. adding or subtracting terms, whether they be main effects, pairwise interactions, or up to 3-way interactions between attributes, will change the aic value for the model. a lower aic criterion indicates a better fit to the data and therefore a better model. to compute the stepwise fit we used the r function'step’. for more information on the stepwise fit algorithm as well as the aic criterion we ask that the readers refer to the'step’ function reference. we ran both of these llm fitting procedures for all experiments separately expecting to find general agreement between coefficients of the corresponding models created. multinomial logistic regression to predict spine type from neighbor types in order to predict spine type we first determined which attributes contributed most to spine type prediction. given the complexity of the multidimensional llm and the various interactions and conditional frequencies that would impinge on this issue, we decided to determine these attributes by analyzing 2-way contingency tables for spine type vs. sd, bo, div, as well as the spine types of the nearest neighbors. this analysis helped us pick attributes that would be useful as the predictors in the multinomial logistic regression explained below. when the response variable of a regression takes binary values "logistic regression" is used. this is an approach which uses a linear combination of the predictor variables to predict the log-odds of a success. since our response variable was spine type and it can take values, we needed to use a "multinomial logistic regression" which attempts to model the probability of any of multiple possible outcomes. we did not use the attributes sd or bo as predictors variables since the results of both the llm analysis and 2-way contingency tables mentioned above told us that these quantities were not as relevant for spine type prediction. therefore our model consisted of spine type as the output variable and the div, 1st, 2nd and 3rd nearest neighbor type along the dendrite as the predictor variables. we tried using only or nearest neighbors, however the results proved inconclusive because the prediction probabilities for each of the types were predominantly close to 1/ if we used more than the nearest neighbors we sometimes ended up spanning a segment of dendrite which we did not consider to be "local", so we decided that nearest neighbors provide the most useful information in the case of this study. the mlr analysis we performed in this paper does disregard the actual inter-spine distances, meaning that if the nearest neighbors are very close or very far apart we still treat them the same. we did this partially because adding the distance variables would complicate the model significantly, but also because we believe that over a large population of spines such as the one we have, these differences in distance will average out and we will still get a general picture of the trends between neighboring spine types. to verify that this was true we computed a histogram showing the distribution of 3rd nearest neighbor distances for each spine, shown in figure although the maximum distance to any 3rd nearest neighbor is extremely high we can see from the histogram as well as the fact that the median 3rd nearest neighbor distance was μm that this distance is clearly an outlier case and that the majority of 3rd nearest neighbor distances lie below μm. suppose the output variable categories are denoted by, corresponding to mushroom, stubby or thin spines, with being the reference category. if y i denotes the observed outcome of the output variable, and x i is the corresponding vector of the neighbor types and div for the ith observation, one regression is run for the logit probability of each category k, with β k representing the vector of regression coefficients in the kth regression. this is done for all but the reference category, whose probability is then obtained by subtracting all other probabilities from one. note that because the predictor variables were spine types, which were nominal as opposed to ordinal variables, the predictor variables x i must be represented with a "dummy coding". this means each neighbor type was represented by predictor variables, where corresponded to mushroom type, corresponded to stubby type and corresponded to thin type. this does not need to be done for the output variable y. with the addition of the div, which does not have to be dummy coded since it is an ordinal variable, this made each x i vector of length the regressions are then written as: p=exp1+exp+exp p=exp1+exp+exp and p=1-p-p=11+exp+exp the parameters are estimated typically by using an iterative procedure such as "iteratively re-weighted least squares" or, more commonly by a numerical approach such as the broyden-fletcher-goldfarb-shanno method. in our case we create an mlr using the command multinom in the r package nnet which uses bfgs by calling the r function optim. it can be seen that logp)=β1xi logp)=β2xi so that the beta coefficients represent the change in the log odds of the dependent variable being in a particular category with respect to the reference category i.e. the thin type, for a unit change of the corresponding independent variable. to check if the models created from all three experiments were in agreement, we ran the mlr separately for each experiment. to satisfy one of the major assumptions of this analysis, namely that the data must be a set of independent observations, we took randomly sampled spines of each type from each experiment to use for the parameter estimation. we chose to select equal proportions of each spine type in order to remove any bias in the model towards the less frequent thin spines, and was the largest number we could justify using since there were only thin spines in experiment we verified that these randomly sampled spines did not lie within μm of the image border so that we were fairly certain their nearest neighbors did not fall outside of the image plane. note that due to the tortuosity of the dendritic structure this did not mean that our sample was necessarily biased towards spines which were proximal to the soma. we did not verify explicitly that the sampled spines were not neighbors of each other, since we assumed that the variation captured by the random sampling was enough to ensure some level of independence. the idea was to aim for an independent set of observations which represented the entire "population" of spines in that experiment. to be clear we used all spines for the llm model and k-function analysis, only the mlr model required random sampling since we were using neighbor information which would have been redundant if we considered every spine. to verify that the prediction of spine type provided by the mlr was better than what we would get purely by their relative abundance i.e. without neighboring spine type information, we computed something similar to a "bayes factor". bayes factor is a method of choosing between two models on the basis of the observed data. in our case, the first prediction model was simply the prior global probability of finding a given spine type based on its frequency in the particular experiment under consideration. the second model was the mlr prediction model using the neighbor type information. we computed p/p and reasoned that values considerably larger than one indicated the neighboring spine type information was helpful in the prediction of the central spine type. linear network k-function as a tool for testing spatial point patterns originally proposed by ripley in 1981, the purpose of the k-function is to estimate whether or not there is clustering or repulsion present in a given spatial point process. the common null hypothesis is that the points within the observation window are distributed as a homogeneous poisson process, which is also termed "completely spatially random" or csr. this means that the density of points does not vary depending on the spatial parameters i.e. x and y in the 2d euclidean case, or the location along the dendritic network in our case. in order to determine if this is a valid null hypothesis for our data, we created q-q plots for individual dendrites which compared the quantiles of the sd values of observed spines to the theoretical quantiles for the csr case. if the two distributions being compared were similar, the points in the q-q plot would approximately lie on the line y = x. in order to create the theoretical quantiles it is necessary to know the values of sd at any location on the given network, not just at the spine locations. once we have this we can partition the network into epsilon small segments and assign each segment a value if it contains a spine and otherwise based on the csr assumptions. we did this using code provided to us by adrian baddeley and gopal nair at the commonwealth scientific and industrial research organization, australia. the k-function computes the expected number of points within a distance t of an arbitrary point p, therefore the empirical value in 2d euclidean space for the csr case will be proportional to the circular area, λπt the proportionality constant λ represents the density of points in the homogeneous poisson case, and can be estimated by finding the total number of points n divided by the total area of the observation window a. ripley’s k-function, which is a function of t, is a very useful tool because it describes the 2 nd order characteristics of the point process at several scales t. if we ignore the edge effects due to the observation window, the observed k^ can be written as: k^=|a|n2∑i∑j≠ii where i stands for the indicator function, and d ij stands for the euclidean distance between two points p i and p j. in the above equation, we see that the expectation is normalized by 1/λ since λ=n|a|, so we infer that theoretically k = πt implies spatial independence of points, or a csr point process. therefore, if k is the theoretical csr value of the function and k^ is the observed function, then k^>k implies clustering between points and k^<k implies repulsion. it is possible to extend this function to multi-type point patterns or to higher dimensional data. since our particular point process consists of spines which lie along the "linear network" of the dendritic tree we were primarily concerned with inter-spine distances along the dendrite as opposed to in euclidean space. therefore we used a version of the k-function developed recently for linear networks by okabe and yamada. this modified version of the k-function takes into account the structure of the linear network on which the point process resides and imitates the euclidean space k-function described above. the linear network k-function is calculated as follows: k^=ℓtn2∑i=1n∑j≠ii where ℓ t is the length of the total network l t. the theoretical csr for this case is described as follows: k=1ℓt∫p∈ltℓpdp where p is a point belonging to the set of all points p = {p …,p n }, and ℓ p is the length of the subset of the network l p where the distance between p and any other point is ≤ t. note that here the distance d ij stands for the linear network distance along the dendrite. accounting for variability in the length ℓ p means the formula takes into account the edge effects due to the observation window inherently, but at the cost of added complexity. the computation of the theoretical linear network k-function requires us to find lpi, the subset of l t where the network distance between a specific point p i and any other point is ≤ t, and ℓpi, the length of that subset, for every point p i. a visualization of the quantities d ij, l t, ℓ t, lpi, and ℓpi is shown in figure lpi, and ℓpi which were used to compute the linear network k-function. here dij is the linear network distance shown by the gray line between points pi and pj. lt is the entirety of the single dendritic network and ℓt is the length of lt. similarly, lpi is subset of the network where the distance between a point pi and any other point is ≤ t and ℓpi is the length of lpi. in this particular example there are spines which fall within lpi and would be counted in determining the empirical function value k^, however point pj falls outside this radius and would therefore not be counted. note that although many biological applications of point processes treat individual observations as replicate patterns coming from the same underlying distribution, we cannot do that using the above definition of the network linear k-function due to the change in linear network structure from dendrite to dendrite. the term "dendrite" here refers to the entire dendritic tree resulting from a single root branch of a neuron. other in-vivo studies focus on clustering of spines which lie on the same unbranched section of the dendrite, however we focus on the entire dendritic tree under the hypothesis that it follows rule-based distributions of spines due to anatomical constraints and integration of the a signal over the entire dendrite. one can infer from figure that since the geometry of the linear network changes from dendrite to dendrite, so do the total lengths of the networks ℓ t, the ranges of possible t-values and the amount of dendritic length that is present within a given distance of any point. we did not simply normalize the lengths of the networks to a scale because it is desirable for the t-axis to retain its real physical values in order to make conclusions about the scale of clustering or repulsion among spines. however, we did desire to compare the linear network k-functions of various dendrites in a meaningful way. for this reason we used a corrected version of the network k-function that intrinsically compensates for the geometry of the network called ang’s correction. the observed k-function then becomes: k^=ℓtn∑i=1n∑j≠iim where m is the number of points of l lying at the exact distance t away from the point i measured by the shortest path. that is, the contribution to the function from each pair of points is weighted by the reciprocal of the number of points that are situated at the same distance from i as j is. as a result, the theoretical csr case is simply k = t for all 0 ≤ t < t. this enables direct comparison of t-values across dendrites, as we will see in the results section. simulations and q-values to test the null hypothesis that the locations of spines on the dendrites were indeed csr, we created a summary statistic which encompasses the difference between the empirical k^ and the theoretical k under csr. the summary or "test statistic" we used, is the max absolute difference over t, viz. d=maxt|k-k^|. one method for obtaining a distribution of d proposed by diggle is to bootstrap the residuals, or differences between the observed and theoretical values. however a more heuristic and intuitive way is to simulate the csr case for each dendrite, compute the k-function for each of these simulations, and find the simulated distribution of our test statistic. we then found the p-value of the observed difference d from this simulated distribution. specifically, we carried out csr simulations for each dendrite by placing uniform points on a line [ ℓ t ], and mapping them to that specific dendrite’s linear network structure. the number of points simulated per dendrite equaled the number of observed spines for that dendrite, thus preserving the overall density λ. this means the same number of spines that existed on each dendrite were randomly placed along the linear network specific to that dendrite. we used these simulations to obtain values of the summary statistic, say d. then the p-value for each dendrite was simply the proportion of simulated values that fell above the observed or experimental value of d, i.e. the rank of this d within the values of d, or nrank/. this p-value approach is similar to the test which rejects the null hypothesis if the graph of the observed k-function lies outside the "point-wise simulation envelope" at any value of t. a simulation envelope is essentially a graphical measure of how far a function can deviate from the theoretical value without being considered significant at a given level. as mentioned above in our case the envelope is calculated by first creating the csr simulations of a point pattern on a given dendritic network with the same observed network intensity, then calculating the linear k-function for each of these simulations. to perform a two-sided significance test at the 10% level, the 5% and 95% percentiles are then calculated based off the lowest and highest linear k-function values per t-value, hence the term "point-wise". plotting these values as a function of t gives one a visual idea of the spread that is produced by chance mechanisms alone. if the observed k-function for a given t-value does not fall outside these percentiles, it is considered insignificant for that t-value at the 10% significance level. we make use of the r package'spatstat’ for obtaining the point-wise simulation envelope. because we have a multitude of hypothesis tests and p-values, to reach a conclusion about the general trend for each div and experiment, we used the concept of false discovery rate. the fdr is defined as π0=#truenulltests#totaltests controlling the overall fdr, or expected proportion of incorrectly rejected null hypotheses termed "false discoveries", is a statistical method commonly used in multiple hypothesis testing which increases the statistical power of each test. what is more general and useful however, is a test-specific fdr measure. this essentially allows us to look at all possible significance thresholds at once, as well as provide each test with a measure of significance that can be easily interpreted. this is accomplished by calculating an analogue of the p-value for each test called a "q-value". a p-value of implies that 5% of all tests will result in false positives, whereas a q-value of implies that 5% of significant tests will result in false positives. since the latter is clearly a far smaller quantity, q-values generally indicate fewer significant tests than p-values for a given significance threshold and provide a far more accurate indication of the level of false positives in the case of multiple hypothesis testing. for q-value estimation we used the qvalue package available from.  data analyzed we performed three biological replicate experiments resulting in a total of neurons from the following time points: div div and div. this provided a rich and complete data set resulting in dendritic branches and spines. example images from each div along with zoomed in dendritic segments where spines and annotations are visible are shown in figure scale bars are shown in red in panels a-c and the yellow rectangular boxes in panels a-c show the region of interest which has been zoomed in on in panels d-f respectively. panels d-f are all at the same resolution. the number of spines per μm, or λ, for each dendrite in different experiments and time points is shown in figure we chose to include this in order to help the reader compare these neuronal culture results with other experimental paradigms with which they may be more familiar. it is clear from the histograms that the distribution of spine density for div is skewed toward lower values as compared to the density for div as expected. the image data as well as spine and trace annotations are made publicly available through the bisque system and the url is given in the section titled "availability of supporting data". we chose bisque over other databases like neuromorpho.org because it allows us to upload multiple layers of annotations as opposed to only the digital reconstruction files. we calculated a 2-way contingency table over all experiments and spine types and obtained table from this table we note the high frequency of mushroom and stubby spines as compared to thin spines, and also the fact that the ratio of types does not remain the same per experiment even though they were indeed biological replicates. in fact, a pearson’s chi-squared test on table shows dependence between the spine type counts and experiment number, χ2 = , p < . we believe that the large experimental variation between spine type proportions and counts in each experiment was a positive result because this meant that statistical agreement across all experiments relating to spine type clustering and density estimation carries heavier weight than if the experiments were more uniform in these quantities, or if we had pooled data from all experiments together. also, if all experiments were unusually homogeneous there could be a possibility that it is a result of our specific culturing, imaging or spine extraction methods rather than a true representation of the underlying biological process. the various biological systems to which these techniques will be applied will certainly have this type of variability. spine type is independent of distance from soma as described in the methods section, we calculated a stepwise-fit of the log-linear model starting with just a constant term, and at each step choosing to add the main effects and possible 2-way interactions between main effects one-by-one if they decreased the corresponding aic value. the captions above tables and show the final models arrived at for each of the experiments as well as their corresponding aic values. the tables indicate the change in the aic value that would occur from adding or omitting each of the terms in the first column. this gives us an idea of how important that term was to the model. the rows of the table are ordered by their overall contribution to the model, i.e. the term in the first column of the first row of each table had the lowest aic value and was therefore the most important to the overall model. if the reader requires further information on the aic criterion or how to interpret this table we ask them to refer to chapter of. despite the fact that they were included in the final stepwise fit model for experiments and the aic values in tables and show that in all experiments the interaction between spine type and soma distance as well as spine type and branch order were the least important in modeling the overall frequency table of occurrences. this implies that the correlation between these quantities was not very high, therefore we reason that it was not necessary to use either sd or bo to predict the spine type in the mlr created in the following section. we also noticed that the term marking the interaction between bo and sd was the most important pairwise term in all stepwise fit models. it is expected that bo and sd are correlated because both necessarily increase as we move away from the cell body. indeed, running a 2-way chi-square test on the contingency table of the discretized versions of these variables showed us high dependence, χ2 = ,p<. we also saw a high level of dependence between div and sd  = , p <  ) and between div and bo  = , p <  ). this was intuitive as well since we expect both bo and sd to generally increase with div. it is possible that the type vs. sd relationship could have also been estimated using a sholl-type analysis where we count the number of each type occurring within concentric circles from the soma and verify that it is constant, however this would not necessarily produce the same results. the crucial difference between our approach and the sholl approach is that in our approach the "distance from soma measures" the actual distance along the centerline of the dendrite instead of the radial distance from the cell center. this is especially important for dendrites with high tortuosity, since the radial distance in those cases will not correspond to the dendritic distance from the cell body. many studies of cultured neurons use sholl analysis, however they use it in its original form for counting dendritic intersections and do not comment on the relation to spine density or type. to our knowledge this is the first study to quantify the spine density vs. distance to the soma in dissociated neuronal cultures. three-way and 4-way interactions are generally known to be weak and difficult to interpret, however in the interest of exploring all possibilities we computed the maximum likelihood fit using all attributes as well as a stepwise fit model which allows for 3-way interactions between attributes. the table presented in additional file results from the llm which models all possible interactions of all attributes, i.e. up to the fourth order. the coefficients presented in the table are those which were significant at the % level, and the corresponding p-values are shown in the last column. the table contains the interactions which were more important to the model, and shows that of these interactions only one between type and either bo or sd, was shown as being significant over all experiments. this verifies once again that neither bo nor sd were highly correlated with the spine type. in addition to this, the stepwise fit models in additional file show that if we did allow 3rd order interactions, the strongest 3rd order correlation over all experiments was that of div, sd and bo, again affirming that all of these quantities should intuitively increase together. spines tend to cluster with other spines of the same type in creating a regression model, we first ascertain that the predictor variables used are not only useful in predicting the output variable, but also that they do not provide redundant information as this can throw off the model fitting process. using all spines in the dataset, we performed a chi-square test on the 2-way contingency tables of spine type versus binned sd and bo, div, and the types of the nearest neighbors as described in the log-linear model section above. due to the aforementioned dependence between the type and experiment number we performed the test separately for each experiment and the results are shown in table from the table we can see that the div and the nearest neighbors showed clear dependency with spine type in all experiments, whereas sd and bo showed independence at the 5% significance level in experiments and respectively. since we expected sd and bo to have a similar relationship with type due to the high correlation mentioned above, and we had found this was not a very strong relationship, we chose to use only div, n n and n as predictors for spine type in the mlr model. n n n the resulting beta coefficients for each of the predictor variables are shown in table here "n1-var1" refers to the beta coefficent of the first dummy variable for the type of the first nearest neighbor; "n1-var2" refers to the second dummy variable, and so on. the "mushroom" row is omitted because it is the reference category and its probability is obtained as shown in eqn. we computed the prediction probabilities for each spine type given each combination of neighbor types for each experiment separately to determine the agreement between experiments. a selected set of results are shown below in tables and the highest probability for each row is marked by an asterisk. note that in these tables all divs in all experiments predicted the spine type to be mushroom when its nearest neighbors were mushroom type, and stubby when the nearest neighbors were stubby type. thin types were the most probable when the three nearest neighbors were thin type in all but experiment div and div the probabilities for cases where all of the nearest neighbors were not of the same type have been omitted for brevity and because they did not show any clear trends. * denotes the highest probability per row. *denotes the highest probability per row. * denotes the highest probability per row. the bayes factor results in table show that the proportional gain in information for the spine type in question was always greater than one for the prediction of a particular type when the neighborhood types were all of that same type. due to the low frequency of thin spines, their corresponding bayes factors were higher than that of other types, meaning that their prediction probabilities benefit more than other types from neighborhood type information. dendritic spine densities are completely spatially random we created q-q plots as described above based on the quantiles of spine counts vs. distance from the soma and found that upon visual inspection almost all dendrites follow the theoretical uniform distribution closely enough to assume that the density of the spines was homogeneous and therefore the csr case was a viable null hypothesis. we selected example dendrites and their q-q plots are shown in figure we randomly selected dendrite from each div and each biological replicate to ensure the diversity of the set. the y = x line is marked in red, and the observed q-q values are marked as black circles. note that because this is a graphical method for comparing two probability distributions there was no p-value or significance level associated. of all the dendrites analyzed, only three of them were considered non-csr at the 5% significance level. figure shows histograms of the p-values of all dendrites separated into each div and experiment number. the 5% significance level is shown by the red vertical line in each case. we then computed the q-values for each dendrite and found that they are all equal to this is not surprising according to the explanation of the q-value above. recall that q-values equal to imply that 100% of the significant tests resulted in false positives, i.e. there were no significant tests. we therefore conclude that regardless of the maturity of the neuron, or the variation over biological replicate experiments, the locations of spines along all of the dendrites we analyzed were completely spatially random. as mentioned above, the k-function is a function of the inter-point distance, t, that we consider around each observed point. the range of t-values is determined by the total length of the network ℓ t, therefore because each dendrite has a different network length it also has a different range of t-values. our chosen summary statistic throws away this information by computing the maximum absolute deviation over all t in order to determine whether that value deviates significantly from the spatially random case. however it may be of interest to determine whether clustering or repulsion between spines occured at specific inter-point distances t. ang’s correction normalizes the k-function such that the theoretical k = t for all t, so we can easily use this as a reference point. figure shows the k-function for the same example dendrites used for the q-q plots of figure each graph shows the observed k^ function, the theoretical k function as well as the two-sided 5% and 95% point-wise simulation envelopes as a function of the radius t. following the description of the point-wise simulation envelope above we calculated these lower and upper envelopes at the 5% and 95% percentiles per t-value in the interest of checking if any t-value fell outside of this range. since the black curves do not leave the gray shaded area for any value of t, the deviation from spatially random was insignificant at the 10% level for every t-value and is in agreement with our previous conclusion using the mad statistic. this observation holds for almost all of the dendrites we inspected visually, with no specific t-value evidencing either repulsion or clustering. we randomly selected dendrite from each div and each biological replicate to ensure the diversity of the set. each graph shows the observed k^ function, the theoretical k function as well as the two-sided 5% and 95% point-wise simulation envelopes as a function of the radius t. we see here that the black curves do not leave the gray shaded area for any value of t, which means that the deviation from spatially random is insignificant at the 10% level for every t-value.  the models used in this work allow spatial prediction of spine types, which has not previously been studied. the conclusions presented here relate to qualities of neurons in dissociated culture. we acknowledge that some of these results will most likely not hold for in vivo settings due to neuronal interactions not modeled here, but maintain that the statistical methods used here will be useful and easily applicable. specifically, we found here that spine type and density are not dependent on the distance from the cell body, and these observations are likely to change for in vitro slices or micro-injection of fixed brain tissue. we also note that we chose not to deconvolve our data because of its high contrast. we acknowledge that this choice may have precluded the image analysis software from detecting some stubby spines among the halo of the bright dendrites, but we do not feel this significantly impacted our results. as a partial compensation for this effect we used neuronstudio’s in-built automatic z-smear compensation, and for more details on this we refer the reader to. although in this study the spine distributions seemed to be completely spatially random it is possible that we will find studies using different neuronal types and treatments where this is not true. in these cases, where spine density may vary with distance from the cell body, it would be interesting to test for inhomogeneous patterns of points such as the hard core strauss process used in. we could also place an exponentially decaying function to model the interaction between spine types within a certain radius or experiment with other pairwise interaction functions such as those used by diggle, gates and stibbard or diggle and gratton in. we find it an interesting result that spines were not spatially clustered when type was disregarded, as shown by the linear network k-function analysis, however spine types do tend to group together as shown by the mlr analysis. we would like to note that these results are not contradictory because they are in fact measuring different quantities. the mlr results tells us that, regardless of their densities along the dendrite, if we have a spine which is of a given type, its nearest neighbors are likely to be of the same type. the k-function, on the other hand, tells us that regardless of type the spines’ locations along the dendritic network are spatially random. these two results provide complementary information and together could aid us in future modeling tasks such as simulation of neuronal growth. for example, we could first place spines uniformly along the dendritic network, and then decide the types of those spines based on the type of information given by the mlr model. as future work we plan to analyze the network cross k-function of the dendritic network, which models the spine distribution as a multi-type point process and therefore provides information about repulsion and clustering of each spine type with each other spine type, modeling both density and type simultaneously. generally previous studies such as have relied on physiology or biochemical markers to validate their neuronal properties. the quantitative morphological features described here provide an additional phenotypic dimension for these analyses. likewise these approaches can be applied to phenotypic analyses of neuronal cultures following over-expression or suppression of specific genes to capture their effect on a complex phenotype. as mentioned in the introduction section, the only other study we are aware of which analyzes clustering of dendritic spines in monkey brains is. the authors of this work study the number of "clustered spines" on each dendritic segment, where a cluster is defined as a group of or more spines. the method used here defines clustering as a statistically significant positive deviation in the linear k-function from the theoretical value of the spatially random linear k-function. we believe our method to be more principled and our results easier to interpret than those of due to the more formal statistical definition of clustering. we chose to use dissociated hippocampal cultures because they are widely used and they allow us to perform an in-depth and automated analysis with larger spine populations than most previous studies. these approaches will be important in assessing features of neurons derived from human induced pluripotent stem cells which have so far not been characterized by detailed morphological features. our paper utilized a highly simplified neuronal culture system to develop the statistical and computational tools for more advanced in vivo studies needed to address the aforementioned bigger biological questions. our overall hypothesis was that we can utilize imaging and statistical analyses to capture features of spine distributions that can be used for testing hypotheses in in-vivo settings. indeed, we have been conservative about hypotheses and findings concerning spine type clustering because any conclusions we might reach on the specifics of spine distribution would be limited to the neuronal culture system we studied. availability of supporting data all the image stacks and neuronstudio annotation files supporting the results of this article are available in the bisque repository, http://bisque.ece.ucsb.edu/client_service/view?resource=http://bisque.ece.ucsb.edu/data_service/dataset/ abbreviations div: day in vitro; llm: log-linear model; sd: distance from the soma to the spine along the centerline of the dendrite; bo: branch order of the dendritic segment on which a spine lies; n n n3: the spine types of the 1st, 2nd, and 3rd nearest neighbors along the dendrite, respectively; mlr: multinomial logistic regression; csr: completely spatially random; fdr: false discovery rate; mad: maximum absolute deviation. competing interests the authors declare that they have no competing interests. authors’ contributions aj designed the study, wrote the software to perform the statistical analyses, interpreted the results and drafted the manuscript. sb cultured and collected the data and drafted the section on cell imaging. kk conceived of the study and helped to draft the manuscript. bm coordinated the effort between departments and reviewed the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 table of 4-way llm coefficients. this table shows the 4-way interaction llm coefficients which are significant at the % level. note that only one interaction between type and either branch order or soma distance is significant in the entire table. this further proves the result that these interactions are not very important to the overall model of frequencies. click here for file additional file 2 aic stepwise models for 3-way llm. this table shows the results of the aic stepwise algorithm using an llm with up to 3-way interactions. the models arrived at by this method are shown in the caption above the table. from this table we can see that if we do allow 3rd order interactions, the strongest 3rd order correlation over all experiments is that of div, sd and bo, which makes sense because all three of these quantities should intuitively increase together. click here for file acknowledgements we would like to acknowledge sruti aiyaswamy for her help with the neuron annotation, and the statlab at ucsb for checking over the statistical analyses included here. this work was supported by the larry l. hillblom foundation, the tau consortium and an nsf award iii- 
although these results may vary with other systems, our primary contribution is the set of statistical tools for morphological modeling of spines which can be used to assess neuronal cultures following gene manipulation such as rnai, and to study induced pluripotent stem cells differentiated to neurons. we fit a log-linear model to the contingency table of spine features such as spine type and distance from the soma to first determine which features were important in modeling the spines, as well as the relationships between such features. our study indicated that in the culture system, dendritic spine densities were "completely spatially random", spine type and distance from the soma were independent quantities, and most importantly, spines had a tendency to cluster with other spines of the same type. a multinomial logistic regression was then used to predict the spine types using the features suggested by the log-linear model, along with neighboring spine information. dendritic spines serve as key computational structures in brain plasticity. our aim in this study was to perform exploratory analyses based on the population distributions of dendritic spines with regard to their morphological characteristics and period of growth in dissociated hippocampal neurons. dendritic spinesrat hippocampal culturelinear network k-functionmorphological modelingspatial statisticspoint processesneuronal growth finally, an important variant of ripley’s k-function applicable to linear networks was used to study the spatial distribution of spines along dendrites.
 dendritic spines serve as key computational structures in brain plasticity. much remains to be learned about their spatial and temporal distribution among neurons. our aim in this study was to perform exploratory analyses based on the population distributions of dendritic spines with regard to their morphological characteristics and period of growth in dissociated hippocampal neurons. we fit a log-linear model to the contingency table of spine features such as spine type and distance from the soma to first determine which features were important in modeling the spines, as well as the relationships between such features. a multinomial logistic regression was then used to predict the spine types using the features suggested by the log-linear model, along with neighboring spine information. finally, an important variant of ripley’s k-function applicable to linear networks was used to study the spatial distribution of spines along dendrites.  our study indicated that in the culture system, dendritic spine densities were "completely spatially random", spine type and distance from the soma were independent quantities, and most importantly, spines had a tendency to cluster with other spines of the same type.  although these results may vary with other systems, our primary contribution is the set of statistical tools for morphological modeling of spines which can be used to assess neuronal cultures following gene manipulation such as rnai, and to study induced pluripotent stem cells differentiated to neurons. dendritic spinesrat hippocampal culturelinear network k-functionmorphological modelingspatial statisticspoint processesneuronal growth 