Classification of Fixed Point Network Dynamics from Multiple Node Timeseries Data

Fixed point networks are dynamic networks encoding stimuli via distinct output patterns. Although, such networks are common in neural systems, their structures are typically unknown or poorly characterized. It is thereby valuable to use a supervised approach for resolving how a network encodes inputs of interest and the superposition of those inputs from sampled multiple node time series. In this paper, we show that accomplishing such a task involves finding a low-dimensional state space from supervised noisy recordings. We demonstrate that while standard methods for dimension reduction are unable to provide optimal separation of fixed points and transient trajectories approaching them, the combination of dimension reduction with selection (clustering) and optimization can successfully provide such functionality. Specifically, we propose two methods: Exclusive Threshold Reduction (ETR) and Optimal Exclusive Threshold Reduction (OETR) for finding a basis for the classification state space. We show that the classification space—constructed through the combination of dimension reduction and optimal separation—can directly facilitate recognition of stimuli, and classify complex inputs (mixtures) into similarity classes. We test our methodology on a benchmark data-set recorded from the olfactory system. We also use the benchmark to compare our results with the state-of-the-art. The comparison shows that our methods are capable to construct classification spaces and perform recognition at a significantly better rate than previously proposed approaches.


I. INTRODUCTION
Robust neural networked systems encode their dynamics by producing attractors in a low-dimensional state space.The attractors represent neural network response to various inputs and reflect particular states of the system.Such networks are common in neuronal systems that process sensory stimuli, command motor systems, or store memory [2], [6], [20].The manifestation of attractors ensures robustness of network performance, such that, for a range of network initializations and stimuli, it exhibits reliable dynamics.The simplest type of attractors are input induced fixed points, triggered by input signals (e.g., step functions) into a subset of network nodes, and as a response after transient dynamics, a subset of the network output nodes produce a steady state pattern.In neuronal networks, these patterns are identified as neural codes [3].The network is considered selective when it is capable to distinguish between stimuli by producing distinct fixed points.In particular, the network will produce similar fixed points for similar stimuli, and distinguishable fixed points for distinct stimuli, where similarity is defined by a metric in a low-dimensional space.Such functionality is the primary principle upon which fixed-point networks incorporate recognition and quantification of mixed stimuli.In addition, these networks are typically 'robust' as they are capable of perform under high dynamic input variability, i.e., low signal to noise ratio (SNR).
A fascinating question is how to infer the lowdimensional state space and fixed points within it from sampled time-series data of the network.It is a particularly relevant problem when it is of interest to characterize the functionality of black-box networks: which connectivity and node dynamics are unknown, and probing network response for various stimuli is the only resort.A structural approach for construction of low-dimensional state space is supervised classification.In such an approach, a set of distinct inputs is being applied to the network independently, and sampling of network activity, i.e., time-series of multiple node dynamics, are being obtained.Thereby, for each input, the fragment of sampled time-series corresponds to a matrix with 'nodes' and 'time' dimensions, and the whole database of sampled responses corresponds to a collection of matrices (see Fig. 1).Classification task for such a collection is to find the most informative lowdimensional state space where the number of distinct fixed points is the number of presented distinct inputs, and that the low-rank space could be used for examining the transient dynamics reaching these fixed points and their superpositions.Such a space would be considered optimally selective when the both the fixed points and the transient dynamics leading to the fixed points are maximally separated (orthogonal).For example, for a collection that constitutes responses to three inputs, we expect to find three basis vectors each representing a single stimulus.
For neuronal sensory networks, the collection would typically correspond to instantaneous firing-rates of neurons (peri-stimulus time-histograms) obtained from multi-neuron recordings when the neural system is stimulated.A particularly intriguing system is the olfactory neuronal network in insects.Within the network, the antennal lobe, primary processing unit for olfactory neural responses, is receiving input from olfactory receptors and is able to respond with output activity which discriminates between odorants and classifies them into similarity classes.Experiments have shown that such classification is associated with behavior, and could adapt over time or after training [14], [15].Furthermore, analysis of multiunit electrophysiological recordings from output (projection) neurons shows that the network employs input induced fixed points to classify the olfactory stimuli [12].
Since the collection consists of matrices, a natural methodology for inference of the low-dimensional space is structuring all sampled response matrices into a single matrix and employing classical multivariate matrix decomposition methods, such as Singular Value Decomposition (SVD), or sparse representations computed by L1 minimization, that will identify dominant orthogonal patterns [18], [19].However, in practice, due to the ambiguity in SVD (i.e., indifference for the structure of the data and sensitivity to low SNR in the inputs), the resulting modes are mixed and non-uniquely discriminate single inputs from the collection.Alternatively, application of dimension reduction methods for each of the sub-matrices individually produces unambiguous lowrank representations -however it introduces a problem of gathering the obtained patterns into a single, joint, meaningful representation.This is a generic problem of finding an optimal low-rank representation of multiple instances data matrix.
To address the problem, here we propose a classification method that leverages the benefits of the two approaches: orthogonality and non-ambiguity with respect to the structure of the data.The method, called exclusive threshold reduction (ETR) operates on the low-rank representations obtained from individual applications of data reduction applied to each sub matrix, and creates an orthogonal basis.Effectively, the method exclusively associates each neuron with one of the stimuli and assigns a weight for the association.We show that such an approach can successfully separate response trajectories to the various stimuli when the response matrix is projected onto the basis.In addition, we formulate an optimization routine, called optimal exclusive threshold reduction (OETR) which allows us to achieve maximal separation of fixed points.For evaluation and demonstration of our methodology we have obtained a benchmark database of multi-unit time series recordings from the antennal lobe projection neurons in the Manduca sexta moth.The recordings were obtained from the neuronal network responses to stimulation of constituents of the Datura flower scent (single molecules and mixtures), a major food resource for the moth.The dataset is divided into a training set for construction of classification space, and a testing set to evaluate our methodology and compare with other approaches.In particular, an important test is determining how well the approach discriminates between distinct stimuli and captures similarity of various mixtures -for example, the method should be able to classify the network response to a stimulus as whether it is behavioral or non-behavioral neural code.We show that the OETR method allows such performance for various number of dimensions of classification space and mixture inputs.

