Functional connectivity changes during a working memory task in rat via NMF analysis

Working memory (WM) is necessary in higher cognition. The brain as a complex network is formed by interconnections among neurons. Connectivity results in neural dynamics to support cognition. The first aim is to investigate connectivity dynamics in medial prefrontal cortex (mPFC) networks during WM. As brain neural activity is sparse, the second aim is to find the intrinsic connectivity property in a feature space. Using multi-channel electrode recording techniques, spikes were simultaneously obtained from mPFC of rats that performed a Y-maze WM task. Continuous time series converted from spikes were embedded in a low-dimensional space by non-negative matrix factorization (NMF). mPFC network in original space was constructed by measuring connections among neurons. And the same network in NMF space was constructed by computing connectivity values between the extracted NMF components. Causal density (Cd) and global efficiency (E) were estimated to present the network property. The results showed that Cd and E significantly peaked in the interval right before the maze choice point in correct trials. However, the increase did not emerge in error trials. Additionally, Cd and E in two spaces displayed similar trends in correct trials. The difference was that the measures in NMF space were significantly greater than those in original space. Our findings indicated that the anticipatory changes in mPFC networks may have an effect on future WM behavioral choices. Moreover, the NMF analysis achieves a better characterization for a brain network.


INTRODUCTION
Working memory (WM) has been defined as a limited capacity system which temporarily maintains and manipulates information "online" during various cognitive tasks (Baddeley, 2003). It supports human thought processes by providing an interface between perception, long-term memory and action. Consequently, the neural activity supporting WM underlies a broad range of advanced cognitive functions (Jonides et al., 2008).
The Y-maze spontaneous alternation task is one of the most widely used WM paradigms. Baeg et al. (2003) have found that firing rate changes of PFC neurons is tied to the reward and reward expectancy during a spatial delayed alternation task on a figure 8-shaped maze, implicating that WM is mediated by the PFC neuronal activities. Benchenane et al. (2010) have suggested that prefrontal cortex (PFC) is involved in rats learning new rules on a Y-maze both at the local field potential and spike levels. There is also evidence that PFC neuronal activity displays anticipatory changes prior to future goal-directed behaviors, suggesting that the PFC is responsible for holding neural representations for WM (Jo et al., 2013).
During a cognitive process, neuronal discharges are not independent but serially correlated with each other (Constantinidis and Goldman-Rakic, 2002). Therefore, a brain can be functionally defined as a network. Studies of functional network organization have become one of the most popular topics in the field of neuroscience (Palva et al., 2010;Sakaki et al., 2013). Several research has shown that cognitive performance was mediated by interactions among distributed, functionally specialized brain networks (White et al., 2011;Sala-Llonch et al., 2012). Conversely, specific network disruptions can cause a range of neurological disorders (Seeley et al., 2009). A powerful technique called Granger causality analysis could be used to extract such connections. Brain connectivity in humans or animals during cognitions has been analyzed using EEG (Anderson et al., 2010), fMRI (Schlösser et al., 2006), MEG (Bressler and Seth, 2011) and LFPs (Bressler et al., 2007), providing satisfactory explanations for neural mechanisms. These studies have opened new thought that the functional connectivity in mPFC networks may play an essential role in the WM. With the rapid development of chronic multi-electrode implantation, spikes can be recorded extracellularly. It is well-known that the spikes have better spatial and temporal resolutions than other signals. Therefore, these advances in the electrode implantation and brain network analyses have allowed us to map brain functional connectivity patterns in vivo.
A previous study has demonstrated that selective interaction among neurons is an essential component in neural information processing (Benchenane et al., 2010). Based on Hebb's theories, only a small fraction of neurons ("cell assembly") participate in a task, selecting and preserving most valuable information to conduct future behaviors. Previous studies have demonstrated that sparse-coding is a convincing hypothesis to explain the way to process neural information (Sakurai, 1999;Schweighofer et al., 2001). Several studies have also reported that brain networks are sparse (Sporns et al., 2005;Hwang et al., 2013). However, most studies have only focused on functional connectivity among all recorded neurons (Lin et al., 2009;Nedungadi et al., 2009;Kim et al., 2011). It is cautioned that spurious effects may be introduced by neurons not in assemblies. From a network perspective, it is necessary to separate present vs. absent connections. Many investigators have attempted to extract the true network features by several mathematics methods. Recent work shows that cell assemblies can be extracted, which is used to represent the intrinsic pattern of neuronal activity (Peyrache et al., 2009;Benchenane et al., 2010). Up to now, non-negative matrix factorization (NMF), principal component analysis (PCA) and independent component analysis (ICA) have been already used to extract the latent structure of fMRI data (Lohmann et al., 2007;Shen and Meyer, 2008;Demirci et al., 2009;Zhou et al., 2009). In our study, NMF is employed to learn the parts-based components of realistic objects with non-negative constraints (Lee and Seung, 1999). The whole perception in brain relies on its partsbased representations (Logothetis and Sheinberg, 1996;Lee and Seung, 1999). And neuronal firing rates are never negative. Based on the theories, NMF has found its way into the neuroscience community (Lohmann et al., 2007;Padilla et al., 2010;Liu et al., 2013). It can be seen that the notion of combining parts to form a whole is intuitive and easy-to-interpreted. Here, NMF was used to remove dimensions that did not contribute to the raw firing rate pattern of mPFC neurons. The necessary information was correctly retained by NMF components. Therefore, the method could not only find intrinsic feature of neural activity, but also reduce computational cost.
In summary, a deeper understanding of the mPFC network changes may have profound implications for revealing the neural mechanism underlying WM. In this present work, we examined connectivity dynamics in mPFC spike networks in both NMF and original spaces, in order to investigate: (1) the association between the network changes and WM behavioral choices; (2) whether the network feature can be better represented in NMF space.

