Morphological analysis of dendrites and spines by hybridization of ridge detection with twin support vector machine

Image acquisition

The cortical neurons used for imaging initiated from the embryonic (E) 18th day rat. Neurons were cultured in vitro until day 22. Finally, Green fluorescent protein (GFP) was transfected using Lipofectamine2000 by Lipofectamine 2000 and neuron was imaged on day 24 by Leica SP5 confocal laser scanning microscopy (CLSM) with a 63X magnification. The size of the image is 1,024 × 1,024 with a resolution of 0.24 µm/pixel at the confocal layer. We added a white line in the figures to indicate the scale bar of 10 micron by 41 pixels.

As the images were captured as a Z-stack series a 3D image stack was mapped onto the ZX, YZ, and XY planes, respectively. However, the slices yielded limited information along the optical (Z) and high computing burden. Therefore, only the projection onto the XY plane was considered.

Image preprocessing

According to the imaging technique (Van der Walt et al., 2014), the photomultiplier tubes (Ge et al., 2016) mainly introduced salt and pepper noise. Therefore, the 2D median filter was used to reduce noise and a partial differential equation (PDE) was used to enhance the image. Furthermore, in order to reduce the disturbance of the backbone extraction we used the opening top-hat and majority operator to eliminate small jiggers.

The original gray image I(x, y) and the structure elements (SE) S(i, j) were then set. The opening top-hat operator (Zhang et al., 2014) is defined as following: (1)${\displaystyle I\circ S=\left(I\mathrm{\Theta }S\right)\oplus S}$ in which I∘S is defined as function I executed the opening operation via SE ⊕ and Θ is defined as the dilation and erosion respectively.

