METHODS
published: 30 November 2016
doi: 10.3389/fnhum.2016.00604
Frontiers in Human Neuroscience | www.frontiersin.org 1 November 2016 | Volume 10 | Article 604
Edited by:
Vladimir Litvak,
UCL Institute of Neurology, UK
Reviewed by:
Paul Fredrick Sowman,
Macquarie University, Australia
Anna Jafarpour,
University of California, Berkeley, USA
*Correspondence:
Edmund C. Lalor
edmund_lalor@urmc.rochester.edu
Received: 13 July 2016
Accepted: 11 November 2016
Published: 30 November 2016
Citation:
Crosse MJ, Di Liberto GM, Bednar A
and Lalor EC (2016) The Multivariate
Temporal Response Function (mTRF)
Toolbox: A MATLAB Toolbox for
Relating Neural Signals to Continuous
Stimuli. Front. Hum. Neurosci. 10:604.
doi: 10.3389/fnhum.2016.00604
The Multivariate Temporal Response
Function (mTRF) Toolbox: A MATLAB
Toolbox for Relating Neural Signals
to Continuous Stimuli
Michael J. Crosse
1, 2
, Giovanni M. Di Liberto
1
, Adam Bednar
1, 3
and Edmund C. Lalor
1, 3
*
1
School of Engineering, Trinity Centre for Bioengineering and Trinity College Institute of Neuroscience, Trinity College Dublin,
Dublin, Ireland,
2
Department of Pediatrics and Department of Neuroscience, Albert Einstein College of Medicine, The Bronx,
NY, USA,
3
Department of Biomedical Engineering and Department of Neuroscience, University of Rochester, Rochester, NY,
USA
Understanding how brains process sensory signals in natural environments is one of
the key goals of twenty-first century neuroscience. While brain imaging and invasive
electrophysiology will play key roles in this endeavor, there is also an important role to
be played by noninvasive, macroscopic techniques with high temporal resolution such
as electro- and magnetoencephalography. But challenges exist in determining how best
to analyze such complex, time-varying neural responses to complex, time-varying and
multivariate natural sensory stimuli. There has been a long history of applying system
identification techniques to relate the firing activity of neurons to complex sensory stimuli
and such techniques are now seeing increased application to EEG and MEG data.
One particular example involves fitting a filter—often referred to as a temporal response
function—that describes a mapping between some feature(s) of a sensory stimulus and
the neural response. Here, we first briefly review the history of these system identification
approaches and describe a specific technique for deriving temporal response functions
known as regularized linear regression. We then introduce a new open-source toolbox
for performing this analysis. We describe how it can be used to derive (multivariate)
temporal response functions describing a mapping between stimulus and response in
both directions. We also explain the importance of regularizing the analysis and how this
regularization can be optimized for a particular dataset. We then outline specifically how
the toolbox implements these analyses and provide several examples of the types of
results that the toolbox can produce. Finally, we consider some of the limitations of the
toolbox and opportunities for future development and application.
Keywords: system identification, reverse correlation, stimulus reconstruction, sensory processing, EEG/MEG
INTRODUCTION
Traditionally, research on the electrophysiology of sensory processing in humans has focused
on the rather special case of brief, isolated stimuli because of the need to time-lock to discrete
sensory events in order to estimate event-related potentials (ERPs;
Handy, 2005; Luck, 2014). The
objective is to estimate the impulse response function of the sensory system under investigation by
Crosse et al. The mTRF Toolbox
convolving the system with a transient, impulse-like stimulus
and averaging over several-hundred time-locked response trials.
This approach has been used extensively to study how the
human brain processes various ecological events, even those that
occur in a continuous, dynamic manner such as human speech
(e.g., Salmelin, 2007; Picton, 2013). However, the type of speech
stimuli used in such ERP studies usually consist of individual
phonemes or syllables and are therefore not entirely reflective
of natural, connected speech which is ongoing and abundant
with lexical complexity. Recent studies have begun to use more
naturalistic, extended speech stimuli by focusing their analysis
on me asuring the phase of neural responses across multiple
repetitions of the same speech segment (Luo and Poeppel, 2007;
Zion-Golumbic et al., 2013). While this approach has revealed
novel and import ant insig hts into the neurophysiology of speech
processing, it does not facilitate characterization of the system’s
response function, and in any case, is an indirect measure of how
the brain entrains to the stimulus over time.
A more direct way to investigate neural entrainment to
continuous stimuli is to mathematically model a function that
describes the way a particular property of the stimulus is
mapped onto neural responses, a technique known as system
identification (SI;
Marmarelis, 2004). While there are several
classes of models that can be implemented for this purpose
(reviewed in Wu et al., 2006), the most straightforward class are
linear time-invariant (LTI) systems. Although the human brain
is neither linear nor time-invariant, these assumptions can be
reasonable in certain cases (e.g., Boynton et al., 1996) and allow
for the syst em to be characterized by its impulse response. An SI
method known as “reverse correlation” has become a common
technique for characterizing LTI systems in neurophysiology
(
Ringach and Shapley, 2004), an approach that has long been
established in both visual and auditory animal electrophysiology
(De Boer and Kuyper, 1968; Marmarelis and Marmarelis,
1978; Coppola, 1979). This technique approximates the impulse
response of the sensory system under investigation, except it does
not require the use of discrete stimuli. While this is somewhat
analogous to calculating an ERP, there are important differences
that must be considered: (1) the response function obtained
by reverse correlation only reflects the response of the system
to specific stimulus parameters defined by the experimenter as
opposed to the entire event, (2) reverse correlation makes the
assumption that the input-output relationship of the system is
linear, unlike time-locked averaging and (3) reverse correlation
converges on a more temporally precise estimate of the systems
impulse response than an ERP (which is susceptible to temporal
smearing). Reverse correlation in its simplest form can be
implemented via a straightforward cross-correlation between the
input and output of an LTI system (Ringach and Shapley, 2004).
While this approach has been used to study how speech is
encoded in human brain activity (Ahissar et al., 2001; Abrams
et al., 2008; Aiken and Picton, 2008), it is better suited to stimuli
modulated by a stochastic process such as Gaussian white noise.
As such, most instances of this approach in animal models have
traditionally used white noise stimuli (
De Boer and Kuyper, 1968;
Marmarelis and Marmarelis, 1978; Coppola, 1979; Eggermont
et al., 1983; Ringach et al., 1997
). This work has even inspired
researchers to investigate how such stochastic signals are encoded
in th e human brain (
Lalor et al., 2006, 2009).
That s aid, the human brain has evolved to process ecologically
relevant stimuli that rarely conform to a white random process.
For example, in the context of human neuroscience research,
a proper understanding of how the brain processes natural
speech would surely require that natural speech is used as a
stimulus in the laboratory, given that neurons respond differently
to more complex stimuli (
Theunissen et al., 2000). As such,
researchers using animal models have shifted their focus toward
studying the brain using more naturalistic stimuli thanks to
the development of SI methods such as “normalized reverse
correlation” (NRC; Theunissen et al., 2001), “ridge regression”
(Machens et al., 2004), and “boosting” (David et al., 2007). Each
of these techniques converge on the same theoretical solution
but use different priors and, critically, give an unbiased impulse
response estimate for non-white stimuli. This has inspired
researchers to characterize the “spectrotemporal receptive fields”
of auditory cortical neurons in various animal models (
Depireux
et al., 2001; Tomita and Eggermont, 2005). As a result, researchers
interested in how human speech is processed have begun
to model response functions describing the linear mapping
between properties of natural speech (such as the envelope or
spectrogram) and population responses in bot h animals (David
et al., 2007; Mesgarani et al., 2008) and humans (Lalor and
Foxe, 2010; Ding and Simon, 2012b). There have been similar
efforts to model response functions relating more natural visual
stimulus properties such as motion to neural responses in
humans (Gonçalves et al., 2014), again inspired by previous
single-unit electrophysiology work (
Jones and Palmer, 1987;
David and Gallant, 2005).
Most of the aforementioned studies have modeled the
stimulus-response mapping function in the forward direction
(i.e., forward modeling). However, this mapping can also be
modeled in the reverse direction (i.e., backward modeling),
offering a complementary way to investigate how stimulus
features are encoded in neural response measures. Unlike
forward models, backward model parameters are not readily
neurophysiologically interpretable (see Haufe et al., 2014),
but can be used to reconstruct or decode stimulus features
from the neural response, a method known as “stimulus
reconstruction.” This approach has several advantages over
forward modeling approaches, especially when recording from
population responses using multi-channel systems such as EEG.
Firstly, because reconstruction project s back to the stimulus
domain, it does not require pre-selection of neural response
channels (
Mesgarani et al., 2009). In fact, inclusion of all response
channels in the backward model is advantageous because the
reconstruction method gives a low weighting to irrelevant
channels whilst allowing the model to capture additional
variance using channels potentially excluded by feature selection
approaches (Pasley et al., 2012). Secondly, backward modeling
can offer increased sensitivity to important signal differences
between response channels that are highly correlated with each
other (as is often the case with EEG). It can do this because the
analysis maps the data from all response channels simultaneously
(i.e., in a multivariate manner) and so it can recognize any
Frontiers in Human Neuroscience | www.frontiersin.org 2 November 2016 | Volume 10 | Article 604
Crosse et al. The mTRF Toolbox
inter-channel correlation in the data (Mesgarani et al., 2009).
In contrast, when performing forward modeling, each analysis
is univariate and thus is ignorant of the data on the other
EEG channels. Thirdly, stimulus features that are not explicitly
encoded in the neural response may be inferred from correlated
input features that are encoded. This prevents the model from
allocating resources to the encoding of redundant stimulus
information (
Barlow, 1972). The stimulus reconstruction method
has previously been used to study both the visual and auditory
system in various a nimal models (Bialek et al., 1991; Rieke et al.,
1995; Stanley et al., 1999). More recently, it has been adopted for
studying speech processing in the human brain using intracranial
and non-invasive electrophysiology (Mesgarani et a l., 2009;
Pasley et al., 2012; Ding and Simon, 2013; Martin et al., 2014;
Crosse et al., 2015a, 2016; O’Sullivan et al., 2015).
While certain research groups now regularly use SI to study
sensory processing in the human brain, the approach has perhaps
not yet been as widely adopted throughout the neuroscience
community as it might because of the (at least perceived)
challenges associated with its implementation. The goal of the
present paper is to introduce a recently-developed SI toolbox
that provides a straightforward and flexible implementation
of regularized linear (ridge) regression (
Machens et al., 2004;
Lalor et al., 2006). We begin by summarizing the mat h ematic s
underlying this technique, continue by providing some concrete
examples of how the toolbox can be used and conclude by
discussing some of its applications and important considerations.
REGULARIZED LINEAR REGRESSION
Forward Models: Temporal Response
Function Estimation
Forward models are sometimes referred to as generative or
encoding models because they describe how the system generates
or encodes information (
Haufe et al., 2014). Here, they will be
referred to as temporal response functions (TRFs; Ding and
Simon, 2012b). There are a number of ways of mathematically
describing how the input to a system relates to its output. One
commonly used approach—and the one that will be described
in this paper—is to assume that the output of the system is
related to the input via a simple linear convolution. In the
context of a sensory system where the output is monitored by
N recording channels, let’s assume t hat the instantaneous neural
response r(t, n), sampled at times t = 1...T and at channel n,
consists of a convolution of the stimulus property, s(t), with an
unknown channel-specific TRF, w(τ , n). The response model can
be represented in discrete time as:
r(t, n) =
X
τ
w(τ , n)s(t − τ ) + ε(t, n), (1)
where ε(t, n) is the residual response at each channel not
explained by the model. Essentially, a TRF can be thought
of as a filter that describes the linear transformation of the
ongoing stimulus to the ongoing neural response. The TRF,
w(τ , n), describes this transformation for a specified range of time
lags, τ , relative to the instantaneous occurrence of the stimulus
feature, s(t).
In the context of speech for example, s(t) could be a measure
of t h e spee ch envelope at each moment in time and r(t, n) could
be t h e corresponding EEG response at channel n. The range of
time lags over which to calculate w(τ , n) migh t be that typically
used to capture the cortical response components of an ERP,
e.g., −100–400 ms. The resulting value of the TRF at −100 ms,
would index the relationship between the speech envelope and
the neural response 100 ms earlier (obviously this should have an
amplitude of zero), whereas the TRF at 1 00 ms would index how a
unit change in the amplitude of the speech envelope would affect
the EEG 100 ms later (
Lalor et al., 2009).
The TRF, w(τ , n), is estimated by minimizing the mean-
squared error (MSE) between the actual neural response, r(t, n),
and that predicted by the convolution, ˆr(t, n):
min ε(t, n) =
X
t
r(t, n) − ˆr(t, n)
2
. (2)
In practice, this is solved using reverse correlation (
De Boer
and Kuyper, 1968), which can be easily implemented using the
following matrix operations:
w =
S
T
S
−1
S
T
r, (3)
where S is the lagged time series of the stimulus property, s, and
is defined as follows:
S =
s(1 − τ
min
) s(−τ
min
) · · · s(1) 0 · · · 0
.
.
.
.
.
. · · ·
.
.
. s(1) · · ·
.
.
.
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · · 0
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · · s(1)
s(T)
.
.
. · · ·
.
.
.
.
.
. · · ·
.
.
.
0 s(T) · · ·
.
.
.
.
.
. · · ·
.
.
.
.
.
. 0 · · ·
.
.
.
.
.
. · · ·
.
.
.
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · ·
.
.
.
0 0 · · · s(T) s(T − 1) · · · s(T − τ
max
)
.
(4)
The values τ
min
and τ
max
represent the minimum and maximum
time lags (in samples) respectively. In S, each time lag is arranged
column-wise and non-zero lags are padded with zeros to ensure
causality (
Mesgarani et al., 2009). The window over which the
TRF is calculated is defined as τ
window
= τ
max
− τ
min
and the
dimensions of S are thus T × τ
window
. To include the constant
term (y-intercept) in the regression model, a column of ones i s
concatenated to th e left of S. In Equation (3 ), variable r is a matrix
containing all the neural response data, with channels arranged
column-wise (i.e., a T × N matrix). The resulting TRF, w, is a
τ
window
× N matrix with each column representing the univariate
mapping from s to the neural response at each channel.
Frontiers in Human Neuroscience | www.frontiersin.org 3 November 2016 | Volume 10 | Article 604
Crosse et al. The mTRF Toolbox
One of the important points here is that t his analysis explicitly
takes into account the autocovariance structure of the stimulus.
In non-white stimuli, such as natural speech, the intensity of
the acoustic signal modulates gradually over time, meaning it
is correlated with itself at non-zero time lags. A simple cross-
correlation of a speech envelope and the corresponding neural
response would result in temporal smearing of the impulse
response function. The solution here (Equation 3) is to divide out
the autocovariance structure of the stimulus from the model such
that it removes the correlation between different time points.
The TRF approach, which does this, is therefore less prone to
temporal smearing than a simple cross-correlation approach.
This is demonstrated in a worked example in the next section.
Regularization
An important consideration when calculating the TRF is that
of regularization, i.e., introducing additional information to
solve any ill-posed estimation problems and prevent overfitting.
The ill-posed estimation problem has to do with inverting the
autocovariance matrix, S
T
S. Matrix inversion is particularly
prone to numerical instability when solved with finite precision.
In other words, small changes in S
T
S (such as rounding errors
due to discretization) could cause large changes in w if the former
is ill-conditioned. In other words, the estimate of w can have very
high variance. This does not usually apply when the stimulus
represents a stochastic process because S
T
S would be full rank
(
Lalor et al., 2006). However, the autocorrelation properties of a
non-white stimulus such as speech means t hat it is more likely to
be singular (i.e., have a determinant of zero). Typically, numerical
treatment of an ill-conditioned matrix involves reducing the
variance of the estimate by adding a bias term or “smoothing
solution.” Specifically, because the overall estimation error is
made up of both a bias term (i.e., the difference between the
estimate’s expected value and its true value) and a variance term,
one can deliberately increase the bias so as to reduce the (high)
variance of the estimate by so much as to decrease the overall
estimation error.
Addition of this smoothing term also solves the other main
issue, that of overfitting. The re verse correlation analysis is
utterly agnostic as to the biological nature of the data that it is
being asked to model. As a result, without re gularizat ion, the
resulting TRF will be optimal i n terms of the particular fitting
criterion (e.g., least squares error) for the specific dataset that
was used for the fitting. And, given that those data will be
“noisy,” the TRF can display biologically implausible properties
such as very high-frequency fluctuations. Using this TRF to then
predict unseen data will likely result in suboptimal performance,
because the high frequency fluctuations will not necessarily
correspond well to the “noise” in the new data. In ot h er words,
the TRF has been “overfit” to the specific dataset used in the
training. Regularization serves to prevent overfitting to such
high-frequency, dataset-specific noise along the low-variance
dimensions (
Theunissen et al., 2001; Mesgarani et al., 2008). It
can do this, for example, by penalizing large differences between
neighboring TRF values, thereby forcing the TRF to be smoother.
This makes the TRF less specific to t h e data that was used to fit it
and can help it generalize better to new unseen data.
In practice, both ill-posed problems and overfitting can be
solved simultaneously by weighting the diagonal of S
T
S before
inversion, a method known as Tikhonov regularization or ridge
regression (
Tikhonov and Arsenin, 1977):
w =
S
T
S + λI
−1
S
T
r, (5)
where I is the identity matrix and λ is the smoothing constant
or “ridge parameter.” The ridge parameter can be adjusted using
cross-validation to maximize the correlation between r(t, n),
and ˆr(t, n) (
David and Gallant, 2005). TRF optimization will be
described in more detail in the next section. While this form
of ridge re gression enforces a smoothness constraint on t h e
resulting model by penalizing TRF values as a function of their
distance from zero, another option is to quadratically penalize
the difference between each two neighboring terms of w (
Lalor
et al., 2006):
w =
S
T
S + λM
−1
S
T
r ,where M =
1 −1
−1 2 −1
−1 2 −1
.
.
.
.
.
.
.
.
.
−1 2 −1
−1 1
.
(6)
Tikhonov regularization (Equation 5) reduces overfitting by
smoothing the TRF estimate in a way that is insensitive to
the amplitude of the signal of interest. However, the quadratic
approach (Equation 6) reduces off-sample error whilst preserving
signal amplitude (
Lalor et al., 2006). As a result, this approach
usually leads to an improved estimate of the system’s response (as
indexed by MSE) compared to Tikhonov regularization.
Multivariate Analysis
The previous section focused on the specific case of relating
a single, univariate input stimulus feature (e.g., t h e envelope
of a speech stimulus) separately to each of multiple recording
channels. However, most complex stimuli in nature are not
processed as simple univariate features. For example, when
auditory speech enters the ear, the signal is transformed into a
spectrogram representation by the cochlea, consisting of multiple
frequency bands which project along the auditory pathway (Yang
et al., 1992). The auditory system maps each of t h ese frequency
bands to t he neural representation measured at the cortical level.
This process can be modeled by a multivariate form of the TRF
(i.e., mTRF).
Indeed, it is possible to define an mTRF that linearly
maps a multivariate stimulus feature to each recording channel
(Theunissen et al., 2000; Depireux et al., 2001). Using the above
example, let s(t, f ) represent the spectrogram of a speech signal
at frequency band f = 1...F. To derive the mTRF, the stimulus
lag matrix, S (Equation 4 ), is simply extended such that every
column is replaced with F columns, each representing a different
frequency band (i.e., a T × Fτ
window
matrix). The resulting
mTRF, w(f , τ , n), will be a Fτ
window
× N matrix but can easily be
Frontiers in Human Neuroscience | www.frontiersin.org 4 November 2016 | Volume 10 | Article 604
Crosse et al. The mTRF Toolbox
unwrapped such that each independent variable is represented as
a separate dimension (i.e., a F × τ
window
× N matrix). Here, the
constant term is included by concatenating F columns to the left
of S.
An important consideration in multivariate TRF analysis
is which method of regularization to use. The quadratic
regularization term in Equation (6) was designed to enforce
a smoothness constraint and maintain SNR along the time
dimension, but not any other. For high λ values, this approach
would cause smearing across frequencies; hence it would not
yield an accurate representation of the TRF in each frequency
band. In this case, it will typically be most a ppropriate to use
the identity matrix for regularization (Equation 5) so as to
avoid enforcing a smoothness constraint across the non-time
dimension of the mTRF—although, in some cases, this may
actually be what is desired.
Backward Models: Stimulus
Reconstruction
The previous sections describe how to forward model the linear
mapping between the stimulus and the neural response. While
this approach can be extended to accommodate multivariate
stimulus features, it is suboptimal in the sense that it treats
each neural response channel as an independent univariate
feature. Backward modeling, on the other hand, derives a re verse
stimulus-response mapping by exploiting all of the available
neural data in a multivariate context. Backward models are
sometimes referred to as discriminative or decoding models,
because they attempt to reverse the data generating process by
decoding the stimulus features from the neural response (
Haufe
et al., 2014). Here, they will simply be referred to as de coders.
Decoders can be modeled in much the same way as TRFs.
Suppose the decoder, g(τ , n), represents the linear mapping from
the neural response, r(t, n), back to the stimulus, s(t). This could
be expressed in discrete time as:
ˆs(t) =
X
n
X
τ
r(t + τ , n)g(τ , n), (7)
where ˆs(t) is the reconstructed stimulus property. Here, the
decoder integrates the neural response over a specified range of
time lags τ . Ideally, these lags will capture the window of neural
data that optimizes reconstruction of the stimulus property.
Typically, the most informative lags for reconstruction are
commensurate with those used to capture the major components
of a forward TRF, except in the reverse direction as the decoder
effectively maps backwards in time. To reverse the lags used in
the earlier TRF example (τ
min
= −100 ms, τ
max
= 400 ms),
the values of τ
min
and τ
max
are swapped but their signs remain
unchanged, i.e., τ
min
= −400 ms, τ
max
= 100 ms.
The decoder, g(τ , n), is estimated by minimizing the MSE
between s(t) and ˆs(t):
min ε(t) =
X
t
s(t) − ˆs(t)
2
. (8)
Analogous to the TRF approach, the decoder is computed using
the following matrix operations:
g =
R
T
R + λI
−1
R
T
s (9)
where R is the lagged time series of the response matrix, r. For
simplicity, we will define R for a single-channel response system:
R =
r(1 − τ
min
,1) r(−τ
min
,1) · · · r(1,1) 0 · · · 0
.
.
.
.
.
. · · ·
.
.
. r(1, 1) · · ·
.
.
.
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · · 0
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · · r(1,1)
r(T,1)
.
.
. · · ·
.
.
.
.
.
. · · ·
.
.
.
0 r(T,1) · · ·
.
.
.
.
.
. · · ·
.
.
.
.
.
. 0 · · ·
.
.
.
.
.
. · · ·
.
.
.
.
.
.
.
.
. · · ·
.
.
.
.
.
. · · ·
.
.
.
0 0 · · · r(T,1) r(T − 1,1) · · · r(T − τ
max
,1)
,
(10)
As before, t h is can be extended to the multivariate case of an N-
channel system by replacing each column of R with N columns
(each representing a separate recording channel). For N channels,
the dimensions of R would be T × Nτ
window
. The constant term
is included by concatenating N columns of ones to the left of
R. In the context of speech, the sti mulus variable, s, represents
either a column-wise vector (e.g., envelope) or a T × F matrix
(e.g., spectrogram). The resulting decoder, g, would be a vector of
Nτ
window
samples or a Nτ
window
× F matrix, respecti vely. While
interpretation of decoder weights is not as straight forward as
that of a TRF, one may wish to separate it s dimensions (e.g.,
N×τ
window
×F) to examine the relative weighting of each channel
at a spec ific time lag. The channel weights represent the amount
of information that each channel provides for reconstruction, i.e.,
highly informative channels rec eive weights of greater magnitude
while channels providing little or no information receive weights
closer to zero.
In Equation (9), Tikhonov regularization is used as it is
assumed that the neural response data i s multivariate. As
mentioned above, any bias from the correlation between the
neural response channels is removed in the reconstruction
approach. In practice, this is achieved by dividing out the
autocovariance structure of the neural response (see Equation
9). As a result, channel weighting becomes much more localized
because inter-channel redundancies are no longer encoded in the
model, giving it an advantage over the forward TRF method and
cross-correlation approaches.
MTRF TOOLBOX: IMPLEMENTATION AND
FUNCTIONALITY
This section outlines how regularized linear regression can be
implemented in MATLAB using the mTRF Toolbox (https://
Frontiers in Human Neuroscience | www.frontiersin.org 5 November 2016 | Volume 10 | Article 604