SUBJECTS
Five male Sprague-Dawley rats weighting 300-350 g participated in this study. They were housed in temperature-controlled and clear plastic cages (25 • C) under a 12:12-h light-dark schedule. All Experiments were conducted during the light phase. The experimental protocol was approved by the Committee on Experimental Animal of Tianjin Medical University.

APPARATUS
The apparatus was a Y-maze ( Figure 1A), consisting of a start arm A and two goal arms B and C. Three arms were identical in size (80 cm long, 16 cm wide and 21 cm high) and other physical features. They were linked together in a symmetrical Y shape (120 • of angular deviation from each other). A starting platform on the arm A was separated by a guillotine door. Two food wells were located at the end of each goal arm, respectively. Meanwhile, there was an infrared detector at the fork of the Y-maze to record a reference time (RT). The RT referred to the time for choice behavior.

BEHAVIORAL EXPERIMENTS
The behavioral task was a spatial alternation Y-maze paradigm. First, all the rats were provided water ad libitum and food restriction for 2 days, maintaining at 85% of free-feeding body weight. Next, the rats explored the Y-maze and consumed reward food over two consecutive days (30 min per day). During habituation, they were forced to perform at least 15 entries to reward arms in 30 min. After that, the rats were trained to learn an alternative choice rule on the Y-maze. They were given one training session of 10 entire trials per day. An entire trial consisted of a free choice run and an alternative choice run. A latency for each trial referred to a period when a rat started from start point to reward site. During the free run, entering either goal arms was all rewarded. After ate, the rats ran back to the start location spontaneously and a guillotine door was closed. After 5 s delay, the door was lifted and the rats started to perform the choice run. They must choose a goal arm not previously entered to get reward. An error was defined by re-entering the arm visited in the last free run. On the way home, the arm that had not been visited was blocked by a guillotine door in order to prevent the rats quickly turning around and running into this arm. When the percent correct reached at least 80% over consecutive 2 days, the rats were considered as acquiring the rule used for the Y-maze task. If the rat has not learned the rule on day 8, it would be knocked out.

