Time Fractional Cable Equation And Applications in Neurophysiology

TL;DR

This paper extends the classical cable equation by incorporating a Caputo time fractional derivative, providing analytical solutions that model anomalous diffusion in neurophysiological structures like spiny dendrites.

Contribution

It introduces a novel fractional derivative approach to the cable equation and derives analytical solutions using Laplace Transform, applicable to neurophysiology.

Findings

Analytical solutions expressed with special functions

Applicable to modeling anomalous diffusion in dendrites

Solutions computed and visualized in Matlab

Abstract

We propose an extension of the cable equation by introducing a Caputo time fractional derivative. The fundamental solutions of the most common boundary problems are derived analitically via Laplace Transform, and result be written in terms of known special functions. This generalization could be useful to describe anomalous diffusion phenomena with leakage as signal conduction in spiny dendrites. The presented solutions are computed in Matlab and plotted.

Equations

$${\lambda }^2\frac{{\partial }^2{V}_m(x,t)}{\partial x^2}-\kappa \frac{\partial V_m(x,t)}{\partial t}-V_m(x,t)=0,$$

$$\frac{{\partial }^2{V}_m(X,T)}{\partial X^2}-\frac{\partial V_m(X,T)}{\partial T}-V_m(X,T)=0.$$

$$\frac{{\partial }^2{V}_m(X,T)}{\partial X^2}=(s+1)V_m(X,T),$$

$${\tilde{V}}_m(X,s)=g(s)e^{-\sqrt{s+1}X}.$$

$${\mathcal{G}}_s(X,T)=\frac{X}{\sqrt{4\pi T^3}}{\mathrm{e}}^{-(\frac{{\mathrm{X}}^2}{4T}+\mathrm{T})},$$

$${\mathcal{H}}_s(X,T)={\int }_0^T{\mathcal{G}}_s(X,T^{\prime })d{T}^{\prime }.$$

$$V_m(X,T)={\int }_0^{T}g(T-T^{\prime }){\mathcal{G}}_s(X,T^{\prime })d{T}^{\prime }.$$

$$I=I_0\delta (T)=\frac{-1}{r_i\lambda }\frac{\partial V_m(X,T)}{\partial X},$$

$$I=I_0\theta (T)=\frac{-1}{r_i\lambda }\frac{\partial V_m(X,T)}{\partial X}.$$

$${\tilde{V}}_m(X,s)=\frac{1}{\sqrt{s+1}}{\mathrm{e}}^{-\sqrt{\mathrm{s}+1}\mathrm{X}},$$

$${\mathcal{G}}_m(X,T)=\frac{1}{\sqrt{\pi T}}{\mathrm{e}}^{-(\frac{{\mathrm{X}}^2}{4T}+\mathrm{T})},$$

$${\mathcal{H}}_m(X,T)={\int }_0^T\frac{1}{\sqrt{\pi T^{\prime }}}{\mathrm{e}}^{-(\frac{{\mathrm{X}}^2}{4{\mathrm{T}}^{\prime }}+{\mathrm{T}}^{\prime })}\mathrm{d}{\mathrm{T}}^{\prime },$$

$$\frac{{\partial }^2{V}_m(X,T)}{\partial X^2}-\frac{{\partial }^{\alpha }{V}_m(X,T)}{\partial T^{\alpha }}-V_m(X,T)=0.$$

$$\frac{\partial V_m(X,T)}{\partial T}={}_{RL}{D}_{0,T}^{1-\alpha }\frac{{\partial }^2{V}_m(X,T)}{\partial X^2}-{}_{RL}{D}_{0,T}^{1-\alpha }(V_m-i_e{r}_m).$$

$$\frac{{\partial }^{\alpha }V(T)}{\partial T^{\alpha }}=-(V(T)-V_L)+I_{inj},$$

$$\frac{{\partial }^2{V}_m^{\ast }(X,T)}{\partial X^2}-\frac{{\partial }^{\alpha }{V}_m^{\ast }(X,T)}{\partial T^{\alpha }}=0.$$

$$\begin{array}[]{ll}V^{*}_{m}(X,T)&=\int_{0}^{T}g(T-T^{\prime})\,\mathcal{G}^{*}_{\alpha,s}(X,T^{\prime})\,dT^{\prime}\,,\\ \\ \mathcal{G}^{*}_{\alpha,s}(X,T)&=\frac{1}{T}W_{-{\alpha/2},0}\left(-X/T^{\alpha/2}\right)\,,\end{array}$$

$$W_{\lambda ,\mu }(z):=\sum _{n=0}^{\infty }\frac{z^n}{n!\Gamma [\lambda n+\mu ]},\lambda >-1,\mu \ge 0.$$

$$t^{\mu -1}{W}_{-\nu ,\mu }\left(-x/t^{\nu }\right)÷s^{-\mu }\exp \left(-x{s}^{\nu }\right),0\le \nu <1,\mu >0.$$

$$F_{\nu }(z)=W_{-\nu ,0}(-z),M_{\nu }(z)=W_{-\nu ,1-\nu }(-z),0<\nu <1,$$

$${\mathcal{G}}_{\alpha ,s}^{\ast }(X,T)=\frac{1}{T}{F}_{\alpha /2}\left(X/T^{\alpha /2}\right)=\frac{\alpha }{2}\frac{X}{T^{\alpha /2+1}}{M}_{\alpha /2}\left(X/T^{\alpha /2}\right),$$

$$(s^{\alpha }+1)\widetilde{V_m}(X,s)-\frac{{\partial }^2\widetilde{V_m}(X,s)}{\partial X^2}=0,$$

$$\widetilde{V_m}(X,s)=\tilde{g}(s)e^{-\sqrt{(s^{\alpha }+1)}\cdot X}.$$

$$\tilde{f}(s)=\varphi (s)\cdot \tilde{F}(\psi (s)),$$

$$f(T)={\int }_0^{\infty }F(\tau )G(\tau ,T)d\tau$$

$$G(\tau ;T)÷\tilde{G}(\tau ,s)=\varphi (s)e^{-\tau \psi (s)}$$

$$\varphi (s)=\tilde{g}(s),\psi (s)=s^{\alpha },$$

$$\tilde{F}(s){|}_X={\mathrm{e}}^{-\mathrm{X}\sqrt{\mathrm{s}+1}}.$$

$$G(\tau ,T)={\int }_0^T\frac{g(T-T^{\prime })}{T^{\prime }}{W}_{-\alpha ,0}(-\tau /T^{\prime \alpha })d{T}^{\prime }$$

$$F(T){|}_X=\frac{X}{\sqrt{4\pi T^3}}{e}^{-(\frac{X^2}{4T}+T)}$$

Full text

**Chaos, Solitons and Fractals (2017): Special Issue on Future Directions in Fractional Calculus. Guest Editors: Mark M. Meerschaert, Bruce J. West, Yong Zhou. Published on line 29 April 2017;

http://dx.doi.org/10.1016/j.chaos.2017.04.043 **

MSC 2010: 26A33, 35A22, 44A10, 92C05.

Introduction

The one dimensional cable model is treated in neurophysiology to model the electrical conduction of non-isopotential excitable cells, we remind as example to the textbooks of Johnston and Wu (1994) [8], Weiss (1996) [29] and Tuckwell (1988) [27], and for some mathematical details to Magin (2006) [12]. In particular it describes the spatial and the temporal dependence of transmembrane potential $V_{m}(x,t)$ along the axial $x$ direction of a cylindrical nerve cell segment. The membrane behaviour is summarized by an electrical circuit with an axial internal resistance $r_{i}$, and a transmembrane capacitance $c_{m}$ and a transmembrane resistance $r_{m}$ in parallel, connecting the inner part to the outside [8]. External axial resistance could be eventually included. Transmembrane potential is generated by ionic concentration gradient across the membrane, and is maintained non null at rest (no current) by a combination of passive and active cell mechanisms. Equivalent models can in fact be derived from the Nernst-Planck equation for electro-diffusive motion of ions , see Qian and Sejnowski (1989)[22].

Cell excitation can be caused by electro stimulation of the membrane. The consequent variation in transmembrane potential is transmitted along the cell segment. The resulting differential equation for the transmembrane potential takes the form of a standard diffusion equation with an extra term to account leakage of ions out of the membrane, it results in a decay of the electric signal in space and in time:

[TABLE]

$\lambda=\sqrt{r_{m}/r_{i}}$ and $\kappa=r_{m}c_{m}$ are the characteristic space and time scales of the process, determined by the values of the membrane resistance and capacitance per unit length of the system, see e.g. [12]. For simplicity in the rest of this work, following [12], we will use the dimensionless scaled variables $X=x/\lambda$ and $T=t/\kappa$, so that we consider the equation

