Synthesis and Analysis of the Fixed-Point Hodgkin–Huxley Neuron Model

: In many tasks related to realistic neurons and neural network simulation, the performance of desktop computers is nowhere near enough. To overcome this obstacle, researchers are developing FPGA-based simulators that naturally use ﬁxed-point arithmetic. In these implementations, little attention is usually paid to the choice of numerical method for the discretization of the continuous neuron model. In our study, the implementation accuracy of a neuron described by simpliﬁed Hodgkin–Huxley equations in ﬁxed-point arithmetic is under investigation. The principle of constructing a ﬁxed-point neuron model with various numerical methods is described. Interspike diagrams and refractory period analysis are used for the experimental study of the synthesized discrete maps of the simpliﬁed Hodgkin–Huxley neuron model. We show that the explicit midpoint method is much better suited to simulate the neuron dynamics on an FPGA than the explicit Euler method which is in common use.


Introduction
The implementation of artificial neural networks (ANNs) in modern electronic devices requires different network topologies depending on the solving tasks, which vary from clustering and classification to pattern recognition. Based on the available training data, the learning model for the network can be supervised or unsupervised. The first principle, relying on labeled datasets, takes place in Convolutional Neural Networks (CNNs), which have recently shown impressive performance in cognitive tasks such as recognition [1,2] and prediction [3,4]. During supervised learning, the error between the input target and the output of a CNN is minimized by adjusting the synaptic weights of the network. High accuracy can only be achieved with a large number of training examples and algorithm cycles, which entails the application of massive-power-consuming computers. This imposes a significant restriction on the application of CNNs in embedded systems, motivating the development of neural networks of the third generation-Spiking Neural Networks (SNNs). SNNs, being unsupervised learners, provide the most realistic natural neural network emulation. Like neurons in living organisms, SNNs encode information as sequences of spikes, the precise timing between which is used to update the synaptic weights. The main advantages of SNNs over the previous-generation ANNs are their computational speed, superior classification abilities, and efficiency in control problems, The circuit dynamics is described by the following differential equations: where V is the potential difference between the neuron's membrane and the environment, R is the recovery variable, I is the input current, C = 0.8 µF/сm2 is the membrane capacity, and τ = 1.9 ms is the recovery time constant. The right-hand side of the first equation is the sum of the input current and Na+ and K+ ion currents. The passive leakage current of the original HH model is absorbed into the Na+ current. The second equation represents the behavior of the recovery variable R, which describes the K+ channel as a memristive element.
To move from a continuous model of dynamical system (1) to the set of investigated ordinary differential equation (ODE) solvers, the following methods of numerical integration were used: the Explicit Euler method (EE), the Semi-Explicit Euler method (SEE), the Explicit Midpoint method (EMP), and the Modified Explicit Midpoint method with a smoothing step (MEMP). The choice of explicit methods of the first and second order was determined by the simplicity of their implementation in integer representation and the visibility of the observed numerical effects.

Floating to Fixed Point Model Conversion
Conversions of ODE solvers with floating points to integer solvers were implemented using the approach described in [20]. First, the minimum and maximum possible values of each state variable of system (1) were determined by preliminary simulation. After that, the largest modulus value was selected among all state variables and system parameters to determine the required number of bits to store the integer part of the fixed-point data type (FXP). It was shown that to store the integer part of the state variables and parameters required not less than six bits. Thus, all state variables and constant coefficients of system (1) were converted to the FXP data type, where one bit is allocated for a sign, seven bits are allocated for storing the integer part, and the remaining bits are allocated for the fractional part. It should be noticed that the number of bits of the integer part was increased by one in order to guarantee an overflow avoidance of the bit grid during calculations. In the research, integer models of 32-bit and 64-bit FXP data types were explored, the fractional parts of which take 24 and 56 bits, respectively. Investigation of models with a longer bit grid is not of interest at the moment, since modern computing platforms do not support arithmetic operations on a hardware level for machine words longer than 64 bits. For the same type, the software implementation of these operations requires additional hardware costs that negate all the benefits of using an integer data type. It would be advantageous to get 16-bit FXP system models, which would allow the use of low-power and less-expensive hardware platforms for implementing large neural networks. The circuit dynamics is described by the following differential equations: where V is the potential difference between the neuron's membrane and the environment, R is the recovery variable, I is the input current, C = 0.8 µF/сm2 is the membrane capacity, and τ = 1.9 ms is the recovery time constant. The right-hand side of the first equation is the sum of the input current and Na+ and K+ ion currents. The passive leakage current of the original HH model is absorbed into the Na+ current. The second equation represents the behavior of the recovery variable R, which describes the K+ channel as a memristive element.
To move from a continuous model of dynamical system (1) to the set of investigated ordinary differential equation (ODE) solvers, the following methods of numerical integration were used: the Explicit Euler method (EE), the Semi-Explicit Euler method (SEE), the Explicit Midpoint method (EMP), and the Modified Explicit Midpoint method with a smoothing step (MEMP). The choice of explicit methods of the first and second order was determined by the simplicity of their implementation in integer representation and the visibility of the observed numerical effects.

