Numerical simulation of fractional Cable equation of spiny neuronal dendrites.

In this article, numerical study for the fractional Cable equation which is fundamental equations for modeling neuronal dynamics is introduced by using weighted average of finite difference methods. The stability analysis of the proposed methods is given by a recently proposed procedure similar to the standard John von Neumann stability analysis. A simple and an accurate stability criterion valid for different discretization schemes of the fractional derivative and arbitrary weight factor is introduced and checked numerically. Numerical results, figures, and comparisons have been presented to confirm the theoretical results and efficiency of the proposed method.

Introduction

The Cable equation is one of the most fundamental equations for modeling neuronal dynamics. Due to its significant deviation from the dynamics of Brownian motion, the anomalous diffusion in biological systems cannot be adequately described by the traditional Nernst–Planck equation or its simplification, the Cable equation. Very recently, a modified Cable equation was introduced for modeling the anomalous diffusion in spiny neuronal dendrites [1]. The resulting governing equation, the so-called fractional Cable equation, which is similar to the traditional Cable equation except that the order of derivative with respect to the space and/or time is fractional.

Also, the proposed fractional Cable equation model is better than the standard integer Cable equation, since the fractional derivative can describe the history of the state in all intervals, for more details see [1,2] and the references cited therein.

The main aim of this work is to solve such this equation numerically by an efficient numerical method, fractional weighted average finite difference method (FWA–FDM).

In recent years, considerable interest in fractional calculus has been stimulated by the applications that this calculus finds in numerical analysis and different areas of physics and engineering, possibly including fractal phenomena. The applications range from control theory to transport problems in fractal structures, from relaxation phenomena in disordered media to anomalous reaction kinetics of subdiffusive reagents [2,3]. Fractional differential equations (FDEs) have been of considerable interest in the literatures, see for example [4-13] and the references cited therein, the topic has received a great deal of attention especially in the fields of viscoelastic materials [14], control theory [15], advection and dispersion of solutes in natural porous or fractured media [16], anomalous diffusion, signal processing and image denoising/filtering [17].

In this section, the definitions of the Riemann–Liouville and the Grünwald–Letnikov fractional derivatives are given as follows:

The Riemann–Liouville derivative of order α of the function y(x) is defined bywhere n is the smallest integer exceeding α and Γ (.) is the Gamma function. If , then (1) coincides with the classical n derivative y((x).

The Grünwald–Letnikov definition for the fractional derivatives of order α > 0 of the function y(x) is defined bywhere means the integer part of and are the normalized Grünwald weights which are defined by .

The Grünwald–Letnikov definition is simply a generalization of the ordinary discretization formula for integer order derivatives. The Riemann–Liouville and the Grünwald–Letnikov approaches coincide under relatively weak conditions; if y(x) is continuous and y′(x) is integrable in the interval [0, x], then for every order 0 < α < 1 both the Riemann–Liouville and the Grünwald–Letnikov derivatives exist and coincide for any value inside the interval [0, x]. This fact of fractional calculus ensures the consistency of both definitions for most physical applications, where the functions are expected to be sufficiently smooth [15,18].

The plan of the paper is as follows: In the second section, some fractional formulae and some discrete versions of the fractional derivative are given. Also, the FWA–FDM is developed. In the third section, we study the stability and the accuracy of the presented method. In section ”Numerical results” numerical solutions and exact analytical solutions of a typical fractional Cable problem are compared. The paper ends with some conclusions in section ”Conclusion and remarks.”

We consider the initial-boundary value problem of the fractional Cable equation which is usually written in the following waywhere 0 < α, β ⩽ 1, μ is a constant and is the fractional derivative defined by the Riemann–Liouville operator of order 1 − γ, where γ = α, β. Under the zero boundary conditionsand the following initial condition

In the last few years, appeared many papers to study this model (3)–(5) [5,19-22], the most of these papers study the ordinary case of such system. In this paper, we study the fractional case and use the FWA–FDM to solve this model.

Finite difference scheme for the fractional Cable equation

In this section, we will use the FWA–FDM to obtain the discretization finite difference formula of the Cable Eq. (3). We use the notations Δt and Δx, at time-step length and space-step length, respectively. The coordinates of the mesh points are x = a + jΔx and t = mΔt, and the values of the solution u(x,t) on these grid points are .

For more details about discretization in fractional calculus see [5].

In the first step, the ordinary differential operators are discretized as follows [23]and

In the second step, the Riemann–Liouville operator is discretized as followswherewhere means the integer part of and for simplicity, we choose h = Δt. There are many choices of the weights [5,15], so the above formula is not unique. Let us denote the generating function of the weights by w(z, α), i.e.,Ifthen (9) gives the backward difference formula of the first order, which is called the Grünwald–Letnikov formula. The coefficients can be evaluated by the recursive formula

For γ = 1 the operator becomes the identity operator so that, the consistency of Eqs. (8) and (9) requires , and for k ⩾ 1, which in turn means that w(z,0) = 1.

Now, we are going to obtain the finite difference scheme of the Cable Eq. (3). To achieve this aim, we evaluate this equation at the intermediate point of the grid

Then, we replace the first order time-derivative by the forward difference formula (6) and replace the second order space-derivative by the weighted average of the three-point centered formula (7) at the times t and twith λ is being the weight factor and is the resulting truncation error. The standard difference formula is given by

Now, by substituting from the difference operators given by (6), (7) and (9), we getPut , and under some simplifications we can obtain the following formwhere

Eq. (17) is the fractional weighted average difference scheme. Fortunately, Eq. (17) is tridiagonal system that can be solved using conjugate gradian method. In the case of λ = 1 and , we have the backward Euler fractional quadrature method and the Crank–Nicholson fractional quadrature methods, respectively, which have been studied, e.g., in [24], but at λ = 0 the scheme is called fully implicit.

Stability analysis

In this section, we use the John von Neumann method to study the stability analysis of the weighted average scheme (17).

The fractional weighted average finite difference scheme (WADS) derived in (17) is stable at under the following stability criterion

By using (18), we can write (17) in the following form

In the fractional John von Neumann stability procedure, the stability of the fractional WADS is decided by putting . Inserting this expression into the weighted average difference scheme (20) we obtainsubstitute by ϕ = (1 − λ)N and divide (21) by e we get

Using the known Euler’s formula we have

Under some simplifications, we can write the above equation in the following form

The stability of the scheme is determined by the behavior of ξ. In the John von Neumann method, the stability analysis is carried out using the amplification factor η defined by

Of course, η depends on m. But, let us assume that, as in [13], η is independent of time. Then, inserting this expression into Eq. (24), one getsdivide by ξ to obtain the following formula of η

The scheme will be stable as long as ∣η∣ ⩽ 1, i.e.,considering the time-independent limit value η = −1 and since , thenFrom the above inequality, we can obtainPut , we findwhich can be written in the formPutone finds that the mode is stable when

Although, M depends on m, it turns out that M tends toward its limit value

In this limit the stability condition isbut from Eqs. (10) and (11) with z = −1 one sees that , so thatsince , replacing by its highest value and since , therefore we find that the sufficient condition for the present method to be stable and this completes the proof of the theorem. □

For , the stability condition (19) can be satisfied under specific values of . We can check this note from the results which presented in Table 1.

Table 1

The absolute error of the numerical solution of Eq. (35).

x	The absolute error
0.1	0.3063 × 10−3
0.2	0.5826 × 10−3
0.3	0.8019 × 10−3
0.4	0.9427 × 10−3
0.5	0.9912 × 10−3
0.6	0.9427 × 10−3
0.7	0.8019 × 10−3
0.8	0.5826 × 10−3
0.9	0.3063 × 10−3

Numerical results

In this section, we present two numerical examples to illustrate the efficiency and the validation of the proposed numerical method when applied to solve numerically the fractional Cable equation.

Consider the following initial-boundary problem of the fractional Cable equationon a finite domain 0 < x < 1, with 0 ⩽ t ⩽ T, 0 < α, β < 1 and the following source termwith the boundary conditions u(0, t) = u(1,t) = 0, and the initial condition u(x, 0) = 0.

The exact solution of Eq. (35) is u(x, t) = t2sin(πx).

The behavior of the exact solution and the numerical solution of the proposed fractional Cable Eq. (35) by means of the FWA–FDM with different values of λ, α, β, Δt, Δx and the final time T are presented in Figs. 1–5.

In Table 1, we presented the behavior of the absolute error between the exact solution and the numerical solution of Eq. (35) at and T = 0.01.

Also, in Table 2, we presented the maximum error of the numerical solution for λ = 0,α = 0.2, β = 0.7, T = 0.1 with different values of Δx and Δt.

Table 2

The maximum absolute error for different values of Δx and Δt.

Δx	Δt	Maximum error
120	130	0.00751
1100	150	0.00716
1150	1100	0.00428
1150	1150	0.00234
1150	1200	0.00095
1200	1250	0.00010

Consider the following initial-boundary problem of the fractional Cable equationwith u(0, t) = u(10, t) = 0 and u(x, 0) = 10δ(x − 5), where δ(x) is the Dirac delta function.

The numerical solutions of this example are presented in Figs. 6–10 for different values of the parameters λ, α, β, Δx, Δt and the final time T.

Conclusion and remarks

This paper presented a class of numerical methods for solving the fractional Cable equations. This class of methods is very close to the weighted average finite difference method. Special attention is given to study the stability of the FWA-FDM. To execute this aim, we have resorted to the kind of fractional John von Neumann stability analysis. From the theoretical study, we can conclude that this procedure is suitable and leads to very good predictions for the stability bounds. The presented stability of the fractional weighed average finite difference scheme depends strongly on the value of the weighting parameter λ. Numerical solutions and exact solutions of the proposed problem are compared and the derived stability condition is checked numerically. From this comparison, we can conclude that the numerical solutions are in excellent agreement with the exact solutions. All computations in this paper are running using Matlab programming 8.

Conflict of interest

The authors have declared no conflict of interest.