Data Mining: Machine Learning and Statistical Techniques

The interdisciplinary field of Data Mining (DM) arises from the confluence of statistics and machine learning (artificial intelligence). It provides a technology that helps to analyze and understand the information contained in a database, and it has been used in a large number of fields or applications. Specifically, the concept DM derives from the similarity between the search for valuable information in databases and mining valuable minerals in a mountain. The idea is that the raw material is the data to analyse, and we use a set of learning algorithms acting as diggers to search for valuable nuggets of information (Bigus, 1996). We offer an applied vision of DM techniques, in order to provide a didactic perspective of the data analysis process of these techniques. We analyze and compare the results from applying machine learning algorithms and statistical techniques, under DM methodology, in searching for knowledge models that show the structures and regularities underlying the data analysed. In this sense, some authors have pointed out that DM consists of “the analysis of (often large) observational datasets to find unsuspected relationships and to summarize the data in novel ways that are both understandable and useful to the data owner” (Hand, Mannila & Smyth, 2001), or, more simply, “the search for valuable information in large volumes of data” (Weiss & Indurkhya, 1998), or “the discovery of interesting, unexpected or valuable structures in large databases” (Hand, 2007). Other authors define DM as “the exploration and analysis of large quantities of data in order to discover meaningful patterns and rules” (Berry & Linoff, 2004). These definitions make it clear that DM is an appropriate process for detecting relationships and patterns in large databases (although we point out that it can also be applied in relatively small databases). In this sense, the concept of Knowledge Discovery in Databases (KDD) has been frequently used in the literature to define this process (Han & Kamber, 2000, 2006; Hand et al., 2001), specifying that DM is a stage of the process, and highlighting the need for a previous stage of integration and collection of data (if we start with large raw databases), and also the stage of cleaning and preparing data (data pre-processing) before building descriptive/predictive models in the DM stage (applying suitable techniques to the analysis requirements). On the other hand, several authors have used the concept of DM (instead of KDD) to refer to the complete process (Bigus, 1996; Two Crows, 1999; Paul, Guatam & Balint, 2002; Kantardzic, 2003; Ye, 2003; Larose, 2005).

Methodology
In this section we aim to offer an integrating view of the use of Data Mining (DM) methodology and techniques through a presentation of the procedures common to these techniques in the process of obtaining predictive models, and through a description of the methodological peculiarities associated to them. Specifically, we focus on a description of techniques which allow us to generate categorical response predictive models (classification models). Table 1 presents a classification of some techniques included in DM according to the nature of the data analyzed. In this sense, we present the techniques available depending on the nature of the predictor variables and the output variable. If the output variable is continuous or categorical we find ourselves dealing with supervised learning models, whereas if there is no output variable we are dealing with unsupervised learning models. Specifically, we will deal with the techniques of Neural networks, Classification trees, k-Nearest Neighbor, Naive Bayes and Logistic regression, since they make it possible to analyze categorical output variables in order to generate classification models. Classification is a task that belongs to the category of supervised learning and it refers to the task of analyzing a set of pre-classified objects in order to learn a model (or function) which can be used to classify unknown data in one of several predefined classes . From the perspective of supervised learning, the analysis technique estimates the model from the knowledge that it has of the behaviour of each of the entries in the selected output variable, in such a way that the supervised techniques itself supervises whether or not the model it is building adjusts to the knowledge it has of the reality. In this sense, the aim of Table 1. Data mining techniques according to the nature of the data (Shmueli et al., 2007) supervised learning is to generate knowledge based models which will help in predicting the behaviour of new data. A common requirement in predictive modelling techniques is the use of a data sample (test data) which is independent of the one used in the construction of the model (training data), with the intention of assessing the model's generalization capacity (assessment of the model).
On the other hand, in unsupervised learning there are no known results to guide the algorithm in obtaining the model, but rather this explores the properties of the data with the aim of identifying behaviour patterns with no knowledge of these "a priori". In this way, the aim of unsupervised learning is to generate knowledge based models with a descriptive, not predictive, intention.

Artificial Neural Networks
Artificial Neural Networks (ANN) are data processing systems whose structure and functioning are inspired by biological neural networks. ANN were developed based on the following guidelines: • Information processing occurs in simple elements called neurons.

•
The neurons transmit signals through established connections. • Each connection (communication link) has an associated weight. • Each neuron applies an activation function (usually non linear) to the total entry of connected neurons received (sum of entries weighted according to the connection weights), thus obtaining an output value which will act as the entry value which will be transmitted to the rest of the network. The fundamental characteristics of ANN are parallel processing, distributed memory and adaptability to the surroundings. The processing unit is the artificial neuron, which receives the entries from the neighbouring neurons and calculates an output value, which is sent to all the remaining neurons. As far as the representation of input and output information is concerned, we can find networks with continuous input and output data, networks with discrete or binary input and output data and networks with continuous input data and discrete output data. An ANN is made up of the sequential order of three basic types of nodes or layers: input nodes, output nodes and intermediate nodes (hidden layer) (Figure 1). The input nodes are in charge of receiving the initial values of the data from each case in order to transmit them to the network. The output nodes receive input and calculate the output value. This set of nodes used by the ANN, together with the activation function, makes it possible for the ANN to easily represent non-linear relationships, which are the most difficult as far as multivariate techniques are concerned. The most used activation functions are: the step function, identity function, sigmoid or logistic function and the hyperbolic tangent. There is a large selection of ANN models. A combination of the topology (number of neurons and hidden layers, and how they are connected), the learning paradigm and the learning algorithm define an ANN model (Bigus, 1996). It can be said that an ANN has three advantages which makes it very attractive in data handling: adaptive learning through examples, robustness in handling redundant and inaccurate information and massive parallelism. The most used method in the practical applications of ANN is the multilayer perceptron, which was made popular by Rumelhart et al. (1986). A multilayer perceptron type of ANN starts with an input layer in which each node or neuron corresponds to a predictor variable. These input neurons connect with each of the neurons making up the hidden layer. The nodes in the hidden layer in turn connect with the neurons in another hidden layer. The output layer is made up of one (binary prediction) or more output neurons. In this sort of architecture, the information is always transmitted from the input layer towards the output layer. The popularity of the multilayer perceptron is mainly due to the fact that it is capable of acting as a universal function approximator. More specifically, a "backpropagation" network which contains at least one hidden layer with enough non-linear units can learn any sort of function or continuous relationship between a group of input variables (discrete and/or continuous) and an output variable (discrete or continuous). This property makes multilayer perceptron networks general, flexible and non-linear tools. A complete description of the mathematical foundations associated with the training stage and the functioning stage of the backpropagation algorithm in multilayer perceptron architecture can be found in Rumelhart et al. (1986).
The usefulness of the multilayer perceptron, lies in its ability to learn virtually any relationship between a set of input and output variables. On the other hand, if we use techniques derived from classical statistics such as linear discriminant analysis, this does not have the capacity of calculating non-linear functions and, therefore, will show a lower performance compared to the multilayer perceptron in classification tasks that involve complex non-linear relationships. When the network is used to classify normally the output layer has as many nodes as the number of classes and the node of the output layer with the highest value offers the estimate of the class which the network makes for a certain input. In the special case of two classes it is common to have a node in the output layer, and the classification between the two classes is carried out by applying a cut off point to the node value. If one of the virtues of ANN is that they allow modelling any sort of functional relationship (linear or non linear) between variables and, therefore, act as universal function approximators, another of the outstanding advantages of this technique, compared to classical modelling techniques, is that it does not impose any sort of restriction with respect to the starting data (type of functional relationship between variables), neither does it usually start from specific assumptions (like the type of distribution the data follow). Another virtue of the technique lies in its capacity to estimate good models even despite the existence of noise in the information analyzed, as occurs when there is a presence of omitted values or outlier values in the distribution of the variables. Hence, it is a robust technique when dealing with problems of noise in the information presented; however, this does not mean that the cleaning criteria of the data matrix should be relaxed. Nevertheless, its extreme flexibility lies in the need to have sufficient training data and that it requires more time for its execution than other techniques (Shmueli et al., 2007). It is worth pointing out that in ANN, as well as the set of training data to build the model and the set of independent data (test data) in order to assess its generalization capacity, a third set of independent data (validation set) is used to avoid overfitting the model (during the learning process) which can cause an excessive number of parameters or weights regarding the problem (Hastie et al., 2001, p. 356). Despite the advantages presented concerning the technique, on the other hand, one of the most important criticisms that have been raised against the use of ANN focuses on the fact that a knowledge of the weights within the network does not in general help the interpretation of the underlying process that the prediction of a certain output value generates. To put it another way, the reproaches against the use of this technique are limited to the difficulty in understanding the nature of the internal representations generated by the network in response to a certain problem. Despite this, this perception of ANN as a complex "black box" is not completely true. In this sense, different attempts at interpreting the weights of the neuronal network have arisen, of which the most widely used is the so-called sensitivity analysis (Montaño & Palmer, 2003), implemented in ANN programmes as recently presented by Palmer et al. (2001), under the name of Sensitivity Neural Network 1.0.