A. Classification from Multiple Node Collection
For supervised classification we consider a set of input nodes {I} which receive a set of stimuli {S}.To monitor network response, we consider a set of output nodes denoted as {F }.In the training stage of supervised classification, m independent stimuli {S i } i=1,...,m are applied to the network and the timeseries of output nodes that produce dynamics attracted to stable fixed points are being recorded.Per each stimulus S i , output nodes timeseries that have been recorded correspond to a matrix [F Si ] of dimensions N × T with rows being the nodes and columns are recording time stamps.Therefore the training stage of a supervised classification produces the collection C = {[F S1 ], ..., [F Si ], ..., [F Sm ]}, which includes multi-node times-series of network responses (see Fig. 1 for an example).The goal of the classification is to find an optimal representation such that the responses to distinct stimuli in the collection are separable.Effectively, the problem is related to finding a separating hyperplane between each fixed point and its S i (here i = 2), are injected into network's input nodes and produce fixed point dynamics in output nodes (n1 and n2 are samples of output nodes).For each stimulus output nodes timeseries are recorded in a response matrix, with dimensions of nodes × time, where each node is a row of the matrix.Response matrices to S 1 and S 2 stimuli are shown as color maps and time series of sampled nodes n1 and n2 (rows 1 and 2 of each response matrix) are shown in insets.As the number of output nodes increases, finding a separating hyperplane between the dynamics represented in the nodes space becomes complex as illustrated in the rightmost plot.associated trajectories that approach/leave it and all other fixed points and their associated trajectories.In practice, the recordings are noisy since incorporate large O(1) variability in time and space of the input.
To quantify the success of various classification methods we have established a benchmark data-set for classification and recognition.The benchmark consists of a collection of responses of neurons within a neurobiological network that exhibits fixed point dynamics.Specifically, to create the benchmark, we have obtained electrophysiological multi-neuron recordings from projection (output) neurons within the antennal lobe neuronal network, the primary olfactory processing unit in insects.It was shown that the response of this neuronal network to odor stimuli is expressed in terms of fixed point dynamics of projection neurons [12].To obtain the collection, spiking activity of multiple projection neurons (N=106) was recorded from the antennal lobe of Manduca sexta moth subject to distinct odor molecules (odorants) and their mixtures extracted from the Datura wrightii flower, a major food resource for the moth.Since realistic response times are of order of few hundreds of milliseconds (200-400 msec), each response was recorded for 1 sec.The recordings were repeated over 5 trials with restoration intervals in between the trials.
The benchmark includes responses to 8 odorants (labeled as: bea, bol, lin, car, ner, far, myr, ger), 1 control odor (labeled as: crl), and 8 mixture odors (labeled as: e2, e3, e3b, p2, p3, p4, p5, p9), see Table in Fig. 2 for more details.Mixtures are labeled as 'behavioral' -identifying that the behavioral response to these stimuli is similar to response to Datura floral scent -, or 'non-behavioral'to denote the mixtures and odorants which do not elicit a significant behavioral response.To work with continuous time series, neural spiking activity was transformed to peri-stimulus time histograms (instantaneous firing rates) which are nonnegative values.
With the established benchmark we first demonstrate that finding a hyperplane that separates the data points is not obvious due to the high-dimensionality and variability (noise) of the data, thus precluding classification.Therefore, methods to directly identify such a separation are not expected to yield a precise classification.To test several standard methods with the benchmark, we first used the Support Vector Machine (SVM) classifier [5], [13] with a binary classification scheme for which the data was separated into two classes: responses triggered by stimuli of interest and responses triggered by other stimuli.SVM is supposed to classify the responses by finding the optimal hyperplane that separates all data points into one of the two classes (see further details on the method in the Appendix).For this binary classification problem we found the performance to be low with average classification accuracies around 50% for various binary classes, as we show in Fig. 3.We find that both precision (percentage of true positives out of instances classified as positive) and recall (percentage of true positives out of instances expected to be positive) errors are ≈ 60% and ≈ 50% respectively, indicating poor performance in the usefulness of the classification (precision) and completeness of classification (recall), comparable to a random guessing strategy for the binary test sets that we have created.
A possible explanation for the failure of SVM on the benchmark is that the points are not linearly separable.We therefore tested nonlinear versions of the classifier, e.g., by using a Gaussian kernel, however they produced similar classification errors.Another hypothesis for low performance of SVM could stem from response dynamics to reside in a lower dimension than both the nodes dimension and the number of data points upon which the separating hyperplane is computed.This creates a situation where the data is being overfit.Another, obstacle could be in the form of the classes being imbalanced, i.e. when SVM classifiers are trained on an imbalanced dataset, they can produce biases towards the majority class and misrepresent the minority class.We therefore tested other techniques designed to deal with class imbalances, such as SMOTING and RUSBoost [16] with the benchmark.However, we did not obtain significant improvement with these methods than SVM classification (Fig. 3).We therefore conclude that direct classification from the complete data-set is difficult and approaches that incorporate pre-processing of the data are necessary.In particular, we identify that finding an appropriate low dimensional state space, where distinct fixed points and their associated dynamics are easily separable, could significantly simplify classification.

B. Classification using Matrix Decomposition Applied to Timeseries Collection
Since classification methods applied to the highdimensional collection C are unable to provide robust representation, an alternative approach is to represent the data in low-dimensional space.The structure of the collection C as a set of matrices suggests that matrix decomposition could be a useful tool for finding such low-dimensional basis.There are several possibilities to create a matrix to be used for decomposition from the collection C. The first possibility is to format the collection into a concatenated matrix consisting of submatrices Data reduction methods can be then applied to the matrix F to decompose it into vectors a k (t i ) (time-dependent coefficients) and g k (j) (spatial patterns/neural codes): Our goal is to find the most informative decomposition where the number of patterns is the number of distinct presented inputs and they span a low rank space of projections in which the dynamics of these patterns, i.e., time-dependent coefficients a k (t i ), and their superpositions could be examined.For example, for a matrix that constitutes the responses to three inputs we expect to get three pattern vectors each representing a single stimulus.Classical multivariate decomposition of the matrix F, e.g., Singular Value Decomposition (SVD) / Principal Component Analysis (PCA) that decompose the matrix F into U ΣV T , or sparse representations computed by L1 minimization, can identify dominant orthogonal patterns [7], [18].
When computing the SVD for the matrix F, the procedure we call SVDCon, the column vectors of orthonormal matrix V T correspond to pattern vectors.The k-th column vector of V T , the vector g k (j) (or denoted as P C k in Fig. 3), is the k-th pattern.Each pattern has associated time-dependent coefficients vector a k (t i ) = σ k α k (t i ), where α k (t i ) is the k-th row vector of the orthonormal matrix U and σ k is the k-th singular value, the k-th element of the diagonal matrix Σ. Singular values indicate the weight of each pattern vector, since they are nonnegative elements ordered as σ 1 ≥ σ 2 ≥ • • • ≥ 0 and since both α k (t i ) and g k (j) are normal vectors.To determine the relative weights of pattern vectors it is possible to define the relative energy yields that the total energy is normalized n k=i E i = 1, and each E k indicates the relative 'percentage' of the energy in the pattern P C k .The distribution of E k values provides an estimate for the effectiveness of the decomposition to find a low dimensional representation.When several first E k values stand out it identifies that the patterns corresponding to these singular values are more dominant than others and the truncation of the remainder modes would maintain a reasonable approximation of the original matrix.
We apply SVDCon to a concatenated matrix F generated from three response matrices,F = [F S1 , F S2 , F S3 ] to three distinct odorant stimuli from the benchmark set.We observe that E k values distribution is such that only the first singular value stands out (E 1 = 0.5) and all other are significantly smaller.Such distribution typically identifies that SVD was not able to capture the variability in data.Furthermore, examination of time dependent coefficients associated with the first three patterns shows that they are overlapping and indistinguishable (Fig. 3).
Effectively, these results indicate that there is discrepancy between the expected three dimensional state space and the outcome of SVDCon.
Another possibility to format the collection to a matrix form is to consider each response matrix F Si separately, i.e., Using this form, SVD can be applied to each matrix, procedure we call SVDSep.In this procedure m sets of decompositions are obtained, where m is the number of distinct stimuli.We apply SVDSep onto three stimuli collection, as for SVDCon, however here it is formatted into three separate response matrices and obtain 3 decompositions.We observe that each decomposition is dominated by 1 pattern (E 1 > 0.5 for each of the three matrices).Therefore collecting dominant pattern vectors from all decompositions into a common set could represent a projection space for the distinct responses.More generally, a procedure following SVDSep would take the first pattern vector P C 1 Si (i.e.most dominant mode) from each decomposition of [F Si ] and store them as column vectors in a library matrix L, which has the form: where there are n rows accounting for every node measured in the network.Notably, however, column vectors of L are taken from separate decompositions and therefore non-orthogonal to each other and do not share the same Euclidean space.Therefore projections onto the column vectors of the library L are not guaranteed to be non-overlapping.Indeed, we observe that that projection of the response matrices onto L for the constructed library L from the benchmark experiment of three distinct odorants, does not produce separable fixed points and associated trajectories (see Fig. 3).The nonorthogonality of the basis warrants development of an approach to transform the library so that the projected fixed points and their trajectories are maximally separable.