[TABLE]

Some interesting quantities to neurophysiology are connected to First kind boundary condition (the Signalling Problem) and Second Kind boundary condition problems. Signalling Problem is interesting to understand how the system evolves when excited at one end with a specific potential profile, Second Kind boundary condition problem is interesting because it can be related to the profile of a current injected across the membrane.

In Signalling Problems the cable is considered of semi-infinite length ($0\leq X<\infty$), initially quiescent for $T<0$ and excited for $T\geq 0$ at the accessible end ($X=0$) with a given input in membrane potential $V_{m}(0,T)=g(T)$. The solution can be derived via the Laplace Transform (LT) approach:

[TABLE]

and the LT of the solution results

[TABLE]

Relevant cases are impulsive input $g(T)=\delta(T)$ and unit step input $g(T)=\theta(T)$ where $\delta(T)$ and $\theta(T)$ denote the Dirac and the Heaviside functions, respectively. The solutions corresponding to these inputs can be obtained by LT inversion [12] and read in our notation

[TABLE]

and

[TABLE]

We refer to $\mathcal{G}_{s}$ to as the fundamental solution or the Green function for the Signalling Problem of the (linear) cable equation in Eq.(2), whereas to $\mathcal{H}_{s}$ to as the step response. As known, the Green function is used in the time convolution integral to represent the solution corresponding to any given input $g(T)$ as follows

[TABLE]

The spatial variance associated to this model is known to evolve linearly in time.

If we consider an impulse or a step current injected at some point $X$ the problem is subjected to the following boundary conditions, specifically

[TABLE]

or

[TABLE]

We consider the adimensional current $I=I_{0}r_{i}\lambda$ and put it to unity for convenience. Appling the impulse in $X=0$ the LT reduces to

[TABLE]

the Green function and the step response function (when a step current is applied in $X=0$) reads, respectively,

[TABLE]

and

[TABLE]

We emphasize that in this standard case the Green function $\mathcal{G}_{m}(X,T)$ is equal to the Green function for the Cauchy problem $\mathcal{G}_{c}(X,T)$, for an infinite cable up to constant coefficients.

The motion of ions along the nerve cells is conditionated by this model, that predicts a mean square displacement of diffusing ions that scales linearly with time. By the way significative deviations from linear behaviour have been measured by experiments, see e.g. Jacobs et al. (1997) [7], Nimchinsky et al. (2002) [19], Santamaria et al. (2006) [23] and Ionescu et al. (2017) [6]. A relevant medical and biological example is the anomalous subdiffusion in neuronal dendritic spines, see Suetsugu and Mehraein (1980) [25], Duan (2003) [2]. Particularly appropriate systems are spiny Purkenje cell dendrites characterized by both spiny and not spiny branches. Spiny branches are in fact characterized by subdiffusive dynamics, while not spiny branches are not. The spatial variance of a diffusing inert tracer (concentration of) in spiny branches evolves as a sub-linear power law of time, and the diffusion with smaller values of the power exponent is associated to higher spine density [23], as spines behave as a trap for the diffusing molecules.

Anomalous subdiffusion can be modelled in several ways introducing some fractional component into the classical cable model. The fractional cable model developed in this section is defined by replacing the first order time derivative in Eq.(2) with a fractional derivative of order $\alpha\in(0,1)$ of Caputo type see e.g. Gorenflo and Mainardi (1997) [3] and Podlubny (1999) [21].

[TABLE]

The solutions of the most relevant boundary problems (Signalling Problem, Cauchy Problem, Second Kind Boundary Problem) are explicitly calculated in integral form containing Wright functions. Thanks to the variability of the parameter $\alpha$, the corresponding solutions are expected to better describe the qualitative behaviour of the membrane potential observed in experiments respect to the standard case $\alpha=1$.

From a mathematical point of view this model is a simple extension to fractional behaviour of the Neuronal Cable Model and it turns to be in some special cases (vanishing initial conditions) equivalent to the equation developed by Henry et al. (2008) [5]:

[TABLE]

