Reliability of neuronal information conveyed by unreliable neuristor-based leaky integrate-and-fire neurons: a model study

We conducted simulations on the neuronal behavior of neuristor-based leaky integrate-and-fire (NLIF) neurons. The phase-plane analysis on the NLIF neuron highlights its spiking dynamics – determined by two nullclines conditional on the variables on the plane. Particular emphasis was placed on the operational noise arising from the variability of the threshold switching behavior in the neuron on each switching event. As a consequence, we found that the NLIF neuron exhibits a Poisson-like noise in spiking, delimiting the reliability of the information conveyed by individual NLIF neurons. To highlight neuronal information coding at a higher level, a population of noisy NLIF neurons was analyzed in regard to probability of successful information decoding given the Poisson-like noise of each neuron. The result demonstrates highly probable success in decoding in spite of large variability – due to the variability of the threshold switching behavior – of individual neurons.

Scientific RepoRts | 5:09776 | DOi: 10.1038/srep09776 To date, a great deal of efforts have been made to realize the different types of artificial hardware neurons, including the LIF neuron, using conventional complementary metal-oxide-semiconductor (CMOS) technologies 16,17 . This CMOS-based approach has been the mainstream approach of neuromorphic (hardware) engineering. A recent emerging research trend in neuromorphic engineering is the increasing adoption of alternative approaches to realize artificial neurons and synapses. These emerging approaches differ from the mainstream approach in that neural functionalities are implemented by introducing functional materials-based elements that could partly replace CMOS-based elements in the former 12 . One of the advantages of these new approaches is that it may enable the circuitry of ANNs to be substantially simplified by using a less number of CMOS elements than the conventional approach. An example of such approaches is a recent breakthrough by Pickett et al., who achieved the LIF neuron by using two pairs of a Mott insulator and a capacitor 18 . Basically, the neuristor concept introduced by Crane 19 in 1962 was employed in the LIF neuron model; hence this LIF neuron is termed as neuristor-based LIF (NLIF) neuron in this study so as to differentiate it from the standard LIF neuron.
Mott insulators are known to undergo temperature-driven insulator-to-metal transitions, so that the conductivity abruptly increases when the lattice temperature exceeds the transition temperature 20 . This transition is reversible, that is, the initial conductivity is recovered when the lattice temperature again falls below the threshold for the reverse transition, i.e., metal-to-insulator transition. As the lattice temperature change is due to Joule heating, the current-voltage (I-V) behavior of the Mott insulator is estimated to be volatile, i.e., threshold switching 18,21,22 . Threshold switching plays a key role in the functioning of the NLIF neuron 18,23 . Other than the Mott insulator, amorphous higher chalcogenides [24][25][26][27] , Si n + /p/n + junctions 28 , and particular transition metal oxides such as NbO x 23,29 are also known to exhibit volatile threshold switching.
Regardless of the type of threshold switch (TS) in the NLIF neuron, the variability of threshold switching behavior cannot be completely avoided 24,26 . It is, therefore, required to determine the effect of this variability on the neuronal behavior of the NLIF neuron that often leads to neuronal noises. The quantitative understanding of the noise is of significant importance when employing the NLIF neuron in both in silico and hardware-based ANNs because some noise properties are acceptable in ANNs insomuch as the noise does not cause serious errors during information processing 30 .
In this study, we conducted simulation on the NLIF neuron in order to identify its neuronal behavior -its noise characteristics and information representation in the presence of noise. Indeed, the employed analysis methods are widely used in characterization of other neuron models 1,3,31,32 , so that one can readily compare the characteristics of the NLIF neuron model with those of other models.
We first attempted to find optimized operational windows for the variables in the NLIF neuron model. It was indeed not an easy task to find the windows owing to the many variables involved simultaneously. We suggest a method to find successful spiking conditions by conducting static and dynamic calculations on the NLIF neuron circuit. The acquired windows should be narrowed down by taking into account the optimal selectivity of individual NLIF neurons for stimulation. The neuronal selectivity is one of the essential functions of neurons, given that they work as information encoders, in particular, in the presence of noise. Next, we allowed the variability of the TS in an individual NLIF neuron under optimal firing conditions. Such variability is most likely seen every switching cycle 24,26 -switching-event-driven variability. Threshold switching events repeatedly occur throughout the external stimulation duration, rendering the variability to feature a noise in the neuron's response. By analyzing the noise property, the relationship between the distribution and the consequent noise is understood and compared with the noise present in biological neurons. A question arising from the analysis on individual NLIF neurons is "Can conveying information, such as encoding and decoding, be achieved in a reliable manner by a population of these individual NLIF neurons?" This question is related to neuronal behavior at a higher dimension, i.e., the group, rather than at the individual neuronal level. The Bayesian decoder was employed so as to examine the reliability of the information conveyed by a population of NLIF neurons. As a result, the reliability was evaluated by means of "uncertainty. "

