VIEW Communicated by Valerie Ventura
Computing Confidence Intervals for Point Process Models
Sridevi V. Sarma
sree@jhu.edu
Department of Biomedical Engineering, Johns Hopkins University, Baltimore,
MD 21218, U.S.A.
David P. Nguyen
dpnguyen@alum.mit.edu
Neuroscience Statistics Research Laboratory, Department of Anesthesia and
Critical Care, Massachusetts General Hospital, Boston, MA 02114, and
Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology,
Cambridge, MA 02139, U.S.A.
Gabriela Czanner
gabriela.czanner@warwick.ac.uk
Warwick Manufacturing Group and Warwick Medical School,
University of Warwick, Coventry CV4 7AL, U.K.
Sylvia Wirth
swirth@isc.cnrs.fr
Institut des Sciences Cognitives, CNRS, 69675 Bron, France
Matthew A. Wilson
mwilson@mit.edu
Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology,
Cambridge, MA 02139, U.S.A.
Wendy Suzuki
wendy@cns.nyu.edu
Center for Neural Science, New York University, New York, NY 10003, U.S.A.
Emery N. Brown
enb@neurostat.mit.edu
Neuroscience Statistics Research Laboratory, Department of Anesthesia,
Critical Care and Pain Medicine, Massachusetts General Hospital, Boston,
MA 02114; Harvard?MIT, Division of Health Sciences and Technology,
Massachusetts Institute of Technology, Cambridge, MA 02139;
and Department of Brain and Cognitive Sciences, Massachusetts Institute
of Technology, Cambridge, MA 02139, U.S.A.
Sridevi Sarma and David Nguyen contributed equally to this view.
Neural Computation 23, 2731?2745 (2011)
c?
2011 Massachusetts Institute of Technology
2732 S. Sarma et al.
Characterizing neural spiking activity as a function of intrinsic and ex-
trinsic factors is important in neuroscience. Point process models are
valuable for capturing such information; however, the process of fully
applying these models is not always obvious. A complete model applica-
tion has four broad steps: specification of the model, estimation of model
parameters given observed data, verification of the model using good-
nessoffit,andcharacterizationofthemodelusingconfidencebounds.Of
thesesteps,onlythefirstthreehavebeenappliedwidelyintheliterature,
suggesting the need to dedicate a discussion to how the time-rescaling
theorem,incombinationwithparametricbootstrapsampling,canbegen-
erallyusedtocomputeconfidenceboundsofpointprocessmodels.Inour
first example, we use a generalized linear model of spiking propensity to
demonstrate that confidence bounds derived from bootstrap simulations
are consistent with those computed from closed-form analytic solutions.
In our second example, we consider an adaptive point process model
of hippocampal place field plasticity for which no analytical confidence
bounds can be derived. We demonstrate how to simulate bootstrap sam-
ples from adaptive point process models, how to use these samples to
generate confidence bounds, and how to statistically test the hypothesis
that neural representations at two time points are significantly differ-
ent. These examples have been designed as useful guides for performing
scientific inference based on point process models.
1 Introduction
Receptivefieldsofneuronsinthebrainchangeinresponsetoenvironmental
stimuli as well as learning. For example, in the cat visual system, retinal
lesionsleadtoreorganizationofcorticaltopography(Pettet&Gilbert,1992).
Peripheral nerve sectioning can substantially alter the receptive fields of
neurons in monkey somatosensory and motor cortices (Kaas, Merzenich,
& Killackey, 1983; Merzenich et al., 1984). Similarly, the directional tuning
of neural receptive fields in monkey motor cortex changes as the animal
learns to compensate for an externally applied force field while moving
a manipulandum (Gandolfo, Li, Benda, Schioppa, & Bizzi, 2000). In the
rat hippocampus, the pyramidal neurons in the CA1 region have spatial
receptive fields. As a rat executes a behavioral task, a given CA1 neuron
fires only in a restricted region of the experimental environment, termed
the cell?s spatial or place receptive field (O?Keefe & Dostrovsky, 1971). Place
fields change in a reliable manner as the animal executes its task (Mehta,
Barnes, & McNaughton, 1997; Mehta, Quirk, & Wilson, 2000). Analysis
of such neural dynamics is crucial for understanding how different brain
regions update neural representations with learning and experience.
These examples highlight the need for models that integrate neuro-
physiological knowledge with internal and external covariates and capture
Computing Confidence Intervals for Point Process Models 2733
time-dependent processes of spiking activity, such as learning, adaptation,
andplasticityonthemillisecondtimescales.Thepointprocessframeworkis
abletomeettheseneeds,asdemonstratedoverawidespectrumofneuronal
data types (Brown, Nguyen, Frank, Wilson, & Solo, 2001; Eden, Frank, Bar-
bieri, Solo, & Brown, 2004; Ergun, Barbieri, Eden, Wilson, & Brown, 2007;
Sarma et al., 2010) because it naturally extends the Poisson process and
can therefore utilize a wealth of existing statistical solutions (Ogata, 1988;
Brown, Barbieri, Ventura, Kass, & Frank, 2002; Daley & Vere-Jones, 2003).
Here, in particular, we focus on a general solution to the problem: What is
the probable range of my model parameters, and how do I perform signifi-
cance testing on my data given a point process model?
We describe a general procedure for computing confidence intervals of
point process models based on the parametric bootstrap, a method that
first fits a parametric statistical model to observed data and then uses that
model to simulate surrogate data (Davison & Hinkley, 1997). The term boot-
strapping generally refers to the process of using surrogate data for the
purposes of estimation and inference. In the following sections, we show
that by using the time-rescaling theorem (Ogata, 1988; Brown et al., 2002;
Daley & Vere-Jones, 2003), a fundamental property of the point process, it is
possible to draw surrogate spike trains that are statistically probable given
the real observed spike train and a point process model of good fit. Surro-
gate spike trains are then pooled and processed to construct bootstrapped
conditional probability densities from which the confidence intervals for
model parameters are derived (Efron & Tibshirani, 1993; Gilks, Richardson,
& Spiegelhalter, 1998; Brown et al., 2002; Gentle, 2003).
We consider two examples where parametric bootstrapping is used to
infer the confidence bounds of estimated model parameters that describe
spike rate. The first is an analysis of a neuron recorded from the medial
temporal lobe of a monkey performing a visual saccade task (Efron & Tib-
shirani, 1993; Gilks et al., 1998; Brown et al., 2002; Gentle, 2003). The model
of the neuron relates intrinsic effects such as the neuron?s own spiking
history and concurrent ensemble activity, as well as the extrinsic effects of
the presented stimuli (Efron & Tibshirani, 1993; Gilks et al., 1998; Brown
et al., 2002; Gentle, 2003). We chose this particular example for its closed-
form solution, which allowed us to compute confidence bounds using a
theoretical maximum likelihood solution and a bootstrap estimation proce-
dure, and demonstrate that the two provide consistently similar confidence
bounds. The second model is motivated by the desire to improve the tem-
poral resolution of models that quantify neuronal dynamics. Using point
process modeling and steepest-descent adaptive filtering, we characterize
experience-dependent changes in one rodent hippocampal place receptive
field on a millisecond timescale (Brown et al., 2001; Frank, Eden, Solo, Wil-
son, & Brown, 2002). Although this adaptive model has many advantages,
it has no closed-form solution and therefore no closed-form solution for
computing confidence intervals. Using the parametric bootstrap sampler,
2734 S. Sarma et al.
we demonstrate how to obtain confidence-bound estimates for this model
and perform hypothesis testing to show that receptive field plasticity is
statistically significant over time.
2 Methods
2.1 Point Processes. Point process models of neural spike train data are
well suited for quantifying interactions among spike trains, environmental
factors,andintrinsicfactors(Brownetal.,2001;Edenetal.,2004;Ergunetal.,
2007;Sarmaetal.,2010).Hereweestablishtheframeworkforbootstrapping
(or simulating) structured spike trains from point process models (we refer
readers to more comprehensive works for theoretical underpinnings Cox &
Isham, 1980; Snyder, Miller, & Snyder, 1991; Daley & Vere-Jones, 2003).
Inapointprocessframework,neuralspikedataarerepresentedbybinary
random events that occur in continuous time and over the observation
interval (0, T]. Values of 1 represent the time of one spike, and values of
0 represent no spiking activity. The spike times are parameterized by the
set u
1
,u
2
,...,u
N
T
where N
T
represents the number of counted spikes in
the observation interval. Given a time t in the observation interval, the
instantaneous rate of spiking is given by the conditional intensity function
(CIF),
?(t | H(t),?(t)) = lim
Delta1?0
P(N(t + Delta1) ? N(t) = 1 | H(t),?(t))
Delta1
, (2.1)
where spiking activity is a function of spike history H(t) and time-varying
model parameters ?(t).Thespikehistorymaybeusedtoimposearefrac-
tory period or oscillatory dynamics, while time-varying parameters allow
the modeling of phenomena such as plasticity and learning on a single-
unit level. In the time interval, (t,t + Delta1), the link between the CIF and the
random number of spikes emitted by a neuron is the probability density
Poisson(?(t | H(t),?(t))Delta1).
Given an experimentally observed spike train, one goal of point process
modeling is to determine the functional form of the CIF and estimate pa-
rameters values that allow the CIF to match observed spiking activity. The
instantaneous log likelihood,
l(?(t) | H(t)) = log[?(t) | H(t),?(t))]
dN(t)
dt
? ?(t | H(t),?(t)), (2.2)
defines an optimization surface that is suitable for parameter estimation
of stationary and adaptive and dynamic models. For the purposes of com-
puter implementation, we approximate the instantaneous log likelihood in
discrete time for t = iDelta1 as
l(?
i
| H
i
) ? log[?(H
i
,?
i
)]n
i
? ?(H
i
,?
i
), (2.3)
Computing Confidence Intervals for Point Process Models 2735
where the interval Delta1 is on the order of milliseconds and smaller than any
observed interspike interval, and n
i
is the number of spikes observed in
the time interval ((i ? 1)Delta1,iDelta1]. For a small Delta1,itcanbeshownthatthedis-
crete instantaneous log likelihood is equivalent to a Bernoulli process of the
form P(n
i
| H
i
) ? (?
i
Delta1)
n
i
(1 ? ?
i
Delta1)
1?n
i
(Truccolo, Eden, Fellows, Donoghue,
& Brown, 2005). Therefore, for adequately small time steps, instead of com-
puting the Poisson density above, we need only compute the following
approximation to characterize spiking activity for time index i:
P(n
i
= 1 | H
i
,?
i
) ? ?(H
i
,?
i
)Delta1. (2.4)
The time-rescaling theorem, a cornerstone of point process theory, states
that if a spike train is generated from (or is consistent with) a CIF, then it is
possible to transform arbitrarily distributed interspike intervals into inde-
pendent and identically distributed (i.i.d.) exponential random variables,
?
k
, by essentially weighting time with the CIF:
bracketleftbigg
?
k
=
integraldisplay
u
k+1
u
k
?(t |?) dt
bracketrightbigg
? Exponential(1). (2.5)
Additionally, the interspike intervals can be transformed into i.i.d. uni-
formly distributed random variables,
[v
k
= 1 ? exp(?
k
)] ? Uniform(0, 1), (2.6)
thus making goodness-of-fit diagnostics, such as the Kilmogorov-Smirnov
(K-S)testandautocorrelationanalyses,astraightforwardprocedure(Ogata,
1988; Brown et al., 2002; Daley & Vere-Jones, 2003).
Our focus now is to demonstrate an equally important, yet less known,
application of the time-rescaling theorem. That is, the theorem may be gen-
erally applied to continuously differentiable point process models to ob-
tain confidence bounds on model parameter estimates (Ogata, 1988; Brown
et al., 2002; Daley & Vere-Jones, 2003). The statistical framework embod-
ied in equations 2.1 to 2.6 is fully compatible with Monte Carlo methods,
thus allowing bootstrapped spike trains to be generated from the nonho-
mogeneous Poisson process that is parameterized by a time-varying CIF,
?(t | H(t),?(t)). These samples may then be used to construct confidence
bounds for any test statistic that is based on the parameters of the CIF, ?(t),
and history H(t).
2.2 BootstrapSamplingofPointProcessModels. Theparametricboot-
strap sampler, a particular type of sampling method, may be used for in-
ference (Davison & Hinkley, 1997). Here, an analytical model is estimated
from the data in order to provide a basis for deriving statistics of interest
analytically or generating simulated data samples with similar statistical
qualities as the original data.
2736 S. Sarma et al.
In the case of neuronal spike trains, a parametric bootstrap sampler
allows confidence bounds to be computed around the instantaneous rate
givenbytheCIFestimate;suchconfidenceboundsallowsustobetterquan-
tify the statistical significance of firing rate comparisons between experi-
mental conditions. In this section, we outline the algorithm for parametric
bootstrap sampling and then apply it to two real data examples.
Thegoalofourbootstrapsampleristoestimatetheunderlyingprobabil-
ity distribution of the CIF parameters, ?(?), in light of the history of spiking
H(?).ByutilizingthetimerescalingtheoremdescribedinBrownetal.(2002),
we can generate b = 1,...,B spike train samples, and for each sample, the
CIF model parameters,
?
?
b
i
, are estimated for a fixed b and 0 < iDelta1 ? T.The
K-S and independence tests are used to validate the bootstrapped model
parameters (Johnson & Kotz, 1969). From the bootstrap samples, we may
generate confidence intervals for some test statistic of interest, ?(?). The
value of B is both large and determined online by convergence criteria
described below.
Let b=1 to start, increment b by 1 with each completion of step 6, and let
b = B denote the last sample to be generated:
1. Let k index be the spike to be generated. Initialize the algorithm by
setting k = 1 and ?u
b
0
= 0.
2. Draw ?v from a uniform distribution on [0, 1], and set ??
b
k
=?log(?v).
3. Find the sampled spike time, ?u
b
k
, as the solution to ??
b
k
=
integraltext
?u
b
k
?u
b
k?1
?(t | H(t),
?
?(t)) dt . It is important to note that the spiking his-
tory, H(t), in the integral expression is accumulated as you sample.
4. If ?u
b
k
> T, then discard ?u
b
k
,letN
b
(T) = k ? 1, and go to step 5; oth-
erwise increment k by 1 and return to step 2. Let the collection of
simulated spike times from each bootstrap iteration be denoted as
?
H
b
(t) ={?u
b
k
fork = 1,2,...,N
b
(t)}.
5. Using the simulated spike history,
?
H
b
(t), reestimate the parameter
vector
?
?
b
(t) using the same estimation procedure used to obtain
?
?(t),
and check for goodness of fit.
6. Let the sample statistic be defined as ??
b
(t) = ?(t |
?
H
b
(t) ?
?
?
b
(t)) for
b ? [1, B].
The topic of convergence deals with how to determine a value for B
that allows the simulation to fully sample from the distribution of ?(t)
for all times of interest. Many solutions have been proposed, especially in
the field of Markov chain Monte Carlo (MCMC) simulation (Gilks et al.,
1998). While we can learn from the MCMC literature, we are not creating
a Markov chain in this simulation because each sample does not depend
on the previous one. Thus, we do not consider the need to define burn-in
periods and problems of being trapped in state spaces of local optima.
Hence, it is possible to monitor convergence in our simulation using
a simple approach. For every tenth bootstrap sample, we computed the
Computing Confidence Intervals for Point Process Models 2737
emprical 2.5% and 97.5% points from our test statistic set, [??
1
(t),...,??
b
(t)],
where time t is chosen to reflect a time of interest in the experiment. When
the dynamic range of each percentile is bounded by a small interval, ?,the
simulation is said to have converged (Billingsley, 1999).
When the choice of ? is not clear, it is possible to use more automated
methods of convergence detection. Many of these approaches require the
simulation of M parallel bootstrap simulations, with the idea that at the
point of convergence, the statistical properties of each simulation will be
the same. One method in particular, the potential scale reduction factor
(PSRF), computes the square root of the ratio of the between-simulation
posterior variance estimate and the within-simulation variance to deter-
mine convergence (Gilks et al., 1998). As the simulation progresses, the
PSRF should begin to decrease and eventually reach 1 at convergence.
3 Results
3.1 ComparisonofBootstrapandAnalyticalConfidenceBounds. The-
oretically,bootstrapestimationgenerallycanbeusedwithanypointprocess
model with an integrable CIF (Brown et al., 2001). However, the application
of this theory is not always straightforward. We first provide an example
using a generalized linear model (GLM) framework (McCullagh & Nelder,
1990), which is becoming increasingly useful in neuroscience. In particular,
using the GLM form for the CIF, we compare our bootstrap-derived confi-
dence bounds with those obtained theoretically using maximum likelihood
estimation.
InaGLM,thelogoftheCIFismodeledasalinearfunctionofparameters
thatmultiplyfunctionsofthecovariatesr(t)thatdescribetheneuralactivity
dependencies, that is,
log(?(t | ?)) =
J
summationdisplay
j=1
?
j
f
j
(r(t),t), (3.1)
where ?
j
? R for each j and f
j
are known basis functions. The model param-
eters are ?
j
for j = 1,2,...,J .
TheGLMisanextensionofthemultiplelinearregressionmodelinwhich
the variable being predicted, in this case spike times, need not be gaussian
(McCullagh & Nelder, 1990). The GLM also provides an efficient computa-
tional scheme for model parameter estimation and a likelihood framework
forconductingstatisticalinferencesbasedontheestimatedmodelanderror
bounds on model parameters (Brown, Barbieri, Eden, & Frank, 2003).
Czanner et al. (2008) set out to examine the patterns of neural activity
observed during the acquisition of new associative memories in monkeys
(seeFigure1).Theyproposedastate-spaceGLM(describedinsection2.1)of
the CIF to capture both intrinsic effects, such as the neuron?s own spiking
2738 S. Sarma et al.
Figure1: (A)Schematicillustrationoflocationsceneassociationtask.(B)Raster
plotfromsceneepochtoendofeyemovementresponse.(C)Estimatedstimulus
coefficients of GLM overlaid with error bounds computed analytically (gray)
and via our bootstrap (black). (D) Kilmogorov-Smirnov plot for GLM model
with temporal history.
history and concurrent ensemble activity, and the extrinsic effect of the
stimulusbothwithinandacrosstrials.Usingtheirmodelappliedtoasingle
neuron spike train, we compute error bounds for a GLM point process
model using simulated and analytical methods.
In Czanner et al. (2008), the discretized CIF is modeled over an interval
of length T as the product of the stimulus intensity and the effect of the
history dependence of the neural process as
?(lDelta1 | ?
S
,?
H
, H
l
) = ?
S
(lDelta1 | ?
S
)?
H
(lDelta1 | ?
H
, H
l
), (3.2)
wherel = 1,2,...,L denotesthediscretetimestepthatcorrespondstocon-
tinuous timet ? ((l ? 1)Delta1,lDelta1], Delta1 = TL
?1
; H
l
defines the spiking history up
totimelDelta1;?
S
= (?
S
1
,?
S
2
,...,?
S
L
)areparametersoftheGLMthatcharacterize
thestimulusortaskeffectonthespikingactivity;and?
H
= (?
H
1
,?
H
2
,...,?
H
J
)
is a constant vector of parameters that describe the dependence of current
spiking activity on the spiking history.
The GLM model of the stimulus effect is
log(?
s
(lDelta1 | ?
s
)) =
L
summationdisplay
r=1
?
s
r
g
r
(lDelta1), (3.3)
Computing Confidence Intervals for Point Process Models 2739
where g
r
(lDelta1) is a sequence of unit pulse functions with equal length, that
is,
g
r
(lDelta1) =
braceleftbigg
1ifl = r
0 otherwise.
(3.4)
Since we compute only error bounds for the stimulus parameters, we
omit the details of the history effects, which are described in Czanner et al.,
(2008).
We use data from one neuron for which a static GLM model provides a
good fit. Twenty-two trials of length 3.2 seconds were used to develop the
modelthatcapturedtheCIFfrom0.3to3.5seconds?thus,thelengthofthe
stimulus parameter vector L = 32. The time bin length Delta1 = 1 msec. The K-S
plot for the model is shown in Figure 1D, which stays within the computed
95% confidence bounds for the degree of agreement using the distribution
of the K-S statistic.
3.1.1 Analytic Approach. Since
?
? is a maximum likelihood (ML) estimate
of ?, it can be shown to asymptotically converge to a gaussian distribution
with mean
?
? (the asymptotic estimate that converges to the true ?) and vari-
ance ?
2
i
= I(
?
?)
?1
, where I(
?
?) is the expected Fisher information evaluated
at
?
? (Brown et al., 2003). Also, since unit pulse functions (g
r
)areused,itis
straightforward to show that the
?
?
S
i
prime
s are mutually independent (Czanner
et al., 2008). Hence, 95% confidence intervals can be computed individually
as
?
?
i
? z
i
such that z
i
= ?
i
Phi1
?1
(0.975).
3.1.2 Bootstrap Approach. Apply the bootstrap algorithm in section 2.2 to
generate K samples of spike trains that comprise one bootstrap sample:
1. We first fix ?
H
and compute
?
?
S
from the bootstrap sample.
2. Repeat the process to generate B = 1000 bootstrapped ML estimates
of ?
S
.
3. Checktheconvergenceofthebootstrapsampledistribution.Ifneeded,
generate more bootstrap samples such that B = B + 1000.
4. Compute the 95% confidence interval from all the bootstrap samples
via the 5th and 95th percentiles of the samples.
a. Given a time step l from our model above, for each bootstrap
sample indexed by s, compute the CIF conditioned on the sample:
?
?
l
(
?
?
S
) ? ?(lDelta1 |
?
?
S
,?
H
, H
l
). The result will be B samples of
?
?
l
(
?
?
S
).
b. Given these B samples, we compute the confidence interval for
each l using percentiles of 2.5 and 97.5.
Figures 1B and 1C illustrate the raw data in the form of raster plots
for each neuron and corresponding stimulus coefficients overlaid with
95% confidence bounds computed both analytically and via bootstrap
simulation. As seen in Figure 1C, the confidence interval obtained by using
2740 S. Sarma et al.
the asymptotic normality of a maximum likelihood estimate matches the
parametric bootstrap confidence interval for that estimate (which is ex-
pected because both procedures rely on the same model). This result helps
to confirm that our parametric bootstrap design will provide the theoreti-
cally comparable confidence bounds.
3.2 Spline-Based Adaptive Point Process Modeling of Place
Cells. Previously we developed and applied adaptive spline-based point
process filters to model the receptive field plasticity of hippocampal place
cells (Brown et al., 2001; Frank et al., 2002). The spline representation is
flexible, requires few assumptions, and lends itself to an adaptive estima-
tion framework. These applications demonstrated that adaptive estimation
of CIF parameters may yield consistent data model agreement; however,
we did not demonstrate how to obtain error bounds for our CIF param-
eter estimates. Here, we briefly reiterate the nonparametric, spline-based
CIF model and steepest-descent estimator, and then apply the bootstrap
method for obtaining error bounds.
3.2.1 Theory. Lettimebediscretizedsuchthatt = iDelta1andletthesubscript
on the function denote the time index for all parameters, for example,
F
i
(r) = F(r
i
) = F(r(iDelta1)). The proposed CIF for hippocampal place cells is
a function that is decomposed into a spatial and temporal component,
?
i
(x,? | H
i
,?
S
i
,?
T
i
) = ?
S
i
(x | ?
S
i
)?
T
i
(?/?
T
i
). (3.5)
The spatial dependency, ?
S
i
(?), depends on x, which is the location of
the rat on a linear track (see section 3.3.2 for the experimental description).
The temporal dependency, ?
T
i
(?), depends on ?
i
= iDelta1 ? t
lastspike
which in-
corporates the history dependence and is defined as the time since the last
spike. For convenience, we state the time dependency of the parameters in
equation 3.5 using the index i.
The two components of our CIF have the same functional form, a
Catmull-Rom or cardinal spline, which we first define generally as ?(?/?)
(Bartels, Beatty, & Barsky, 1987), where ? will later be replaced by x or ?.
A priori, we define the range supported by the spline, (?
1
,?
J
], and J +
2 uniformly spaced, time-varying control points, ?
i
={?
i,0
,?
i,1
,...,?
i,J +1
},
located at the fixed locations {?
0
,?
1
,...,?
J +1
}. The spline model is
?(? | ?)=[?(?)
3
?(?)
2
?(?)1]
?
?
?
?
?
?0.51.5 ?1.50.5
1 ?2.52?0.5
?0.500.50
0100
?
?
?
?
?
?
?
?
?
?
?
i,j?1
?
i,j
?
i,j+1
?
i,j+2
?
?
?
?
?
,
(3.6)
Computing Confidence Intervals for Point Process Models 2741
for ? ? (?
j
,?
j+1
] and where ?(?) = (? ? ?
j
)/(?
j+1
? ?
j
). At any point in
the continuous range of (?
1
,?
J
], the spline is fully determined by the four
nearest control points.
The goal is to quantify the plasticity of the place-receptive field by se-
quentially estimating the vector ?
i
for each incremental time step iDelta1 =
Delta1,2Delta1,...toT (Brown et al., 2001). The steepest-descent solution attempts
to move
?
?
i
toward the solution that will maximize the instantaneous log
likelihood in equation 2.3 at each time step:
?
?
i
=
?
?
i?1
? ?
?l
t
(?)
??
vextendsingle
vextendsingle
vextendsingle
vextendsingle
?=
?
?
i?1
. (3.7)
Theparameter ? isanadaptivegainparameterthatdeterminesthespeed
ofadaptation.ThefullderivationofthisadaptivefiltercanbefoundinFrank
et al. (2002).
3.2.2 Application. In the following example, a rat runs back and forth
on a 1D track that is 330 cm in length. The activity of a pyramidal neu-
ron from the CA1 region of the dorsal hippocampus is observed during
behavior. In Figure 2A, the trajectory of the rat is shown in space-time by
the gray line, and the action potentials are shown by the black dots. We
focus this example on the condition where the rat is running toward the
0 location.
The spatial spline contains 33 controls points that span from ?11 cm to
342 cm. The temporal spline is defined on the log
10
scale with 35 control
points ranging from ?4 to 3 (equivalently 0.1 msec to 1000 sec). We initial-
izedthealgorithmasDelta1 = 16.7ms,?
S
= ?
T
= 60,?
S
0
= 0,and?
T
0
= 0.Werun
the adaptive filter forward and backward in time like a temporal smoother
until the values of
?
?
S
i
and
?
?
T
i
do not change significantly over successive
passes. We then validate the model using the K-S statistic (see Figure 2B).
We illustrate the utility of the point process sampler by computing the 95%
confidence bounds for
?
?
S
i
(?)and
?
?
T
i
(?) at time points t
1
= 9840 sec and t
2
= 10,714 sec by obtaining B = 500 bootstrap samples of
?
?
S
and
?
?
T
(see
Figures 2C and 2D). Convergence criteria of the variance of the 5th and
95th percentiles were met for B = 500.
In order to state if the change in the neuron?s receptive field is signifi-
cant,wechoseourteststatisticstobe ?
S
=
integraltext
x=342
x=11
[?
S
t
2
(x |
?
?
S
) ? ?
S
t
1
(x |
?
?
S
)]dx
and ?
T
=
integraltext
?
?
=3
?
?
=?4
[?
T
t
2
(?
?
|
?
?
T
) ? ?
T
t
1
(?
?
|
?
?
T
)]d?
?
, where
?
?
S
or
?
?
T
denotes a
bootstrap sample and ?
?
= log
10
(?). Under a null hypothesis where the
receptive has not changed between t
1
and t
2
, the distributions for ?
S
and ?
T
should be symmetric and centered at 0. In Figures 2E and 2F,
we find that the null hypotheses are not supported, with both p-values
equaling 0.
2742 S. Sarma et al.
Figure 2: (A) Rat position overlaid with spike events. Each dot marks an ac-
tion potential of the neuron. The two lines indicate points in time where the
error bounds of the conditional intensity function were produced. (B) The
Kilmogorov-Smirnov test demonstrating the adaptive filter model fit is ap-
propriate. (C) Conditional spatial intensity function for each of the two time
stamps.Theshadedregionscorrespondtothe95%confidencebounds.(D)Con-
ditionaltemporalintensityfunctionforeachofthetwotimestamps.Theshaded
regions correspond to the 95% confidence bounds. (E) Probability distribution
for the difference in the test statistic, ?
s
, between two times of interest shown in
C. (F) Probability distribution for the difference in the test statistic, ?
T
, between
two times of interest shown in D.
Computing Confidence Intervals for Point Process Models 2743
4 Conclusions and Future Work
We have discussed a general parametric bootstrap method for computing
confidenceintervalsofpointprocessfiltersthatmayormaynothaveclosed-
form expressions. Hypothesis tests, which are crucial in the process of sci-
entific inference, greatly rely on the accuracy of confidence intervals. Thus,
we demonstrated in our first example that confidence intervals computed
by bootstrapping are consistent with those that are theoretically derived. In
the second example, we showed that the parametric bootstrap also allowed
us to compute confidence bounds for an iterative filter and determine if
significant changes occurred in neuronal coding.
The time-rescaling theorem, an elegant property from point process the-
ory, is the basis for useful diagnostic measures such as the Q-Q test and
Kolmogorov-Smirnov test, as well as for the inversion sampler used here to
derive bootstrap samples from original data. It allows us to take any cross-
validatedpointprocessmodelthatisintegrableandsimulaterandomspike
trains that have statistically similar features to the observed data. It is im-
portant to note that while we have proposed to bootstrap samples using
inversion sampling, other methods of random sampling, such as rejection
sampling (Efron & Tibshirani, 1993; Gilks et al., 1998; Brown et al., 2002;
Gentle, 2003), may be used with possibly easier implementation. In addi-
tion, it is worth noting that the accuracy of the confidence intervals derived
by parametric bootstrapping is predicated on the proper validation of the
point process model by goodness-of-fit measures such as the K-S test and
tests of random independence.
The method presented here increases the value of an existing point
process model by allowing comparisons to be made within a data set
without having to alter the functional form of the CIF. In addition, de-
pending on the extent of the data set and model, it is entirely pos-
sible to perform multiple hypothesis tests under the same parametric
bootstrap simulation. Thus, with a well-designed analysis approach us-
ing point process modeling and parametric bootstrap, it may be pos-
sible to obtain additional scientific inferences that are statistically well
supported.
Acknowledgments
We thank Uri Eden for valuable discussions. Support was provided by
Burrough?s Wellcome Fund and L?Oreal for Women in Science fellow-
ship to S.V.S, NIH MH59733 and NIH Pioneer Award DP1 OD003646-01
to E.N.B, NIH DA015644 to E.N.B. and W.A.S., NIH MH58847 to W.A.S., a
McKnightFoundationgranttoW.A.S.,andaJohnMerckScholarsAwardto
W.A.S.
2744 S. Sarma et al.
References
Bartels, R. H., Beatty, J. C., & Barsky, B. A. (1987). An introduction to splines for use in
computer graphics and geometric modeling. San Francisco: Morgan Kaufmann.
Billingsley, P. (1999). Convergence of probability measures. Hoboken, NJ: Wiley.
Brown, E. N., Barbieri, R., Eden, U. T., & Frank, L. M. (2003). Likelihood methods for
neural data analysis. In J. Feng (Ed.), Computational neuroscience: A comprehensive
approach (pp. 253?286). Boca Raton, FL: CRC.
Brown, E. N., Barbieri, R., Ventura, V., Kass, R. E., & Frank, L. M. (2002). The time-
rescaling theorem and its application to neural spike train data analysis. Neural
Computation, 14(2), 325?346.
Brown, E. N., Nguyen, D. P., Frank, L. M., Wilson, M. A., & Solo, V. (2001). An
analysis of neural receptive field plasticity by point process adaptive filtering.
Proc. Natl. Acad. Sci. U.S.A., 98(21), 12261?12266.
Cox, D. R., & Isham, V. (1980). Point processes. London: Chapman and Hall.
Czanner, G., Eden, U. T., Wirth, S., Yanike, M., Suzuki, W. A., & Brown, E. N.
(2008). Analysis of between-trial and within-trial neural spiking dynamics. J.
Neurophysiol., 99(5), 2672?2693.
Daley, D. J., & Vere-Jones, D. (2003). An introduction to the theory of point processes.
New York: Springer.
Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application.
Cambridge: Cambridge University Press.
Eden, U. T., Frank, L. M., Barbieri, R., Solo, V., & Brown, E. N. (2004). Dynamic anal-
ysis of neural encoding by point process adaptive filtering. Neural Computation,
16(5), 971?998.
Efron,B.,&Tibshirani,R.(1993).Anintroductiontothebootstrap.NewYork:Chapman
& Hall.
Ergun,A.,Barbieri,R.,Eden,U.T.,Wilson,M.A.,&Brown,E.N.(2007).Construction
of point process adaptive filter algorithms for neural systems using sequential
Monte Carlo methods. IEEE Trans. Biomed. Eng., 54(3), 419?428.
Frank, L. M., Eden, U. T., Solo, V., Wilson, M. A., & Brown, E. N. (2002). Contrasting
patternsofreceptivefieldplasticityinthehippocampusandtheentorhinalcortex:
An adaptive filtering approach. J. Neurosci., 22(9), 3817?3830.
Gandolfo, F., Li, C., Benda, B. J., Schioppa, C. P., & Bizzi, E. (2000). Cortical correlates
of learning in monkeys adapting to a new dynamical environment. Proc. Natl.
Acad. Sci. U.S.A., 97(5), 2259?2263.
Gentle, J. E. (2003). Random number generation and Monte Carlo methods.NewYork:
Springer-Verlag.
Gilks, W. R., Richardson, S., & Spiegelhalter, D. J. (1998). Markov chain Monte Carlo in
practice. Boca Raton, FL: Chapman & Hall.
Johnson, N. L., & Kotz, S. (1969). Distributions in statistics. Hoboken, NJ: Wiley.
Kaas, J. H., Merzenich, M. M., & Killackey, H. P. (1983). The reorganization of so-
matosensory cortex following peripheral-nerve damage in adult and developing
mammals. Annual Review of Neuroscience, 6, 325?356.
McCullagh, P., & Nelder, J. A. (1990). Generalized linear models (2nd ed.). Boca Raton,
FL: CRC.
Computing Confidence Intervals for Point Process Models 2745
Mehta, M. R., Barnes, C. A., & McNaughton, B. L. (1997). Experience-dependent,
asymmetric expansion of hippocampal place fields. Proc. Natl. Acad. Sci. U.S.A.,
94(16), 8918?8921.
Mehta, M. R., Quirk, M. C., & Wilson, M. A. (2000). Experience-dependent asym-
metric shape of hippocampal receptive fields. Neuron, 25(3), 707?715.
Merzenich, M. M., Nelson, R. J., Stryker, M. P., Cynader, M. S., Schoppmann, A., &
Zook, J. M. (1984). Somatosensory cortical map changes following digit amputa-
tion in adult monkeys. J. Comp. Neurol., 224(4), 591?605.
O?Keefe, J., & Dostrovsky, J. (1971). The hippocampus as a spatial map: Prelim-
inary evidence from unit activity in the freely-moving rat. Brain Res., 34(1), 171?
175.
Ogata, Y. (1988). Statistical-models for earthquake occurrences and residual analysis
for point-processes. Journal of the American Statistical Association, 83(401), 9?27.
Pettet, M. W., & Gilbert, C. D. (1992). Dynamic changes in receptive-field size in cat
primary visual cortex. Proc. Natl. Acad. Sci. U.S.A., 89(17), 8366?8370.
Sarma, S. V., Eden, U. T., Cheng, M. L., Williams, Z. M., Hu, R., Eskandar, E.,
et al. (2010). Using point process models to compare neural spiking activity
in the subthalamic nucleus of Parkinson?s patients and a healthy primate. IEEE
Trans. Biomed. Eng., 57(6), 1297?1305.
Snyder, D. L., Miller, M. I., & Snyder, D. L. (1991). Random point processes in time and
space. New York: Springer-Verlag.
Truccolo, W., Eden, U. T., Fellows, M. R., Donoghue, J. P., & Brown, E. N. (2005). A
point process framework for relating neural spiking activity to spiking history,
neural ensemble, and extrinsiccovariate effects.J.Neurophysiol.,93(2), 1074?1089.
Received November 10, 2009; accepted April 26, 2011.
This article has been cited by: