Phase diagrams and dynamics of a computationally efficient map-based neuron model

We introduce a new map-based neuron model derived from the dynamical perceptron family that has the best compromise between computational efficiency, analytical tractability, reduced parameter space and many dynamical behaviors. We calculate bifurcation and phase diagrams analytically and computationally that underpins a rich repertoire of autonomous and excitable dynamical behaviors. We report the existence of a new regime of cardiac spikes corresponding to nonchaotic aperiodic behavior. We compare the features of our model to standard neuron models currently available in the literature.


Introduction
Modeling the brain is not a simple task. Usually scientists opt for simple models that provide insights about the original phenomenon one is trying to study. On the other hand, these models may lack some fundamental features that would play an important role in the considered phenomenon, specially when studying complex systems like the brain. Nowadays supercomputers have allowed the Neuroscience community to propose large-scale models for either brain functions or brain electrophysiology [1][2][3][4][5].
Most of these models represent each neuron by a point in space that performs every computation expected to happen through the whole extended body of a real neuron. Herz et al. [6] have grouped single-neuron models in five broad categories (from detailed compartmental models to statistical simple models), but all of the considered models are generally described by ordinary differential equations (ODE). The more complex the model, the more details on how is the shape and dynamics (both in space and time) of the action potential propagating in the neuron membrane. Although very simple, single-compartment neurons are capable of performing some fundamental tasks of information processing in single-neurons [6].
Map-based neuron models have emerged in the last two decades as a new class of models that have simpler equations and are computationally efficient [7,8]. Maps are discrete-time equations of continuous state variables, and map-based models describe the variation in time a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Models
The KT is a two-dimensional model defined by Eqs (1a) and (1b) below. The KTz is a tridimentional model defined by Eqs (1a)-(1c). In both models the gain function f(u) is a hyperbolic tangent. To obtain the KTLog and KTzLog models we approximate the hyperbolic tangent by its first order Taylor expansion, the logistic function f(u) = u/(1 + |u|): yðt þ 1Þ ¼ xðtÞ; zðt þ 1Þ ¼ ð1 À dÞzðtÞ À lðxðtÞ À x R Þ; ð1cÞ In Eqs (1a)-(1c) x(t) is the actual membrane potential (in arbitrary units) of the cell at time t (measured in arbitrary time steps, ts), y(t) is a recovery variable, z(t) is the slow total ionic current (which may generate bursts and cardiac spikes) and I(t) is an arbitrary potential generated by external currents (due to synapses and/or electrodes). The K and T parameters control the fast spiking dynamics, the parameter δ is the inverse recovery time of z(t) and controls the refractory period, and λ and x R adjust the slow spiking and bursting dynamics. Particularly, λ controls the damping of oscillations whereas x R controls bursting duration. H is a bias membrane potential and t is the discrete time step. Every variable/parameter is in arbitrary scale.
The membrane potential of the cell represented by x(t) in Eq (1a) may display a wide variety of autonomous and excitable behaviors. Some examples are cardiac spikes, bursting, slow or fast spiking, types I and II excitability, rebound spikes and bursts, nerve blocking effect, accommodation, and many others. These and more behaviors will be discussed in Section 3. The following subsections are dedicated to explore the bifurcations of the model. We first present the phase diagrams to provide an overview of the dynamics. Following, we explicitly provide closed forms for the fixed points (FP) and their eigenvalues to trace specific bifurcation diagrams.
We separate the presentation of the model in two cases in order to simplify its understanding: (I) the two-dimensional KTLog model: δ = λ = z(0) = 0; and (II) the tridimensional KTzLog model: H = 0 and any δ, λ and z(0). Case I is the fast subsystem and case II is equivalent to taking case I and transforming its H parameter into the dynamical slow variable z(t), H (case I) z(t), with arbitrary δ, λ and x R . Parameters K > 0 and T > 0 are assumed in this paper. The external current I(t) 6 ¼ 0 is used to study the excitable behaviors of the model, but is kept zero for tracing all the bifurcation and phase diagrams.

Phase diagrams
A phase diagram contains the stability limits of the model attractors in the parameter space. It gives a general picture of the frontiers between stable attractors of the model, be they FP or oscillatory attractors (OA). It is worth noticing that both cases of our map present three different types of OA: periodic, quasi-periodic and chaotic. Chaotic attractors are going to be discussed in Section 3.
A FP (x Ã , y Ã , z Ã ) is the point that does not change with map iteration, such that xðt þ 1Þ; yðt þ 1Þ; zðt þ 1Þ ½ ¼ xðtÞ; yðtÞ; zðtÞ A FP may also be referred to as a 1-cycle, because the map iteration exactly repeats after every 1 ts. Any periodic Q-cycle (Q > 1) may be determined similarly by . However, the Q-cycle regions in the phase diagram of case I are computationally calculated due to the variety of possible OA in the model. For a given FP (x Ã , y Ã , z Ã ), the stability limit is obtained using the condition max jL 1 j; jL 2 j; jL 3 j f g ¼ 1; where Λ i are the eigenvalues of the FP. Eq (3) yields an equation for x Ã over the stability limit. Therefore, the stability limit curve may be obtained by applying the x Ã thus calculated to Eq (2) and isolating the parameter of interest as function of the other parameters. Notice that for case I there are only two eigenvalues in Eq (3). Also, both y Ã and z Ã are linear functions of x Ã so it is sufficient to calculate x Ã to fully determine the FP. A given OA corresponds to a neuronal oscillatory behavior such as fast spiking (FS), bursting (BS) or cardiac spiking (CS). We use the interspike interval (ISI) distribution, PðISIÞ, to further distinguish between oscillatory behaviors for case II. We give details about the ISI distribution in the supporting S1 File. We also calculate the maximum amplitude, A, of the oscillations, where x(t) is a given OA, in order to identify pure subthreshold oscillations (SO) phases.

