A discrete Huber-Braun neuron model: from nodal properties to network performance

Many of the well-known neuron models are continuous time systems with complex mathematical definitions. Literatures have shown that a discrete mathematical model can effectively replicate the complete dynamical behaviour of a neuron with much reduced complexity. Hence, we propose a new discrete neuron model derived from the Huber-Braun neuron with two additional slow and subthreshold currents alongside the ion channel currents. We have also introduced temperature dependent ion channels to study its effects on the firing pattern of the neuron. With bifurcation and Lyapunov exponents we showed the chaotic and periodic regions of the discrete model. Further to study the complexity of the neuron model, we have used the sample entropy algorithm. Though the individual neuron analysis gives us an idea about the dynamical properties, it’s the collective behaviour which decides the overall behavioural pattern of the neuron. Hence, we investigate the spatiotemporal behaviour of the discrete neuron model in single- and two-layer network. We have considered obstacle as an important factor which changes the excitability of the neurons in the network. Literatures have shown that spiral waves can play a positive role in breaking through quiescent areas of the brain as a pacemaker by creating a coherence resonance behaviour. Hence, we are interested in studying the induced spiral waves in the network. In this condition when an obstacle is introduced the wave propagation is disturbed and we could see multiple wave re-entry and spiral waves. In a two-layer network when the obstacle is considered only in one layer and stimulus applied to the layer having the obstacle, the wave re-entry is seen in both the layer though the other layer is not exposed to obstacle. But when both the layers are inserted with an obstacle and stimuli also applied to the layers, they behave like independent layers with no coupling effect. In a two-layer network, stimulus play an important role in spatiotemporal dynamics of the network.


