Some Fractional Stochastic Models for Neuronal Activity with Different Time-Scales and Correlated Inputs

: In order to describe neuronal dynamics on different time-scales, we propose a stochastic model based on two coupled fractional stochastic differential equations, with different fractional orders. For the specified choice of involved functions and parameters, we provide three specific models, with/without leakage, with fractional/non-fractional correlated inputs. We give explicit expressions of the process representing the voltage variation in the neuronal membrane. Expectation values and covariances are given and compared. Numerical evaluations of the average behaviors of involved processes are presented and discussed.


Introduction
The fractional calculus ( [1,2]) used in the framework of stochastic modeling (see [3] and references therein) allows a twofold advantage: a mathematical generalization of existing stochastic processes and the construction of more refined models, including much more detailed phenomenological evidence.Specifically, in neuronal modeling, the introduction of fractional non local operators in stochastic equations allows the description of confluent dynamics on different time scales.This is the idea on which this contribution is written: to provide a fractional stochastic model that generalizes the previous ones, but it is also suitable to describe dynamics such as those involved in the neuronal activity subject to interactions and inputs of different nature.
In the scientific literature of recent decades, we can see a very large number of contributions devoted to refine existing stochastic neuronal models, among them the well-known Leaky Integrate-and-Fire (LIF) model, in order to enrich the description with further phenomenological aspects.The purpose of such refinements is, essentially, that to include adaptation behaviors, correlated inputs, correlation between inputs, evolution on different time-scales, etc. (see, for instance, [4][5][6][7]).The stochastic LIF model, essentially based on a stochastic Langevin equation ( [8]), has become unsatisfactory: it only includes a deltacorrelated white noise, whereas correlated inputs (see, for instance, [9]), in some cases modeled by colored noises ( [10]), appear more satisfactory to explain the high variability in the neuronal response to inputs and adaptation phenomena ( [9,11]).Furthermore, its intrinsic Markov property, or equivalently, memory-less property, makes it not suitable to model adaptation dynamics, or any other firing activity affected by memory effects.Moreover, only one time-scale can be considered, even if different time-scales seem to be mandatory when it is required to include the effects of the activity of ionic channels and/or synaptic inputs on the neuronal firing.In some previous contributions (see [12][13][14][15] and references therein), some specialized LIF models were provided in order to satisfy such kind of requirements.
Inspired by papers [10,16], we proposed in [15] a LIF model with a colored noise by using a time-non-homogeneous Gauss-Markov (GM) process.Approximations of the first passage time (FPT), whose probability density is not available in closed form, for such a process through a constant level (the firing threshold) was obtained and discussed with some others such as those in [10,17].Successively, some further investigations about firing rates were provided, and specifically in [18], a modified LIF-type stochastic model was considered with a non-delta correlated input process in place of the white noise.Mathematical tools and simulation approaches were adopted in [18] to perform an analysis of that LIF-type model with specified correlated inputs.The theory of integrated GM processes ( [12] ) was evoked in the study of such kinds of models.Indeed, the neuronal voltage was modeled by means of the integral over time of a correlated input, an Ornstein-Uhlenbeck (OU) process.In [12], advances in the study of fractional integrals applied to GM processes were made in order to devise more general models, also in the context of the neuronal modeling.Finally, with the aim to put together the strategy to include the correlated input ( [15,18]) and the use of fractional integrals ( [12]), in [14], we proposed a first neuronal model based on the fractional integral of a correlated OU process.The joint effect of these two approaches, separately considered in previous papers, was discussed in [14].
The motivation of the present paper is that to provide a neuronal model improved respect to that presented in [14] and oriented to gain a mathematical generalization.From this improved version new specialized neuronal models can be derived.Here, our specific purpose is to study a system of two coupled fractional stochastic equations by which it is possible to describe the dynamics of the neuronal membrane voltage by means of the solution process subject to an input modeled by another process, possibly evolving on a different time-scale.
Specifically, the neuronal activity is modeled by the process denoted by V(t) as the α-fractional time integral of the input process, denoted by η(t), evolving on a different more or less finer time-scale tuned by a fractional order β.The values of the fractional orders α and β can be chosen on the basis of neuro-physiological conditions and evidence.
The novelty of the proposed model is in the use of fractional derivatives and integrals in place of the classical ones for modeling both the voltage of the neuronal membrane and a correlated input.Indeed, due to the non-locality property, the use of fractional operators is suggested by the need to model a sort of "memory"in the evolution of the neuronal membrane voltage and in the input itself.Models including memory effects are essentially based on fractional Langevin equations, or Langevin equations with colored noise or correlated inputs ( [10,17,18]).Here, in particular, we consider fractional Langevin equations and a correlated input.
The present paper is organized as follows: in Section 2, we define the general model based on two fractional coupled stochastic equations.Modeling justifications and biological interpretation are discussed in Section 2.1.The theoretical assumptions for the existence of solutions are specified in Section 2.2, which is dedicated to the mathematical setting for the general fractional model ( [2,19]).Then, in successive sections, we derive specific models for different choices of the involved functions and the fractional orders α and β.Specifically, we consider three cases: Model 1 (without leakage) in Section 3, Model 2 (with leakage) in Section 4, and finally Model 3, which is the general FNM model, in Section 5.For each case, we give the explicit form of the stochastic process V(t) solving the proposed equations and the corresponding mean and covariance functions.In Section 6, we provide numerical evaluations of the mean of the process V(t) for the three considered models.In particular, we show and compare the average behavior of the solution processes subject to different inputs and for different values of fractional orders.Finally, in Section 7, we write some conclusions and indicate possible strategies for future simulations of sample paths of such processes and their FPT, useful to estimate neuronal firing times.