Case I.
The neuron has three parameters (K, T and H). We then trace diagrams K × T for H = 0 and then T × H for many K. We use the stability limit condition given in Eq (3) in order to derive the phase boundaries. The FP and its eigenvalues are given in the next subsection [Eq (15)] together with the description of each bifurcation.
There are three analytically determined stability limits for H = 0 (depicted in Fig 1A): T ¼ K; for K > 0:5; ð6Þ where Eq (5) separates a region of one stable fixed point (1FP) from two stable fixed points (2FP) via a supercritical pitchfork bifurcation, Eq (6) separates 1FP from an oscillatory attractor region (OA) by a supercritical Neimark-Sacker bifurcation and Eq (7) separates 2FP from OA region through a subcritical Neimark-Sacker bifurcation. The curve given by Eq (7) is the rightmost dashed line of the bistability (BI) region in Fig 1A. These bifurcations are going to be further discussed in the next section. The OA!2FP stability limit has been determined by map iteration and is the leftmost dashed curve of the BI region in Fig 1A. The BI region contains at least two kinds of stable attractors: OAs and FP. The OA region is filled up with Arnold tongues, each corresponding to a different periodic OA. The labels for these tongues in Fig 1A are winding numbers w P/Q, where Q is the attractor period (i.e. the amount of time steps the map takes to exactly repeat itself) and P is the amount of complete oscillations inside a period Q. The OAs in this region are born in a supercritical Neimark-Sacker bifurcation (the line T = K c ) over which |Λ + | = |Λ − | = 1. Thus, we may determine the labels simply by w The tongues are thinner than in the original hyperbolic tangent model [10] and they correspond to the periodic attractors in the bifurcation diagrams of the next subsection. The overall structure of the diagram is preserved by the logistic function approximation. Details about the algorithm used to determine the tongues are given by Tragtenberg & Yokoi [20].
Parameter H is a third axis going out of the paper in Fig 1A. Then for H 6 ¼ 0, we have three qualitatively distinct phase diagrams (see Fig 1B-1D) that correspond to the stability limits of Eqs (5), (6) and (7): K < 0.5; 0.5 K < 1 and K > 1. Vertical dashed lines in Fig 1A indicate K values used to trace the other three phase diagrams in the figure. The solid lines in Fig 1B-1D are the FP stability limits: derived from Eq (2), where x Ã are the FP values calculated from Eq (3): ; for K < 0:5; for K ! 0:5; We have two distinct solutions for Eqs (8) and (9) because the eigenvalues change from real to complex on K = 0.5. Eq (9) holds only over the bifurcation line, Eq (8). The full expression for the FP will be derived in the next subsection.
The curves defined by Eq (8) are vertical slices of the three-dimensional parameter space (perpendicular to K-axis and parallel to T-axis in Fig 1A). The dashed lines in Fig 1C and 1D are the stability limits of the OA phase and have been determined by map iteration. The OA phase is also full of Arnold tongues which are not shown for simplicity. The phase diagrams in Fig 1 give us a general picture of all the bifurcations of the model (some of which will be further detailed in the next subsection).
There is again a striking similarity between our approximation and the original model by Kinouchi & Tragtenberg [10] with two important differences: first, the two multi-critical points (MCP) from the original model coalesced into only one MCP (Fig 1C), sharpening the top-most part of the curves into a cusp-like form; second, the original supercritical Neimark-Sacker bifurcation has now turned into subcritical for K > 1 ( Fig 1D) and only the beak of this curve corresponds to a supercritical Neimark-Sacker bifurcation.
Comparing the diagrams in Fig 1, we conclude that the richest region in terms of dynamical behavior is within 0.5 < K < 1 because it presents all the phases of the model. Thus, we will study the 3-dimensional model only inside this range, even though other as interesting behaviors might be found for other values of K.
2.1.2 Case II. The variable z brings three new parameters to analyze (δ, λ and x R ). We trace the x R × T diagram for fixed δ = λ = 0.001 and K = 0.6 applying the interspike interval (ISI) method described in Methods section to determine the phases. It is possible to trace diagrams using the ISI method to any combination of parameters, but we will present only this case in order to compare with the diagram of the hyperbolic tangent model presented by Kuva et al. [11].   A and compare it with the same method applied to the tanh KTz. Panel B can also be compared to the diagram created using only map iteration in [11].
The FP stability limit may be easily calculated via an adiabatic approximation [considering that z(t) varies very slowly compared to x(t)]. It consists in taking H = z Ã = (x R −x Ã )λ/δ in case I and solving for the stability limits [Eqs (8) and (9)] in the plane x R × T [29]. We obtain (for 0 < δ = λ ( 1): Eq (10) is plotted as a dashed line in Fig 2A. This curve is way simpler than the one obtained with the hyperbolic tangent model (dashed line in Fig 2B) [29]. In principle, the stability limits of the KTzLog model could be precisely calculated for any δ and λ range, due to its analytical eigenvalues (see next subsection). The curve given by Eq (10) approximates [OðdÞ] the line of a supercritical Neimark-Sacker bifurcation between the FP and subthreshold oscillations (small amplitude OAs) in the plane x R × T.
Both diagrams present the same autonomous OA phases: cardiac spiking (CS), fast spiking (FS), bursting (BS) and subthreshold oscillations (SO). These behaviors are depicted in Fig 7I, 7G, 7J and 7M, respectivelly. The SO region is very narrow for the logistic case and is better identified by the amplitude of oscillations presented together with the bifurcation diagrams. Slow bursting (SB) has been separated from regular bursting using an arbitrary threshold only to emphasize that the interburst interval (IBI) diverges quickly near the bifurcation, as will be discussed ahead.
The ISI method provides an automated way to separate phases and also sheds light on structures and bifurcations which are very difficult to find via map iteration. Two features were identified using this method: the aperiodic cardiac spiking (ACS) region and the dust in the phase diagram between CS and BS. Both are also present in the original tanh version, but had not been previously identified using map iteration [11].
The ACS region is composed of aperiodic non-chaotic cardiac spiking attractors (see Fig A, panel D, in S1 File). The periods of these OAs are distributed within a broad skewed bellshaped distribution. The OA over the dust region between CS and BS is characterized by constant switching between CS and BS behavior (see Fig 7L and right-hand plot of Fig 6C). The dust region comprises very chaotic behaviors and might have a fractal shrimp structure in the plane x R × T, similarly to the parameter space of other oscillators [30,31]. This will be investigated in forthcoming works.

Case I.
Making use of the logistic function along with the FP definition, it is possible to write Eqs (1a) and (1b) as: where p T + |p 0 | and p 0 (1 − K)x Ã + H.
The eigenvalue, Λ, equation is: with solutions: The dependence of Λ on x Ã AE and H is implicitly inside p. These eigenvalues assume a complex value for K > 0.5, when H = 0, allowing for the appearance of OAs.
We plot the bifurcation diagrams ( Fig 3A-3H) versus each of the three parameters of the map, K, T or H. As usual, in these diagrams the stable FPs (|Λ ± | < 1) are plotted with solid lines whereas the unstable FPs (|Λ + | > 1 and |Λ − | > 1) are plotted with dotted lines. Saddle FPs (at least one of the eigenvalues has a different sign than the others) are plotted with dashed lines.
The gray regions in Fig 3A-3H correspond to OAs obtained by map iteration. The heights of the colored areas give the maximum amplitude of the OAs, calculated with Eq (4). These regions are filled up with many periodic attractors, or Q-cycles, windows. Most of these windows are too narrow to be perceived, so we do not show them. In Section 2.1, we trace the phase diagram of the model's case (I) with H = 0, clearly showing its periodic attractors.
The dark gray region in these bifurcation diagrams represent co-stability of OAs with at least one stable FP. These regions' internal borders then present subcritical Neimark-Sacker bifurcations (also known as Andronov-Hopf bifurcations for continuous time systems).
The bifurcation of x Ã over T parameter in Fig 3A and 3B changes from a pitchfork-like bifurcation in A for H = 0 to a cusp-like catastrophe in B when |H| > 0. In the case of H < 0 (H > 0), the bottom (top) branch gets disconnected from the top (bottom) branch in panel B. The behavior depicted in this figure is valid for 0.5 < K < 1. For 0 < K < 0.5, the OAs region vanishes, whereas for K > 1 the 2FP region vanishes. The behavior of the map is symmetrical in H. This fact is evident in the phase diagrams of Section 2.1. Fig 3C-3F shows the cusp catastrophe surface projection on the x Ã × H plan (two saddlenode bifurcations between unstable and saddle equilibria). The hysteresis curve in C becomes a single stable FP in F with increasing T, passing through a region of stable OA. In these panels, the FP loses stability via a subcritical Neimark-Sacker bifurcation, but the OAs that emerge in Fig 3C are only unstable and bifurcate by touching the saddle equilibrium as |H| increases. This figure is valid for 0.5 < K < 1. The OA regions vanish for K < 0.5 and, on the other hand, the 2FP region vanishes for K > 1. The value of the the parameters K, T and H at which the two FPs coalesce into one is easily obtained analytically using the stability diagrams which we shall present in Section 2.1. Fig 3G and 3H shows the bifurcation of x Ã according to K parameter for H = 0 and two different T. As T is increased, the OAs region decouples from the two stable FPs (2FP) region, changing from a subcritical [ Fig 3G] to a supercritical [ Fig 3H] Neimark-Sacker bifurcation for T ! 0.5. For H 6 ¼ 0, the symmetry is broken similarly to the bifurcation over T in Fig 3B, i.e. either the top (H > 0) or bottom (H < 0) FP survives and the other vanishes in a Saddle-Node bifurcation. Also, OAs exist only within a certain range of H close to zero and any |H| > 0 turns the supercritical Neimark-Sack bifurcation into subcritical. Again, there is symmetry in H ! −H.

Case II.
The FP equation for case II is easily derived from Eq (12), one should just replace H by z Ã = (x R −x Ã )λ/δ, resulting in the equation  (13) and (15) give stable FPs (-), unstable FPs (Á Á Á) and saddle FPs (---); and Eq (4)  where we define α λ/δ. The solution is simply: The FPs given in Eq (17) exist only if they satisfy the same conditions listed for Case I, replacing H by z Ã in the third condition, namely The Jacobian eigenvalues polynomial is given by: where the coefficients B, C and D are given: We plot the bifurcation diagrams of the model in Fig 4A-4F, with FPs given by Eq (17) and FPs stability given by Eq (19). Solid lines are stable FPs and dashed lines are saddle FPs. OAs again reside in the colored areas.
For Fig 4A-4C, the colored area contains the points x(t) of the actual OA attractor for each T iterated for 10 6 ts. If we had iterated the map for more time steps, the attractors would have more densely distributed points for each T. For Fig 4D-4F, the colored areas' heights give the amplitude of the OAs, Eq (4). BS lies in the dark gray region and FS or CS lie in the light gray region. The horizontal gaps are periodic attractor windows (corresponding to Arnold tongues in the phase diagrams of the previous subsection).
The z slow oscillation makes the 2FP region from the 2-dimensional model disappear, making it a single FP region. The values of x R , λ and δ determine the steady value of z Ã = (x R −x Ã )λ/ δ, which is equivalent to H in case I. Thus, the described behavior is symmetrical in x R , although the FPs would be positive.
In Fig 4A-4C, the FP loses stability via supercritical Neimark-Sacker bifurcation to, initially, SO. The 2FP region that is present in case I generates the CS region for case II. If one considers z varying very slowly (a quasi-static approximation), then a CS may be regarded as a slow switching between two stable FP. This statement can be verified by comparing Figs 3A, 3B and 4A-4C and noticing the high density of attractor points in CS region for case II where there were previously the two FPs of case I. The high density of the attractor close to the case I FP highlights the presence of a Saddle-Node bifurcation nearby (see Fig 3C) and is known as attractor ruins or ghosts [32,33].
In Fig 4D-4F, the FP also loses stability via a supercritical Neimark-Sacker bifurcation to a SO regime. The SO grows in amplitude according to a power law as K is increased and give rise to the BS region. The BS region shrinks as T increases. The OA amplitudes decrease as T increases. The decreasing of amplitude of OAs may also be explained by looking to case I: For H = 0 and T = K (K > 0.5) there is a supercritical Neimark-Sacker bifurcation in case I (see Figs 1A and 3A), therefore for |x R λ/δ| ≳ 0 and T ≲ K the amplitude of OAs should be small. In panel D, particularly, there is a large periodic window in the right-hand side of the OA region, displaying a 5-cycle attractor.
The FS and BS regions get extinguished as |x R | increases (see phase diagram in Fig 2A). The bifurcations' description presented so far for case II is valid for 0.5 < K < 1. For 0 < K < 0.5, there is only FP and CS behavior whereas for K > 1 there is no CS, but FS and FP. The same pattern repeats for x R > 0, but the oscillations occur around the symmetrical positive FP. The general picture of the bifurcations for all K values is better observed in Fig 1. We calculate ISIs in order to reveal the inner structure and the bifurcations of the OAs. We can clearly separate regimes of FP, FS, BS, CS and ACS (Figs 5C, 5D and 6). As far as we know, an ACS region has never been reported before. We may find narrow regions of chaotic Dynamics of an efficient map-based neuron model attractors inside the aperiodic regime, although the neuron is not generally chaotic inside the ACS region (see Fig 11C ahead). The aperiodic behavior appears due to failed attempts of spikes (see Fig 7N).
The ISI bifurcation diagrams are shown in Fig 5 (top blue curves). The bottom black points correspond to the amplitude, A, of the OAs as function of T (top panels) and x R (bottom panels). The amplitude highlights the presence of SOs, which appear through a supercritical Neimark-Sacker bifurcation of the FP. Except for the transition between SO and ACS, the amplitude continuously varies for all the considered bifurcations in Fig 5; the apparent discontinuity is due to the a fast variation in A (within a range of 10 −10 of T or x R ). The appearance of large amplitude OA and the qualitative change in oscillations are marked with vertical dotted lines. Fig 5A-5C correspond to approximately the same parameters as those in Fig 4A-4C, respectively.
The FS regime coexists with CS regime in a small range of parameters T and x R . Thus, there is a first-order-like transition between these two regimes (see Fig 5A). The transition from BS to FS, on the other hand, happens through increasing duration of a single burst of spikes, Δt BS , similarly to what happens for the leech heart interneuron modeled by Shilnikov & Cymbalyuk [34] using a conductance-based approach. Fig 5E shows the ISI profile of this transition whereas Fig 6A shows Δt BS close to the transition point x BS R % À 0:075 for the same set of parameters. Notice that as the interburst interval (largest ISI) stays nearly steady as x R ! x BS R , burst duration increases boundlessly according to hDt BS i $ jx R À x BS R j À 1=2 characterizing a Homoclinic bifurcation of a Saddle-Node periodic orbit [34]. The border between CS and BS ( Fig 5B) is filled up with chaotic attractors due to the loss of stability of cardiac spikes (see e.g. Fig 11A and 11B). In Fig 6C, we show the behavior of the membrane potential, x(t) versus t from Eq (1a), as the transition CS!BS is approached. For  (4), has discontinuities at every bifurcation. These discontinuities tend to disappear as |x R |!0 because then the model approaches a supercritical Neimark-Sacker bifurcation. The amplitude also highlights the presence of subthreshold oscillations in the transition from FP to BS or to ACS. https://doi.org/10.1371/journal.pone.0174621.g005 Dynamics of an efficient map-based neuron model small T (left plot in Fig 6C), the cardiac spikes are very regular resembling the FitzHugh-Nagumo model [35,36]. For increasing T, the OA loses stability and displays small oscillations even during the slow decrease of the action potential (middle plot). These small oscillations are caused by the subcritial Neimark-Sacker bifurcation of the fast subsystem (Fig 3C), similarly to what happens in early afterdepolarizations of the action potentials of conductance-based models [28,37]. The BS region starts after a very chaotic region of mixed behavior of cardiac spikes and bursts (right plot) taking place for 0.24 ≲ T ≲ 0.25 (the dust region in Fig 2A).
The transition ACS!SO happens through an infinite period bifurcation (Fig 5C and 5D).
Here, x CS R is the bifurcation point. The same scaling law holds for hISIi versus the T parameter, around the bifurcation point T CS . The variance of ISI diverges close to the transition point ðT CS ; x CS R Þ as a power law VarðISIÞ $ jx R À x CS R j À 1:8 (inset of Fig 6B).
The bifurcation that separates BS and SO (Fig 5B and 5E) is similar to that separating ACS and SO. The largest ISI (or the interburst interval, IBI) diverges according to hIBIi $ jx R À x SO R j À 1=2 and the smaller ISI remains constant. The same power law is found over the T parameter.

Fig 6. Infinite period bifurcation, blue-sky catastrophe and CS!BS transition. Panel A:
Burst duration, Δt BS , as a function of jx R À x BS R j, where x BS R is the bifurcation point in which fast spiking begins through a blue-sky catastrophe [34]. The inverse square root is only valid very close to the bifurcation. Panel B: Period of oscillations, ISI j , as function of jx R À x CS R j, where x CS R is the bifurcation point given by Eq (10) in the quasi-static approximation. Notice how hISIi goes to infinity as the system gets closer to the bifurcation according to the inverse square root. This figure presents the same data as in Fig 5D. A small beak near jx R À x CS R j % 10 À 1 separates aperiodic (left) from periodic (right) cardiac spiking. Panel B inset: variance of ISI, a power law with slope −1.8. Panel C: Transition from CS to BS for increasing T and fixed K = 0.6, δ = λ = 0.001 and x R = −0.2. Notice as slow oscillations lose stability going through a behavior of mixed cardiac spikes and bursting. The shape of action potential in panel C, middle, resembles that of a early depolarized action potential [28]. https://doi.org/10.1371/journal.pone.0174621.g006 Dynamics of an efficient map-based neuron model

Behaviors
This section is dedicated to show that this map exhibits many dynamical behaviors of excitable systems, like neuronal cells and muscle cells, especially the action potential of cardiac cells [38].
In this work, the solution x(t) of the map, Eq (1a), is the actual membrane potential of an excitable cell. Figs 7 and 8 will have x (in arbitrary units) in the vertical axis (plotted with circles and lines) and t (in time steps) in the horizontal axis. The solution of Eq (1a) is the set of points plotted in these figures; the lines are there only as a guide to the eyes. Excitable behavior is accompanied by the curve of the input current, I(t) (in arbitrary units), right below the membrane potential curve. The solution x(t) corresponds only to circles, but lines are drawn connecting the points in order to guide the eyes of the reader. Parameters for every behavior in both figures are given in Table A in S1 File.

Case I
The bidimensional model can only present three autonomous behaviors: resting state, fast spiking and subthreshold oscillations. Fast spiking in this model is sometimes chaotic. The excitable behaviors in turn are more diverse and are presented in Fig 7A-7F. The excitable behaviors are generally found near the subcritical Neimark-Sacker bifurcation.
The observed excitable behaviors are: relaxation to FP, transient oscillations, tonic spiking, nerve blocking, bistability and excitable spiking. All of them are also present in the original hyperbolic tangent model [10]. However, the amplitude of the logistic model oscillations is somewhat smaller than the original, because the logistic function varies slower for |T| ! 0 than the hyperbolic tangent.  KT model presents true excitability because, in most bifurcations, a saddle equilibrium separates the basins of attractions of both FPs (see bifurcation diagrams in Section 2). Each curve in Fig 7F is the result of the neuron receiving a small delta pulse stimulus in the first time step of the simulation. The excitability threshold is easily found by simply taking the distance (in the parameter space) between the simulated neuron and the stability limit of the FP calculated in Section 2.

Case II
In addition to the autonomous behaviors of case I, KTz neuron also presents: slow spiking, fast and slow cardiac spiking, aperiodic cardiac spiking, fast and slow bursting and chaotic bursting (mixed cardiac spiking and bursting found in the boundary between these two regimes, as shown in both Figs 7L and 6C center plot).
Cardiac spikes are fully described only by other two differential equation models, namely: the Luo & Rudy model [39] and the FitzHugh-Nagumo model [35,36]. The first map-based cardiac spike was obtained using the hyperbolic tangent KTz model [8,11]. The logistic simplification allows us to precisely visualize the mechanism of generating this kind of spikes in maps, differently from other approaches that involve solving the complex differential equations of Luo & Rudy using Euler method with unity step [40].
The existence of the dynamical slow variable z generates a large number of slow-fast dynamical behaviors. They are depicted in Fig 8. We quickly explain these behaviors following Izhikevich [21] in order to be able to compare our map-based approach computational efficiency with other models in Section 4.
Tonic spiking: Neurons usually present a quiescent behavior until excited by a continuous current, when they start firing spikes periodically (Fig 8A). The train of spikes is present only while stimulated by the external current. This behavior is found in almost every known neuron.  Table A in S1 File. https://doi.org/10.1371/journal.pone.0174621.g008 Dynamics of an efficient map-based neuron model Phasic spiking: The neuron fires only once at the beginning of stimulation then becomes quiescent. This is known as phasic spiking (Fig 8B) and is present in general excitable elements.
Tonic bursting: KTz neurons are also able to present tonic bursting (Fig 8C) while the continuous current is active. This feature is important since it is related to the production of the gamma oscillations in the brain [21].
Phasic bursting: Similarly to phasic spiking, the neuron reacts to the stimuli by a unique burst of spikes (Fig 8D).
Mixed mode: Some neurons can exhibit a mixed behavior: they fire a phasic burst spike then start tonic spiking (Fig 8E). Our case II neuron is capable of representing this spiking mode.
Class 1 excitability: Some neurons present a variable spiking frequency that depends on the magnitude of the injected current [33]. These neurons are said to display a Class 1 excitability ( Fig 8F) and are able to code its input in an unique frequency output. Such behavior may be found in cortical pyramidal neurons and is typically caused by a Saddle-Node on an invariant circle bifurcation [33].
Class 2 excitability: There is another class of excitability in which neurons are not capable of firing arbitrarily high frequency spikes (Fig 8G). This behavior is mostly independent of the intensity of the applied current and is typically found in cortical inhibitory interneurons. These spikes are generally caused either by a Saddle-Node bifurcation or by a Neimark-Sacker bifurcation [33].
Subthreshold oscillations: The presence of oscillations that are insufficient to trigger a spike is common in a large variety of neurons (Fig 8H). These small oscillations may happen autonomously or after a spike, as pictured in this panel.
Resonator: Some neurons can act as resonators. They react only to stimuli in which the frequency is the same as their subthreshold oscillations frequency. In Fig 8I the neuron does not spike when the stimulus frequency is too high.
Integrator: In a neuron without subthreshold oscillations, the higher the input frequency, the higher the probability of spiking. This kind of neuron is called integrator (Fig 8J).
Rebound spike: Some neurons can fire after an inhibitory input is received. This phenomenon is called post-inhibitory, or rebound, spike (Fig 8K).
Rebound burst: Instead of single-peak spike some neurons can fire a burst of spikes after receiving inhibitory stimuli (Fig 8L).
Threshold variability: Our model is capable of reproducing the phenomenon of threshold variability (Fig 8M). A stimulus which otherwise would not be able to generate a spike, indeed generates a spike after a small perturbation of the local membrane potential.
Bistability: Some neocortical neurons can present two different stable modes. This can be seen in Fig 8N, in which a stable FP and a stable OA coexist. An external stimulus changes the relative positions of the nullclines momentarily and leaves the neuron either in the FP or in the OA basin of attraction.
Accommodation: When receiving a slowly increasing continuous input, certain neurons fail to fire. This indicates an accommodation to the stimulus, as shown in Fig 8O. If the same neuron is stimulated by a small and brief input it is able to fire a spike.

Chaotic attractors
We use the Eckmann-Ruelle method [41] to calculate the largest Lyapunov exponent, λ L , and find strange attractors. We can find chaotic regions in parameter space for the KTLog model, although these regions are narrower and displaced when compared to the KT model. Typical Lyapunov exponents in the logistic approximation are also smaller than the hyperbolic tangent model. A typical chaotic attractor for the case I and H = 0 is in Fig 9. Its capacity dimension is 1.35 (2) and the Lyapunov dimension [42] is 1.158(2), thus categorizing it as a strange attractor.
We computationally check the chaotic behavior by iterating the map for two sets of initial conditions very close to each other (*10 −8 ) yielding x 1 (t) and x 2 (t) and Δx(t) = |x 2 (t)−x 1 (t)|. Only a few time steps are needed for the separation between the curves to become visible (Fig  10 top panels). We plot Δx versus t in Fig 10 bottom  We also applied this method to determine the largest Lyapunov exponent for a vast region of parameters of case II. We show λ L as function of T in Fig 11A-11C for the same set of parameters of Fig 5A-5C, respectively.
In fact, Fig 11 shows the Lyapunov exponent over the horizontal dotted lines of the x R × T diagram in Fig 2. Notice that chaotic behavior is very strong on the frontier between CS and BS and inside BS region. Also, positive Lyapunov exponent has been found on the edges of ACS region (Fig 11C). Dynamics of an efficient map-based neuron model

Computational efficiency
Usually, the performance of a program depends on many uncontrollable variables, such as code style, language specific implementations, background memory operations and usage, function calls, concurrent processing due to operating system threads, etc. Thus, there is a diversity of measures for computational efficiency. Izhikevich [21] defines it through the amount of floating-point operations (FLOPs) the neuron model needs in order to evolve 1 ms (model time) of its dynamics, considering that a spike takes around 1 ms to rise and fall.
The intrinsic dynamics of the model also matters for performance: after a spike, each model reaches the FP after different number of time steps. For instance, the KTzLog model takes only 136 ts (roughly 1,632 FLOPs) to reach the FP whereas the Izhikevich model takes 1,158 ts (roughly 15,054 FLOPs). These convergence times clearly depend on the specific phase portrait of the model for the chosen parameters. Thus, we follow a similar approach to Girardi-Schappo et al. [8] and measure the CPU cycles needed for a single time step of the model to be   Table 1 for the more popular models. A CPU cycle is the fundamental operation cycle of a computer and thus it takes into consideration all the processing operations. It consists of retrieving a program instruction from the memory, determining the needed actions to process the instruction and executing these actions [43].
The number of FLOPs/ts of both logistic and hyperbolic KTz models are 12 FLOPs/ts and 31 FLOPs/ts, respectively. We show them in Fig 12 together with other standard models and their quantity of biological features. We estimated the quantity of biological features of our neuron model through Figs 7 and 8. Izhikevich [21] defines a set of 22 biological features in order to measure "biological relevance" of the model. However, notice that some features are not accounted for in the proposed list, such as cardiac spikes, nerve blocking and transient oscillations. Out of that list, we found 15 biological features in both logistic and hyperbolic KTz plus these three aforementioned behaviors. It does not mean that the remaining 7 features of Izhikevich's list are not present in KTz family, because there is a lot more to explore in our proposed neuron map. Table 1. Computational efficiency of some standard single-neuron models. The ts to FP column is the number of time steps taken to approach the FP with a tolerance of 10 −8 . The ms/ts column shows the proposed scale factors for the time step of each model, assuming a spikes takes about 1 ms. The number of FLOPs for the conductance based model is estimated assuming a 4 th order Runge-Kutta ODE solver with step 0.001 ms. Each model's average CPU cycles has been performed over 5,000 realizations of 3,000 ms considering the lowest conversion factor in the table for the considered model. Dynamics of an efficient map-based neuron model

Model
Map-based models are seven to twenty times more efficient than ordinary differential equation (ODE) models if we consider the CPU cycles in Table 1. Even though Rulkov, Izhikevich and our model have similar number of FLOPs/ts, our model time step takes about one and a half times more CPU cycles to be evaluated than the other two models. These extra cycles are needed in order to deal with possible overflow operations for large absolute arguments of the logistic function during map evaluation. Notice however that the logistic approximation optimized both the FP convergence time (in ts) and the quantity of CPU cycles per time step in relation to the hyperbolic model. And, more interestingly, the logistic approximation kept all the dynamical behaviors of the hyperbolic KTz.
We also consider two typical network cases: linear and mean-field. In both cases, there are N neurons and the synaptic current over an element i is given by a gap junction model [8]: where the sum runs over the neighbors of i and G > 0 is an excitatory conductance that is adjusted for each case slightly above the necessary threshold to propagate activity. The coupling scheme for each map-based model is described in details elsewhere [8]. The linear network is directional (with a total of N − 1 synapses), neurons are adjusted in an excitable regime and only the left-most neuron is initially excited (i = 1). The signal propagates until it eventually reaches neuron i = N. The mean-field network is a complete graph without self-interactions [total of N(N − 1) synapses], and random initial conditions, so the neurons (adjusted in a bursting regime) tend to synchronize after some transient activity. The conductance-based model is either the standard Hodgkin-Huxley model for the linear network or the leech heart interneuron model given by [34] due to its bursting regime. Fig 13 shows the CPU cycles/ts for each network of each model. Here, the network time step is evaluated as follows: 1) calculate the signal of every synapse; 2) for each neuron, sum up its input signal; 3) calculate the new membrane potential of each neuron considering its input current.
Both networks have similar trends as N is increased. The number of CPU cycles/ts increases proportionally to the number of synapses in the network for large N (notice the black solid Dynamics of an efficient map-based neuron model curve in both panels of Fig 13). The conductance-based models demand significantly more computational power. Except for the tanh KTz, the KTzLog, Izhikevich and Rulkov maps have equivalent efficiency. The large fluctuations for small N appear because the background operating system processes have a high impact on the CPU cycles of fast jobs. It is worth noticing that the precise number of CPU cycles strongly depends on the processor, memory architecture and operating system used to compile and run the program. However, the statistical trend should still hold for different computers. We used Microsoft(R) Windows(TM) 10, version 10.0.14393 Build 14393, on a x64 Intel(R) Core(TM) i7-4500U CPU @ 1.80GHz, 2401 Mhz, 2 Cores (4 Logical Processors), using 16 GB of DDR3-RAM.