Introduction
Investigating the local dynamics of neurons is vital in understanding their oscillatory behaviour (Song et al. 2020) but will not be useful in analysing their collective dynamics. To mention, phenomena such as synchronization, chimera states, Spiral waves etc., cannot be studied with such low population neuron models. Hence in the recent years the research throng has shifted towards the neuronal network modelling for biological realism. The challenge starts from modelling the action potentials such as frequency, inter spike intervals (ISI), bursting or frequency of bursting. A map-based neuron is identified with effective in computationally feasible and tailored for investigating phenomenological behaviours (Izhikevich and Hoppensteadt 2004;Kinouchi and Tragtenberg 1996;Rulkov 2002;Shilnikov and Rulkov 2004). The dynamical map can be defined as dynamical systems with discrete space and discrete time, but with continuous state variables. In (Kinouchi and Tragtenberg 1996) an approach is presented with simplification of continuous time differential equation models into discrete time and keeping the desired properties. The early work on Map-based models (Karthikeyan et al. 2021;Rajagopal et al. 2021a, b, c) enumerated with advantages such as independency from integration and integration step adjustments. But the existence of piecewise linear function affects the performance particularly phase plane analysis (Kaneko 1993;Tanaka et al. 2006). Formulation of network topology for combining the discrete time elements with synapse model is the other challenge. Among various types, diffusive coupling models (Ibarz et al. 2011;Rulkov 2002) found its attraction because it uses synaptic current through a time dependent (or not) conductance. In (De Vries 2012;Rulkov 2002) Rulkov and De vries investigated the synchronization and existence of spatiotemporal chaos with neuron coupling strength. (Rulkov 2002;Shilnikov and Rulkov 2004;Tanaka et al. 2006) extended the model to study the chaotic phase synchronization and burst synchronization and proved that discrete maps are superior to study the network behaviour of neuronal system. The major challenge is, system stability not matching with the map stability and results with restriction to study the delay nature and temperature effects on the neuron synapse and network (McGahan and Keener 2020; Budzinski et al. 2019;Wang et al. 2015).
In (Braun et al. 1994) Huber Braun derived Temperature dependent Hodgkin-Huxley-type neurons and presented various dynamic nature of the neuronal system under different temperature conditions, which other models fails to show. The model ensembled with two slow, subthreshold currents along with the leak and channel currents and henceforth called Huber Braun (HB) neurons. The model can expose different spiking patterns and bursting and chaotic nature (Bazhenov et al 2005;Mikhail et al. 2007;Shilnikov and Rulkov 2004). The experimental evidences show that the transmission of signals between neurons in single layer is considerably different from neurons in different layers (Izhikevich 2000;Tan et al. 2020;Kamimura 2019;Wang et al. 2019), it extends the scope of study from single layer into multilayer network. The conductance is often changed when the ion channels are blocked. The ratio of working ion channels of potassium (and calcium) to the total ion channels of potassium (and calcium) can be used to partially blocking the ion channels of a fraction of neurons (local poisoned area). The unexcited or partially excited neurons can be identified as defects and presence of bunch of defected neurons is denoted as obstacle (Lenk et al. 2010). These obstacles potentially influence the propagation of travelling waves and contributes re-entrant of wavefronts and generation of spiral waves etc., (Starobin et al. 1996;Rostami et al. 2018a, b).
The advent of irregular patterns in electrically coupled networks are studied in many literatures (Kotini and Anninos 1997;Wei and Luo 2007;Ibarz et al. 2007;de Pontes et al. 2008;Kaneko 1993;Wang et al. 2008;Tanaka et al. 2006;Zou et al. 2009) but many not considered the presence of obstacles and very few discussed influences of the obstacle size and orientation (Rajagopal et al. 2019a, b;Majumder et al. 2018;Rostami et al. 2018a, b;Prado wt al. 2014;Wang et al. 2018;Wu et al. 2019;Feng et al. 2019;Yuan et al. 2017). Even those are not considered multilayer networking strategy.
A rotating spiral wave with its frequency higher than the frequency of a pacemaker often disturbs the regularity of heartbeat and thus causing arrhythmias and cardiac death. More than a spiral seed the danger imposed by multiple broken spiral waves are more because of its multiple frequency of rotation which increases the turbulent activity of the heart tissues and thus causing atrial fibrillation (Rostami et al. 2018a, b). This unsynchronised mechanical contrition of the heart muscles causes the irregular pumping of the blood. Such situations can often arise because of the heterogeneity and obstacles (Starobin et al. 1996) but this theory was not completely proven in the literatures as the relation between cardiac arrhythmias and reduced excitability is still a topic of debate. Using a discrete neuron model such complexities are resolved as the resetting happens in every step without loss of the actual behaviour of the neuron.
Motivated by this in the present study we propose a discrete HB neuron model and investigating the behavioural patterns of the discrete neuron such as periodic and chaotic bursting. Further to study the network behaviour we used a 2D lattice of the discrete HD neurons.
The manuscript is organized as follows Sect. 2. Describes the model of Discrete Huber-Braun neurons with their parameter values. Section 3. Discussing about the dynamical analysis of the model including bifurcation diagram, Lyapunov spectrum and complexity measure Sect. 4. Deals with formulation of network model and the effects of obstacle on the spatiotemporal behaviour. Finally the conclusion section is provided.

Discrete Huber-Braun neuron (DHB) model
The continuous time model of a Huber-Braun type neuron which is the modification of the well-known Hodgkin-Huxley models derived by including two slow, subthreshold currents along with the leak and channel currents was discussed in (Finke et al. 2008). But many literatures have shown that discrete neuron maps are valid phenomenological neuron models (.Braun et al.). Normally discretisation can be applied with respect to space, time and field in neuron models. We are interested in applying discretisation in time and propose a new discrete neuron model derived from the continuous time Huber-Braun type neurons. The new neuron map can be called as Discrete Huber-Braun neuron (DHB) model whose mathematical model is defined by.
Vðk þ 1Þ ¼ ÀI l ðkÞ À I d ðkÞ À I r ðkÞ À I sd ðkÞ À I sr ðkÞ C M a r ðk þ 1Þ ¼ /ðTÞ ða r1 ðkÞ À a r ðkÞÞ s l a sd ðk þ 1Þ ¼ /ðTÞ ða sd1 ðkÞ À a sd ðkÞÞ s l a sr ðk þ 1Þ ¼ /ðTÞ ðÀgI sd ðkÞ À ka sr ðkÞÞ s sr ð1Þ where the currents are defined by the relations with I l defining the leak current, I d and I r representing depolarising and repolarising currents, I sd and I sr representing the sub-threshold depolarising and repolarising currents respectively.
In all the current and function definitions the suffix d denotes depolarising and r denotes repolarising. The slow depolarising and repolarising parameters and functions are defined with a suffix sd and sr respectively. The steady state activation functions are defined by the relations, In order to define the temperature dependencies for the neuron model (1), we use the Q 10 laws (Finke et al. 2008) defined by the mathematical relations.
The system parameters are adopted from (Finke et al. 2008) and defined as.
Activation time constant : s r ¼ 2 ms; s sd ¼ 10 ms; Set the initial conditions

