Neuronal Dynamics: from Complexity to Simplicity

ABSTRACT


A. INTRODUCTION
The neuroscientists employ the mathematical models mainly to understand the brain activities from a cellular to the network level. Various models have been introduced to explain the neuron dynamics namely the Hodgkin-Huxley model (Hodgkin & Huxley, 1952), FitzHugh-Nagumo systems (FitzHugh, 1961), including Leaky Integrate-and-Fire (LIF) neuron (Gerstner, Kistler, Naud, & Paninski, 2014).
The four-dimensional nonlinear Hodgkin-Huxley model describes quantitatively the initiation and propagation of action potential. It helps to understand the brain functionalities such as ion channels modulation, the circumstances which control both rate and timing of action potential, and the neural coding and information processing (Catterall, Raman, Robinson, Sejnowski, & Paulsen, 2012). The changes on ion channel function like the opening and closing of sodium and potassium gating variables are the key concepts introduced by Hodgkin and Huxley in their model. The ion channels response may be affected through a neuromodulation in the axon while the gates are very sensitive to the incoming stimulus even if it is weak (Burke & Bender, 2019) (Mathie, Kennard, & Veale, 2003). As a consequence, the neuron shows a periodic firing activity when the stimulus exceeds a critical value.
Fitzhugh and Nagumo proposed a two-dimensional nonlinear system to explain the basic excitability properties as exhibited by the Hodgkin-Huxley model (Izhikevich & FitzHugh, FitzHugh-Nagumo model, 2006) (Sherwood, 2013). In principle, the Fitzhugh-Nagumo model is the extension of the Van Der Pol oscillator (see equations 2-3 in Ref. (Kanamaru, 2007)): they introduced the variable input current in the Van Der Pol fast equation to mimics the experimental injection of a current into the cell membrane. The cubic and line nullclines intersect at one equilibrium point and its stability is determined by a constant current. If no current is applied, the equilibrium is stable i.e., the neuron remains in the resting membrane potential and no action potential is generated. If the current is weak, the neuron is in the subthreshold regime and again no action potential is generated. Here, the trajectories decay quickly to the resting potential from a particular initial condition. On the other hand, if the current is sufficiently strong the equilibrium becomes unstable and the limit-cycle oscillation is born (Hopf-bifurcation) i.e., the neuron fires a spike periodically.
In contrast to the Hodgkin-Huxley and Fitzhugh-Nagumo models, the LIF neuron neglects many physiological details of the neuron such as the activation of ion channels, the shape of the action potential etc. The model is one-dimensional linear differential equation describing the evolution of membrane potential accompanied by the reset constraint (Gerstner, Kistler, Naud, & Paninski, 2014). If the membrane potential passing the threshold, it suddenly jumps to the reset value. One can expect a regular spiking behavior in the LIF neuron as those displayed by the two other models, if the equilibrium point is larger than the threshold.
It is relatively easy to simulate a single neuron, yet the task is not always trivial since it relays on the assumptions being considered. On the other hand, one must deal with the computational costs when simulating a large-scale network (Brette, et al., 2007). Fortunately, the emergence of periodic oscillation in the real neuron paves a way to treat the state variables, e.g., membrane potential as a phase. We may eliminate the electro-chemical aspects of the neuron and then limit ourselves only to the position of membrane potential along the periodic orbit, or equivalently, the phase on the circle (Izhikevich & Ermentrout, 2008). The resulting model of pulse-coupled phase oscillators is relatively simple and computationally efficient (Politi & Rosenblum, 2015). Despite of its simplicity, the model is indeed able to produce the rich dynamical regimes, ranging from synchronization transition (Strogatz, 2000), self-partial synchronization (Barabash, Belykh, Osipov, & Belykh, 2021), including irregular collective behaviour (Afifurrahman, Ullner, & Politi, 2021).
The phase reduction is a commonly used method to reduce the multidimensional dynamical systems into one-dimensional phase equation (Nakao, 2016). In this study, we discuss the phase reduction process to acquire the phase representation for LIF neurons under the following conditions: (1) the neurons are all identical and show periodic firing; and (2) weak synaptic interactions. The first condition depends on the appropriate parameters selected for the neuron's model being studied, while the second one ensures a small deviation of the neuron from its limit-cycle (Winfree, 1967). The LIF model is selected as it has been repeatedly investigated in the literatures and widely used both for analytical and numerical studies (Abbott, 1999) (Brunel & Van Rossum, 2007) (Coombes, Thul, & Wedgwood, 2012) (Godinez, Sossa, & Montero, 2017). Our research systematically verifies the equivalence of phase oscillator and integrate-and-fire models studied in (Politi & Rosenblum, 2015).

B. METHODS
The research procedures consist of three steps: (1) models review; (2) phase reduction; and (3) simulation. More specifically, we review the mathematical model for the LIF neuron and discuss its behaviour in the presence of constant input. We then discuss the phase model and derive the corresponding phase representation for the LIF neuron. Next, we explain the effect of pulses to the neuron's dynamics and transform the LIF neuron into the phase model by assuming a weak pulse input. Accordingly, we extend the phase reduction approach to obtain a simplified version of LIF neuronal networks, namely pulse-coupled phase oscillators. At the end, we verify the existence of collective behaviour (e.g., synchronization) through numerical studies. We employ the Euler forward integration scheme to simulate our model of pulsecoupled phase oscillators and use a statistical measure to test the emergence of synchronization. The research procedures are summarized in Figure 1.  The LIF model is defined by a linear differential equation describing the membrane potential u(t) when it is stimulated by an input current. If u(t) reaches a threshold , the neuron is said to fire a spike at time ts and u(t) is instantaneously reset to ur at t = ts+.
combined with a reset condition for = + where is the membrane time constant and defined as = .
In the absence of any input i.e. I(t)=0, the membrane potential u(t) decays to the resting membrane potential as → ∞. When the passive membrane potential is stimulated by a constant current I(t) = I0, its trajectory is obtained by integrating Eq. (1) from the initial value u(0) = ur. The solution for t > 0 is given by From Eq.
(3) we can clearly see that ( ) = 0 when → ∞. The neuron emits a spike, if the equilibrium point RI0 is larger than , otherwise there is no spike emitted. The time needed by the neuron to reach the threshold can be found by solving where ts is the period of oscillation. In case the neuron has a refractory time tr, i.e., a time interval when the membrane potential stays constant and equal to the reset value after the spike has fired, the Eq. (5) becomes The firing rate, , of single neuron is given by 314 | JTAM (Jurnal Teori dan Aplikasi Matematika) | Vol. 7, No. 2, April 2023, pp. 310-323 = 1 .
b. Phase Oscillator Model According to Eq. (4) the neuron fires a spike periodically if the membrane potential exceeds a threshold (see also Figure 2). Here, we discuss a simplest mathematical model to represent the periodic oscillation. It approximates the behaviour of nonlinear systems with a multidimensional state variable by using only a single-phase variable. Suppose that is nonlinear dynamical system having a stable periodic orbit ( ) with a period T0 and x0 is the initial condition. The only relevant parameter needed to characterise the oscillator is the position along the closed curve. This is naturally a phase-like variable and can be measured as the amount of time elapsed from the crossing of a given threshold.
The phase is normally defined on the interval [0,2 ). Here, we rescale the time in such way that every point on the oscillation can be uniquely described by the phase in the unit interval, i.e., Φ ∈ [0,1). More precisely, the phase is defined as the time multiplied by the frequency (Nakao, 2016), i.e.
Taking a derivative of Eq. (9) with respect to time gives which defines the phase model for the system (8) Thus, by definition of the phase in Eq. (9), we have where T0 is defined in Eq. (6). Figure 2 illustrates the time series of membrane potential u(t) and the corresponding phase Φ( ) for LIF model driven by a constant current. In this simulation we choose a set of parameters: = 1Ω, 0 = 20 , = 0 , = 15 , = 10 , = 0 . As seen, the neuron fires a spike in regular manner with a constant period. The dashed line represents the spike trains. The main point is that there one can see du/dt is not constant, while dΦ( )/dt it is, as shown in Figure 2.