ELECTROPHYSIOLOGICAL RECORDING AND DATA PREPROCESSING
Following 2 weeks recovery, the rats were again conducted into performed the Y-maze task. At the same time, neuronal activity was recorded. The RT for each trial was detected by an infrared sensor mounted on the Y-maze. The latencies for all trials were also recorded. To be clear, the trials mentioned in following paper referred to the alternative choice trials which could be used to measure the WM performance. A total 114 trials were analyzed (n correct in total = 96, n error in total = 18. Rat 1: n correct = 20, n error = 4; Rat 2: n correct = 13, n error = 2; Rat 3: n correct = 18, n error = 3; Rat 4: n correct = 27, n error = 5; Rat 5: n correct = 18, n error = 4) in the present paper. Besides, the baseline activity was taken during the inter-trial intervals. We FIGURE 1 | Working memory behavior. (A) A working memory task on a Y-maze (overhead view). An entire trial contained a free choice run (left) and an alternative choice run (right). Left: Arm A is a start arm. Reward sites are located at the ends of two goal arm B and C. RT was recorded by an infrared detector (red triangle). Right: A black line represents a correct choice and a gray line represents an error choice. (B) As training processed, the average rate of correct choice gradually increased (n = 5 rats). (C) As training continued, the average latency (n = 5 rats) gradually decreased. (D) The latency distribution (from the guillotine door to the choice point, n Rat1 = 24 trials, n Rat2 = 15 trials, n Rat3 = 21 trials, n Rat4 = 32 trials, n Rat5 = 22 trials) for each rat. (E) Average latency (from the guillotine door to the reward site, n Rat1 = 24 trials, n Rat2 = 15 trials, n Rat3 = 21 trials, n Rat4 = 32 trials, n Rat5 = 22 trials) over entire trials for each rat.
compared the certain duration of mPFC neuronal activity (duration: 0.25 s pre-measure peak and 0.25 s post-measure peak) in correct trials with the baseline activity.
During the recording courses, the rats were connected to the Cerebus Acquisition system (Cyberkinetics, Foxborough, MA) via a rotating and lightweight cable. They were allowed to move freely on the Y-maze. Wide-band neural signals through the multi-electrode array were fed into a head stage amplifier with 5000 × gain and band-pass filtered between 250 and 7.5 kHz. The obtained signals whose amplitudes exceeded pre-set voltage thresholds were recorded, digitized at 30 kHz and stored for offline spike sorting. Pre-set voltage thresholds for spike detection were done on negative amplitudes for each electrode independently. The threshold level was set as -7.5 times the standard deviation of raw signals (Reyes-Puerta et al., 2014). Any waveforms that crossed the corresponding threshold were detected and stored with time stamps. Since neural activity per channel was actually synthesized by several neurons nearby an electrode tip, Off-Line Sorter software (Plexon) was applied to extract multiple spike trains associated with different neurons by a semiautomatic clustering selection method (Lewicki, 1998) followed by manual clustering.

HISTOLOGY
After electro-physiological recordings, all the rats were sacrificed for histological analyses. The animals anaesthetized with sodium pentobarbital (40 mg/kg) were perfused intracardiacally with saline followed by 10% formalin saline. Then the brains were removed from the calvaria and trimming. Serial, 150 μm thick, coronal sections were carried out. Tracks and lesion sites were observed under a light microscope.

CONSTRUCTION OF A CONTINUOUS TIME SERIES FROM A SPIKE TRAIN
Due to the discrete nature, spikes were not suitable for the subsequent connectivity analysis. A spike train must be converted into a continuous instantaneous firing rate time series (Zhu et al., 2005). The new series could encode the same temporal firing pattern of neuron as original spikes.
Specifically, the mathematical procedure consisted of three steps (Zhu et al., 2005). First, the instantaneous firing rate could be described by the following function: where t i − 1 , t i and t i + 1 were considered as three successive spike times. Next, a small time interval δT was chosen to obtain the instantaneous integrated rate as follow: Note that δT was much smaller than the mean interspike interval T. According to the recommendation (Zhu et al., 2005), δ was set to 0.25. Finally, a continuous instantaneous firing rate series was yielded by smoothing out above functions in time (by Lagrange interpolation methods). In this way, each trial was represented by a matrix which consisted of m neurons of n instantaneous firing rates.