Concluding remarks
We studied a model of action potential generation using difference equations obtained through the first order Taylor approximation of the hyperbolic tangent KTz neuron model [11]. We presented detailed bifurcation diagrams and fixed point stability diagrams for our new model, the so-called KTzLog neuron. We have shown that the KTzLog neuron reproduces many neuronal excitatory and autonomous behaviors observed experimentally. We have compared our model's computational performance to other widely used models and concluded that its efficiency is comparable to that of the most efficient neuron models, both in isolation or in a network.
We highlight the presence of cardiac spikes within a large parameter region of the KTzLog model. This behavior is only possible because we do not follow an integrate-and-fire procedure for spike generation. Instead, the action potential of our model has its own slow-fast dynamics. Other authors have attempted cardiac spikes modeling, either by discretizing Luo & Rudy model [40], or by proposing piecewise-continuous functions [44]. Or yet, some authors proposed the modeling of the action potential duration through a difference equation [45]. However, only the KTz map family displays cardiac spikes with a reduced set of parameters on top of a simple continuous sigmoid dynamics. This system was fundamentally conceived as a map [8] and the logistic approximation allows us to analytically study bifurcations and stability of orbits without loss of any dynamical behavior. This fact may allow us to better understand phenomena such as early afterdepolarization and cardiac arrhythmia: we have identified a bifurcation that is known to cause early afterdepolarized action potentials in ODE neurons [37].
An aperiodic cardiac spike phase was found in a large parameter region of both versions of KTz model (logistic and tanh). Surprisingly, the Lyapunov exponents inside this phase are negative, indicating a non-chaotic behavior. We are not aware of another work that has shown the existence of this non-chaotic and aperiodic behavior in excitable systems. We will investigate these subjects in forthcoming works.
Mesbah et al. [16] have proposed a uni-dimensional map-based neuron model that displays a few excitable behaviors. However, such behaviors are obtained by applying external currents to different parameters, making it very hard for neurons to be coupled in order to produce an heterogeneous neuronal network (as the brain is expected to be). On the other hand, external input currents and synaptic coupling is easily attainable by adding a current term to the KTzLog model equation similarly to what is done with conductance based models. An example of synaptic map was proposed by Kuva et al. [11], and used by Girardi-Schappo et al. [22] to study neuronal avalanches. Thus, it is straightforward to build biologically motivated networks of KTz maps in order to study higher brain functions, such as cortical processing, synchronization phenomena and so on.
KTzLog map has only five parameters (or six, if H is considered a polarizing current for the three-dimensional model) displaying almost as many excitable behaviors as the Izhikevich [14] model (which has nine parameters). In this sense, the KTz family, and specially its logistic variation, provides minimal models to study neuronal bifurcations. Also, the logistic case displays all of its bifurcations in their normal forms, making the KTzLog model one of the canonical neuron models of choice for studying unknown mechanisms underlying dynamical phenomena in neurons (such as early afterdepolarization).
Modeling complex functions of the nervous system demands an extensive gathering of information about the system one wants to mimic. During this process, a certain degree of simplification is unavoidable. Neurons are the basic processing units of the brain, but they are generally the first elements of the simulation that are reduced. However, these simplifications are sometimes useful, as they may result in a more reliable model in the context of the function one is trying to understand [6,46]. On the other hand, modelers must use models that allow a certain degree of controllability of its neurons' dynamics in order to be able to make inferences about the simulated experimental setup. In any case, the chosen neuron model should be capable to provide simultaneously a good computational performance and enough features to make the model reliable and easy to use. We have shown that the KTzLog model meet these essential requirements.
Furthermore, future work with our model include understanding the chaotic attractors and their coexistence with periodic attractors, the shrimp structure in the x R × T diagram, compartmental modeling, cardiac spike bifurcations and cardiac tissue modeling.
Supporting information S1 File. Supporting information: Phase diagrams and dynamics of a computationally efficient map-based neuron model. Details about the ISI method used to determine the OA in this paper and about the parameters of the model for each behavior depicted in Figs 7 and 8. Fig A, Typical ISI distributions. The four different types of ISI distribution PðISIÞ are displayed in panels A to D (top) with their corresponding map iteration (bottom). Solid and dashed lines are only there to guide the eyes. Panel A: Fast spiking (FS)-a single well defined peak in PðISIÞ such that hISIi < ISI th . Panel B: Bursting (BS)-two peaks are generally present in PðISIÞ for BS behavior; if chaotic bursting is present, both peaks will be broadened; slow bursting phase has the large ISI larger than an arbitrary threshold. Panel C: Periodic cardiac spiking (CS)-a single peak in PðISIÞ such that hISIi > ISI th . Panel D: Aperiodic cardiac spiking (ACS)-a single peak broad distribution PðISIÞ shaped similarly to a lognormal curve. Fixed point (FP) has no ISI. Table A