Pulse Input
Let's consider two identical LIF neurons where 1 ( ) and 2 ( ) are used to describe the membrane potential for neuron 1 and 2, respectively. The evolution equation of membrane potentials for the two neurons is given by and combined with the reset conditions 1 ≥ ⟹ 1 = 2 ≥ ⟹ 2 = .
The variables 1 ( ) and 2 ( ) characterize the synaptic inputs for neuron 1 and 2, respectively and defined as 2 ( ) = ∑ ( − 1 ) + 0 The synapse is excitatory (inhibitory) if > 0( < 0). The sums run over nth pulse emitted by the neuron 1 and 2 at a time 1 and 2 , respectively. Each pulse being emitted by a single neuron is described by the function s(t). The shape of a single pulse can take many forms including delta , exponential, and alpha pulses (Gerstner, Kistler, Naud, & Paninski, 2014) (Afifurrahman, Ullner, & Politi, Collective dynamics in the presence of finite-width pulses, 2021), as shown in Figure 3. with the inhibitory -pulse input. Figure 3 depicts the time trace of membrane potential for the two LIF neurons interacting through delta pulses. The initial conditions for the membrane potential 1,2 (0) are selected from uniform and random distributions. The same parameter set is chosen as in isolated neuron driven by a constant current I0 (see Figure 2).
A general remark concerns on the oscillatory behaviour of the neuron 1 and 2 under the action of a constant current I0. When = 0, the membrane potential for neuron 1 and 2 evolves independently and have equal free-running period T. We notice for the mutually inhibitory synapse; the firing of the presynaptic neuron kicks the membrane potential of the postsynaptic neuron backward by a fixed amount leading to a delay on the firing of the postsynaptic neuron.
For mutually excitatory synapse, the firing of the presynaptic neuron kicks the membrane potential of the postsynaptic neuron forward which accelerates the firing rate of the postsynaptic neuron. Interestingly, after a finite amount of time the phases of two neuron locked and they start to fire simultaneously as a result of such mutual interaction. This issue has been initially studied rigorously by Peskin (1975), and later on generalized by Mirollo and Strogatz (Mirollo & Strogatz, 1990) for any two pulse-coupled oscillators with delta pulses.

Phase Description for LIF Neuron: Weak Pulse Input
According to the Figure 3, the neuron reacts to the incoming delta pulse by shifting its trajectory in the opposite direction of the motion due to inhibitory synapse. If the synapse is excitatory, the trajectory is shifted in the same direction as the motion, causing an acceleration on the firing rate. The response of the oscillator to the external stimulus (e.g., pulse) is encoded in a phase-dependent function, called phase response curve (Canavier, 2006).
First of all, let's recall the phase representation of a stable periodic orbit Φ⁄ = where is the frequency of oscillation. If the oscillator weakly interacts with its environment (weak interaction means that the periodic orbit is not substantially affected by the external force), then the phase model can be written as (Pikovsky, Rosenblum, & Kurths, 2001) where is the phase of forcing, Q is a periodic function of both arguments, and is the coupling strength. In many cases, Q can be expressed as where Γ(Φ) is the PRC and ( ) is the forcing function. Let's give a look again at the LIF neuron. Assume that the membrane potential evolves according to with the reset condition and constant current I0. A delta pulse arrives at time tp and it kicks the membrane potential u by a fixed amount .
If is sufficiently small (typically ≪ 1), the Eq. (18) can be reduced to the phase model as follows The expression inside the square bracket (R.H.S) determines the PRC. In order to obtain the phase-dependent function, it is necessary to substitute u(t) given in Eq. (3) into Γ( ( )) = (− ( ) + 0 ) ⁄ and these yields 318 | JTAM (Jurnal Teori dan Aplikasi Matematika) | Vol. 7, No. 2, April 2023, pp. 310-323 while T0 is defined in Eq. (6). The internal complexity of LIF neuron is now encoded in the function Γ(Φ). Notice that there are alternative ways to derive the PRC for the LIF model (see (Abbott & Vreeswijk, Asynchronous states in networks of pulse-coupled oscillators, 1993) (Lewis & Rinzel, 2003) for further references).

Networks of Pulse-Coupled Phase Oscillators
Now we extend the phase reduction approach to N system sizes. First, consider a population of N network of identical LIF neurons. The evolution of the membrane potential of the jth neuron (j = 1, 2, …, N) is given by with the reset condition ≥ ⟹ = .
Analogous to the case of two neurons, the variable ( ) can be defined as where I0 is the input current. The connectivity from the presynaptic neuron k to the postsynaptic neuron j is defined by a matrix G where Gj,k = 1 if → is active, and Gj,k = 0 otherwise. The first sum on the R.H.S is performed over all presynaptic neurons (k = 1, 2, …, N), while the second sum runs over nth pulse emitted by the kth neuron at times . The scaling factor 1/N is important to keep the expected input finite when N goes to infinity. In principle, the transformation of Eq. (20) to the phase oscillators can be done by following the same procedure as described in the section 1.3. Doing so, we obtain where Γ(Φ ) has already defined in Eq. (19). We refer to the model as pulse-coupled phase oscillators. In comparison to the original LIF model in Eq. (20), each neuron is now characterized by the phase Φ ≤ 1. If Φ passed a threshold Φ ℎ = 1, it is reset to the value Φ = 0 and enters a refractory time tr. At the same time, the nth pulse is emitted and sent to the neuron j at time .