The General Fractional Model
Here, we denote by FNM the fractional neuronal model, the new model based on the following coupled stochastic differential Equations (SDEs): for t > 0, α, β ∈ (0, 1] and σ, γ (positive) intensity parameters of the stochastic terms η and dW,, respectively.(In Section 2.1, we provide the meaning of all involved functions and parameters).Here, D α and D β denote the fractional Caputo derivative ( [1,3]), that, for all f ∈ C 1 , is defined as follows: with ν = α, β ∈ (0, 1],, respectively.Some mathematical notes: The fractional differential Equation (1) can also be considered by assuming the existence of the process V(t) such that: where, for a f continuous function, I α ( f ) is the fractional Riemann-Liouville (RL) integral ( [20]) of order α defined as follows with the Gamma Euler function (Γ), i.e., Γ(z) = +∞ 0 t z−1 e −t dt , z > 0. (Recall that we use with the same meaning the double notation I α ( f )(t) = I α ( f (t)), even if the first one denotes mainly the integral operator I α ( f )).
The above fractional integrals evaluated on stochastic processes have to be intended in pathwise interpretation.And, note that the SDE (1) is well defined due to the fact that the Caputo fractional derivative D α is the right inverse of the RL integral of α order (see [20]), i.e., for f (0) = 0, justified by considering that D α ( f ) = I 1−α ( f ′ ) and by taking into account that (see [3] for details) Similarly, we assume the existence of an input process η(t) such that: Here, the stochastic process I β (dW(t)) is defined as the Lévy fractional Brownian motion, i.e., we assume that where W is a standard Brownian motion and the integral has to be interpreted in the Itô sense.We remark that only if β ∈ 1 2 , 1 , for fixed t > 0, the function (t − τ) β−1 is in L 2 (0, t).Thus, the above process is well-defined and it is a Gaussian process starting from 0 at time 0 with probability 1.
Furthermore, the integral forms ( 4) and ( 6) allow to apply Theorem 4 of [21], by which we obtain Gaussian solutions.
Finally, we note that a more general model could have been considered if in place of (1) we put the following SDE with noise terms involving W and W independent Brownian motions but, in this paper, we consider the special case of FNM model ( 1) and (2), i.e., we exclude the noise in the first equation in such a way the dynamics of V(t) is the fractional integral (not perturbed by any noise) of the correlated input η(t).We postpone the study of the FNM model (8) to a future work.