which can be derived from a modified Nernst-Planck equation, with diffusion constant replaced by fractional derivatives of Riemann-Liouville type. Other studies consider similar approaches, see e.g. Langlands et al. (2009) [9], Langlands et al. (2011) [10], Liu et al. (2011) [11], Moaddy et al. [16], often concentrating on the Initial Value Problem (Cauchy Problem). The introduction of a Caputo time derivative to model the voltage response in neurons has already been proposed by Teka et al. (2014) [26] to model spiking adaptation for a homogeneous membrane patch, where the space derivatives vanish, named fractional leaky integrate-and-fire model:

[TABLE]

where an external injected current $I_{inj}$ is included.

Beside the apparent simplicity our approach allows to reproduce at least qualitatively the main characteristics observed in experiments [19], [7],[2],[25][23]. Fractional calculus is often used to catch in a parsimonious mathematical description some underlying complex behaviour. We remind that Caputo’s fractional derivative is a non-local operator and for this reason, as pointed out in [26], it could be also introduced to explain behaviours like multiple timescale dynamics and memory effects, related to the complexity of the medium. Further generalizations of this model as well the introduction of external injected current and a more directly relation with ions motion, could be analysed in future, to refine the biological relevance of the model.

Solution of the Signalling Problem via Laplace Transform

The solution of the Signalling Problem has been derived via LT by Vitali and Mainardi [28], however the inversion of the LT solution for Eq. (13) requires special efforts because of the term $V_{m}(X,T)$.

When this term is not present, the resulting equation is the well known time fractional diffusion equation:

[TABLE]

for which the solutions of the corresponding Cauchy and Signalling Problems have been derived in the 1990’s by Mainardi in terms of two auxiliary Wright functions (of the second type) [13, 14]. Specifically for the Signalling Problem, the general solution there provided in integral convolution form reads

[TABLE]

where $\mathcal{G}^{*}_{\alpha,s}(X,T)$ denotes the Green function of the Signalling Problem of the fractional time diffusion equation (Eq.(16)) and $W_{-{\alpha/2},0}(\cdot)$ is a particular case of the transcendental function known as Wright function

[TABLE]

This function, entire in the complex plane, is discussed extensively in the Appendix F of Mainardi’s book (2010) [15] where the interested reader can find the following relevant LT pairs, rigorously derived by Stankovic (1970) [24]:

[TABLE]

Here we have adopted an obvious notation to denote the juxtaposition of a locally integrable function of time $t$ with its LT in $s$ with $x$ a positive parameter. It is worth to recall the distinction of the Wright functions in first type ($\lambda\geq 0$) and second type ($-1<\lambda\leq 0$) and, among the latter ones, the relevance of the two auxiliary functions introduced in [13]:

[TABLE]

inter-related as $F_{\nu}(z)=\nu zM_{\nu}(z)$. Indeed the relevance of both the Wright functions has been outlined by several authors in diffusion and stochastic processes. Particular attention is due to the $M$-Wright function (also referred to as the Mainardi function in the book of Podlubny (1999) [21]) that, since for $\nu=1/2$ reduces to $\exp{(-z^{2}/4)}/\sqrt{\pi}$, is considered a suitable generalization of the Gaussian density, see Pagnini (2013) [20] and references therein.

Then the Green function for the Signalling Problem of the time fractional diffusion equation (Eq.(16)) can be written in the original form provided by Mainardi (1996) [13] as

[TABLE]

where the superscript $*$ is added to distinguish the time fractional diffusion equation from our fractional cable equation, both depending on the order $\alpha\in(0,1)$.