C. Optimal Exculsive Threshold Reduction Method
To find separation between fixed points associated with distinct stimuli, we propose a simple method that will obtain an orthogonal set of vectors from the matrix L. The method, called exclusive threshold reduction (ETR), operates on the nodes dimension of L (rows): it selects the maximal value in each row and sets to zero all other elements in that row.The procedure is performed on each row, and once completed, leaves a new matrix with N nonzero elements out of N 2 elements, effectively associating each node with a stimulus.For example, if for node n i the element in the third column is the maximal element then node n i is associated with S 3 and the weight for the association is the value of the element, see top row of Fig. 4 for graphical illustration of the method.Therefore, the ETR method both associates nodes/neurons with stimuli, the different axes in Fig. 4, and assigns weights to each association.By definition, ETR guarantees an orthogonal set of vectors, i.e., an orthogonal matrix O of the form where τ is the threshold value.If τ = 0, as in many applications, the generated matrix O is a basis for the nodes space and a projection space for trajectories associated with stimuli.The matrix defines a mapping from the high-to low-dimensional system.Such a mapping was used to recover network connectivity in conjunction with proper orthogonal decomposition of population model equations for the antennal lobe neuronal network [17].
The mapping is supposed to group node responses and capture the exclusive features associated with each stimulus.When we test the projection of three odorant benchmark matrices onto the matrix O, producing timedependent trajectories in stimulus space, we observe that the trajectories are separate from each other.We also locate the fixed points, which coordinates are defined as the outcome of the product L T O.We find their location to be close to the stimulus axis of each trajectory (Fig. 4A).These experiments indicate that the ETR method is effective in mapping distinct trajectories to their own axes, and can facilitate a distinct classification space.
To test the scalability of the approach we apply ETR on various subsets of responses to odorants in the benchmark set.We start with 2 odorants and increment by one the number of included matrices to 8 odorants.Application of ETR on the latter produces 8 dimensional space (eight column matrix O).In these experiments we further confirm the scalable performance of the approach and its ability to produce separate fixed points and trajectories (as shown in Fig. 4A).
While ETR turns out to be valuable in separation of trajectories, it does not ensure generically, especially when the number of matrices in the collection grows, that the fixed points are optimally separated and hence noisy trajectories are attracted to the fixed points.To optimally separate the fixed points we propose to introduce additional weights contained in the diagonal matrix D The formulation assigns weights that scale the weight of each individual node obtained by ETR to ensure that the fixed points are orthonormal (Fig. 4B) and thus satisfies the assumption of distinct stimuli used in supervised classification.With such scaling fixed points coordinates are computed as L T D w O and in order for them to be orthonormal are expected to be exclusive in the associated axis of each stimulus.We therefore require that the coordinates of the fixed points will be represented by the identity matrix, i.e., L T D w O = I.To solve for the weights, we formulate the following convex optimization problem We denote the generalized approach of solving the optimization problem above, in conjunction with the ETR method, as Optimal Exclusive Threshold Reduction (OETR) method.It is expected to produce optimally separable projected trajectories for multiple distinct stimuli and hence the vectors of D w O are optimal axes of the classification state space.To solve the optimization problem, computational packages for convex optimization can be used.For example, we have used the CVX package implemented in MATLAB [8].Application of OETR to the three distinct stimuli benchmark set indeed produces more optimally separable fixed points and their associated trajectories as shown in Fig. 4B.Next, we compare OETR with other approaches, such as the Independent Component Analysis (ICA), to exclude that the benchmark set is too simple problem of separation.The ICA method that we apply is information based algorithm (Infomax) particularly designed to obtain separable signals from a collection of inseparable signals, such as the library L, see Methods for more details on our ICA implementation and based on [9]- [11].With ICA, our results show that projections on the obtained vectors remain to be overlapping and do not produce efficient classification (Fig. 4C).