Some Biological and Modeling Justifications
The adoption of the two above coupled fractional equations allows to use a α-fractional stochastic process V(t) for modeling the neuronal evolution of the membrane voltage and a β-fractional stochastic correlated input η(t) integrating a current term I.In this sense, we intend that the model can describe neuronal activities on different time-scales modulated by parameters α and β.
In order to clarify the notation, and to show how this model born as a generalization of the previous ones, we firstly recall the essentials on the neuronal models under consideration.In particular, we recall that the SDE on which the classical stochastic LIF model is based, which is the following stochastic version of the Langevin equation: The solution process V(t) stands for the voltage of the neuronal membrane, whereas C m is for the membrane capacitance, g L is for the leak conductance, V L is a specified "resting" value of potential, dW(t) is the white noise (with W(t) a standard Brownian motion (BM)), σ is the intensity (or amplitude) of the noise.Such terminology, proper of electrical circuits, is justified by the representation of the neuronal membrane similar a capacitor with capacitance C m in parallel with a resistor (see, for instance, [6,8]).It is possible to have the so-called leakage, because the charge inside the neuronal cell, in absence of the electrical input current, slowly leak through the neuronal membrane towards to the resting potential value V L , or value after leakage.The leak conductance g L is such that 1/g L stands for the leak resistance; moreover, the ratio C m g L denotes the membrane time constant.
The occurrence of a neuronal spike, or equivalently a large depolarization of the voltage when this attains a critical value V th , is modeled by the FPT of the process V(t) through the threshold value V th .The firing activity corresponds to the emission of sequences of spikes with different rates, and it is the basis of the information transmission between neurons.This mechanism is triggered by electro-chemical variations due to entrance and/or exit of potassium and sodium ions through ionic channels (gates) of the neuronal membrane, but also due to the injections of external currents or synaptic inputs originated by surrounding neurons.It appears clear that these dynamics occur on different time scales.Finally, the membrane voltage, which acts as an integrator of different inputs, can be modeled by means of the solution of a linear differential equation.
More specifically, Equation ( 2) is a fractional version of the stochastic LIF model ( 9).Here, we adopted (2), and its solution process η, for modeling a correlated input that is fractionally integrated by the process V(t), solving Equation (1).
We again remark that the use of fractional operators is made to include memory in the resultant model.Moreover, the choice of different fractional orders, for the voltage V(t) and the input η(t), is made for distinguishing the different time-scales of the above-explained dynamics.
Furthermore, in the present model, the η process can be view as a synaptic input coming from other neurons and it depends on an additional current term denoted by I. Indeed, I models a synaptic or an injected current embodied in the correlated input η(t).It is possible to perform different mathematical choices for I: beyond the cases of excitatory or inhibitory currents corresponding to positive or negative I terms, we can choose a constant, a periodical function, and a stochastic process [7,18,22], each one with a proper biological justification.Indeed, the form of I can reflect different physiological circumstances; anyway, here, even if we consider a constant I, the suggested above dynamics corresponds to the activity (modeled by V(t)) of a neuron affected by multiple interactions originated from the activity of neighborhood neurons, of correlated inputs, and from codified currents on different time-scales, all summarized in η(t).
A first step in such modeling direction is made in [14], where the proposed model is essentially constructed on the two SDEs: There, the stochastic correlated input is modeled by the Gaussian process η(t), and a deterministic input current I is involved in its drift.In [14], the functions F and G was suitably specified.In the case, F(t, V L , η) is independent on V(t), the solution process V(t) is essentially the fractional integral of an OU-type process and its FPT was estimated by simulations in order to provide approximations of firing times in [14].Hence, it is now clear that the FNM model ( 1) and ( 2) constitutes the generalization of ( 10) and (11).Indeed, if in (2) we assume β = 1, the fractional derivative D β reduces to the classical derivative and the input process η(t) belongs to the class of Ornstein-Uhlenbecktype processes that are properly Gauss-Markov correlated processes.In such a case, η(t) is the same considered in (11).
In this way, the membrane potential V(t) is the α-fractional Caputo integral of the correlated input η(t).
Here, the purpose is to study the fractional models derived as particular cases of the proposed fractional FNM model based on the stochastic pair (V(t), η(t)), solutions of ( 4) and (6), respectively.Note that the correlated input process η(t) drives the dynamics of membrane voltage V(t) in different ways, depending on the choice of involved functions f , g.In addition, we also justify the adoption of Equations ( 1) and (2), or equivalently, the coupled processes (4)-( 6), remarking that the fractional order α of the derivative in (1) can allow us to describe dynamics of V(t) on a time-scale different from that of the input η(t), regulated by β.We expect that the possibility to estimate from the data the values of α and β makes actually applicable the model and adherent to the real phenomenon.