EXTRACTED COMPONENTS IN NMF SPACE
X is an m-by-n input matrix. m is the number of neurons. n is the number of instantaneous firing rates. NMF was performed on the input matrix to decompose it into non-negative basis matrix (W) and non-negative encoding matrix (H), as follows: where FRO is known as Frobenius norm and e is residual error matrix (Paatero and Tapper, 1994;Lee and Seung, 1999). W is a representation of NMF space. H is a low-dimensional representation of the input data (Padilla et al., 2010). In NMF analysis, the estimation of rank r for factorization is a critical matter. It is required that at least: The NMF computation was iterated by varying the parameter r from 1 to m (the number of neurons). The number whose variance accounted for (VAF) was just greater than 90% was selected as the rank r (Torres-Oviedo et al., 2006). The VAF is expressed as follows: For each trial, there was a NMF space. The VAF and rank r for each one were calculated. The rank r varied between animals, since different numbers of neurons were recorded as well as neurons with different contributions. There were also different ranks constructed from the same animal in different trials, as different neurons were either active or silent in different trials, or patterns of activity changed. In order to perform statistical analysis on data, maximum rank r for each rat was determined as the number of NMF components which keeps the VAF for each trial more than 90%. That means the neural activity pattern was well accounted for by extracted NMF components. Thus, the rank r was consistent within an animal, across trials. According to the final rank, we recalculated the corresponding VAF for each trial.

Network construction
Granger causal connectivity analysis (GCCA) is a powerful approach for assessing functional connectivity among neural signals, revealing how one time series might exert influence on another (Zhou et al., 2013). If the past of one series does not contribute to the future of another one, the Granger causality (GC) between them is 0. When applied to a brain, this method has a potential to map the connectivity pattern of a brain network. A network is defined as a set of nodes linked by edges. In original space, nodes represent a set of neurons recorded in mPFC. In NMF space, nodes represent a group of NMF components. Edges represent the Granger causality values between their respective nodes. A 0.5-s moving window (moving step: 0.125 s) was applied to time series of every trials. Here we give a brief mathematical description about the Granger causality analysis (Seth, 2010). First, the multivariate time series can be described as follows: Next, a signal X is represented by a multivariate autoregressive model as follows: Frontiers in Behavioral Neuroscience www.frontiersin.org January 2015 | Volume 9 | Article 2 | 4 where p is the model order specified by the Bayesian information criterion, A(n) is the model coefficients, and E(t) contains the prediction errors for all time series, and ξ m (t) is the prediction error for x m . The E(t) covariance matrix of the unrestricted model is defined as follows: In a restricted model omitting variable x j , E(t) is represented as follows: Thus, the E R (t) covariance matrix of the restricted model omitting x j is: The Granger causality from variable x j to variable x i , conditioned on all other variables is given by: Having constructed a Granger causality matrix, it is important to assess the statistical significance. Only significant connections were retained by F-tests (Seth, 2010;Frye et al., 2012). The tests were corrected for multiple comparisons by means of the false discovery rate (Benjamini and Hochberg, 1995). Connections which were not statistically significant were set to zeros.

Network properties
In the present study, we focus on the global properties in mPFC networks during the WM task. Both causal density (Cd) and global efficiency (E) can provide overall descriptions of a network. Cd is a measure of causal interactions involving dynamical integration (neurons coordinate with each other) and differentiation (neurons contribute in different ways) among network nodes. Given a causal network identified by GCCA, Cd is defined as (Seth et al., 2006): where Num GC is the total number of significant connections, N is the number of nodes. The values of Cd are bounded between 0 and 1. A neural network only containing independent neurons would have low Cd. By contrast, a network with high Cd indicates that nodes globally coordinate with each other and are dynamically distinct (Seth, 2010) to accomplish a task at the same time.
E is a measure of the capacity for nodes to propagate information in parallel across a network (Hwang et al., 2013). In general, the integration feature of a network is based on paths, which represent potential routes of information flow between pairs of nodes (Rubinov and Sporns, 2010). The path length is the sum of edge lengths along this path. Every edge length is determined by calculating the inverse of connection strength. The shortest path length between node i and j is the path with minimum total length (Yan et al., 2011). E could be defined mathematically as Bernhardt et al. (2011): where L ij is the shortest path length between node i and j, N is the number of nodes. Previous studies had demonstrated that E was more suitable for sparser networks (Achard and Bullmore, 2007;Estrada and Hatano, 2008;Rubinov and Sporns, 2010). As the mPFC network in our study is sparse, inversing path length simplifies the numerical issues in the computation of E in sparser networks. Moreover, E measures the capacity to transfer parallel information. Since the brain supports large amounts of parallel information processing, E is a superior measure of network integration. A network having a high E suggests the existence of strong associations between nodes. Information is communicated with a high level of efficiency throughout the network.

