Numerical Approximation for Nonlinear Noisy Leaky Integrate-and-Fire Neuronal Model

Dipty Sharma 1,*, Paramjeet Singh 1, Ravi P. Agarwal 2 and Mehmet Emir Koksal 3 1 School of Mathematics, Thapar Institute of Engineering & Technology, Patiala 147004, India; paramjeet.singh@thapar.edu 2 Department of Mathematics, Texas A& M University-Kingsville, Kingsville, TX 78363, USA; agarwal@tamuk.edu 3 Department of Mathematics, Ondokuz Mayis University, Atakum, Samsun 55139, Turkey; mekoksal@omu.edu.tr * Correspondence: dipty.sharma@thapar.edu or diptysharma16@gmail.com


Introduction
The large-scale neural network models in computational neuroscience have become familiar.The classical description of these (excitatory-inhibitory) neural network models is based on the deterministic/stochastic system.One of the most common models is known as the noisy leaky integrate-and-fire (NLIF) neuron model in which the behavior of the whole population of neurons is encoded in a stochastic differential equation (SDE) for the time evolution of membrane potential of a single neuron representative of the network.The dynamics of a single neuron is given by ( [1][2][3][4][5][6]) where V(t) represents the membrane potential of a single neuron and τ m is the time relaxation of the membrane potential in the absence of any communication.The communication of a single neuron with the network is modeled by the synaptic input current, I(t).The form of I(t) is a stochastic process, given by ( [7]) Here, each spike is treated as a delta function, and if a spike occurs at time t = t 0 , it is denoted by δ(t − t 0 ).The terms t n E m and t n I m in above equation represent the time of mth-spike receiving from nth-presynaptic neuron for excitatory and inhibitory neurons, respectively.Moreover, the terms N E and N I are the total number of presynaptic neurons, where J E and J I are the strength of the synapses for excitatory and inhibitory neurons, respectively.Since the above form of synaptic input current is the discrete Poisson process, it becomes very difficult for further investigation.In addition, the researchers have used the diffusion approximation in which the synaptic input current I(t) is approximated by a continuous in time Ornstein-Uhlenbeck-type stochastic process as given by Initially, it is assumed that every neuron generates spikes according to a stationary Poisson process with constant probability of generating a spike per unit time r, and it is also assumed that all these processes are independent between neurons, because of these assumptions, the mean value of the current, indicated by µ c , is given by br = (N E J E − N I J I )r and its variance, σ 2 c = (N E J 2 E + N I J 2 I )r.Here, we have to depict the likelihood of firing per unit time of the Poissonian spike train r and it is thus recognized as the firing rate, which should be computed as r = I ext + N(t), where N(t) is the mean firing rate of the network.On the other side, B t is the standard Brownian motion in above equation.
The next important factor in the modeling is that the neurons generate a spike only when its membrane potential V(t) arrives at a certain voltage, known as threshold V F , and instantly reset toward a resetting potential V R < V F and sends a signal over the network.
Comprising the continuous form of I(t) in SDE model (1), we obtain where τ m = 1.We suppose that the voltage of a neuron arrives at threshold level at time t − o , i.e., V(t − o ) = V F and after that the voltage arrives suddenly at resting potential, i.e., V(t Furthermore, one can write the associated FPE with source term by using Ito's rule [8], for the evolution of probability density function p(v, t) ≥ 0 of finding neurons at a voltage v ∈ (−∞, V F ], with time (t ≥ 0) where 2 , and N(t) is the mean firing rate of the network which is computed as the flux of neuron at the firing voltage.The source term of the Equation (2) comes from the fact that when the neurons generate spikes and send the signals over the network, their voltage immediately reset to the reset potential V R .At the relaxation time, no neuron have the firing voltage, for this reason, the initial and boundary conditions are given by We can easily verify the conservation of the total number of neurons in above equation.For this purpose, we need to describe the mean firing rate for the network.Since the mean firing rate is the flux of neurons at V F , the value of N(t) is given by N(t) = −a(N(t)) ∂p(V F ,t) ∂v .Integrating (3) across the voltage domain and using the above boundary conditions, we obtain the required condition.
Equation (2) represents the evolution of probability density function and therefore If we translate the new voltage variable by considering In [8], the authors provide theoretical and numerical analysis using finite difference approximation.For the references of mathematical aspects of nonlinear NLIF models, we refer [9][10][11][12][13].We are concern about finding the value of the unknown p(v, t) using alternative approach.The problem (2)-( 4) cannot be solved analytically, because of its complexity arising from nonlinearity and having a source term [14][15][16].Therefore, numerical methods are generally used, for example, finite difference method (FDM) is used to find the approximate solution of the governing equation [8].However, FDM has some disadvantages, for instance, the singularity in the delta source term, as in the above mentioned equation, makes the solution divergent.In order to use the FDM appropriately, the governing equation must be modified, due to the fact that the procedure becomes complicated.
Hence, in the present work, we propose a formulation based on finite element approximation to find the solution of governing equation.FEM is one of the powerful numerical methods for the solutions of problems that describe the real-life situations.Moreover, the characteristics of FEM tackle the singularity problem in an effective manner [17][18][19][20][21][22][23][24][25][26][27].The applicability of FEM regarding this model problem is demonstrated in the final section.
Description of the present paper is as follows: We consider the NLIF model described by Equations ( 2) and (3).In Section 2, we develop the numerical approximation based on the finite element approach.The analysis for the stability is provided in Section 3. We report some numerical examples from [8] and discuss the solution behavior graphically in Section 4. In the last section, we conclude the work done in this research article.

Preliminaries
Here, we state some basic definitions and auxiliary results, which will be used throughout the manuscript.As we are studying a nonlinear version of the FPE, we start with the notion of weak solution.Definition 1.We say that a pair of non-negative functions (p, N) with p ∈ L ∞ (R + ; ) is a weak solution of (2) and (3) if for any test function Here, the space L p (Ω), 1 ≤ p < ∞, refers to the space of functions such that f p is integrable in Ω, while L ∞ corresponds to the space of bounded functions in Ω.The set of infinitely differentiable functions in Ω is denoted by C ∞ (Ω) used as test functions in the notion of weak solution.The blow-up of solution and a priori estimates are given in [8].We here just state the results.

Theorem 1. (Blow-up) Assume that the drift and diffusion coefficients satisfy
for all −∞ < v ≤ V F and all N ≥ 0, and let us consider the average-excitatory network where b > 0.
is close enough to e µV F , then there are no global-in-time weak solutions to ( 2)-( 4).
Lemma 1 (A priori estimates).Assume h(v, N) = −v + bN, a(N) = a 0 + a 1 N on the drift and diffusion coefficients and that (p, N) is a global-in-time solution of ( 2)-( 4) in the sense of Definition 1 fast decaying at −∞, then the following a priori estimates hold for all T > 0: In a latest work [8], it was demonstrated that the problem ( 2)-( 4) can produce a finite time blow-up solution for excitatory networks b > 0 when the initial data is concentrated near sufficient to the threshold voltage.This result was obtained by giving no information about the behavior at the blow-up time.In a recent work [12], we state the theorem gives a characterization of this blow-up time when it occurs for b > 0.
There exist a classical solution of ( 2)-( 4) in the time interval [0, T * ) with T * > 0. The maximal existence time T * > 0 can be characterized as Moreover, when b ≤ 0 we have that T * = ∞, while for b > 0 there exist classical solutions which blow up at a finite time T * and consequently have diverging mean firing rate as t ↑ T * .
We now state the main result on steady states from [8].

1.
Under either the conditions b > 0 and 2a then there exists at least one steady state solution to (2)-(4).

2.
If both then there are at least two steady states to solution to (2)-(4).

3.
There is no steady state to (2)-( 4) under the high connectivity condition'

Finite Element Approximation
We construct the numerical approximation of the problem given in ( 2) and ( 3) in two ways: First we use FEM for space discretization that provides a system of ordinary differential equations, which is then solved by Euler's backward difference for time.Spatial discretization involves the construction of a weak formulation of problem over a given domain Ω = [v 0 , v n ] with specified boundary conditions at v = v 0 and v = v n .Weak formulation of the problem (2) and ( 3) is obtained by multiplying the equation with some test function w(v) and integrating over Ω, Ω w(v) ∂p ∂t In the present study, the mean firing rate N(t) is approximated by using backward FDM.Performing the integration by parts in above equation, we get the following equation: This resulting integral ( 8) is called weak formulation because it allows approximate function with less continuity (or differentiability) than the strong form given in Equation ( 2).Once we obtained the weak formulation, next step is to discretize the weak form for the easy representation and to capture the local effects more precisely.Weak form discretization consists of dividing the entire domain into set of elements, then developing the finite element model by seeking the approximation of a solution over a typical element.This discretization is tackled by taking n non-overlapping elements say D i = [v i , v i+1 ] for i = 1, 2, . . ., n with step size h given by : The unknown function p(v, t) must be approximated in a manner so that continuity or differentiability demands by weak formulation can be met.Since the weak formulation contains the first order derivative of p, any function with non-zero first derivative would be a candidate for approximation.Thus, semi discretization consists of finding where pj are the nodal values and ψ j are the basis functions given by There are many choices of weight function w(v) to be used.Particular choice for the weight function w(v) in Galerkin approach is the same as the choice of basis function ψ j (v).Thus, substituting weight function w = ψ l (x), l = i, i + 1 and approximation for the solution defined by Equation (10) in weak formulation obtained in Equation ( 9) leads to the following equation: On simplification of the above, we get as follows Solving Equation ( 12), we get the system of ordinary differential equations for p = ( p0 , p1 ) T , which can be expressed in matrix notation given by where , where By assembling the contribution from all elements, we get the following system for the global nodal vector p = [p 0 , p 1 , . . ., p n ] T on the entire domain where f is the column vector with all entries are zero except at reset potential V R .Ordinary differential Equation ( 14) requires implicit and stable time-stepping method to avoid extremely small time-step.
Firstly we discretize the time domain [0, T] into m subintervals with time step ∆t.We use the Euler's backward difference in time and get the following system from Equation ( 14): this algebraic system (15) can be solved for p j .

Stability Analysis of the Scheme
Fourier method is a very flexible tool for analyzing stability developed by von Neumann.In this method, initial data is demonstrated in terms of finite Fourier series and we examine the growth of individual Fourier component.After assembling the reduced system (13) and using b 1 = b N(t k−1 ) and b 2 = a(N(t k−1 )), we get the finite element difference-differential equation at the i−th node given by where Let us denote the error as e k−1 j at the (k − 1)th stage for nodal value j to satisfy Equation ( 16) At each time level, error can be expanded as e k j = A ξ k e iβ j h .Thus, substituting the error in Equation ( 17), we get as follows: Simplification of the above equation, we get the following Performing some algebraic manipulation, we get Whenever Q ≥ 0, we get |ξ| ≤ 1, hence, the scheme is conditionally stable.

Numerical Experiments
In this section, we present some numerical examples to demonstrate the behavior of the solutions of the nonlinear NLIF model.The performance of the developed scheme is tested by comparing our results with the existing scheme in the literature.Consider the system (2) and (3) with initial data as follows where change in mean v 0 and variance σ 2 0 describe different scenario of a solution.We notice that the behavior of the solution depends upon the value of excitatory (b > 0) and inhibitory (b < 0) average network.First we take constant diffusion coefficient a(N) = a 0 and find the effect on solution with change in value of b.
In Figure 1a, we find that after some time the solution p(v, t) goes to steady state by taking excitatory case i.e., b = 0.5 > 0 small enough.(20) with mean v 0 = 0 and variance σ 2 0 = 0.25.The system is considered for the excitatory case i.e., by taking b = 0.5 > 0 with activity dependent noise a(N(t)) = a 0 = 1. Figure 1a depicts the numerical solution p(v, t) using FEM at different time level using initial data (20); Figure 1b shows the comparison of the existing scheme and FEM for numerical Solution p(v, t) at t = 1.5.
The approximate solution p(v, t) for v ∈ [−4, 2] at different time levels t > 0 is plotted in Figure 1 with a reset potential V R = 1.From Figure 1a, we see that height of impulse decreases as time increases and after some time it reaches a steady state.In Figure 1b, we perform numerical approximation based upon FEM and compare the results obtained in [8] at a final time t = 1.5.
The evolution of firing rate N(t) with the time t > 0, is plotted in Figure 2. We find that firing rate has different range with change in b > 0 excitatory case as well as the inhibitory case b < 0. In Figure 2a, we consider the case when the initial data is centered at v 0 = 0 with b = 0.5.We observe that the solution reaches a steady state.We also take the cases when the initial data is centered at v 0 = −1 and different values of b = 3, 1.5, −1.5, to find different phenomena based on these values in Figure 2.
The errors for the approximate solution simulated in Figures 1, 3 and 4 are plotted in Figures 5-7 respectively.Numerical values of these errors and CPU time (MATLAB and Statistics Toolbox Release 2013a, The MathWorks, Inc., Natick, MA, United States) of the two methods are shown in Tables 1 and 2.
where p * is the numerical solution at finest grid N * .For the numerical experiments, we find the errors with the finest grid being N * = 320.Our test outcomes show that the error defined above is sure a monotone decreasing function as N increases i.e., N = 20, 40, 80, 160.In Figure 3a, we performed numerical approximation based upon FEM and compared the results obtained in [8] at a time t = 0.0408.The evolution of firing rate N(t) with the time t > 0, is plotted in Figure 3b, which describes the blow-up situation when the initial data is concentrated around v 0 = 1.5 with b = 1.5 > 0. Tables 3 and 4 gives the different error values at the final time t = 0.0408 for the same data that is graphically represented in Figure 6.From Figure 4a, it is clear that when initial data is concentrated enough near the threshold point V F , solution blows up at a finite time t = 0.0025, which is earlier than the phenomena described in Figure 3b.This happens because initial data is concentrated enough near the threshold point i.e., v 0 = 1.83, with b = 0.5 > 0 small enough.For different error values and CPU time for the data graphically represented in Figure 4, see Tables 5 and 6 .In Figure 8, we treat the cases for a(N) = a 0 + a 1 N(t), a 0 , a 1 > 0 type activity dependence noise.Figure 8a shows that by taking b = 0.5 and a(N(t)) = 0.5 + N(t)/8 solution goes to steady state.This further indicates that solution goes to a steady state earlier than the solution behavior provided in Figure 8b by taking a(N(t)) = 0.4 + N(t)/100.Figure 8c shows blow-up situation of a solution by taking b > 1 and a(N(t)) = 0.5 + N(t)/8.From Figure 8d, we find that by reducing the noise factor a(N(t)) = 0.4 + N(t)/100, solution goes to steady state.
The behavior of the solution by taking noise factor a(N(t)) = 1 + N(t)/100 for both the cases of excitatory and inhibitory are represented in Figure 9. From Figure 10a, we find that solution blows up in a finite time and right figure indicate the situation of steady state by reducing the noise factor.20) with mean v 0 = 1.5 and variance σ 2 0 = 0.005 and with activity dependent noise a(N(t)) = 1 + N(t)/100.Left: The system is considered for the excitatory case i.e., by taking b = 0.5 > 0; Right: The system is considered for the inhibitory case i.e., by taking b = −0.5 < 0.

