Brain plasticity depends on the functional and structural reorganization of the synapses. The majority of the excitatory synapses are located on dendritic spines, which are small membranous protrusions localized on the surface of neuronal dendrites. The important feature of dendritic spines is their structural variability, which ranges from long, filopodia spines to short stubby and mushroom-shaped spines. Dendritic spines are typically composed of a head that is connected to the dendrite by a neck. The size of the spine head is proportional to the postsynaptic density area and correlates with the postsynaptic receptor content and synaptic strength (Nusser et al., 1998; Kharazia and Weinberg, 1999; Takumi et al., 1999). The length of the dendritic spine neck is correlated with the postsynaptic potential (Araya et al., 2006; Tønnesen et al., 2014). Thus, the dendritic spine shape has been accepted for determining the strength of the synaptic connections and is thought to underlie the processes of information coding and memory storage in the brain. Furthermore, alterations in dendritic spine shape, size, and density are associated with a number of brain disorders (DeKosky and Scheff, 1990; Irwin et al., 2001; Selemon et al., 2006; Knobloch and Mansuy, 2008; Sweet et al., 2008; Hutsler and Zhang, 2010; Penzes et al., 2011; Levenga and Willemsen, 2012).

The morphology of spines can change in an activity-dependent manner. The structural plasticity of dendritic spines is related to synaptic function, as the morphological modifications of pre-existing spines as well as the formation or loss of synapses accompany the learning and memory processes (Xu et al., 2009; Yang et al., 2009; for reviews see Holtmaat and Svoboda, 2009; Caroni et al., 2012). The cellular models of synaptic plasticity, long-term potentiation (LTP) and long-term depression (LTD) associate the synaptic strength with spine enlargement and spine shrinkage, respectively (Yuste and Bonhoeffer, 2001; Holtmaat and Svoboda, 2009; Kasai et al., 2010).

Understanding dendritic spine shape taxonomy and shape transitions upon synaptic potentiation is of great importance. The common approach in analysis of dendritic spine morphological changes is to categorize the spines into subpopulations based on whether they are stubby, mushroom, thin, or filopodia shaped (Su et al., 2014). Of importance, the existing categorization of dendritic spine shapes (Su et al., 2014) does not provide a clear definition of each group. Moreover, the literature lacks methods allowing an automated assigment between dendritic spine shapes. Additionally, a recent report suggests the existence of dendritic spine shapes continuum rather than specific spine subclasses (Loewenstein et al., 2015). Thus, the existing classifications have to be revisited and a new automatic classification method with clear mathematical rules should be derived. To address this issue, we employed a new methodological approach with potential applicability in the studies of dendritic spine shape taxonomy and transitions in time. Our clustering-based approach permits analysis of dendritic spine dynamics in multi-dimensional feature space by reducing the complexity subpopulations. To test this method, similarly to previous works, we potentiated the synapses with cLTP stimulation that produces a long-lasting increase in network activity and mimics several aspects of LTP, including synaptic receptor incorporation to the dendritic spine membrane. The morphology of single dendritic spines was assessed using time-lapse imaging of living neurons. In the rest of the paper, we refer to a population of spines stimulated by cLTP as ACTIVE, and the non-treated spines are denoted as CONTROL.

The paper is organized as follows. In Section 2, we describe the process of data gathering and data representation and the statistical approach to analysis of spine shapes. First, we analyze the basic characteristics of features in the populations named ACTIVE and CONTROL and conclude that before a meaningful comparison can be performed, the populations need to be normalized. In Section 3, we develop a simple but meaningful numerical representations of spines and we provide an approach to dendritic spine taxonomy construction and models of shape transitions together with statistical tests for model comparisons. For taxonomy development, we propose a clustering-based approach that does not depend on subjective decisions of experts and can accommodate arbitrary numerical features. Later, we introduce a corresponding probabilistic model of spine transitions between the clusters in time. We also propose a bootstrap-based approach and two statistical tests that are applied for the purpose of the comparison of models built for different populations of spines. Finally, in Section 4, we present an example of method application. We summarize our work in Section 5.

In this section, we describe the statistical analyses of mammalian dendritic cell populations ACTIVE and CONTROL. A comparison of the descriptor distributions showed that an initial data preprocessing was necessary, which we performed by carefully choosing subsets of the spines from both populations². Finally, we introduce the algorithm for spine representation dimensionality reduction.

