Quantization of Map-Based Neuronal Model for Embedded Simulations of Neurobiological Networks in Real-Time

: The discreet-time (map-based) approach to modeling nonlinear dynamics of spiking and spiking-bursting activity of neurons has demonstrated its very high efficiency in simulations of neuro-biologically realistic behavior both in large-scale network models for brain activity studies and in real-time operation of Central Pattern Generator network models for biomimetic robotics. This paper studies the next step in improving the model computational efficiency that includes quantization of model variables and makes the network models suitable for embedded solutions. We modify a map-based neuron model to enable simulations using only integer arithmetic and demonstrate a significant reduction of computation time in an embedded system using readily available, inexpensive ARM Cortex L4 microprocessors.


Introduction
The functionality of a brain and other neuronal systems of living organisms rely on the synaptic coupling organization among populations of neurons forming a function specific network of neurons with proper pathways of electrical activity in forms of short electrical spikes known as action potentials. The spikes generated in a neuron propagate along its axon towards synaptic inputs of the other neurons making the network achieve its functional goals.
Dynamical properties of spiking electrical activity of neurons, conditioned by the intrinsic properties of neurons or evoked in the response to external stimuli or other synaptic inputs, provide critical elements in shaping proper neuronal activity and, therefore, the functionality of neurobiological networks. This spiking activity propagating along the synaptic pathways is a key element of information processing in the brain, classification of sensory inputs, decisionmaking and motor control (Shepherd, 2004). Studies of the mechanisms behind the formation of network spatio-temporal patterns of spiking activity rely on neuronal models that can capture essential characteristics of spiking activity for specific types of neurons (Dayan and Abbott, 2001).
Numerous recent studies focus on the development of artificial (engineering) systems that use elements of spiking activity to create a network that could perform some elements of brain functionality, sensory processing, motor control (Dura-Bernal et al., 2014) and other functions. They mainly rely on the development and implementation of a spiking neuron model and synapses that replicate the dynamics of neurons as close to real biological neurons as possible that fits in computational hardware. Two main directions of these studies are the use of analog computing simulating the network with specially designed VLSI integrated circuits (Indiveri et al., 2006;Silver et al., 2007) and numerical simulation of the model with specially designed microchips, GPUs FPGAs or DSPs (Yavuz et al., 2016;Zbrzeski et al., 2016;Cheung et al., 2016;Brette et al., 2007 and references therein).
This paper focuses on the design, implementation and validation of a neuronal model that could be used in real-time simulations of relatively large networks using off-the-shelf 32-bit microprocessors. We base our design on the map-based approach to modeling spiking neuronal activity in discrete times. Such models are very computationally efficient, realistically replicate the dynamics of spiking activity and have been successfully used in many studies of neuronal activity in large-scale cortical networks, olfactory systems and Central Pattern Generator (CPG) networks. The mapbased model of swimming motor control of the CPG network of a Lamprey was implemented using a TI DSP and operated in real-time on a Lamprey-based robot (Westphal et al., 2011).

Map-Based Modeling of Neuronal Activity
Consider a model of spiking and spiking-bursting activity designed in the form of a two-dimensional map and generating neuronal activity states at discrete moments of time. This model is very efficient in numerical simulation since its time step of consecutive iterations can be made to be on the order of spike duration and the dynamics of the spiking does not depend on the step size, which is usually fixed at 0.5 or 1.0 ms, depending on the network model. The equation of such a map can be written in the form proposed in (Rulkov, 2002): (1) where, variables x n and y n are fast and slow variables describing the state of the neuron at discrete moments of time indicated by n. β n and σ n are input variables that model synaptic inputs and injected currents and other stimuli. The intrinsic dynamics of the map, Equation 1 is controlled by parameters µ, σ and α. The spiking in the model is due to the specific shape of the nonlinear function which can be written in the following form: Note, Equation 2 slightly differs from the function proposed originally in (Rulkov, 2002). As shown in (Komarov et al., 2016) it is more efficient for large-scale network simulations, simplifies the equation and improves the spike shape consistency.