D. Recognition Metrics and Classifiers based on the Classification State Space
With the classification space constructed using ETR and OETR we consider how it can be utilized for recognition and classification of novel stimuli.These stimuli include mixtures of odorants upon which the axes of the state space were constructed and responses to odorant stimuli on trial by trial basis, where each trial is a novel stochastic realization of the trajectory than the one used for space construction.Since the fixed points and their associated trajectories are well separated in the classification space, we propose the convex hull metric defined by the fixed point and the associated trajectory.A simplification of the convex hull metric (and more robust) is an m-dimensional hyperellipse centered at the fixed point where m is the dimension of the classification state space.The metric s t is a pointwise metric which provides for each point x 1 , ..., x m of the tested projected trajectory, sampled at time t, whether it lies inside the hyperellipse (q ≤ 0 : s = 1) or outside (q > 0 : s = 0).Integration of the pointwise metric over a time interval of interest, typically the time of the response, e.g., for the benchmark here it is the time of the stimulus being ON (500 ms) provides a total score S of the trajectory being located within the hyperellipse and indicates the similarity of the tested trajectory with the stimulus and its associated fixed point and hyperellipse.Normalization of the total score S over the total number points provides a recognition ratio score Rec.Such a metric is simple to implement and is specially designed to harness the optimal separation of fixed points in the classification space and expected to provide consistent scores for testing similarity between different trajectories in the classification space by being robust to different stochastic realizations of the trajectories (due to noise and other perturbation effects).Other metrics, measuring the time scales of convergence to the fixed point or focusing on specific properties of the trajectories are not expected to be robust for low signal-to-noise ratio (SNR) dynamics (see Fig. 5).Such properties of the dynamics are typical to multi-unit recordings from neuronal networks.For example, in the benchmark set SNR is less than 3 and stimuli trajectories exhibit short timescales and do not necessarily converge to the fixed points.Rather, they only approach their vicinity.
To investigate the performance of recognition and classification using the hyperellipse metric we computed the similarity of various stimuli (17 stimuli from the benchmark set) with B1 mixture from the benchmark consisting of the odors: Benzaldehyde (S1); Benzyl Alcohol (S2); and Linalool (S3).This mixture was shown to include the dominant constituents of the Datura scent -a important flower nectar source of the moths -and moths stimulated with it elicit a 'behavioral' response similar to stimulation by the full Datura scent, and similarly to mixtures B2 and B3 [14].The goal of the classification is examine the mixture trajectories and determine which are similar to B1.Strong similarity indicates that the mixture is classified as 'behavioral' as well.The goal of recognition is to decide instantaneously given a single trial whether the stimulus is B1.Notably, behavioral experiments show that similarity in responses to these different mixtures was not solely due to the concentration of the constituents making up the mixture, and therefore classification methods from neural dynamics are in need.For example, experiments show that when Linalool (component with the least concentration) is taken out of B1 mixture the new mixture becomes non-behavioral.On the other hand, adding other constituents from the Datura scent replacing Linalool does not restore the ability of the mixture to elicit behavior.
To compare between different approaches we constructed classification spaces using four methods: SVD-Sep, ICA, ETR and OETR for the odorants S1-S8 (see table in Fig. 3).We also varied the dimension m of the classification space from m = 1(S1) incrementally to m = 8(S1, .., S8).For each method and dimension of Fig. 5: Classification and Recognition Employing the Classification Space.A: ETR projection of response that corresponds to B1 (blue) and its associated 3D hyperellipse (black sphere) that is used for recognition of B1.All points that fall inside the hyperellipse are marked as "recognition evidence" (s t = 1).For reference we also show the projections of responses to B1 constituents, S1, S2, S3, in gray.the classification space we computed the score Rec for each trajectory (5 trials of 17 stimuli = 85 trajectories) with respect to hyperellipse that corresponds to B1 stimulus (sphere with fixed radius).We show our results in Fig. 5.In classification, for each dimension of the space, Rec scores were computed for all stimuli and trials.The average Rec score was then computed over the trials of each stimulus.The distribution of average scores { Rec } for 17 stimuli is then normalized by the maximal max{ Rec } value.A decision line for binary classification of 'behavioral' (B class) or 'nonbehavioral' (S or E classes) is set as the midpoint between the mean of the distribution and maximal value: d = ( { Rec } + 1)/2 (depicted as blue line in the inset showing distributions for m = 7 in Fig. 5).Values above d are classified as behavioral while values below d are classified as non-behavioral.In Fig. 5 the bars that correspond to erroneous classification are marked with red color and correctly classified bars are marked with green color.We show the total accuracies of the four methods (precision × recall) for each dimension as a bar plot in the middle plot of Fig. 5 and also show the distributions for m = 7 below it.
Classification results demonstrate that ETR and OETR are significant improvements over other dimensional reduction techniques, which in turn are significant improvement over methods that do not employ dimension reduction.All methods perform poorly when the classification space is of m = 1 or m = 2 dimensions.Indeed, since B1 consists of 3 independent stimuli, the results indicate that there is not enough data for discrimination based on less dimensions than the constitutents of the tested stimulus.As the dimension increases OETR achieves perfect accuracies for m ≥ 3. ETR achieves 100% accuracy for m = 3 and then stabilizes on lower value of ≈ 70% with high precision but lower recall rates.ICA achieves 100% accuracy for m = 4 but when the dimension increases both precision and recall rates drop to values as for m = 1 and m = 2 dimensions.SVDSep produces constantly low accuracies for all dimensions.From the distributions, as shown for m = 7, we get insights on the differences in the performance of the methods.In particular, two stimuli are used for validation: B1, which is expected to produce a high Rec score, and E6 (control), which is expected to produce a low Rec score.We observe that the scores for B1 are high across methods, however in SVDSep B1 does not cross the threshold d line.For E6, only ETR and OETR successfully assign a significantly lower score to E6.As expected, individual odorants S1,..,S8 and nonbehavioral mixtures E1,...E6 are far from threshold in OETR method and mostly far from threshold in ETR.The form of the distributions transforms from nearly uniform (for SVDSep) to spiked distribution with high values for B stimuli and low for other (for OETR).The latter form allows for the line d to easily separate the B class from other classes.Indeed we observe that while ETR is able to perform with high precision, it is unable to recall one of the behavioral stimuli (B2).Next, we proceed and test whether the methods are sensitive to the choice of the hyperellipse.In particular, we vary the radius of B1 sphere in the range of 0.5 -0.85 for dimension m = 8 and compute classification accuracies.We observe a clear separation between SVDSep/ICA and ETR/OETR groups of methods.We find that recall rates are stable across methods, however precision rates are sensitive.Precision rates of SVDSep and ICA are less than 50% overall and indicate that these methods are not robust.By contrast, ETR and OETR are more robust.ETR achieves 100% precision for a range of 0.19 radii and OETR achieves 100% precision in a much wider range of 0.32.
For recognition task, we compute the Rec score for m = 8 classification space for each trial (85 trials).If the score crosses a threshold of 70% of the Rec score of a target averaged trajectory (i.e. for more than 400 msec the trajectories are in the same hyperellipse) the tested trajectory is recognized as associated with the target stimulus.Using this method, we test recognition for B1 hyperellipse (sphere of radius 0.65 in our case).As in classification, our results show that OETR is more accurate than other methods.Recall and precision rates of OETR and ETR are higher than the rates of SVDSep and ICA.Interestingly, when we test for recognition of B1 stimulus only (i.e.we expect to recognize only 5 noisy trajectories of B1 stimulus as positive) precision rates do not reach 100% in any of the methods; there are other trajectories marked positively as B1 although not B1.These trajectories turn out to be mainly from B class.We show that by expanding the class of recognition and to allow for any stimulus in B class to be marked as B1.In such a case, we obtain 100% precision for OETR and ETR.These results demonstrate that noisy B trajectories appear extremely similar to each other in our recognition scheme.OETR is hence consistent with experimental observations and metrics which indicate that B class stimuli are behaviorally indistinguishable.

III. CONCLUSION
Classification of fixed point networks responses from sampled network activity is a challenging fundamental problem, especially when the inputs into a network are highly variable and noisy.We have addressed this problem by proposing novel supervised classification methods which find a low-dimensional representation from a collection of matrices of multiple node time-series data from the network.To test our methods we have established a benchmark database of multichannel timeseries recordings from a real neural fixed point network: the antennal lobe in the Manduca sexta moth.We have shown that the OETR method that we have introduced allowed us to create a classification space that separates fixed points and their transient trajectories with high accuracy.Furthermore it allows to represent and classify noisy mixtures of stimuli.In contrast, we have shown that traditional methods such as SVM and Boosting that do not rely on dimension reduction and classical matrix decomposition and reduction methods such as SVD and ICA do not succeed at classifying fixed point networks dynamics with comparable performance to OETR.
Recordings from the olfactory network are valuable benchmark for testing classification.The network is known for its robustness in discrimination of scents composed of numerous distinct molecules in a highly dynamic environment.It has been established that the network employs fixed points for odor representation that are being read and processed by higher layer olfactory processing units (mushroom body in insects).In this respect, the ETR and OETR methods that we have proposed are simple to implement since they rely on dominant mean pattern (first SVD mode) associated with each distinct odorant, maximum rule and linear weights.Such components are natural to neural networks.Thereby application of the proposed methods could help in understanding the processes that higher layer units employ to instantaneously read information from fixed point networks.Here we focused on fixed point networks.These networks are fundamental and simplest attractor networks.It is thereby plausible that our methodology will lay the foundations for future work to extend the methodology to other more complex attractor networks such as limit cycle networks, lines of fixed points networks that are also ubiquitous networks in neural systems.