Simulations a. Setup
Our object of study is N pulse-coupled phase oscillators defined in Eq. (22). We perform the Euler algorithm in C programming languages with the integration time step ∆ = 10 −3 . The table below highlights all parameters being used for the simulation, as shown in Table 1.  Mirollo & Strogatz, 1990) > 0 (Mirollo & Strogatz, 1990) , = 1 for all ≠ (Mirollo & Strogatz, 1990) , = 1 for all ≠ (Mirollo & Strogatz, 1990) ( ) = ( ) pulse (Mirollo & Strogatz, 1990) ( ) = ( ) pulse (Mirollo & Strogatz, 1990) Notice that, the phase-response curve Γ(Φ ) is the same for all neurons and its value depends on the parameters imposed in the original LIF neuron. Furthermore, we assume that the neurons interact to its neighbourhoods through pulses. b. Statistical Measure The degree of synchrony in a large population of neurons is measured based on the fluctuations of the global variable (e.g., the average phase). The formula is given by (Golomb D. , 2007).
where ⟨∘⟩ represents an ensemble average, while the over-bar is a time average. The numerator is the variance of the ensemble average while the denominator is the ensemble mean of the single-neuron's variance. If the neurons in the population fire in asynchronized manner, then = 0. However, if all those neurons fire in complete synchrony, exactly at the same time, then = 1. If the neurons are firing synchronously but not in a complete synchronization, then the neurons are said to be partially synchronized. Here, the order parameter is strictly larger than 0 and smaller than 1. c. Numerical Analysis We simulate a network of N = 100 phase oscillators with all-to-all connectivity, i.e., , = 1 for every ≠ while the coupling strength is set to = 0.2. The initial conditions are selected randomly from a uniform distribution within the unit interval Φ (0) ∈ [0,1]. The first two panels in Figure 4 show the time series of Φ 1 and Φ 3 out of 100 oscillators, respectively. It is clear that the two oscillators fire a spike exactly at the same time and follow periodic dynamics with a period approximately ≈ 13.9 . The value fits the formula given in the Eq (6) for LIF neuron. On top of that, the ensemble average 〈Φ〉( ) in the lower panel shows the same scenario as the two other cases, meaning that all oscillators are identical, as shown in Figure 4. In order to be more quantitative, we compute ( ) for a longer time span (see Figure 5). There one can see the degree of synchronization settles around ( → ∞) ≈ 1, which means that all neurons are completely synchronized. The finding is in line with the synchronous behaviour studied originally in a single population of LIF neurons by (Mirollo & Strogatz, 1990). A detail stability analysis of synchronous states in the generic model of phase oscillators is investigated by (Afifurrahman, Ullner, & Politi, Stability of synchronous states in sparse neuronal networks, 2020). In this study, however, we limit ourselves on testing the emergence of synchronization to confirm the validity of the model being studied. In addition, the pulse-coupled phase oscillators (22) can be examined by imposing a scaling → ∞. In fact, a true synchronous state is independent of the system size N (Golomb, Hansel, & Mato, 2001). The black and red curves in Figure 5 correspond to N = 100 and N = 200, consecutively. As seen, both synchrony measures are relatively close to 1, which indicates that the synchronous state is preserved, as shown I Figure 5. Interestingly, the degree of synchronization ( ) converges to 1 as N grows (the red curve is on the top of black one) suggesting the existence of a large number but finite N at which = 1. This observation is consistent with the law of large numbers where the precision is proportional to the system sizes (Weisstein, 2023).

D. CONCLUSION AND SUGGESTIONS
In this paper, we have described the implementation of phase reduction principles to obtain the phase oscillators model for LIF neurons. First, we consider simple case where a single neuron is injected by a constant current. We construct the variable time as a function of membrane potential, i.e., ( ) and use the definition of the phase to obtain Φ( ). In the second case, we consider the dynamics of two identical neurons interacting through pulses. While assuming the two neurons are weakly coupled, we derive the exponential shape of phase response curve for LIF neuron. Next, we extend these formulations for any N networks of LIF neurons and derive the corresponding phase oscillators model.
In the last step, we perform the numerical simulation to validate the phase oscillators model of Eq. (22) and compare the results with the original LIF neurons by confirming the existence and persistence of synchronous periodic dynamics. We conclude that both models are indeed equivalent and recommend the phase oscillators as an alternative model to investigate a synchronization in the populations of neurons. However, one must take into account both time and computational costs when doing the simulation for a large system size. Finally, the question of how fit the model with the experimental data is worth to explore as this allows us for better understanding of the brain activity.