Dynamical properties
In this section, we focus on the complex dynamics of the system with system parameters g l , g d and T, where bifurcation diagram, Lyapunov exponents (LEs) and sample entropy algorithm are employed. It should be noted that the Jacobian matrix of the neural model is obtained used the matlab function ''J = jacobian(.)'', and the QR decomposition method is applied to estimate Lyapunov exponents.
In the figures, since LE3 and LE4 are both smaller than zero, we just illustrate the Lyapunov exponents LE1 and LE2 for simplification. The initial conditions are given as

Bifurcation and Lyapunov exponents
Let g l =0.1 and g d =1.5, vary the parameter T from 8°C to 28°C with step size of 0.0401. The obtained bifurcation diagram and Lyapunov exponents are shown in Fig. 1. It shows that there is a large periodic window when T [ [12.14, 21.1]. In the left and right side of this periodic window, the system is chaotic. However, the system is nonchaotic when T \ 8.9°C. It also shows that there are many small periodic windows with the variations of T in the chaotic regions. Secondly, the Lyapunov exponents-based chaos diagrams are obtained in different parameter planes. Let g l =0.1, T varies from 8 to 28 with step size of 0.1 and g d varies from 1 to 2 with step size of 0.005. The analysis result in the parameter plane T-g d is shown in Fig. 2a. Let g d =1.5, T varies from 8°C to 28°C with step size of 0.1 and g l varies from 0 to 1.5 with step size of 0.0075. The analysis result in the parameter plane T À g l is shown in Fig. 2b. Let T = 11, g d varies from 1 to 2 with step size of 0.005 and g l varies from 0 to 1.5 with step size of 0.0075. The analysis result in the parameter plane T À g d is shown in Fig. 2c. In Fig. 2, the threshold value for chaos is set as 0.01. Thus we just show the regions with largest Lyapunov exponents than 0.01 and keep the rest states including periodic and divergent states with white color. As shown in those chaos diagrams, the system has wide regions for chaos. Meanwhile, there are also regions with Lyapunov exponents smaller than zero, and the white regions mean that the system is divergent. It shows that the system has rich dynamics with the variation of parameters g l , g d and T. (Bifurcation and Lyapunov exponents for the parameters g l & g d provided in the supplementary material).

Complexity measure
In this paper, complexity measure of time series generated by the neural is carried out based on the sample entropy (SampEn) algorithm.
The physical significance of SampEn is not difficulty to find out. It measures the probability of the generation of new patterns in the time series. The larger the measure value is, the greater the probability of the generation of new patterns. At present complexity analysis plays an important role for dynamics of different chaotic systems (Ma et al. 2021;Natiq et al. 2019;Politi 2017). Complexity of the system with the variation of different parameters are shown in Fig. 3. The results shows high complexity with chaotic state while show low complexity and zero complexity for the other states. Obviously, for those periodic windows, there shows low complexity measure results. When g d takes values between 1.4 and 1.63, the complexity increases with the increase of g d , and when T [ 21.1°C, the complexity of the system increases with the increase of T. However, compared with dynamic analysis results, SampEn can build a clear version regarding the complexity of the system. Thus, we analyzed complexity of the neural model in the parameter planes.
We use the same parameter settings and initial conditions to simulate the system for complexity measure. The results are shown in Fig. 4 where we draw the SampEn based chaos diagrams using those measured values which are larger than 0.1. According to the above analysis results, it shows that the neural system can also have high

Spiral waves in the DHBN model
The major challenge of a map-based neuron model is to correlate the behaviour with a biophysical neuron characteristic. Replicating the network behaviour of a continuous time neuron model with a map-based neuron is considered complex as phenomenon like the spiral waves and multiarmed spiral waves are not very well documented by literatures (Huang et al. 2010). Hence a discrete neuron model capable of replicating the behaviour of an ODE neuron model in both local and collective sense is very rarely discussed. The DHBN model proposed in this work has effectively showed the characteristics of a biophysical neuron (Braun et al. 2000) and we now are interested in investigating the network behaviour of the DHBN model.
For this we arrange the neuron in a network and depending on the number of layers in the network we subdivide the discussion in two cases. In the first case we investigate using a single lattice array of DHBN model and in the second case we consider two-layer network with each layer being the lattice array of the DHBN model. Further we expand the investigation by considering two separate scenarios which could affect the excitability of the neurons and thus could contribute to re-entry and spiral waves. In the first scenario we considered different sizes and orientations of obstacles and show their impact on the wave propagation and re-entry. The mathematical model of the DHBN network can be defined as.
The stimuli are applied to the network with the condition h is ; h jq ¼ 1 when i ¼ s; j ¼ q . The stimulus is applied from the left side of the network by setting s ¼ 100; q ¼ 1 . The stimuli settings are A ¼ 0:1; x ¼ 0:01 and we capture the spatiotemporal behaviour of the network at the end of 6000 s.
When the stimuli are applied to the network, the tissues are highly excited and we could observe simple travelling waves from left to right boundary and no signs of wave reentry. Varying the stimulus amplitude, frequency and the system parameters doesn't change the excitability of the neurons. For certain amplitude and frequency values of the stimulus the entire network goes unstable. Hence, we are interested in investigations different scenarios which will influence the wave re-entry and spiral waves in the network.
In our previous work (Rostami et al. 2018a, b), we have showed that obstacles in the network can initiate wave reentry and can contribute to the spiral wave formations. Here we consider that the network has some defects defined by the obstacles on the path of the wave propagation. Effects of obstacle on the spatiotemporal behaviour During absence of obstacles the network tissues are in high excitability and hence plane waves propagate in the media. Now we insert an obstacle in the network and investigate the effect of the size and orientation of the obstacle on the wave front. We consider the system parameters as in (5) and the stimulus setting are same. For the orientations of the obstacle, we have considered parallel and perpendicular obstacles with reference to the wave front. Different sizes of obstacles having heights H 2 60; 20 ½ and width W 2 20; 10; 5; 1 ½ with the measuring unit being the nodes. First, we considered parallel obstacles with different combinations of H and W as shown in Fig. 5. When the obstacle width is less than 5 nodes, there is no wave re-entry in the media and we could note just a phase shift in the wave propagation as in Fig. 5a for the [H, W] = {60, 1}. When the width is increased to 5 and keeping the height as 60, the obstacle creates wave re-entry and spiral wave seeds in the right boundary of the obstacle. Further increasing the width to 10 two spiral pools are created at the top and bottom ends of the obstacle and we could also note that the spiral waves are formed inside the obstacle too. Broken multiple spiral waves can be generated by further increasing the width to 20 as in Fig. 5d. Now we decrease the height to 1 and width to 5 and no phenomenal changes in the excitability of the neurons and hence travelling plane waves are noted. Now an increase in just the width of the obstacle to 10 making it a square block could induce wave re-entry and spiral waves shown in Fig. 5f.
Considering the obstacles placed perpendicular to the wave front, we considered different sizes of obstacles with width W 2 f20; 10; 5g and length L 2 f60; 20; 10; 5g. Firstly, we consider an obstacle of 10 9 5 and the excitation is applied to the network. The snapshots shown in Fig. 6a confirms that the obstacle breaks the wave at the interaction points while no re-entry or spiral wave phenomenon observed. By increasing the width to 20 and length kept at 10, the spiral wave seeds are seen towards the upper and lower part of the obstacle as in Fig. 6c. For an obstacle of width 20 and length 60, symmetric spiral waves originate at the right side of the obstacle and rotates with a larger arm radius. But comparing Fig. 6 with 6, we could comment that the obstacles parallel to the wave front are more effective in the excitability and wave re-entry compared to the perpendicular obstacles. Now the question arises of what will be effect of obstacles in multilayer networks. Will the obstacle in a layer influence the excitability in the other layers too? If so under what circumstances such excitability changes happens. To answer all these questions, we construct the mathematical model of the two-layer DHBN network as in (7). are used for the analysis with obstacles greater than five node width we could note that wave re-entry and spiral waves are formed near the boundary opposite to the wave entry side L;ij ¼ ÀI k L L;ij À I k D L;ij À I k R L;ij À I k SD L;ij À I k The inter layer coupling term is defined by D LU ðV L;ij À V U;ij Þ for the upper layer and D 21 ðV U;ij À V L;ij Þ for the lower layer. The intra layer coupling is like the one used in single layer network (7) with r U ; r L defining the respective coupling strengths. We have considered stimulus whose amplitudes are A U and A L for upper and lower layers respectively and by setting the amplitude of the stimulus in the respective layer to '0' we could show the effect of coupling effects between layers.
In the first discussion we consider the obstacles to be parallel to the wave front and stimulus is applied only to the upper layer by setting A L ¼ 0 . We have considered similar sizes of obstacles as considered for the single layer discussion and the obstacle is inserted only in the upper layer. Though there is no obstacle in the lower layer, the excitability of tissues is controlled by the first layer through the coupling as there is no stimulus in the lower layer. Hence, we could note that both the layers operate in synchrony and the obstacle effects seen in the upper layer is transferred to the lower layer and thus both the layers showing similar re-entry and spiral wave phenomenon. These observations can be easily seen in the snapshots shown in Fig. 7. There is an identifiable difference between the upper layer and lower layer snapshots with the lower layer displaying a much lesser difference between the media and excited tissues while the upper layer shows a definitive larger difference. This is because the excited state of neurons in the lower layer is of lesser amplitude as the coupling though is capable of achieving coherent behaviour between the layers, they are not effective in coupling the exact excitability conditions to the lower layer.
When stimulus is applied to only upper layer the obstacle effects are coupled to the lower layer even though there is no physical obstacle present. Now we consider stimulus in both the layers and use the same orientation and size of obstacle as in the previous case. Now we could show that though the coupling strength or the stimulus amplitude is note changed as from Fig. 7, the obstacle effects in the lower layer are not displayed. In other words, now the lower layer tissues are excited by the applied stimulus rather than the upper layer effects as shown in Fig. 8. Thus, the effects of coupling bond between the layers are now eliminated by applying stimulus to the respective layers. Also, we could note that the difference between the excited amplitude and the media in the lower layer are significantly higher and thus overcoming the blur effect seen in Fig. 7. Though not presented here, the same phenomenon can be observed when we consider the lower layer obstacle and excitation with the upper layer not exposed to obstacles or stimulus.

Conclusion
A discrete neuron map derived by modifying the continuous time Huber-Braun neuron model is proposed in this paper. The novelty of the model is the inclusion of two slow, subthreshold currents along with the leak and channel currents. We have considered the temperature effect on various current channels of the neurons and depending on the temperature the neuron shows periodic spiking and chaotic bursting behaviour. The bifurcation studies reveal the different periodic and chaotic regime of the DHBN model. To clearly show the chaotic and periodic regime in a 2D parametric space, we have used the Lyapunov exponents-based chaos diagram. A sample entropy algorithm is proposed to measure the complexity of the DHBN model which is correlated with the Lyapunov exponent chaos diagram. Though nodal dynamics will make us understand the complete behavioural pattern of the neuron, it is the collective dynamics which plays a significant role in biophysical models. Hence, we constructed a simple 2D network of DHBN modes and investigated the network behaviour by applying a periodic force. The network shows simple plane wave propagation with no sign of wave reentry on spiral waves. Now our interest is to investigate the  60,20] conditions which can induce wave re-entry. We introduced an obstacle in the network and decided to investigate the effects. For this we have chosen two orientations and six different sizes of obstacle. The observation shows that wave re-entry and spiral waves are formed near the boundary opposite to the wave entry side. Also, we have shown the minimum dimensions of the obstacles which can cause wave re-entry and induce spiral waves. To check the obstacle effects on multilayer networks, we constructed a two-layer DHBN network. Firstly, we showed that though the obstacle is considered only in one layer, the stimulus applied to the layers play an important role is the wave reentry in the layers. This is proved by setting only stimulus and obstacles to the first layer, but the second layer also displays stimulus wave re-entry like the first layer. This coherent behaviour is removed when we apply stimulus to the second layer also.