Floating to Fixed Point Model Conversion
Conversions of ODE solvers with floating points to integer solvers were implemented using the approach described in [20]. First, the minimum and maximum possible values of each state variable of system (1) were determined by preliminary simulation. After that, the largest modulus value was selected among all state variables and system parameters to determine the required number of bits to store the integer part of the fixed-point data type (FXP). It was shown that to store the integer part of the state variables and parameters required not less than six bits. Thus, all state variables and constant coefficients of system (1) were converted to the FXP data type, where one bit is allocated for a sign, seven bits are allocated for storing the integer part, and the remaining bits are allocated for the fractional part. It should be noticed that the number of bits of the integer part was increased by one in order to guarantee an overflow avoidance of the bit grid during calculations. In the research, integer models of 32-bit and 64-bit FXP data types were explored, the fractional parts of which take 24 and 56 bits, respectively. Investigation of models with a longer bit grid is not of interest at the moment, since modern computing platforms do not support arithmetic operations on a hardware level for machine words longer than 64 bits. For the same type, the software implementation of these Electronics 2020, 9, 434 4 of 16 operations requires additional hardware costs that negate all the benefits of using an integer data type. It would be advantageous to get 16-bit FXP system models, which would allow the use of low-power and less-expensive hardware platforms for implementing large neural networks. However, preliminary simulation showed that system (1) cannot be adequately represented in a case of such limited bit length.
After converting all variables and constants to the FXP data type, adequate simulation is still not guaranteed. In the process of calculations, the values of state variables can overflow the bit grid. This situation is most probable while calculating the (17.81 + 47.71V + 32.63V 2 ) part of the first equation of system (1). This was considered while generating the FXP solvers and compensated by organizing the correct order of calculation. For example, the algorithm suitable for the integer implementation of the ODE solver of system (1), constructed on the basis of the Euler method, looks as follows: where h is the constant integration step size, and i is the solution time.

Accuracy of the Fixed Point Simulation
Let us estimate the accuracy of various implementations of system (1) on the basis of the difference between the values of the state variables of the 32-bit and 64-bit FXP models in comparison to the floating-point model (double, DBL) in the time domain. Figure 2 provides charts of absolute error for the EMP-based solvers.
part of the first equation of system (1). This was considered while generating the FXP solvers and compensated by organizing the correct order of calculation. For example, the algorithm suitable for the integer implementation of the ODE solver of system (1), constructed on the basis of the Euler method, looks as follows: where h is the constant integration step size, and i is the solution time.

Accuracy of the Fixed Point Simulation
Let us estimate the accuracy of various implementations of system (1) on the basis of the difference between the values of the state variables of the 32-bit and 64-bit FXP models in comparison to the floating-point model (double, DBL) in the time domain. Figure 2 provides charts of absolute error for the EMP-based solvers.   Figure 3 shows that error accumulation leads to a complete divergence between the trajectories of the two models. For the 32-bit EMP-based solver, it occurs after 30 ms of simulation, while the 64-bit solver switches to a different operational mode after 400 ms. It should be noted that the process of distancing the trajectory of the integer solution from the trajectory of the solution of the original algorithm is inevitable not only due to the limitations of the bit grid but also because of the different order of arithmetic operations in the solvers. However, the stability of the solution is retained, as was confirmed in a series of experiments.

