1. Introduction

Parkinson’s disease (PD) is a common neurodegenerative disease that is characterized by the progressive loss of dopaminergic (DAergic) neurons in substantia nigra pars compacta (SNpc), which results in both motor and cognitive deficits with cardinal symptoms of bradykinesia, rigidity, tremor, and postural instability [1]. Many studies have been done to understand the pathogenesis of PD in experimental cell and animal models [2,3]. However, its complete pathological mechanism remains unknown, thus making it difficult to develop efficient treatments to target the original cause of PD. It has been shown that people who are intoxicated with 1-methyl-4-phenyl-1, 2, 3, 6-tetrahydropyridine (MPTP) develop a syndrome that is nearly identical to PD [4]. Since this discovery, MPTP has become one of the most used neurotoxins that can induce PD-related DAergic neuronal death in mice and non-human primates. An active metabolite of MPTP, 1-Methyl-4-phenyl-pyridium (MPP⁺), is actively transported inside DAergic neurons by dopamine and noradrenaline transporters, resulting in the prevention of oxidative phosphorylation at the level of complex I and an increase in the steady state levels of mitochondrial superoxide [5].

Superoxide is the precursor of most reactive oxygen species (ROS) and a mediator in oxidative chain reactions. ROS lead to severe harmful effects to the cells when they are over-accumulated. This phenomenon is known as oxidative stress. Oxidative stress results in the generation and aggregation of misfolded proteins, which can also lead to excessive ROS and neurotoxicity [6]. In PD, the death of DAergic neurons can be caused by the over-accumulation and posttranslational modification of a-synuclein by oxidative stress [7,8]. This indicates that oxidative stress might be one of the primary mechanisms behind the onset and progression of DAergic neurodegeneration [9]. Therefore, it is important to identify molecular key players, such as genes, proteins, and metabolites, in the downstream processes of oxidative stress for drug target discovery in PD.

A growing number of studies on postmortem PD samples suggest that the activation of common upstream components of PCD pathways may lead to the co-existence of different morphological types of cell death [10]. Co-existence of different morphological types of cell death suggests involvedness of molecular pathways linked to PCD with diverse morphology of cell death [10]. Significant advances have been made in experimental studies of PD, especially through the use of pheochromocytoma (e.g., PC12) and human neuroblastoma (SH-SY5Y, SH-EP) cells. These cells mimic many aspects of the DAergic neuronal death that is observed in PD when treated with neurotoxins such as MPP⁺ and 6-hydroxydopamine (6-OHDA) [11]. Our previous work investigated genome-wide gene expression changes in two types of human neuroblastoma cells (SH-SY5Y [N-type] and SH-EP [S-type]) that were treated with MPP⁺ [12,13]. Neuroblastoma cells have been grouped into three major types depending on morphologies and biochemical properties: neuroblastic (N-type) cells, substrate-adherent (S-type) cells and malignant neuroblastoma stem cells, called I-type (intermediate type) which expresses features of both N- and S-type cells [14]. Both SH-SY5Y and SH-EP cell lines were derivatives of the cell line, human neuroblastoma SK-N-SH.

In this study, parkinsonian neurotoxin MPP⁺ response genes and their interactions in MPP⁺-treated SH-SHSY and SH-EP cells were compared with a linear Gaussian state space model (SSM) [15], which allows one to estimate an internal state vector representing the expression level of gene modules from time-series microarray data. Gene modules might represent groups of co-expressed genes that have similar biological functions from gene expression data of short time courses. There were only 14% common MPP⁺ response genes between N-type SH-SY5Y and S-type SH-EP cells, and gene networks that were estimated from time-series microarray data showed cell-type specificity. This indicates that MPP⁺-induced neuronal cell death might be dependent on cell type. In addition, our results might be helpful for developing potential neuroprotective strategies for PD, as well as understanding the co-existence of apoptotic and nonapoptotic DAergic cell death in postmortem PD patients.

2. Materials and Method

2.1. Microarray data

Time-series microarray data were obtained from our previous reports [12,13] and consisted of eight time points, including control (before MPP⁺ treatment) and 0.5, 1.5, 3, 6, 9, 12, and 24 h after MPP⁺ treatment for SH-SY5Y and SH-EP cells. Microarray data included three replicates and two replicates for each time point in SH-SY5Y and SH-EP cells, respectively. All microarray data were produced with the Illumina bead array (human HT-12 expression v.4 bead array), which includes the probes for 47, 231 genes with well-established or provisional annotation. Raw expression data were log2 transformed and normalized with quantile normalization. After the filtration of probe sets that were redundantly mapped to a single gene and probes that had no corresponding gene name, i.e., gene symbol, the remaining 32, 421 probe sets all had one-to-one mapping to human gene symbols.