STATISTICAL ANALYSIS
In this study, all data was expressed as mean ± standard error (SEM). Statistical analysis was carried out using One-or Two-Way ANOVA for repeated measures of the network features. T-test was used to compare the measures in the WM and baseline states. The significance level was set at P < 0.05.

NMF-GCCA ALGORITHM
In order to be more clearly, the general procedure for researching the mPFC network changes in original and NMF spaces was summarized as follows: (1) Obtain spike trains in rat mPFC during a WM task and convert each discrete train into a continuous time series.

NEURONAL ACTIVITY IN mPFC DURING WORKING MEMORY TASKS
Spikes in mPFC were recorded while rats performed WM tasks. Under a light microscope, the electrode tracks and lesion sites were unequivocally identified for 5 rats. The region of an implanted electrode array is provided in Figures 2A-D. Figure 3A displays representative spike activity in a correct trial. 15-channel spike trains were sorted, obtaining 24 single-unit spike trains ( Figure 3B). Based on sorting results, the relationship between channel numbers and neuron numbers is shown in Figure 3C. For example, two neurons (neuron 13 and 14) were isolated from the Channel 7 ( Figure 3D). Under this method, in total of 117 neurons (Rat 1: n = 17; Rat 2: n = 24; Rat 3: n = 32; Rat 4: n = 23; Rat 5: n = 21) were collected for analysis. The firing rates of correct trials were averaged using a series of sliding windows (width: 0.5 s, moving step: 0.125 s; Rat 1: n = 20 trials; Rat 2: n = 13 trials; Rat 3: n = 18 trials; Rat 4: n = 27 trials; Rat 5: n = 18 trials). One-Way ANOVA shows main effect of time on average firing rates (Rat 1-Rat 5, p < 0.05), indicating that the values significantly peaked before rats arriving at the maze choice point.  Figure 4A shows single-unit spikes during a WM trial. The spike trains were then converted into continuous-time series via the instantaneous firing rate method. Two representative signals (Neuron 8 and 16) are shown in Figure 4B. The results suggest that continuous series can describe the temporal firing patterns of neurons. Then the resulting time series were embedded in NMF space. Figure 4C depicts that VAF varied with rank r. By using the criterion previously described, the data shown in Figure 4A has a rank of 7. According to the selected rank, 7 NMF components were identified from this trial ( Figure 4D).

CAUSAL DENSITY DYNAMICS OF mPFC NETWORKS DURING CORRECT TRIALS IN ORIGINAL AND NMF SPACES
Similar to the mean connection weights, Causal density (Cd) increased to peaks in the interval right before rats' correct behavioral choices. The distinct increases were observed in both original and NMF spaces for 5 rats (Figure 6, from left to right; Two-Way ANOVA, factors: NMF / original space (S) and time on maze (T); Rat 1, n = 20 trials, main effect: S: p < 0.01, T: p < 0.05, interaction p < 0.05; Rat 2, n = 13 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.05; Rat 3, n = 18 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.01; Rat 4, n = 27 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.01; Rat 5, n = 18 trials, main effect: S: p < 0.01, T: p < 0.05, interaction p < 0.05). Additionally, Cd in NMF space was significantly higher than that in original space. The results indicated that the changes of Cd in correct trials were represented more obvious by NMF analysis.