Nonlinear Map of Spiking
Spike generation in model Equation 1 is provided by the nonlinear dynamics of fast variable x n and can be considered independently from the rest of the model using the following one-dimensional map: where, u is treated as a control parameter. Figure 1 illustrates the mechanism of spike formation and positions of stable and unstable fixed points x sp and x up , respectively, indicated by green and blue filled dots. By definition, fixed points, in a one-dimensional map, are located at the intersection of the first segment of the function Equation 2 with the diagonal (x n+1 = x n ), Fig. 1. Therefore, when the parameter u increases the fixed points move closer together, merge and then disappear when u crosses bifurcation value . This value is found using the fixed point equation x * = α(1-x * ) −1 + u, when it has a single solution (merged solution) in domain x * < -0.5. Such bifurcation scenario of the fast subsystem fixed points is a key element of many simple models of spiking neurons, such as nonlinear integrate-and-fire and the Izhikevich model. It is important for our study to emphasize that near the bifurcation state the dynamics of the system is very sensitive to the resolution set for variable x n because changes of the variable in that region of phase space can be arbitrarily small. Simulations with insufficient resolution can treat such small changes as zero and artificially stall the trajectory.
The one-dimensional map, Equation 3, demonstrates only the solutions in the form of steady state or tonic spiking and are typically used within the range of parameters 3 < α < 6 and u cr -0.3 < u < u cr +0.2. This is sufficient to mimic the dynamics of active sodium (Na+) and potassium (K+) pumps in generating action potentials (spikes). To expand the dynamics of this model to describe behavior of various types of neurons the parameter u is substituted with a variable or a set of variables that capture dynamics of relatively slow ionic channels of the neuronal membrane. An example of such a model is Equation1, which has been used in many studies as a reduced model of a Regular Spiking Cell in a cortical network  and a Spiking-Bursting neuron in modeling CPG networks (Ayers and Rulkov, 2008).
The detailed parametric analysis of nonlinear dynamics of Equation 1 with β n = 0 and σ n = 0 can be found elsewhere (Shilnikov and Rulkov, 2003). An example of the spiking-bursting oscillations that occurs in the model without any external inputs is shown in Fig. 2. The phase portrait in Fig. 2a shows the trajectory (black dots connected with a line) and manifolds of slow motions (blue curves). The corresponding waveforms of fast x n and slow y n variables are shown in Fig. 2b and 2c. The response dynamics of the model is controlled with two input variables β n and σ n . The use of two input variables significantly expands the possibilities in shaping the desired response behavior of the models without adding new dynamical variables .

Modeling of Response Spiking Behavior
In an experimental study, the type of the measured neuron is usually identified by analyzing its response characteristics to an electrical stimulus, such as a rectangular pulse of current, injected into the neuron through a microelectrode. In the model Equation 1 any such received stimulus is described by input variables β n and σ n , which can written as a sum of all inputs: are weights used to tune the model to produce a desired response behavior to fit a particular type of neuron.
Examples of tuning this model to regular spiking cortical cells can be found in (Bazhenov et al., 2005;Rulkov and Bazhenov, 2008). Effects of spike deceleration and after-hyperpolarization rebound bursts controlled by parameters β D and σ D are shown in Fig. 3, where α = 3.9, µ = 0.001, σ = 0.0, β D = 0.2, σ D = 1.0 and duration of the D n I pulse of 400 iterations that starts from initial state The mechanism of shaping response behavior can be examined by following the map trajectory in the phase portraits, Fig. 3, where blue curves x sp (y) and x up (y) correspond to the stable and unstable manifolds of slow motions, respectively. The green fixed point indicates the baseline state of the neuron while the red fixed point corresponds to the perturbed state during pulse action. The red dotted arrow indicates fast motion at pulse onset and solid red-a slow drift during pulse action. Similar behavior at the termination of the pulse is indicated with green arrows.