2.2. Identification of MPP⁺ response genes

Identification of the MPP⁺ response genes from time-series microarray data was carried out by a software package called EDGE [16]. To use this program, the relative expression level at each time point was required. Thus, the expression value of each time point was divided by that of control and was log2 transformed. For example, let $z_{i}^{j, k}$ be the relative expression level (log2 ratio) of gene i on replicate k, where there are i = 1, 2, …, 32, 421 probe sets and k = 1, 2, 3 for SH-SY5Y cells and k = 1, 2 for SH-EP cells. For each replicate, there were j = 1, …, 7 time points (which correspond to 0.5, 1.5, 3, 6, 9, 12, 24 h after MPP⁺ treatment). The relative expression values are modeled by

$z_{i}^{j, k}={{\mu }_{i}}({{t}^{(j)}})+\xi _{i}^{j, k}$

where ${{\mu }_{i}}({{t}^{(j)}})$ is average (unknown) expression value for gene i evaluated at time ${{t}^{(j)}}$ and the average time curve ${{\mu }_{i}}(t)$ for gene i can be written in terms of a p-dimensional linear basis $[{{s}_{1}}(t), ..., {{s}_{p}}(t)]$ :

${{\mu }_{i}}(t)={{\beta }_{0, i}}+{{\beta }_{i, 1}}{{s}_{1}}(t)+{{\beta }_{i, 2}}{{s}_{2}}+...+{{\beta }_{i, p}}{{s}_{p}}(t)$

where ${{\beta }_{0, i}}$ is the intercept term, p is the same for all genes (the value of dimension p in this study was 3), and $\xi _{i}^{j, k}$ is the random deviation from this curve and is modeled as independent random variables with mean zero and gene dependent $\sigma _{i}^{2}$. To assess differentially expressed genes by MPP⁺ treatment, i.e., MPP⁺ response genes, the goodness-of-model fit under null hypothesis ${{H}_{0, i}}$ that ${{\mu }_{i}}(t)=0$ was compared to that under the alternative hypothesis ${{H}_{1, i}}$ formulated with the general parameterization ${{\mu }_{i}}(t)={{\beta }_{0, i}}+{{\beta }_{i, 1}}{{s}_{1}}(t)+{{\beta }_{i, 2}}{{s}_{2}}+{{\beta }_{i, 3}}{{s}_{3}}(t)$ with some non-zero coefficients, by calculating for gene i a F statistic similar to the one used in ANOVA:

$${{F}_{i}}=\frac{SS_{i}^{0}-SS_{i}^{1}}{SS_{i}^{1}}$$

where $SS_{i}^{0}$ is the sum of squares of the residuals obtained from the null model, and $SS_{i}^{1}$ is the sum of squares of the residuals obtained from the alternative model. The top 1, 000 probes from the q-value rank were selected as MPP⁺ response genes, and the cut-off q-value was 0.01.

2.3. Gene Ontology (GO) analysis

For the estimation of GO terms that are significantly enriched in selected gene list, such as MPP⁺ response genes, a web-based program called GOrilla was employed [17]. This program uses hypergeometric distribution to identify significant GO terms. Given a total number of genes N with B of these genes associated with a particular GO term and n of these genes in the target set, the probability that b or more genes from the target set are associated with the given GO term is given by the hypergeometric tail:

$\Pr \text{ob(}X\ge b)=\sum\limits_{i=b}^{\min (n, B)}{\frac{\left( \begin{matrix} n \\ i \\ \end{matrix} \right)\left( \begin{matrix} N-n \\ B-i \\ \end{matrix} \right)}{\left( \begin{matrix} N \\ B \\ \end{matrix} \right)}}$

The GOrilla program has recognized 22, 466 genes out of 31, 421 gene symbols, and only 18, 002 of these genes are associated with a GO term.

2.4. Identification of gene modules and module gene interactions with the SSM

To find gene modules from MPP⁺ response genes using time-series gene expression data, a linear Gaussian SSM proposed by Hirose et al. [15] was employed. In the SSM, state variables are used to represent unobservable biological entities, such as transcription factors or protein complexes, and temporal gene expression profiles are mapped to the state variables to model the dynamics and dependencies of the state variables. Briefly, let ${{y}_{t}}\in {{R}^{p}}, \text{ }t\in {{T}_{obs}}\subseteq T$ be a series of vectors containing observed expression levels of p genes at the t-th time points. In this study, the value of p and t was 1, 000 and 8, respectively. A sequence of the observation vectors, ${{Y}_{{{T}_{obs}}}}=\left\{ {{y}_{t}} \right\}, \text{ }n\in {{T}_{obs}}$, is modeled by supposing that ${{y}_{t}}$ is generated from the k-dimensional hidden-state variable denoted by ${{x}_{t}}$. In this study, the dimension of state variable was chosen as four (i.e., k = 4), based on the BIC suggested by Yamaguchi et al. [18], for the SSM of SH-SY5Y and SH-EP cells. The SSM is shown as follows.

