ON THE RETURN PROCESS WITH REFRACTORINESS FOR A NON-HOMOGENEOUS ORNSTEIN-UHLENBECK NEURONAL MODEL

An Ornstein-Uhlenbeck diffusion process is considered as a model for the membrane potential activity of a single neuron. We assume that the neuron is subject to a sequence of inhibitory and excitatory post-synaptic potentials that occur with time-dependent rates. The resulting process is characterized by time-dependent drift. For this model, we construct the return process describing the membrane potential. It is a non homogeneous OrnsteinUhlenbeck process with jumps on which the effect of random refractoriness is introduced. An asymptotic analysis of the process modeling the number of firings and the distribution of interspike intervals is performed under the assumption of exponential distribution for the firing time. Some numerical evaluations are performed to provide quantitative information on the role of the parameters.

1. Introduction and background.In 1964, in [7] Gernstein and Mandelbrot proposed a model of neuronal activity based on the Wiener process.These authors demonstrated that with a suitable choice of parameters, the histograms of the interspike intervals, experimentally recorded, could be plotted with a good degree of approximation with the average of the first passage time (FPT) for a temporally homogeneous Wiener process.Since then various other models based on diffusion stochastic processes have been proposed to describe the evolution of the neuronal membrane potential.In particular, to take account of the exponential decay exhibited by the membrane potential in the absence of input of any type, in 1971 Capocelli and Ricciardi in [3] proposed a model based on the Ornstein-Uhlenbeck (OU) process.This model has been used widely to describe the activity of a single neuron (see, for instance [2], [9], [10], [14]).More recently attention has been paid to the estimation of the parameters involved in this model (cf.[4], [5], [11]).Several ways exist to derive this model, one of these consists of assuming that the neuron is subject to a sequence of inhibitory and excitatory postsynaptic potentials (PSP's) with constant amplitude that occur according to the Poisson's law.Further, it is assumed that, in the absence of input, the membrane potential decays exponentially to a resting value, called resting potential, with a time constant which we denote by θ.When this constant diverges, the OU model yields to the Wiener model.To describe the spikes trains a return process is built on a non homogeneous OU process in which the effect of random refractoriness is introduced.In this regard we recall that the first attempt to study the effect of refractoriness in a point process was made in [12] and in [15].
In the present work we assume that the inputs, while remaining a constant amplitude, are characterized by time-dependent rates, meaning that some external stimulations are induced on the neuron; so that the involved Poisson process is not homogeneous.In Section 2 the model, based on a non-homogeneous OU process, is introduced.A comparison between the obtained OU model and the corresponding time-homogeneous process is done analyzing the trajectories of the two processes and considering the relative entropy of distributions characterizing the two models.Particular attention is paid to the FPT random variable because it represents the "theoretical counterpart"of the neuronal firing time, so that the FPT's probability density function (pdf) describes the pdf of the firing time.In this regard it should be noted that for the OU process the FPT's pdf is not known in closed form if not for thresholds that are not of particular interest in the neuronal context, nonetheless, for the FPT pdf of the OU process is possible to make use of an asymptotic behavior of exponential type (cf.[8]).To study the train of spikes, in Section 3, we build the return process.It is a continuous process with jumps.The number of firings and the distribution of interspike intervals are studied under the assumption of exponential distribution for the firing time.In Section 4 we introduce random downtimes which delay spikes, simulating the effect of refractoriness.A theoretical and numerical analysis of the return process in the presence of constant and exponential refractoriness is performed.
2. The model.To construct the model, we assume that the neuronal membrane potential is subject to a sequence of inhibitory and excitatory postsynaptic potentials characterized by constant magnitude occurring with time-dependent rates: , where A i (t), A e (t) are positive function of time and σ 2 > 0.Moreover, in the absence of inputs the membrane potential decays to the resting potential with a time constant θ > 0. So, making use of a standard procedure (cf., for instance, [13]), it can be proved that the evolution of the neuronal membrane potential is described via a diffusion process {X(t), t ≥ 0} defined in R whose infinitesimal moments are related to the rates.In particular, the drift and infinitesimal variance of X(t) are respectively, with Note that when θ diverges X(t) becomes a Wiener process with drift µ(t).Our study is led on the model characterized by a generic function m(t).Moreover, to give a quantitative information on the evolution of the membrane potential, we focus on the case m(t) = A sin t because this situation reflects some oscillatory effects of the environment acting on the neuron.In general, X(t) is solution of the following stochastic equation: where B(t) is a standard Wiener process.Eq. ( 2) describes the evolution of the membrane potential.
To analyze the effect of the time dependent drift, in the following we denote by X(t) the process obtained from X(t) when m(t) = 0. Hence, X(t) is described by the following equation: Of course X(t) is a time homogeneous OU process with drift Ã1 (x) = −x/θ + µ and infinitesimal variance Ã2 = σ 2 .In Figure 1 the sample paths of X(t) and of X(t) with periodic m(t) are compared.In particular, we have chosen µ = −14 mV/ms, θ = 5 ms, σ = 1 mV/ms 1/2 , X(0) = X(0) = −70 mV and we have considered various amplitudes of m(t): m(t) = sin t (magenta line), m(t) = 5 sin t (black line) and m(t) = 10 sin t (red line).The sample path of X(t) is flat when it is compared to the others; hence the introduction of m(t) makes the process more fluctuating.Moreover, by increasing A (magnitude of the environmental fluctuations) the sample paths of the process become more and more oscillating.Note that X(t) can be obtained from X(t) via a transformation as shown in the following Proposition.Proposition 1.The process with is a homogeneous OU process characterized by drift and infinitesimal variance Proof.Let v[X(t), t] = X(t), from the Ito's lemma we have that X(t) satisfies the following stochastic equation where It follows that is the Ito's equation for the diffusion process X(t) characterized by infinitesimal moments (6).

2.1.
Transition probability density function.The transition pdf of X(t) is a normal density:  with mean and variance respectively.Note that if m(t) is such that exists and it is finite, then the steady state density of X(t) is: Alternatively, when m(t) is a periodic function with period Q it is possible to consider lim n→∞ f (x, t + nQ|y, τ ) = w(x, t); if this limit exists w(x, t) plays a role analogous to steady state density.
To analyze the influence of m(t) on the transition pdf, we consider the relative entropy between f (x, t|y, τ ) and f (x, t|y, τ ): where with M (t|y, τ ) = y e −(t−τ )/θ + µθ[1 − e −(t−τ )/θ ], represents the transition pdf of X(t).The relative entropy, although it is not symmetrical, is used as a measure of the distance between f e f .Making use of ( 9) and ( 14), from (13) we have In Figure 2  2.2.First passage time problem.Let S ∈ R be a state of process X(t) representing the firing threshold.Let be the FPT through S and let g(S, t|y, τ ) = dP (T y ≤ t) dt be FPT pdf.The random variable T y describes the time of occurrence of neuronal spike and g(S, t|y, τ ) is the theoretical counterpart of the firing pdf for the neuron.Making use of the Proposition 1 we can study the FPT problem of X(t) through S via the FPT problem of X(t) through In particular, denoting by g[ S(t), t|ỹ, τ ] the FPT pdf of X(t) from ỹ = y + d(τ ) through S(t), we have that g(S, t|y, τ ) = g[ S(t), t|ỹ, τ ]; so we focus on g[ S(t), t|ỹ, τ ].
Unfortunately, the FPT pdf g is known analytically only in particular cases that aren't of interest in the present context.However, numerical approximations for FPT pdf can be obtained via appropriate algorithms (cf., for example, [1]).Furthermore, since X(t) admits steady state density under large assumptions, one can prove that the FPT pdf exhibits an exponential behavior for large times when the boundary is far from the starting point (cf.[8]).
In particular, two cases can be distinguished: 1. if S(t) admits limit and S = lim t→∞ S(t), then for large times one has: where where with U (t) = lim n→∞ S(t + n Q) and w(x) defined in (17).In the present context m(t) = A sin(t), so that U (t) = S + d(t) with d(t) given in (8), hence one has: In Figure 4 the approximation of the FPT pdf g[ S(t), t|τ ] obtained via ( 19) is plotted for S(t) = −60 + d(t) with m(t) = sin(t) for µ = −14, θ = 5 and σ = 2, while in Figure 5 we have chosen σ = 3.  3. The return process.In this Section the return process constructed on X(t) is studied to describe the train of spikes characterizing the neuronal activity.Note that in the Case 1. of Section 2 the firing threshold is asymptotically constant so the firing time mean is time independent.This situation has been widely studied in [6], so that in the following we will analyze the return process constructed on X(t) making use of the exponential approximation of the FPT pdf (19).
3.1.Description of the process.Let Z(t) be the return process constructed from X(t) in the following way.Starting at Z(0) = X(0) = η, the process goes on as described by (4) until S(t), defined in (16), is reached for the first time.After this time the process is instantaneously reset to η and then evolves as described by (4) until S(t) is reached again.
The process Z(t) consists of recurrent cycles I 1 , I 2 , . . . of random duration described by I 1 , I 2 , . ... The random variable I i represents the i-th interspike interval (ISI); indeed I 1 is the waiting time for the first firing and, for i = 2, 3, . .., I i measures the time interval elapsing between the (i − 1)-th and the i-th firing.
Each sample path of Z(t) is solution of the following stochastic equation: where B(t) is a standard Wiener process, P (t) is a non-homogeneous Poisson process with intensity λ(t) given in (20) and the term S(t) − η represents the amplitude of the jumps.In this context, each return occurs simultaneously with a jump and it represents a neuronal spike.Figures 6-9 show the sample-paths of Z(t) assuming that the process returns on the starting point η when the boundary S(t) is reached.In Figure 6, we have σ = 1 and S(t) = −50 + d(t) with d(t) given in (8), so the distance between the starting point η = −70 and the threshold is approximatively 20; moreover η is different from the equilibrium point of the process being µθ = −50.In Figure 7 the boundary has been approached to η = −70, simultaneously the width of the environmental oscillations has been increased, consequently the frequency of the jumps increases.In Figures 8 and 9

Analysis of interspike intervals (ISI).
To analyze the sequence of the spikes we denote with T = ( T1 , T2 , . . ., Tn ) the random vector that represents the instants of time in which single firings occur.Note that each Ti is a FPT, moreover the variables Ti conditioned from Ti−1 = t i−1 are independent and distributed according to g[ S(t), t|t i−1 ].Now we study the joint density of T.
For t 0 < t 1 < t 2 < . . .< t n it follows:  From ( 22) the marginal density and the distribution function of the i-th element Ti of T can be determined.Indeed one has and, after k < i integrations, one has:  hence, it results: From (23) the marginal distribution function of Ti can be determined.In particular one has: placing x = Λ t0 (t) so that dx = λ(t)dt, it follows: We use (24) to determine the probability of occurrence of k firings up to time t.
To this aim, we denote by M (t) the stochastic process that counts the number of firings records in (t 0 , t).It results that M (t) is a Poisson process with intensity Λ t0 (t).Indeed, To analyze the ISI distribution we assume that T0 = t 0 and for n = 0, 1, . . .we denote with I n+1 = Tn+1 − Tn the random variable describing the duration of the (n + 1)-th ISI, i.e. the duration of the time interval between the n-th and (n + 1)-th firing.One has: so, assuming that (19) holds, it follows that Therefore, I n+1 conditioned from Tn = t n (n = 0, 1, . ..) is distributed according to an exponential law with parameter Λ tn (t n + x) and its pdf is 4. The effect of refractoriness.The refractoriness is the time interval of variable duration that follows a spike during which the neuron is incapable of responding to input signals.We introduce refractoriness periods in the return process so that the interspike's interval, starting from the second one, can be considered as consisting of the sum of two terms: the first one represents the refractory period following the firing, the second term describes the time for firing from the state η.Therefore we construct a new process Z R (t) describing the evolution of membrane potential in the presence of refractoriness as follows.Starting at Z R (t 0 ) = X(t 0 ) = η, a firing takes place when X(t) attains for the first time the firing threshold S(t), after which the neuron is unable to fire again for a period of refractoriness of random duration.At the end of this period, Z R (t) is instantaneously reset to η.The subsequent evolution of the process then goes on as described by X(t), until the boundary is again reached.A new firing then occurs, followed by the period of refractoriness, and so on.In Figures 10 and 11 sample paths of Z R (t) are showed.The process Z R (t) consists of recurrent cycles F 0 , R 1 , F 1 , R 2 , . . .each of random duration.The duration of cycle F i is represented by the random variable F i described by the FPT pdf of X(t) through S(t) starting from η.Moreover, for i = 1, 2, . .., the refractory period R i is represented by the random variable R i that we assume to have pdf φ ti (t) depending on the time t i at which the last spike is occurred.In particular, for i = 0, 1, . .., the duration F i of F i denotes the time interval elapsing between the i-th reset of the membrane potential at the value η and the (i + 1)-th FPT from η to S(t).Instead, for i = 1, 2, . .., R i indicates the duration of the i-th refractory period.We note that the random variables F i aren't iid because they depend on the instant of the last reset as well as R 1 , R 2 , . . .that depend on the last firing time.
The random variables that describe the ISI are given by: To study the ISI pdf's we denote by T = ( T 1 , T 2 , . . ., T n ) the random vector that represents the instants of time in which single firings occur during the evolution of Z R (t).Note that the variable be the ISI distribution conditioned by occurrence of the last firing at the time t n .
Proposition 2. For n = 1, 2, . . . the ISI conditional distribution is: and the ISI conditional pdf results: Proof.To obtain the distribution F In+1| Tn (x|t n ) we note that x represents the width of the interval (t n , t n + x).So, after the instant of the n-th firing, occurred at time t n , there is a refractory period that can have width at most x and ends at time r.
In the remaining interval (r, t n + x) another firing occurs.It follows that: from which, recalling (20), Eq. (30) immediately follows.By taking the derivative with respect to x, we get: that leds to (31).
In the following we consider two types of refractoriness.In the first case the refractoriness is constant and its duration is 1/ξ, whereas in the second case we consider a refractoriness period of random duration characterized by exponential distribution with parameter ξ, so that its mean duration is the same as that of the constant case.4.1.Constant refractory period.We assume that the refractoriness period is constant and of duration 1/ξ with ξ > 0. Since in this case the refractoriness is not random we can describe its pdf via the Dirac delta function denoted by δ(•).So, assuming that the last spike occurs at the instant τ , one has: Thanks to Proposition 2, the ISI distribution can be evaluated.In particular, from (30), recalling (33) and making use of the properties of the Dirac delta function, it follows that: with Λ τ (t) defined in (20).Hence, one has: where H(x) is the Heaviside unit step function.Furthermore, the ISI pdf is: In Figure 12 the ISI pdf in the presence of constant refractoriness is plotted.We have chosen ξ = 0.1 (red line) which represents a refractoriness of 10 milliseconds and ξ = 1 (black line).Note that when the refractoriness period is longer (red line) the ISI pdf assumes higher values at times greater than 10 milliseconds.Indeed, for a fixed time x ≥ 10, the spikes are more frequent for shorter refractoriness, consequently the distribution of the ISI decreases when ξ increases.The ISI pdf without refractoriness f In+1| Tn and the ISI pdf in the presence of constant refractoriness f In+1| Tn can be compared.Indeed, in Figure 13 are plotted f In+1| Tn and f In+1| Tn for refractoriness of duration 1/ξ = 10 ms, assuming that the last spike occurs at the same time t n : note that the ISI pdf with refractoriness is not shifted with respect to that without refractoriness, as you would expect.This is due to the assumption that the two densities are evaluated for the same instant t n .To point out this property we consider the ISI pdf without refractoriness for t n + 1/ξ as showed in Figure 14 where f In+1| Tn (x|t n ) and f In+1| Tn (x|t n + 1 ξ ) are compared.4.2.Exponential refractory period.We suppose that the refractory period is a random variable having an exponential distribution of parameter ξ.Hence, assuming that the last firing occurs at time τ one has: To determine the ISI distribution we make use of the Proposition 2. In particular, from (30), recalling (37), one has: from which it follows: where Λ τ (t) is defined in (20).Moreover, the ISI density in the presence of exponential refractoriness is In Figure 15 the ISI density f In+1| Tn (x|t n ) in the presence of exponential refractoriness, obtained from Eq. ( 40), is plotted for ξ = 0.1 (red line) and for ξ = 1 (black line).Note that for small amplitudes (small values of x) the ISI pdf for ξ = 1 exceeds the density for ξ = 0.1 and then this behavior reverses, differently from the case of constant refractoriness (cf. Figure 12).Indeed, since the refractoriness mean is 1/ξ ms, one has that for small values of x, corresponding to small ISI durations, it is more likely that the neuron with smaller mean refractoriness fires.As in the case of constant refractoriness we compare the ISI densities with and without refractoriness.In Figure 16 are plotted f In+1| Tn (x|t n ) (blue line) and f In+1| Tn (x|t n ) (red line) with ξ = 0.1.Note that, also in this case, the curves are not shifted, whereas this property is verified when we compare f In+1| Tn (x|t n + 1 ξ ) and f In+1| Tn (x|t n ) at least for great values of x as we can see in Figure 17.

5.
Conclusions.In this paper we have analyzed the membrane potential activity of a neuron subject to some external stimulations.To this aim, we have considered an OU process characterized by a time-dependent drift in which there appears a function m(t) that represents the additional external input acting on the neuron.The result has been a non-homogeneous OU model.We have considered a periodic function m(t), as this situation reflects some oscillatory effects of the environment acting on the neuron.To analyze the ISI distribution a return process has been considered.In this process we have introduced random downtimes which delay  spikes, simulating the effect of refractoriness.The expression of the ISI distribution has been obtained.This distribution is conditioned by the time in which the last fire occurs.We have focused on constant and exponential refractoriness characterized by the same mean value.Some similarities between the ISI pdf with refractoriness and without refractoriness have been observed.In particular, our analysis has showed that the ISI pdf in the presence of refractoriness is shifted with respect to the ISI pdf in the absence of refractoriness provided the latter is suitably conditioned.This observation supports the proposed model.
Future research may investigate the behavior of the model with different external inputs (different choices of the function m(t)) as well as different refractoriness distributions.

Figure 3 .
Figure 3.For the same parameters of Figure 2, D(f | f ) is plotted with m(t) = sin t (magenta line), m(t) = 5 sin t (black line) and m(t) = 10 sin t (red line).
for µ = −14, θ = 5, σ = 2 we plot the time functions f (−65, t| − 70, 0) (blue line) and f (−65, t| − 70, 0) with m(t) = A sin t by choosing A = 1 (magenta line), A = 2 (black line) and A = 3 (red line).The function f (−65, t| − 70, 0) fluctuates around f (−65, t| − 70, 0) and the width of the fluctuations increases as A increases.To analyze the distance between f and f , in Figure 2 the relative entropy (15) is plotted for the same parameters of Figure2.Note that the relative entropy vanishes at times in which f = f , moreover it increases as A increases in according to the behavior shown in Figure2.
we have chosen η = µθ = −70 and S(t) such that | S(t) − η| ≈ 10.Note that in Figures 8 the boundary is rarely reached although σ = 3; increasing σ the frequency of the jumps increases as shown in Figures 9 where σ = 4.

Figure 10 .
Figure 10.A sample path of ZR(t) is plotted for the same parameters of Figure 6, with a constant refractoriness of 5 millsec.

Figure 11 .
Figure 11.A sample path of ZR(t) is plotted for the same parameters of Figure 8, with a constant refractoriness of 5 millsec.

Figure 12 .
Figure 12.For the same parameters of Figure8and tn= 5 ms the ISI pdf's in the presence of constant refractoriness, given in (36), are plotted for ξ = 1 (black line) and ξ = 0.1 (red line).

Figure 13 .
Figure 13.The ISI pdf's f I n+1 | Tn (x|tn) (blue line) and f I n+1 | Tn (x|tn) with ξ = 0.1ms (red line) are plotted for the same parameters of Figure 12.

Figure 15 .
Figure 15.The ISI pdf given in Eq. (40) with ξ = 1 (black line) and with ξ = 0.1 (red line) are plotted for the same parameters of Figure 8 and for tn= 5 ms.