changes in the development of technologies that we use every day. Sometimes the novelty and the complexity of those frameworks make that takes some time to be fully exploited in other related fields of science such as Computational Neuroscience, or different areas of Engineering.
That enforces the need to develop good tutorials or premiers to explain how to implement the algorithms that we need to ask the scientific questions of interest, and also to be able to solve different problems using the newly available tools.
One relevant problem is how to build the models for the study of dynamical systems and how to extract the relevant information to be able to answer the scientific questions of interest.
Neural Networks are algorithms that allow us to model different systems. According to Universal Approximation Theorem, a neural network with one hidden layer containing a sufficient but finite number of neurons can approximate any continuous function to a reasonable accuracy under certain conditions for activation functions [5] . In particular, RNNs allows to model dynamical systems, since a dynamical system can be approximated by continuous-time RNN [6] .
The problem of training neural networks to perform different tasks is relevant across various disciplines that go beyond Machine Learning. In particular, RNN are of great interest in different scientific communities. They are widely used in Neurosciences to describe the behavior of the cortex, an area that presents great recurrence in its connections [7] . These systems also have great relevance concerning control systems and other areas such as electronics [8] [9] [10] .
In the present work, the modeling of a complex system is presented, using Tensorflow and Keras. The reason for this selection is that these new scientific libraries are open source, and its use is rapidly growing. Good documentation can be found for software development [11] [12] [13] , and also new tools such as Google Colaboratory allow implementing and testing directly online network implementation, with the focus on ML.
In this work, a RNN was chosen because it is relevant in various fields. It was trained to perform a time-series processing task, that was inspired in Neuroscience tasks [14] . The implementation of the network, the training, and the tools are presented, as well as different forms to obtain the information that allows describing the system appropriately.
Training a RNN to perform temporal tasks has many difficulties and can be done through various paradigms. Here it is proposed to approach the problem through supervised learning. The entire procedure is described in detail, using open-source software. Gradient descendant minimization was used, to take advantage of different optimized implementations of the current algorithms available.
The Flip Flop task was chosen as a case-example, in fact, a 3-bit memory, which is a task composed of Flip Flops. The parametrization of the task is as described in [15] , but also revisited here. All steps will be described in detail, from the parameterization of the task to the description of the dynamics of the trained networks. This example is used to show how the problem of training networks can be studied using these new computing tools applied in any temporal task in general and also discussing the limitations that networks have and the alternatives to solve them.
The rest of the paper is organized as follows. First, in Section 2, the description of the dynamics of the RNN and the discretization is presented. In Section 3, the task parametrization is shown. Then, in Section 4, the training protocol is shown. In Section 5, the results and different analyses of the network and tools are discussed. Finally, in Section 6 final remarks are presented.
The dynamics of the units in the RNN model is described by equation 1 [16] , where units have index i, with i = 1, 2..., n.
In Equation 1, the matrix elements w Rec ij are the synaptic connection strengths of the matrix W Rec and w in ij the matrix elements of W in from the input units. x j are the components of the vector X of the input signal. τ represents the time constant of the system and σ is a non-linear activation function.
The network is fully connected, and matrices have weights given by a certain parametrization of interest. In our case, a normal distribution with zero mean and variance 1 N . The network has three layers: the input, the recurrent hidden layer, and the output layer. The readout output in terms of the matrix elements w out ij from W out is described by Equation 2 .
For this work, it was considered σ() = tanh() and τ = 1, without loss of generality. The model is discretized using the Euler method. A simple schema of the model is presented in Figure 1 , in this case, with three inputs and three outputs corresponding to the inputs and memory states of the 3-bit Flip Flop.
As described in the Section 1, model is implemented in Python using Keras and Tensorflow [1, 2] , which allows making use of all current algorithms and optimizations developed and maintained by a massive research community.
In vector form, the equations 1 and 2 can be written as: 
and respectively:
The system represented by Equation 1 is approximated using Euler's method with a step time δt. A value of τ = 1 was considered. Then the dynamics of the discrete-time RNN is done by means of:
The value of the time step is δt = 1 mS,for the time evolution. Then, from equation 5, the activity of the recurrent units at the next time step is given by Equation 6 .
Tensorflow has a recurrent layer implemented to represent Equation 6, where it is possible to chose the initialization of parameters, units and activation function.
tf.keras.layers.SimpleRNN( units, activation="tanh", kernel_initializer="glorot_uniform", recurrent_initializer="orthogonal",**kwargs)
Code for a Recurrent layer defined in Tensor Flow.
To define a Recurrent Network with the topology of Figure 1 , it is possible to build a sequential model with layers such as:
Code for the sequential model defined for Figure 1 .
Where input shape = (N one, 3) is the shape of the vector input, activation = tanh corresponds to the definition of activation function, and Dense is the output layer. In this way, we have completed the first step that consists of the model definition in terms of the code.
Once the task was chosen, the requirements must be translated into the algorithm that allows generating the training set. To parameterize the task, the following criteria is considered for the present work:
-The possible states of the Flip Flop are represented in such a way that a positive pulse represents a set and a negative pulse represents a reset. -The state of the output will change, corresponding to the input command.
-A certain delay was considered after the falling edge of the input signal.
The training data set consists of time-series with pulses of fixed duration that represent SET and RESET. Those signals can be activated randomly and are separated by a random time interval. In all time-series, a certain noise level has been superimposed on the input. Each input elicitate a target output according to the Flip Flop rule: if we have a set signal or positive pulse, the output is high-state. If we have a reset signal or negative pulse output is lowstate, otherwise, the output remains in the previous state.
The number of inputs in the network corresponds to the number of bits that can be stored. A Flip Flop is a one-bit of memory. In this case, we have registers formed by three Flip Flops (a 3 bits memory).
To generate the full training data set, it is necessary to generate sets of tensors of size sample size with the input time series of length time s eries l engh, for each of the three inputs and outputs. To do that efficiently, we can use Numpy arrays [17] . With present work, sample code is provided in present work to generate the Flip Flop data set. Three of components of the full set x train − y train are shown in Figure 2 . The target output, y train, was simulated with a time delay answer of 20 mS. Each row (and color) corresponds to one of the inputs and each column to a different sample. In this way, each training sample consists of a Numpy array. Training data set pairs defined as Numpy arrays. 
The training methods for neural networks can be unsupervised or supervised. We will focus on applying a supervised method.
Studies, where one kind of gradient descent method is applied, stands out in the literature. An example is the paradigm of Reservoir Computing paradigm with liquid-or echo-state networks [18] , where the modifications of the network weights are made in the weights of the output layer, W out .
Other outstanding methods, were developed by Sussillo and Abbot. They have developed a method called FORCE that allows them to reproduce complex output patterns, including human motion-captured data [19] . Modifications to the algorithm have also been applied successfully in various applications [20] .
The other method used for the estimation of the gradient in RNNs is called Backpropagation Through Time (BPTT), and then some optimization method for minimizing the gradient. Given the recent advances in the implementation of this method with the Open-source libraries previously mentioned, this is the method explored here. There are also other back propagation-based methods been published recently, for example in [21] , authors propose to use fractional calculus to improve the conventional BPTT.
In this work, supervised learning was used, with standard backpropagation through time implementing an Adaptive SGD training method provided by Keras framework.
First, recurrent weights were obtained using a random normal distribution with the orthogonal condition on the matrix. Noisy square pulse signals, as presented in Figure, were used at the inputs.
The appropriate loss function to train the model is the mean square error between the target function and the output of the network. It is defined as:
where Z target j (t) is the desired target function and Z j (t) is the actual output. The training sample consisted of more than 15000 samples.
In terms of the code, the previously defined steps are: compile the model and train it. In the first step, the loss function and the optimizer algorithm are selected. In the training stage, different detail on the training data set, epochs, and other characteristics of training. The parameters of a neural network are the weights of the connections. In this case, these parameters are learned during the training stage. The hyperparameters are parameters of your neural network that can not be learned via gradient descent or some other training method. These include learning rate, number of layers, number of neurons in a given layer.
Tuning the hyperparameters means the process of choosing the best values of them. Typically this is done by evaluating the performance of the network on a validation set. Then, we have to change the hyperparameters and reevaluate, choosing the values that give the best performance on the validation set.
How do we choose these values? Often there is good standard initialization related to each particular task of interest, here an example of the criteria is provided for the Flip Flop case in Table 4 . In our example the parameters are also given in the table.
Another aspect to consider is the regularization of the model. Regularization refers to train our model well enough that it can generalize over data it hasn't seen before.
To summarize, in the training stage is where we have to consider the size of the network, data set, noise, regularization terms that are appropriate for the considered task, in our case, the Flip Flop. Table 1 Model's parameters and criteria for the network's implementation and training.
Finally, after training, we can obtain a set of RNN that can perform the tasks. Now we search the aspects to analyze regarding the network's collective behavior. Also, we chose a way for the visualization of the model obtained and tools that allow us to extract the relevant information. For example, it is possible to visualize the connectivity matrix (recurrent weight matrix), as it is shown on the left side of Figure 3 . The columns represent the output connection of the i-neuron, and the rows are the input connection. The color bar represents the intensity of the connection considered. With such visualization, it is not clear which is the relevant information, apart from the fact that after training, most of the weights remain close to zero. In case that we have imposed constraints on the connections, such as Dales' Law, or any particular constraint they will be visible in this stage and this representation will be more useful.
Different Linear Algebra operations are available in the Numpy that are optimized to be used with the array structures. If we perform, for example, a decomposition of W rec in their eigenvectors and eigenvalues, we can obtain the eigenvalue distribution as it is shown in the right side of Figure 3 .
from numpy import linalg as LA eigenvalues, eigenvectors= LA.eig(Matrix)
How to use linear algebra library from Numpy for eigenvalue decomposition.
We can visualize that, except for a small group of eigenvalues that migrated out of the unit circle, the rest remain on the unitary circle (orthogonal initial condition). This situation is repeated in all the obtained simulations. Eigenvalues outside the unitary circle seem to be related to the behavior (or modes) observed for the different stimuli. Other additional information, related to the connectivity matrix, could also be obtained. Other possible studies are related to the response of the network when applying the different stimuli.
Since we have a large number of units, and for each and activity vector, dimensionally reduction methods are appropriate to analyze such behavior, and it has been used widely in different works related to large-scale neural recordings [22, 23] .
Scikit-learn [24] is, for example, a python open source library based on Numpy that allows us to perform dimensionality reduction, feature extraction, and normalization among other efficient methods for predictive data analysis. Such decomposition could be, for example, Principal Component Analysis (PCA) or Single Value Decomposition (SVD). from sklearn.decomposition import PCA from sklearn.decomposition import TruncatedSVD How to import scikit learn libraries to perform single value decomposition and principal component analysis.
We can use such tools to extract relevant features of our system. It is well known that the different memory states in a 3-bit memory distribute in the vertex of a cube-like form in the space state [15] . It was shown by then when they explored the hypothesis that fixed points, both stable and unstable, and the linearized dynamics around them, can reveal aspects of how RNNs implement their computations.
Here, a data set was built for testing. It generates eight different memory states, as it is shown on the left side of Figure 4 . Time series of 600 mS were considered to generate all the possible different memory states of the 3-bit memory by choosing the correct commutation for the inputs. Input and output response is shown.
The testing set is injected into the network (right upper panel of the Figure) and then we can analyze the activity of the units using SVD applied on the activity vector H(t). The behavior of the system was represented into the three axes of the greatest variance. The bottom right part of Figure 4 shows the activity in the reduced state space (3-dimensional) . Each vertex corresponding to each memory state is marked in different colors.
As it is well known, different variations of the realizations, in terms of weight distribution and dynamical behavior are possible when training networks for the same task. Such is the case is shown in Figure 5 , where four different realizations for the trained networks elicited with the same testing data set are shown. Vertex is distributed in different positions, and a cube-like structure is rotated in different spatial directions. We can study and classify the behavior of the obtained systems by comparing different network realizations.
Additional analysis could be considered depending on the interest in the aspect to be studied. Here a minimal analysis was proposed describing in detail the steps, visualization tool, criteria, and implementation. The code is provided for training and analysis but also to be a framework open for different tasks or additional studies.
In this work, all steps to build and analyze an RNN have been presented. We started from the model description in terms of the equations, discretization, and code implementation. Then we described the task parametrization and network training protocol. We also presented a set of tools to analyze the results using open-source scientific libraries making use of the different visualization tools that allow extracting relevant features.
We use the Flip Flop task as an example, but other relevant tasks could be considered. The use of the new open-source scientific tools that are designed and maintained for large communities, such as the tools used here, allows enhancing research. This is why we are currently using tools that are more transparent in terms of code and documentation because they are open to being modified and improved for thousands of users.
The author declare that she has no conflict of interest. 
Code is provided in an open repository after paper publication and with a MIT License. https : //github.com/katejarne/3 − bit − F F − tutorial