GAUSS-DIFFUSION PROCESSES FOR MODELING THE DYNAMICS OF A COUPLE OF INTERACTING NEURONS

With the aim to describe the interaction between a couple of neurons a stochastic model is proposed and formalized. In such a model, maintaining statements of the Leaky Integrate-and-Fire framework, we include a random component in the synaptic current, whose role is to modify the equilibrium point of the membrane potential of one of the two neurons and when a spike of the other one occurs it is turned on. The initial and after spike reset positions do not allow to identify the inter-spike intervals with the corresponding first passage times. However, we are able to apply some well-known results for the first passage time problem for the Ornstein-Uhlenbeck process in order to obtain (i) an approximation of the probability density function of the interspike intervals in one-way-type interaction and (ii) an approximation of the tail of the probability density function of the inter-spike intervals in the mutual interaction. Such an approximation is admissible for small instantaneous firing rates of both neurons.

1. Introduction.In the last decade an increasing attention has been dedicated to the explanation of the dynamics of couples of interacting neurons [6,7,16,17], of little and large networks of neurons [14,18,19] by using different stochastic models.Some of these models are essentially based on Stein equations [6,7], jump-diffusion processes [18] or Gaussian processes [1] subject to Poissonian inputs originated by the surrounding neurons.More recently a different approach is presented in [15] in which the dynamics of pairs of neurons is described through copulas.
Inspired by these papers, by combining previous investigations on the dynamics of pairs of neurons [6,7] and by some analytical results related to the first passage time problem for the linear time-inhomogeneous Ornstein-Uhlenbeck (OU) process [3,4,8,12,13], which we refer as the generalized OU process, we now analyze the response of the model sketched in Fig. 1, that seems to be able to describe the interactions between the two neurons.
With reference to the left side of Fig. 1, the synaptic current, responsible for the variation of the membrane potential of a neuron (say Neuron 1), is subject to the influence of a constant voltage generator.Such a device turns on as a consequence of a post-synaptic potential due to one of the two neurons (say Neuron 2) and stays on until an action potential of Neuron 1 occurs.The same phenomenon occurs for Neuron 2. Nothing prevents that only one of the two neurons is effectively influenced by the other: it happens, for example, when one of the generators produces a zero voltage (as shown, with respect to the Neuron 1, in the right side of Fig. 1).
2. The model.Denoting by V S a deterministic level representing the neuronal potential threshold of both neurons and by V i = V i (t) : t ≥ 0 the sub-threshold membrane potential time course of Neuron i (i = 1, 2), 1 we consider the following stochastic differential equations: 1) Whenever the membrane potential of one of them achieves V S , it is usual to say that the neuron fires or a spike occurs; after a spike the membrane potential of the corresponding neuron is reset to a preassigned value, say v 0 , below the potential threshold.
In Eqs.(1), W 1 = W 1 (t) : t ≥ 0 and W 2 = W 2 (t) : t ≥ 0 are two independent standard Brownian motions, while B i,3−i = B i,3−i (t) : t ≥ 0 is a stochastic process whose value at time t depends on σ-algebra σ({V 1 (s), V 2 (s) : 0 ≤ s ≤ t}).Appropriate regularity conditions on a i (t), B i,3−i and σ 2 i (t) are required.In the sequel, in order to describe the time evolution of the membrane potential of two interacting neurons, as far as for Leaky Integrate-and-Fire (LIF) models for a single neuron, the following functions are considered: As usual, θ > 0, ρ and µ represent the membrane decay time constant, the membrane resting potential and a constant current (due to an excitatory-inhibitory synaptic balance or/and to an external stimulus), respectively.With the aim to specify the process I i,3−i = I i,3−i (t) : t ≥ 0 we have to introduce some additional settings.Let T i,0 = 0, we denote by {T i,n } n∈N0 the sequence of spike instants of the Neuron i and by the set (eventually empty) of its firing times belonging to a generic interval between two successive spikes of the Neuron 3 − i.It is now useful to introduce the sequence {L i,n } n∈N : Therefore, {I i,3−i (t) : t ≥ 0} describes a synaptic current (see, for instance, in [11,16,17]) obeying: where, In Eqs.(4), k 1 and k 2 represent the intensities of the two neurons interaction:3 excitatory (inhibitory) in the case of positive (negative) values; while zero intensity represents the case in which the corresponding neuron is not affected by the spikes of the other one.Furthermore, α is the decay time constant of the synaptic current.Finally, as for V i (T + i,n ) = v 0 , even the synaptic current is subject to a reset to a preassigned value, say i 0 , after each spike, i.e.I i,3−i (T + i,n ) = i 0 (n ∈ N).With the aim to make the contribution of the second term on the right hand side of Eqs.(4) dependent on the time elapsed since the last spike of the same neuron, we take into account the following (not continuous) solution Fig. 2 is a qualitative illustration of the features of the above described model.Some useful remarks are listed here.
1.With regards to positions in (5), one recognizes that 2. An interval between two successive spikes of the same neuron is called interspike interval (ISI).The dynamics and the reset positions described above do not allow us to identify each ISI of Neuron i with the first passage time (FPT) of the process V i = V i (t) : t ≥ T i,n , starting at v 0 at time T i,n = i i,n , through a threshold V S > v 0 : in fact, when the Neuron i fires the neuronal potential of the Neuron 3 − i continues its evolution without any reset.
v S ,1(t ) In the middle, the synaptic current acting on Neuron 1.On top, a path v 1 (t) of the process V 1 and its firing threshold.Note that, after a spike of Neuron 2 the synaptic current decreases, reducing the oscillation ) = 0 and, for t > t 1,1 , one has h 1 (t) = 0 and h 2 (t) = 1 or h 1 (t) = 1 and h 2 (t) = 0.
3. We consider the values of the sequences {L i,n } n∈N defined in (3) as the "reference times": in such times there is a change of the course of the synaptic currents of Neuron 3 − i. 4. To take in account the case of a time dependent threshold, in the present model we have to introduce a prototype function V S (t) by which a stochastic threshold {V S, i (t) : t ≥ 0} can be defined as follows: The extension to the case of different values of v 0 , i 0 , θ, α, ρ and µ and different V S (t) for the two neurons is easy to accomplish.In the sequel, we consider them coincident in order to make easier the required calculations.
3. First passage time for the OU process.Here we consider one single neuron, by describing its sub-threshold membrane potential, {V (t) : t ≥ t 0 },4 within the framework described above.To this purpose, we keep unchanged meaning and symbols (if necessary without the subscript used as a reference of the considered neuron) of the quantities introduced before.Assuming the absence of the synaptic current (i 0 = 0), the equation for the sub-threshold membrane potential is the following with V (t 0 ) = v 0 .It is well known that the (Gaussian-Diffusion) OU process V = V (t) : t ≥ t 0 with mean and variance respectively, is the (unique) solution of Eq. ( 7).Hence, m V (t|v 0 , t 0 ) and C V (τ, t|v 0 , t 0 ) completely specify the distribution of any order; in particular, the transition probability density function coincides with that of a Gaussian random variable with as mean and variance, respectively.Let S(t) (t ≥ t 0 ) be the firing threshold of the considered neuron, in the sequel assumed as a continuous function.Then, the first passage time of {V (t) : is an (fair) absolutely continuous random variable, for which we denote by its probability density function. 5or S(t) ∈ C 2 ([t 0 , +∞[), in this specific context, the following function has a particular relevance: Indeed, the following integral equation holds (see, [2]): Note that, by choosing S(t) = (ρ + µθ) + c 1 e (t−t0)/θ + c 2 e −(t−t0)/θ , with c 1 and c 2 arbitrary constants, one has a null kernel of the Eq. ( 8) for all t 0 ≤ τ ≤ t.In this case, the probability density function of the FPT is equal to the known term in Eq. ( 8).However, this result is not relevant in the considered biological context.Eq. ( 8) is very important, since via a numerical procedure a good approximation of the required function can be obtained by it (see, [8]).The main drawback, due to computational complexity of the considered quadrature rules, is related to the determination of the tail of g V [S(t), t|v 0 , t 0 ] when S(t) m V (t|v 0 , t 0 ) for all t > t 0 .In this situation, if lim t→+∞ S(t) = S (asymptotically constant threshold), the following result, adapted by [10], can be fruitfully used.
for t − t 0 max{β, θ} one has: 4. One-way interaction.In the present section we intend to determine the distribution of the inter-spike interval T i of the two neurons subject to the dynamics sketched on the right side of Fig. 1.For this purpose, in Eqs.(6), it is sufficient to require k 1 = 0 and k 2 = 0.
In such a case, the spikes of the Neuron 2 do not affect the membrane potential of the Neuron 1 and the ISI's distribution of Neuron 1 is equal to the first passage time distribution of the corresponding membrane potential through the firing threshold.This cannot be stated for the ISI's distribution of the Neuron 2. In order to evaluate such a distribution it is necessary to set the appropriate initial conditions.Specifically, we have to observe the dynamics until both neurons fire and then impose the initial conditions in these times.To fix ideas, we choose two times, 0 < t 1 < t 2 , such that V 1 (t 1 ) = v 0 , V 2 (t 2 ) = v 0 and no spike of Neuron 1 occurs between t 1 and t 2 .
4.1.Distribution of T 1 .Here, we refer to X 1 = X 1 (t) : t ≥ t 1 as the OU process of Section 3, having infinitesimal variance equal to σ 2 1 and initial condition X 1 (t 1 ) = v 0 .For t > t 1 , Eq. ( 1), written for i = 1, is linear with a degenerate-type initial condition.This remark allows us to state that V 1 , conditioned to stay in v 0 at time t 1 , is a Gauss-Diffusion process having mean value i 0 e −(s−t1)/α e (s−t1)/θ ds = m X 1 (t|v 0 , t 1 ) + m 1 (t|t 1 ), (9) where Furthermore, since the covariance of {V 1 (t) : t ≥ t 1 } is determined only by means of −1/θ and σ 2 1 , we can say that it is equal to C X 1 (τ, t|v 0 , t 1 ).It follows that therefore, given that m 1 (t 1 |t 1 ) = 0 ⇒ X 1 (t 1 ) = V 1 (t 1 ) = v 0 , one has: Finally, it results: Note that, with characteristic time constant max{α, θ}, we have so, let be if Theorem 3.1 applied to the threshold S(t) = V S − m 1 (t|t 1 ) provides, for t − t 1 max{α, θ}, the following approximation: 4.2.Distribution of T 2 .Now, we refer to X 2 = X 2 (t) : t ≥ t 2 as the OU process of Section 3, having infinitesimal variance equal to σ 2 2 and initial condition X 2 (t 2 ) = v 0 , with t 1 < t 2 and no other spike of Neuron 1 occurs in [t 1 , t 2 ].Here with T (1) 2 we denote the first ISI of Neuron 2 after t 2 .Then, by conditioning with respect to F t := σ({V 1 (s), V 2 (s) : 0 ≤ s ≤ t}) the process H 1 = H 1 (t) : 0 ≤ s ≤ t is known and Eq. ( 1), written for i = 2, is linear with a degenerate-type initial condition.Therefore, under F t2 , the process V 2 , conditioned to stay in v 0 at time t 2 , is a Gauss-Diffusion process with mean value where, Since the process H 1 is measurable, one has: 1 − e −(s−t2)/α e (s−t2)/θ P 1,2 (s|t 1 , t 2 ) ds, where, So that {V 2 (t) : t ≥ t 2 } is a Gauss-Diffusion with mean and covariance C X 2 (τ, t|v 0 , t 2 ).It follows that Thereby, since m 2 (t and Since Then, by setting max{α, θ}, the following approximation: Note that ( 18) and ( 19) are challenging to use due to the numerical evaluations of P 1,2 (s|t 1 , t 2 ) in m 2,1 (t|t 1 , t 2 ) and h V 2 (in which the calculation of l 2 is required).Finally, we explicitly observe that the probability density function of T 2 is the average of g V 2 (V S , t|v 0 , t 2 ) with respect to the probability density function of the random variable describing the time interval elapsing from the last occurrence of a Neuron 1 spike and t 2 .

4.
3. An approximation of the distribution of T 2 .Let now suppose These statements allow us to apply Eq. ( 14) for t ≥ t 2 : the asymptotic regimen of the membrane potential of Neuron 1 is guaranteed with respect to states and times by the first and the second one, respectively.Thus, for the function P 1,2 (s|t 1 , t 2 ) defined in (16), for s ≥ t 2 , one has: By virtue of this approximation in Eq. ( 15), we obtain: Therefore, by setting we can say that the process V 2 = V 2 (t) : t ≥ t 2 is a Gauss-Diffusion process that, after t 2 , approximates V 2 and Finally, one has: We emphasize that the function m 2,1 , given in analytical form in (21), does not depend on t 1 , so that: and the approximation (22) also holds for the probability density function of T 2 .In Figs. 3 and 4 the case of Neuron 2 with a special inhibitory synapse is considered; the values of the involved parameters are indicated in the corresponding caption.We obtained the simulated ISIs by virtue of Eqs.(1): it is evident that spikes of Neuron 1 are able to decrease the firing activity of Neuron 2. Finally, on the right side of Fig. 4, the agreement between the ISIs' histogram of Neuron 2 with the approximated density function seems to be very satisfactory.5. Mutual interaction.Here we limit ourselves to present a result holding in an appropriate asymptotic regimen for both neurons; extensions to a more general case will be object of our future studies.Being in the framework described in the previous section, the idea that we want to develop begins with the already highlighted remark that in Eq. ( 21) the initial time t 1 disappears.Since, Then, by setting ), for t − t 2 max{α, θ}, provides the following approximation: Now, from Eq. (24) it appears that h V 2 is totally unrelated to the evolution of the membrane potential of Neuron 1. Accordingly, setting by exchanging the role of Neuron 1 with the one of Neuron 2, if for t − t 1 max{α, θ}, we obtain: In Fig.
6. Summary.The idea developed in the present article consists in modeling the interaction between two neurons by using two generalized coupled stochastic Leaky Integrated-and-Fire equations, each one describing the stochastic evolution of the neuron membrane voltage (for biological background see, for instance, [5]).The coupling is realized by including a function for the synaptic current that jumps when the other neuron fires.We are able to determine two Gauss-Diffusion processes suitable to describe the above dynamics by obtaining their mean and covariance function applying classical rules for calculation of the conditional expected value of the processes.The Eq. (15) shows that the mean of one of such processes involves the distribution of the FPT of the other one.
By proceeding along this line, one can determine an approximation of the FPT distribution of each process, by solving a system of non singular second-type Volterra integral equations via a numerical procedure.
We explicitly note that, in this interactive scheme, it is not possible to identify the ISI of each neuron with the FPT of the respective above stochastic process: the reason essentially stays in the fact that the reset of one of the two neurons occurs (every time) after a different time spent from the last firing time of the other neuron.The complete solution of the problem could be obtained by averaging the FPT distribution with the distribution of this random time.Our future research will be devoted to deal with such a problem.
However under the hypothesis that one of the two membrane potentials (e.g. that of the Neuron 1) stays in a particular asymptotic regime (in short: its equilibrium level has to be quite below the firing threshold) we obtain in the present paper an approximation of the FPT distribution of the membrane potential of the Neuron 2 (see Eq. ( 22)) that does not depend on the last firing time of the Neuron 1: in such a case the distribution of the ISIs of the Neuron 2 can be identified with the determined FPT distribution (see Eq. ( 23)).
Finally, we compare our numerical approximation of the FPT probability density function with the histogram of ISIs simulated by applying the numerical Euler scheme to the involved stochastic differential equations (1).A quite satisfactory agreement between the results is shown in several graphics for the case of one-way interaction and in the case of mutual interaction between the two neurons.

Figure 1 .
Figure 1.Schemes of two interacting neurons.Inside each box, the arrows represent fluxes of synaptic current.The vertical one takes in account all the other synapses, while the horizontal one the synapses of the interacting neuron.

Figure 2 .
Figure 2. On the left: Qualitative description of the model for Neuron 1; lowercase letters indicate the sample paths of the involved stochastic processes.On time axis, firing times of each neuron and reference times of (5) are shown.On bottom, plot of the function which indicates the firing times of Neuron 2.In the middle, the synaptic current acting on Neuron 1.On top, a path v 1 (t) of the process V 1 and its firing threshold.Note that, after a spike of Neuron 2 the synaptic current decreases, reducing the oscillation point of function v 1 (t) (k 1 < 0), h 2 (t + 1,. ) = 0, and h 2 (t + 2,. ) = 1.On the right: Same description with reference to Neuron 2 (k 2 < 0).Note that h 1 (t + 1,. ) = 1, h 1 (t + 2,. ) = 0 and, for t > t 1,1 , one has h 1 (t) = 0 and h 2 (t) = 1 or h 1 (t) = 1 and h 2 (t) = 0.

Figure 4 .
Figure 4. ISIs' histogram and probability density function of FPT for Neuron 2. On the left: With no interaction (k 2 = 0).On the right: The interaction with Neuron 1 is considered (k 2 = −1).Other parameters are chosen as in Fig. 3. Continuous lines depict numerical evaluation obtained by virtue of Eq. (22).
5, with the choice of parameters indicated in the caption, tails of probability density function of FPT for both neurons and corresponding ISIs' histograms obtained via simulation of Eqs.(1), are shown.The agreement seems to be quite satisfactory.