Resonant Spike Generation
The simplified HH neuron model can demonstrate a stationary mode, damped subthreshold oscillations, and spike generation mode. Therefore, considering Izhikevich's classification [21], the HH model is a bistable resonator.
In Figure 4 the represented dynamical maps demonstrate the behavior of system (1), which was implemented using the solvers under investigation, when applying a pulse signal to its input with amplitude I = 4.5 µA and different values of period T and pulse width w. The period T ranged from 20 to 40 ms with step ∆T = 0.05, and the pulse width ranged from 10 to 30 ms with step ∆w = 0.05. The white areas on the diagrams of Figure 4 correspond to the stationary mode and the mode of subthreshold oscillations, while blue areas correspond to the spiking mode for the considered 32-bit and 64-bit FXP and DBL solvers built using the following methods: EE, SEE, EMP, and MEMP. The integration step size for all experiments was the same h = 0.01 ms.
The results showed that the excitability regions mainly form diagonal bands. The number of bands depends on the integration method used. The solvers based on the IE and SEE methods create the smallest excitability area: just three bands. The EE method leads to the highest number of bands in the excitability area, while the second-order methods EMP and MEMP give intermediate results.
The transition to an integer data representation reduces the excitability area for each of the considered solvers. Table 1 provides a numerical assessment of this phenomenon. The number of lost excitability areas for the 32-bit and 64-bit FXP solvers relative to the DBL solvers for given simulation parameters was calculated as a percentage.   Figure 3 shows that error accumulation leads to a complete divergence between the trajectories of the two models. For the 32-bit EMP-based solver, it occurs after 30 ms of simulation, while the 64-bit solver switches to a different operational mode after 400 ms. It should be noted that the process of distancing the trajectory of the integer solution from the trajectory of the solution of the original algorithm is inevitable not only due to the limitations of the bit grid but also because of the different order of arithmetic operations in the solvers. However, the stability of the solution is retained, as was confirmed in a series of experiments.

Resonant Spike Generation
The simplified HH neuron model can demonstrate a stationary mode, damped subthreshold oscillations, and spike generation mode. Therefore, considering Izhikevich's classification [21], the HH model is a bistable resonator.
In Figure 4 the represented dynamical maps demonstrate the behavior of system (1), which was implemented using the solvers under investigation, when applying a pulse signal to its input with amplitude I = 4.5 µA and different values of period T and pulse width w. The period T ranged from 20 to 40 ms with step ∆T = 0.05, and the pulse width ranged from 10 to 30 ms with step ∆w = 0.05. The white areas on the diagrams of Figure 4 correspond to the stationary mode and the mode of subthreshold oscillations, while blue areas correspond to the spiking mode for the considered 32-bit and 64-bit FXP and DBL solvers built using the following methods: EE, SEE, EMP, and MEMP. The integration step size for all experiments was the same h = 0.01 ms.   The results showed that the excitability regions mainly form diagonal bands. The number of bands depends on the integration method used. The solvers based on the IE and SEE methods create the smallest excitability area: just three bands. The EE method leads to the highest number of bands in the excitability area, while the second-order methods EMP and MEMP give intermediate results.
The transition to an integer data representation reduces the excitability area for each of the considered solvers. Table 1 provides a numerical assessment of this phenomenon. The number of lost excitability areas for the 32-bit and 64-bit FXP solvers relative to the DBL solvers for given simulation parameters was calculated as a percentage. As one can see from Table 1, the FXP solvers of system (1) based on the EE method most strongly lost sensitivity to the pulse signal. However, these losses were calculated in thousandths of a percent.

