ON A SPIKE TRAIN PROBABILITY MODEL WITH INTERACTING NEURAL UNITS

. We investigate an extension of the spike train stochastic model based on the conditional intensity, in which the recovery function includes an interaction between several excitatory neural units. Such function is proposed as depending both on the time elapsed since the last spike and on the last spiking unit. Our approach, being somewhat related to the competing risks model, allows to obtain the general form of the interspike distribution and of the probability of consecutive spikes from the same unit. Various results are ﬁnally presented in the two cases when the free ﬁring rate function (i) is constant, and (ii) has a sinusoidal form.

1. Introduction.Since the seminal papers by Gerstein and Mandelbrot [19] and Stein [32], many efforts have been directed to the formulation of stochastic models for single neuron's activity aimed to describe the relevant features of the behaviour exhibited by neural cells.We mention the contributions by Ricciardi [27] and Ricciardi et al. [30], and the bibliography therein, as a reference to mathematical models and methods on this subject.
Various researches have been carried out by the authors of this paper on the construction and analysis of models, based on stochastic processes and aimed to describe dynamic systems of interest in different fields.Their research activity has been performed continuously thanks to the precious guidance and support of Professor Luigi M. Ricciardi, to whose unforgettable memory this paper is gratefully dedicated.
Among the numerous investigations performed in biomathematics under his advice and supervision (mainly in neuronal modeling, population dynamics, subcellular stochastic modeling) we recall the following themes: -the characterization of the time course of the neuronal membrane potential as an instantaneous return process (Ricciardi et al. [29]), -the description of neuronal units subject to time-dependent inputs via Gauss-Markov processes (Di Crescenzo et al. [12]), -analysis of the interaction between neuronal units of Stein type based on Monte-Carlo simulations (Di Crescenzo et al. [18]), -stochastic modeling of the evolution of a multi-species population, where competition is regulated by colonization, death and replacement of individuals (Di Crescenzo et al. [13]), -analysis of birth-death processes and time-non-homogeneous Markov processes in the presence of catastrophes (Di Crescenzo et al. [14], [15]), -the study of stochastic processes suitable to describe the displacements performed by single myosin heads along actin filaments during the rising phases (Buonocore et al. [4], [5]).
Along the lines traced by some of the above contributions, in this paper we discuss a suitable extension of a spike train stochastic model to neuronal networks with interacting units.
In several investigations the synaptic inputs that carry the stochastic component of the neuronal activity is modeled by Poisson processes with a fixed spike rate (see Amit and Brunel [1], Bernander et al. [3], Softky and Koch [31], for instance).We recall that the customary assumption based on Poisson processes allows the approximation of the synaptic input of a typical neuron by a stationary uncorrelated Gaussian process due to the superposition of a large number of incoming spikes (hence a sum of many Poisson processes) of either excitatory as well as inhibitory type (see Ricciardi [27]).However, models based on homogeneous Poisson processes fail to capture the relevant feature of the neural activity consisting in the refractory period.See, for instance, Hampel and Lansky [20] for an investigation on parametric and nonparametric refractory period estimation methods.The refractory period is sometimes modeled by means of a dead time, i.e. the time interval following every firing during which the neuron cannot fire again.This leads to a delayed Poisson process, obtained by a step change to the rate of a Poisson process (see Deger et al. [11], Johnson [21], Ricciardi [28]).
Aiming to include the neuronal refractory period and to describe properties of spike trains, another approach has been adopted recently by various authors.It is based on the assumption that the inhomogeneous Poisson process describing the number of neuronal firings has a conditional intensity function expressed as product of the free firing rate function and a suitable recovery function.
We purpose to investigate the spike train model based on the conditional intensity, where the recovery function is aimed not only to include the refractory period, but also to devise the interaction between several excitatory neural units.This is performed via a suitable choice of the monotone recovery function, which is increasing when describes the effect of excitatory neurons and is decreasing when models the refractory period.This scheme allows studying various statistics related to the firing activity, by following an approach analogous to the competing risks model (see Di Crescenzo and Longobardi [16]).In the homogeneous case it is shown that the overall activity of the network exhibits exponentially distributed interspike intervals.In addition, it seems that other suitable choices of the recovery function yield further dynamics, such as the bi-exponential and periodic behaviors investigated by Mazzoni et al. [24].This is the plan of the paper: In Section 2 we describe the background on the conditional intensity function model.Section 3 presents a suitable extension of this model to the case of a network formed by a fixed number of units, in which the recovery function depends both on the time elapsed since the last spike and on the last spiking unit.A comprehensive discussion on this model is also given, with attention to the conditional random variables describing the time length between consecutive spikes.A connection with the competing risks model is also pinpointed.Section 4 is devoted to investigate the model in detail.We determine the general form of the interspike distribution and of the probability of consecutive spikes from the same unit.Explicit expressions are thus obtained in the special case of constant free firing rate function, when the interspike distribution is shown to be exponential.We also consider the case when the free firing rate is of sinusoidal type.The spike intertimes density is thus given in closed form, whereas the mean and the variance are obtained and shown for some suitable instances by means of numerical computations.

