The Mixture of Autoregressive Hidden Markov Models of Morphology for Dentritic Spines During Activation Process

1. Introduction

Dendritic spines are short protrusions that harbor excitatory synapses, which are believed to play a major role in neuronal plasticity and integration through their structural reorganization (Yuste and Denk, 1995; Segal, 2002; Xu et al., 2009). Many physiological and pathological phenomena rely on brain plasticity, including learning and memory, epileptogenesis, drug addiction, and postinjury recovery. Their peculiar shape suggests that spines can serve as an autonomous postsynaptic compartment that isolates chemical and electrical signaling. How neuronal activity modifies the morphology of the spine and how these modifications affect synaptic transmission and plasticity are intriguing issues. 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. Indeed, the induction of long-term potentiation (LTP) or long-term depression is associated with the enlargement or shrinkage of the spine, respectively (Bosch and Hayashi, 2012).

Classification of dendritic spines, according to their shape and size, is a popular strategy to evaluate maturation and pathological changes of neurons. Spines are often classified into four morphological classes: filipodia, thin, stubby, and mushroom (Peters and Kaiserman-Abramof, 1970). It must be noted, however, that the existing categorization of dendritic spine shapes does not provide a clear definition of each group. Some researchers prefer to add two additional categories: branched and cup-shaped spines (Harris and Kater, 1994; Hering and Sheng, 2001). An important question is whether the categorization of spines according to their morphology represents a rigid classification of distinct entities or tentative labeling of transient spine states. Live-cell imaging studies clearly show that spines are very dynamic and undergo reversible transformation between thin and mushroom morphologies in a time scale of minutes to hours (Fischer et al., 1998).

Another important question is how the changes of the shapes of the dendritic spines in time correlate with learning and memory process—is there a correlation between changes of the class of the dendritic spines and changes of the enlargement or shrinkage of the spine and changes with the synaptic strength. Another very important question is how to build a model based on the population of dendritic spines including their classification changes, synaptic strength, learning, and memory process. These questions can help us find answers to better understand the process of learning and memory and how the synaptic plasticity process works. Several neuromorphological tracing tools are available to segment and classify spines in two-dimensional (2D) and/or three-dimensional (3D). The major problem is determining which descriptors in 2D or 3D give the most accurate classification (by accurate classification we mean stable classification of dendritic spines, highest prediction accuracy, and highest biological relevance). Because of the dynamic changes of dendritic spine shapes (and, therefore, its class) in time, we should use time-dependent statistical models to answer for this problem. One of the most important models used in statistics and machine learning is hidden Markov model (HMM), which was described by Rabiner (1989) and Rabiner and Juang (1993). Figure 1A shows a sequence of states (S) and observations (O) for first-order HMM (where sequence of states is S = $S_1,\dots ,S_T$ and sequence of observations is $O=O_1,\dots ,O_T$).

A first-order HMM instantiates two simplifying assumptions. First, as with a first-order Markov chain, the probability of a particular state depends only on the previous state:

Second, the probability of an output observation O_i depends only on the state that produced the observation S_i and not on any other states or any other observations:

Unfortunately, HMM can be poor at capturing dependency between observations because of the statistical assumptions it makes. In many situations, the assumption that observation has a common effect in another observation in the future simplifies the design of the Bayesian network: directed arcs should flow forward in time. The extensions of the HMM called autoregressive HMM (ARHMM) and mixture ARHMM (MARHMM) (Fig. 1B, C) are used especially for time series data (Rezaei Tabar et al., 2019). The ARHMM can be calculated as the forward-directed ARHMM encourages correlation among observations by adding direct dependencies between them. By adding direct dependencies between observations, samples drawn from forward-directed ARHMM are thus smoother than samples from an HMM, usually making a better generative model in time series problems (Stanculescu et al., 2014; Rezaei Tabar et al., 2019). The forward-directed ARHMM is a Bayesian network and obeys the following two conditional independence relations:

Using these conditional independence relations, the joint distribution of a sequence of states and observations can be included as follows:

The backward-directed ARHMM is also known as a Bayesian network. According to D-separation concept in Bayesian network, the backward-directed ARHMM also has the following two conditional independence relations:

The same way, forward-directed ARHMM also has the following two conditional independence relations:

Also we can define the forward $\left(f_k\left(t\right)\right)$ and backward $\left(b_k\left(t\right)\right)$ algorithms that are to find out a recursive way to represent the variable sequence in both models. The forward probability represents the probability of the observation sequence up to time t and the state k at time t, given model $\lambda 1$ (or $\lambda 2$) as the following formulas:

For some stochastic processes, it may be more appropriate to use an ARHMM directed backward in time (Rezaei Tabar et al., 2019). The backward probability represents the probability of the partial observation sequence from $t+1$ to the end, given state k at time t as follows:

In Equations (7)–(10), $\lambda 1$ and $\lambda 2$ are set of parameters where $\lambda$ = (the transition matrix, the initial emission matrix, the emission matrix, and vector of initial state).

Note that the index L − R is timestamp from left to right for representing the forward-directed ARHMM and R − L is timestamp from right to left for representing the backward-directed ARHMM.

For MARHMM, equation describing the model is as follows:

where $P_1\left(O|{\lambda }_1\right)$ and $P_2\left(O|{\lambda }_2\right)$ are the probability of the observation sequences given the forward-directed ARHMM and backward-directed ARHMM, respectively, and ${\alpha }_1$ and ${\alpha }_2$ are mixing weights such that ${\alpha }_1,{\alpha }_2\ge 0$ and ${\alpha }_1+{\alpha }_2=1$. $P\left(O|\lambda \right)$ can be computed by obtained the Baum–Welch algorithm.