Chaotic Spike Generation
A series of computational experiments were carried out for the case when dynamical system (1) goes into a chaotic mode of generating action potentials in response to a signal I = I 0 + A sin(2πωt). The interspike interval histograms were used as an analysis tool. They represent the dependence between the number of intervals n among spikes and their duration t. In [22], it was concluded that the application of interspike interval histograms to solvers based on first-order methods is impractical because they do not maintain the chaotic mode of the system. The diagrams for the DBL, 32-bit FXP, and 64-bit FXP solvers built using the methods EMP and MEMP are represented in Figure 5. Simulations were conducted on the time interval T = 105 ms with integration step h = 0.01 for all solvers under consideration. Figure 5 shows that the total number and duration of time intervals differ for all neuron models under consideration. However, it cannot be said that 32-bit FXP solvers are inferior to 64-bit FXP solvers in this series of experiments. It should be also noticed that differences in data type more strongly affect the solvers based on the EMP method: for the FXP64 EMP solver, intervals longer than 22 ms were not detected, although such intervals were observed on the prototype of the DBL solver.

Refractory Period
In this part of the research, we investigated the refractory period of the neuron model solvers. This period is the response time of a neuron model when two short-term current pulses are sent to its input. A common representation of input and output signals is provided in Figure 6. The interspike interval histograms were used as an analysis tool. They represent the dependence between the number of intervals n among spikes and their duration t. In [22], it was concluded that the application of interspike interval histograms to solvers based on first-order methods is impractical because they do not maintain the chaotic mode of the system. The diagrams for the DBL, 32-bit FXP, and 64-bit FXP solvers built using the methods EMP and MEMP are represented in  Figure 5 shows that the total number and duration of time intervals differ for all neuron models under consideration. However, it cannot be said that 32-bit FXP solvers are inferior to 64-bit FXP solvers in this series of experiments. It should be also noticed that differences in data type more strongly affect the solvers based on the EMP method: for the FXP64 EMP solver, intervals longer than 22 ms were not detected, although such intervals were observed on the prototype of the DBL solver.

Refractory Period
In this part of the research, we investigated the refractory period of the neuron model solvers. This period is the response time of a neuron model when two short-term current pulses are sent to its input. A common representation of input and output signals is provided in Figure 6. The time interval between the first and the second pulses ranged from 50 to 70 ms. The pulse width was 1 ms. Simulations were conducted on the time interval T = 100 ms with integration step h = 0.005 ms. The results provided in Figure 6 were obtained for the input pulse amplitude A = 20 µA. Figure 7а demonstrates that the occurrence times of spikes for system (1) almost coincide. The experiments showed that the data type of the EE, SEE, EMP, and MEMP solvers affects this characteristic significantly less than the numerical integration method that forms the basis of the solver. The maximum values of discrepancy in the time of occurrence of the second spike for the FXP The time interval between the first and the second pulses ranged from 50 to 70 ms. The pulse width was 1 ms. Simulations were conducted on the time interval T = 100 ms with integration step h = 0.005 ms. The results provided in Figure 6 were obtained for the input pulse amplitude A = 20 µA. Figure 7a demonstrates that the occurrence times of spikes for system (1) almost coincide. The experiments showed that the data type of the EE, SEE, EMP, and MEMP solvers affects this characteristic significantly less than the numerical integration method that forms the basis of the solver. The maximum values of discrepancy in the time of occurrence of the second spike for the FXP solvers compared to the DBL solvers based on the same numerical integration method are given in Table 2. Figure 7b represents the differences between the DBL solvers on the basis of EE, SEE, EMP, and MEMP. The difference in time of occurrence of the second spike here reached tenths of a millisecond. The time interval between the first and the second pulses ranged from 50 to 70 ms. The pulse width was 1 ms. Simulations were conducted on the time interval T = 100 ms with integration step h = 0.005 ms. The results provided in Figure 6 were obtained for the input pulse amplitude A = 20 µA. Figure 7а demonstrates that the occurrence times of spikes for system (1) almost coincide. The experiments showed that the data type of the EE, SEE, EMP, and MEMP solvers affects this characteristic significantly less than the numerical integration method that forms the basis of the solver. The maximum values of discrepancy in the time of occurrence of the second spike for the FXP solvers compared to the DBL solvers based on the same numerical integration method are given in Table 2. Figure 7b represents the differences between the DBL solvers on the basis of EE, SEE, EMP, and MEMP. The difference in time of occurrence of the second spike here reached tenths of a millisecond.   Table 2 demonstrates that all FXP solvers of system (1) give a similar discrepancy in occurrence time of the system (1) response regardless of the accuracy order of the integration method on which they are based.
In the case of input pulse amplitude A = 12 µA and the same simulation parameters, we obtained the results represented in Figure 8. and FXP solvers while sending two pulse signals with amplitude A = 20 µA. In the case of this amplitude, disappearance of the second spike at some intervals of time delay of the second input pulse was observed. The differences in behavior for the considered solvers of system (1) are retained. This can be seen in Figure 8b. However, compared to the previous experiment, the solvers on the basis of EE, SEE, and EMP match up to hundredths of a millisecond, and just the MEMP-based solver demonstrates a significantly different time of occurrence of the second spike on the system (1) output, calculated in tenths of a millisecond.

