{'fraction_non_alphanumeric': 0.05285921411544983, 'fraction_numerical': 0.008158651365099648, 'mean_word_length': 4.061442786069652, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On stationary distributions of stochastic neural networks

The paper deals with non-linear Poisson neuron network models with bounded memory dynamics, that can include both Hebbian learning mechanisms and refractory periods. The state of a network is described by the times elapsed since its neurons fired within the post-synaptic transfer kernel memory span, and the current strengths of synaptic connections, the state spaces of our models being hierarchies of finite-dimensional components. We establish ergodicity of the stochastic processes describing the behaviour of the networks and prove the existence of continuously differentiable stationary distribution densities (with respect to the Lebesgue measures of corresponding dimensionality) on the components of the state space and find upper bounds for them. For the density components, we derive a system of differential equations that can be solved in a few simplest cases only. Approaches to approximate computation of the stationary density are discussed. One is to reduce the dimensionality of the problem by modifying the network so that each neuron cannot fire if the number of spikes it emitted within the post-synaptic transfer kernel memory span reaches a given threshold. We show that the stationary distribution of this `truncated' network converges to that of the unrestricted one as the threshold increases, and that the convergence is at a super-exponential rate. A complementary approach uses discrete Markov chain approximations to the network process. We derive linear systems for the stationary distributions of these Markov chains and prove that these distributions converge weakly to the stationary laws for the original processes.


Introduction
Neurons are electrically excitable cells whose main function is to process and transmit information. They connect to each other to form neural networks that constitute core components of the nervous system, and so building and studying mathematical models of such networks is of key interest. To justify the modelling approach used in this paper (as described in detail in Section 2), we will briefly describe the mechanism enabling neurons to communicate with each other.
The anatomy of a neuron involves three distinct parts with different electrical activity functions: dendrites that form a tree and contain post-synaptic receptors (inputs), the cell body (soma) that integrates the input currents coming from the dendrites, and a long-limbed axon that terminates with pre-synaptic buttons (outputs). A typical feature of the neuronal electrical activity is the propagation of membrane depolarisation. The membrane of a resting neuron is polarised. Brief high-amplitude depolarisations that propagate from the soma along the axon are called action potentials (or spikes) and have a characteristic shape. When a spike reaches an axonal termination that "connects" to a post-synaptic neuron, neurotransmitters are released into the extra cellular space and excite receptors on the post-synaptic neuron (usually on dendrites). This generates a local variation of the membrane potential in that neuron, which propagates towards the soma. The soma can be seen as a spatio-temporal integrator of these post-synaptic potentials (PSP) to generate an output spike. The soma potential often remains close to the resting value for a few milliseconds after firing an action potential, which is referred to as the refractory period (through which the neuron cannot fire again). These basic elements of the neuronal information processing actually depend upon many different mechanisms at the molecular level, such as ionic concentrations, density of ion channels, axonal myelination, and types of neurotransmitters (for a review, we refer the reader to [1]).
There exists extensive literature on mathematical modelling of both individual neurons and neural networks. We refer the reader interested in neurophysiological principles of neuron and brain operation to [1,27] and more advanced expositions in [31,32] of the circuitry of the brain. A detailed (but accessible and rather non-technical) discussion of brain networks, covering structural, functional, and effective connectivity and their respective dynamics, is presented in [33] (the book also contains extensive bibliography of the relevant research work in the complex network theory). A detailed overview of the computational modelling of nervous systems from the molecular and cellular level, including mathematical modelling of adaptation and learning, is given in [13]. Monograph [23] is a systematic study of the relationship of electrophysiology, nonlinear dynamics, and computational properties of neurons. One can also mention here [25,15] and refer to [17] for a recent review of the literature in the area.
As the duration of an action potential is relatively short (usually less than 1 ms for sodium-based action potentials), for modelling purposes spikes are often considered to be instantaneous. Hence a neuron can be modelled using a point process whose intensity depends on the past activity of the neuron, its incoming synaptic stimulation or other mechanisms. The use of point processes, such as Hawkes processes, for modelling the spiking activity of neurons dates back to papers [9,10,12] and made it possible to study analytically the neuronal response to various input stimulations (for a review, see [15]). Applications of particular neural non-linear point process models to real data can be found in [28,29,30,34,35].
In Section 2 we present descriptions of two network models we are dealing with in this paper. Both are "assembled" of non-linear Poisson neurons that can be viewed as extended versions of self-exciting Hawkes point processes, the difference between the two models being that the former has constant strengths of synaptic connections between neurons, whereas in the latter, to model Hebbian learning, we allow the strengths to change depending on the order in which the connected neurons are firing. Instead of using the formalism of point processes (as e.g. in [6,8]), we choose an alternative description in terms of multivariate Markov processes whose states represent networks' spiking histories, with state spaces being products of hierarchies of simplices. This approach proves to be rather convenient and allows one to demonstrate ergodicity of the network processes under rather general conditions (Section 3) and, moreover, to study the stationary distributions thereof. We show that the stationary distribution of the network has a smooth density with respect to a natural measure on the state space and give upper bounds for the components of that density on different components of the state space of the process.
In Section 4 we discuss a way to reduce the dimensionality of the model and approximate its stationary distribution with more tractable objects. The approach is based on "truncating" the original process by "forbidding" neurons to fire once they have fired a given number n of spikes recently (within the "memory window" of the neuron). Thus modified process will still be Markovian and ergodic, its stationary distribution confined to a space of lower dimensionality and approximating that of the original process at a super-exponential rate in n. In fact, such dynamics do make physical sense when the existence of the refractory period is taken into account, but one can further simplify the model by choosing an even lower threshold n.
Section 5 deals with the problem of computing the stationary distributions of the networks. We derive systems of differential equations for the stationary distributions (unfortunately, they seem to be tractable in the simplest cases only, that are discussed as examples). Moreover, we prove that the stationary distribution of our network process can be approximated by those of discrete Markov chains constructed as discretised (in both time and space) versions of the process. Computing the stationary distributions for the chains is more feasible, as it only requires solving systems of (a large number of) linear algebraic equations.

Network dynamics and its description by Markov processes
We consider a model neural network consisting of N neurons and M external sources. Both external sources and neurons can fire spikes, which are assumed to be generated by a random mechanism. External sources are assumed to fire according to independent Poisson processes with constant ratesρ k , k ∈ {1, 2, . . . , M } (all the quantities related to external sources will be labelled with hats, and all the processes in the paper will be assumed to be right-continuous), whereas a neuron's instantaneous firing rate is determined by the value of the activation function of the so-called synaptic influx. For neuron i, the latter is the sum of all PSPs generated by spikes arriving to the ith neuron's synapses from external sources and other neurons, and also background activity. More precisely, assuming that { T k,n } n∈Z are times at which external source k ∈ {1, . . . , M } fired, {T j,n } n∈Z are times at which neuron j ∈ {1, . . . , N } fired, andǫ ik (t) and ǫ ij (t) are post-synaptic response kernel functions describing the effects on neuron i potential from accepting spikes through synapses connecting source k to neuron i and neuron j to neuron i, respectively, the total time t synaptic influx for neuron i is given by where v i = const represents the background activity for neuron i, the synaptic weights W ik (t) and W ij (t) can be positive (excitatory synapse) or negative (inhibitory synapse) and, to reflect brain plasticity (e.g. to model Hebbian learning), they can depend on time as well. If, at time t, there is no synaptic connection of external source k to neuron i, we simply have W ik (t) = 0, and likewise for network neurons' connections.
The kernelsǫ ik (t) ≥ 0 and ǫ ij (t) ≥ 0 are assumed to vanish outside a compact interval: for any k ≤ M and i, j ≤ N, which ensures causality and also means that the direct effect of any given spike on a neuron completely disappears within a finite time Θ (for real life neurons, the order of magnitude of Θ is 10 2 ms).
In a previous series of papers by one of the authors ( [18,20]), the case where all kernels were identical to some function ǫ, but incorporating individual synaptic delaysd ik and d ij :ǫ was considered. The delays account for both the axonal propagation of action potential up to the synaptic site, and for the diffusion time of the neurotransmitters in the synaptic cleft. In the present paper, we allow each synapse to have individual properties.
The effect of the synaptic influx on the behaviour of neuron i is expressed via an activation function ς i (·), which is assumed to be continuous non-decreasing (and usually "S-shaped"), with Namely, denoting by T i (t) the time of the last spike fired by neuron i prior to time t, by F t the σ-algebra generated by the evolution of our system up to time t and setting ∆(t) := (t, t + ∆) for ∆ > 0, we have, as ∆ → 0, where we used a left-continuous function r(·) ∈ [0, 1] to model the existence of the so-called absolute refractory period, i.e. the time period during which a just fired neuron is unable to fire again. One can take e.g.
the indicator function of the complement of the interval (0, δ AR ] (for real life neurons, δ AR is about 1 ms). Whatever the shape of r, we always assume that r(s) = 1 for s ≥ Θ.
In addition to (1), we assume that P more than one neuron fires during ∆(t) which basically means that, given the past history F t , the instantaneous firing of different neurons is driven by independent random mechanisms. Note that the widely studied classical Hawkes process (see [22,5,7,8,16]) corresponds to the identity activation function ς i in (1), and that the positivity of ς i means that neurons can fire spikes in the absence of any external stimulation. The use of non-linear bounded functions ς i is motivated by the experimentally observed saturation of the neuronal firing rate when its excitation increase. Observe also that the temporal spread of the synaptic responses (modelled by ǫ ij ) induces specific temporal correlations between the neuronal spike trains, which can be evaluated for the case of linear activation functions ς i [22,19,20].
We will consider two types of models that differ in their assumptions concerning synaptic weights: Model I assumes that all the synaptic weights are constant: W ik (t) ≡ W ik = const and W ij (t) ≡ W ij = const for any k ≤ M and i, j ≤ N .
Model II assumes a Hebbian learning mechanism in the form of spike-timing dependent plasticity (STDP) : if, within a short enough time interval, there are spikes at both pre-synaptic and post-synaptic sides of a connection, this can change the weight of the connection. The weight increases if the post-synaptic spike follows the pre-synaptic one (reinforcement of the synapse), and decreases otherwise (depression of the synapse).
For simplicity we assume that, for each of the connections, the synaptic weight can assume finitely many values: for a common finite L, W ik (t) ∈ G ik := {ĝ ik (1) ≤ĝ ik (2) ≤ · · · ≤ĝ ik (L)}, W ij (t) ∈ G ij := {g ij (1) ≤ g ij (2) ≤ · · · ≤ g ij (L)}, and the following discrete approximation of the STDP mechanisms discussed e.g. in [11]. For any i, j ∈ {1, . . . , N } and m ∈ {1, . . . , L}, we have a collection of points In real life situations, the length δ LW := max{u ij (m, m), |u ij (m, m + 1)|} of the "learning window" is about 10 2 ms. We assume that δ LW < Θ. Now suppose that, for a given time t, one has W ij (t−) = g ij (m) and either t = T i (t) or t = T j (t) (i.e. one of the neurons i, j fired at time t). Then we put , . . . , L}, while otherwise the value of the weight remains unchanged. A similar rule applies to the weights W ik . Note that the above mechanism allows one to model the emergence of new synaptic connections as well. Altogether, our network models provide a certain degree of biologically realism together with a mathematical framework that allows a tractable analysis.
Having described the rules governing of our neural network, we will now present a Markov process model for it. Observe that, at time t, the knowledge of all the current synaptic weights and the times of all the spikes fired in the network within the time interval (t − Θ, t] is all the information from the past and present that one needs to uniquely specify the probability distribution of the future evolution of the system. Therefore, to obtain a Markovian description of the network, we denote bŷ ν k (t) the number of spikes fired by external source k in the time window (t − Θ, t], k ≤ M . Ifν k (t) = 0, then we say that source k was at the state , which is the time till the last spike fired by k prior to the "present" time t disappears from the moving window (s − Θ, s], s ≥ t, and then we denote by X k,2 (t) : , the time till the second last spike fired by k prior to time t disappears from the moving window (s − Θ, s], and so on, so that in this case we always have s. (as having two spikes at exactly the same time is a zero probability event).
Likewise, the state of neuron i is described by the vector being the number of spikes fired by i during (t − Θ, t]. Now the complete history of spikes within the time window (t − Θ, t] is described by the vector For Model I, this vector will completely specify the state of the network. The state space for the process Z will be taken to be is the union of simplices is a face of the (n + 1)dimensional one, E (n+1) , n ≥ 0. We will endow S with the product σ-algebra For Model II, we need to specify in addition the state of the synaptic connections. This can be done by using the matrices in which to non-existent connections there will correspond zero entries. The new process Z * := ( X; X; W ; W ) will have the state space endowed with the natural product σ-algebra that we will denote by S * .
Model I dynamics. Suppose we are given an initial condition Z(0) ∈ S. Following the earlier description of the dynamics of our network, in the case of Model I the process Z is a piece-wise deterministic (linear) Markov process, which evolves for t > 0 as follows.
(i) Inside time intervals free of jumps, one has, for any k ≤ M, i ≤ N , n > 1, This means that all the non-zero components of the process decay at the unit rate, and when the "first visible in the window" spike of, say, neuron i that occurred at time T i,n "disappears" from the moving time window (t−Θ, t] at time t ′ = T i,n +Θ, the number of positive components of X i drops by one: the instantaneous firing rate for source k isρ k , and for neuron i it is given by (iii) When source k fires (say, at time t ′ = T k,m ), the only change in the state of Z is in the component X k :ν a.s. (as it is impossible to simultaneously "lose" a spike in the time window and acquire a new one), and the new values of the components are: Likewise, a spike fired by neuron i will mean similar changes in the component X i .
It is quite straightforward to write down the generator of the process Z, following the above description. The vector field specifying the dynamics of the process between jumps is piece-wise linear, and it changes its direction when, for one of the components X k ∈ E or X i ∈ E, the respective integral curve running inside E (n) , n > 1, hits the face E (n−1) of that simplex and then continues inside that lower dimensional simplex. The domain of the generator will consist of all bounded functions S → R that are path-continuous and differentiable for that vector field (cf. [24]).
Model II dynamics. The trajectories of Z * will also be piece-wise deterministic (linear), with its first two components following (3) and the last two remaining unchanged between successive jumps. Jumps occur at the times when either external sources or neurons fire spikes, and, given the current state of the process is z * = (x; x;ŵ; w) withŵ = (ŵ ik ) and w = (w ij ), the instantaneous firing intensities areρ k for source k and, instead of (4), for neuron i. When a spike is fired, the change in the components X and X is exactly the same as for Model I (see part (iii) of the description of its dynamics above), whereas the synaptic weight W ik can change when the spike was fired either by source k or by neuron i. As δ LW < Θ, the state of Z * just prior to the spike completely specifies to what value the synaptic weight should change, according to the learning rules listed in the description of Model II. Similarly for the weights W ij that can change when the spike is fired by either of the neurons i and j.
Thus we see that, in the case of Model II, Z * is a well-defined Markov continuous time process completely describing the dynamics of the system. Its generator will differ from the one for Z by the presence of terms related to jumps in the synaptic weights' values.

Ergodicity and the properties of stationary distributions
In this section, we establish strong ergodicity of the Markov process Z * describing the dynamics of Model II. As Model I is a special case of the latter, this means that the process Z is also ergodic. The latter fact is actually an immediate consequence of Theorem 5 in [6] on stability of multivariate point processes with bounded memory dynamics (see also Theorem 6 in [6] for stability of a nonlinear multivariate Hawkes process with PSP transfer kernels having unbounded supports, and [26]). However, even in the case of Model I, our Markov process framework allows us to come up with much shorter and simpler proof of stability.
Theorem 1. Under the stated assumptions for Model II, the process Z * is strongly ergodic: it has a unique stationary distribution π on (S * , S * ) such that Moreover, the convergence is exponentially fast.
It would be most interesting to know the properties of the stationary distribution π. One basic fact that we can easily establish is that the distribution π S of the first two components of ( X(∞); X(∞); W (∞), W (∞)) ∼ π on (S, S ) has a density w.r.t. some naturally chosen measure. Moreover, we can obtain upper bounds for the density.
More precisely, that natural measure on (S, S ) is taken to be the product measure µ M +N , where µ n being the n-dimensional Lebesgue measure and B n being Borel subsets of the respective n-dimensional simplices that can be identified with We use the convention that µ 0 is just the unit mass at 0.
is the product of the respective E 0 -sets. To simplify the formulation of the next theorem, we will slightly abuse notation by identifying the sets E (m;n) with E (m;n) 0 and so considering the latter as the components of the state space S (so that the components of µ are actually given on finite-dimensional spaces).
where Σρ := kρ k , Λ := max{max kρk , max i ς i }, Σ m := k≤M m k , and Σ n := i≤N n i . The density function ψ m,n is continuously differentiable in the interior of E (m,n) 0 and has finite limits on its boundary. Remark 1. As will be easily seen from the proof of Theorem 2, if we assume that the function r has form (2) (and so is not continuously differentiable), the assertion of the theorem will remain true with the only amendment that the density components ψ m,n will be continuously differentiable inside their supports in the spaces of the respective dimensionalities.
Proof of Theorem 1. It is obvious that Z * is aperiodic and stochastically continuous. So it suffices to show that the Markov "skeleton" chain Y = {Y n := Z * (Θn), n = 0, 1, 2, . . .}, has a recurrent state whose first hitting time distribution tail decays exponentially fast uniformly in the chain's initial state Z * (0) (see e.g. Theorem 18.1 in [2]).
First we will use a standard argument to show that the tail of τ := inf{n > 0 : Z(Θn) = (0; 0)} (i.e. the first value n such that there were no spikes in (Θ(n − 1), Θn]) admits such a bound. Indeed, setting for convenience P z * (·) := P(·|Z * (0) = z * ), we have, for any z * ∈ S * and t > 0, where Σ ς := i ς i , and (· · · ) represents the argument of R * i along the trajectory of Z * on [0, t] that started at z * and experienced no jumps. Now setting A n := {no spikes in (Θ(n − 1), Θn]} we obtain, using recursively the Markov property and bound (8), that, for n ≥ 0, Next we observe that ( W (Θn); W (Θn)) is clearly an indecomposable aperiodic finite Markov chain, and hence it is ergodic. Take any fixed state (ŵ ′ ; w ′ ) of this chain; as it is well known, for any initial condition, the first hitting time of that state has an exponentially fast decaying distribution tail. Hence it is obvious that the state (0; 0;ŵ ′ ; w ′ ) ∈ S * will be positive recurrent for the chain Y , and that the tail of the first hitting time of that state by Y will admit a geometrically fast vanishing upper bound uniform in the initial condition of the state. The theorem is proved.
Proof of Theorem 2. First we will establish existence of density for transition probabilities, and then infer the desired result from that fact.
Suppose our process started at point Z * (0) = v * and, at time Θ, was at a Clearly, the function p(v * , z) is continuously differentiable in z in the interior of E (m,n) , has finite limits on its boundary, and admits an upper bound of the form Next, in view of (10), for any B ∈ S , we can use Fubini's theorem to write This means that π S does have density ψ w.r.t. µ, and (11) implies that ψ admits the desired upper bound.
That ψ is continuously differentiable in the relative interiors of the components of its supporting space follows from representation (10), the last relation in (12) and the assumption that ς i and the kernel functions ǫ ij are all continuously differentiable, the ǫ's vanishing outside [0, Θ]. Theorem 2 is proved.

Approximation of π by finite-dimensional distributions
In this section we will be dealing with the simpler Model I. Even for that model, the state space is an infinite hierarchy of multidimensional simplices, so working with non-trivial distributions on it is not easy. The natural question in such a situation is whether one can find an appropriate approximation to the distribution in question, together with an approximation error bound. For our model, a tempting approach to finding such approximations is to consider "truncated" processes Z n in which none of the neurons is "allowed" to fire more than the fixed number n ≥ 1 times within any given time interval of length Θ. In fact, if the model assumes existence of absolute refractory periods of positive length by stipulating, say, that (2) holds true, then that condition will automatically be satisfied (note, however, that one can still apply truncation with n < Θ/δ AR to reduce dimensionality).
The only difference in the dynamics of the process Z n compared to those of Z is that neurons' firing intensities will now be given by (cf. (4)). It is obvious that Z n will also be an ergodic Markov process. Denote its stationary distribution on (S, S ) by π n , while for the stationary distribution of Z (on the same measurable space) we will re-use notation π.
Proof. We will use coupling. Assume that Π is a Poisson random field of unit intensity on R + × R, given on some probability space, and construct a process {(Z(t), Z n (t))} t≥0 with the state space S × S, whose components follow the original and "truncated" dynamics, respectively, start at a common state Z(0) = Z n (0) ∈ S, and are driven by the field Π via the following simple mechanism. Introduce intervals Clearly, (Z, Z n ) is a well-defined Markov process, and its components follow the desired dynamics. Note also that the process will be ergodic, like each of its components (this is obvious e.g. from Theorem 1). Now denote by (Z(∞), Z n (∞)) a random element of S × S whose distribution coincides with the stationary distribution of (Z, Z n ). Using the standard argument, it is easily seen that π(B) − π n (B) ≤ P(Z(∞) = Z n (∞)) =: P n , B ∈ S .
To bound P n , denote by |T K |, K ∈ N, the total length of the set and observe that, from the ergodicity of (Z, Z n ), one has As we are interested in the stationary distribution, we can assume w.l.o.g. that the common starting point of Z and Z n has no components x i in E (m) , m ≥ n. Then the trajectories Z(t) and Z n (t), having originated at the same point, will coincide with each other till the time T ′ when one of the values X i (t), i = 1, . . . , N, enters E (n) . Then the respective neuron i will stay silent in Z n at least till the time when the number of spikes produced by i and "visible" in the time window (t − Θ, t] drops below n, while in Z the respective neuron will still be able to fire. Thus, past that time point T ′ , the trajectories Z(t) and Z n (t) can diverge. They will have to meet again, though, and the latest that will occur is at the end of the next "silent interval" of length Θ, which, in its turn, occurs no later than at the time inf t > T ′ : To make use of the above argument to obtain an upper bound for P n , introduce two random sequences, {κ m } m≥0 and {γ m } m≥0 , as follows. Letting for brevity θ j := Θj, set Now if Z(θ j−1 ) = Z n (θ j−1 ) but, for some t ∈ (θ j−1 , θ j ], one has Z(t) = Z n (t), then, for some i ≤ N, at least one of the following two relations must hold: if none of the two holds then, in any time interval of length Θ within (θ j−2 , θ j ], neuron i will have fewer than n spikes). Thus the values κ m "mark" time intervals where Z and Z n may split, whereas γ m "mark" those intervals following κ m where Z and Z n must merge (provided that they have split indeed). Set Clearly, for t > 0 and an arbitrary fixed ε > 0, First we will bound the last term. Observe that are i.i.d. Bernoulli random variables with success probability Since V i,j− has the Poisson distribution with parameter λ i := Θς i /2, one can use Taylor's formula for the exponential series (with remainder in Lagrange form) and then Stirling's formula to write Therefore Assuming that p n < 1 (otherwise the bound in the theorem will become trivial), and that δ := ε − p n ≡ ε − Eχ j > 0, we obtain by virtue of Theorem 10 from Chapter 5 of [3]. Now we will turn to the first term on the RHS of (15). From the definitions of our random variables, it is obvious that |T K | ≤ Θ m≤H K (γ m − κ m ), and so From the strong Markov property it follows that η m := γ m −κ m are i.i.d. geometric random variables, with P( Clearly, Eη 1 = 1/q and ϕ(a) := Ee aη 1 < ∞ for a < − ln(1 − q), Therefore, assuming w.l.o.g. that εK is integer, we have by the exponential Chebyshev's inequality that One can see from (17) that, choosing t = t K := ΘεK(1 + h)/q for an arbitrary fixed h > 0, we will have, for small enough a, the bound Q ≤ e −εcK for some c = c(a, h).
From here, (15) and (16) we obtain the bound Clearly, K P(|T K | > t K ) < ∞, and so, by Borel-Cantelli lemma, with probability one we have |T K | ≤ t K for all large enough K. Therefore, As this holds for any ε > p n and h > 0, we conclude that P n ≤ p n /q, which completes the proof of the theorem.

Computing the stationary distribution
It is not difficult to derive differential equations (and boundary conditions for them) that the components ψ m,n of the stationary density of Z will satisfy in the case of Model I. They may be derived from the general relation where A is the infinitesimal generator of the process, f a function from a suitable subset of the domain of A, and, as before, Z(∞) ∼ π. It may be easier, however, to obtain them via a direct argument, making use of our Theorem 2 (of which the conditions will be assumed satisfied in this section unless we explicitly state otherwise).
To show how to do that, we will first consider the simple case of a network with one external source and one neuron (with feedback). Suppose that the neuron has an absolute refractory period (so that the state space is actually finitedimensional). For simplicity, we assume throughout this section that Θ = 1 (which clearly does not restrict generality).
In this case, the state space of the process is just E × E, so that each state (x; x) = (x 1 ,x 2 , . . . ,x m , 0, 0, . . . ; x 1 , x 2 , . . . , x n , 0, 0, . . .) (note that here we suppress the unnecessary first subscript indicating the number of the external source or neuron; likewise,ρ will denote hereρ 1 etc.) belongs to one of the components E (m) × E (n) , m, n ≥ 0. The respective density components we will denote by ψ m,n .
Subtracting from both sides the integral of ψ m,n (y + δ) over I z one can then easily verify that the relation implies that the following differential equation must be satisfied: for θ ∈ (0, 1 − max{x 1 , x 1 }), ∂ ∂θ ψ m,n (z + θ) = (ρ + R(z + θ))ψ m,n (z + θ) − ψ m+1,n ((x + θ, 0; x + θ)) − ψ m,n+1 ((x + θ; x + θ, 0)). (20) Of course, the equation will hold along the whole interval formed by the intersection of E , where x * := (x 2 , x 3 , . . . , x n ), caused by a new spike fired by the neuron. Using a probabilistic argument similar to the one above, it is easy to see that the following must hold: A similar equation will hold in the case wherex 1 > x 1 , but then the coefficient of ψ m−1,n on the right hand side of the respective relation will simply beρ. The case where only one of m, n is zero is treated similarly.
Solving equations of the form (20) with boundary conditions (21), complemented by (19) and the condition that m,n ψ m,n d(µ m ⊗ µ n ) = 1, is hardly possible except for the simplest cases. One such case is considered in the following example.
Example 1. Consider the case of a single neuron with feedback and no external sources. Moreover, assume that the firing rate function has the property R(x) = 0 for all x ∈ E (n) , n ≥ 2, so that there cannot be more than two spikes in any given time interval of length Θ = 1 (say, due to the length of the absolutely refractory period exceeding 1/2). Thus the state space of the system is just E (0) × E (1) × E (2) (which we again can and will identify with E 0 ), the density components being ψ n , n = 0, 1, 2 (for n > 2, all ψ n ≡ 0).
Together with (26) this completely specifies the density function ψ (computing ψ 0 is trivial).
Example 2. One can also obtain a closed form solution in the case of a single neuron with feedback and no absolutely refractory period, but under the special assumption that the neuron's PSP kernel is exponential: ǫ(t) = e −αt 1(t ≥ 0) for some α > 0. Of course, this violates the basic assumption that ǫ has a finite support and is smooth, but such a simplified mathematical model could still provide a useful approximation to more realistic ones. It is not hard to see that in this case the dynamics of the system can be described by a univariate Markov process Y (t) := n ǫ(t − T n ) > 0. Assuming without loss of generality that α = 1, one can see that the process Y is driven by the Ornstein-Uhlenbeck type equation t > 0, Z(t) being a pure jump process with unit jumps and instantaneous jump rate γ(Y (t)), where γ(y) := ς(v + W y) and W is the weight of the "self-connection" of our neuron. The infinitesimal generator A of the process Y is clearly its domain's description can be found e.g. in Proposition 4 in [21]. It is not hard to see that Y is ergodic (see e.g. [4]) and so has a unique stationary distribution that we will again denote by π. Substituting (28) into (18) Routine calculation now leads to where we used notation x + := max{x, 0} for the positive part of x. As this equation holds for a large enough class of functions f (see e.g. [21]) and b is continuous and locally bounded, we conclude that π has a locally bounded and continuous density ψ on (0, ∞), satisfying This equation can be solved recursively, on intervals J n := (n, n + 1), n 0. Straightforward calculations show that the stationary density is given by ψ(y) := ϕ n (y), y ∈ J n , n = 0, 1, 2, . . . , where the ϕ n are found recursively as ϕ 0 (y) = ψ(1) exp where, for f defined on J n−1 with a finite f (n−), we denote by Φ n [f ](y), y ∈ J n , the solution φ of the problem The only unknown constant ψ(1) is just the normalizing factor that is to be determined from ∞ 0 ψ(y)dy = 1. At the integer points the density ψ can be defined by continuity.
In the general case, one can only hope to compute approximations to the stationary distribution of the network. One way to do that is to discretise the state space and approximate the differential equations for the density components ψ m,n discussed at the beginning of this section with respective difference equations, and then to solve the latter. However, although the existence of solution to the original system of differential equation follows from Theorem 2, establishing its uniqueness and also the convergence of the solutions to the systems of approximating difference equations presents a challenge. We will follow an alternative approach by first approximating the original stochastic process with a sequence of finite ergodic Markov chains in discrete time, and then proving convergence of their stationary distributions to the desired stationary distribution of Z. In the rest of the section, we will be dealing with our Model I, but one can easily see that analogous results hold for Model II as well. The only reason why we restrict ourselves to Model I here is that the formulation of results for the more general model is much more cumbersome.
The next theorem provides a way for numerical calculation of the stationary distribution π of our original process Z. Endow S with the topology of componentwise convergence and introduce the following notation. For ξ = (ξ 1 , ξ 2 , . . . , ξ q ) ∈ q E, denote by are transition probabilities from those states to ζ.
Moreover, as q → ∞, the distributions q π • q F −1 converge weakly to the stationary distribution π of Z .
Proof. That the chain q Z is ergodic is obvious since it is finite, irreducible and aperiodic. The system of equations that q π is claimed to satisfy is just an explicit form of the usual matrix equation q π = q π q P for stationary probabilities, q P being the transition probabilities matrix of our chain. So we only need to prove the last claim of the theorem. Recall that we used Z n to denote a "truncated version" of the process Z, of which the components X i cannot take values in spaces of dimensionality higher than n (see Section 4). Here we will use the same notation for a similarly "truncated" versions where the components X k are likewise constrained. It is easy to see that the assertions of Theorem 3 remains true in this case as well (with a different value for C).
Denote by q Z n a similarly truncated version of the chain q Z and observe that a complete analog of Theorem 3, with the same bound as in (14) (of which the right-hand side does not depend on q), will hold true for that process as well.
Now fix an arbitrary ε > 0 and choose n so large that the right-hand side of (14) is less than ε. That means that the stationary distributions π and π n of the processes Z and Z n , respectively, will be ε-close in total variation, and the same will apply to the stationary distributions q π and q π n of the processes q Z and q Z n , too, so that sup q>0 sup B∈S π(B) − π n (B) + sup B⊂qS | q π(B) − q π n (B) < 2ε.
This observation implies that it suffices to prove the claim of Theorem 4 for the truncated processes Z n and q Z n that take values in the finite-dimensional space R K , K := n(M + N ),. To simplify notation, we will suppress the superscript n in the next two paragraphs, so that Z will mean there Z n etc.
To prove convergence of the stationary distributions, first assume that Z(0) = q Z(0) = 0 ∈ R K and then observe that, as q → ∞, the distributions of the processes { q Z(⌊qt⌋)} t≥0 weakly converge to that of {Z(t)} t≥0 in the Skorokhod space D R K [0, ∞) (see e.g. Section 5 in Chapter 2 in [14]; by ⌊x⌋ we denote the integral part of x). This can be seen, for instance, from Theorem 2.6 in Chapter 4 in [14] (in fact, the purpose of the "truncation" that we did above as the first step in the proof was to make the state space of the processes locally compact, which is one of the conditions of the theorem). Indeed, extend the domain of the transition operator q T of the chain q Z to all bounded measurable functions f defined on the state space of Z (= Z n ∈ R K ) by setting q T f (z) := E f ( q Z(s + 1))| q Z(s) = h⌊qz⌋ , where ⌊y⌋ denotes the vector whose components are equal to the integral parts of the respective components of y, and let T (t)f (z) := E f (Z(u + t))|Z(u) = z , u, t > 0, be the transition semigroup of Z. Then the conditions of the above-mentioned theorem from [14] will be met provided that we show that, for any continuous function f on R K and any t > 0, one has lim q→∞ sup z | q T ⌊qt⌋ f (z) − T (t)f (z)| = 0.
Because of the semigroup property, it suffices to prove that convergence for t ∈ [0, 1] only, which is not hard to do using representations of the form (10). Indeed, assume that t = 1 (recall that we assumed that Θ = 1 here); the argument in the case t < 1 will be similar, but we will need to integrate over subspaces then, which makes everything even more cumbersome. Partition the component E (m;n) 0 of the domain of integration of T (t)f (z) into cubes of edge length h with vertices on the grid (hZ) Σm+Σn (we use here notation from Theorem 2 and ignore incomplete cubes, i.e. the ones that intersect the "skew" faces of E (m;n) 0 , as their contribution to the integrals will be asymptotically negligible as q → ∞). Fixing one of these cubes, we observe that the probability of the arrival of the chain q Z starting at the point x to the "left bottom" vertex of the cube after ⌊qt⌋ steps will be given by h Σm+Σn times a product approximating the quantity p(x, z) similar to p(v * , z) from (10) (recall that we are dealing with Model I here, so that we do not need the "extended" state variable z * ). Thus q T ⌊qt⌋ f (x) will essentially be an integral sum approximating the integral T (t)f (x), and as the function f is continuous, it is a simple technical exercise to show that (34) holds true.
The last step in the proof is to observe that a bound of the form (9) will hold uniformly in q, n ∈ N for the processes {Z n (t)} t≥0 and { q Z n (⌊qt⌋)} t≥0 as well (resurrecting now the superscripts n ). Therefore there exists a t ε < ∞ such that (recall that we assumed zero initial conditions for all the processes Z n and q Z n ) one has sup q>0 sup B∈B(R K ) P(Z n (t ε ) ∈ B) − π n (B) + |P( q Z n (qt ε ) ∈ B) − q π n (B) < ε. Now the desired assertion follows from (33) and the weak convergence of the distributions of q Z n (qt ε ) to that of Z n (t ε ) as q → ∞. Theorem 4 is proved.