Results
Optimal firing conditions of individual NLIF neurons. In a single NLIF neuron circuit, standard circuit elements such as resistors (R 1 , R 2 , and R L ), capacitors (C 1 and C 2 ), and TSs (S 1 and S 2 ) are in use, as shown in Fig. 1a. 18 When it comes to a network of NLIF neurons, V 2 in Fig. 1a is relayed to a neighboring neuron through a synapse, so that V 2 works as the output voltage, corresponding to membrane potential. To generate a positive spike V 2 , V dc2 and V dc1 need to be positive and negative, respectively. For simplicity, it is assumed that V dc2 = −V dc1 = V d . The two dc voltage sources (V dc1 and V dc2 ) effectively supply power -enabling active operation -and determine the spiking dynamics including the spike's height and the undershoot level following the spike. The spiking dynamics will be explained in detail later. The key component in the NLIF neuron is the TS that performs monostable switching. The monostability of a TS can be understood from the schematic of current-voltage (I-V) hysteresis of the TS illustrated in Fig. 1b. The behavior of the TS is described by four parameters: R on , R off , V on , and V off , which denote the on-and off-state resistances and threshold voltages for off-to-on and on-to-off switching, respectively.
The assumption of linear I-V behaviors in both the states allows constant R on and R off in a given operational voltage range. For simplicity, switches S 1 and S 2 are assumed to be identical. The NLIF neuron fires spikes only when switch S 2 flickers at a given input current (I in ). Note that the term "flicker" means completing a threshold switching cycle, for instance, that along the arrows in Fig. 1b. To meet this requirement, the five standard circuit components (R 1 , R 2 , R L , C 1 , and C 2 ), the four operational parameters of the TS (R on , R off , V on , and V off ), and the dc voltage (V d ), i.e., ten variables in total, should be optimized. Owing to the difficulty in optimizing such a large number of variables, it is required to rule out several variables, in particular the operational parameters of the TS, that are most likely estimated from available experimental data 24,26,27,29 . In this calculation, R off and V on were set to 1 Mohm and 1 V, respectively, so that only two variables (R on and V off ) of the TS remain. They were converted to the following normalized variables: R off /R on and V off /V on . These threshold switching parameters are summarized in Table 1. Note that these ratios are often employed in characterizing resistive switching devices. A further reduction in the number of variables was made by setting R d and V d to 1 Gohm and 0.9 V, respectively. R L works as a voltage divider in this single neuron; it is desired to be large. V d needs to be close to, but smaller than, V on so as to turn on switch S 2 with a small input current I in ; 90 percent of V on , i.e., 0.9 V, was taken as V d .
To arrive at a condition of flickering switch S 2 at a given I in , time-independent calculations were performed with capacitors C 1 and C 2 ruled out (see Fig. 2a). The calculations provided I in and R 2 windows for spiking at given R off /R on and V off /V on values. The condition drawn from these static calculations is a "prerequisite" for successful spiking in the time domain. This is because the capacitors only determine the rate of voltage redistribution in the NLIF neuron upon switching of S 1 and S 2 , and the voltages across the two switches will eventually reach V on . Meeting the four requirements, shown in Fig. 2b and described below, allows S 2 to flicker. Note that on-switching of switch S 1 is a necessary condition for that of switch S 2 , but off-switching of switch S 1 is unnecessary for that of switch S 2 . Requirement i: setting R off for both switches in the circuit results in a voltage across switch S 1 (|V 1 +V d |) that is larger than V on , leading to the off-to-on switching of switch S 1 , given the aforementioned necessary condition for on-switching of switch S 2 . Requirement ii: setting R on and R off for switches S 1 and S 2 , respectively, results in a voltage across switch S 2 (|V 2 -V d |) that is larger than V on , leading to off-to-on switching of switch S 2 . Requirement iii: setting R off and R on for switches S 1 and S 2 , respectively, allows on-to-off switching of switch S 2 by decreasing |V 2 -V d | below V off . Requirement iv: setting R on for both switches allows on-to-off switching of switch S 2 regardless of on-or off-switching of switch S 1 , given that off-switching of switch S 1 is not a necessary condition for that of switch S 2 . Satisfying these requirements, R 2 windows with respect to input current I in and combinations of R off /R on (5, 10, 20, and 50) and V off /V on (0.3, 0.5, and 0.7) values were obtained as indicated using the grey zones in Fig. 2c. The white zones correspond to the failure of spiking. Insomuch as a current rather than a voltage is applied, R 1 and R 2 are independent variables as can be seen in Fig. 2c. It is noticed that the higher R off /R on and the lower V off /V on ratio are, the wider R 2  window is. For the following calculations, we chose moderate parameters of the TS (R off /R on = 20 and V off / V on = 0.5) and R 1 and R 2 of 100 kohm. As mentioned earlier, the windows drawn from the static calculations serve as necessary, rather than sufficient, conditions for successful firing of the NLIF neuron in a time domain; therefore, capacitors C 1 and C 2 need be optimized as well. The only concern in spiking in due course would be the sequential on-switching events of switches S 1 and S 2 when both are in the off-state, i.e., aforementioned Requirements i and ii are met in consecutive order. The major role of capacitors C 1 and C 2 in spiking is time-dependent redistribution of V 1 and V 2 upon switching of S 1 and S 2 . The capacitors determine the rate of the redistribution. That is, the higher the capacitance, the lower the rate. To satisfy the above-mentioned requirements, the time required for the evolution of V 2 -eventually leading to |V 2 -V d |≥V on , i.e., on-switching of switch S 2 upon the on-switching of switch S 1 -should be sufficiently short to hinder the off-switching of switch S 1 in the meantime. Otherwise, an increase in |V 2 -V d | in due course, owing to the on-switching of switch S 1 , would be abruptly diminished before V on is reached. Note that it was assumed that the switching times of S 1 and S 2 are sufficiently short to have a negligible impact on the time-dependent voltage redistribution. In addition, regarding the off-switching of switch S 2 , the static calculations basically assume the instability of the on-state of switch S 2 (see Requirements iii and iv), and thus off-switching occurs regardless of capacitances of C 1 and C 2 .
Given the above-mentioned requirements, a capacitance window for spiking at an input current of 1 μA is obtained as shown in Fig. 3a. The input current profile with respect to time is plotted in Fig. 3d. The maximum capacitance of C 2 for successful spiking at a given capacitance of C 1 tends to increase monotonically with that of C 1 . A higher capacitance of C 1 allows a longer discharging time of C 1 ; the discharging arises from the on-switching of switch S 1 and continues as far as |V 1 +V d | >V off , i.e., off-switching of switch S 1 . A higher capacitance of C 2 allows the charging time of C 2 to be longer; the charging arises from an increase in |V 2 | (V 2 < 0) occurring upon the prior on-switching of switch S 1 and continues until |V 2 -V d | ≥V on , i.e., on-switching of switch S 2 . Thus, a higher capacitance of C 1 enables the capacitance range of C 2 to widen, leading to the formation of the capacitance window shown in Fig. 3a. To highlight the capacitance dependence, four C 1 and C 2 pairs, denoted by α, β, γ, and δ in Fig. 3a, were sampled and for each pair a "membrane potential, " i.e., V 2 , the profile with respect to time was evaluated for the given input current I in . The results plotted are shown in Fig. 3b and 3c. In case of δ, the charging period of capacitor C 2 is longer than the discharging period of capacitor C 1 , and thus switch S 1 recovers its off-state before a transition of switch S 2 into the on-state. Therefore, no spiking is observed (see Fig. 3c).
To closely look at the evolution of V 1 , V 2 , R S1 (R of S 1 ), and R S2 (R of S 2 ) for case β, their time-dependent behaviors are zoomed in in Fig. 3e and 3f.
The spiking dynamics is described by the membrane potential V 2 and the auxiliary variable V 1 as follows: The NLIF neuron model is similar to the LIF neuron model regarding such that capacitor C 2 integrates potential and fires a spike when the threshold for the on-switching of S 2 is reached. However, a difference lies in the auxiliary variable V 1 in Eqs. (1) and (2). Given these two variables, the spiking dynamics can be mapped onto a V 2 -V 1 phase-plane, which is analogous to two-dimensional Hodgkin-Huxley neuron models such as FitzHugh-Nagumo model 33,34 . Eqs. (1) and (2) express V 1 -and V 2 -nullclines as the following linear equations: respectively. R S1 and R S2 are history-and voltage-dependent, and thus so are these nullclines. When R S1 = R off and R S1 = R on , the V 1 -nullcline is given by For case β, the evolution of V 1 , V 2 , R S1 , and R S2 is zoomed in in (e,f) to identify the self-consistent relation between them. The black and the red dashed line denote thresholds for on-and off-switching, respectively.
for case of R S2 = R off and R S2 = R on , respectively. When I in = 0, V 1 and V 2 stay at a stable fixed point (V 2 = 0.044, V 1 = -0.042) that is indicated in Fig. 4a. The parameters used in the phase analysis are listed in Table 2, which correspond to case β shown in Fig. 3. Upon the application of a constant current, this V 1 -nullcline in Eq. (5) shifts upwards by so that the fixed point moves far in the above-threshold region as shown in Fig. 4b. Consequently, the (V 2 , V 1 ) trajectory moves towards the fixed point in due course as seen in Fig. 4b until the on-switching condition for emerges at this moment, leading the trajectory to a new fixed point -again outside the sub-threshold region (see Fig. 4c). On its way, the off-switching condition for S 1 , (5)) is thus recovered, changing the direction again (see Fig. 4d). The path undergoes another change when it reaches the on-switching condition for S 2 , V 2 ≤ − |V on | + V dc2 = −0.1 V, as a consequence of emergence of the above-threshold V 2 -nullcline given by Eq. (8) (Fig. 4e). Encountering the off-switching condition for S 2 , ), so that the path heads to a new fixed point as shown in Fig. 4f. The subsequent spiking dynamics follows the limit cycle that is indicated using a grey line in Fig. 4f. The spiking dynamics of the NLIF neuron differs from that of the Hodgkin-Huxley neuron mainly in the fact that a stable fixed point varies upon V 1 and V 2 on the phase-plane, due to the V 1 -and V 2 -nullclines conditional on V 1 and V 2 .
Notably, the limit cycle is confined in the area on the phase-plane as seen in Fig. 4. That is, the spike's height and the level of the following undershoot are determined by V dc1 , V dc2 , V on , and V off , so that they are important parameters in spike's shape design.
Neuronal selectivity of the NLIF neuron. The circuit parameters of the NLIF neuron are required to be further optimized by taking into account the neuronal selectivity of the NLIF neuron for stimulation. As an information encoder, the NLIF neuron should be able to represent "distinguishable" responses to different stimuli. Neuronal responses are typically parameterized by the spiking rate or the spike number -also known as activity -in a given time period. Regarding the neuronal selectivity, the NLIF neuron needs to vary its firing rate upon input current I in . In particular, the firing rate and the input current are expected to be in a one-to-one correspondence relationship. Otherwise, one can hardly estimate the stimulus by counting the number of spikes, implying difficulty in "decoding" neuronal information. This difficulty in decoding consequently reduces the amount of information conveyed by the neuron 35 . In general, a neuronal encoding process is described by a = G[I in (s)], where a and s denote activity and stimulus, respectively. In this study, we define neuronal "activity" denoting the number of spikes in a time period of 30 ms. The function G in a biological neuron is nonlinear and exhibits a threshold value for activation, i.e., firing, of the neuron. This function is often referred to as the neuronal response function in which the activity is determined by input current I in .
For cases of aforementioned α, β, and γ, the neuronal response functions were simulated in the input current I in range (0-1.0 μA) and they are plotted in Fig. 5a. The neuronal response functions tend to increase monotonically with input current I in as long as the current is larger than a threshold of approximately 0.25 μA. This threshold results from a threshold voltage for the on-switching of switches S 1 and S 2 (V on ). All these functions appear to fulfill the aforementioned requirements for successful neuronal encoding. Nevertheless, a higher da/dI in value is favorable considering the fact that it reduces the uncertainty in discrimination when "noisy" neuronal responses are decoded. This issue will be revisited later when dealing with the noisy behavior of the NLIF neuron. Thus, case β appears to be most favorable and further discussion in this study will be narrowed down on this particular case. The corresponding parameters are listed in Table 2.
In biological neurons, the input current I in is understood to be determined by stimulus s, that is, I in is a function of stimulus s. Each individual neuron has a preferred stimulus s p , at which I in injected into a given neuron becomes the maximum. Note that this corresponds to the case of controlled neurophysiology experiments. The neuronal in-vivo function is, however, by and large triggered by an incident spike train(s) that provides a time-varying, rather than constant, synaptic current 1 . For simplicity, only one-dimensional stimuli are of concern in this study. For instance, stimulus s can be a one-dimensional visual stimulus such as the orientation angle of a light bar for the primary visual cortex 36,37 and a wind direction for the cricket cercal system 38 . For convenience, the orientation of a light bar is considered as a one-dimensional stimulus in this study. It is assumed that the input current is described by a Gaussian function whose maximum is placed at preferred stimulus s p as follows:    bell-shaped tuning curve appears consistent with that obtained for typical biological neurons. The response of the NLIF neuron shows its maximum at an orientation of 0°, corresponding to its preferred orientation, and tails around the preferred orientation. That is, stimuli within an orientation range of approximately −40° to 40° are able to activate the NLIF neuron although they are not exactly coincident with the preferred orientation. As a matter of fact, this imperfect-looking tuning curve enables a population of neurons with the limited number of preferred orientations to encode continuous, i.e., analog, information 31 . If neurons represented delta-function-like tuning curves, then an infinite number of such neurons would be required for encoding analog orientation information.
Noisy NLIF neuron. The tuning curve shown in Fig. 5c is of a perfectly working NLIF neuron. The orientation information encoded in the neuron can be decoded without uncertainty. Now, an arising question is "how large is the impact of imperfect behavior of the NLIF neuron on neuronal encoding and decoding?" Imperfect behavior is most likely caused by the variability of switching parameters of switches S 1 and S 2 (R on , R off , V on , and V off ). For actual experimental NLIF neurons, variations in such parameters cannot be avoided. Thus, an attempt to determine the quantitative uncertainty in processing neuronal information, caused by such variations, was made by evaluating the neuronal encoding and decoding processes with varying switching parameters. Firing each spike in a spike train involves off → on → off switching of each S 1 and S 2 in a consecutive order. Given that the switching parameters in an experimental switch, in general, varies on each switching cycle 24,26 , it is rather natural to assign different switching parameters to each switch immediately after each switching event. That is, repeating a random update on the parameters lasts throughout the entire spiking period. In this regard, the switching-event-driven randomness leads to time-varying variability, and thus "noise" rather than heterogeneity 39 . Such variations change the inter-spike interval (ISI) while a constant I in is applied for a given time period. Spiking behavior involving the variation is shown in Fig. 6. For this simulation, R on and R off of switches S 1 and S 2 were randomly sampled using Gaussian PDFs -centered at 50 kohm and 1 Mohm, respectively, with various deviations. After each of on-and off-switching events, a new resistance was assigned to the switches. The switching-event-driven update therefore lets R on and R off fluctuate in time as shown in Fig. 6b and d, implying a noise. It should be noted that the nullclines, Eqs. (5)- (8), are determined by R on and R off of S 1 and S 2 , and thus their variation essentially alters the nullclines and the corresponding fixed point. The random update, therefore, alters the trajectory of (V 2 , V 1 ) on the phase-plane. The spiking dynamics on the phase-plane and in the corresponding time-domain is shown in Fig. 6f and 6a, respectively. In this regard, the noise of the NLIF neuron is distinguished from other noise models for LIF models, e.g. diffusive noise given white noise and/or noisy synaptic current 40 , and for conductance-based model such as Hodgkin-Huxley neuron 41,42 .
The observed noise characteristics were quantitatively analyzed by examining a relationship between the mean and the variance of the spike number (activity) for a time period of 30 ms (see Fig. 6g). This and B≈1, where σ n 2 and n denote the variance and the mean of the activity, respectively 43 . This type of noise is referred to as the Poisson noise because the spike generation takes after the Poisson process 1 . The Poisson noise results in a straight line whose slope is unity as indicated with the dashed line in Fig. 6g. The relationship in Fig. 6g approximately features linearitysimilar to the Poisson noise -albeit steeper than unity (ca. 1.3). The NLIF neuron therefore exhibits a Poisson-like, rather than perfect Poisson, noise and the variance is by and large larger than the Poisson noise at given mean activities. Interestingly, for the 5 percent deviation case, the variance is much smaller than that of the Poisson neuron at mean activities of approximately 30. This is attributed to the activity limit by the capacitors' charging and discharging times that restrict the integration time for spiking (see Figure S1 in Supplementary Information).
The stochastic characteristics of this seemingly Poisson-like noise were further confirmed by analyzing the distribution of ISIs and the autocorrelation of the spikes in a given spike train; the results are shown in Fig. 6h and i, respectively. In Fig. 6h, the distribution is better fitted to a Gamma, rather than exponential, function regarding the effective refractory time caused by the finite recharging time of mainly capacitor C 1 1 . The evaluated autocorrelation data in Fig. 6i represent typical delta-function-like distribution, suggesting no correlation between the spikes. These noise analyses, therefore, identify Poisson-like noise characteristics of the observed noise.
To achieve successful spiking, the standard deviations of the switching parameters should be confined within particular ranges. Mostly, failure of spiking takes place when switch S 2 becomes stuck to its on-state and high membrane potential (V 2 ) is maintained. Typical examples of successful spiking and failure cases are seen in Fig. 7a and c, respectively. Insomuch as switch S 2 keeps its on-state in case of failure, the membrane potential remains high, so that no further switching of switch S 1 occurs. The current pulse duration was set to 30 ms. The number of successful spiking events was evaluated with separately varying the standard deviation of each switching parameter to determine the tolerance limit of each parameter (see Fig. 7e-h). We also varied input current I in (0.4, 0.6, 0.8, and 1.0 μA) in order to identify its effect on success in spiking. It turned out that the tolerance limit for R on and R off reaches up to approximately 30 percent, whereas the limit for V off is less than 20 percent as shown in Fig. 7e-h. It should be noted that the current pulse duration affects the probability of occurrence of an "R on -stuck" event. This happens because the longer the pulse lasts, the more spikes are likely to be emitted and the more likely that after one of these spikes switching parameters are set to a configuration that does not support further firing.
The switching parameters were chosen in a random manner by employing a Gaussian probability density function (PDF) with particular standard deviations. The means of these distributions were placed at the values used in the calculation of the perfect tuning curve in Fig. 5c (R on : 50 kohm, R off : 1 Mohm, V on : 1 V, and V off : 0.5 V, the other parameters are shown in Table 2). The encoding process of a noisy individual NLIF neuron was evaluated by calculating its tuning curve based on statistics. R on and R off simultaneously varied at different standard deviations and the mean activity at each orientation was obtained out of 100 trials. Given the very limited tolerance for V off variation as shown in Fig. 7h, no variation in V off was taken into account. Variation in V on was also ruled out in light of its negligible effect on successful operation probability (see Fig. 7f). The calculated tuning curves for four different deviations (5, 10, 20,   Fig. 8. The maximum activity tends to decrease with increasing standard deviation. This rapid decrease in activity results from the fact that higher deviation renders switching parameters more likely to settle into a configuration that does not support further spiking. Unexpected spiking also takes place, in particular, at orientations outside the active range (−40° -40°). Unlike the ideal tuning curve in Fig. 5c, a one-to-one correspondence relationship between orientation and activity is no longer satisfied. Therefore, a significant difficulty in decoding the neuronal information arises, consequently reducing the amount of information conveyed by the noisy neuron 35 . Representation of a population of NLIF neurons. The noise in the individual NLIF neuron seems to be an obstacle to appropriate neuronal information processing because of the difficulty in decoding caused by the noise. Fortunately, neuronal information processing in the brain does not strongly rely on individual neurons; instead, the task is in general performed by a population of individual neurons 2,31,35 . Nevertheless, the neuronal noise can still contaminate the population response. Some types of correlations between neurons in a population are known to reduce errors to some extent 35 , but this does not seem to be the general case. A possible answer to the question "how do brains as groups of unreliable (noisy) neurons work reliably?" is that populations of neurons may encode and decode "probability distributions" rather than particular values [44][45][46] . In other words, encoding and decoding are viewed as processes retuning probability distributions over all possible values: response and stimulus distributions for encoding and decoding, respectively. Especially, decoding is most likely based on a statistical inference process, in particular, Bayesian inference 32,45,47,48 . In fact, some psychophysical evidence for Bayesian inference have been found in, for instance, contrast-depending velocity perception 45,49 . Given the role of the NLIF neuron in either hardware-based or in silico systems, it is then an important task to examine the NLIF neuron as a Bayesian decoder, quantitatively evaluating probability distributions over the orientation at given degrees of variability of the switching parameters.
According where s and r denote stimulus and response, respectively. The notation P . Thus, we can evaluate the probability of population representation of a particular pattern → r when subject to a given stimulus s. We snapshotted → r patterns of the population of NLIF neurons with resistance deviations of 5, 10, 20, and 30 percent at a stimulus of 0°, and the patterns are plotted in Fig. 9a, b, c, and d, respectively. The different preferred orientations of the population let a few neurons preferring stimuli in the vicinity of 0° be activated despite the noise complicating the patterns.
Finally, the aforementioned Bayesian decoding was done for the patterns, leading the posterior PDFs shown in Fig. 9e, f, g, and h, respectively. The posterior PDFs are more or less noisy showing data scattering; the larger the deviation, the larger the data scattering. This data scattering is also a matter of population size, i.e., the more neurons in the population, the less decoding error. Given a large increase in calculation time with increasing the number of neurons in the population, we placed 20 neurons in the population; however, the Bayesian decoding of larger population sizes definitely enables correct estimation. The calculated posterior PDFs were fitted using Gaussian PDFs so as to evaluate the center and standard deviation of each posterior PDF. As can be seen in Fig. 9, the center of each PDF is found to be placed around 0°. This revealed that this Bayesian decoder most likely give a correct answer and a correct inference will be made if made by means of the Bayes' rule despite the present Poisson-like noise. Nevertheless, note that statistics cannot be free from error in any cases so that the Bayesian decoding can give a wrong answer at times. Besides, it turns out that the standard deviation in the decoding becomes larger as increasing the variability of R on and R off of the TSs (see Fig. 9i).
The Bayesian decoding results were compared with the case of populations of Poisson neurons. The likelihood PDF as well as the posterior PDF of a population of independent Poisson neurons is given by a closed-form expression; an increase in the number of Poisson neurons in the population leads to a posterior PDF of a Gaussian form 32 . At a given stimulus, the likelihood PDF of each independent Poisson neuron was analytically calculated with similar activity as that of the NLIF neuron at the same stimulus. As a result, the posterior PDF of a population of Poisson neurons with the four resistance deviations could be obtained from the population response patterns → r shown in Fig. 9a, b, c, and d; the calculated PDFs are plotted using dashed lines in Fig. 9e, f, g, and h, respectively. In fact, the 20 Poisson neurons already provide a Gaussian PDF as shown in the figures. In comparison with the population of Poisson neurons, it is noticed that the Poisson-like NLIF neurons represent smaller maxima and larger deviations than the Poisson neurons under the same condition except the case of 5 percent deviation of TSs' resistance (see Fig. 9a). A difference in the standard deviation of the posterior PDF between the NLIF neurons and the Poisson neurons is observed in Fig. 9i. The larger uncertainty deviation of the Bayesian decoding for the NLIF neurons arises from the larger deviation of activity of the NLIF neuron than that of a Poisson neuron as shown in Fig. 6g. Likewise, the larger maximum of posterior PDF of the NLIF neurons in Fig. 9e than the Poisson neurons is understood in terms of the smaller variance of activity at 5 percent deviation of TSs' resistance at the high mean activities, shown in Fig. 6g. The smaller variance is attributed to the activity limit by charging and discharging of the capacitors.

Discussion
The NLIF neuron studied in this work can serve as a prototypical in silico neuron model exhibiting a Poisson-like noise. The circuitry is simple and perhaps easy enough to be implemented in large-scale ANNs. In particular, as a result of this study, it is understood that the variability of the TSs' resistance leads to such a Poisson-like noise that the noise behavior of this prototypical in silico neuron needs to be under control and appropriately designed to meet the noise behavior required for ANNs built for specific purposes.
Nevertheless, when it comes to hardware realization of such NLIF neurons, there are several practical obstacles that should be overcome to realize the goal. What is of significant importance in the Poisson-like NLIF neuronal behavior is the minimum variability of V off of the TS in the NLIF neuron. As discussed earlier, the tolerance limit of V off is merely a few percent unlike that of the other switching parameters, i.e., R on , R off , and V on . Thus, "reliability" of this unreliable neuron requires meeting this stringent requirement for ensuring reliable operation. Apart from this restriction, other requirements discussed earlier may be satisfied by appropriate choices of TS materials, systems, and their design. Another important issue that potentially hinders practical use of this type of neuron is the long-term reliability of switches S 1 and S 2 , which are subject to the relatively high dc-voltage stress (V dc1 and V dc2 ). The dc voltages allow active operation of the NLIF neuron, working as effective power suppliers. Regarding the limit cycle confined in the area -(−|V on |+V dc2 ≤ V 2 ≤ −|V off | + V dc2 and |V off | + V dc1 ≤ V 1 ≤ |V on | + V dc1 ) on the phase-plane (see Fig. 4) -dc voltages close to, but smaller than, V on need to be applied to switches S 1 and S 2 , and thus the consequent electrical stress most likely affects the switches adversely. Eventually, it most likely leads to dielectric breakdown when a dielectric layer is in use as a TS material. This issue is also directly related to a high power consumption problem. The constant application of dc voltages during the lifetime of the NLIF neuron gives rise to severe power consumption, which is definitely against one of the inherent advantages of neuromorphic systems over standard digital systems, i.e., low power consumption. Therefore, addressing these significant problems properly accelerates practical use of such NLIF neurons in hardware-based neuromorphic systems.
The most crucial conclusion drawn from this study is that the potential variability of behavior of the TS is allowed up to a certain level as long as the Bayesian decoder is able to discriminate the encoded information correctly. In addition, the uncertainty, i.e., standard deviation, of the posterior PDF shrinks when introducing a larger number of NLIF neurons in the population. In general, the statistical accuracy of a survey increases with the number of samples. Thus, the uncertainty of posterior and likelihood of individual NLIF neurons is compensated by the increase in accuracy. An increase in the number of neurons in the population, therefore, tolerates a larger variability of switching parameters of the TSs. Nevertheless, confining the variability within a tolerance range is still of significant importance, especially confining that of V on .

Methods
Neuronal response function calculation. Eqs. (1) and (2)    The superscripts of V 1 and V 2 denote the ith node in the time domain. Eqs. (11) and (12) were numerically solved and R S1 and R S2 were updated when the evaluated potential reaches the threshold switching conditions.