The Mathematical Setting for Fractional SDEs and Solutions
Mathematical results about the resolution of fractional SDEs are scarce and fragmentary.Some papers, including [23,24], are mainly dedicated to the existence and uniqueness of the solutions.
We found (1) and ( 2) really appropriate for solving SDEs with the method of variation in the constant formula for Caputo fractional differential equations, as described in [25].Indeed, the assumption of the standard Lipschitz condition on the coefficient of considered SDEs is satisfied.
Remark 1.We remark that if B(t) is a stochastic process with sample paths almost surely (a.s.) Hölder continuous (as those we will consider), a stochastic version of FDe can be defined similarly to Equation (12) in path-wise sense.Note that (1) of proposed model FNM is a stochastic version of the fractional differential equation as (12) with f (t, V(t), V L ) ≡ A and ση(t) = B(t).
In the next theorem, we provide the explicit form of the solution of FDe (12).
Theorem 1.The solution of (12) in [0, T] is: The proof can be found in [27].Note that in case of a stochastic version of FDe (12) with B(t), the stochastic process of a complete filtered probability space (Ω, F , {F t } t∈[0,∞) , P) with a.s.Hölder continuous paths, the explicit solution given in (18) is the corresponding stochastic process solution.
The definition of a Caputo-fractional SDE follows.
In [24], it was proved that for any given X 0 , there exists a unique solution, whose explicit form was found with a variation constant formula.We recall such solution in the next theorem.
Theorem 2. The solution process of the fractional SDE (14) The proof can be found in [25].
Really useful for our application purposes is also the following corollary.
Corollary 1.Consider the fraction SDE: The explicit expression of the solution of ( 19) is: Note that the case α = 1 can also be covered by previous theoretical results, taking into account that in such a case the fractional derivative and integral reduce to the corresponding classical ones.Consequently, it is easy to prove that the above results become exactly equal to those well-known in the ordinary and stochastic differential calculus.Remark 2. Under all above assumptions we choose the coefficients of fractional SDE here considered, in such a way solutions exist, are unique, and their explicit expressions are also available.The validity of above mathematical results in a compact set [0, T] is perfectly adequate for such kind of models, because the dynamics have to be described until the finite time instant T corresponding to the firing time of the neuron, object of future investigations.
Moreover, about the equations on which the FNM model relies, we remark that (1) is a stochastic version of a FDe as (12) with stochastic term B(t) ≡ η(t), whereas (2) is a typical fractional SDE as defined in (14).
In next sections, by specifying the involved functions, we identify the coefficients of equations and we provide the explicit expressions of the corresponding solutions.

Model 1:
The FNM Model without Leakage with a Correlated Input: α ∈ (0, 1), β = 1 Beyond the choice of fractional orders, we have to choose the functions f and g involved in the FNM model ( 1) and (2).The case of was extensively studied in [14] as an integrate and fire (IF) neuronal model without leakage.
Here, specifically, we consider that model as a particular case of FNM model ( 1) and (2) with α ∈ (0, 1), β = 1, essentially based on the two following equations: There, the meaning of the involved parameters g L , C m , V L is the same previously specified in Section 2.1, whereas the addition parameter τ is related to a characteristic time involved in the correlation time of the input process η(t), (see [15] for details).
In this first instance, we note that Equation ( 21) includes the input process η(t) in place of a typically noise term.
In particular, for the specified involved functions, we note that Equation ( 21) is a stochastic version of a FDe as (12) with A = 0, B(t) = g L V L C m + ση(t), whereas ( 22) is a typical a typical SDE, that can be viewed as the fraction one defined in (19) for α = 1 and with A ≡ − 1 τ , B(t) ≡ I τ , σ(t) ≡ γ.