Applying the LT to Eq.(13 with the boundary conditions required by the Signalling Problem, that is $V_{m}(X,0^{+})=0$, $V_{m}(0,T)=g(T)$, we have:

[TABLE]

which is a second order equation in the variable $X$ with solution:

[TABLE]

Because of the shift constant in the square root of the LT in Eq.(23), the inversion is no longer straightforward with the Wright functions as it is in the time fractional diffusion equation (Eq.(16)). Consequently, we have overcome this difficulty recurring to the application of the Efros theorem, see e.g. the book by Graf (2004) [4] that generalizes the well known convolution theorem for LTs. For sake of convenience let us hereafter recall this theorem, usually not so well-known in the literature. The Efros theorem states that if we can write a LT $\widetilde{f}(s)$ as:

[TABLE]

where the function $\widetilde{F}(s)$ has a known inverse LT $F(T)$, the inverse LT can be written in the form:

[TABLE]

where:

[TABLE]

In Eq.(23), LT solution of our Signalling Problem, we thus have:

[TABLE]

and

[TABLE]

Then, having $\widetilde{G}(\tau,s)=\widetilde{g}(s)\,\rm{e}^{-\tau s^{\alpha}}$, thanks to the standard convolution theorem of LTs, we obtain:

[TABLE]

where $W_{-\alpha,0}$ is the F-Wright function, and

[TABLE]

is the solution in Eq.(5) of the standard cable equation (Eq.(2)).

Then, the general solution for the Signalling Problem can be written in terms of known functions:

[TABLE]

Substituting $g(T)=\delta(T)$ in the general solution in Eq.(31) we obtain the Green function for the fractional model (Eq.(13)), shown in Fig.(1):

[TABLE]

When $g(T)=\theta(T)$ we obtain the step response of our fractional cable equation :

[TABLE]

where $\mathcal{H}^{*}_{2\alpha,s}(\tau,T)$ is the step response function for the time fractional diffusion equation. After some manipulations including the change of variable $z=\tau/{T^{\prime}}^{\alpha}$ and integrating by parts after using the recurrence relation of Wright functions:

[TABLE]

and the relation between the auxiliary functions: $F_{\nu}(z)=\nu z\,M_{\nu}(z)$ we may rewrite the step-response solution as:

[TABLE]

where $\mathcal{H}_{\alpha,s}(X,T)$ is the step response function for the standard cable model and $\mathcal{G}^{*}_{2\alpha,c}(\tau,T)$ is the fundamental solution of the time fractional diffusion equation for the Cauchy Problem. The same expression can easier be derived by direct application of the Efros theorem and is plotted in Fig.2.

The Green function for the Cauchy Problem

Consider an infinite cable with boundary conditions $V_{m}(\pm\infty,T)=0$ and initial condition $V_{m}(X,0)=f(X)$. The general solution of the Cauchy problem is related to the Green function $\mathcal{G}_{\alpha,c}(X,T)$ through the following relation:

[TABLE]

$\mathcal{G}_{\alpha,c}(X,T)$ can be derived via Laplace Transform:

[TABLE]

boundary conditions imposes:

[TABLE]

Imposing $\widetilde{\mathcal{G}}_{\alpha,c}(0^{-},s)=\widetilde{\mathcal{G}}_{\alpha,c}(0^{+},s)$ leads to $c_{1}(s)=c_{2}(s)$. Integrating Eq.(13) over $X$ from $0^{-}$ to $0^{+}$ we have:

[TABLE]

the coefficients result:

[TABLE]

the resulting LT of the Green function reads:

[TABLE]

The inversion can be easily performed for $X>0$, thanks again to the Efros theorem, and extended by symmetry respect to the $X$-axes for $X<0$.

Let’s consider $\phi(s)=\frac{1}{s^{1-\alpha}}$, $\psi(s)=s^{\alpha}$, following the theorem we may set $G(\tau,s)=\frac{1}{s^{1-\alpha}}e^{-\tau s^{\alpha}}$ and $F(X,s)=\frac{1}{2\sqrt{s+1}}e^{-X\sqrt{s+1}}$, that have known inverse LT:

[TABLE]

and

[TABLE]

The inverse LT for the Green function is plotted in Fig.3 and reads:

[TABLE]

Response to injected current

An interesting biological problem is to consider an injected current in the system. Transmembrane potential is related to the transmembrane current through the relation $-I=\frac{\partial^{2}V_{m}(X.T)}{\partial X^{2}}$, where the minus sign is due to the direction of the current, in this case flowing inside the cell. Let’s consider a singular point injected current in $X=0$, it takes the form $I(X,T)=I_{0}\delta(X)f(T)$. Integrating from $0^{-}$ to $0^{+}$ we obtain the relation

[TABLE]

We recall the LT for the semi-infinite cable for an initially undisturbed cable:

[TABLE]

At the boundary condition we have:

[TABLE]

if we consider an impulse injection of current in $X=0$ we have I_{0}\delta(T)=-\frac{\partial V_{m}(X,T)}{\partial X}\big{|}_{X=0^{+}}. Applying this condition to the LT we obtain:

[TABLE]

leading to the following Laplace Transformed solution:

[TABLE]

According to the previous derivations it is then straightforward that the inverse LT takes the form:

[TABLE]

represented in Fig.4.

For a generic boundary $I_{0}\widetilde{f}(s)$ we obtain:

[TABLE]

The general solution becomes:

[TABLE]

The solution is symmetric respect to $X$, the problem can be then extended to the infinite cable introducing a factor $1/2$: $\mathcal{G}^{\infty}_{\alpha,m}(X,T)=\frac{1}{2}\mathcal{G}_{\alpha,m}(X,T)$

The extension to the infinite cable case admits also the following generalization, current injection in $X_{0}\neq 0$ is equivalent to shift the cable of the same value $X_{0}$, then:

[TABLE]

When the injected current is a step function we obtain the following LT solution:

[TABLE]

considering $\phi(s)=\frac{1}{s^{1-\alpha}}$, $\psi(s)=s^{\alpha}$ we have $G(\tau,s)=\frac{1}{s^{1-\alpha}}e^{-\tau s^{\alpha}}$ and $F(X,s)=\frac{1}{s\sqrt{s+1}}e^{-X\sqrt{s+1}}$, that have known inverse LT Eq.(55) can be simplyfied to:

[TABLE]

which is shown in Fig.5.

Conclusions

The cable model, fractional or linear, is used to describe subthreshold potentials, or passive potentials, associated to dendritic processes in neurons. The travelling potential is summed up in the center of the cell, called soma, and an action potential is produced when a threshold is exceeded. Anomalous regimes of diffusion can then have a deep impact on the communication strength.

Diffusion results more anomalous, i.e. the fractional exponent $\alpha$ decreases, with increasing spine density [23]. Decreasing spine density is characteristic of aging [7],[2], pathologies as neurological disorders [19] and Down’s sindrome [25], then subdiffusive regimes are in some sense associated to a healthy condition. It has been suggested that increasing spine density should serve to compensate time delay of postsynaptic potentials along dendrites and to reduce their long time temporal attenuation [5].

Looking at our plotted solutions for the fractional cable equation when an impulsive potential is applied at the accessible end it can be noted from Fig.1 that peak high decreases more rapidly with decreasing $\alpha$ at early times, viceversa is less suppressed at longer times, and the cross over time increases with decreasing $\alpha$. Looking at the potential versus time it can also be noted that potential functions associated to lower $\alpha$ last for longer time at appreciable intensity and arrive faster at early times with respect to the normal diffusion case ($\alpha=1$). By the way, when a constant potential is applied at the accessible end we note from Fig.2 that the exponential suppression of the potential along the dendrite is reduced for high $X$ values with respect to normal diffusion. Instead for small $X$ the potential results just slightly more suppressed in the sub-diffusion process. These behaviours can be noticed also for the other cases in Fig.3 and Fig.4 - 5.

From a mathematical point of view the Efros theorem extends the concept of convolution as an integral form that is consistent with a subordination-type integral. However such integral form does not necessary connote a subordinated process, as it has been shown in literature for the generalized grey Brownian motion (ggBm) by Mura and Pagnini (2008) [18] and in a more extended way by Molina-Garcia et al. (2016) [17], but could also be interpreted as a consequence of the random nature of the media in which particles are diffusing. This model can be also read as a generalization of time fractional diffusion processes where mass is not conserved due to leakage. This approach naturally recover the solution for the time fractional case in the limit in which the leakage is put to zero in the integral forms.

In conclusion the presented fractional cable model satisfies the main biological features of the dendritic cell Signalling Problem. With respect to models solved as Cauchy problem, our approach could include specific time dependent boundary conditions, which will allow to reconstruct with accuracy the expected signal at the soma if the model will result capable to predict real data behaviour. Furthermore the solutions can be computed directly, i.e. calculating the integral associated, as well as by Laplace Transform inversion, see e.g. Abate and Ward (2006) [1] without any remarkable issue. Further generalizations of the model as the inclusion of external currents in the equation and comparison with data will be investigated in future works.

Aknowledgements

The work of F. M. has been carried out in the framework of the activities of the National Group of Mathematical Physics (INdAM-GNFM). All the authors acknowledge support by the Italian Ministry of Education and the Interdepartmental Center ”Luigi Galvani” for integrated studies of Bioinformatics, Biophysics and Biocomplexity of the University of Bologna. The authors would also like to thank PhD. G. Pagnini (BCAM, Bilbao, Spain) for valuable comments and discussion and the anonymous referees for their constructive remarks and suggestions that helped to improve the manuscript.