Nonlinear Event-Related Responses
in
fMRI
Karl
J.
Friston, Oliver Josephs, Geraint Rees,
Robert Turner
This paper presents an approach to characterizing evoked
hemodynamic responses in fMRl based on nonlinear system
identification, in particular the use of Volterra series. The
approach employed enables one to estimate Volterra kernels
that describe the relationship between stimulus presentation
and the hemodynamic responses that ensue. Volterra series
are essentially high-order extensions of linear convolution or
“smoothing.” These kernels, therefore, represent a nonlinear
characterization of the hemodynamic response function that
can model the responses to stimuli in different contexts (in
this work, different rates of word presentation) and interac-
tions among stimuli. The nonlinear components of the re-
sponses were shown to be statistically significant, and the
kernel estimates were validated using an independent event-
related fMRl experiment. One important manifestation
of
these nonlinear effects is a modulation of stimulus-specific
responses by preceding stimuli that are proximate in time.
This means that responses at high-stimulus presentation
rates saturate and, in some instances, show an inverted
U
behavior. This behavior appears
to
be specific to
BOLD
ef-
fects (as distinct from evoked changes in cerebral blood flow)
and may represent a hemodynamic “refractoriness.” The aim
of this paper is to describe the theory and techniques upon
which these conclusions were based and to discuss the im-
plications for experimental design and analysis.
Key words: nonlinear system identification; functional neuro-
imaging; fMRI; hemodynamic response function; Volterra se-
ries.
INTRODUCTION
This paper
is
about evoked hemodynamic responses in
functional magnetic resonance imaging
(fMRI)
and how
the measured blood oxygen level dependent
(BOLD)
ef-
fects can be related to underlying neuronal activity. In
particular we investigate the nonlinear nature of this
response, the significance of the nonlinear components,
and how they affect the design and interpretation of fMRI
experiments.
In Friston
et al.
(I),
we presented a model of observed
hemodynamic responses, in fMRI time series, that obtain
when the underlying neuronal activity (inferred
on
the
basis
of
changing task conditions)
is
convolved
or
smoothed with a
hemodynamic response function.
This
model was subsequently elaborated in the context of the
general linear model
(2-4).
The general linear model is
variously employed in linear system identification, from
MRM
3941-52
(1998)
From the Wellcome Department
of
Cognitive Neurology, Institute
of
Neu-
rology, Queen Square, London, United Kingdom.
Address correspondence to: Karl
J.
Friston, The Wellcome Department
of
Cognitive Neurology, Institute
of
Neurology, Queen Square, London, UK
WCl N
3BG.
Received January
21,
1997; revised
June
17, 1997; accepted June 77,
1997.
This work was supported by the Wellcome Trust.
Copyright
0
1998 by Williams
&
Wilkins
All rights of reproduction in any form reserved.
0740-31 94/98 $3.00
41
a
signal-processing perspective,
or
multiple linear regres-
sion
or
AnCova in statistics. These approaches to mod-
eling and characterizing
fMRI
time series are predicated
on the assumption that the relationship between evoked
neuronal activity and the ensuing hemodynamic re-
sponse can be approximated by a linear convolution
using a fixed and time-invariant hemodynamic response
function
(1,
5).
This assumption has been evaluated by
comparing estimates
of
the hemodynamic response func-
tion using different stimuli
(6)
and linear system identi-
fication. In this paper, we present a nonlinear character-
ization of the hemodynamic response function using
nonlinear system identification and explicitly assess
both the significance and behavior of the nonlinear com-
ponents (where they exist).
We have employed a parametric experimental design
using evoked responses to words presented aurally at
varying frequencies. This variation allowed
us
to exam-
ine the “stability” of the hemodynamic response to single
words and the interactions among stimuli when pre-
sented close together. These interactions represent non-
linear effects that we were able to estimate and make
statistical inferences about. The importance of this work
relates to understanding the nonlinear relationship be-
tween evoked responses and sensory
or
behavioral pa-
rameters (such as presentation rate) and the implied con-
straints imposed upon experimental design and analysis.
For example, we were able to resolve the apparent dis-
crepancy between linear increases in blood
flow
in re-
sponse to increasing word presentation rates
(7)
and the
nonlinear dependency
of
the
BOLD
signal
(8).
From a
data analysis perspective, the framework described in
this paper can be seen as a generalization of linear ap-
proaches that characterize hemodynamic responses,
evoked by single events, in terms of basis functions of
peri-stimulus time
(9).
This paper is divided into three sections. The first
section describes the theoretical background to analyzing
nonlinear
or
dynamic systems using Volterra series as a
general model relating changes in neuronal activity
(or
stimuli) to hemodynamic responses. By using a second-
order expansion, we were able to reformulate the prob-
lem in terms of the general linear model and therein
provide for both parameter estimation and statistical in-
ference about the effects observed. The second section
uses the results of the first section to analyze two fMRI
experiments, comprising two single-subject aural-stimu-
lation paradigms. In the first,
epoch-related
experiment
blocks,
or
epochs, of words were presented at different
rates. On the basis of this experiment, we were able to
estimate a high-order
or
nonlinear hemodynamic re-
sponse function and use it to predict the response that
would have been evoked by a single word. In the second,
event-related
experiment, words were presented in iso-
lation. This allowed
us
to validate the estimated re-
sponses to single words from the first study in terms of
42
Friston
et
al.
empirically determined event-related responses from the
second. The third section uses the nonlinear model of
evoked responses of the previous sections to look in
detail at the nonlinear interactions.
By
performing
“vir-
tual” experiments on the model, we show that these
interactions can be thought of in terms of
a
hemodynamic
“refractoriness” in which a prior stimulus modulates the
response to a subsequent stimulus,
if
it occurs within a
second
or
so.
This modulation represents an interaction,
over time, between the responses to successive stimuli
and results in reduced responsiveness at high-stimulus
frequencies.
THEORETICAL BACKGROUND
Nonlinear System Identification
Neuronal and neurophysiological dynamics are inher-
ently nonlinear and lend themselves to modeling by non-
linear dynamic systems. However, due to the complexity
of biological systems, it is difficult to find analytic equa-
tions that describe them adequately [although see
Vazquez and No11
(10)
for a compelling example]. An
alternative
is
to take a very general model and obtain the
specific parameters that enable the model to describe the
system in question
(11).
A
common example of this func-
tional approach to system identification is the use of
Volterra series. The Volterra series is an extension of the
Taylor series representation to cover dynamic systems
and has the general form
y(t)
=
ho
+
. . .
and soon
fit)
is
the output, in this case the hemodynamic response
or
fh4RI
signal, and
u(t)
the input, in this case neuronal
activity as indexed by the stimulus rate employed in
our
experiments.
hl’(T1,.
. . .
7,)
is the nth order Volterra ker-
nel. It can be shown that these series can represent any
analytic time-invariant system
(1 1).
The Volterra series
has been described as a “power series with memory” [see
Chapters
2
and
3
of
Bendat
(12)
for a fuller discussion].
The problem of characterizing the relationship between
the stimulus function
or
neuronal activity
u(t)
and the
hemodynamic response
fit)
reduces to estimating the
kernel coefficients
h”.
For
long time series, with rela-
tively noiseless data, a number of approaches could be
tried. Wray and Green
(11)
describe a technique using
time-delay neuronal networks. In
our
experience, the
nature of
fMRI
data does not permit t.he use of such
techniques,
so
we have adopted a standard least squares
approach. This has the advantage
of
providing
for
statis-
tical inference using the general linear model (see be-
low). To do this, one must first linearize the problem.
Consider the second-order approximation to the above
expansion with finite “memory”
T:
y(t)
=
ha
h2(T1,T2)’
U(t
-
T1)
*
U(t
-
72)
*
dTldT2
+
II;
Note that the integrals start at zero. This reflects the fact
that
our
system is “causal” in the sense that neuronal
changes precede hemodynamic responses. In this formu-
lation, the first-order coefficients
hl(t)
correspond to the
[linear] hemodynamic response function as described in
Friston
et
al.
(1).
The new terms depend on the second-
order coefficients
h2
and are the primary focus of this
paper.
The second step in making the estimates of
ho, hl,
and
h2
more tractable, for noisy data like
fMRI,
is
to
expand
the kernels in terms of a small number
P
of temporal
basis functions This allows
us
to estimate the
coefficients of this expansion using standard least
squares:
97).
PP
Now define a new set of response variables
x,(t)
that
represent the original time series
u(t)
convolved with the
ith basis function
xi(t)
=
bi(~,)
-
u(t
~
T~)~T~
I
Substituting this expression into Eq.
[I]
and including an
explicit error term
e(t)
gives
P
PP
This is simply a general linear model with response
variable
fit),
the observed time series, and explanatory
variables
1,
xj(t),
and
xi(t).xj(t)
at the [discrete] times at
which they are observed. These explanatory variables
(convolved time series
of
neuronal activity
or
stimulus
presentation rate) constitute the columns of the
design
Nonlinear Even t-Rela
ted
Responses
-0.05
43
~
'.'
matrix.
The unknown parameters are
8,
g',
and
8
from
which the kernel coefficients
ho, hl,
and
h2
are derived,
using Eq.
[2].
Having reformulated the problem in this
way, we can now use standard analysis procedures de-
veloped for serially correlated
fMRI
time series that em-
ploy the general linear model
(2,
4).
These procedures
provide parameter estimates (i.e., estimates of the basis
function coefficients and, implicitly, the kernels them-
selves) and statistical parametric maps
(SPMs)
testing the
significance
of
a hemodynamic response at each and
every voxel. In this paper, we will
use
SPMs of the F
statistic
(SPM{FJ)
that test the joint contribution of effects
considered of interest (the remaining effects,
or
columns
of the design matrix, are called confounds). Below we
will present
SPM{FJs
testing for the significance of the
first- and second-order coefficients
hl
and
h2
and
SPM{F]s that test for the nonlinear effects
h2
alone by
treating the first-order effects
hl
as confounds. These
analyses involve using design matrices with and without
the effects of interest and assessing the reduction in error
with the F statistic.
In this work, we used only three basis functions, i.e.,
P
=
3
(Fig.
1).
These were gamma density functions
peaking during the early, intermediate, and late compo-
nents of the anticipated hemodynamic response. The
choice of these functions was motivated by prior knowl-
edge about the form of the [linear] hernodynamic re-
sponse function. This form is usually well approximated
by a linear combination of two
or
more gamma density
functions.
A
special case of this is the Poisson form
adopted in Friston
et al.
(1)
that corresponds to a single
gamma density with equal mean and variance. Clearly
the choice of basis functions is dictated by the nature of
the data and the amount of temporal detail that one
wants to model.
In
some instances (e.g., multislice acqui-
sition), there are differences in the times that one voxel
time series
is
acquired in relation to another. To accom-
modate these slight shifts in time, we often supplement
the basis functions with their temporal derivatives (Fig.
1).
The role of these derivatives can be seen intuitively by
basis
functions
0.25,
I
i
-0.1
'
I
0
5
10
15
20
25
30
time
(seconds)
FIG.
1.
Basis functions
b,(t)
(solid lines) and their derivatives (bro-
ken lines) used
in
the expansion of the Volterra kernels
h'(7,)
and
h2(7,,
7.J.
These are gamma density functions
with
mean and
variance
2'
(i
=
2,
3,
and
4).
These gamma density functions can
be thought
of
as a set of time-scaled Poisson functions because
their mean and variance are equal.
noting that adding
(or
subtracting) the temporal deriva-
tive shifts the basis function backwards
(or
forwards) in
time.
In
this paper, derivatives were only used in the
analysis of the event-related study where temporal ef-
fects were more acute.
THE fMRl EXPERIMENTS
Experimental Design and Data Acquisition
In this section, we apply the theory presented above to
fMRI time series obtained from a single normal male
subject during passive listening
to
words presented alone
or
continuously at different rates. The data were acquired
at
2
Tesla using a Magnetom VISION (Siemens, Erlangen)
whole body MRI system, equipped with a head volume
coil. Contiguous multislice T,*-weighted
fMRI
images
were obtained with a gradient echo-planar sequence
us-
ing an axial slice orientation
(TE
=
40
ms,
TR
=
1.7
s,
64
X
64
X
16
voxels
[19.2
X
19.2
X
4.8
cm]). After
discarding initial scans (to allow for magnetic saturation
effects), each time series comprised 1200 (first study) and
1000
(second study) volume images with 3-mm isotropic
voxels. In the first, epoch-
or
rate-related experiment, the
subject listened to monosyllabic or bisyllabic concrete
nouns (i.e., dog, radio, mountain, gate) presented at five
different rates
(10,
15,
30,
60, and
90
worddmin)
for
epochs of 34
s
(20
scans), intercalated with periods
of
rest. The five presentation rates were successively re-
peated according to a Latin Square design. In the second,
event-related study, the subject listened to (nonrepeat-
ing) nouns presented once every
34
s.
Data Preprocessing
The data were analyzed with SPM96 (Wellcome Depart-
ment of Cognitive Neurology, http://www.fil.ion.
ucl.ac.uk/spm). The time series were realigned, corrected
for movement-related effects, and spatially normalized
into the standard space of Talairach and Tournoux
(13)
using the subject's coregistered structural
Tl
scan (14,
15).
The data were spatially smoothed with a 5-mm iso-
tropic Gaussian kernel and temporally smoothed with a
,8
-s
Gaussian kernel. Because we also smoothed the de-
sign matrix, the temporal smoothing does not affect the
kernel or response function estimates
(2).
Epoch-Related Responses
The data were analyzed using a design matrix that in-
cluded the explanatory variables (convolved time series)
in Eq.
[3].
The basis functions employed in this analysis
were a series of gamma density functions as shown in
Fig.
1
(solid lines). The stimulus function
u(t),
the sup-
posed neuronal activity, was simply the word presenta-
tion rate at which the scan was acquired. We also used
more comprehensive forms for
u(t)
that involved model-
ing each word individually, but the results were very
similar to the simpler snalysis presented here. The re-
sulting SPM(FJ, reflecting the significance of an evoked
response (or more formally, testing the null hypothesis
that all
h1
and
h2
were jointly zero), is shown in Fig.
2
44
Friston
et
al.
Design Matrix
SPM{F}
-
Rate experiment
SPM{F}
-
Event-related
200
400
600
800
1000
1200
10
20
30
40
50
Design Matrix
100
200
300
400
500
600
700
800
900
1000
along with the design matrix
used. The left-hand side of the
upper design matrix com-
prises the explanatory variable
xi(t)
and xi(t).xi(t). The remain-
ing columns contain the con-
stant (used to estimate
go)
and
other effects designated as
confounds (low-frequency ar-
tifacts, global effects, and
so
on). This
SPM(F}
has been
thresholded
[F
=
32,
P
<
0.001
corrected for multiple compar-
isons (IS)] and shows wide-
spread responses in bilateral
temporal regions with the
most significant effects evi-
dent in the periauditory re-
gions.
The estimated kernels
ho,
h’,
and
h2
for a voxel in the
left superior temporal gyrus
(-56, -28,12 mm) are shown in
Fig.
3.
As
might be expected,
the first-order kernel resem-
bles the hemodynamic re-
sponse functions identified
using linear analyses such as
least squares deconvolution
[e.g., Fig. 6 in Ref.
(l)]
or
linear
regression
[
e.g., Fig.
5
in Boyn-
ton
et
al.
(6)]. Of note is the
protracted undershoot that
lasts for about 16
s.
The sec-
ond-order kernel is shown be-
low and is remarkable for the
pronounced negativity on the
lower left, flanked by smaller
positive lobes. This negativity
suggests that
if
neuronal activ-
ity has been high in the past
few seconds, then the hemo-
20
40
60
dynamic response will be sup-
FIG.
2.
Top left: SPM{F} testing for the significance of the
first-
and second-order kernel coeffi-
cients
(b’
and
b2)
in
the
first
(rate) experiment.
This
is
a maximum intensity projection
of
a statistical
process of the
F
ratio, following a multiple regression analysis at each voxel. The format
is
standard
and provides three orthogonal projections
in
the standard space conforming
to
that described
in
Talairach and Tournoux
(13).
The grey scale
is
arbitrary, and the
SPM{F)
has been thresholded at
32
(f
<
0.001 corrected). Top
right:
The design matrix used
in
the analysis. The design matrix
comprises the explanatory variables
in
the general linear model.
It
has one row for each of the
1200
scans and one column for each explanatory variable
or
effect modeled. The left-hand columns
contain the explanatory variables of interest
x,(t)
and
x,(t).x,(t),
where
x,(t)
is
word presentation rate
u(t)
convolved
with
the
basis functions
b,(t)
in
Fig.
1.
The remaining columns contain covariates or
effects of no interest designated as confounds. These include (left to right) a constant
term
(b’),
periodic (discrete cosine set) functions of time, to remove low-frequency artifacts and
drifts,
global
or whole brain activity
G(t),
and interactions between global effects and those of interest
G(t).x,(t)
and
G(t).x,(t).x,(t).
The latter confounds remove effects that have no regional specificity. Lower left:
SPM{F) as above
but
for the second event-related experiment.
This
SPM{F} has been thresholded
at
F
=
16
(f
<
0.001
corrected). Note the similarity between the two SPM{F}s, despite the fact that
they
derive from different experimental designs and completely independent data. Lower right:
The
design matrix employed. The effects of interest are the first-order terms
x,(t)
derived
by
convolving
u(t)
with
all
six
functions
in
Fig.
1.
u(t)
was
in
this
instance one when a single word was presented
and
zero elsewhere. The confounds are as described above.
pressed. The positive lobes
suggest that this suppression
is ameliorated
if
the underly-
ing neuronal activity is
sus-
tained, i.e., is high in the re-
cent
(4
s)
and more distant
past
(8
s).
There are two fur-
ther important points to note.
First, the second-order kernel
is symmetrical. This will al-
ways
be
the case because the
contribution of
u(t
-
Tl).u(f
-
T~)
to the response is exactly the
same as
u(t
-
T~).u(~
-
T~).
The
second unpredicted and more
intriguing observation
is
that
the second-order kernel is very
similar to the “product”
of
the
first-order kernel times itself
Nonlinear Event-Related Responses
0th
order kernel
(hO)
=
146
1st
order kernel
{hi}
time
(seconds}
FIG.
3.
The Volterra kernels
ho,
h’,
and
h2
based on parameter
estimates from a voxel in the left superior temporal gyrus at
-56,
-28, 12 mm. These kernels can be thought of as a characterization
of the second-order hemodynamic response function. The first-
order kernel (upper panel) represents the (first-order) component
usually presented in linear analyses. The second-order kernel
(lower panel)
is
presented in image format. The color scale is
arbitrary; white
is
positive and black is negative. The insert on the
right represents
[-h1(T,).h1(T2)],
the second-order kernel that would
be predicted by a simple model that involved convolution with
h1
and then some nonlinear scalar function.
(insert in Fig.
3).
In other words,
h2(~l,~2)
is roughly
proportional to
hl(i-J.hl(~J.
We return to the implica-
tions of this in a subsequent section.
The basic form of the second-order kernel estimates
was very similar for all the voxels in the periauditory
region. The nature
of
these nonlinear effects will be dem-
onstrated more intuitively in the final section of this
paper. Using these kernel estimates, we can estimate
responses to any temporal pattern of words presented by
using
Eq.
[I]
and any suitable function
u(t).
In the first
instance, we present the estimated responses to the stim-
uli actually used: The predicted and observed responses
for the first
17
min of the time series are shown in Fig.
4
(upper panel). Predicted responses to
34-s
epochs at two
presentation rates
(30
and 60 words/min) and the ad-
justed responses observed
are
shown in the lower panel
45
of Fig.
4
(adjusted data is simply the original data after
the confounding effects have been removed). The agree-
ment is evident.
Event-Related Responses
By specifying
a
stimulus function
u(t)
that models the
occurrence of
a
single word, we can use
Eq.
[I] and the
kernel estimates in Fig.
3
to simulate the hemodynamic
response
of
this brain region to single word. This simu-
lated event-related response is shown in Fig.
5
(upper
panel). One observes a peak at about
4
s
followed by a
protracted undershoot lasting for about 16
s.
This re-
sponse
is
“simulated” using a model whose parameters
were determined without ever presenting single words in
isolation
(i.e,,
the Volterra series model based on the rate
experiment).
A
validation of the model can be effected in
terms of the empirically determined event-related re-
sponse
to
actual single words using the second experi-
ment.
The same analysis described above was applied to the
event-related, single word experiment. In this instance,
by virtue of the fact that the words were presented very
sparsely, there is no opportunity for the responses to
fitted and observed response
I
143
0
200
400
600
aoo
time
(seconds)
responses at
30
and 60 wpm
143
i
0
10
20
30 40
50
60
time
(seconds)
FIG.
4.
Top panel: Fitted or predicted (broken line) and observed
(solid line) responses at the same voxel as in Fig.
3,
over the first
1024
s.
The observed responses here are adjusted such that those
effects that can be modeled by the confounds have been re-
moved. Lower panel: Predicted and observed responses for ep-
ochs
of
30
(light grey) and
60
(dark grey) wordshin (wpm) super-
imposed upon each other.