The input process η(t)
We firstly specify that η(t), form the well-known SDE (22), is a non-stationary Gauss-Markov (GM) Ornstein-Uhlenbeck (OU) process.Its mean function is: and covariance Furthermore, it can be written in terms of a standard BM process W(t) with the transformed time r(t), as follows: where r(t) = k(t) h(t) is the increasing function, and functions k(s), h(t) are such that the covariance (24) can be factorized as c(s, t) = k(s)h(t).Hence, from (24), we have Finally, an explicit form of η(t) ( [12,14]) is: from which it is possible to understand how the stochastic input process η(t) depends, in particular, on its initial value η 1,0 and on the constant input current term I.
The voltage process V 1 (t) From ( 21) and the last form of process η(t) we can also specify the form of the process V 1 (t).Indeed, first we re-write (21) by using (27) or, equivalently, and finally, by applying Theorem 1 to Equation ( 21), we can explicitly write: Furthermore, in this, case the process V 1 (t) is Gaussian (see, for instance, [21]) with mean where is the Mittag-Leffler function ( [26]).Its covariance can be derived by considering that the last term in RHS of ( 30) is the fractional integral of an OU process and from [12], we can specify its covariance as follows: In [12] it was proven that with Remark 3. We note that the process V 1 preserves also in its covariance the effect of the input η(t) under the fractional integration.This appears evident by writing (32) equivalently as The analysis carried out constitutes the proof of the following proposition.
Proposition 1.The voltage process V 1 (t) solution of (21) with the correlated input process η(t) solution of ( 22) is a Gaussian process with mean (31) and covariance (32).
In [12], a simulation strategy, also applied in [14], for such a process and its FPT can be found.

Model 2:
The FNM Model with Leakage and a Correlated Input: α ∈ (0, 1), β = 1 As suggested in [14], if in the model ( 10) and (11), that corresponds to FNM model ( 1) and ( 2) for α ∈ (0, 1) and (by inserting the so-called leakage, as in a LIF model), and we have to address the following systems of equations: In particular, for the specified involved functions, we note that Equation ( 35) is a stochastic version of a FDe as (12) C m , whereas (36) is a typical SDE, that can be viewed as the fraction one defined in (19) for α = 1 and with A ≡ − 1 τ , B(t) ≡ I τ , σ(t) ≡ σ/τ.Furthermore, note that in order to gain the integral representations of ( 35) and (36) as in ( 4) and ( 6), we also have to set V 2,0 = 0. Remark 4.More generally, the agreement between a biological circumstance and the choice of initial values V(0), η(0) can be found when it is required to model the firing activity of neurons with or without reset of the neuronal voltage and of the input after each firing time.This phenomenon is typically called the case of endogenous/exogenous setting: specifically, the endogenous setting is related to the reset of voltage and of the input after a spike (a firing time), otherwise the exogenous case is considered.
Hence, the integral form of equations in ( 35) and ( 36) is: and where I is the Lebesgue measure integral and I(dW)(t) ≡ W(t), in such a way we recognize η(t), the solution of the SDE (38), to be an OU process, essentially the same of the previous section, i.e., a Gauss-Markov OU-type process with mean ( 23) and covariance (24).
The voltage process V 2 (t) Referring to [1,14,26,28], and by applying Theorem 1 to Equation (35), the process V 2 (t) has the following expression: Equivalently, by recognizing from (42) that we can also write: Hence, V 2 (t) is Gaussian process (see, Theorem 4 in [21]) that, for V 2,0 = 0, has the following mean where m(t) is the mean of η(t) as specified in (23).Equation (40) shows that the V 2 (t) process involves fractional integrals of a Mittag-Leffler function and of the process η(t).In particular, recalling the explicit expression of η(t) in (27), and setting V 2,0 = 0 in (40), we have that the stochastic term of V 2 (t) can be specifically written as From (43), we distinguish the stochastic terms that contributes to the calculus of the covariance.Indeed, recalling r(t) = 1 τ (e 2t/τ − 1) and setting W r (t) = W 1 τ (e 2t/τ − 1) and the covariance of V 2 (t) can be finally written (in terms of the Brownian motion) as Remark 5.It is fundamental to compare the covariances of the processes V 1 and V 2 in order to understand the insertion of leakage what effect has determined on the voltage process.
First, of all, we note that: and comparing (32) with (45), we understand that the additional function E α,α − g L (t−•) α C m is included in the fractional integral involved in the covariance of V 2 respect to that of V 1 .Then, we can write the covariance (in terms of η(t) process) that put in evidence how the covariance of the voltage process V 2 is affected by the input process η(t).Differently for V 1 , the covariance of the process V 2 does not involve simply the fractional integral of η(t) but it preserves the effect of the input η(t) temperated by the inserted function E α,α in the fractional integrals involved in its covariance.
The analysis carried out in this section constitutes the proof of the following proposition.