Dynamics of Synaptic Inputs
A spike generated by a neuron propagates towards the other neurons that receive and respond to it due to the synaptic connections located in the dendrite of the receiving (postsynaptic) neuron. There are electrical synapses, working as a conductor between neighboring neurons and chemical synapses where a spike propagating over an axon can reach dendrites of rather remote neurons forming a complex network structure. The dynamics and network organization of chemical synapses play a key role in the formation of neurobiological network functions. The basic elements of synapse dynamics is the modulation of membrane conductance at synapse contact evoked by a received presynaptic spike, which is characterized by elevation and relaxation times (Destexhe et al., 1994). There are known several distinct types of synapses characterized by different timing rates and conductance. We adopt a model that captures the response to all presynaptic spikes of the same type of synaptic connections in the postsynaptic neuron using the following simple set of equations. The synaptic current syn n I computed using synapse conductance g n and variable x n of the postsynaptic neuron: where, parameter x RP defines reverse potential and, therefore, the synapse type (excitatory or inhibitory). The temporal modulation of conductance is modeled as the combination of two exponential curves:

Quantization of Model for Integer Computation
The main goal of this study is to modify the models of neurons and synapses, considered above, for computation using readily available inexpensive off-the-shelf hardware utilizing only integer operations and verify model performance. To allow such computation we need to rescale all variables and parameters to provide sufficient resolution that supports proper dynamical behavior of the model. The method of such scaling in the case of the Izhikevich model (Izhikevich, 2003) has been discussed in detail in (Jin et al., 2008). In our case the map-based model can be rewritten as: The values of scaling parameters p x and p y are selected to be as large as possible assuming that all computation still fits within 32-bit operations of the microprocessor. Note that in Equation 4 we shifted slow variable, y n by the value

Experimental Studies of Model Performance
To study the efficiency of models and demonstrate their computation in real time we built a simple computer using an STM32F427 microprocessor. We utilized an old empty Printed Circuit Board (PCB), which we had designed before for another project and rewired it to fit our study of neuron simulations. The populated STM computer PCB with an array of 4 12-bit DACs, serial DATA link and power control we built for this study is shown in Fig. 4.
The STM32F427 was set to run with a 180MHz clock. It also supports 8 Analog Input Channels with internal 12-bit ADCs. The microprocessor workflow is designed to operate in one of two modes selected by a command switch in the code. In the first mode, we call it real-time continuous operation, the microprocessor runs all computations for each iteration within controlled 500µs time intervals taking input signals from Analog Inputs and showing the neuronal responses on the Analog Outputs. In this case only 4 selected variables are available at the Analog Outputs. The other "two state" outputs, such as spikes, are made available at digital output pins that can be used to control external actuators, generate motor activity controlled by a simulated CPG model or other tasks. The use of serial DATA in this case is very limited. Figure 5 and 6 present such a study for three implementations of the same neuronal model using double Equation 1, int32_t Equation 5 and float Equation 1 types for variables and parameters.

Simulations of Map-Based Networks
To deal with simulations of large networks with various topologies of synaptic coupling we developed code that can be compiled for simulations on a PC or for real-time operation in the microprocessor using only a minor change in the parameter definitions. In the PC case, the network topology and all parameters are defined in a text file, which is read by the simulation program. In the case of embedded implementations, an additional program uses the same text file to generate a header file that is included in the compilation of the microprocessor based code. In both cases the network description is translated into data structures in the PC or microprocessor memory. The network simulation code represents models for neurons, synapses and interconnection topology and delays as structures. This provides convenient flexibility with minor computational overhead.
The synapse model includes code for spike propagation delays between pre and post-synaptic neurons. Each spike delay line is implemented with an int32_t variable and propagation is accomplished by a binary shift applied to the variable every iteration of the network. Time delay is achieved by reading the proper bit in the presynaptic delay variable as it is defined in the synapse connectivity of the postsynaptic neuron. A block diagram illustrating the structure of synapse simulation code in the case of two types of synapses, Excitatory (E) and Inhibitory (I), is presented in Fig. 8. To simplify the network diagram description the argument n of variable x i (n) represents the iteration step (time) and subscript i, the neuron index. The strengths of synaptic connections between neurons are defined using weights _ synType post pre W applied to received spikes in the summation block (in the text below we call it Presynaptic Inputs). Note that weights can be made to vary in time during simulations to model various mechanisms of synaptic plasticity.