2.
A spike train probability model.A customary believe in neuroscience is based on the hypothesis that the neural coding adopted by the brain to handle information is based on the neuronal spike (the number of spikes in the time unit), or on the temporal occurrence of spikes (the sequence of spikes).Within both paradigms, since spikes have very short duration, point processes or counting processes are commonly used as probability models of spike trains.
The occurrence of neuronal spikes is often described by the inhomogeneous Poisson process.It is a continuous-time stochastic process {N (t); t ≥ 0}, with state space the set of non-negative integers, where N (t) denotes the number of spikes of a single neural unit occurring in [0, t] (see, for instance, Burkitt [7] and [8] for comprehensive reviews of the integrate-and-fire neuron model, where the stochastic synaptic inputs are described as a temporally homogeneous and inhomogeneous Poisson process).The intensity function of the inhomogeneous Poisson process is defined as follows: It represents the intensity of occurrence of a spike at time t in a single neural unit.Various choices of λ(t) have been proposed in the past.In the simplest case it is constant in t, this leading to a homogeneous Poisson process.Function ( 1) is useful to describe various quantities of interest.For instance, let τ j be the j-th spike time (j = 1, 2, . ..) of a single unit; denote by Λ(t) = t 0 λ(s)ds the mean function of N (t), and assume that Λ(t) < +∞ for any finite t ≥ 0, with lim t→+∞ Λ(t) = +∞; then the probability density function of τ j is: 1) is based on the assumption that the following conditional intensity function exists:

A customary extension of definition (
where 0 < τ 1 < τ 2 < . . .< τ N (t) is the sequence of spike times occurring in [0, t].Function (2) thus describes the intensity of occurrence of a new spike at time t conditional on the spike times occurred in [0, t].
In order to describe specific properties of spike trains, such as the neuronal refractory period, various authors follow an approach based on the assumption that λ(t | τ 1 , τ 2 , . . ., τ N (t) ) is expressed as product of two suitable functions (see, for instance, Berry and Meister [6], Johnson and Swami [22], Kass and Ventura [23], Miller [25]), i.e.
In Eq. ( 3), s(•) and r(•) are suitable non-negative functions, s being known as the free firing rate function and r as the recovery function.Recently, Chan and Loh [9] investigated this model with reference to template matching of multiple spike trains, and to maximum likelihood estimators of the free firing rate and recovery functions.
We notice that model ( 3) is Markovian because the conditional intensity of spikes is assumed to depend only on the present time t and on the duration t − τ N (t) since the last spike.

3.
A model for interacting neural units.We aim to study the model described by Eq. ( 3) in a more general case that includes interaction among units.Indeed, we consider a network of d excitatory neural units, say . ., N d (t) be counting processes, where N i (t) describes the number of spikes of unit U i in [0, t], for 1 ≤ i ≤ d.Moreover, we denote by τ i,k the k-th spike time, k = 1, 2, . .., of unit U i , for 1 ≤ i ≤ d.The sequence of overall spike times of the network occurring in [0, t] will be denoted as where the counting process counts the total number of spikes occurring in [0, t].For k = 1, 2, . . .and 1 ≤ i ≤ d, we set In analogy with the model expressed by (3), the conditional intensity function of the unit U i , for 1 ≤ i ≤ d, is assumed to have the following form, for t ≥ 0: where G t collects all information related to the activity up to time t, i.e.
Function s(t) is non-negative and such that +∞ τ s(t) dt = +∞ for any τ > 0. As for model (3), it is named free firing rate function, since it describes the spiking intensity of the network's units due to external inputs, and in absence of firing activity.From Eq. ( 7) we note that if means that the occurrence of the first spike is uniform over the d units.Moreover, we have: In the general setting s(t) is a time-varying function, which allows for the description of stimuli with varying amplitudes such as modulated inputs.Again, function r i (•; •) is non-negative, and is called the recovery function of unit U i .Its main role is the inclusion in the model of the refractory period of U i , and also of the effect of the spiking activity of the other network units.
Remark 1. Due to Eq. ( 7) the intensity function of N i (t) does not depend on i when N (t) = 0, whereas it depends on the counting process (5) through τ •N (t) and Z N (t) , when N (t) ≥ 1.The firing activity of the i-th neural unit is thus governed by the last spiking time, τ •N (t) , and by the last spiking unit of the network, Z N (t) .Moreover, N 1 (t), N 2 (t), . . ., N d (t) are conditionally independent processes, in the sense that the distribution of each of such counting processes depends on the remaining d − 1 processes only through the sum (5).
From now on we suppose that the recovery function appearing in the right-handside of ( 7) is given by: for all 1 ≤ i ≤ d and 1 ≤ j ≤ d, where: (i) coefficients c i,j are such that (ii) u(t) is a non-negative continuous function, decreasing for all t ∈ [0, +∞), with u(0) = 1 and lim t→+∞ u(t) = 0.
We point out that the above assumptions concerning Eq. ( 9) yield the following features of the model: • Coefficients c i,j measure the strength of the spiking activity of U j on the network units.Conditioning on Z N (t) = j, thus being U j the last spiking unit before t, we have: (a) If i = j then c j,j = −1; this describes the auto-inhibition due to a neuron spike, i.e. the effect of the refractory period.(b) If i = j then the coefficients c i,j are strictly positive, this yielding a full interaction (of excitatory type) among the network's units.In some sense, they give a measure of the synaptic strength from U j (the presynaptic neuron) to U i (the postsynaptic neuron).
• Function u(•) describes the effect over time of the spiking activity on the network units.When t is close to last spiking time τ •N (t) , the last spiking neuron, U j , is less likely to process the stimuli arriving according to the free firing rate function s(•), since this being in agreement with the effect of the refractory period.Moreover, for all t and 1 ≤ j ≤ d we have • All other units U i , i = j, receive a stimulus from the last spiking neuron U j , the strength of the stimulus being regulated by c i,j .In this case, for all t and i = j, it is • The effect of the last spike tends to vanish as time proceeds; indeed, for all 1 ≤ i ≤ d and 1 Note that an accurate choice of the recovery function r i (t − τ •N (t) ; Z N (t) = j) should treat the cases i = j and i = j as different since they arise from distinct physical situations.When i = j we deal with the auto-inhibition of a neuron due to spikes, and then the modeling of the refractory period should include time-delay effects in function u(•).On the contrary, when i = j we deal with the interaction between different neurons, and thus such delay is not required.Nevertheless, in order to make the model mathematically treatable, the cases i = j and i = j have been unified in the right-hand-side of Eq. ( 9).On the other hand the condition (11) implies that, within the present model, spikes close in time from the same neuron are very unlikely.
Remark 2. In order to assess the plausibility of the above assumptions in a model of neural spike trains, we point out that the mean interspike intervals (of the superposed spike trains) should be larger than the characteristic time scale of the recovery function.In a broad sense, the model is physiologically plausible when the recovery function (9) decreases rapidly as t increases.
Recalling Remark 1, the first spike occurs according to the free firing rate s(t) (see Eq. ( 8)), so that τ •1 has distribution function Moreover, the probability that the first spike is generated by unit U i is uniform, since P (Z 6) and (7).We now introduce the random vectors X where, in agreement with (7), is a non-negative random variable having hazard rate s(t) r i (t−τ •N (t) ; Z N (t) = j).Assuming that the k-th spike of the network was generated by unit U j at time τ describes the time length between τ •k and the next spike, conditional on the event that the latter spike is generated by unit U i , 1 ≤ i ≤ d.From the above assumptions it follows that the spiking process is regenerative, in the sense that the distribution of X does not depend on k.Hence, we shall write X (τ ) i,j when it is not necessary to specify the index k.Moreover, as soon as a spike occurs, the firing activity restarts afresh according to the scheme described by Eqs. ( 7) and (9).We notice that the components of vector (12) are not observable, whereas the following random variables are observable: for 1 ≤ j ≤ d.Clearly, T (τ ) j denotes the time length between a spike discharged at time τ by unit U j and the next spike produced in the network, the unit producing such spike being described by δ (τ ) j .On the ground of Eqs. ( 7) and ( 9), the distribution function of X (τ ) i,j is given by In the following we shall denote by the probability that a spike of unit U j , occured at time τ , is followed by a spike of the same unit.
We remark that the above framework can be viewed as referring to the classical "competing risks model".The latter deals with failure times subject to multiple causes of failure, and deserves interest in various fields such as survival analysis and reliability theory.In the present case the roles of failures and of failure causes are played, respectively, by the observed spikes and by the firing network units.General properties of the competing risks model can be found for instance in Crowder [10], whereas recent results on such model related to ageing notions and shock models are given in Di Crescenzo and Longobardi [16] and [17], respectively.4. Analysis of the model.Aiming to give a deeper description of the model introduced in the previous section, we first consider the simple case where the network is composed of d = 2 units.Due to (10), for d = 2 and i, j = 1, 2 we have so that Eq. ( 9) becomes for i, j = 1, 2. Recalling ( 12) and ( 13), now we deal with the random vectors whose components are not observable.On the contrary, the random variables T (τ ) j and δ (τ ) j , j = 1, 2, defined in (13), are observable.Since the matrix ||c i,j || in this case is symmetric (cf.( 16)), we can introduce two random variables X (τ ) − and X (τ ) + , by renaming the components of the random vector (17) as follows: 1 Hence, from the given assumption it is not hard to prove that X − and X + are non-negative independent random variables, where X A sample of activity of a network with d = 2 units.time length between a spike occurring at time τ and the next spike, conditional on the event that the latter spike is due to the same unit (resp., the other unit).An example of activity of a network with d = 2 units is shown in Figure 1 where, for instance, X Recalling (14), the complementary distribution functions and the probability density functions of variables ( 18) can be expressed respectively as for t ≥ 0.Moreover, due to (15), and since X − and X + are independent, when d = 2 the probability that a spike of a generic unit, occured at time τ , is followed by a spike of the same unit is given by We are now able to provide the expressions of ( 20) and of the distribution function of the observable random variable Note that T (τ ) describes the intertime between a spike occurring at time τ and the subsequent spike.A relevant role is played by the free firing rate function s(•) and by the auxiliary function u(•) appearing in the recovery function (9).
Since u(•) is a non-negative function, from (22) we have q (τ ) ≤ 1/2.Thus it is more likely that consecutive spikes are displayed by different units rather than the same unit.The function φ τ (t), defined in Eq. ( 24), is named cumulative firing rate.
The analysis of the model in the case of a network of d units can be performed by taking into account that, similarly to (20), the probability ( 15) is given by where f i,j (t | τ ) and F i,j (t | τ ) denote respectively the probability density and the complementary distribution function of X (τ ) i,j , for i, j = 1, 2, . . ., d. Due to (10), the terms in the right-hand-side of (25) do not depend on j, and thus we are now able to give the following extension of Proposition 1.
Hereafter, in Sections 4.1 and 4.2 we consider two special cases in which s(t) is constant and of sinusoidal type.

4.1.
Constant free firing rate.In this section we discuss the homogeneous case, in which the external inputs arrive to the network's units according to a constant intensity.We thus assume that the free firing rate is constant, so that s(t) = λ for all t ≥ 0, (28) with λ > 0. We point out that in this case the distribution functions given in (14) do not depend on τ , and thus can be expressed as follows: where If r = 1 the following expression of q holds: where c is defined in (32), and Γ(•, •) is the upper incomplete gamma function.
For both cases treated above, Figure 2 shows some plots of q as function of c, with various choices of r, and for d = 2.
We point out that Proposition 3 states that the interspike intervals described by T are exponentially distributed.This is significantly different from the distribution functions specified in (29).

4.2.
Sinusoidal free firing rate.Several papers on neuronal activity focus on modulated stimuli described by periodic inputs.For instance we recall Tateno et al. [33], where the problem of finding the period of the oscillation in an oscillator driven by a period input is studied by means of a first-passage-time approach, and Yoshino et al. [34], where the effect of periodic pulse trains on oscillatory regimes neuronal membranes is investigated.More recent researches studied the behaviour of the leaky integrate-and-fire model driven by a sinusoidal current or slowly fluctuating signal (see, for instance, Barbi et al. [2], Picchini et al. [26]).
Aiming to include the presence of periodic external stimuli in model (7), in this section we consider the inhomogeneous case in which the time-varying free firing rate is given by where |A| ≤ λ and P > 0. Hence, due to (27) the density of the spike intertimes T (τ ) for a network of d units is where, due to (24), the cumulative firing rate is Figure 3 displays some plots of density (34) for some choices of the involved parameters.It shows that the multimodality of such density reflects the periodicity of the free firing rate (33).Figure 4 gives the mean M = E[T (τ ) ] and the variance In this case a closed-form expression of probability q (τ ) seems not available.However, it can be numerically evaluated by making use of Proposition 1. See Figure 5 for some plots of q (τ ) when u(t) = e −t , t ≥ 0. In particular, the oscillating behaviour of q (τ ) with respect to τ is evident for large values of A (see the right panel of Figure 5).product of the free firing rate function and a suitable recovery function.We have proposed an extension dealing with a neural network composed of d excitatory units, in which the recovery function of each unit depends both on the time elapsed since the last spike and on the last spiking unit.Our approach, which is somewhat related to the competing risks model, leads to the general form of the interspike distribution and of the probability of consecutive spikes from the same unit.Explicit results have been found when the free firing rate function is constant.We also considered the case when the free firing rate is sinusoidal, for which the density, the mean and the variance of the spike intertimes is investigated by means of numerical evaluations.In both cases we studied the probability that a spike of a generic unit, occured at a fixed time, is followed by a spike of the same unit.

5 .Figure 5 .
Figure 5. Plots of q (τ ) in the sinusoidal free firing rate case as a function of A (left panel) and of τ (right panel), for u(t) = e −t , t ≥ 0, with d = 2, λ = 1 and P = 2.