Decision Trees
Decision trees (DT) are sequential partitions of a set of data that maximise the differences of a dependent variable (response or output variable). They offer a concise way of defining groups that are consistent in their attributes but which vary in terms of the dependent variable. D T a r e m a d e u p o f n o d e s ( i n p u t v a r i a b l e s ) , branches (groups of entries in the input variables) and leaves or leaf nodes (values of the output variable). The construction of a DT is based on the principle of "divide and conquer": through a supervised learning algorithm, successive divisions of the multivariable space are carried out in order to maximise the distance between groups in each division (that is, carry out partitions that discriminate). The division process finalizes when all the entries of a branch have the same value in the output variable (pure leaf node), giving rise to the complete model (maximum specified). The further down the input variables are in the tree, the less important they are in the output classification (and the less generalization they allow, due to the decrease in the number of inputs in the descending branches).
To avoid overfitting the model, the tree can be pruned by eliminating the branches with few or scarcely significant entries. As a result, if we start from the complete model, after the tree pruning this will gain in generalization capacity (assessed with test data), at the expense of reducing the degree of purity of its leaves (Larose, 2005). There are different learning algorithms designed to obtain DT models (see Table 2). The learning algorithm determines the following aspects: • Specific compatibility with the type of variables: nature of the input variables and the output variable. • Assessment procedure of the distance between groups in each division: division criteria.

•
Restrictions can be placed on the number of branches each node can be divided into.

•
Pruning parameters [pre-pruning / post-pruning]: minimum number of entries per node or branch, critical value of the division, performance difference between the extended and reduced tree. Pre-pruning implies using stopping criteria during the construction of the tree, whereas post-pruning applies the pruning parameters to the whole tree. The most used algorithms (Table 2) are CART (Classification And Regression Trees), CHAID (Chi-Squared Automatic Interaction Detection), QUEST (Quick, Unbiased, Efficient Statistical Tree) and C4.5 / C5.0. The CART algorithm was designed by Breiman et al. (1984) and it generates binary decision trees, where each node is divided exactly into two branches. In this way, if the input variable is nominal and has more than two categories, it groups different categories in one branch. If the input variable is nominal or continuous, it still generates two branches, associating a set of values limited by the operators to each one of them "less than or equal to" or "greater" than a certain value. The CART algorithm makes it possible to introduce nominal, ordinal and continuous input data into the model. The output variable of the model may likewise be nominal, ordinal or continuous. The CHAID algorithm (Kass, 1980) was originally designed to handle only categorical variables. Nevertheless, nowadays it makes it possible to handle nominal and ordinal categorical output data and continuous variables. The tree construction process is based on the calculation of the significance of a statistical contrast as a criterion in order to decide the hierarchy of importance of the predictor variables, and to establish clusters of similar values (statistically homogeneous) with respect to the output variable, keeping all the values that turn out to be heterogeneous (distinct) unaltered. Similar values are melted in one category, forming part of one branch of the tree. The statistical test used depends on the level of measurement of the output variable. If the aforementioned variable is continuous, the F test is used. If the output variable is categorical, the Chi-square test is used.
A differential characteristic between the CART and CHAID algorithms is that the latter allows the division of each node into more than one branch; therefore it tends to create much wider trees than the binary development methods. The QUEST algorithm (Loh & Shih, 1997) can be used if the output is nominal-categorical (allows the creation of classification trees). The tree construction process is also based on the calculation of the significance of a statistical contrast. For each input variable, if this is a nominal categorical variable, it calculates the critical level of a Pearson Chi-square independence contrast between the input variable and the output variable. If the input variable is ordinal or continuous, it uses the F test. The C5.0 algorithm (Quinlan, 1997) only admits categorical output variables. Input variables may be categorical or continuous. This algorithm is the result of the evolution of algorithm C4.5 (Quinlan, 1993) designed by the same author and which has as a nucleus the ID3 version (Quinlan, 1986). The ID3 algorithm is based on the concept of information gain to select the best attribute.  Table 2. Comparative between learning algorithms for decision trees (amplified version of Gervilla et al., 2009) One of the most outstanding advantages of DT is their descriptive nature, which allows us to easily understand and interpret the decisions made by the model, as we have access to the rules that are used in the predictive task (an aspect that is not taken into consideration in other machine learning techniques, such as ANN). Thus, DT allow the graphic representation of a series of rules concernin g t h e d e c i s i o n t h a t m u s t b e m a d e i n t h e assignment of an output value for a certain entry, offering a friendly, intuitive explanation of the results. On the other hand, the decision rules provided by a tree model have a predictive value (not only descriptive) from the moment in which their accuracy is assessed from independent data (test data) to the ones used in the construction of the model (training data). Another attractive characteristic of DT is that they are intrinsically robust to missing values, as they handle them without having to impute values or eliminate observations. Nevertheless, DT do have some weaknesses (Shmueli et al., 2007): they are sensitive to small changes in the data and, unlike the models that assume a particular relationship between the response and the prediction (e.g. a linear relationships like a linear regression), DT are not linear and are not parametrical. This allows for a wide range of relationships between the predictors and the response, but it can be a weakness: given the fact that the partitions are carried out on unique predictors rather than on combinations of predictors, DT probably omit relationships between predictors, particularly linear structures such as linear or logistic regression models. Another drawback in the construction of DT is the problem of overfitting the model, that is to say, the DT includes not only the real patterns or structures present in the data, but also part of the "noise". To reduce this problem as much as possible there are several strategies:

Type of prediction
• Strategies that slow down the growth of the tree before it reaches the perfect classification of the examples in the training set (for instance, the CHAID algorithm).

•
Strategies that make it possible for the tree to grow completely and afterwards carry out some pruning (for instance, the CART and C4.5 algorithms). These latter have been shown to be more efficient than the former. One final disadvantage of DT is that they need a large set of data in order to build a good classifier. Nevertheless, Breiman & Cutler (2004) have introduced "Random Forests", which deal with these limitations. The basic idea is to create multiple DT from the data (and thus obtain the "forest") and to combine their result to obtain a better classifier (see Breiman, 2001).

k-Nearest Neighbor
When we run into new situations, we human beings are guided by memories of similar situations we have experienced in the past. This is the basis of the k-Nearest Neighbor (k-NN) technique. That is, the k-NN technique is based on the concept of similarity. Moreover, this technique constructs a classification method without making assumptions concerning the shape of the function that relates the dependent variable with the independent variables. The aim is to identify in a dynamic way k observations in the training data that are similar to a new observation that we want to classify. In this way, k similar (neighbouring) observations are used to classify the observation specifically in a class (see Figure 2). More specifically, k-NN looks for observations in the training data that are similar or near to the observation that has to be classified, based on the values of the independent variables (attributes). Then, depending on the classes of these nearby observations, it assigns a class to the observation that it wishes to classify, taking the majority vote of the neighbours to determine the class. In other words, it counts the number of cases in each class and assigns the new case to the one that most of its neighbours belongs to (Two Crows, 1999). Even though the method may appear "naive", it is capable of competing with the other more sophisticated classification methods. Hence, where the linear model is rigid, k-NN is extremely flexible. The performance of this technique for data of the same size depends on k, and on the measurement used to determine which observations are nearby. As a result, in the application of the technique we must take into account how many neighbours are to be considered (k value), how we measure the distance, how we combine the information for more than one observation and whether all the neighbours should bear the same weight (see Berry & Linoff, 2004). As has been said, the k-NN technique classifies an unknown example in the most common class by using its k nearby neighbours. It assumes that all the examples correspond to points in an n-dimensional space. A neighbour is considered nearby if it has the least distance in the n-dimensional space of attributes . If we set k=1, the unknown example is classified in the class of the nearest neighbour in the training data. Although there is no formula to choose the k number, it is worth noting that if we choose a small k value, the classification may possibly be too affected by outlier values or unusual observations. On the other hand, choosing a not very small k value will tend to damp any idiosyncratic behaviour learnt from the training data. Nevertheless, if we choose a k value that is too large, locally interesting behaviour will be overlooked. You can let the data help to solve this problem by following a procedure of cross validation. That is, we can try several k values with different training sets chosen at random and choose the k value that minimizes the classification error. With respect to the way of measuring the distance, the most common distance function is the Euclidean distance (1), where x and y represent the m values of the attributes of two cases.
Although, alternatively, the Manhattan distance may also be used (2) The Euclidean distance has three drawbacks: • The distance depends on the units chosen to measure the variables.

•
It does not take into account the variability of the different variables.

•
It ignores the correlation between variables. One solution is to use a measurement called statistical distance (or Mahalanobis distance). By way of advantages of the k-NN technique we can highlight Mitchell, 1997;Nisbet et al., 2009): first of all, it does not simplify the distribution of objects in space in a set of comprehensible characteristics; instead, the training set is stored completely as a description of this distribution. Furthermore, the k-NN method is intuitive, easy to implement and effective in practice. It can construct a different approximation to the target function for each new example to be classified, which is advantageous when the target function is very complex, but it can be described by a collection of less complex local approximations. Finally, this technique builds a classification method without carrying out assumptions concerning the shape of the function that relates the dependent variable (classification variable) with the independent variables (attributes).
On the other hand, the most important disadvantage is that k-NN is very sensitive to the presence of irrelevant parameters. Other disadvantages: • The time to find nearby neighbours in a large training set may be very high.

•
The number of observations that are needed in the training data increases exponentially with the number of dimensions (variables).

Naive Bayes
Bayesian methods use the Bayes rule or formula (3) (based on Bayes' theorem), which expresses a very powerful framework to combine information from the sample with expert opinion (prior probability) so as to produce an up-dated expert opinion (posterior probability) (Giudici, 2003). Specifically, the Naive Bayes technique (NB) is a very powerful classification technique and is one of the most widely used ones, due to its computationally simple process (Hand & Yu, 2001). As it is based on Bayes' theorem, it can predict the probability of a given case belonging to a certain class. Its computational simplicity is due to the assumption known as class conditional independence (this assumes that the effect of an attribute value on a certain class is independent of the values of the other attributes), and in this sense it is considered a "naive" classifier (Han & Kamber, 2000. This classifier predicts that a case A will belong to class C i which has the highest X conditioned posterior probability (set of attributes of the case in the predictor variables). Bayes' theorem allows us to define the Bayes' formula (3) that this posterior probability provides; and since P(X) is constant for all the classes, we only need to maximize P(X|C i )P(C i ) in the classification process.
From a set of training data, P(C i ) how many times each class C i occurs in these data can be estimated. To reduce the computational cost of estimating P(X|C i ) for all possible x k (predictor variables), the classifier uses precisely the "naive" assumption that the attributes used to describe X are conditionally independent from each other given class C i . This conditional independence can be found in the expression (4), where the m value indicates the number of predictor variables that participate in the classification.
The studies that compare classification algorithms (e.g. Michie et al., 1994) have often shown that NB is comparable in its working with ANN and DT, and in fact exceeds these sophisticated classifiers if the attributes are conditionally independent given the class. Recent theoretical analyses have shown why NB is so robust (Domingos & Pazzani, 1997;Rish, 2001). The appeal of NB classifiers lies in their simplicity, computational efficiency and good performance in classification. What is more, NB can easily handle unknown values or missing values. Nevertheless, it has three important drawbacks (Shmueli et al., 2007): first of all, it requires a large number of cases in order to obtain good results; secondly, if a prediction category is not present in the training data, the technique assumes that a new case with this category in the predictor has zero probability; finally, even though we obtain a good performance if the aim is to classify or order cases according to their probability of belonging to a certain class, this method offers very biased results when the aim is to estimate the probability of belonging to a class. Above and beyond these drawbacks, the NB technique is straightforward to use, it adjusts to the data and is easy to interpret. In addition, it requires only one exploration of the data. This simplicity, parsimony and interpretability has led it to enjoy widespread popularity, especially in the literature of machine learning .

Logistic Regression
Linear regression is used to approach the relationship between a continuous response variable and a set of predictor variables. However, when the response variable is categorical, linear regression is not appropriate. Logistic regression (LR) is a generalized linear model. It is used mainly to predict binary variables (with values like yes/no or 0/1). Thus, LR techniques may be used to classify a new observation, whose group is unknown, in one of the groups, based on the values of the predictor variables. As with linear regression, the classification depends on the linear combination of the attributes. The logistic function (5) In LR we follow two steps: the first step consists of estimating the probability of belonging to each group and in the second step we use a cut off point with these probabilities to classify each case in one of the groups. The parameters of the model are estimated with the method of maximum likelihood through a process of successive iterations. For a more detailed explanation of the LR technique, consult Larose (2006). Lastly, it is worth commenting that LR can produce stable results with relatively few data (Harris-Jones & Haines, 1998). On the other hand, the fact that traditional regression is so widely accepted, easily implemented and generally understood makes it even more attractive (Hill et al., 2004). What is more, LR shows behaviour analogous to a diagnostic test.

Applying data mining techniques
In this section we go further, from the applied point of view, into the integration within DM methodology of the techniques described in the previous section. We use the Weka platform to meet this aim. In the first section, we give a brief description of the functionalities integrated in the Weka interface, with the aim of providing the reader with a presentation of the tool and to specify its virtues. In the second section, we propose a case study in order to compare with Weka the performance of predictive models obtained with the techniques indicated.

Weka interface
Weka (Waikato Environment for Knowledge Analysis) is a data mining platform distributed under public license GNU-GPL: it is free software that can be freely used, copied, studied, modified and distributed and it is protected from appropriation attempts that would restrict these user liberties. Bearing in mind its characteristics, we find that it is a tool which, first of all, has an interactive interface which contains four user-machine interaction modalities (Fig. 3): • Explorer: is the most used mode and the most descriptive.  Within the Explorer modality, it provides extensive support to the overall process of data mining (Fig. 4): 1. Access to databases, exploration and selection of data and data processing: • Preprocess: functionality aimed at importation, transformation (application of filters) and data extraction. • Visualize: functionality aimed at the visualization of data using graphic techniques. 2. Predictive and descriptive modelling. Compiles a wide range of data mining procedures for the obtention of knowledge models: • Classify (classification and regression): predictive modelling (supervised learning). • Cluster (grouping) and Associate (association rules): descriptive modelling (unsupervised learning). • Select attributes: selection of predictive attributes. Lastly, it is worth highlighting the possibility of system extensibility: it allows the user to modify Weka by integrating new functionality developed in Java code, using its structure and object oriented functional design. This represents the main advantage as opposed to other closed code data mining platforms (commercial programmes). In Witten & Frank (2005) you can find a detailed description of the different modalities of interaction with Weka.

Case study
Once we have presented the methodological basis of the techniques involved in this work, in this section we aim to contribute comparative elements of the information provided by these techniques in an applied context. These comparative elements refer to the predictive power (accuracy) of the knowledge models generated and the descriptive components that add informative value to the decision making processes (classifications). Thus we aim to provide a more integrating view, if possible, of DM methodology, since we provide common assessment parameters in order to compare the results obtained. Nevertheless, from the analysis of these results we do not aim to reach substantial conclusions related to the context from which the data used comes, but rather our intention is to divulge to the readers a series of methodological tools that allow us to detect knowledge patterns in a way that is practically automatic and, on the other hand, to make it easier to interpret the descriptive elements associated with the assessment of the models obtained. From an initial sample of 9300 young people aged between 14 and 18 years, in which information concerning variables intervening in the consumption of addictive substances was collected, we selected a sub-sample of 2526 young people. We are interested in studying the relationship between the consumption / non consumption of cannabis (output variable) among the people surveyed and the reasons the subjects surveyed have for consuming or not consuming drugs (input variables). Specifically, we collected fifteen possible reasons (variables) for consuming drugs and eleven reasons (variables) for not consuming; the possible response to each of these variables is dichotomous (yes/no). On the other hand, if in the initial sample (complete) the percentage of consumers of cannabis is nearly 18%, as opposed to 78% of non consumers, the sub-sample selected shows a greater balance between consumers/non consumers (44.4%/55.6%); this equilibrium (or balancing) is justified by methodological motives, as there must be a similar number of entries in each of the output categories (consumes/does not consume), so that they can be equally represented in the modelling stage (detection of classificatory patterns). In a session with Explorer, first of all we loaded (Open file…) in the Preprocess section the data to be analyzed, whose structure (database) has been adapted to the Weka format: Arff file. Once the data file is open, it is possible to explore the variables they contain (Fig. 5). It is also possible to read data in the CSV (comma delimited) format from Weka, although it is not possible to import databases in the more widely used formats such as Excel, Access, SPSS, etc. However, there is the possibility of converting these other more common formats into the native Arff format from the data mining platform RapidMiner (free, open code software, which also allows the use of the algorithms included in Weka). For instance, if the data source is in SPSS format, we can indicate whether we are interested in extracting the names of the variables and/or their labels in another format (Arff, in this case), and whether we are interested in using the labels of the values (option be default) instead of the numerical values. From the section Classify, we access the different modelling techniques integrated in Weka (Fig. 6). For instance, we can select the J48 classifier (weka.classifiers.trees.J48); this classifier uses the C4.5 algorithm (Quinlan, 1993) to generate a classification tree which is in agreement with a series of parameters determined by the algorithm (to edit them, click on the classifier) and other parameters determined by data mining methodology (in Test options). In the example (Fig. 6) we have indicated that the J48 classifier uses 70%of the sample (training data) to create the model, and the rest as test data. The output variable is also indicated (predicted variable), which by default is the last variable in the database. The Start button allows us to generate the model and access (in Classifier output) the model's assessment results. It can be observed that the model has correctly classified 599 of the 758 test patterns (79%), with a larger percentage of hits in the category Does not consume (83.4%) than in the category Consumes (73.9%). It is possible to access the graphic representation of the classification tree (Fig. 7) through the options of the contextual menu (Visualize tree) of the model generated (in the Result list).

Fig. 6. Evaluation of the selected decision tree model on test data
What is more, the descriptive nature of decision trees allows us to easily derive the classification rules used by the model, with information from their support (percentage of cases in which the antecedent and consequent of the rule -correct prediction -are found with respect to the total subjects) and their confidence (percentage of cases in which the antecedent and consequent of the rule are found with respect to the number of subjects in which the antecedent of the rule is found). In Table 3 we show some examples. If we focus on the rule associated to leaf 1 (Leaf1 in Fig. 7), it can be seen that out of the 1768 subjects who took part in the generation of the model (training stage), 550 (31.3%) of them considerantecedents of the rule -their pleasurable nature a reason for consuming drugs (rcons3) and do not consider the fact that their friends consume a reason for consuming (rcons9) and at the same time -consequent of the rule -they are consumers of cannabis. Hence, if we want to classify new subjects in the output variable (consumes / does not consume), in the case in which their values in the output variables coincide with the model's rule, this would indicate that they are consumers of cannabis, with a confidence of 82.3% in the decision adopted. The confidence of the rule indicates to us, therefore, that there are 17.7% of subjects (118 cases) who fulfil the antecedent of the rule, but not the consequent (they are not consumers of cannabis).  Table 3. Some classification rules from the selected decision tree As mentioned above, with Weka we have generated other classifying models from the techniques described in the previous section, using the same training (70%) and test (30%) sub-samples in all of them. Specifically, to generate a neural network we used the MultilayerPerceptron classifier (weka.classifiers.functions.MultilayerPerceptron) (it uses the backpropagation algorithm), keeping 20% of the training data as validation data (validationSetSize parameter of the classifier). In order to obtain the k-Nearest Neighbor model, we used the IBk classifier (weka.classifiers.lazy.IBk), with the Euclidean distance function. Lastly, we generated a Naive Bayes model (weka.classifiers.bayes.NaiveBayes) and a logistic regression model (weka.classifiers.functions.Logistic). Table 4 shows the comparison of the correct classifications (accuracy) between models, through the respective confusion matrices (bivariate contingency tables in which the real classification categories are crossed with the entries in the categories estimated by the model). It is shown that in all the models selected there is a relatively high accuracy (ranging between 76.6% for the IBk model and 80.1% for the Logistic model). If we analyze the accuracy according to the classificatory category, in all the models the percentage of correct classifications in the category Does not consume is greater (ranging between 81.4% for the NaiveBayes model and 90.5% for the IBk model), whereas in the Consumes category the range of correct classifications moves between 60.5% for the IBk model and 75.9% for the Logistic model. Since the same test sub-sample was used in all the models selected, it is possible to carry out a comparison of the coincidences in the classification of the cases (  We can also analyze the confusion matrix between the cases that converge on the classification predicted by the five models (576 subjects) and the real classification category (  Table 6. Confusion matrix between prediction agreement and actual classification

www.intechopen.com
Weka also allows us to access the individual prediction data if the option Output predictions is activated using the button More options… (see Fig. 6); these predictions are included in the window of results (Classifier output) together with the model's assessment information. In this way, we can compare the predicted classification and confidence level (probability) in this classification, case by case and for each model. We extracted the individual predictions obtained in each of the five models to a data base outside Weka (Fig. 8), so as to be able to compare the degree of agreement and the confidence level in the joint classificatory decision (0: does not consume; 1: consumes) for a certain case. For instance, we can compare the predicted classification in cases 1 and 6; in both cases, the five models converge on the classification (consumes), even though the confidence level for this decision is lower (on average) in subject 1 (p=0.805) than in subject 6 (p=0.925). Weka also allows us to reassess a given model with a new set of data independent from the one used in its construction. Through the contextual menu of the Result list section of Explorer (see Fig. 6), the model can be loaded in memory (Load model option) if it had been saved (Save model option) before closing the Weka session in which it was built. Once the model has been loaded in the memory, a new set of data must be chosen through the Supplied test set option, and finally select Re-evaluate model on current test set in the contextual menu of the model. This option can also serve in new cases to find out the classification proposal of already validated predictive models (by activating the Output predictions option). Finally, it is possible to access the Weka Experimenter model in order to configure experiments (simulations) with different predictive model techniques applied to the same set of data (Fig. 9), and thus be able to assess and compare the capacity of the techniques to generate good predictive models. We proposed an experiment to compare the five classificatory techniques under study in relation to the ZeroR classifier (weka.classifiers.rules.ZeroR). This is the simplest classifier, which predicts the mode (the most repeated value) for a categorical output variable (and the mean for a numerical output variable). In the experiment, we indicate the division of the sample into training (70%) and test (30%) data and 10 repetitions (models) were specified for each technique.
The comparison was carried out from the mean percentage values of the cases correctly classified in the ten models generated for each technique and their standard deviations (hidden by default, but which can be shown if the corresponding option is chosen). As well as the percentage of correct cases (Percent_correct) it is possible to choose other indices of interest in the Comparison field (see Fig. 9), such as the area under ROC curve (Area_under_ROC). The values of these indices for each of the specified repetitions are stored in an Arff data file, whose variables can be analyzed in Weka itself or in other data analysis programmes (with prior extraction to a compatible data format). By pressing the button Experiment (Fig. 9), the mean values and standard deviations of the chosen comparison index are used to study (from a statistical t test) whether there are significant differences between the classificatory accuracy of the models obtained from the first of the techniques defined in the experiment (in this case, ZeroR) and the classificatory accuracy of each of the remaining techniques used. As can be observed, the ZeroR classifier indicates us that the mean percentage of correct classifications (with test data) is 55.6%, which obviously corresponds to the percentage of subjects in the sample who are non consumers (most frequent category). The rest of the classifiers show a much higher mean percentage of classifications, whose difference with respect to this reference value is, in all cases, statistically significant.

Discussion
When faced with the question, what is the best algorithm for classification? There is evidently no general answer that can help us to know prior to an analysis of the data which technique or algorithm I should apply in order to obtain the best classificatory model. In this sense, Nisbet et al. (2009, p. 256) indicate us that if different classificatory algorithms are used, we will discover that the best algorithm for classifying a set of data may not work well in another set of data; in other words, different techniques or algorithms have a better functioning in different data sets, and in this sense, they claim that "using a diversity of algorithms is best". They even establish an analogy between the process of creation of knowledge models and the process of sculpting a statue, which they call "the art of data mining" (Nisbet et al., 2009, p. 46). In the literature you can find DM definitions which point in this same direction, for instance: "data mining is the art of discovering meaningful patterns in data" (Pyle, 2003). Weiss & Indurkhya (1998, p. 21), reflect on this question and ask the question, "Data Mining: Art or Science?", which they answer, "no universal best approach is describable for data mining; making good decisions is part art, part science"; in this sense, these authors combine the art and science of DM in their work: they use science when it is known and effective, and offer guidance in practical issues that are not easily quantifiable.
Our scientific contribution to DM has been precisely to present arguments that point in the direction of convincing applied researchers of the usefulness of using different techniques or algorithms in the search for knowledge models that will help in decision making concerning a given problem. Researchers should take on the role of designer in this task: they should design the data selection, cleaning and preparation processes, as well as the model obtention and validation processes through the wide repertory of associated techniques, algorithms and parameters that are at their disposal. Precisely, the Weka platform offers an ideal space to combine the art and science of DM in an effective way, as demonstrated in the previous section of the chapter.


extfieldsofstudy:
Computer Science
Medicine

s2fieldsofstudy:
Computer Science
Biology

text:
Emergent Properties of Interacting Populations of Spiking Neurons

Dynamic neuronal networks are a key paradigm of increasing importance in brain research, concerned with the functional analysis of biological neuronal networks and, at the same time, with the synthesis of artificial brain-like systems. In this context, neuronal network models serve as mathematical tools to understand the function of brains, but they might as well develop into future tools for enhancing certain functions of our nervous system. Here, we present and discuss our recent achievements in developing multiplicative point processes into a viable mathematical framework for spiking network modeling. The perspective is that the dynamic behavior of these neuronal networks is faithfully reflected by a set of non-linear rate equations, describing all interactions on the population level. These equations are similar in structure to Lotka-Volterra equations, well known by their use in modeling predator-prey relations in population biology, but abundant applications to economic theory have also been described. We present a number of biologically relevant examples for spiking network function, which can be studied with the help of the aforementioned correspondence between spike trains and specific systems of non-linear coupled ordinary differential equations. We claim that, enabled by the use of multiplicative point processes, we can make essential contributions to a more thorough understanding of the dynamical properties of interacting neuronal populations.

INTRODUCTION
Dynamic neuronal networks represent an important new paradigm in neuroscience. The field has increasing impact in fundamental brain research, and it is relevant for some new branches of neural engineering. It is concerned with the analysis of biological neuronal networks and, at the same time, with the synthesis of artificial brain-like systems. Neuronal network models serve as mathematical tools to understand the function of brains, and they may develop into future tools to enhance brain function. Both analysis and synthesis of neuronal networks rely on a thorough understanding of the dynamical properties of neural populations. Currently, this understanding is episodic and elusive, and even the most basic questions about non-linear networks impose serious mathematical challenges.
In a previous paper  we introduced a framework that allows the concise analysis of dynamic activity patterns in structured networks of spiking neurons. Here, we exploit the new theory to analyze small meta-networks of neural populations ("networks of networks"), with a particular focus on functional properties that are considered relevant for biological brains.
Neural spike trains exhibit fluctuations of random appearance, which span different spatial and temporal scales, and which have both intrinsic and extrinsic causes (Marom, 2010). On the one hand, isolated neurons fire with limited reliability, depending on the circumstances of their stimulation (Mainen and Sejnowski, 1995). On the other hand, ongoing activity of unknown origin strongly modulates the activity of cortical networks that are part of an active brain (Arieli et al., 1996). It was also shown theoretically that balanced excitation and inhibition in large recurrent networks is bound to generate a very complex dynamics, resulting in strong intrinsic fluctuations of activity (van Vreeswijk and Sompolinsky, 1998;Jahnke et al., 2009). All this suggests the use of stochastic methods to model neuronal populations. In turn, it is possible that biological networks exploit such random effects to boost their computational properties, as, e.g., in stochastic resonance (Wiesenfeld and Moss, 1995;McDonnell and Abbott, 2009).
A useful framework to define and study stochastic dynamics for spiking neurons and networks is that of a Wiener cascade (Wiener, 1958;Rotter, 1996;Herz et al., 2006), which is related to the concept of escape noise or threshold noise (Plesser and Gerstner, 2000). Cascade models are characterized by linear input integration, and a non-linear transfer function, which specifies how the linear state variable (typically the membrane potential) is mapped onto the instantaneous firing rate of the neuron, see Figure 1.
For exponential transfer functions, specifically, networks of cascade units have a predictable mean-field dynamics . The corresponding coupled non-linear differential equations for the ensemble firing rates turn out to be a competitive system of Lotka-Volterra type. A similar approach has been developed by Buice et al. (2010) starting from certain field equations (Wilson and Cowan, 1972). The long history of FIGURE 1 | Schematic of cascade models. In a neural model with escape noise, the state variable (membrane potential) V is mapped to the instantaneous firing rate r via the transfer function f. The firing rate is used to generate a Poissonian spike train. Thus, the probability of firing a spike in a short interval (t, t + ε) is given by f (V (t ))ε. Escape noise models have the advantage that the spike generation itself is stochastic, making the statistical analysis of such neurons more accessible than deterministic integrate-and-fire neurons with additive noise.
Lotka-Volterra equations in the mathematical literature, starting with Lotka (1910) and Volterra (1926), opens very interesting possibilities for the applications pursued here. In particular, it is possible to reversely engineer any prescribed Lotka-Volterra type competitive system, and to specify the connectivity pattern of the network of neuronal populations that have the given rate equations as a mean-field description. The suggested framework has two main advantages as compared to other models of networks of spiking neurons.
The first advantage is of practical nature: the essential behavior of spiking networks can be studied mathematically by analyzing the associated rate equations. Surprisingly, the rate equations provide an accurate description of the population dynamics not only in the case where the system converges to an attractive fixed point, but also if the rates undergo highly transient dynamics, even in the chaotic case. An important criterion for our choice of the examples was to demonstrate this feature in different situations. Moreover, the known mathematical properties of the non-linear rate equations can be exploited to construct networks of neural populations with prescribed behavior. In the sequel, we will give several examples for this modeling strategy. The second criterion for the choice of the examples was their relevance for better understanding biological spiking networks in the brain, and for designing artificial spiking networks in the emerging field of neural engineering, in which Lotka-Volterra equations already found application (Asai et al., 1999).
The second advantage offered by our framework is more of conceptual nature. Multiplicatively interacting point processes, i.e., cascade models with exponential transfer function, offer a framework to easily implement so-called "scaled-rate" models (Marom, 2010) for networks of spiking neurons. Marom suggests that One good reason to bother with formulation of an abstract model is the hope that it will lead up to a mathematical construct that dramatically reduces the dimensionality of the problem at hand.
Multiplicatively interacting point processes offer an amazingly rich behavior already in low-dimensional cases, thus providing experimentalists and theoreticians with a powerful and innovative set of tools.

COMPARISON TO EXISTING POPULATION RATE MODELS
Rate-based descriptions of neural networks have a long tradition in neuroscience. To embed our work into this context we briefly review the most important approaches.
Temporal coarse-graining was introduced by Wilson and Cowan (1972) in their seminal paper. In this approach, the dynamic variables are given by the expected number of nonrefractory cells receiving input above threshold in a (small) region. Averaging over times shorter than the time-scale of the population dynamics (which is assumed to be in the range of the refractory period) leads to non-linear rate equations. These equations take the form of a linear equation, complemented by a sigmoidal coupling function. These equations have recently been generalized to take fluctuations and higher-order moments into account (Buice and Cowan, 2007;Buice et al., 2010), as well as time-dependent inputs (Ledoux and Brunel, 2011). Already in the original paper by Wilson and Cowan arguments have been given that guarantee that fast dynamics on the time-scale of the refractory period does not affect the population dynamics. Different rate equations have been derived by resorting to ensemble averaging rather than temporal coarse-graining for the case of spiking neurons with refractoriness (Deger et al., 2010). Indeed, the refractory period was found to have a strong effect on the population dynamics.
A very successful, and by now classical theory for studying rate dynamics in homogeneous networks of integrate-and-fire neurons is based on a Fokker-Planck approach (Abbott and van Vreeswijk, 1993;Fusi and Mattia, 1999;Brunel, 2000;Knight et al., 2000). Some results about correlations arising from network interactions have recently been obtained, though (Ostojic et al., 2009). We will not follow this approach here, since the use of second-order partial differential equations makes the study of network dynamics difficult.
Ensemble averaging is the preferred approach for studying time-evolution of activity in neural networks. Ensemble averaging can be used to study neural dynamics both on the single cell and on the population level (Kriener et al., 2008;Deger et al., 2010;Pernice et al., 2011). In this manuscript, in particular, we assume that the instantaneous firing rate of a neuron is a function of its membrane potential, and we use ensemble averaging to obtain the average firing rate of a homogeneous population of independent units. For the transfer function, a common choice for population models is a threshold-linear non-linearity (Wilson and Cowan, 1972;Ledoux and Brunel, 2011). Such models have been used, for example, for the investigation of orientation selectivity in the early visual system (Carandini and Ringach, 1997;Ernst et al., 2001). Recently, a theory for the analytic derivation of transfer functions in the case of leaky integrate-and-fire neurons has been proposed (Ostojic and Brunel, 2011). In our manuscript, we chose an exponential transfer function (Jolivet et al., 2006). Both the reset and refractoriness are incorporated in the self-inhibition of neurons. The rate equations arising from this choice are of Lotka-Volterra type and have already been used in modeling neural activity (McCarley and Hobson, 1975; Frontiers in Computational Neuroscience www.frontiersin.org and Tanaka, 1997;Billock et al., 2001;Rabinovich et al., 2006;Cardanobile and Rotter, 2010). In some cases, Lotka-Volterra equations display the same attractor landscape as linear integrators complemented by a threshold-linear non-linearity (Hahnloser et al., 2003). An important feature of our method is the correspondence between stability properties of the rate equation and of the underlying point process model. This correspondence goes far beyond the matching of attractive fixed points, and includes more complex attractors and associated bifurcations. The main aim of the present work is to demonstrate the correspondence between the point process model and the associated rate equations in different situations. We study feed-forward networks, but also recurrent networks both in the case of single (attractive or repulsive) fixed points, as well as in the multistable case. Finally, we show that stable limit cycles are also correctly mapped from the rate equation to the point process model.

RATE EQUATIONS
Here, we briefly recapitulate the model introduced in Cardanobile and Rotter (2010). In this model, a neuronal population is characterized by its mean membrane potential V (t ). The mean membrane potential is assumed to obey a perfect integrator dynamics (Tuckwell, 1989). (1) represent the spike train of population j, and α ij the net synaptic coupling from neurons in the population j to neurons in the population i. At this step, we have assumed that the populations are statistically homogeneous in the sense that the distribution of synaptic strengths from neurons of population j to neurons in population i is narrow. Action potentials are emitted according to a stochastic mechanism, as in escape noise models for single neuron dynamics (Gerstner and Kistler, 2002).
Such models rely on the assumption that the instantaneous firing rate λ(t ) is only a function of the instantaneous value of the internal state variable V (t ), which is the mean membrane potential in our case.
The interpretation of escape noise models is the following: the membrane potential as well as the spike threshold are subject to spontaneous fluctuations on very short time scales. These fluctuations are both due to internal noise (on the level of synapses and on the level of membrane) and to fluctuating input. Therefore, the probability of observing a spike in a small interval (t, t + t ) is a non-linear function of the membrane potential V (t ) at time t. This transfer function depends on the state of both the neuron and the network. A good choice for the transfer function for neurons is an exponential function (Rotter, 1996;Carandini, 2004;Jolivet et al., 2006). We have verified through simulations (data not shown) that this choice is also suitable for populations of integrate-and-fire neurons.
Under these assumptions, we can write, up to normalization constants Inserting Eq. 2 into Eq. 1, rearranging, taking the ensemble mean, and ignoring covariances leads to the following non-linear dynamic equationṡ (3) We refer to them as the rate equations of the system. They reflect the ensemble behavior of a homogeneous population of spiking neurons. We will use these equations exactly with this interpretation in mind. We would like to stress that the same type of dynamic equations has previously been employed in the neuroscientific literature (McCarley and Hobson, 1975;Fukai and Tanaka, 1997;Billock et al., 2001;Rabinovich et al., 2006), although without the biological motivation of the model offered here. We use the term homogeneous loosely, without a specific statistical framework in mind. Although the exact matching of spiking units to rate equations is a very interesting issue, the focus of this paper is to show that the rate equations offer a precise description of the expected behavior of multiplicatively interacting point process also in the time-dependent regime and to exploit the possibilities offered by this description.
In our firing rate equation, we have disregarded the leak term in the voltage dynamics. The question of whether leak terms can be consistently included into the equation is subject of current research.
The focus of the present work is to demonstrate with neuroscientifically relevant examples that the mapping from the spiking model to Lotka-Volterra equations holds for all stable regimes. The examples included here are multistable systems, and systems with stable limit cycles. We have verified the mapping also for chaotic attractors and Hopf bifurcations, but we have excluded these examples from the present manuscript for the sake of readability. The arguments given in this and our previous manuscript, and the examples of microcircuit design discussed here demonstrate unambiguously that the framework of Lotka-Volterra equations can indeed establish a solid connection between spiking dynamics of neural networks and their mean-field description.

MONTE CARLO SIMULATIONS
All simulations of networks of spiking neurons were implemented in the programming language Python (van Rossum, 1995), the scripts are available upon request. We employed time-driven solvers, based on a fixed step size. Time steps were chosen between 0.5 and 5 ms, depending on the expected spike rates. The goal was to keep the probability of missing a spike as small as possible. In all cases, we have checked that our results are robust against a further decrease of step size (data not shown).
For the simulations, we resorted to a natural reformulation of the model in discrete time. Let S i (t ) again denote the spike train of neuron i, now with the convention that S i (t ) = 1 if neuron i fired a spike at time t, and S i (t ) = 0 otherwise. Due to the exponential transfer function, firing rates couple multiplicatively to incoming spikes such that where dt is the step size of the simulation. In each time step (i) a spike in neuron i is generated with probability 1 − e −λ i (t )dt , and then (ii) all spikes are propagated, and the rates are updated according to Eq. 4.
Simulations are quite effective, since only the state variables of those neurons must be updated that receive spikes. In addition, Python supports vectorized computations of using the numpy (Oliphant, 2006) module. For the for the estimation of the probability of observing an unstable network (Figure 4) and for the simulation of the spiking central pattern generator (Figure 5) we used an event-based solver for increased precision, based on the Doob-Gillespie algorithm (Doob, 1945;Gillespie, 1976).

POWER-LAW RELAXATION IN FEED-FORWARD SYSTEMS
A specific strength of our approach is the possibility of determining the stationary response rate analytically, for a constant stimulus. To appreciate this result, it is instructive to first look at a very elementary "network," composed of one input population and one output population (Figure 4). We assume that the input population fires at a constant rate λ in ≥ 0. The strength of the connection is α syn , which can be either positive (excitation) or negative (inhibition). We further assume that the output population is self-inhibiting with strength α self < 0, implementing some sort of membrane potential reset after having produced a spike. This population model describes a set of unconnected, independently firing neurons. Equivalently, it can be interpreted as the expected behavior of a single neuron, inferred from multiple independent observations (trials). This arrangement typically leads to an output spike train that is slightly more regular than a Poisson process (Softky and Koch, 1993;Shadlen and Newsome, 1998), and that has weak negative serial correlations between adjacent inter-spike intervals (Nawrot et al., 2008).
The system (3) of rate equations simplifies here tȯ This non-linear ordinary differential equation governs the evolution of λ out (t ). Some elementary conclusions can be drawn right away: the output rate λ out (t ) remains always non-negative, provided the initial value is chosen non-negative. The two fixed points, which characterize the stationary solutions of Eq. 5, are 0 and λ * = − α syn α self λ in . It is easy to verify that the more positive of the two fixed points is always globally attractive, and the other one is always globally repelling. That is, arbitrary non-negative initial output rates relax back to max {0, λ * }. Considering α syn λ in as the total signed input, a threshold-linear output characteristic ("halfwave rectifier") results from this behavior, cf. also (Wilson and Cowan, 1972).
Apart from this simple non-linear response to stationary inputs, a feed-forward multiplicative unit also exhibits a power-law transient in its relaxation behavior, see Figure 2. Power-law relaxation to equilibrium has been described in real neurons (Lundstrom et al., 2008). In this paper fractional differentiation was used to model neuronal integration properties, with the aim describe slow time-scale dynamic properties.

INPUT-SELECTIVE RESPONSE
If a unit has more than one input population, its total signed input is the sum of all input rates, weighted by the corresponding coupling constants. This suggests that the equilibrium dynamics of a multi-layer feed-forward network can be obtained by the following simple step-by-step procedure: first, one computes the equilibrium firing rates of the populations at the lowest level, induced by their respective inputs. The equilibrium rates of populations at higher levels are then obtained by the same method, treating the activity of all lower level populations as input.
As an example to illustrate the potential of such simple feedforward networks, we demonstrate how to implement a system that is selective for inputs of strength within a prescribed range. This can be achieved with an arrangement of four populations, as depicted in Figure 3. Following the algorithm explained above, the output rate is given by where we used the notation [λ] + : = max{λ,0}. Furthermore, λ I is a constant level of inhibition, and λ in is the variable excitatory input.
For very strong input λ in , the total excitation to the output unit is given only by the paths including the input, i.e., by (α 41 − α 31 α 43 α 33 )λ in . So, if the term in the brackets is positive, the output activity will increase for increasing input, whereas if it is negative, the output activity will be 0 for strong inputs. Since we want to construct a network selective for a certain range of inputs, we impose that such that the output vanishes for high input rates. On the other hand, given that the input vanishes for high rates, the output is maximal at the point where the input rate to unit 3 overcomes the global inhibitory signal to unit 3, i.e., when Inserting this value in Eq. 6 we obtain a maximal output of We now assume that λ max out > 0. Since α 31 > α 32 then λ max out > 0 if and only if The only negative term in the sum is −α 42 α 32 . So we obtain a bound on the global inhibition level given by Under the conditions Eq. 7 and Eq. 8, the network is sensitive for rates in a certain range, given by the onset of the two threshold-linear functions. An example for a specific set of network parameters is given in Figure 3.
All computations in this network are performed on the level of mean firing rates of the neurons. Therefore, it is important to assess the behavior of the corresponding spiking network, e.g., its performance in a discrimination task. Results are given in Figure 3 and show, using ROC analysis, that different inputs can be well distinguished based on spike counts that are extracted from observations of finite duration.

EXCITATION-INHIBITION INTERPLAY
The situation for a recurrent network is typically more involved. The attractive regimes are not necessarily corresponding to stable fixed points. Depending on the number interacting units, competitive Lotka-Volterra systems can display a variety of different phenomena: (i) continuous manifold of periodic trajectories in the classical two-dimensional case (Arnol'd, 1984); (ii) periodic orbits with increasing periods, single and multiple limit cycles in three-dimensions (May and Leonard, 1975;Hofbauer and So, 1994;Rabinovich et al., 2006); (iii) genuine chaotic dynamics in four dimensions (Wang and Xiao, 2010).
Here, we focus on one specific question of both theoretical and practical relevance: do the stability properties predicted by the Frontiers in Computational Neuroscience www.frontiersin.org rate Eq. 3 also determine the stability of the spiking system? Similar investigations have, in fact, been carried out for networks of integrate-and-fire neurons (Brunel, 2000;Mattia and Del Giudice, 2002;Ledoux and Brunel, 2011) as well as for binary neurons Sompolinsky, 1996, 1998). We obtain a condition similar to the saddle-node bifurcation caused by a rate instability in LIF networks (Mattia and Del Giudice, 2004;Ledoux and Brunel, 2011).
To start our investigation, we recall that a sufficient condition for the dynamic stability of a non-linear equation can be derived from the existence of a so-called Lyapunov function for the system. Under certain conditions, such a function can be constructed for our networks. Specifically, if A r denotes the recurrent part of the coupling matrix, then the system is stable, if all eigenvalues of A r have a negative real part. With a slight abuse of terminology we will refer to such (not necessarily symmetric) matrices as negative definite. In this case, the total firing rate of all populations taken together is a Lyapunov function of the system, thus guaranteeing stability. This gives a handy sufficient (but not a necessary) condition for bounded rates in the system .
In two-dimensional systems describing mixed networks of homogeneous excitatory and inhibitory populations, dynamic stability can be assessed independently of a Lyapunov function. Such a system has dynamic variables λ E (t ) and λ I (t ), respectively, which are coupled by the two equationṡ where α EE , α IE > 0 and α II , α EI < 0. For such a system, the conditions for stability can be stated exactly, here made plausible using a simple heuristic argument. If the system was feed-forward and in equilibrium, the stationary rates of both populations would be λ E = − α EI α EE λ I and λ I = − α IE α II λ E , respectively. Inserting one expression into the other, we obtain This is a self-consistency condition. The system has a transition when or, equivalently, when We tested numerically the validity of this heuristic condition. Figure 4 shows simulation results that confirm the heuristic arguments made here. The condition given here is, in fact, a dissipativity condition which guarantees that the rate decays to 0 in absence of external input. Accordingly, if external input is present, the condition guarantees that the expected rate settles to an equilibrium value. It is interesting to mention that a linear system with the same coupling matrix has an additional stability condition. In fact, the condition given in Eq. 12 is a condition on the determinant of the connection matrix. The system is stable if and only if the determinant is larger than 0. For linear systems, one needs that additional constraints that the trace must be negative, which is not needed here.

SPIKING CENTRAL PATTERN GENERATOR
Central pattern generators have been studied on an abstract level by many authors, for a review see Rabinovich et al. (2006). Lotka-Volterra type equations have also been used as specific models (Venaille et al., 2005). However, no conclusive argument has been FIGURE 4 | Stability properties. Shown is a Monte Carlo experiment that illustrates the condition for dynamic network stability derived in the main text. 200 Networks were initialized, with random couplings drawn from an exponential distribution with mean 1. For each network, a simulation was started and stopped whenever the total population rate would leave the range (10 −10 , 10 200 ). If the total rate was small, the network was classified as stable, otherwise as unstable. For each network 1000 trials were performed. Plotted is the probability of being classified as unstable against the value of . This probability has a sharp transition from 0 to 1 at the value η = 0.

Frontiers in Computational Neuroscience
www.frontiersin.org given why these systems would describe the mean-field dynamics of a neuronal network. In our framework, in contrast, using Eq. 3 and the exponential transfer function, one can map the parameters of the Lotka-Volterra oscillators onto the parameters of a spiking network, which is known to have the same mean-field dynamics.
The main goal of this section is to demonstrate that the spiking networks behave according to the predictions of the rate equations also the latter have stable limit cycles, rather than attractive fixed points.
There are two fundamental alternatives to generate periodic patterns with Lotka-Volterra equations. Already the classical twodimensional Lotka-Volterra system used in population dynamics predicts periodic orbits under certain conditions. More specifically, in that case a whole family of periodic orbits exists, surrounding a stable fixed point. This can be proven by using the fact that the system possesses a constant of motion with closed contour lines. As a consequence, oscillations can occur at different distances from the fixed point, and with possibly different frequencies (Arnol'd, 1984). Thus, due to intrinsic noise of the spiking system, these oscillations will not generate periodic patterns with a stable, precise frequency. Multiplicative point processes indeed exhibit periodic activity patterns with drifting frequencies, when connections are chosen according to the classical Lotka-Volterra system. However, we do not address this type of oscillator in this work.
Alternatively, it has been observed that higher dimensional Lotka-Volterra systems can have stable limit cycles (May and Leonard, 1975). In this case, the noise produced by the system is counterbalanced by the attractivity of the limit cycle. As a consequence, such oscillations have a stable frequency. Furthermore, since the limit cycle is attractive, the spiking network described by the rate equation is attracted to a periodic activity regime.
Mapping the parameters of such Lotka-Volterra systems back to its spiking counterpart allows to construct a spiking system generating precise periodic patterns of spiking. Figure 5 depicts an example of such behavior.
It is actually possible to design spiking networks that implement a given central pattern generator using spikes as signals. We have observed that, although strictly the firing rate description holds only for the mean-field model, periodic behavior is also exhibited in single trials, due to the absence of a competing stable regime. However, the spiking oscillator does not oscillate at exactly the same frequency as the rate equation. This is due to the intrinsically generated noise which continuously perturbs the system, continuously inducing transient behavior. Nevertheless, the oscillation frequencies of the rate equation and of the spiking counterpart are monotonically related, as it is shown in Figure 5, bottom panel.

A LINEAR CLASSIFIER BASED ON WINNER-TAKE-ALL DYNAMICS
Networks with mutual inhibitory couplings that exhibit some kind of winner-take-all dynamics have long attracted the interest of neuroscientists in different contexts (Fukai and Tanaka, 1997;Wang, 2002). Fully connected inhibitory networks of multiplicative processes implement such winner-takes-all behavior, provided the total inhibitory output of each population exceeds its self-inhibition. It can be shown that the probability of winning the competition is largest for the population which receives the strongest total input. According to the rate Eq. 3, the input to the ith population is given by where λ j denotes the firing rate of the jth population in the input layer, comprising K populations in total. Denote by the vector of the synaptic strengths of the K input populations projecting to the ith unit. The total synaptic input to population i can thus be written as Frontiers in Computational Neuroscience www.frontiersin.org Here, v,w denotes the standard scalar product between two vectors v and w. The population i win receiving the strongest input maximizes α i ,λ . Therefore, the highest probability of winning corresponds to the population for which the synaptic input vector is best-matching the input firing rates. In this sense, a completely connected inhibitory network of multiplicative processes performs a linear classification of input vectors.
The situation considered here is similar to the model described by Wang (2002). The network performs a decision (Figure 6), converging with high probability to the configuration where only the population with the best-matching input is firing. Additionally, the outcome of the (random) decision can be manipulated by adding a "bias" to one population (Figure 7). The performance of the classifier degrades if the input rate vectors are almost collinear, and it improves if they are close to orthogonal. This behavior imposes a limit to the number of patterns that can be linearly discriminated, but it also suggests that performance can be boosted by increasing the dimensionality of the input vector. Our results bear a similarity with the forbidden set theory developed by Hahnloser et al. (2003).
The main difference is that we are concerned with the mapping of Lotka-Volterra equations to spiking systems.

DISCUSSION AND OUTLOOK
In our study, we have demonstrated that multiplicative point processes (MPPs) provide a powerful new tool for the study of biological neuronal networks and their functional properties. We have constructed a number of previously unknown spiking network models with functional properties that can be predicted from the dynamics of the associated firing rate equations. In fact, the availability of a higher-level description of the spiking network dynamics in terms of population firing rates opens many new possibilities in terms of mathematical analysis, but also in terms of "engineering" spiking networks with desired properties. In particular, MPPs have the following features: (i) They can be fitted to networks of homogeneous neural populations in terms of their first-order statistics (firing rates). The probabilities for deciding between the different alternatives depend on the network parameters, on the initial states, and on the input levels.
FIGURE 7 | Providing a bias to winner-takes-all networks. Additional spikes are provided at the beginning of the simulation to one of the populations of a symmetric inhibitory network. Shown is the relative frequency that the primed alternative A * is selected. Left panel: relative frequency for selecting A * , displayed as a function of the number of additional spikes provided to one population at the beginning of the simulation. Right panel: reliability of the decision as a function of the background noise for a fixed number (n = 5) of additional spikes. Noise level was changed by increasing the synaptic strength while simultaneously decreasing the input firing rates.
Frontiers in Computational Neuroscience www.frontiersin.org (ii) Their expected behavior (ensemble mean) can be precisely characterized in terms of coupled ordinary differential equations. (iii) The mean-field rate equations give also useful information also about single trial spike statistics.
These features make our framework ideal for use in fundamental brain research, as well as in engineering applications. Specifically, a working theory for interacting neural populations could help in understanding the functional dynamics of interconnected brain areas. Inspired by considerations based on the rate description, research in modeling of the basal ganglia network has already been conducted by the authors (Kumar et al., 2011).

BIOPHYSICAL EXTENSIONS
In this manuscript, we have addressed the construction of neural circuits with prescribed properties, using a first-order non-linear rate description. There are several directions in which the theory could be extended to a more detailed and more precise biophysical theory. In previous work, Toyoizumi et al. (2009) derived rate equations for leaky integrate-and-fire neurons based on a Markov process model combined with a noisy threshold. The Markov approximation is often combined with a Gaussian approximation of the membrane potentials. This method is mainly limited by the complexity of the resulting rate equations. For multiplicatively interacting point processes, there are two possible directions in which the model could be extended. On the one hand, one could strive for more realistic models which include additional dynamic phenomena observed in the activity of neurons. Promising steps have been already performed by including refractoriness into dynamical models (Deger et al., 2010). The differential equations derived in this context could be easily adapted to describe also networks, for the price of increased mathematical complexity. Synaptic transmission delays can be also incorporated into the description, leading to systems of delay Lotka-Volterra equations. Such equations have been extensively studied in the mathematical literature (Bomze, 1983(Bomze, , 1995Tainaka, 1988), also in a stochastic context (Mao et al., 2003;Cattiaux and Méléard, 2010) and with applications to neural networks (Yi and Tan, 2002).
To include more complex neuronal dynamics into the model one can use a relation connecting the derivative of the meanfield and the conditional derivative of a stochastic process (Nelson, 1987, Chap. 9) Here E 0 denotes the best prediction for λ(t ) at time 0, the "mean-field." Based on this relation, seems possible, at least in principle, to derive rate equations for any given neuronal dynamics, provided that the rate only couples to the membrane potential. The precise translation of this formalism into useful rate equations depends, of course, on the details of the neuron model under consideration.

HIGHER-ORDER MOMENTS AND CORRELATIONS
Another interesting question is whether the non-Gaussian approach we have adopted in this work can be adapted to study higher-order statistics of neural networks. Buice et al. (2010) introduced a path integral formulation of neuronal interactions with the aim of deriving expansions for higher-order moments. One limit of their approach is the intrinsic complexity of the path integral formalism.
For multiplicative processes one could attempt a similar treatment exploiting the following idea. Assuming that unit i has rate λ i (0) at time 0, the rate at time t is given by where N j (t ) denotes the number of spikes from unit j until time t. Differentiating, this expression leads tȯ SinceṄ j (t ) is nothing but the spike train produced by unit j, taking expectations and ignoring covariances leads to Eq. 3. In fact, this is an alternative derivation of those equations. The advantage of this derivation is that now it becomes apparent that similar higher-order equations can be derived using the formula After differentiations, mixed third-order moments appear in the derivative of the rate product λ i (t )λ j (t ). Unfortunately, it is not clear what are the terms that can be ignored. Further research will be needed to derive rate equations for higher moments.