EE Solver SEE Solver EMP Solver MEMP Solver
The FXP solvers give a greater discrepancy with DBL solvers in terms of response time with the input signal amplitude A = 12 µA. Similar to in the previous experiment, this difference is less significant than the difference between solvers based on various integration methods. This is represented in Table 3.  In the case of this amplitude, disappearance of the second spike at some intervals of time delay of the second input pulse was observed. The differences in behavior for the considered solvers of system (1) are retained. This can be seen in Figure 8b. However, compared to the previous experiment, the solvers on the basis of EE, SEE, and EMP match up to hundredths of a millisecond, and just the MEMP-based solver demonstrates a significantly different time of occurrence of the second spike on the system (1) output, calculated in tenths of a millisecond.

Hysteresis
The FXP solvers give a greater discrepancy with DBL solvers in terms of response time with the input signal amplitude A = 12 µA. Similar to in the previous experiment, this difference is less significant than the difference between solvers based on various integration methods. This is represented in Table 3.

Hysteresis
In this part of the study, we explored the response of various neuron models to the half-period of the triangular wave. Such analysis allows for estimating the difference in hysteresis for system (1) solvers. Figure 9 illustrates an example of the input and output signals. The simulation was performed on the time interval T = 400 ms with integration step h = 0.005 ms. The period of the input signal was Т = 800 ms, and the amplitude of the input signal was A = 40 µA.
In this part of the study, we explored the response of various neuron models to the half-period of the triangular wave. Such analysis allows for estimating the difference in hysteresis for system (1) solvers. Figure 9 illustrates an example of the input and output signals. The simulation was performed on the time interval T = 400 ms with integration step h = 0.005 ms. The period of the input signal was Т = 800 ms, and the amplitude of the input signal was A = 40 µA. In order to evaluate the differences between the solvers in detail, the interspike interval histograms of system (1) were plotted. An example of these diagrams is shown in Figure 10. The time in milliseconds is plotted on the y axis and the number of interspike intervals on the x axis. The processing of the results obtained during these two experiments is summarized in Tables 4 and 5.  FXP32  25  25  22  30  FXP64  0  7  7  8   Table 5. The largest interspike interval errors between FXP and DBL solvers, ms. In order to evaluate the differences between the solvers in detail, the interspike interval histograms of system (1) were plotted. An example of these diagrams is shown in Figure 10. The time in milliseconds is plotted on the y axis and the number of interspike intervals on the x axis.

