Image Anal Stereol 2008;27:47-52 Original Research Paper MODELLING AND SIMULATION OF A NEUROPHYSIOLOGICAL EXPERIMENT BY SPATIO-TEMPORAL POINT PROCESSES Viktor Benesˇ and Blazˇena Frcalova´ Charles University, Faculty of Mathematics and Physics, Department of Probability and Mathematical Statistics, Sokolovska´ 83, 18675 Praha 8, Czech Republic e-mail: benesv@karlin.mff.cuni.cz, frcalova@karlin.mff.cuni.cz (Accepted February 14, 2008) ABSTRACT We present a stochastic model of an experiment monitoring the spiking activity of a place cell of hippocampus of an experimental animal moving in an arena. Doubly stochastic spatio-temporal point process is used to model and quantify overdispersion. Stochastic intensity is modelled by a Le´vy based random field while the animal path is simplified to a discrete random walk. In a simulation study first a method suggested previously is used. Then it is shown that a solution of the filtering problem yields the desired inference to the random intensity. Two approaches are suggested and the new one based on finite point process density is applied. Using Markov chain Monte Carlo we obtain numerical results from the simulated model. The methodology is discussed. Keywords: filtering, overdispersion, spatio-temporal point process, spike train. INTRODUCTION The spiking activity of the place cells of hippocampus varies with the position of an experimental rat in its arena (Fenton and Muller, 1998). In this process overdispersion takes place, which means the variability is higher than that of a Poisson process. In this case a doubly stochastic Poisson process is a suitable model, see La´nsky´ and Vaillant (2000). In the present paper the aim is to develop stochastic modelling and present new computational methods for the estimation of the spiking intensity characteristics. A simulation study enables to demonstrate the suggested approaches. The stochastic point process theory (Daley and Vere-Jones, 1988) enables to model the neurophysiological experiment in at least two ways. A spatio-temporal point process of events in three-dimensional space considers two dimensions for the arena area and simultaneously the third dimension for the time. A slightly different model is a temporal marked point process with spatial marks. Here times of spikes are marked by the instantenous location of the rat. We will study the spatio-temporal point processes here, see Schoenberg et al. (2002) for more references.Another approach based on the conditional intensity modelling of a temporal point process was developed in Eden et al. (2004). The Poisson process in the d-dimensional Euclidean space Rd is such that in each bounded measurable set B ? Rd the number of events has Poisson distribution with mean lm(B), where lm is the intensity measure. Moreover the numbers of events in disjoint sets are stochastically independent. Overdispersion is modelled by the Cox process X in Rd (Moller and Waagepetersen, 2003) with random driving measure Am such that conditionally X\Am = X is a Poisson process with intensity measure A. For the Cox process, denoting X(B) the random counting measure (number of events in B), we have the overdispersion formula, varX(B) var[Am(B)] ~EXB = E[hm{B)\ ¦ which is a quantity equal to 1 just for the Poisson process (deterministic Am) and bigger else. Cox processes belong to a class of doubly stochastic point processes. In our situation we need a more complex model involving a random movement of the experimental rat. We will develop such a model, then the overdispersion will be first quantified analogously to La´nsky and Vaillant (2000). Further we propose a new approach to the estimation of intensity characteristics from the observed rat path with time monitored spikes. It is based on nonlinear filtering (Fishman and Snyder, 1976) using original techniques based on Markov chain Monte Carlo for point processes. Simulation study of a simple model is added with quantification of overdispersion and intensity characteristics. 47 Benesˇ V et al: Spatio-temporal point processes STOCHASTIC MODELS SPACE-TIME INTENSITY MODELS INTENSITY AND OVERDISPERSION Consider a bounded arena A c R2 and a stochastic process Y = {Yt, t > 0} in A which describes the position y of a rat at time t. Realisation of Y is almost surely a rectifiable curve y = \yt, t > 0}. Further consider a driving random intensity function A(u,t), u G A, t > 0 where the driving measure is Am(D) = JDA(u,t)dudt for a measurable set D c AxR+. Then we define a doubly stochastic Poisson process X with driving measure Am and path Y so that conditionally given A = A and Yt = yt, the number of spikes on y from time 0 to T is Poisson distributed with mean 0l(yt,t)dt. Following La´nsky and Vaillant (2000), the mean number of spikes during [0,T] fired in B C A given intensity A and path Y is E(X(B) | A, Y) = / lB(yt)A(yt,t)dt, 0 and considering further in this subsection all quantities conditioned by Y (omitting this symbol in conditioning) it holds var(X(B)) =E(var(X(B) | A)) +var(E(X(B) | A)) ( / = E(X(B))+var lB(yt)A(yt,t)dt 0 The Fano factor, (Fano, 1947), is then vari J0 lB(ytjA(yt,ty d t J F=1 + - 0T 1B(yt)EL(yt,t)dt (1) which is the event number variance to mean ratio, equal to 1 for homogeneous Poisson process. In experiments A is often homogeneous in time and inhomogeneous in space, i.e., E(A(x,t)) = \ix and var(A(x,t)) = o2, x e A, do not depend on t. Assume that A is divided in l boxes where \ix is piecewise constant and denote /xi intensity mean in i-th box. Consider experimental data in a form (Fenton and Muller, 1998) (tij, nij), i = 1,..., l; j = 1,..., ki, where tij is the duration of j-th stay of the rat in i-th box and nij is the number of spikes during this stay. A natural estimator of the expected intensity /xi is then ßi Hj nij 2" j tij A broad class of spatio-temporal models of the driving intensity function A is that used for Levy driven Cox processes (LCP), cf. Prokesˇova et al. (2006). Here A(x,0= / f{x, tyiy,s)Z[d(y,s)), At{x) (3) where measurable At(x) e R3 is called an ambit set, it determines the part of domain of Z which influences the behaviour of A at the point (x,t). Further f(xt\{y,s) are deterministic functions and Z is a Levy basis, an independently scattered random measure, infinitely divisible, see Barndorff-Nielsen and Schmiegel (2004). It must be chosen so that A is almost surely nonnegative and locally integrable. The class of LCP has useful properties, e.g. if X is LCP then XC is LCP where XC(B) = X(B x [0,T]) is the cumulative spatial point process. A special case of LCP is OUCP (Cox process driven by Ornstein-Uhlenbeck type process), see Lechnerova´ et al. (2008), with driving intensity t A(u,t)= / es-tZ(Bs-t(u)xds) —oo (4) (2) Here Bs(u) = {x gR2; x(x ,u) < ~s}, s < 0, ^ is a metric, Z a Levy basis. The class of OUCP is related to shot noise Cox processes (Moller and Waagepetersen, 2003) and its generating functional is obtained in a closed form. Further we will consider the driving intensity (Eq. 4) of Ornstein-Uhlenbeck type with Poisson basis Z (Barndorff-Nielsen and Schmiegel, 2004). Poisson basis is a Levy basis such that for measurable bounded sets B c R3 the random variable Z{B) has Poisson distribution, the random measure Z has a density which is a jump process. A numerical evaluation of the model is thus relatively simple which suggests to investigate it. SIMULATION AND COMPUTATION SIMULATION OF THE MODEL To simulate the model we will evaluate A in a bounded window W = A x [0, T]. To do it we need to consider a larger region. The shape of this region is derived from the form of sets Bs in Eq. 4, where — °° in the integral limit is substituted by a finite time. This is a good approximation because of the exponential term in Eq. 4. Start with a discretization step A > 0 in space and time. T _ 48 Image Anal Stereol 2008;27:47-52 Typically W is a parallelepiped divided in cubes, each (i, j,k)-th cube is represented by its central point. We put in Eq. 4, V) Bs(x1,x2) = [x1+s,x1-s]x[x2 + s,x2-s],s<0 to get an approximation k+n i+r-k (5) A(i,j,k) = ^e^A L I Z{r,l,q), (6) r=k l=i-r+k q=j-r+k where Z(r,l,q) are independent Poisson distributed random variables with mean equal to a constant multiple of the volume of a cube. If this constant is a > 0 we have the n-th term in Eq. 6 of mean value (2n-1)2aA3e nA (7) which should be small relative to evaluated A values. To simulate an inhomogeneous A we may vary the parameter of the Poisson distribution in cubes. A simple model of the rat path Y is considered on a discrete square arenaA = {1,...,m}2 putting Yt a symmetric random walk on A with reflecting walls. The animal starts from the position (m+1/2,m+1/2), m odd, with constant speed in each time interval moves a step A in a random direction on the square grid with equal probability 1/4. Y is here independent of A. For the evaluation of Fano factor in Eq. 1, we get an approximation T M __ / IB (yt ) A {yt, t) dt~AY lB(yq)A(yq, q) 0 q=1 (8) where yq is the location at time q. Empirical mean, variance is substituted in Eq. 1, obtained from N realisations of A, respectively. Finally we simulate the spatio-temporal doubly stochastic Poisson point process X of spikes. Given A. the number of spikes during time q, q + A is Poisson distributed with mean AA(yq,q). NUMERICAL EXAMPLE We present numerical results of the simulation for the grid size m = 9. The grid A is subdivided onto boxes Aij, i,j= 1,..., 3, see the thick lines in Fig. 1a. Further parameters are N = 20 (number of realisations of A), M = 1000 (number of time steps), A = 0.1, n = 65 in Eq. 6. We thus obtain the term (Eq. 7) equal to 0.025a. We choose a equal to 100 in the lower left box A31 and equal to 1 elsewhere to obtain an inhomogeneous A. Thus the n-th term in Eq. 6 has mean value 0.02636, which is small in comparison to A values in Fig. 1b. L_L. d) /\ A rm ;v JL °Wii \» ^ /v HP Fig. 1.3x3 grid of boxes, a) simulated random walk with spikes, b) a simulated inhomogeneous intensity, c) Fano factors, d) histogram of the numbers of spikes, e) histogram of time units spent in a box, f) estimated expected intensity. The data and results (estimation of Fano factor and conditional intensity) are in Fig. 1. The histograms in d) and e) are naturally of a similar shape resulting from a). The Fano factors estimated in c) are clearly above 1. Using Eq. 2 the estimated expected intensity in f) corresponds to the theoretical intensity in b). Unfortunately for further characteristics estimation, e.g., for the variance of the driving intensity, there are no simple estimators available. Therefore we proceed using the filtering techniques. FILTERING BACKGROUND AND DISCRETIZATION Another rigorous approach to the problem of the inference of the driving intensity A and its characteristics is the filtering (Fishman and Snyder, 1976; Lechnerova´ et al., 2008). Given a realisation of a spatio-temporal point process X of spikes and given a trajectory of the path Y, the solution of the nonlinear filtering problem is the conditional expectation E [A \X, Y] which is not explicitly available. Here the Bayes formula for probability densities is useful f&\x,y)~f(x\k,y)f(k\y), since from the definition of the doubly stochastic Poisson process f(x\X,y) is a density of an inhomogeneous Poisson process with intensity function A. We will assume in the following that A,Y are stochastically independent. Then, f(A|x)ocf(x|A,y)f(A) (9) 49 Benesˇ V et al: Spatio-temporal point processes and the aim is to simulate a sample from the density f{X\x) which enables to solve the filtering problem and estimate empirically any characteristics of A including the overdispersion. Simulation is possible using Markov chain Monte Carlo (MCMC) techniques by means of one of the following approaches: the use of a) discrete grid, b) point process densities w.r.t. unit Poisson process. Consider a parallelepiped W C R3. In the approach a) we consider a spatio-temporal grid of cubes with centroids {c n k} C W as above and denote Xijk = A(ci jk). Let ijk be the counts of spikes in each cube. Then for the joint conditional densities it holds f{{Kjk}\{nijk}) x f({ni jk} | {hk})f({hk}) « Ilffni jk I Kk)Y[f{{kjk, ij} I {Ai j,k-1,i,j}) . ijk k In the right hand side of the last formula there is a product of marginal Poisson probabilities f(ni jk Xijk) (because of independence property of the Cox process given intensity) and a second product of transition densities (from the Markov property in time of a spatio-temporal OU process, see Barndorff-Nielsen and Schmiegel, 2004). In this second product {Aijk,i ,j} is the set of all Xijk where i,j vary along the index set, but k is fixed (while in {Ai jk} all indexes vary). The use of MCMC requires a closed formula for the transition densities, which may have atoms and are hard to evaluate . The algorithm was demonstrated for a different model of spatio-temporal log-Gaussian Cox process in Brix and Diggle (2001). The approximation works well if the intensity varies slowly w.r.t. the grid step size which is not the case for the jump model (Eq. 4). Therefore we look for another method. POINT PROCESS APPROACH TO FILTERING Because of the jump character of the model (Eq. 4) with Poisson basis, the second approach to filtering, based on the point process of jumps of Z is available. Let W = A x [0,T], denote the data x = {t7} where each %j reflects time and location of a spike on y. In Eq. 9, we have now with intensity function X. For the driving intensity function we have from Eq. 4 a representation, f(x|A,y)ocexp V 0 X{yt,t)dt TT A(Ti Tiex with the proportionality constant independent of both x and X. It is a density (w.r.t. a unit Poisson process) of an inhomogeneous spatio-temporal Poisson process A(v,0 = Xetj-tlfl,_t(v)(s7) tj 0, which represents A. The Metropolis-Hastings birth-death algorithm (Moller and Waagepetersen, 2003) can be used to simulate the MCMC chain {E,j,tjY, l = 0,...,J, where l is the index of iteration and J is the total number of iterations. The distribution of the chain tends to the desired conditional distribution (Eq. 9). Using ergodicity properties of the chain we can estimate statistical characteristics of A. T 50 Image Anal Stereol 2008;27:47-52 ESTIMATION OF CONDITIONAL INTENSITY CHARACTERISTICS Using the MCMC iterations of the auxiliary point process we evaluate the l-th iteration L(l)(v,t) from Eq. 10. Let 0 < K < J be the burn-in of the chain, k = J -K, and let a single realisation of Y,X be given. Analogously to Eq. 1, we pay attention here to time averaged characteristics, but any quantity derived from L can be obtained. For the discrete random walk with M time steps let 1 M-1 L(l)(v) = 1a L« (v, qA) M q=0 we get the conditional expectation of intensity estimated as, EL(v) 1 a L ^W(v (11) l=K+1 and the conditional variance of intensity estimated as varL(v) 1 k-1 J 2 a ( L^ > (v) - E L(v) J . (12) l=K+1 From N independent realisations of Y,X by evaluating corresponding N MCMC chains as above and averaging over realisations we can obtain estimates of unconditional quantities E[L(v,t)], var[L(v,t)]. NUMERICAL RESULTS AND CONCLUSIONS We simulated an experiment Y,X with homogeneous intensity L and used the point process approach to the filtering. A grid of 2 x 2 boxes was used with 20 time steps, m = 2,A = 1,a = 1 and BJx 1,x 2) = fx1 + us, x 1 - us1 X \x 2 + us,x2 - us1, s < 0, where u = 5. The MCMC chain of time averaged conditional intensity is in Fig. 3 together with its characteristics. We can observe an approximately Gaussian character in b) of the chain in a) and its slow mixing (decrease of autocorrelation function) in c). Denote Eij the mean in Eq. 11 and s ii the standard deviation (square root of Eq. 12). The results obtained on the grid are Eij sij 0.746 0.758 0.768 0.714 0.097 0.099 0.101 0.104 (13) (14) 200 400 600 800 b) c) 0.5 0.7 0.9 1.1 0 5 15 25 Fig. 3. 85×103 iterations of MCMC chain in the first of 4 boxes, a) the conditional expectation of intensity and its b) histogram, c) autocorrelation function. The horizontal axes in a) and c) correspond to number of iterations times 100. The horizontal axis in b) and the vertical axis in a) correspond to conditional expectation of intensity values. The blue line level in a) is the time averaged realisation of Eq. 10 at the corresponding grid point, used for the simulation of data X. The red line is the resulting level E11 = 0.746 in Eq. 13. The bounds in a) correspond to E11 ±2s11 where s11 = 0.097 is the standard deviation in Eq. 14. In the paper we discussed the stochastic modelling and simulation of an action consisting of random events which appear in time and space, with a similarity to a neurophysiological experiment. We developed the filtering approach to the estimation of conditional characteristics of the driving intensity of the spatio-temporal Cox process model. From the two approaches to filtering the one based on discretization is analytically and algorithmically demanding, since the transition distribution is hard to evaluate. On the other hand the presented point process approach is only computationally demanding, which does not make serious problems. An ultimate goal is a real data evalution where the questions of parameter estimation and degree of fit of a suitable model will arise. However these problems can be solved within the presented Bayesian MCMC methodology, in principle. a) 0 _ _ _ _ 51 Benesˇ V et al: Spatio-temporal point processes ACKNOWLEDGEMENT The research was supported by the projects ˇ MSM0021620839,IAA101120604 andGACR201/05/ H007. REFERENCES Barndorff-Nielsen O, Schmiegel J (2004). Le´vy based tempo-spatial modelling; with applications to turbulence. Usp Mat Nauk 159:63–90. Brix A, Diggle P (2001). Spatio-temporal prediction for log-Gaussian Cox processes. J Royal Statist Soc B 63:823– 41. Daley DJ, Vere-Jones D (1988). An introduction to the theory of point processes. New York: Springer. Eden UT, Frank LM, Barbieri R, Solo V, Brown EN (2004). Dynamic analyses of information encoding by neural ensembles. Neural Comput 19, 5:971–98. Fano U (1947). Ionization yields radiations. II. The fluctuations of the number of ions. Phys Rev 72:26–9. Fenton AA, Muller RU (1998). Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proc Nat Acad Sci USA 95:3182–7. Fishman PM, Snyder D (1976). The statistical analysis of space-time point processes. IEEE Trans Inf Theory 22:257–74. La´nsky´ P, Vaillant J (2000). Stochastic model of the overdispersion in the place cell discharge. Biosystems 58:27–32. Lechnerova´ R, Helisova´ K, Benesˇ V (2008). Cox point processes driven by Ornstein-Uhlenbeck type processes. Methodol Comput Appl Probab, in print. DOI: 10.1007/s11009-007-9055-1 Moller J, Waagepetersen R (2003). Statistics and simulations of spatial point processes. Singapore: World Sci. Prokesˇova´ M, Hellmund G, Vedel Jensen E (2006). On spatio-temporal Le´vy based Cox processes. In: Lechnerova´ R, Saxl I, Benesˇ V, eds. S4G. Proceedings of the International Conference on Spatial Statistics, ˇ Stereology, Stochastic Geometry. Praha: JCMF, 112–7. Schoenberg FP, Brillinger DR, Guttorp PM (2002). Point processes, spatial-temporal. In: El-Shaarawi A, Piegorsch W, eds. Encyclopedia of Environmetrics, vol. 3. New York: Wiley, 1573–7. 52