GLOBAL EFFICIENCY CHANGES OF mPFC NETWORKS DURING CORRECT TRIALS IN ORIGINAL AND NMF SPACES
Global efficiency (E) changes of mPFC networks were examined in both original and NMF spaces for each rat (Figure 7). During correct trials, we found the trends of increased E occurring before the choice point in two spaces (Figure 7, from left to right; Two-Way ANOVA, factors: NMF / original space (S) and time on maze (T); Rat 1, n = 20 trials, main effect: S: p < 0.01, T: p < 0.05, interaction p < 0.05; Rat 2, n = 13 trials, main effect: S: p < 0.01, T: p < 0.05, interaction p < 0.05; Rat 3, n = 18 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.01; Rat 4, n = 27 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.01; Rat 5, n = 18 trials, main effect: S: p < 0.01, T: p < 0.01, interaction p < 0.05). The statistical results indicate that there were profound modifications in E of mPFC networks during correct trials. And like the Cd, E in NMF space was significantly greater than that in original space, indicating that the WM-related mPFC network dynamics was highlighted in NMF space.

mPFC NETWORK PROPERTIES IN CERTAIN WM AND BASELINE STATES
The above results demonstrated that the Cd and E in mPFC networks significantly increased before the maze choice point during correct trials. In order to examine whether the increases were meaningful changes, we compared the measures in certain periods (duration: 0.25 s pre-measure peak and 0.25 s postmeasure peak, n = 96 trials, 5 rats) in correct trials to measures in the baseline states (n = 91 trials, 5 rats). As shown in Figure 8, t-test reveals that there are significant differences between the certain WM and baseline groups ( Figure 8A, Cd: p < 0.01; Figure 8B, E: p < 0.01).

mPFC NETWORK PROPERTIES IN CORRECT AND ERROR WM TRIALS
Since Cd and E in mPFC networks changed dynamically during correct trials, we attempted to ascertain the network properties