Map-Based Network Model
As an example, that is useful for validating the dynamics and testing the performance of integer computation in an embedded microprocessor environment, we consider a simple network consisting of two neurons demonstrating the effects of spike adaptation and rebound bursts coupled with chemical inhibitory synapses as shown in Fig. 9.
The baseline state of each neuron is set to be in the silent condition. The injected rectangular current excites the first neuron, which starts spiking and inhibits the second neuron. The second neuron demonstrates the effect of a rebound burst after the inhibition from the first neuron is terminated. Then a pulse of injected current is applied to the second neuron and the process repeats in the opposite sequence.
An example of such simulations computed using floating point operations is presented in Fig. 10a Fig. 10c and 10d. It uses scaling multipliers p x = p y = 16384 and p s = 1000. Panels Fig. 10a and 10c show response waveforms of the neuronal activity and panels b and d waveforms of the synaptic conductance activated by the presynaptic spikes. One can see that both simulations exhibit very similar dynamical behavior including very fine elements of neuronal excitability. Note that in both cases the first rebound burst spikes evoke rebound spikes in the opposite neuron. This indicates that the subthreshold dynamics of both neuronal models are the same, despite the difference in the computation method. At the same time in the case of the scaled model code for integer operations the computation time is about ten times more efficient than in the case of floating point including all the overhead components associated with the structure architecture of the code.
Measurements of computation time in the microprocessor for various sizes of network and connection topology performed for the cases of float and int32_t were used to derive approximate estimations of time needed for computation of each element of the simulated network when it contains 5 or more neurons. The results of this study are summarized in Table 2. For example if the network has 50 neurons and each neuron has two types of synaptic inputs (e.g., excitatory and inhibitory) with 20 presynaptic inputs of each type, each time iteration of the network in the microprocessor will require 50*(1µs+2*0.7µs+2*20*0.076µs) = 272µs for int32_t and 50*(10µs+2*4.4µs+2*20*0.095µs) =1130µs for float. The results of computational time analysis demonstrate that quantization of a map-based model for dynamics of neurons and synaptic conductance give a significant improvement in simulation speed. The use of a linear model of synaptic conductance, Equation 4, enables one to sum presynaptic inputs with corresponding weights very efficiently in both methods of simulation, because the floating-point Multiply-And-Accumulate (MAC) loop is optimized in STM32F4xx microprocessors.

Conclusion
The presented results demonstrate that a transition from floating point to integer implementation of the map-based models can significantly speed up computation of large networks of neurons. In the reported study we used neuronal and synaptic models scaled and coded using 32-bit integer variables and parameters. This increased the computational efficiency of the model Equation 1 and 2 about ten times and the overall network simulation, about five times. This improvement is very important and beneficial for real time implementation of neurobiological network models using low power off-the-shelf microprocessors. All computations presented in this study were done using an STM32 F427 running at a clock rate of 180MHz and consuming about 60mA at 4.0V. The size of the simulated network did not have a significant effect on the power consumption. This power consumption figure can be reduced by either reducing the processor clock rate and/or placing the processor in low power mode after it completes an iteration and subsequently resuming full power at the start of the next step interval.
The presented results of quantization are focused on the use of 32-bit integers, where the map-model shows a very good behavioral match with the original model Equation 1 and 2. Further decrease of resolution may be possible and has been reported in similar types of models, see for example (Jin et al., 2008). However, in our studies we observed that use of quantized models with resolution of 16-bits and below sometimes cause the onset of artificial instability and unexpected oscillations in sub threshold dynamics. The dynamical mechanisms and conditions of such distortions needs to be studied in more detail as the use of lower resolution could be beneficial for FPGA and other similar implementations.
The considered neuronal map-based models and the developed software and hardware implementation reveal a power efficient, readily available and inexpensive real-time engine for medium size neurobiological network simulations. These network models are good candidates for use in the design of wearable electronics for controllers in bioinspired prosthetic, sensors and robotic devices.