IV. METHODS AND PROCEDURES
The benchmark collection and code that implements the various methods is available online on GitHub: https://github.com/shlizee/fpnets-classification.
Below we include further information on the algorithms that we have applied and compared with ETR and OETR methods.

A. Class Imbalances
If the error or misclassification percentage is the only metric taken, learning algorithms like SVM appear to work well, but this is not an accurate result.The reason the error is so low, is because the first class is not classified at all and the second class is perfectly classified.Therefore where tp is true positive, fp is false positive, and fn is false negative.Another tactic, is re-sampling the data set.This includes, oversampling the class with fewer instances to try and achieve an equal balance between the classes, and under sampling the data to try and reduce the instances in the larger class.However, over-or under sampling can cause biases in the data as either the model is overfit, as in the case of over fitting, or, in the case of under sampling, information is lost.Another technique to improve classification performances when dealing with class imbalances, is boosting.Boosting can improve the performance of weak classifiers, such as decision trees, by re-sampling the training data by its assigned weights.Boosting techniques constantly tweak weak learners until they converge towards a strong learner.

B. SVM Binary Classification from Raw Data
In our application of SVM approach, a hyperplane H is created to divide the classes into their proper labels (-1 or 1) and is designed to give the largest margin between these classes.The hyperplane is then constructed as, w T x i + b = 0, where 1 w is the distance to H + and H − and w is the vector of weights for the features in x.In order to obtain the optimal hyperlpane, 1 w must be maximized, or w must be minimized.Thus, the Maximum Likelihood Estimator (MLE) for SVMs is found by minimizing 1 N N i=1 max{0, 1 − y i w T x i } + λ w 2 , where y is the class label and λ is the penalty on the weights w.

C. RUSBoosting
RUSBoosting is a hybrid method of random undersampling and boosting, which is designed to improve performance when classifying imbalanced data.RUS-Boost randomly removes instances of the majority class to try and equalize it with the minority class.RUSBoost has been shown to work well, as well as, perform simply and speedily.Seiffert et al. [16] demonstrated that RUSBoosting shares similiar performances to SMOTE-Boosting, which is a common algorithm for managing class imbalances, but does so faster and more simply.Given a training set S of examples (x 1 , y 1 )...(x n , y n ) with a minority class y τ ∈ Y, |Y | = 2.A new set D is generated from random undersampling, and a weak learner (i.e. a decision tree) is called to create a hypothesis h t .The pseudo loss, for S and D t , is then calculated using t = (i,y):yi =y D t (i)(1 − h t (x i , y i ) + h t (x i , y i )).
The weight update parameter is then derived with: α t = t 1− t .Normalize this and the final hypothesis H(x) = argmax y∈Y T t=1 h t (x, y)log 1 αt .