Model 3:
The General FNM Model with Leakage and a Fractional Correlated Input: α, β ∈ (0, 1) Finally, we consider fractional SDEs of the model FNM to describe both dynamics with leakage for V(t) and η(t).Indeed, in (1), we set f (t, and in (2), we set g(t, η(t), I) = − η(t)−I τ .For this case, we refer to processes V 3 and y in place of the general ones V and η; hence, we have to address the following systems of equations: In order to apply the mathematical results of Section 2.2, we note that Equation ( 47) is a stochastic version of a FDe (12) C m , whereas (48) is a typical β-fractional SDE as (19) and with A ≡ − 1 τ , B(t) ≡ I τ , σ(t) ≡ σ/τ.Furthermore, note that in order to also exploit the integral representations of ( 47) and (48) as in ( 4) and ( 6), we have also to set V 3,0 = 0 and y 0 = 0. Recall conditions of α, β in order to apply theoretical results of Section 2.2.

The fractional input process y(t)
As in [25], and by applying Theorem 2 and, more specifically, Corollary 1 to Equation (48), the solution process y(t) has the following explicit form: It is a Gaussian process ( [21], Theorem 4) with the following mean: The calculus of its covariance leads to: where the stochastic fractional integral I β is that defined in (17).

The voltage process V 3 (t)
To Equation (47), we can apply Theorem 1 and we obtain the following expression for the process V 3 (t): In particular, for V 3,0 = 0 and, by substituting (49) with y 0 = 0, in place of y(s) in (52), we can specify: The mean is: Finally, the covariance is where the stochastic fractional integral I β is that defined in (17).The covariance of V 3 (t) is given in the above form (55) in order to make it comparable with the covariance (51) of y(t).
The analysis carried out in this section constitutes the proof of the following proposition.
Proposition 3. The voltage process V 3 (t) solution of (47) with fractional correlated input process y(t), solution of (48), with specified initial conditions, has mean (54) and covariance given in (55).

Numerical Results and Comparisons
In this section, we show the plots of mean functions of the input process η(t) and of the processes V 1 (t), V 2 (t) and V 3 (t) for specified values of parameters.In particular, we compare the different average behaviors of the membrane processes in correspondence of positive and negative input current terms I and for different resting voltages V L .The comparisons are also provided in order to understand how different fractional orders of α and β affect the evolution of the modeling processes.
We remark that, by preserving the choice of V i,0 = 0, for i = 1, 2, 3, and η 0 = 0 in an endogenous setting of neuronal voltage and input (see, for instance, [18,29]), we assign the values of all involved parameters as specified in Table 1.Note that such values are suitably re-scaled respect to the physiological ones (see [18]).Fractional order for the input β = 0.8 Comparing Figure 1 with Figures 2-4 In left side of Figure 1, we plot the mean of the input process as in (23) with η 0 = 0, τ = 5, for two values of current term I.The red curve corresponds to an inhibitory current term I = −0.01,and the blue curve corresponds to an excitatory current term I = 0.01.Take into account the increasing (decreasing) behavior on time of the input correlated process η(t) on the consequent change in the behavior of V 1 (t) (in Figure 2) and V 2 (t) (in Figure 3).
The effect of η(t) (of Figure 1 Left) on V 1 (t) and V 2 (t): In contrast to the intuition, the negative current term I = −0.01 in the mean of η(t) (the red curve) causes a positive depolarization (a peak) in the profiles of V 1 (t) and V 2 (t), as shown in Figure 2 Left and Figure 3 Left, respectively, for any considered values of fractional order α.A symmetrical behavior can be observed in the profiles of of V 1 (t) and V 2 (t), as shown in Figure 2 Right and Figure 3 Right, respectively, with a well (an hyper-polarization) in place of the peak, in correspondence of the positive current term I = 0.01 in the mean of η(t) (the blue curve in Figure 1 Left).
These results are also a consequence of the chosen value of the resting potential V L , with opposite sign respect to the sign of the current term I; see the captions of Figures 2  and 3, in which are specified the values of involved parameters, for details.
In the right side of Figure 1, we plot the mean of η(t) as in (50) with η 0 = 0, τ = 5.The values of β are indicated in the figure.For the positive current term I = 0.01 in the mean (50) of η(t), we plot the curves on top of Figure 1 Right.For the negative current term I = −0.01 in the mean (50) of η(t), we plot the curves on bottom of Figure 1 Right.
The effect of η(t) (of Figure 1 Right) with β = 0.8 on V 3 (t): Again, in contrast with intuition, the negative current term I = −0.01 in the mean of η(t) (the bottom curves) causes a positive depolarization (a peak, even if less sharp respect those in V 1 (t) and V 2 (t)) in the profile of V 3 (t), as shown in Figure 4 Left, for any considered values of fractional order α.Symmetrically, the positive current term I = 0.01 in the mean of η(t) (the top curves) causes a hyper-polarization (a negative well) in the profile of V 3 (t), as shown in Figure 4 Right, for any considered values of fractional order α.Also, in this case, the chosen value of the resting potential V L played a determinant role in the behavior of the membrane process V 3 (t).
Finally, we deduce that the effect of η(t), corresponding to different stimuli modeled by I, on membrane processes V i (t), for i = 1, 2, 3, is more-or-less evident, and it can induce depolarization or hyper-polarization answers, when combined with some parameters of membrane processes, and in particular with the values of the resting potential.
Note that the choice of characteristic time τ = 5 for the η(t) process is suggested by the comparison with the characteristic time C m /g L = 10 of the membrane process V i (t).Indeed, we suppose that a less value for τ could have made less evident the effects on the membrane processes.Keeping this idea in mind, compare the time axis of Figure 1 with those of remaining figures: in such a way, it is possible to determine and compare the variation periods of these functions.This comment is given in order to underline that it is also possible to detect the persistence of the effect of input η(t) on the membrane process.The effect of different fractional orders α: The main comment on these last Figures is related to the effect on the profiles of such means by varying the value of the fractional order α.We note that in all cases, i.e., in Figures 2-4, the same ordering is observable.We refer to the ordering between the means of the same voltage process as α varies.The fact that in each case (and in each corresponding figure) the value 0.8 of α determines the higher mean (or the lower one, in case of a negative sign) is justifiable by the use of the fractional integral of order 0.8 able to preserve much more 'łmemory"of the input effect and for much more time respect to cases corresponding to less fractional orders.
On the basis of this fact, we also adopted the value 0.8 for the fractional order β for the η(t) process considered in Figure 4.

Concluding Remarks
In this contribution, a fractional model is proposed for neuronal dynamics essentially based on two coupled stochastic fractional differential Equations ( 1) and (2).In order to explain how the model can describe the biological phenomenon, some justifications are given in Section 2.
Note that we sometimes disregarded to specify the physical dimension of all involved functions and processes, but we remark that such specification, essential to apply the model, can be found in Table 1 for the re-scaled values, or in [18] for the original values.
We study the model on the basis of the mathematical setting given in Section 2.2; the explicit expressions of the processes, the means and covariances in the three different considered cases are given in Sections 3-5.These results show the rule played of the current term I, the correlated input η(t), perturbed by a white noise, on the voltage process V(t).We also provide figures with plots of means for η(t) and V(t) processes.The numerical results are studied and commented in detail.Hence, we try to clarify how some memory effects can be preserved and detected when such kind of models is applied.
Furthermore, this contribution is also particularly useful to design simulation algorithms with the purpose to obtain estimations for firing times by means of first passage time of these processes through the level V th , the firing threshold.Our aim will be specialize simulation algorithms, in some cases already available, based on the generation of sample paths of Gaussian processes (see, for instance, [12,30]), but also to adopt and specialize other algorithms based on Euler-Maruyama schemes for stochastic fractional differential equation (see, for instance, [31]).Successive investigations about numerical and simulation approximations of FPTs for this model is mandatory and will be object of a future work.

Table 1 .
Assigned values of the parameters for numerical evaluations.