Conclusions
In this article, we proposed a finite element method to find the approximate solution of the nonlinear NLIF model.The performance of the proposed method is validated by comparing with an existence scheme in the literature.The approximate solutions determined by Galerkin finite element method have same accuracy as achieved by high-order finite difference scheme (WENO-FDM).The proposed scheme takes less computational time as compared to WENO-FDM.The reason behind that the existing scheme contains many computational factors such as smoothness indicator functions and non-negativity weights etc.Moreover, we also included the role of both excitatory and inhibitory impulses in the model equation.The stability analysis of the proposed scheme is discussed which shows that the scheme is conditionally stable.The behavior of the solution is plotted by taking some test examples.The results reveal that the continuous Galerkin FEM is better than the WENO-FDM for simulating dynamics of large-scale neuronal networks in the brain.

Figure 1 .
The approximate solution p(v, t) for initial data

Figure 3 .
Left: The approximate solution p(v, t), Right: Firing rate N(t) for initial data given in Equation(20) with mean v 0 = 1.5 and variance σ 2 0 = 0.005.The system is considered for the excitatory case i.e., by taking b = 1.5 > 0 with activity dependent noise a(N(t)) = a 0 = 1.(a): Comparison of the existing scheme and FEM for the numerical Solution p(v, t); (b): Time evolution of the firing rate N(t).