D. ICA
We are comparing our methods with the Infomax ICA algorithm which given observed components, in our case the matrix L T , finds a linear transformation A from L T to S, i.e. S = AL T .The matrix S is decomposed into independent components that minimize mutual information.The assumption of the ICA model is that the components in L are statistically independent.We are using the implementation of the infomax based on [4], in which a maximum likelihood estimation is found by using the log-likelihood in the form LH = T t=1 n i=1 logf i (w T i S(t)) + T log|det(A)| where f i is the density functions of the columns of L T .This is equivalent to the Infomax Principle, M 2 = H(φ 1 (w T 1 x), ..., φ n (w T n x)) where φ i are nonlinear scalar functions.To solve the Infomax algorithm we used Python software package [].

Fig. 1 :
Fig. 1: Supervised Classification of Fixed Point Network Responses.Left to right: Distinct input signals, stimuli S i (here i = 2), are injected into network's input nodes and produce fixed point dynamics in output nodes (n1 and n2 are samples of output nodes).For each stimulus output nodes timeseries are recorded in a response matrix, with dimensions of nodes × time, where each node is a row of the matrix.Response matrices to S 1 and S 2 stimuli are shown as color maps and time series of sampled nodes n1 and n2 (rows 1 and 2 of each response matrix) are shown in insets.As the number of output nodes increases, finding a separating hyperplane between the dynamics represented in the nodes space becomes complex as illustrated in the rightmost plot.

Fig. 2 :
Fig. 2: Benchmark data set of supervised recordings and superpositions of stimuli from the olfactory network.Recordings of extracellular neural responses of 106 neurons for 8 individual stimuli (odorants) that constitute the Datura scent, labeled as S1,...,S8.In addition, responses to mixed stimuli labeled as behavioral B1,..,B3, non-behavioral E1,...,E5 and control stimulus E6 were recorded.For each stimulus responses were recorded over 5 distinct trials.

Fig. 3 :
Fig. 3: Possible Methods for Classification From Multiple Node Timeseries Collection Top row: Application of SVM and RUSBoosting classifiers to perform binary classification directly on the collection C in the benchmark set.Precision and recall percentages are computed and shown in the right graph.Middle row: Treatment of the collection as a concatenated matrix, Eq. 1 -SVDCon method.Projection of the responses onto low-dimensional space spanned by P C 1 − P C 3 obtained from SVD on concatenated matrix shown on the right.Bottom row: Treatment of the collection as a set of individual matrices, Eq. 3 -SVDSep method, and projection onto the first PC vector of each matrix SVD shown on the right.

Fig. 4 :
Fig. 4: Exclusive Threshold Reduction (ETR) and Optimization (OETR).A: The exclusive threshold reduction (ETR) method applies a maximizing rule (top left) to each neuron/node in the matrix L and produces the matrix O where L T O defines coordinates of fixed points (bottom left).Projection of distinct stimuli trajectories onto the matrix O resulting in a phase space where each stimulus trajectory is correlated with its associated axis (top right; compare with C).B: Optimal exclusive threshold reduction (OETR) re-weights the nodes by solving a convex optimization problem associated where D w is optimized to find the optimal weighting of L T D w O = I (top left).The coordinates of fixed points are then defined as L T D w O (bottom left).This further separates the transient trajectories and associated fixed points (top right).C: projection onto basis obtained by ICA method.

B 1 :
Fig.5: Classification and Recognition Employing the Classification Space.A: ETR projection of response that corresponds to B1 (blue) and its associated 3D hyperellipse (black sphere) that is used for recognition of B1.All points that fall inside the hyperellipse are marked as "recognition evidence" (s t = 1).For reference we also show the projections of responses to B1 constituents, S1, S2, S3, in gray.B 1 : Classification accuracy into behavioral/nonbehavioral classes for OETR (red), ETR (yellow), ICA (light blue), SVDSep (navy) for increasing dimension of the Classification Space.B 2 : Breakdown of normalized Rec scores per per stimulus (for 17 stimuli) for each method and m = 7.The horizontal blue line indicates the threshold line d for identification whether the response is similar to B1/behavioral.The color of the bars symbolize correctness of classification: correctly classified stimuli are marked by green, and incorrectly classified stimuli are marked by red.C: Precision of the each method when the radius of the classifying hyperellipse (sphere) is varied.D: Recognition precisions using B1 sphere for instantenously identifying B1 stimuli responses (right) and when when false positives from B class are allowed (left).
, different tactics are needed to work with imbalanced data.The first tactic, is to change the performance metric considered.Two informative measures are precision (the positive predictive value, or, what fraction of the labeled class actually are in the class) and recall (the sensitivity, or, what fraction of the relevant class was retrieved).Precision (P ) and recall (R) are defined as follows: P = tp tp+f p , R = tp tp+f n