DISCUSSION
This present study explored mPFC network changes during a WM task, employing multi-channel spikes. We embedded spike trains into a NMF space, followed by network analysis in this transformed space. The principle findings here are (1) Cd and E of mPFC networks changed dynamically while the rats correctly performed WM tasks; (2) before the rats arrived at the maze choice point, Cd and E had already augmented to peaks; (3) mPFC network properties were highlighted by NMF analysis. Furthermore, mPFC network measures of certain WM states in correct trials are significantly higher than measures in baseline states. And we did not observe increased network measures in error trials. The prefrontal cortex is at the top of the cortical perceptionaction cycle. WM is supported by this high-order brain area. In the Y-maze task, rats must keep reward location in free choice run "in mind" for short periods, which is used to execute subsequent alternative choices. With current multi-unit recording techniques, we can simultaneously record spikes from many neurons and track their activity changes directly. Previous studies have demonstrated that neurons in mPFC maintain WM-related information through their synaptic connections (Hains et al., 2009;Paspalas et al., 2013). The connections enable mPFC networks to efficiently process WM-related information. By studying the functional properties of mPFC networks, we may gain insight into the true nature of mPFC that underlies WM. For recent years, neurophysiologists have studied functional interactions among neurons during special cognitions (Lin et al., 2009;Stevenson et al., 2009;Hirabayashi et al., 2010;Masud and Borisyuk, 2011). Cd of certain WM states in correct trials (duration: 0.25 s pre-measure peak and 0.25 s post-measure peak, n = 96 trials, 5 rats) is significantly greater than that in baseline states (n = 91 trials, 5 rats, t-test, * * p < 0.01). (B) E of certain WM states in correct trials (duration: 0.25 s pre-measure peak and 0.25 s post-measure peak, n = 96 trials, 5 rats) is significantly higher than that in baseline states (n = 91 trials, 5 rats, t-test, * * p < 0.01). In our study, well trained rats made choices at the maze choice point, on the basis of a prediction about which goal arm was the most likely rewarded one. The results show that Cd and E of mPFC networks change dynamically in correct trials. When rats started from the start point, the measures were gradually rising over time. Their peaks all appeared before the choice point. Then the measures begun to fall after rats entered the goal arms. The statistical results suggested that Cd or E in the certain durations (0.25 s pre-measure peak and 0.25 s post-measure peak) of correct trials was significantly higher than the measures in baseline states, indicating that there were real changes in Cd and E. Additionally, the significant increase pattern was not found during runs toward the choice point in error trials. mPFC network dynamics was different between correct and error trials, suggesting an association between the network states and correct behavioral choices. It is inferred that increases in mPFC network measures might have an influence as rats assessed the possible results of choices, consistent with a previous WM study (Rossetti and Carboni, 2005).
Methods for accurately estimating functional connectivity may be important for many neuroscience issues. Brain neural activity is sparse which is the point of the feature space transformation. Low-dimensional representations chosen from high dimensional data can preserve or highlight some interested features in the data (Cunningham and Yu, 2014). NMF is one of the dimensionality reduction methods. Previous studies had demonstrated that NMF could extract the dominant features in fMRI (Lohmann et al., 2007), EMG (Chiovetto et al., 2013), SPECT (Padilla et al., 2010), ECG (Ghoraani et al., 2011) database and odor profiles (Castro et al., 2013). In our study, the network analysis in NMF space presented here approaches an efficient description of mPFC network properties. The results in Figure 5A show that mPFC networks in original space are naturally sparse, in that many Granger causality values are exactly equal to zero. It is possible that some dimensions are redundant, resulting in the latent neuronal activity being confined to a low-dimensional slice of original space. NMF was used to remove dimensions that did not contribute to the raw firing rate pattern of mPFC neurons. The necessary information was correctly retained by NMF components. In Figure 5B, the average connections of mPFC networks are significantly greater than those in Figure 5A, indicating that featured mPFC activity may densely occupy the NMF space. Therefore, the functional connectivity among NMF components may offer a more accurate description of mPFC network features. The dynamics of Cd (or E) in two spaces tracking one another means that NMF components preserve the necessary information in original data. Higher Cd (or E) means that the network properties in NMF space are more apparent. The results suggested that NMF analysis faithfully captured the patterns of mPFC networks and made them more obvious. Our approach enables us to study network dynamics from a novel perspective that are not observable in original space, opening a possibility of better identifying the network characteristics.
Numerous researches have shown that mPFC is linked with WM. In the present study, food-restricted rats were rewarded if they chose the goal arm not entered in the previous free run. In this process, mPFC recalled and integrated task-relevant information to guide behavior appropriately. Previous electrophysiological evidence has suggested that lesion or inactivation of mPFC often caused increased error responding and longer retention delays in rats performing delayed spatial alternation tasks (Wang and Cai, 2006;Yoon et al., 2008;Horst and Laubach, 2009). From a molecular and cellular perspective, Kim and colleagues also have revealed that mPFC is responsible for WM by using mice with mPFC-selective PLC-β1-knockdown . When carried out a series of behavioral tests, they found that the mice with impaired mPFC were incapable of acquiring the task. However, rats did not depend only on the WM to perform trials during the Y-maze task. The other process might be involved. For instance, rats should inhibit immature response to the previously rewarded arm to get reward in the alternative choice run. Therefore, whether the observed mPFC activity also presumably reflects and drives many aspects of behavior (such as performance monitoring and cognitive control) deserves discussion. Previous study has suggested that the activity of a significant percentage of the mPFC neurons did represent and maintain prior task-relevant information to guide appropriate WM task behavior (Jones and Wilson, 2005). Moreover, a well-controlled study has demonstrated that elevated delay-period activity in mouse mPFC appeared to play a more important role in memory retention than inhibitory control, motor selection, and decision-making (Liu et al., 2014). These studies, together with the results reported here, support that the recorded mPFC activity closely linked with WM during Y-maze tasks.

CONCLUSIONS
In summary, we have studied changes of mPFC spike networks during WM by means of NMF analysis. The findings point to a possible mechanism wherein the significant increases in mPFC network measures might conduct rats' subsequent choice behaviors. The anticipatory changes can be predictors of rat's future WM performance. Further results indicated that by considering the sparse feature of neural activity, the NMF method can provide a better characterization of network dynamics with less computational cost.

ACKNOWLEDGMENTS
Xin Tian designed the experiment; Jing Wei, Wenwen Bai and Tiaotiao Liu carried out the experiments; Jing Wei analyzed the data; Jing Wei and Xin Tian wrote the manuscript. All authors read and approved the final manuscript. This work was supported by the National Natural Science Foundation of China (91132722 and 61375113).