The ARHMM and MARHMM parameters are estimated by the expectation-maximization (EM) algorithm (Rezaei Tabar et al., 2019). The algorithm arises in many computational biology applications that involve probabilistic models, and it also enables parameter estimation in probabilistic models with incomplete data. The EM algorithm computes probabilities for each possible completion of the missing data, using the current parameters. These probabilities are used to create a weighted training set consisting of all possible completions of the data. Finally, a modified version of maximum likelihood estimation that deals with weighted training examples provides new parameter estimates (McLachlan and Krishnan, 2007). By using weighted training examples rather than choosing the single best completion, the EM algorithm accounts for the confidence of the model in each completion of the data (Chuong and Batzoglou, 2008). In the rest of the article, we refer to a population of spines stimulated by LTP as dynamic data.

2. Data Preparation

Details about data preparation are described in Bokota et al. (2016) in section “Data preparation and analysis.” Here we describe only the differences in chemical substances that were used in experiments in dynamic data sets gathered at three time points. In this experiment, the chemical LTP was induced through the application of bath in a mixture of 50 μM forskolin, 50 μM picrotoxin, and 0.1 μM rolipram [each dissolved in dimethyl sulfoxide (DMSO)] in a maintenance medium. Subsequently, experiment was obtained from three points in time (timestamps): control (timestamp t₀), 10 minutes after stimulation (timestamp $t_{10}$), and 40 minutes after stimulation (timestamp $t_{40}$).

3. Methods

Details about all four descriptors that were used to classify dendritic spines are given in Supplementary Data S1. In Supplementary Data S2, details about how we prepare the data and how the vectors and matrices should look like to run the code are described. Classification of the dendritic spines in all three timestamps was done by four descriptors. Two descriptors were used for 2D classification and two were used for 3D classification. A Viterbi algorithm that was used in machine learning is also described in Supplementary Data (Supplementary Data S1).

In Matlab (The MathWorks, 2018), we implemented the ARHMM and MARHMM algorithms, whereas Viterbi analysis was conducted in R in version 3.4 in HMM package (Development, 2010). We took dendritic spines marked by SpineMagick descriptor (Ruszczycki et al., 2012) and marked the same dendritic spines using 2D programs (2dSpAn and Spinetools) and 3D programs (3dSpAn and Neurolucida 360). Depending on the descriptor we were working on, dendritic spines were classified as stubby, filipodia, mushroom, thin, and long (combining filipodia and thin into one class) and “not existing” (i.e., after stimulation we observed that these spines did not exist anymore—this class represents spines that were absorbed by the neuron). Finally, we proceed to run the descriptor using the following parameters: data, level, transition probability matrix, initial probability vector, number of iterations, and dependence matrix As an output, the algorithm produces an estimated transition probability matrix, estimated emission probability matrix, and logarithm of probability (logarithm of likelihood) at each iteration. Also, we use methods such as principal component analysis (Jolliffe, 2002), fuzzy partition coefficient, and hierarchical clustering by k-means the same way as by Bokota et al. (2016) and Urban et al. (2019)—the results and comments are described in Supplementary Data S3.

4. Results

In this section, we describe the sensitivity analysis results of MARHMM and machine learning for various combinations of the parameters: initial probability vector, transition probability matrix, and dependence probability matrix, different $\alpha$ values. We also show how the changes in synaptic strength over time depend on the changes in the area surface (for 2D descriptors) or volume (for 3D descriptors). Also, we show that differences in transition probabilities between states in time depend on 2D or 3D classification.

5. Conclusion

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 assures 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 phenomenon is related to dendritic spine enlargement in response to stimulation (Bokota et al., 2016). The presented methods can be used not only in dendritic spines for neurobiological experiments but also in other biological branches [genetics (Kröger et al., 2017), molecular biology (Qing, 2017), immunology (Hill et al., 2017), hematology (Efficace et al., 2017)], or even in other science problems such as in medicine (diagnostic, epidemiologic) (Porcelli and Guidi, 2015; DasMahapatra et al., 2017), psychology (psychology diagnostic) (Porcelli and Guidi, 2015), or sociology (relationships between people) (Wong et al, 2017). In each instance in which we want to know about the next step, we can use the forward-directed ARHMM, if we want to know about the previous step, we can use backward-directed ARHMM, and finally the MARHMM can be used for considering both previous and next steps. Also results from MARHMM used in machine learning can be used in developing a new spice classification descriptor. Also based on results from MARHMM, we can choose the best descriptor, and by using other tools we try to build probabilistic models.

Supplementary Material

Ethics Statement

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 the EU directive on animal experimentation.

Authors' Contributions

P.U. and V.R.T. conceived the project and performed the analysis based on data from experiments. P.U., M.D., G.B., and V.R.T. prepared the figures and wrote the article. N.D. aided in 3D segmentation of dendritic spines using the 3DSpAn descriptor. D.P. and S.B. were supervising the project and conducted extensive discussions during our analysis. All authors read and approved the final article.

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

This study was supported by the Polish National Science Center (Grant No. 2014/15/B/ST6/05082), Foundation for Polish Science (TEAM to D.P.), and by the grant from the Department of Science and Technology, India, under Indo-Polish/Polish-Indo project no.: DST/INT/POL/P-36/2016. The study was cosupported by Grant No. 1U54DK107967-01 “Nucleome Positioning System for Spatiotemporal Genome Organization and Regulation” within 4DNucleome NIH program. S.B. was funded by Department of Biotechnology grant (Grant No. BT/PR16356/BID/7/596/2016).

Supplementary Material

Supplementary Data S1

Supplementary Data S2

Supplementary Data S3

Supplementary Data S4