Figure 4 .
Left: The approximate solution p(v, t), Right: Firing rate N(t) for initial data given in Equation(20) with mean v 0 = 1.83 and variance σ 2 0 = 0.003.The system is considered for the excitatory case i.e., by taking b = 0.5 > 0 with activity dependent noise a(N(t)) = a 0 = 1.(a): Comparison of existing scheme and FEM for numerical Solution p(v, t); (b): Time evolution of the firing rate N(t).

Figure 5 .
Figure 5. Error for the approximate solution p(v, t) plotted in Figure 1.

Figure 6 .
Figure 6.Error for the approximate solution p(v, t) plotted in Figure 3a.

Figure 7 .
Figure 7. Error for the approximate solution p(v, t) plotted in Figure 4a.

Figure 9 .
The approximate solution p(v, t) for initial data given in Equation (

Table 1 .
Error table using FEM for the approximate solution p(v, t) graphically represented in Figure1, at a final time t = 1.5 with the finest grid being N * = 320.

Table 2 .
Error table using finite difference-WENO scheme for the approximate solution p(v, t) graphically represented in Figure1, at a final time t = 1.5 with the finest grid being N * = 320.Errors of the numerical solution p(v, t) are calculated in different norms • 1 , • 2 , • ∞ , norms which are defined as follows

Table 3 .
Error table using FEM for the approximate solution p(v, t) graphically represented in Figure3, at a final time t = 0.0408 with the finest grid being N * = 320.

Table 4 .
Error table using WENO-FDM for the approximate solution p(v, t) graphically represented in Figure3, at a final time t = 0.0408 with the finest grid being N * = 320.

Table 5 .
Error table using FEM for the approximate solution p(v, t) graphically represented in Figure4, at a final time t = 0.00255 with the finest grid being N * = 320.

Table 6 .
Error table using WENO-FDM for the approximate solution p(v, t) graphically represented in Figure4, at a final time t = 0.00255 with the finest grid being N * = 320.