$$\begin{align} & {{x}_{t}}=F{{x}_{t-1}}+{{v}_{t}}, \text{ }t\in T\text{ (system model)} \\ & {{y}_{t}}=H{{x}_{t}}+{{w}_{t}}, \text{ }t\in {{T}_{obs}}\text{ (observation model)} \\ \end{align}$$

where F is the state transition matrix $(k\times k),$ H is the observation matrix $(p\times k),$ and ${{v}_{t}}\tilde{\ }{{N}_{k}}({{0}_{k}}, Q)\text{ and }{{w}_{t}}\tilde{\ }{{N}_{p}}({{0}_{p}}, R)$ ${{v}_{t}}\tilde{\ }{{N}_{k}}({{0}_{k}}, Q)\text{ and }{{w}_{t}}\tilde{\ }{{N}_{p}}({{0}_{p}}, R)$ are the system noise and the observation noise, respectively. The initial state vector ${{x}_{0}}$ is assumed to be a Gaussian random vector with mean vector ${{\mu }_{0}}$ and covariance matrix ${{\sum }_{0}}$, i.e., ${{x}_{0}}\tilde{\ }{{N}_{k}}({{\mu }_{0}}, {{\sum }_{0}})$. Parameters including H, F, R, Q, ${{\mu }_{0}}$, and the state vector in the model were estimated using the SiGN-SSM software [19] under the following constraints:

$$\begin{align} & \cdot \text{ }Q={{I}_{k}} \\ & \cdot \text{ }{{H}^{'}}{{R}^{-1}}H=\Lambda \equiv \text{diag(}{{\lambda }_{1}}, ..., {{\lambda }_{k}})\text{ where }{{\lambda }_{1}}>{{\lambda }_{2}}>\cdots >{{\lambda }_{k}} \\ & \cdot \text{ An arbitrary sign condition is imposed on the elements of the first row of }H \\ \end{align}$$

After the estimation of parameters, gene expression vectors $({{y}_{t}})$ can be mapped onto the state space ${{R}^{k}}$ with the projection matrix $D\text{ (}k\times p\text{)}$ as follows:

$$\begin{align} & {{x}_{t}}=\text{ }D{{R}^{-1/2}}({{y}_{t}}-{{w}_{t}}), \text{ }t\in {{T}_{obs}} \\ & D={{\Lambda }^{-1}}{{H}^{T}}{{R}^{-1/2}} \\ \end{align}$$

The gene module is identified by using the above projection matrix. That is, the element in the i-th row and j-th column of projection matrix $D\text{ }(k\times k)$ represents the effects of the j-th gene on the i-th module and was designated as ${{d}_{ij}}$. If genes have large $\left| {{d}_{ij}} \right|$ for each module (each row of matrix D), then groups of genes, i.e., modules, can be obtained. In this study, four gene modules including 100 genes were obtained by selection of the highest and lowest ranked 50 genes in each row of projection matrix D. As the signs of ${{d}_{ij}}$, which corresponds to the highest and lowest ranked 50 genes in each row of projection matrix D, were positive and negative, respectively, each module could also be divided into two sub-modules (i.e., ${{M}_{i+}}\text{ and }{{M}_{i-}}$).

The interactions or regulation structures at gene levels were obtained by the conversion of the SSM into a parsimonious representation of the first-order vector autoregressive form, as shown below:

$${{R}^{-1/2}}({{y}_{n}}-{{w}_{n}})=\Psi {{R}^{-1/2}}({{y}_{n-1}}-{{w}_{n-1}})+{{R}^{-1/2}}H{{v}_{n}}$$

where the autoregressive coefficient matrix is given by $\Psi \equiv {{D}^{T}}\Lambda FD$. Matrix $\Psi \text{ (}p\times p\text{)}$ represents temporal interactions in gene level in the following manner: once the k modules ${{x}_{n-1}}=D{{R}^{-1/2}}({{y}_{n-1}}-{{w}_{n-1}})$ are given at time n−1, the current module are calculated with the scaled-system matrix $\Lambda F$ as $\Lambda F{{x}_{n-1}}=\Lambda FD{{R}^{-1/2}}({{y}_{n-1}}-{{w}_{n-1}})$.

2.5. Construction of gene regulation networks and visualization

The Interaction matrix ${{\Psi }_{m}}$ between module genes was obtained by extraction of the elements corresponding to module genes in the matrix $\Psi$ estimated in the SSM. As the total number of module genes was 331 in SH-SY5Y and 318 in SH-EP cells, the size of ${{\Psi }_{m}}$ was 331×331 and 318×318 for SH-SY5Y and SH-EP cells, respectively. Gene networks were constructed by the selection of elements satisfying $\left| {{\Psi }_{m, ij}} \right|$ > 0.012, where the edge density was nearly constant for both cases of SH-SY5Y and SH-EP cells. The estimated networks were visualized by using Cytoscape v3.2.1 (http://www.cytoscape.org), and communities inside the networks were determined with Cytoscape plug-in ModuLand [20].

3. Results and Discussion

3.1. Comparison of MPP⁺ response genes between SH-SY5Y and SH-EP cells

Time-series microarray data for MPP⁺-treated neuroblastoma SH-SY5Y and SH-EP cells were obtained from previous studies [12,13]. These microarray data consisted of control (time 0) and seven time points (0.5, 1.5, 3, 6, 9, 12, and 24 h) after MPP⁺ treatment, and each time point was repeated three times for SH-SY5Y cells and two times for SH-EP cells. Prior to the comparison of SH-SY5Y and SH-EP cells, differentially expressed genes (DEGs) during the time course were identified with the Extraction of Differential Gene Expression (EDGE) program [16], which introduces the time variable through a gene expression curve expanded over the polynomial or B-spline basis with the coefficients. The most significant curves were selected by q-values using a false discovery ratio (FDR)-like procedure. From the rank by q-value, the top thousands of genes were selected from 32, 421 genes for both types of cells. These selected genes were considered to be MPP⁺ response genes, and the q-value cutoff was 0.01. The number of common MPP⁺ response genes between SH-SY5Y and SH-EP cells was 140, and the remaining 1, 720 genes were cell type specific. This indicates that the molecular response of MPP⁺ toxicity might differ between SH-SY5Y and SH-EP cells. On the other hand, common MPP⁺ response genes between SH-SY5Y and SH-EP cells might be involved in core mechanisms for MPP⁺ toxicity.

To examine geneontology (GO) terms that are significantly overrepresented in the 140 MPP⁺ response genes between SH-SY5Y and SH-EP cells, the web-based program GOrilla [17] was employed. The results are shown in Table 1. Most of the significantly enriched GO biological process (BP) terms showed involvement in the positive regulation of transcription in response to endoplasmic reticulum (ER) stress. This indicates that ER stress might play a key role in MPP⁺-induced neuronal cell death. In addition, the significant enrichment of CHOP (CCAAT/enhancer-binding protein)-ATF4 (activating transcription factor 4) complex in GO cellular component (CC) terms suggests that the expression of CHOP-ATF4-targeted genes might be changed by MPP⁺ treatment.

Most of the common MPP⁺ response genes showed similar expression patterns. However, the gene expression of 13 common MPP⁺ response genes represented clear contrasts at 12 and 24 h after MPP⁺ treatment (Table 2). Genes such as ETS2, ZHX2, H1F0, C10orf35, and CNBP showed down-regulation and up-regulation in SH-SY5Y and SH-EP cells, respectively. In contrast, the opposite regulation in both types of cells was observed for CTDSPL, ABCB6, DACT3, P704P, TMBIM4, and ODZ4. In SH-SY5Y cells, down-regulation of H1F0 might induce defects in nucleosome structure, whereas up-regulation of TMBIM4 may inhibit MPP⁺ toxicity-induced apoptosis. In SH-EP cells, down-regulation of ABCB6 may decrease the adenosine triphosphate-dependent uptake of toxic molecules, such as MPP⁺, into mitochondria, whereas up-regulation of transcription factors, such as ETS2, ZHX2, and CNBP, might contribute to changes in the expression of other MPP⁺ response genes. The existence of the differential regulation of common MPP⁺ response genes indicates that cellular downstream mechanisms toward the same neurotoxin (i.e., MPP⁺) might be cell type dependent. Both SH-SY5Y and SH-EP cells are derived via sub-cloning of SK-N-SH cells, but they are morphologically different. SH-SY5Y is an N-type cell line that expresses neurotransmitters and various neuronal cell-surface markers, whereas SH-EP is an epithelial substrate-adherent Schwannian (S-type) cell line that expresses proteins that are characteristic of Schwann cells and lack neuronal markers [14]. The physiological and genomic differences between these two types of cells might be responsible for the differential regulation of the same genes toward MPP⁺ toxicity.

3.2. Identification of gene modules in MPP⁺-treated SH-SY5Y and SH-EP cells with SSM

Modularity is recognized as a design principle in biological systems, as it has been observed in protein-protein interactions, metabolic networks, and transcriptional regulation networks [21]. The identification of gene modules involved in MPP⁺-induced neuronal cell death might open a way to systems approach to PD. Here, the time-series microarray data of MPP⁺-treated SH-SY5Y and SH-EP cells were used to find gene modules, which represent groups of transcriptionally co-expressed genes with similar biological functions. In this study, a linear Gaussian SSM, which was suggested by Hirose et al. [15], was employed for the identification of gene modules using time-dependent gene expression values of MPP⁺ response genes in SH-SY5Y and SH-EP cells. In the SSM, a gene module is considered as a set of genes that are regulated in concert by a shared program that governs their behavior [15]. To find gene modules for MPP⁺-treated SH-SY5Y and SH-EP cells, a sequence of 1, 000-dimensional observation vectors $\{{{y}_{1}}, {{y}_{2}}, ..., {{y}_{8}}\}$ (i.e., gene expression values at 8 time points for 1, 000 MPP⁺ response genes) was modeled by assuming that at each time step, ${{y}_{t}}$ was generated from a k-dimensional hidden-state variable that was denoted by ${{x}_{t}}$. Furthermore, the sequence $\{{{x}_{1}}, {{x}_{2}}, ..., {{x}_{8}}\}$ defined a first-order Markov process. The dimension of state vector k was determined was chosen by minimizing Bayesian information criterion (BIC). Here, k = 4 was used for SH-SY5Y and SH-EP cells. The SSM can be described by the two following equations:

$\begin{align} & {{x}_{t}}=F{{x}_{t-1}}+{{v}_{t}}\text{, }t\in \{1, 2, ..., 8\}\text{ (system model)} \\ & {{y}_{t}}=H{{x}_{t}}+{{w}_{t}}\text{, }t\in \{1, 2, ..., 8\}\text{ (observation model)} \\ \end{align}$

where F is the state transition matrix (4×4), H is the observation matrix (1000×4), ${{v}_{t}}\sim {{N}_{4}}({{0}_{4}}, Q)$ is the system noise, and ${{w}_{t}}\sim {{N}_{1000}}({{0}_{1000}}, R)$ is the observation noise. The initial state vector ${{x}_{0}}$ is assumed to be a Gaussian random vector with mean vector ${{\mu }_{0}}$ and covariance matrix $\sum\nolimits_{0}{, }\text{ }i.e., {{x}_{0}}\sim {{N}_{k}}({{\mu }_{0}}, \sum\nolimits_{0}{)}$. Parameters, including H, F, R, Q, ${{\mu }_{0}}$, and the state vector ${{x}_{t}}$, were estimated using maximum likelihood estimation with the EM algorithm [15,19]. As the number of time points was too small compared with that of MPP⁺ response genes, all replicates of time-series expression data were incorporated into the parameter estimation. All the parameters in the SSM were estimated with the SiGN-SMM software [19], which estimates a SSM for analyzing short-time and/or replicated time-series gene expre ssion profiles.

During the parameter estimation process, the state vectors were constructed such that they were likely to follow the first-order Markov process:

${{x}_{t}}=D{{R}^{-1/2}}({{y}_{t}}-{{w}_{t}}), \text{ }t\in \{1, 2, ..., 8\}$

where D is the projection matrix (4×1000) and defined by ${{(H'{{R}^{-1}}H)}^{-1}}H'{{R}^{-1/2}}$. The element in the i-th row and j-th column of matrix D represents effects of the j-th gene on the i-th module and was designated as ${{d}_{ij}}$. Gene modules, which are relevant to temporal gene expression patterns, were identified by selecting genes that have large $\left| {{d}_{ij}} \right|$ in each row of matrix D. That is, the genes at each module were selected according to the following method: the j-th gene was associated to the i-th positive module ${{M}_{i+}}$ or the i-th negative module ${{M}_{i-}}$ if the element ${{d}_{ij}}$ in the i-th row of matrix D was ranked in the highest or lowest 50, respectively. As the number of rows in matrix D is “4”, 4 modules were identified, and each module included 2 sub-modules. The expression profiles of module genes in SH-SY5Y and SH-EP cells are shown in Figure 1.

The genes located at the same sub-modules represented very similar expression patterns as well as good reproducibility for replicates. The aggregation of genes with similar expression patterns indicates the high relevance of the genes in each module. Table 3 shows over-represented GO BP terms in each module of SH-SY5Y and SH-EP cells. The relevance of module genes for the GO BP was more apparent in SH-EP cells than in SH-SY5Y cells. For example, in SH-EP cells, the genes in represented enrichment of GO terms related to nucleic acid transport, whereas those in showed over-representation of GO terms related to cell cycle. This indicates that genes belonging to the same module might be involved in similar biological processes. Also, each module might play a role in the cellular response to MPP⁺ treatment. Table S1 shows gene information for each module.

The estimated matrix F (4×4) included the temporal relationship between modules, as shown in the system model. That is, the element in the i-th row and j-th column of matrix F represented the effects of the j-th module on the i-th module. These interactions between modules are shown in Figure 2. In both SH-SY5Y and SH-EP cells, all modules showed a self-loop edge with positive regulation, and their strengths were stronger than those among other modules. Positive auto-regulation in gene regulation networks is often mediated through signal molecules and may cause bi-stability [22]. Thus, the positive auto-regulation of modules that include cellular damage-related genes might result in cell death. Weak regulation among modules indicates that genes in each module might be barely affected by genes that belong to the other module. On the other hand, strong self-regulation of each module suggests that genes in each module might be highly associated.

3.3. Estimation of module genes’ regulation structure in MPP⁺-treated SH-SH5Y and SH-EP cells from the SSM

The SSM that was employed for the estimation of gene modules can be converted to the first-order autoregressive form [15], as shown below:

$\begin{align} & {{R}^{-1/2}}({{y}_{t}}-{{w}_{t}})=\Psi {{R}^{-1/2}}({{y}_{t-1}}-{{w}_{t-1}})+{{R}^{-1/2}}{{H}_{t}}{{v}_{t}} \\ & \Psi ={{D}^{T}}\Lambda FD, \text{ }\Lambda \text{=}{{H}^{T}}{{R}^{-1}}H \\ \end{align}$

where all parameter matrixes, including $R(p\times p), \text{ }D(k\times p), \text{ }\Lambda (k\times k), \text{ }F(k\times k), \text{ and }H(p\times k)$ (here, p—the number of MPP⁺ response genes, and k—dimension of state variable, are 1000 and 4, respectively), have already been calculated during the gene module estimation process with the SSM in the previous section. The coefficient matrix $\Psi (1000\times 1000)$ represents the temporal structure of gene-level networks in the following way: once the k modules ${{x}_{t-1}}=D{{R}^{-1/2}}({{y}_{t-1}}-{{w}_{t-1}})$ are given at time t-1, the current modules are generated with the scaled-system matrix $\Lambda F$ as $\Lambda F{{x}_{t-1}}=\Lambda FD{{R}^{-1/2}}({{y}_{t-1}}-{{w}_{t-1}})$. In addition, the current modules regulate expression values of p genes in ${{y}_{t}}$ with the reconstruction matrix ${{D}^{T}}$ as ${{D}^{T}}\Lambda F{{x}_{t-1}}=\Psi {{R}^{-1/2}}$. Thus, matrix D, F and Ψ in the SSM correspond to transcriptional modulation, module-module interactions, and gene-gene interactions, respectively. The element in the i-th row and j-th column of matrix Ψ represents the effects of the j-th gene on the i-th gene and was designated as ${{\Psi }_{ij}}$. The absolute value of ${{\Psi }_{ij}}$ represents the regulation strength, and the sign of ${{\Psi }_{ij}}$ indicates activation (+, positive regulation) or inhibition (−, negative regulation). The regulation structure between module genes might be important for understanding the underlying cellular mechanism of MPP⁺-induced neuronal cell death. Thus, instead of all MPP⁺ response genes, we focused on interactions between module genes that were identified in the previous section. In SH-SY5Y and SH-EP cells, the total number of modules was 4, and each module included 100 MPP⁺ response genes. However, the total number of module genes was 331 and 318 in SH-SY5Y and SH-EP cells, respectively. This was due to the existence of genes belonging to multiple modules. To investigate interactions between module genes, the module gene interaction matrix ${{\Psi }_{m}}$ was constructed by selection of ${{\Psi }_{ij}}$ corresponding to module genes.

Thus, the size of ${{\Psi }_{m}}$ was 331×331 and 318×318 for SH-SY5Y and SH-EP cells, respectively. The topology of gene interaction or regulation network is generally sparse, i.e., a small constant number of edges per node, much smaller than the total number of nodes. To construct sparse networks from the module gene interaction matrix ${{\Psi }_{m}}$, a threshold for interaction value can be used. Because an almost constant average edge density of estimated networks from the module gene interaction matrix ${{\Psi }_{m}}$ was observed at near a threshold of 0.012 for both SH-SY5Y and SH-EP cells, 0.012 was used as the threshold. That is, the only elements satisfying $\left| {{\Psi }_{m}} \right|>0.012$ were selected for the construction of the sparse gene network. Estimated gene networks with a threshold of 0.012 for SH-SY5Y cells included 249 edges and 87 nodes, and had an average edge density of 2.86 (Figure 3). The color of edge in the figure represents regulation type (i.e., red for activation and blue for inhibition), and its width indicates regulation strength. Three communities were identified with Cytoscape plug-in ModuLand [20] from the networks of SH-SY5Y cells, and the nodes (genes) in the same community are shown with the same color. Genes including DCN (decorin) and HIST1H2BK (histone cluster 1) showed high node degrees in the leftmost community, whereas C5orf40 (fibronectin type III domain containing 9) and SGSM1 (small G protein signaling modulator 1) represented high node degrees in the middle and rightmost community, respectively (Figure 3). The relationship of hub genes, such as DCN, HIST1H2BK, and C5orf40, to the cellular structure indicates that MPP⁺ toxicity in SH-SY5Y cells might induce instability of the cellular structure.

The regulation structure between genes with node degrees above four is shown in Figure 4. Strong positive auto-regulation was observed at hub genes, such as DCN, HIST1H2BK, and C5orf40. A zinc finger protein 266 (ZNF266) showed positive association with two hub genes, including DCN and MGC27121, but showed negative association with HIST1H2BK. ZNF266 was located in a cluster of similar genes encoding zinc finger proteins on chromosome 19 and may be involved in gene regulation. At 24 h after MPP⁺ treatment, DCN was upregulated (fold change = 2.4), whereas HIST1H2BK was down-regulated (fold change = 0.30). These regulation patterns showed a good agreement with estimated regulation structure between these genes. SGSM1 showed no association with other genes, as none of the genes that were directly connected to SGSM1 had more than four edges. In the community with purple nodes (Figure 3), all genes except for TMEM79 (transmembrane protein 79) showed negative association with SGSM1. This regulation structure showed a good agreement with the time-dependent expression patterns of genes in the community (Table S2). That is, the expression values of TMEM79 and SGSM1 showed time-dependent increases, whereas those of BCKDK (branched cahin ketoacid dehydrogenase), CDK6 (cyclin dependent kinase 6), ENPP2 (ectonucleotide pyrophosphatase/phosphodiesterase 2), and ICA1 (islet cell autoantigen 1) represented time-dependent decrease. BCKDK is found in the mitochondrion, where it phosphorylates and inactivates BCKD. Thus, down-regulation of BCKDK might be related to mitochondrial stress.

In the case of SH-EP cells with a threshold of 0.012, the estimated networks included 478 edges and 99 nodes and had an average edge density of 4.83 (Figure 5). In these networks, three communities were identified by ModuLand [20], and each community node is shown as a different color. Eighty-eight percent of genes in the networks appeared in the leftmost community. Twenty-three genes had node degrees above 10, and their regulation structure is shown in Figure 6. Except for genes including MSH6 (mutS homolog 6) and RBCK1 (RanBP-type and C3HC4-type zinc finger containing 1), all hub genes came from the same community, including cyan node. The topmost gene MSH6 encodes a member of the DNA mismatch repair MutS family, which helps in the recognition of mismatched nucleotides prior to their mismatch [23]. This gene showed negative association with RBCK1, whose gene product has a function of E3 ubiquitin-protein ligase. E3 ubiquitin-protein ligase accepts ubiquitin from specific E2 ubiquitin-conjugating enzymes and transfers it to its substrates [24]. At 24 h after MPP⁺ treatment, MSH6 and RBCK1 showed down-regulation (fold change = 0.47) and up-regulation (fold change = 2.0), respectively. The negative association of MSH6 and RBCK1 showed a good agreement with their regulation patterns. In addition, the down-regulation of MSH6 may induce DNA damage and increase the number of misfolded proteins, which result in the up-regulation of ubiquitination-related genes, such as RBCK1. This indicates that MPP⁺ might induce the acceleration of DNA damage and ER stress by misfolded proteins in SH-EP cells. Positive auto-regulation was observed at all hub genes except for MTHFD2 (methylenetetrahydrofolate dehydrogenase [NADP⁺-dependent] 2) and ZCCHC8 (zinc finger, CCHC domain containing 8). Positive auto-regulation may contribute to the continuous increase or decrease in gene expression patterns. The fold changes of zinc finger genes, such as ZCCHC8, ZNF772 (zinc finger protein 772) and ZNF26 (zinc finger protein 26), were high at 12 and 24 h after MPP⁺ treatment (Table S2).

This indicates that regulated or targeted genes by these zinc finger proteins might play a role in MPP⁺ toxicity. The hub gene TRIB3 (tribbles pseudokinase 3), which is involved in DDIT3/CHOP-dependent cell death during ER stress, showed very high fold changes at 24 h after MPP⁺ treatment, whereas the expression of hub gene NINJ2 (ninjurin 2), which may have a role in nerve regeneration after nerve injury, showed a time-dependent decrease. It is interesting to note that RNA genes, such as RNU86 (U86 small nucleolar RNA), EPB41L4A (EPB41L4A antisense RNA 1), SNHG1 (small nucleolar RNA gene 1), and C7orf40 (small nucleolar RNA host gene 15), were observed in the hub gene networks, and most of their expression showed time-dependent increases. Two hub genes, including CTH (cystathionine gamma-lyase) and CARS (cysteinyl-tRNA synthetase), showed up-regulation from 9 h after MPP⁺ treatment (Table S2). The cells under oxidative stress can require high amounts of cysteine, because the synthesis of glutathione might depend on its availability, which may be the most important intracellular defense against oxidative stress. Thus, the increasing expression values of CTH and CARS in MPP⁺-treated SH-EP cells might indicate that MPP⁺-triggered oxidative stress might increase with time.

4. Conclusions

Growing evidence suggests that the mitochondrion-dependent programmed cell death (PCD) pathways may lead to different morphological types of cell death in the brain of PD patients [10]. This suggests that PCD pathways that are triggered by the inhibition of complex I might result in cell death via various downstream processes, such as mitochondrion-mediated caspase-dependent apoptosis or caspase-independent non-apoptotic cell death. The factors that affect various downstream processes of mitochondrial complex I inhibition are unknown. Therefore, various downstream processes of mitochondrial complex I inhibition should be considered as potential neuroprotective strategies for PD. Because the inhibition of mitochondrial complex I by MPP⁺ reproduces several PD-linked cellular alterations, such as increased production of ROS, oxidative damage to lipids, DNA, and proteins [25,26,27], MPP⁺-treated cells have been employed as a PD model. We compared MPP⁺ responses in two types of human neuroblastoma cells (SH-SY5Y and SH-EP), which were derived from the common cell line SK-N-SH, but their morphologies and biochemical properties are different. SH-SY5Y cells belong to the N-type which has neuroblast characteristics, such as neurotransmitter biosynthetic enzymatic activities and catecholamine uptake. In contrast, SH-EP cells are similar to immature glial cells (S-type), because they are not neuronal but do produce type I and III collagens [28].

SH-SY5Y and SH-EP cells shared only 14% of their own MPP⁺ response genes. GO analysis for these common MPP⁺ response genes showed significance in ER stress by misfolded proteins. This indicates that ER stress by misfolded proteins might play a key role in MPP⁺-induced neuronal cell death. However, the mechanism that relates complex I inhibition to ER stress by misfolded proteins might be different, as 86% of their own MPP⁺ response genes displayed cell-type specificity. To identify groups of influential genes, i.e., gene modules, during MPP⁺-induced neuronal cell death, the SSM consisting of system and observation model was employed. Time-dependent expression values of MPP⁺ response genes were used for the estimation of SSM parameters, which provides a basis of module and gene interactions. Gene expression patterns within the same gene module were quite similar, which suggests the coordinate regulation of genes in the same module (Figure 1). As module genes have influential effects on the expression of other genes, their interaction structure might provide an insight into MPP⁺-induced neuronal cell death. Interactions in gene levels were estimated using parameters of the SSM, which were calculated during the module identification process. Gene networks for SH-SY5Y and SH-EP cells were constructed based on the gene-level interaction matrix. In the SH-SY5Y networks, a hub gene, HIST1H2BK, was severely down-regulated at 24 h after MPP⁺ treatment (Table S2). This indicates that the MPP⁺-mediated inhibition of mitochondrial complex I might affect nucleosome structure. Nucleosomes are consisting of four major classes of histones (H3, H4, H2A, and H2B) [29] and build the basic units of eukaryotic chromatin. Packaged DNA into chromatin also plays as a key regulator with transcriptions factors [30]. Two hub genes, MSH6 and RBCK1, were down-regulated and up-regulated at 24 h after MPP⁺ treatment (Table S2) in the SH-EP networks, respectively.

This suggests that MPP⁺-mediated inhibition of mitochondrial complex I might induce ER stress by misfolded proteins, which might be produced via DNA mutation. The comparison of two networks indicates that the MPP⁺-mediated inhibition of mitochondrial complex I might induce neuronal cell death via various cellular downstream processes rather than a single determined process. Therefore, nucleosome structure defects and ER stress by misfolded proteins might be key downstream processes for MPP⁺ treated SH-SY5Y and SH-EP cells, respectively. This result might explain why postmortem PD patients exhibit variations in cell death types, such as mitochondrion-mediated caspase-dependent apoptosis and caspase-independent nonapoptotic cell death. In addition, hub genes in SH-SY5Y and SH-EP networks might be considered as potential targets for PD treatment or therapy.

Acknowledgments

This study was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (MEST) (No.2012R1A1A2005622), and was also supported by Dong Yang University in the year 2014.

Conflicts of Interest

The author declares no conflict of interest.

Supplementary

Table S1. Module genes of SH-SY5Y and SH-EP cells, as estimated by the state space model. Table S2. Fold changes of node genes shown in Figures 3 (networks of SH-SY5Y cells) and 5 (networks of SH-EP cells) at each time point after MPP⁺ treatment. Provided as separate files.