The majority operator is defined as: (2)${\displaystyle P=\left\{\begin{array}{c}1,\text{more than}n\text{positive pixels in its 3 by 3 window}\hfill \\ 0,\text{otherwise}\hfill \end{array}\right.}$ indicating that a pixel is considered to be part of the major line if there are more than n positive pixels in its 3*3 window. Otherwise, the pixel’s noise information must be relegated to 0.

Ridge based backbone extraction

Along the pipeline of spine analysis the backbone detection is a critical step, underlining the importance of finding an effective and robust backbone detection algorithm. In this section we tried to use curvilinear structure detectors to find the dendritic backbone. One of the most well-known curvilinear structure detection methods is to use the principal curves proposed by Kegl & Krzyzak (2002) and Kegl et al. (2000), in which the curve is defined as a line passing through the middle of the data in a certain sense. However, this method was restrictive in that it assumed there would be single lines without self-intersection. Considering that the dendrite has multiple branches, we followed the detector that extracts the curvilinear structure from the data of ridge curves by the density estimation—so-called the ridge-based curvilinear structure detector. The multiple dendritic branches were detected on the assumption that the branches were separated by low density areas. The critical step then becomes the density estimation from the original image data. Here, we employed the nonparametric estimation directly from the samples represented by the observed variable X within a specific compact domain ΩϵR_d . First, the image was normalized into binary image. We set T = 1 for the foreground clutter and T = 0 for the background clutter.

The probabilities (3)${\displaystyle D\left(T=1\right)=\eta ,\text{and}D\left(T=0\right)=1-\eta }$ in which, η ∈ [0, 1]. We defined the generating function ${\left\{f_i\right\}}_{i=1}^n$, $i=\left\{1,2,3,3\dots n\right\}$ with unit length of parameters, where n represents the number of the dendrite. ${\left\{f_i\right\}}_{i=1}^n$ depends on two parameters i ∈ {1, 2, 3, ...., n}, θ ∈ D_i where D_i is represented by the data samples along each dendrite. The output of X depends on function f_i and noise ε : N_d(0, σ²). Density (4)${\displaystyle P\left(I=i\right)=w_i}$ in which w_i > 0, ∑w_i = 1. (5)${\displaystyle D_x\left(x|T=1,I=i,\mathrm{\Theta }=\theta \right)=\frac{1}{{\left(\sqrt{2\pi }\sigma \right)}^d}exp\left(-\frac{{∥x-f_i\left(\theta \right)∥}^2}{2{\sigma }^2}\right).}$

We employ marginalized and the joint density to determine total density, which depends solely on the observed sample point X. (6)${\displaystyle D_{x,T,I,}\left(x,T,i,\theta \right)=D_x\left(x|T=1,I=i,\mathrm{\Theta }=\theta \right)D_{T,I,\mathrm{\Theta }}\left(1,i,\theta \right)=D_x\left(x|T=1,I=i,\mathrm{\Theta }=\theta \right)D_{\mathrm{\Theta }}\left(\theta |I=i\right)D_{T,I}\left(1,i\right)=D_x\left(x|T=1,I=i,\mathrm{\Theta }=\theta \right)D_{\mathrm{\Theta }}\left(\theta |I=i\right)D\left(I=i\right)D\left(T=1\right)}$ and (7)${\displaystyle D_{x,T,I,}\left(x,0,i,\theta \right)=D_x\left(x|T=0\right)D\left(T=0\right)}$ in which, D_x,T,I,(x, 0, i, θ), and Θ represent the joint density and axis along the specified dendrite, respectively. When we summed the joint density over the domains of the discrete random variable I and T (subsequently integrated over the domain of the continuous variable Θ), marginal density (Alghalith, 2016) that provides the probabilities of various values of the variables in the subset without reference to the values of the other variables could be determined as: (8)${\displaystyle d_x\left(x\right)=\frac{\rho }{{\left(\sqrt{2\pi }\sigma \right)}^d}\sum _{i=1}^n{w}_i{\int }_{D_i}exp\left(-\frac{{∥x-f_i\left(\theta \right)∥}^2}{2{\sigma }^2}\right)D_{\mathrm{\Theta }}\left(\theta |I=i\right)d\theta +\frac{1-\eta }{V\left(\mathrm{\Omega }\right)}}$ in which V(Ω) is the volume of domain Ω.

The kernel density estimation was then used, which could be determined from a set of samples $Y={\left\{y_i\right\}}_{i=1}^N\subset R^2$ of the probability density D: (9)${\displaystyle \hat{D}=\frac{1}{N}\sum _{i=1}^N{K}_H\left(x-y_i\right)}$ in which K_H is the Gaussian function (Cecchin et al., 2015) with a symmetric and positive kernel bandwidth in the domain R² → [0, ∞], and H ∈ R^2×2.

We assumed that the ridge was present in the area of high probability density. Thus, the dendrites were thought to be separated by a low density area. The ridge points (Micheal, Vani & Sanjeevi, 2014) were determined based on the principal that modes lie on a ridge curve. Accordingly, we first found modes of $\hat{D}$, and constructed the ridge point sets passing through these modes by using them as initial points. The algorithm is listed in the following steps:

Step 1: The modes of $\hat{D}$ are obtained by establishing the local maximum from each sample point.

Step 2: Step1 is iteratively applied to find the ridge curve components.

Step 3: The ridge curve segments are traced based on the following Eq. (10) adapted to the initial value as (10)${\displaystyle \frac{\mathrm{d}}{\mathrm{d}\theta }\left(P\left(x\left(\theta \right)\right){\nabla }^{2}D\left(x\left(\theta \right)\right)\frac{\nabla D\left(x\left(\theta \right)\right)}{\parallel \nabla D\left(x\left(\theta \right)\right)\parallel }\right)=0,x\left(0\right)=x_0.}$

Step 4: A given point is projected onto the n-dimensional ridge set of $\hat{D}$, which is the density estimate.

Step 5: The ridge curve sets are separated by a low density area with a threshold of the density value ϑ.

The spurs caused by the connected spines were excluded based on the observation that there was only one exact main ridge curve structure passing through an intersection point. The proposed method spitted any other ridge curve components passing through such a point.

Boundary location based on the density invariance

The local line direction was determined via the 2nd directional derivatives of the image (Zamani et al., 2012). Along the two directions perpendicular to the local line, the sharp difference of the density is thought to be the boundary of the dendrite based on the gray scale image.

The density of each pixel is given I(p) of point (x, y) in the original image, and α is the predefined pixel intensity value: (11)${\displaystyle \text{if}\left\{\begin{array}{c}I\left(p\right)\ge \alpha ,p\text{belongs to the line pixel}\hfill \\ I\left(p\right)<\alpha ,p\text{does not belong to the line pixel}.\hfill \end{array}\right.}$

Hu moment invariants based feature extraction

In order to establish biological significance, it is essential to characterize the shape of the spinal structures. Spines are historically classified into three categories based on the shape of specific structures, including the spinal neck and head. Therefore, shape descriptors are critical for spinal classification and quantitative analysis. In this paper we employ a shape descriptor designated as the image moment. For a 2D Neuron image I(x, y), we define the raw moment of order (p + q) as (12)${\displaystyle M_{pq}=\sum _x\sum _y{x}^p{y}^{q}I\left(x,y\right)}$ where p, q = 0, 1, 2, ….

Central moment: In practice, the central moments μ are usually utilized to replace the raw moment in Eq. (13) (13)${\displaystyle {\mu }_{pq}=\sum _x\sum _y{\left(x-\bar{x}\right)}^p{\left(y-\bar{y}\right)}^{q}I\left(x,y\right)}$ (14)${\displaystyle \bar{x}=\frac{M_{10}}{M_{00}},\bar{y}=\frac{M_{01}}{M_{00}}.}$

Central moments are translational-invariant.

Normalized central moment: When dividing the corresponding central moment by the properly scaled (00)th moment, Central moments can be extended to be both scale and translation invariant, The division was called normalized central moment. (15)${\displaystyle {\eta }_{pq}=\frac{{\mu }_{pq}}{{\mu }_{00}^{\left(\frac{p+q}{2}+1\right)}}.}$

Hu moment invariants: In order to allow for invariable rotation the above moments were reformulated. There are two different methods for producing rotation moment invariants described by Hu & Mingkuei (1962). The first method is principal axes. However, when images do not have unique axes that are rotationally symmetrical, the first can be broken down. The second method is Hu moment invariants (HuMI). Hu derived these equations from algebraic invariants and applied them to the moment generating function under the condition of a rotation transformation. They consist of a set of nonlinear centralized moment (NCM) equations as (16)${\displaystyle \begin{array}{ccc}{H}_1\hfill & =\hfill & {\eta }_{20}+{\eta }_{02}\hfill \\ H_2\hfill & =\hfill & {\left({\eta }_{20}-{\eta }_{02}\right)}^2+4{\eta }_{11}^2\hfill \\ H_3\hfill & =\hfill & {\left({\eta }_{30}-3{\eta }_{12}\right)}^2+{\left(3{\eta }_{21}-{\eta }_{03}\right)}^2\hfill \\ H_4\hfill & =\hfill & {\left({\eta }_{30}+{\eta }_{12}\right)}^2+{\left({\eta }_{21}+{\eta }_{03}\right)}^2\hfill \\ H_5\hfill & =\hfill & \left({\eta }_{30}-3{\eta }_{12}\right)\left({\eta }_{30}+{\eta }_{12}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-3{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]\hfill \\ \hfill & \hfill & +\left(3{\eta }_{21}-{\eta }_{03}\right)\left({\eta }_{21}+{\eta }_{03}\right)\left[3{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]\hfill \\ H_6\hfill & =\hfill & \left({\eta }_{20}-{\eta }_{02}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]+4{\eta }_{11}\left({\eta }_{30}+{\eta }_{12}\right)\left({\eta }_{21}+{\eta }_{03}\right)\hfill \\ H_7\hfill & =\hfill & \left(3{\eta }_{21}-{\eta }_{03}\right)\left({\eta }_{30}+{\eta }_{12}\right)\left[{\left({\eta }_{30}+{\eta }_{12}\right)}^2-3{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right]\hfill \\ \hfill & \hfill & -\left({\eta }_{30}-3{\eta }_{12}\right)\left({\eta }_{21}+{\eta }_{03}\right)\left[3{\left({\eta }_{30}+{\eta }_{12}\right)}^2-{\left({\eta }_{21}+{\eta }_{03}\right)}^2\right].\hfill \end{array}}$

It has been clearly observed that HuMI consist of absolute orthogonal (i.e., rotation) moment invariants which can be used for rotation, translation and scale invariant pattern recognition (Noronha & Nayak, 2013; Xiang et al., 2014; Žunic & Žunic, 2013).

Twin SVM based spine classification

We did not use artificial neural network (Hernandez-Serna & Jimenez-Segura, 2014) since it may suffer from overfitting. The principle of the standard SVM involves setting up two parallel planes so that each is nearest to one of two datasets with the two planes being as far apart as possible (Kim et al., 2015; Modinos et al., 2013; Wu & Zhang, 2012; Zhang & Wang, 2015). The planes are described as (17)${\displaystyle w^{T}x+b=0}$ which lie midway between the bounding planes provided by (18)${\displaystyle w^{T}x+b=1\text{and}{w}^{T}x+b=-1.}$

The goal of the standard SVM is to maximize the margin in the form as (19)${\displaystyle \begin{array}{c}\underset{w,b}{Min}\frac{1}{2}{w}^{T}w\hfill \\ \text{subject to}{A}_{i}w\ge 1-b\text{for}{y}_i=1\text{and}{A}_{i}w\ge 1-b\text{for}{y}_i=-1.\hfill \end{array}}$

Jayadeva, Khemchandani & Chandra (2007) proposed a Twin SVM (TSVM) that generated a pair of nonparallel planes for data classification, instead of two parallel planes via the standard SVM. However, in TSVM it was necessary to solve two quadratic programming problems, with one in SVM. Reports indicated that TSVM had better performance with regards to for classification and run time than did standard SVM.

The TSVM classifier is implemented by solving the following two quadratic programming problems (20)${\displaystyle \begin{array}{c}\underset{w^{\left(1\right)},b^{\left(1\right)},q}{Min}\frac{1}{2}{\left(A{w}^{\left(1\right)}+e_{{}^1}{b}^{\left(1\right)}\right)}^T\left(A{w}^{\left(1\right)}+e_1{b}^{\left(1\right)}\right)+c_1{e}_2^{T}q\hfill \\ \text{subject to}-\left(B{w}^{\left(1\right)}+e_2{b}^{\left(1\right)}\right)+q\ge e_2,q\ge 0\hfill \end{array}}$ (21)${\displaystyle \begin{array}{c}\underset{w^{\left(2\right)},b^{\left(2\right)}q}{Min}\frac{1}{2}{\left(B{w}^{\left(2\right)}+e_2{b}^{\left(2\right)}\right)}^T\left(B{w}^{\left(2\right)}+e_2{b}^{\left(2\right)}\right)+c_2{e}_1^{T}q\hfill \\ \text{subject to}-\left(A{w}^{\left(2\right)}+e_1{b}^{\left(2\right)}\right)+q\ge e_1,q\ge 0\hfill \end{array}}$ in which c₁, c₂ > 0 are parameters, and e₁, e₂ are vectors of those with appropriate dimensions.

The TSVM sets up a hyper-planes for each class and classifies a sample according to which hyper-plane any given point is closest to relative to Euclidean distance. TSVM attempts to find the minimal squared distance based on Eq. (21) , which requires keeping the hyper-plane near to the points of one class. Simultaneously, a hyper-plane is required to be at a certain distance from the other class. Should the hyper-plane fall within predefined distances, a set of error variables is used to measure the error in the form of Eq. (19) . In order to improve the classification accuracy we employed the kernel nonlinear classifier Eqs. (22)–(23) to generate surfaces instead of planes as Eqs. (20)–(21). (22)${\displaystyle \begin{array}{c}{Min}_{U^{\left(2\right)},b^{\left(2\right)},q}\frac{1}{2}{∥\left(K\left(B,C^T\right)u^{\left(2\right)}+e_2{b}^{\left(2\right)}\right)∥}^2+c_1{e}_1^{T}q\hfill \\ \text{subject to}\left(K\left(A,C^T\right)u^{\left(2\right)}+e_1{b}^{\left(2\right)}\right)+q\ge e_1,q\ge 0\hfill \end{array}}$ (23)${\displaystyle \begin{array}{c}{Max}_{\gamma }{e}_1^T\gamma -\frac{1}{2}{\gamma }^{T}L{\left(N^{T}N\right)}^{-1}{L}^T\gamma \hfill \\ \text{subject to}0\le \gamma \le c_2\hfill \end{array}}$ in which $L=\left[K\left(A,C^T\right),e_1\right],N=\left[K\left(B,C^T\right),e_2\right]$.

The advantage of the TSVM is that it can be up to four times faster than the traditional SVM depending on computation size. In addition, TSVM has better generation than the traditional SVM.

Backbone extraction and dendrite location

Figure 3 shows different cases of the backbone extraction based on the ridge detection and boundary location via the density difference. The purple color indicates the backbone, and the red color indicates the boundary. Figure 3A show the dendrite location result with multiple branches. Figures 3C and 3D show a single dendrite with one branch. Figures 3E and 3F show a single dendrite with another one branch dendrite. Figures 3G–3J respectively provide one single dendrite. Figure 3 proves that the backbone extraction based on the ridge detection and boundary location based on the density difference is effective.

Spine analysis

Figures 4A–4E illustrate the spine detection results of the sub-images for the condition in which that the dendrite backbone and boundary are marked. The pseudo spines were excluded based on the TSVM. After the dendrite and backbone were detected, our first step was to apply the Hu moment. We extracted the features of the remaining connected components and classified them as either true spines or pseudo-spines. Figure 4 shows two Examples of the pseudo spines, marked by a purple circle.

Table 1 shows the quantities data corresponding to Figs. 4A–4E. As we can find from Table 1, the total length of the dendrite, the number, total length and the area of the spines are recorded for the statistics analysis to reach the biology meaning.

Figure 6 demonstrates that the area and length have similar distributions. For the normalization, L represents for largest value for the smallest area and R for the true value of the spine It is provided as: (24)${\displaystyle N=\frac{R}{L-S}.}$

Finally, The TSVM was utilized to classify the positive spines into predefined categories-MushRoom, Stubby and Thin-using the seven extracted Hu features. Table 2 shows the classification result for the pictures in Fig. 4. The average computation of the classification is 0.0012 s for 20 spines. The classification accuracy can reach 99%. Table 3 shows the Spine classification statistical results, we can find that the “MushRoom” type takes the vast majority of all the spines in our dataset followed by the “Thin:” types.

Comparison of detection rate

We compared our proposed methods with Su’s algorithm and the manual mark. We randomly selected six images which include about 2,000 spines from the dataset and tested the accuracy. Figures 7 and 8 show two image examples to demonstrate the result of spine detection based on different methods, including via manual mark, Su’s method and our proposed RTSVM.

Table 4 shows the number of detected spines based on above three methods. For the first image (ROI_1), the detected number is 33, 25, and 29; for the second (ROI_2) is 28, 22, 23, respectively. It should be noted, however, that the manual method is vulnerable to mistaking as spines pixels belonging to the dendrite. Thus, it will have a much higher false alarm rate than the other methods. Conversely, Su’s method missed some positive spines due to the blurred dendrite boundary.

Comparison of the classification

We employed traditional support vector machine (SVM), generalized eigenvector proximal support vector machine (GESVM) and Twin support vector machine (TSVM) for the spine classification. Table 5 shows the classification result based on different SVMs for ROI_ 1 and ROI_ 2. We can find from Table 5 that GEPSVM and traditional SVM obtain similar classification results and are different from the results of TSVM. TSVM has better performance in the term of classification accuracy when compared to the ground truth (MushRoom 21, Stubby 12, Thin 19), which is provided by five experts in the neuroscience field who manually classified the spines into three predefined categories including “Mushroom,” “Stubby” and “Thin” according to the human vision by using Photoshop. For the conflict of the manual marking, the minority was supposed to subordinate to the major. Meanwhile, the computation time is much faster via TSVM than other two methods, which is 0.0015 s for the 52 spines instead of 0.007 s and 0.0054 s respectively by SVM and GESVM.

Meanwhile, we compared our classification results based on the TSVM and Su’s algorithm. The result is shown in Table 6. We can find that our method obtain better classification result than Su’s algorithm (Su et al., 2014) based on the ground truth.