In this section, we present PCA and apply two clustering methods to construct the spine shape taxonomy in an unsupervised way. Further, we build the probabilistic model of shape changes in time. Finally, the bootstrap analysis is presented to statistically evaluate the differences between both the resting and potentiated populations.

Here, we provide an example of an application of the method. To obtain the taxonomy of spine shapes, we applied cmeans and hierarchical clustering to DESCRIPTORS^PCA for ACTIVE ∪ CONTROL for t₀ and t₁. To select the proper values of the parameters, we used WSS plots with “knee” shapes (see Figure S3). We obtained k = 10 for hierarchical and k = 8, m = 4 for cmeans clustering. According to WSS measures, these values ensure a good balance between the complexity of the results; i.e., the number of clusters and quality of cluster fitness.

Figure 2A presents the results of hierarchical clustering calculated for ACTIVE ∪ CONTROL according to the procedure described in Section 3.2. Each spine is represented with a single point, and the colors represent the cluster memberships. For each cluster, we identified three representative spines lying nearest to the cluster center. Representative spines are shown in Figure 2B. The obtained clusters express the universal taxonomy of shapes that will be later employed for the computation of Shape Transition Model for ACTIVE, CONTROL, ACTIVE300, and CONTROL300. Representative spines can be used for visual inspection and biological interpretation.

Apart from hierarchical clustering, we also consider cmeans clustering. Table 3 presents the comparison of the prediction error E for both methods. Values were obtained using 10-fold cross-validation. For the purposes of cross-validation, the dataset was randomly separated into ten (10) subsets. Numbers from the same column but in different rows should not be compared. Different clustering methods result in different shape clusters that have different members and thus are incomparable. Although the errors E for different methods have different ranges and cannot be compared, different models with the same method can be compared.

The Shape Transition Model is compared with three baselines. The first baseline is the majority vote model, where all spines from a particular cluster move to a single destination cluster that is selected as the most popular choice. The second baseline is the model, where we assume that all spines remain in the initial clusters, i.e., weights in t₁ are the same as in t₀. Finally, the third baseline assumes random values for the probability P. For both clustering methods, the Shape Transition Model has the smallest error E and predicts spine behavior the best.

Transition graphs of the Shape Transition Model for ACTIVE and CONTROL are shown in Figures 3A,B. Each cluster of shapes is represented by an oval. Initial sizes, i.e., weights of clusters (for hierarchical clustering equivalent to number of spines), are listed. Edges representing transitions are labeled with probability P. They are filtered out, and only transitions (probabilities) greater than 20% of the initial weight are visible.

Only five clusters (numbers 1–5) are well represented in the data. Clusters 1, 2, and 4 are the most dense. Clusters 3 and 5 are interpreted as peripheral. Finally, clusters 6–10 have only a few spines. For transitions from clusters 6–10, high errors were obtained. For example, SE for P(C₉ → C₁₀|C₉) is equal to 66%. Conclusions concerning clusters 6–10 are not reliable. Analogous plots of clustering results and transition graphs for cmeans are presented in Figures S4, S5A,B.

Graphs presented in Figures 3A,B should not be compared because they are computed for populations of different characteristic at t₀. Alternatively, Figure 4 presents a comparison of the transition graphs for CONTROL300 and ACTIVE300 for hierarchical clustering (the exact values of the probabilities can be found in Table S1). A similar analysis for cmeans is presented in Figure S6, and the values of the transitions in percent can be found in Table S2.

In the case of the CONTROL300 and ACTIVE300 subsets (Figure 4), only clusters 1, 2, and 4 contain enough spines to produce credible conclusions. For CONTROL300, cluster 1 has a slightly stronger inertia than for ACTIVE300 (91 vs. 87% spines remained in the same cluster). For cluster 2, the situation is the opposite: 41% of spines from cluster 2 for CONTROL300 remain in cluster 2 compared with 67% for ACTIVE300. For both populations, a large transition of spines from cluster 2 to cluster 1 is observable. However, for CONTROL300, it is present for 52% of the spines, whereas for ACTIVE300, it is present only for 28%. Another difference is visible for transitions from cluster 4. For CONTROL300, 73% of spines move to cluster 1 and 27% to cluster 2. For ACTIVE300, only 54% of spines move to cluster 1, and the rest move to clusters 2–5. Unfortunately, none of the observed differences is significant when the errors are taken into consideration. Therefore, to identify such differences, the models must be compared as a whole.