Data type EE Solver SEE Solver EMP Solver MEMP Solver
In this part of the study, we explored the response of various neuron models to the half-period of the triangular wave. Such analysis allows for estimating the difference in hysteresis for system (1) solvers. Figure 9 illustrates an example of the input and output signals. The simulation was performed on the time interval T = 400 ms with integration step h = 0.005 ms. The period of the input signal was Т = 800 ms, and the amplitude of the input signal was A = 40 µA. In order to evaluate the differences between the solvers in detail, the interspike interval histograms of system (1) were plotted. An example of these diagrams is shown in Figure 10. The time in milliseconds is plotted on the y axis and the number of interspike intervals on the x axis. The processing of the results obtained during these two experiments is summarized in Tables 4 and 5.  FXP32  25  25  22  30  FXP64  0  7  7  8   Table 5. The largest interspike interval errors between FXP and DBL solvers, ms. The processing of the results obtained during these two experiments is summarized in Tables 4 and 5.  According to the results of the hysteresis analysis, we can conclude that the MEMP method is the worst for fixed-point implementation. The EE method reproduces hysteresis properly even on 32 bits, but it poorly reproduces interspike intervals. The EMP method provides minor error compared to DBL for both bit lengths and can therefore be recommended for practical use.

Variation of the Integration Step Size
In discrete solvers of chaotic systems, the integration step size affects the system mode defining the rate of accumulation of rounding error. It leads to the emergence of a chaotic mode in points where it was not in the original system, and to the disappearance of this mode in the points where it exists in the prototype. To estimate the divergence of the DBL and FXP models using this criterion we built bifurcation diagrams that use the discretization step as a nonlinearity parameter. The charts of such dependency are called h-diagrams [23,24]. h-diagrams of the DBL and FXP solvers built on the basis of EMP are presented in Figure 11. Data type EE Solver SEE Solver EMP Solver MEMP Solver FXP32 0.117 0.003 0.05 0.12 FXP64 0.001 0.097 0.053 0.03 According to the results of the hysteresis analysis, we can conclude that the MEMP method is the worst for fixed-point implementation. The EE method reproduces hysteresis properly even on 32 bits, but it poorly reproduces interspike intervals. The EMP method provides minor error compared to DBL for both bit lengths and can therefore be recommended for practical use.

Variation of the Integration Step Size
In discrete solvers of chaotic systems, the integration step size affects the system mode defining the rate of accumulation of rounding error. It leads to the emergence of a chaotic mode in points where it was not in the original system, and to the disappearance of this mode in the points where it exists in the prototype. To estimate the divergence of the DBL and FXP models using this criterion we built bifurcation diagrams that use the discretization step as a nonlinearity parameter. The charts of such dependency are called h-diagrams [23,24]. h-diagrams of the DBL and FXP solvers built on the basis of EMP are presented in Figure 11. . Differences between the DBL solver and 64-bit FXP solver are insignificant and almost invisible. However, different from them, the 32-bit solver demonstrates chaotic behavior at small integration step sizes. The models based on the EE, SEE, and MEMP methods give close results. A more detailed demonstration of the difference in dependency of chaotic mode on integration step size for the researched models can be seen in Figure  12. The DBL and FXP solvers demonstrate chaotic model behavior for parameters specified in system (1) and values of integration step size h ∈ [0.001; 0.03]. Differences between the DBL solver and 64-bit FXP solver are insignificant and almost invisible. However, different from them, the 32-bit solver demonstrates chaotic behavior at small integration step sizes. The models based on the EE, SEE, and MEMP methods give close results. A more detailed demonstration of the difference in dependency of chaotic mode on integration step size for the researched models can be seen in Figure 12.  Figure 12 shows the dependency of the number of time intervals between spikes on integration step sizes. It can be noticed that methods of the second accuracy order (EMP, MEMP), compared to methods of the first order (EE, SEE), depend less on the integration step size and provide more stable solutions within the entire selected integration step size range. However, they are more sensitive to the data type used. On small integration step sizes ] 007 . 0 ; , a significant difference between the DBL solver and 32-bit FXP solver can be observed. This is not present in solvers built on the basis of first-order methods.