Table 4 presents p-values of RDC and SMD statistics (Equations 5 and 6) used for a comparison of models for ACTIVE300 and CONTROL300. Results below 0.05 are marked in bold font. Detailed plots of the statistical distributions using kernel estimation are shown in Figures S7, S8. For hierarchical clustering, only SMD shows a significant difference between ACTIVE300 and CONTROL300. This statistic compares transitions of spines between shapes, which is well captured by hierarchical clustering. The RDC statistics depend only on changes of distributions, and hierarchical clustering enforces that each spine belongs to only one shape cluster at the particular time point, which may noticeably affect the overall distributions. In contrast, the distributions are well captured by cmeans clustering, where each spine is an arbitrary mixture of shapes and RDC shows a significant difference. Different clustering methods are sensitive to different properties of the data. The selection of the right clustering method and appropriate test depends on the characteristic of the data that is of interest to the researcher.

The majority of excitatory synapses in the brain are located on dendritic spines. These highly dynamic and plastic structures undergo constant morphological changes in different physiological and pathological processes (Kasai et al., 2010). The structure of the dendritic spines is tightly correlated with their function and reflects the synapse properties. Synapse strengthening or weakening along with dendritic spine formation and elimination assure correct processing and storage of the incoming information in the neuronal network. This plastic nature of the dendritic spines allows them to undergo activity-dependent structural modifications, which are thought to underlie learning and memory formation. At the cellular level, the most extensively studied aspect of this phenomena is related to dendritic spine enlargement in response to stimulation.

In this study, we show a computational method that permits statistical analysis of the impact of an externally applied stimulation on the dendritic spine structural dynamics. We applied statistical tests and examined a population consisting of 923 dendritic spines. We used two dissociated neuronal cell cultures and compared the dendritic spine volume and shape changes between two populations at two different states, unstimulated (CONTROL) and LTP-stimulated (ACTIVE), and at two time points (with a 10-min time interval). We preprocessed the datasets and reduced the dendritic spine number to 300 for each analyzed group.

We provided a probabilistic model for dendritic spine population dynamics. First, the resting state model was constructed (Figure 3A). Then, the probabilistic null model for active neurons was built (Figure 3B). We showed that LTP treatment induced transition of filopodia-like spines (cluster 4) into mushroom-shaped spines (cluster 2). For the first time, we provided exact transition probabilities for this morphological transformation (from cluster 4 to cluster 2, the transition probability was found to be 0.27 ± 0.11). Our result show that in ACTIVE300 in comparision to CONTROL300 there is significant group of growing neurons which partially supports the previous studies (Szepesi et al., 2014) that report chemical LTP-induced spine enlargement in dissociated cultures.

Finally, we compared the models for balanced populations (Figure 4). We found differences between active and non-active neurons. Unfortunately, none of the observed differences between the models was significant when particular transitions between shape clusters were considered. Such absence of significant differences can be explained by too small number of spines in analyzed data. Large errors predominated the differences between values in cells of appropriate transition matrices. However, statistically significant differences were detected when models of populations were compared. Different clustering algorithms showed statistically significant differences between the two analyzed groups (ACTIVE300, CONTROL300). Crisp clustering captured the difference in shapes transitions well, whereas fuzzy clustering captured the difference in changes of shape cluster distributions.

All experimental procedures were carried out in accordance with the Ethical Committee on Animal Research of the Nencki Institute, based on the Polish Act on Animal Welfare and other national laws that are in full agreement with EU directive on animal experimentation.

MM and MR performed biological experiments and did manually segmentation of spines from confocal microscope images. TK, MŁ, GB, and DP designed the computational method and performed statistical analysis. GB contributed in the data analysis and visualization. All authors read and approved the final manuscript.

TK and MŁ were partially supported by research fellowships within “Information technologies: research and their interdisciplinary applications” agreement POKL.04.01.01-00-051/10-00. MM and MR were supported by the grants no. 6651/B/P01/2011/40 and UMO-2015/17/B/NZ3/00557 from National Science Centre Poland. DP, GB, TK, and MŁ was supported by the Polish National Science Centre (Grant numbers 2013/09/B/NZ2/00121 and 2014/15/B/ST6/05082) and COST BM1405 and BM1408 EU actions.