Discussion and Conclusions
Software models can be not applicable in many practical and scientific tasks because of their low performance and negative influence of the operating systems in real-time applications. In the last decade, scholars have intensively worked on high-speed realistic implementation of neuromorphic system in hardware. Meanwhile, these studies focused mainly on the interaction of neurons and the overall network architecture, while comparatively little attention has been paid to numerical solvers needed to synthesize finite-difference models reproducing single neuron dynamics. For example, in recent work by Pani et. al. [2], a 32-bit FPGA implementation of the Izhikevich neuron using the explicit Euler method was shown. A good correspondence between the fixed-point model and the floating-point simulation was established; this may confuse followers, as for the various neuron models, different network architecture and different bit-length FPGA simulation results may not be that close to the software simulation results. In our work, we focused on the dynamics of individual neurons and studied it in the case of various numerical methods and bit lengths to help developers correctly choose a numerical implementation.  Figure 12 shows the dependency of the number of time intervals between spikes on integration step sizes. It can be noticed that methods of the second accuracy order (EMP, MEMP), compared to methods of the first order (EE, SEE), depend less on the integration step size and provide more stable solutions within the entire selected integration step size range. However, they are more sensitive to the data type used. On small integration step sizes h ∈ [0.001; 0.007], a significant difference between the DBL solver and 32-bit FXP solver can be observed. This is not present in solvers built on the basis of first-order methods.

Discussion and Conclusions
Software models can be not applicable in many practical and scientific tasks because of their low performance and negative influence of the operating systems in real-time applications. In the last decade, scholars have intensively worked on high-speed realistic implementation of neuromorphic system in hardware. Meanwhile, these studies focused mainly on the interaction of neurons and the overall network architecture, while comparatively little attention has been paid to numerical solvers needed to synthesize finite-difference models reproducing single neuron dynamics. For example, in recent work by Pani et al. [2], a 32-bit FPGA implementation of the Izhikevich neuron using the explicit Euler method was shown. A good correspondence between the fixed-point model and the floating-point simulation was established; this may confuse followers, as for the various neuron models, different network architecture and different bit-length FPGA simulation results may not be that close to the software simulation results. In our work, we focused on the dynamics of individual neurons and studied it in the case of various numerical methods and bit lengths to help developers correctly choose a numerical implementation.
Our present study is devoted to the possibility of a realistic spiking neuron model implementation (the simplified Hodgkin-Huxley model) using fixed-point arithmetic with bit lengths of 32 (FXP32) and 64 (FXP64) bits, which is relevant when implementing neuromorphic systems on FPGAs. For numerical simulation of the neuron dynamics, four methods were chosen: the Explicit Euler method (EE), the Semi-Explicit Euler method (SEE), the Explicit Midpoint method (EMP), and the Modified Explicit Midpoint method with a smoothing step (MEMP). These methods were tested on various problems often found in the scientific literature on natural and artificial neuron studies. Therefore, the obtained results can be considered valid for a wide class of problems associated with neuron simulation.
To summarize the results of all the experiments, Table 6 is given. In the left column for each bit length, the method providing the smallest simulation error compared to the double (DBL) data type is placed. Methods listed in the right-hand columns are the runners-up. The following conclusions can be drawn from the results of the study: • With a 32-bit word length, the EMP method having second-order algebraic accuracy turned out to be the best solver. Among the first-order methods, the SEE method is preferable. The MEMP method seems to be worse for a low-bit implementation due to the higher number of arithmetical operations.

•
With a bit length of 64 bits, the MEMP method turned out to be the most accurate. The other three methods competed well with each other for accuracy in various tests. Taking into account its second-order algebraic accuracy, the EMP method is preferable. • Figure 12 shows that when using the EMP and MEMP methods with different integration steps, the dynamics of the neuron remains relatively unchanged, while the EE and SEE methods strongly affect the dynamics. The difference in the number of time intervals at various steps can reach 2 times, so the use of the Euler method and its modifications in neural network emulators cannot be recommended.

•
Summarizing the two last points, in the general case, one should choose the EMP method to reproduce the dynamics of neurons in fixed-point arithmetic.
In future studies, we will check intermediate bit lengths (such as 40 bits) to establish a smoother dependence of the method preference according to data type. The implementation of a neuron and neuromorphic system model on FPGA is also of interest. Funding: The reported study was supported by RFBR, research project No. 19-07-00496.

Conflicts of Interest:
The authors declare no conflict of interest.