# Online asynchronous detection of error-related potentials in participants
with spinal cord injury by adapting a pre-trained generic classifier

Catarina Lopes-Dias1, Andreea I. Sburlea1, Katharina Breitegger2, Daniela
Wyss2, Harald Drescher2, Renate Wildburger2 and Gernot R. Müller-Putz1,3 1
Institute of Neural Engineering, Graz University of Technology, Graz, Austria
2 AUVA-Rehabilitationsklinik Tobelbad, Tobelbad, Austria  3 BioTechMed, Graz,
Austria<EMAIL_ADDRESS>
###### Abstract

A brain-computer interface (BCI) user awareness of an error is associated with
a cortical signature named error-related potential (ErrP). The incorporation
of ErrPs’ detection in BCIs can improve BCIs’ performance.

Objective: This work is three-folded. First, we investigate if an ErrP
classifier is transferable from able-bodied participants to participants with
spinal cord injury (SCI). Second, we test this generic ErrP classifier with
SCI and control participants, in an online experiment without offline
calibration. Third, we investigate the morphology of ErrPs in both groups of
participants.

Approach: We used previously recorded electroencephalographic (EEG) data from
able-bodied participants to train an ErrP classifier. We tested the classifier
asynchronously, in an online experiment with 16 new participants: 8
participants with SCI and 8 able-bodied control participants. The experiment
had no offline calibration and participants received feedback regarding the
ErrPs’ detection from its start.  For a matter of fluidity of the experiment,
the feedback regarding false positive ErrP detections was not presented to the
participants but these detections were taken into account in the evaluation of
the classifier. The generic classifier was not trained with the user’s brain
signals. Still, its performance was optimized during the online experiment
with the use of personalized decision thresholds. The classifier’s performance
was evaluated using trial-based metrics, which consider the asynchronous
detection of ErrPs during the entire trials’ duration.

Main results: Participants with SCI presented a non-homogenous ErrP
morphology, and four of them did not present clear ErrP signals. The generic
classifier performed above chance level in participants with clear ErrP
signals, independently of the SCI (11 out of 16 participants). Three out of
the five participants that obtained chance level results with the generic
classifier would have not benefited from the use of a personalized classifier.

Significance: This work shows the feasibility of transferring an ErrP
classifier from able-bodied participants to participants with SCI, for
asynchronous detection of ErrPs in an online experiment without offline
calibration, which provided immediate feedback to the users.

Keywords:  error-related potential, asynchronous classification, generic
classifier, online experiment, spinal cord injury, end-users, brain-computer
interface

##  1 Introduction

Brain-computer interfaces (BCIs) can assist people with severe motor
impairments to operate external devices by converting their modulated brain
activity into the control of these devices [1, 2, 3]. Although being a
promising technology, most BCIs are still error-prone, and the frequent
occurrence of errors can spoil the experience of the BCI user. The user’s
awareness of an unintended response from the device that he/she is controlling
is associated with a neural signature known as error-related potential (ErrP)
[4, 5].

ErrPs are associated with conflict monitoring and error processing [6] and can
be measured using non-invasive techniques, such electroencephalography (EEG),
that are often used for BCIs’ control. Therefore, ErrPs can be used to improve
BCIs’ performance either in a corrective manner, by allowing corrective
actions, or in an adaptive manner, by reducing the possibility of future
errors [7, 8, 9, 10, 11]. The real-time detection of ErrPs is pertinent in
BCIs used by persons with motor impairments and also in applications targeting
healthy users [12, 13, 14, 15]. The incorporation of ErrPs’ detection in a BCI
promotes a smoother interaction with its user. Nevertheless, this
incorporation is not widely investigated.

The use of ErrPs in discrete BCIs, which are controlled in discrete steps, is
well established in healthy participants [4, 16, 17, 18, 19, 20, 21, 22, 10]
and has also been marginally tested in potential end-users of BCIs [23].
Still, BCIs are developing in the direction of offering users continuous
control of an external device - continuous BCIs [24, 25, 26, 27, 28]. The
incorporation of ErrPs in such BCIs requires an asynchronous detection of
ErrPs, since the user can realise at any moment, during the control of the
device, that an error has occurred. The asynchronous detection of ErrPs has
been studied in healthy participants, both in offline scenarios [29, 30, 31,
32, 13, 33, 34] and more recently in online scenarios [35].

A possible explanation for the limited use of ErrPs in BCIs can be linked with
most BCIs relying on personalised classifiers, which are constructed with the
user’s brain signals. Since a considerable amount of data is necessary to
reliably train the classifier, personalised classifiers commonly require a
long calibration period before the user can receive feedback of its own brain
signals. In this manner, combining ErrPs with other controlling signals would
imply collecting calibration data for all the different signals, increasing
even more the calibration period. Alternatively, using an ErrP classifier that
would not require calibration with the user could encourage the integration of
ErrPs with other control signals when constructing BCIs. This could be
achieved by either transferring an ErrP classifier across different tasks or
across different participants. Both options have been tested in discrete
tasks, in offline conditions [36, 23, 37, 38, 39, 40, 41, 42, 43] and in
online conditions [20]. Recently, the asynchronous detection of ErrPs with a
generic classifier has been studied in the context of a continuous task, in
offline conditions [44] and in pseudo-online conditions [45].

Very few works addressed the study of ErrPs in potential BCI end-users and the
existing studies are mainly conducted offline. Keyl and colleagues
characterized the morphology of ErrPs of spinal cord injured participants and
compared it with able-bodied control participants [46]. The ErrP morphology
was comparable in the two groups but the ErrPs of the participants in the SCI
group showed smaller peak amplitudes. Kumar and colleagues studied ErrPs
during post-stroke rehabilitation movements[47]. In this work, individual
participants did not display very clear ErrP patterns. Spüler and colleagues
studied ErrPs in six participants with amyotrophic lateral sclerosis (ALS) in
an online experiment, and showed that the incorporation of ErrPs improves the
BCI performance [23]. This work also analysed, offline, the transfer of an
ErrP classifier from ALS participants to able-bodied control participants.

Our study has three main aims. First, we test the feasibility of transferring
an ErrP classifier for asynchronous classification from able-bodied
participants to potential end-users of BCIs, in particular participants with a
high spinal cord injury (SCI). Second, we test the feasibility of using a
generic ErrP classifier asynchronously in an online experiment in which both
participants with SCI and control participants took part. Third, we
investigate the morphology of ErrPs both in participants with SCI and in
control participants.

In the work presented here, we recorded EEG from both participants with SCI
and control participants while testing asynchronously a generic ErrP
classifier in a closed-loop online experiment. The generic classifier has been
trained with the EEG data of 15 able-bodied participants from a previous study
of ours and was not retrained during the experiment [35]. This allowed us to
create an online experiment with no offline calibration period, in which
participants received immediate feedback of their brain signals from the very
beginning of the experiment onwards.

##  2 Methods

###  2.1 Participants

Sixteen volunteers participated in the experiment, eight of which had a spinal
cord injury. The age of the participants with SCI was $37.5\pm 9.7$ years
(mean $\pm$ std). The remaining participants were able-bodied control
participants. Each participant with SCI was matched with a control participant
of the same sex and a maximum age difference of 5 years. The control
participants were $35.9\pm 10.8$ years old (mean $\pm$ std).

All participants with SCI had a spinal cord injury between levels C4 and Th2.
Table 1 summarizes the demographical and clinical data of the participants
with SCI: age, sex, neurological level of injury (NLI) and ASIA impairment
score (AIS).

####  2.1.1 Inclusion and exclusion criteria

All participants had to be of age between 18 and 65 years. Given that the
experimental paradigm required a preserved arm function, all participants with
SCI had to have the injury at level C4 or lower. Participants with SCI were
excluded if they were artificially ventilated or had major spasms due to
possible interference with the EEG measurement. Control participants were
required to be able-bodied and with no history of neurological diseases.

###  2.2 Ethical approval and measurements

This study was approved by the local ethics committee of the Medical
University of Graz (ethical approval number 31-501 ex 18/19) and by the
Allgemeine Unfallversicherungsanstalt (AUVA) ethical committee. All
participants read and signed an informed consent form before the start of the
experiment and were paid for their participation. The EEG measurements of the
participants with SCI took place at AUVA Rehabilitation Clinic Tobelbad and
the EEG measurements of the control participants took place at Graz University
of Technology.

Participant | Age | Sex | NLI | AIS | Time since injury
---|---|---|---|---|---
P1 | 24 | F | C5 | B |  $>$ 10 years
P2 | 29 | M | C7 | C |  $>$ 10 years
P3 | 33 | M | Th2 | D |  $>$ 9 years
P4 | 36 | F | C7 | B |  $>$ 10 years
P5 | 37 | M | C4 | B |  $>$10 years
P6 | 39 | M | C6 | B |  $>$ 10 years
P7 | 48 | M | C4 | B | 6 - 12 months
P8 | 54 | M | C4 | B |  $>$ 1 year
Table 1: Summary of the demographical and clinical data of the participants
with SCI.

###  2.3 Hardware and electrodes’ layout

We recorded EEG data with a sampling rate of $500\text{\,}\mathrm{Hz}$ using
BrainAmp amplifiers and ActiCap caps (Brain Products, Munich, Germany) with 61
active electrodes positioned in a 10-10 layout, as detailed in Figure 1 of the
supplementary material. The ground electrode was placed on AFz and the
reference electrode was placed on the right mastoid. Additionally, we used 3
EOG electrodes that were placed above the nasion and below the outer canthi of
the eyes.

###  2.4 Experimental setup

Figure 1: Top: Experimental setup. Participants sat in front of a table,
attached which was a robotic arm. The participants controlled the robotic arm
during the trials using their hand. Bottom: The experimental protocol
displayed on the monitor. During the pre-trial period, the white square
represented the target of the coming trial. The small rectangle located
centrally, on the bottom part of the screen, represented the home position of
the participant’s hand. A trial started when the participants moved their hand
to its home position. The participants were instructed to move the robot to
the target square during the trials. After the trial (post-trial period), the
target changed colour, indicating whether or not it was reached, and the robot
automatically returned to its home position.

Similarly to the experimental setup described in [35], participants sat in
front of a table, on top of which was a wooden 4-sided box, with open sides
towards the participant and the tabletop, as depicted in Figure 1 (top). On
the ceiling of the box was a Leap Motion device (Leap Motion, San Francisco,
United States) that tracked the participants’ right or left hand, according to
their preferred hand. The position of the Leap Motion on the ceiling of the
box was adjusted to the handedness of each participant. The participants kept
their hand inside the wooden box. On the right side of the participants,
attached to the table, we placed a robotic arm (Jaco Assistive robotic arm -
Kinova Robotics, Bonn, Germany). Differently from [35], a screen monitor was
lying on the wooden box, centred in relation to the robotic arm. The monitor
was slightly inclined, with a 15-degree angle, to offer the participants a
better view of the screen. This change in relation to [35] was introduced to
minimize head and eye movements during the experiment.

###  2.5 Controlling the robotic arm

During the trials, participants could control the robotic arm on a horizontal
plane by moving their preferred hand on the tabletop. To reduce the range of
the participants’ movements, we considered the robot’s hand displacement to be
three times larger than the participants’ hand displacement.

Many participants with SCI had a very closed fist, due to hand spasticity
caused by their injury, and this impaired their hand’s recognition by the Leap
Motion. When this occurred, we inserted a small object in the participants’
hand in order to sustain the hand in a more open position and facilitate its
tracking.

###  2.6 Experiment overview

Before the experiment, we recorded one block in which the participant
performed eye movements [48, 49]. The experiment then consisted of 8 blocks of
30 trials each. $30\text{\,}\mathrm{\char 37\relax}$ of the trials of each
block were _error trials_ (9 trials). The remaining $70\text{\,}\mathrm{\char
37\relax}$ of the trials were _correct trials_ (21 trials). The sequence of
correct and error trials within each block was randomly generated using a
uniform distribution. We defined a maximum of 2 consecutive error trials in
each block and repeated the randomization procedure until the sequence of
trials satisfied this condition. Similarly, the trials of each block were
equally split between the right and the left targets. The sequence of targets
within each block were randomly assigned using a uniform distribution. We
defined a maximum of 3 consecutive trials with the same target in each block
and repeated the randomization procedure until the targets’ sequence satisfied
this requirement.

All the 8 blocks were online blocks: we used a _generic ErrP classifier_ in an
asynchronous manner to give participants real-time feedback of the _ErrP
detections_ during the experiment. For a matter of fluidity of the experiment,
we decided not to give participants feedback of false positive ErrP
detections, i.e., of the _ErrP detections_ that happened when no error had
occurred. This decision assured that all participants experienced the same
number of errors, which aimed to create a comparable expectation regarding the
occurrence of errors across participants. False ErrP detections can occur both
in correct and error trials and were considered when evaluating the
classifier. The details regarding the generic classifier are described in the
section _Generic ErrP classifier_.

###  2.7 Experimental protocol

During the pre-trial period, the monitor displayed two squares, on the top
part of the screen, each with a 14 cm side. As depicted in Figure 1, one of
the squares was filled in white and the other square had no fill. The filled
square represented the target of the coming trial. The centres of the squares
were 35 cm apart and their midpoint was located 30 cm in front of the home
position of the robot’s hand. On the bottom part of the screen was a
rectangle, representing the home position of the participant’s hand. The
position of the participant’s hand in relation to its home position was
depicted by a dot on the screen.

Participants could decide when to start a new trial and could rest for as long
as they needed in between trials. A trial started when the dot entered the
rectangle. This ensured that the participant’s hand was at a similar position
at the beginning of each trial. Participants were instructed to, when they
felt ready to start a new trial, position the dot representing their hand
below the home position’s rectangle, fixate their gaze on the target and
finally enter the rectangle from the bottom. This last step ensured a forward
movement of the robot. Participants were also asked to keep their gaze fixed
at the target during the entire trial in order to prevent eye movements.

The aim of each trial was to move the robot’s hand from its home position to
the target square. During the trials, only the two squares were displayed on
the screen: the white square representing the target and the square with no
fill. A trial ended when the robot’s hand was above the target or after 6
seconds (time out), in case the target has not been reached. After the end of
the trial (post-trial period), the target’s colour changed from white to
either green or red, for $1.2\text{\,}\mathrm{s}$, indicating whether or not
the target was reached, respectively. This feedback was always in line with
the robot’s behaviour. Then, the screen turned black, the robot automatically
returned to its home position and a new pre-trial period would start.

####  2.7.1 Error trials

In these trials, the paradigm triggered an error, during the movement of the
robot towards the target. The error consisted in interrupting the
participant’s control of the robot and adding a $5\text{\,}\mathrm{cm}$
upwards displacement to the robot’s hand. The participants perceived the error
by noticing the robot stopping and lifting and by realizing that the control
of the robot was lost. The errors occurred randomly, when the robot’s hand was
within 6 to 15 $\text{\,}\mathrm{cm}$ from its home position, in the forward
direction. For every error trial, this distance was drawn from a continuous
uniform distribution. In participants with SCI, the error onset occurred, on
average, $1.36\pm 0.14$ $\text{\,}\mathrm{s}$ after the start of the error
trial (mean $\pm$ std). In control participants, the error onset occurred, on
average, $1.30\pm 0.07$ $\text{\,}\mathrm{s}$ after the start of the error
trial (mean $\pm$ std).

We used the generic ErrP classifier in an asynchronous manner to give
participants feedback of the _ErrP detections_ occurring after the error
onset. Figure 2 illustrates all the possible interactions between the
participants and the robot during error trials, taking into account the
generic ErrP classifier feedback. If no ErrP was detected after the error
onset, the robot remained still for the rest of the trial. In this situation,
the total duration of the trial was 6 seconds and afterwards the target square
turned red. Differently, if an ErrP was detected by the classifier after the
error onset, the robot’s hand lowered $5\text{\,}\mathrm{cm}$ and the
participants regained its control. The downward movement informed the
participants of the ErrP detection and consequent regain of control. Since
participants instinctively stopped their hand movement when noticing the
error, they were instructed to reinitiate the movement and move the robot’s
hand to the selected target when regaining control of the robot. To
accommodate the extra movement, we added 6 seconds to the maximal trial
duration, once the first ErrP detection after the error onset occurred. If the
robot reached the target, after the error onset, the target square turned
green. Participants did not receive feedback of the false positive detections
occurring during the error trials, i.e., of the _ErrP detections_ occurring
before the error onset. Prior to the experiment, participants were informed
that errors would occur and were shown the characteristic robot movement
associated with error occurrence, i.e., the robot stopping and lifting.

####  2.7.2 Correct trials

In these trials, the paradigm did not trigger any error. Participants did not
receive feedback of the false positive ErrP detections occurring during the
correct trials. Figure 2 illustrates all the possible interactions between a
participant and the robot during correct trials. Correct trials lasted, on
average,  $2.11\;\pm$ $0.17\text{\,}\mathrm{s}$ for participants with SCI and
$2.05\;\pm$ $0.13\text{\,}\mathrm{s}$ for the control participants (mean $\pm$
std). All participants reached the target in more than
$99.4\text{\,}\mathrm{\char 37\relax}$ of the correct trials.

###  2.8 Data processing

Eye movements and blinks were removed online from the EEG data, using the
subspace subtraction algorithm [48, 49] and the eye movement data recorded
right before the start of the experiment. For the online detection of ErrPs
with the generic classifier, the EEG data were bandpass filtered between 1 and
10 Hz with a causal Butterworth filter of order 4. For the offline
electrophysiological analysis presented here, the EEG data were bandpass
filtered between 1 and 10 Hz with a noncausal Butterworth filter of order 4.

###  2.9 Defining events

In the error trials, we defined the error onset as the moment in which the
robot started its upwards displacement once the participant’s lost its
control. Prior to the experiment, we calculated the robot’s delay on 100
uncorrected errors, i.e., the time difference between the error marker and the
robot upwards displacement. This resulted in an average delay of $0.225\;\pm$
$0.005\text{\,}\mathrm{s}$ (mean $\pm$ std). Since the robot’s delay was
rather stable, we added the average delay to each recorded error marker in
order to obtain the error onset.

Correct trials had no clear onset. Therefore, to obtain comparable onsets in
correct and error trials for the electrophysiological analysis, we defined a
virtual onset for the correct trials at a time point in which errors could
occur in the error trials. For every participant, we defined the virtual onset
for his/her correct trials as the average time difference between the error
onsets and the start of the corresponding trials. For the participants with
SCI, the correct onset occurred, on average, $1.36\;\pm$
$0.14\text{\,}\mathrm{s}$ after the start of the correct trials (mean $\pm$
std). For the control participants, the correct onset occurred, on average,
$1.30\;\pm$ $0.07\text{\,}\mathrm{s}$ after the start of the correct trials
(mean $\pm$ std).

###  2.10 Generic ErrP classifier

We built a generic error-related potential classifier using the EEG data from
15 able-bodied participants of a previous study for ours [35]. None of these
previous participants took part in the experiment described here. The EEG data
from those participants were filtered between 1 and 10 Hz using a causal
Butterworth filter of order 4. Eye movements were removed from the data using
the subspace subtraction algorithm [48].

For each participant from [35], we used the 8 calibration runs of the dataset
and extracted an epoch with $450\text{\,}\mathrm{ms}$ from every trial. In the
error trials, the selected epoch started $300\text{\,}\mathrm{ms}$ after the
error onset. In the correct trials, the selected epoch started
$300\text{\,}\mathrm{ms}$ after the virtual onset. Hence, our initial features
were the amplitudes of the 61 EEG electrodes at all the time points of the
$450\text{\,}\mathrm{ms}$ of each epoch.

In order to remove outlier epochs, we first applied principal component
analysis (PCA) on the initial features and kept the PCA components that
explained $99\text{\,}\mathrm{\char 37\relax}$ of the data variability. Then,
we removed $1\text{\,}\mathrm{\char 37\relax}$ of the correct epochs and
$1\text{\,}\mathrm{\char 37\relax}$ of the error epochs as outliers. The
rejection criterion was based on a large Mahalanobis distance of the rejected
epochs within each class type (error or correct) in the PCA space. After this
step, 2475 correct epochs and 1059 error epochs were kept.

Finally, we repeated the PCA step on the initial feature space, after
discarding the outlier epochs, and kept as features the PCA components that
preserved $99\text{\,}\mathrm{\char 37\relax}$ of the data variability. This
step resulted in 412 PCA components. These components were then used as
features to train a shrinkage-LDA classifier with two classes: error and
correct [50]. The linear scores of the classifier were transformed into
probabilities using a softmax function. The PCA components preserved most of
the activity of the original space, as depicted in Figure 2 of the
supplementary material. Figure 3 of the supplementary material depicts the
classifier pattern, obtained by applying the discriminant feature analysis
(DFA) method to the training matrix with 3534 epochs and 412 features [51].
The generic classifier remained unchanged during the entire experiment. In
[45], we showed that the generic ErrP classifier offers a comparable
performance to a personalized ErrP classifier for the asynchronous detection
of ErrPs. Therefore, we chose not to retrain the classifier with the
participants’ own data.

###  2.11 ErrP detection

Similarly to the classifier developed in [35], the generic classifier
developed here was constructed to be used and evaluated in an asynchronous
manner. In the online experiment, the incoming EEG signals were analysed in
real-time by the ErrP classifier, which received as input an EEG window of
$450\text{\,}\mathrm{ms}$. Consecutive analyzed windows had a leap of
$18\text{\,}\mathrm{ms}$. The classifier’s evaluation of each window resulted
in the probability of the analysed window belonging to either class (correct
or error). Hence, the classifier produced a probability output every
$18\text{\,}\mathrm{ms}$, during the entire duration of each block. We defined
an _ErrP detection_ when two consecutive windows had a probability of
belonging to the error class above a certain threshold $\tau$. In [44], we
evaluated offline the asynchronous ErrP detection with the generic classifier
and tested the effect of varying the decision threshold. From [44] we
concluded that the combination of the generic ErrP classifier with a
personalized decision threshold leads to the achievement of better
performances. Hence, in this online experiment, we decided to apply this
strategy. The procedure to determine the personalized thresholds is described
in the section _Tailoring the decision threshold of the generic classifier to
each participant_.

Figure 2: Experimental protocol and metrics. Graphical representation of the
trial structure, of the interaction between the participant and the robot
during the trials, and of the trial-based metrics used for the evaluation of
the classifier. All the occurrences that are not labelled nor detailed,
inherit the corresponding description from the preceding node.

###  2.12 Metrics to evaluate the ErrP classifier

To evaluate the performance of the generic classifier, we considered the trial
structure of the experiment and the asynchronous nature of the decoding. The
proposed metrics assess a trial as successful or unsuccessful, based on the
asynchronous detection of ErrPs over the entire trial’s duration. This
strategy has been applied to the study of asynchronous detection of ErrPs and
other event-related potentials, in several other works [29, 52, 31, 30, 32,
53, 33, 54, 34, 44, 35, 45]. Figure 2 presents a graphical representation of
the metrics proposed here. Correct trials were labelled negative and error
trials were labelled positive.

####  2.12.1 True negative trials

We defined the true negative trials (TN trials) as the correct trials in which
no _ErrP detection_ occurred during the entire trial duration. For the
classifier’s evaluation we considered the true negative rate (TNR): the
fraction of correct trials that are TN trials, i.e., that have no _ErrP
detections_ 111 The metrics TNR and TPR used here address the asynchronous
detection of ErrPs in a trial-based scenario and are not directly comparable
with the TPR and TPR definitions commonly used in time-locked classification..

####  2.12.2 True positive trials

We defined the true positive trials (TP trials) as the error trials in which
no ErrP detection occurred before the error onset and at least one ErrP
detection occurred within the $1.5\text{\,}\mathrm{s}$ after the error onset.
For the classifier’s evaluation we considered the true positive rate (TPR):
the fraction of error trials that were TP trials. An additional metric, _ErrP
detection rate_ (EDR), considering only the ErrP detections within the 1.5 s
after the error onset, is defined in Figure 5 of the supplementary material,
where its relation with the TPR is detailed  ‡.

####  2.12.3 Chance level

To calculate the chance level for TNR and TPR we performed several
classifications with a classifier in which the training labels were randomly
permuted (500 times for the evaluation of the online detection with the
generic classifier and 50 times for the evaluation of the offline cross-
validation with a personalized classifier). Furthermore, we used permutation
based $p$-values to present the significance of the classification results
obtained with the generic ErrP classifier [55, 56].

###  2.13 False activation rate

The false activation rate (FAR) is the percentage of 1-second long intervals
that are contaminated with at least one false positive ErrP detection [57].
For this evaluation, we considered the entire duration of correct trials and
the period before the error onset in error trials. These periods were divided
into 1-second long intervals and these intervals were evaluated for the
presence of false positive ErrP detections.

###  2.14 Tailoring the decision threshold of the generic classifier to each
participant

In [44] we evaluated offline the asynchronous detection of ErrPs with a
generic classifier similar to the one described here. There, we observed that
the decision threshold $\tau$ that maximized the group performance was
$\tau=0.7$. Moreover, we also concluded that in order to optimize the
individual performance with the generic classifier, participants benefit from
the use of a personalized threshold. Therefore, in this experiment, we decided
to initiate the generic classifier with $\tau=0.7$ in the first block. This
enabled us to skip the offline calibration and allowed us to give participants
immediate feedback of their ErrP detections. Afterwards, we tailored $\tau$ to
each participant. After each of the first 3 blocks, we performed offline an
asynchronous classification with the generic ErrP classifier on all the
available data and tested thresholds between 0 and 1 in steps of 0.025. For
each of the 41 thresholds analysed, we calculated the corresponding TPR and
TNR. The TNR and TPR curves were further smoothed using a moving average with
7 samples. The smoothed curves were named smooth TPR and smooth TNR. For every
participant, we chose the threshold that maximized the product of the smooth
TPR and the smooth TNR. This was considered the threshold that maximized
performance and it was used in the next block. From block 4 onwards, the
generic ErrP classifier was combined with the threshold $\tau$ obtained after
the third block. The generic ErrP classifier was not retrained with the
participants’ data and only the decision threshold was updated based on the
data.

###  2.15 Evaluation of the generic ErrP classifier

We stopped tailoring $\tau$ to each participant after the third block because
we wanted to collect a substantial amount of data in unchanged conditions.
From blocks 4 to 8, all participants used the generic classifier with a fixed
but personalized threshold. Therefore, we only use the data from blocks 4 to 8
to evaluate the performance of the generic classifier, ensuring comparable
conditions across the participants.

###  2.16 Personalized ErrP classifier

In order to evaluate, offline, the performance of a personalized classifier,
we performed 10 times a 5-fold cross-validation in the entire dataset of each
participant, where a personalized classifier was tested in an asynchronous
manner in each fold. There, we also tested all thresholds from 0 to 1 in steps
of 0.025. For every participant, we obtained, in each fold, a TPR and a TNR
for every threshold tested. For every participant, we averaged the TPR and TNR
of the 50 iterations in the cross-validation, obtaining an average TPR and an
average TNR per participant. Finally, we selected the threshold that maximized
the product of the average TPR and the average TNR, for every participant. The
evaluation of the personalized classifier followed the metrics defined in the
section _Metrics to evaluate the ErrP classifier_.

##  3 Results

###  3.1 Neurophysiology

The electrophysiological results presented here comprise the entire recorded
dataset. Figure 3 shows the grand average correct and error signals at channel
FCz (green and red lines, respectively) for participants with SCI and control
participants. The green and red shaded areas depict the
$95\text{\,}\mathrm{\char 37\relax}$ confidence interval for the grand average
signals. The vertical line at $t=$ $0\text{\,}\mathrm{s}$ depicts the error
onset of the error trials and the virtual onset of the correct trials. For the
participants with SCI, the grand average error signal displays a negativity,
with peak amplitude of $-2.4\text{\,}\mathrm{\SIUnitSymbolMicro V}$ at time
$t=$ $0.154\text{\,}\mathrm{s}$ after the error onset, followed by a
positivity, with peak amplitude of $3.8\text{\,}\mathrm{\SIUnitSymbolMicro V}$
at time $t=$ $0.332\text{\,}\mathrm{s}$. For the control participants, the
grand average error signal displays a negativity, with peak amplitude of
$-5.5\text{\,}\mathrm{\SIUnitSymbolMicro V}$ at time $t=$
$0.176\text{\,}\mathrm{s}$ after the error onset, followed by a positivity,
with peak amplitude of $5.8\text{\,}\mathrm{\SIUnitSymbolMicro V}$ at time
$t=$ $0.334\text{\,}\mathrm{s}$. The grand average correct signal displays no
particular peaks, both in participants with SCI and control participants.
Figure 3 displays also the topographic plots of the grand average correct and
error signals at the time points of the peaks of the grand average error
signal.

Figure 3: Grand average correct and error signals at channel FCz (green and
red solid lines, respectively) for participants with SCI and control
participants. The shaded areas represent the $95\text{\,}\mathrm{\char
37\relax}$ confidence interval of the grand average curves. The vertical black
line at $t=$ $0\text{\,}\mathrm{s}$ represents the error onset of the error
trials and the virtual onset of the correct trials. The figure displays also
the topographic plots of the grand-average correct and error signals at the
time points of the peaks in the grand average error signal. Figure 4: Average
correct and error signals at channel FCz of every participant (green and red
lines, respectively). The shaded areas represent the $95\text{\,}\mathrm{\char
37\relax}$ confidence interval of the average signals. The black line at $t=$
$0\text{\,}\mathrm{s}$ represents the error onset of the error trials and the
virtual onset of the correct trials. The grey regions indicate the time points
in which correct and error signals were statistically different (Wilcoxon
ranksum tests, Bonferroni corrected, with $\alpha=0.01$).

As the morphology of the error signals was not homogeneous across
participants, we found it relevant to also present the electrophysiological
results of the individual participants. Figure 4 displays the average correct
and error signals at channel FCz (green and red lines, respectively) of every
participant. The green and red shaded areas depict the
$95\text{\,}\mathrm{\char 37\relax}$ confidence interval for the average
signals. The grey areas indicate the time points in which correct and error
signals were statistically different (Wilcoxon ranksum tests, Bonferroni
corrected, with $\alpha=0.01$). Figures 6 and 7 of the supplementary material
depict the topographic plots of the average correct and error signals of every
participant at different time points.

###  3.2 Adaptation of the classifier’s threshold in the first three
experimental blocks

Figure 5:  Optimization of the decision threshold used with the generic ErrP
classifier. Top: Evolution of the decision threshold: Initial threshold
($\tau=0.7$) and the calculated thresholds after each of the first 3 blocks,
for every participant. Bottom: TNR and TPR obtained offline, after the third
block (dashed green and red lines, respectively) and the corresponding smooth
curves (green and red solid lines). The blue line represents the threshold
that maximizes the product of the smooth curves, which is represented with a
black dotted line.

This experiment required no offline calibration and the participants received
feedback regarding their ErrP detections from its very beginning. This was
possible by combining the generic ErrP classifier with a generic decision
threshold ($\tau=0.7$) for the first experimental block. Still, we used the
first three experimental blocks to reach a fixed personalized decision
threshold. After each of the first three blocks, we updated the decision
threshold $\tau$ in order to maximize the participant’s performance. Hence,
participants used a generic classifier combined with a personalized decision
threshold from block 2 onwards. Figure 5 (top) depicts the initial threshold
($\tau=0.7$) and the calculated thresholds after each of the first 3 blocks,
for every participant. At the end of block 3, the average threshold was
$\tau=0.68$ for the participants with SCI and $\tau=0.59$ for the control
participants. Figure 5 (bottom) shows the TNR and TPR obtained offline after
block 3, for all the tested thresholds (green and red dashed lines,
respectively). It also shows the smooth TNR and smooth TPR, obtained with a
moving average (green and red solid lines). The black dotted line depicts the
product of these smooth curves and the blue vertical line indicates the
threshold that maximizes it. This is the decision threshold used for every
participant from blocks 4 to 8.

###  3.3 Evaluation of the online asynchronous classification with a generic
ErrP classifier

Figure 6: Evaluation of the generic ErrP classifier. Top: Classification
results in terms of true positive rate (TPR) and true negative rate (TNR). The
circles on the individual bars represent the chance-level of the corresponding
metrics. Bottom: Permutation based $p$-value regarding the significance of the
results obtained with the classifier.

To evaluate the asynchronous classification results obtained with the generic
ErrP classifier during the experiment, we only considered the data of the last
five blocks of the experiment, i.e., from blocks 4 to 8, since no parameters
were changed during these blocks.

Figure 6 (top) depicts the classification results obtained with the generic
classifier in terms of true positive rate (TPR) and true negative rate (TNR).
For participants with SCI, we obtained an average TPR of
$46.9\text{\,}\mathrm{\char 37\relax}$ and an average TNR of
$71.9\text{\,}\mathrm{\char 37\relax}$. For control participants, we obtained
an average TPR of $56.4\text{\,}\mathrm{\char 37\relax}$ and an average TNR of
$77.9\text{\,}\mathrm{\char 37\relax}$. The circles on the individual bars
represent the chance level of the corresponding metrics. The chance level
results for each participant were obtained by averaging the classification
results of 500 classifiers in which the training labels were randomly permuted
and by considering the final participant-specific threshold, as depicted in
Figure 4 of the supplementary material. Figure 6 (bottom) presents the
permutation based $p$-values regarding the significance of the classification
results [55, 56]. Figure 5 of the supplementary material depicts the
comparison between the metrics TPR and EDR.  Table 1 of the supplementary
material presents the false activation rate (FAR) in correct and error trials.

Figure 7 illustrates the online asynchronous detection of ErrPs and the
trials’ offline evaluation for participant C1. The dark grey areas represent
the trials and the white marks within them represent the ErrP detections. The
narrow rectangles colour code the trials’ offline evaluation. In these
rectangles, trials successfully classified (true positive trials and true
negative trials) are coded in white and trials with false positive ErrP
detections are coded in grey. The error trials with no ErrP detection are
coded in black.

Figure 7: Online detection of ErrPs and trials’ offline evaluation for
participant C1. Left: Error trials, aligned to the error onset (black vertical
line). Right: Correct trials, aligned to the start of the trial. The dark grey
areas represent the trials and the white marks within them represent the ErrP
detections. The narrow rectangles colour code the trials’ offline evaluation.
In these rectangles, trials successfully classified (true positive trials and
true negative trials) are coded in white and trials with false positive ErrP
detections are coded in light grey. Error trials with no ErrP detections are
coded in black.

###  3.4 Offline evaluation of the asynchronous ErrP classification with a
personalized classifier

To evaluate offline the asynchronous classification results with a
personalized classifier, we considered the 8 experimental blocks and performed
10 times a 5-fold cross-validation. As this evaluation was done offline, we
tested thresholds from 0 to 1 with a leap of 0.025 and the results obtained
are shown in function of the threshold $\tau$.

Figure 8 depicts the grand average TNR and TPR (green and red solid lines,
respectively) as well as the grand average chance level for TNR and TPR (green
and red dashed lines, respectively) in function of the threshold. The shaded
areas represent the $95\text{\,}\mathrm{\char 37\relax}$ confidence intervals
of the grand average curves. The chance level curves were obtained by
performing 10 times a 5-fold cross-validation with 50 classifiers in which the
labels of the training trials were randomly permuted.

Figure 8: Evaluation of the personalized ErrP classifier. Grand average TNR
and TPR (green and red solid lines, respectively) and grand average chance
level TNR and TPR (green and red dashed lines, respectively) in function of
the threshold. The shaded areas indicate the $95\text{\,}\mathrm{\char
37\relax}$ confidence interval for the grand average curves.

Figure 9 depicts, for every participant, the average TNR and TPR (green and
red solid lines, respectively) and the chance level TNR and TPR (green and red
dashed lines, respectively). The blue vertical line indicates the threshold
that maximizes the individual performance with the personalized ErrP
classifier.

Figure 9: Evaluation of the personalized ErrP classifier. Single subject
average TNR and TPR (green and red solid lines, respectively) and chance level
TNR and TPR in function of the threshold (green and red dashed lines,
respectively). The shaded areas indicate the $95\text{\,}\mathrm{\char
37\relax}$ confidence interval for the average curves. The blue vertical line
indicates the threshold that maximizes the individual performance. Figure 10:
Evaluation of the personalized ErrP classifier. Average TNR and TPR calculated
from the cross-validation procedure, with the optimal personalized threshold,
for every participant and their average. The small dots on each bar indicate
the chance level with the considered threshold, for every participant.

Figure 10 depicts the average TNR and TPR obtained in the cross-validation,
when using the optimal personalized decision threshold for every participant
(green and red bars, respectively). The small circles on the bars indicate the
chance level obtained for every participant with the considered threshold. For
participants with SCI, the grand average TNR was $77.9\text{\,}\mathrm{\char
37\relax}$ and the grand average TPR was $55.0\text{\,}\mathrm{\char
37\relax}$. For control participants, the grand average TNR was
$86.1\text{\,}\mathrm{\char 37\relax}$ and the grand average TPR was
$71.5\text{\,}\mathrm{\char 37\relax}$.

##  4 Discussion

In this work we investigated the transfer of a generic ErrP classifier from
able-bodied participants to participants with SCI. The classifier was
developed using the data from able-bodied participants from a previous
experiment of ours [35] and was tested asynchronously in a closed-loop online
experiment in which participants with SCI and able-bodied control participants
took part. Using the classifier asynchronously, the entire trials were
evaluated and not only a time-locked window. The online experiment required no
offline calibration period and the participants received feedback of the ErrP
detections immediately from the start of the experiment onwards. Additionally,
we also analysed the morphology of error-related potentials in participants
with SCI and in able-bodied control participants.

The grand average correct signal displayed, as expected, no particular
potential both in participants with SCI and in control participants. The
correct epochs correspond to the period in which the participants were
continuously controlling the robot and were not associated with any specific
event. The grand average error signal was associated with a fronto-central
activity both in participants with SCI and control participants. The peaks of
the grand average error signal were less pronounced in participants with SCI
than in control participants, as visible in Figure 3. This matches the results
described in [46]. Nevertheless, the electrophysiological patterns of
participants with SCI were rather heterogeneous and half of the participants
with SCI did not display the characteristic error-related activity
(participants P4, P5, P7 and P8). The remaining participants with SCI revealed
patterns comparable to control participants. Therefore, we believe that in our
study the decrease in peak amplitudes observed in the grand average error
signal of participants with SCI is not directly related with the injury, but
rather a consequence of heterogeneity of the signals in the population with
SCI. Several studies reported the effect of psychological factors, such as
depression and anxiety, on error-related potentials [58, 59]. The population
with SCI is particularly vulnerable to emotional disorders and higher levels
of distress [60, 61, 62]. Nevertheless, individual differences are large [62].
To consider a psychiatric evaluation and medication of the participants would
have been valuable for the current work and should be considered in future
studies involving a population with SCI [63]. Interestingly, the error signals
of control participants were also less homogeneous than in our previous
studies with a similar experimental protocol [34, 35]. Several studies showed
that ageing affects error processing and consequently the error-related
potentials, hence we hypothesise that the higher variability observed in the
signals of control participants of this study is related to the wider age
range of the participants in comparison with previous studies [64, 65, 66].

In order to interpret the classification results we focus on the TPR. This
metric considers an interval after the error onset and the period before the
error onset. Hence, it translates not only the classifier’s ability to decode
ErrPs after the occurrence of an error but also its ability to not detect
ErrPs when no error occurs. The TNR only captures the classifier’s ability to
not detect ErrPs when no error occurs. It is still a meaningful metric but the
TNR’s outcome can be artificially increased by the use of a high decision
threshold, as depicted in Figure 4 of the supplementary material. The
classification results of the generic classifier were, on average, lower in
participants with SCI than in control participants. Only half of the
participants with SCI obtained a TPR above chance level. These participants
were the ones that displayed clear error patterns. In the control
participants, seven out of eight participants obtained a TPR above chance
level. The remaining participant (participant C5) did not obtain a TPR above
chance level and did not display a very clear error signal. Summarizing, all
participants that displayed clear ErrP patterns in the electrophysiological
analysis obtained above chance level results with the generic classifier,
independently of the group (SCI or control). It would be rather interesting to
further investigate the factors that affect the error patterns, independently
of the spinal cord injury. These results support that using a generic ErrP
classifier is a valuable option to give immediate feedback to the
participants. Moreover, it indicates that ErrPs are transferable across
participants, and that the transfer can be applied to distinct populations,
such as participants with SCI.

With the generic classifier developed, participants received real-time
feedback of the ErrP detections from the beginning of the experiment. Still,
the first three blocks of the experiment were used to update the threshold
applied to the generic classifier. We made this choice because we had
previously shown that some participants strongly benefit from combining the
generic ErrP classifier with a personalized decision threshold [44]. For most
participants, the threshold was relatively stable after the first block. This
supports the use of a personalized threshold with the generic classifier, as
suggested in [44]. In a real-world online application, the occurrence of
errors can not be easily assessed, since it is determined by a subjective
perception of the BCI user. Such a constraint hinders an objective evaluation
of any ErrP classifier, unless the participants can use a motor-based strategy
to report the occurrence of errors. Still, our approach could be applied to a
real-world asynchronous situation in which the occurrence of errors is
unknown. Nevertheless, in order to establish a personalized decision
threshold, our approach would need, beforehand, a short online application in
which the occurrence of errors is known. Such application could be the
equivalent of one of our experimental blocks, which contained 9 errors and
lasted less than 5 minutes.

In our experiment, we only gave participants feedback of the ErrPs detected
after the error onset. This aimed to assure that participants experienced the
same number of errors and had comparable expectations regarding the occurrence
of errors. Providing participants with feedback of the false positive ErrP
detections would have brought our experiment closer to a real-world
application at the cost of putting participants in dissimilar circumstances,
given that false positive ErrP detections could affect their behaviour and the
generation of ErrPs. For instance, participants with many false positive ErrP
detections would certainly be affected by the feedback in a negative manner.
Either by losing engagement or by disregarding the feedback. Such participants
would no longer perceive the errors as meaningful and relevant and this could
alter their ErrPs.

When testing offline, the asynchronous classification with a personalized
classifier, two participants with SCI (participants P4 and P5) and one control
participant (participant C5) obtained chance level TPR results. This indicates
that the signals of these participants were not sufficiently different to
build a personalized classifier and these participants obtained chance level
results with both generic and personalized classifiers.

The classification results obtained with the personalized classifier are not
directly comparable with the results obtained with the generic classifier
because the classifiers were evaluated on different datasets. In a real-world
scenario, we could provide participants immediate feedback of their brain
signals using a generic classifier, while collecting data to train a
personalized classifier. Simultaneously, we could compare, at regular
intervals, the performance of the personalized and generic classifiers and
swap the generic classifier for a personalized classifier, once the latter
would grant significantly better performance.

##  5 Conclusion

Our work shows that a generic ErrP classifier can be used, asynchronously and
online, by participants with SCI and able-bodied participants. Moreover, the
generic ErrP classifier is transferable from an able-bodied population to a
population with SCI. The developed classifier required no previous calibration
with the participant and granted immediate feedback of the ErrP detections.
Therefore, our findings can help to widespread the incorporation of ErrPs in
BCIs for different types of users.

###  5.1 Acknowledgements

This work was supported by Horizon 2020 ERC Consolidator Grant 681231 ‘Feel
Your Reach’. The authors would like to thank Reinmar J. Kobler, Lea
Hehenberger and Joana Pereira for the fruitful discussions. The authors would
also like to thank the staff and patients of the AUVA-Rehabilitationsklinik
Tobelbad for the warm welcome.

##  6 References

## References

  * [1] C. Brunner, N. Birbaumer, B. Blankertz, C. Guger, A. Kübler, D. Mattia, J. del R Millán, F. Miralles, A. Nijholt, E. Opisso, N. Ramsey, P. Salomon, and G. R. Müller-Putz, “BNCI Horizon 2020: towards a roadmap for the BCI community,” Brain-Computer Interfaces, vol. 2, no. 1, pp. 1–10, 2015.
  * [2] J. d. R. Millán, R. Rupp, G. Müller-Putz, R. Murray-Smith, C. Giugliemma, M. Tangermann, C. Vidaurre, F. Cincotti, A. Kübler, R. Leeb, C. Neuper, K. Müller, and D. Mattia, “Combining brain–computer interfaces and assistive technologies: State-of-the-art and challenges,” Frontiers in Neuroscience, vol. 4, p. 161, 2010.
  * [3] J. R. Wolpaw, N. Birbaumer, D. J. McFarland, G. Pfurtscheller, and T. M. Vaughan, “Brain–computer interfaces for communication and control,” Clinical Neurophysiology, vol. 113, no. 6, pp. 767 – 791, 2002.
  * [4] G. Schalk, J. R. Wolpaw, D. J. McFarland, and G. Pfurtscheller, “EEG-based communication: presence of an error potential,” Clinical Neurophysiology, vol. 111, no. 12, pp. 2138 – 2144, 2000.
  * [5] P. W. Ferrez and J. d. R. Millán, “You are wrong!: Automatic detection of interaction errors from brain waves,” in Proceedings of the 19th International Joint Conference on Artificial Intelligence, IJCAI’05, pp. 1413–1418, 2005.
  * [6] N. Yeung, M. Botvinick, and J. D Cohen, “The neural basis of error detection: Conflict monitoring and the error-related negativity,” Psychological review, vol. 111, pp. 931–59, 11 2004.
  * [7] R. Chavarriaga, A. Sobolewski, and J. d. R. Millán, “Errare machinale est: the use of error-related potentials in brain-machine interfaces,” Frontiers in Neuroscience, vol. 8, p. 208, 2014.
  * [8] A. Llera, M. A. van Gerven, V. Gómez, O. Jensen, and H. J. Kappen, “On the use of interaction error potentials for adaptive brain computer interfaces,” Neural Networks, vol. 24, no. 10, pp. 1120–1127, 2011.
  * [9] A. Llera, V. Gómez, and H. J. Kappen, “Adaptive classification on brain-computer interfaces using reinforcement signals,” Neural Computation, vol. 24, no. 11, pp. 2900–2923, 2012.
  * [10] R. Yousefi, A. R. Sereshkeh, and T. Chau, “Online detection of error-related potentials in multi-class cognitive task-based bcis,” Brain-Computer Interfaces, vol. 6, no. 1-2, pp. 1–12, 2019.
  * [11] X. Artusi, I. K. Niazi, M.-F. Lucas, and D. Farina, “Accuracy of a BCI based on movement-related and error potentials,” in 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 3688–3691, Aug 2011.
  * [12] H. Zhang, R. Chavarriaga, Z. Khaliliardali, L. Gheorghe, I. Iturrate, and J. d. R. M. , “EEG-based decoding of error-related brain activity in a real-world driving task,” Journal of Neural Engineering, vol. 12, p. 066028, nov 2015.
  * [13] M. Spüler and C. Niethammer, “Error-related potentials during continuous feedback: using EEG to detect errors of different type and severity,” Frontiers in Human Neuroscience, vol. 9, p. 155, 2015.
  * [14] A. Kreilinger, H. Hiebel, and G. R. Müller-Putz, “Single versus multiple events error potential detection in a BCI-controlled car game with continuous and discrete feedback,” IEEE Transactions on Biomedical Engineering, vol. 63, no. 3, pp. 519–529, 2016.
  * [15] H. Si-Mohammed, C. Lopes-Dias, M. Duarte, F. Argelaguet, C. Jeunet, G. Casiez, G. R. Müller-Putz, A. Lécuyer, and R. Scherer, “Detecting system errors in virtual reality using eeg through error-related potentials,” in 2020 IEEE Conference on Virtual Reality and 3D User Interfaces (VR), pp. 653–661, 2020.
  * [16] L. C. Parra, C. D. Spence, A. D. Gerson, and P. Sajda, “Response error correction - a demonstration of improved human-machine performance using real-time EEG monitoring,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 11, pp. 173–177, June 2003.
  * [17] A. Kreilinger, C. Neuper, and G. R. Müller-Putz, “Error potential detection during continuous movement of an artificial arm controlled by brain–computer interface,” Medical & Biological Engineering & Computing, vol. 50, no. 3, pp. 223–230, 2012.
  * [18] T. Zander, L. Krol, N. Birbaumer, and K. Gramann, “Neuroadaptive technology enables implicit cursor control based on medial prefrontal cortex activity,” in Proceedings of the National Academy of Sciences of the United States of America, vol. 113(12), pp. 14898–14903, Dez 2016.
  * [19] A. F. Salazar-Gomez, J. DelPreto, S. Gil, F. H. Guenther, and D. Rus, “Correcting robot mistakes in real time using EEG signals,” 2017 IEEE International Conference on Robotics and Automation (ICRA), pp. 6570–6577, May 2017.
  * [20] S.-K. Kim, E. A. Kirchner, A. Stefes, and F. Kirchner, “Intrinsic interactive reinforcement learning - using error-related potentials for real world human-robot interaction,” Scientific Reports (Sci Rep), vol. 7: 17562, 12 2017.
  * [21] M. Mousavi, A. S. Koerner, Q. Zhang, E. Noh, and V. R. de Sa, “Improving motor imagery BCI with user response to feedback,” Brain-Computer Interfaces, vol. 4, no. 1-2, pp. 74–86, 2017.
  * [22] S. K. Ehrlich and G. Cheng, “Human-agent co-adaptation using error-related potentials,” Journal of Neural Engineering, vol. 15, p. 066014, sep 2018\.
  * [23] M. Spüler, M. Bensch, S. Kleih, W. Rosenstiel, M. Bogdan, and A. Kübler, “Online use of error-related potentials in healthy users and people with severe motor impairment increases performance of a P300-BCI,” Clinical Neurophysiology, vol. 123, no. 7, pp. 1328 – 1337, 2012.
  * [24] R. J. Kobler, A. I. Sburlea, and G. R. Müller-Putz, “Tuning characteristics of low-frequency eeg to positions and velocities in visuomotor and oculomotor tracking tasks,” Scientific Reports, vol. 8, no. 1, p. 17713, 2018.
  * [25] P. Ofner and G. R. Müller-Putz, “Decoding of velocities and positions of 3D arm movement from EEG,” in 2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 6406–6409, Aug 2012\.
  * [26] Y. Inoue, H. Mao, S. B. Suway, J. Orellana, and A. B. Schwartz, “Decoding arm speed during reaching,” Nature Communications, vol. 9, no. 1, p. 5243, 2018\.
  * [27] A. A. Nooh, J. Yunus, and S. M. Daud, “A review of asynchronous electroencephalogram-based brain computer interface systems,” in International Conference on Biomedical Engineering and Technology, 2011.
  * [28] B. Z. Allison, C. Brunner, C. Altstätter, I. C. Wagner, S. Grissmann, and C. Neuper, “A hybrid ERD/SSVEP BCI for continuous simultaneous two dimensional cursor control,” Journal of Neuroscience Methods, vol. 209, no. 2, pp. 299 – 307, 2012.
  * [29] J. Omedes, I. Iturrate, and L. Montesano, “Detection of event-less error related potentials,” in Proceedings of IROS 2013 Workshop on Neuroscience and Robotics ”Towards a robot-enabled, neuroscience-guided healthy society”, 2013.
  * [30] J. Omedes, I. Iturrate, and L. Montesano, “Asynchronous detection of error potentials,” in Proceedings of the 6th Brain-Computer Interface Conference 2014, 2014.
  * [31] J. Omedes, I. Iturrate, J. Minguez, and L. Montesano, “Analysis and asynchronous detection of gradually unfolding errors during monitoring tasks,” Journal of Neural Engineering, vol. 12, no. 5, p. 056001, 2015\.
  * [32] J. Omedes, I. Iturrate, R. Chavarriaga, and L. Montesano, “Asynchronous decoding of error potentials during the monitoring of a reaching task,” in 2015 IEEE International Conference on Systems, Man, and Cybernetics (SMC2015), pp. 3116–3121, Oct 2015.
  * [33] C. Lopes-Dias, A. I. Sburlea, and G. R. Müller-Putz, “Error-related potentials with masked and unmasked onset during continuous control and feedback,” in 7th Graz Brain-Computer Interface Conference 2017, pp. 320–332, September 2017.
  * [34] C. Lopes-Dias, A. I. Sburlea, and G. R. Müller-Putz, “Masked and unmasked error-related potentials during continuous control and feedback,” Journal of Neural Engineering, vol. 15, p. 036031, apr 2018.
  * [35] C. Lopes-Dias, A. I. Sburlea, and G. R. Müller-Putz, “Online asynchronous decoding of error-related potentials during the continuous control of a robot,” Scientific Reports, vol. 9, no. 1, p. 17596, 2019.
  * [36] I. Iturrate, L. Montesano, R. Chavarriaga, J. del R. Millán, and J. Minguez, “Minimizing calibration time using inter-subject information of single-trial recognition of error potentials in brain-computer interfaces,” in 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 6369–6372, 2011.
  * [37] I. Iturrate, R. Chavarriaga, L. Montesano, J. Minguez, and J. d. R. Millán, “Latency correction of error potentials between different experiments reduces calibration time for single-trial classification,” in 2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pp. 3288–3291, Aug 2012.
  * [38] I. Iturrate, R. Chavarriaga, L. Montesano, J. Minguez, and J. Millán, “Latency correction of event-related potentials between different experimental protocols,” Journal of Neural Engineering, vol. 11, p. 036005, apr 2014.
  * [39] S. K. Kim and E. A. Kirchner, “Handling few training data: Classifier transfer between different types of error-related potentials,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 24, pp. 320–332, March 2016.
  * [40] S. K. Kim and E. A. Kirchner, “Classifier transferability in the detection of error related potentials from observation to interaction,” in 2013 IEEE International Conference on Systems, Man, and Cybernetics, pp. 3360–3365, Oct 2013.
  * [41] S. Bhattacharyya, A. Konar, D. N. Tibarewala, and M. Hayashibe, “A generic transferable EEG decoder for online detection of error potential in target selection,” Frontiers in Neuroscience, vol. 11, p. 226, 2017.
  * [42] S. K. Ehrlich and G. Cheng, “A feasibility study for validating robot actions using EEG-based error-related potentials,” International Journal of Social Robotics, Nov 2018.
  * [43] F. M. Schönleitner, L. Otter, S. K. Ehrlich, and G. Cheng, “Calibration-free error-related potential decoding with adaptive subject-independent models: A comparative study,” IEEE Transactions on Medical Robotics and Bionics, vol. 2, no. 3, pp. 399–409, 2020.
  * [44] C. Lopes-Dias, A. I. Sburlea, and G. R. Müller-Putz, “Asynchronous detection of error-related potentials using a generic classifier,” in 8th International Brain Computer Interface Conference 2019, Sep 2019.
  * [45] C. Lopes-Dias, A. I. Sburlea, and G. R. Müller-Putz, “Generic error-related potential classifier offers a comparable performance to a personalized classifier,” in EMBC, 2020.
  * [46] P. Keyl, M. Schneiders, C. Schuld, S. Franz, M. Hommelsen, N. Weidner, and R. Rupp, “Differences in characteristics of error-related potentials between individuals with spinal cord injury and age- and sex-matched able-bodied controls,” Frontiers in Neurology, vol. 9, p. 1192, 2019.
  * [47] A. Kumar, Q. Fang, J. Fu, E. Pirogova, and X. Gu, “Error-related neural responses recorded by electroencephalography during post-stroke rehabilitation movements,” Frontiers in Neurorobotics, vol. 13, p. 107, 2019.
  * [48] R. J. Kobler, A. I. Sburlea, and G. R. Müller-Putz, “A comparison of ocular artifact removal methods for block design based electroencephalography experiments,” in Proceedings of the 7th Graz Brain-Computer Interface Conference 2017, Sep 2017.
  * [49] R. J. Kobler, “Corneo-retinal-dipole and eyelid-related eye artifacts can be corrected offline and online in electroencephalographic and magneticencephalographic signals,” NeuroImage, 2020.
  * [50] B. Blankertz, S. Lemm, M. Treder, S. Haufe, and K.-R. Müller, “Single-trial analysis and classification of ERP components — a tutorial,” NeuroImage, vol. 56, no. 2, pp. 814 – 825, 2011.  Multivariate Decoding and Brain Reading.
  * [51] G. A. Giraldi, P. S. Rodrigues, E. C. Kitani, J. R. Sato, and C. E. Thomaz, “Statistical learning approaches for discriminant features selection,” Journal of the Brazilian Computer Society, vol. 14, pp. 7–22, Jun 2008.
  * [52] E. López-Larraz, L. Montesano, Á. Gil-Agudo, and J. Minguez, “Continuous decoding of movement intention of upper limb self-initiated analytic movements from pre-movement eeg correlates,” Journal of NeuroEngineering and Rehabilitation, vol. 11, p. 153, Nov 2014.
  * [53] A. I. Sburlea, L. Montesano, and J. Minguez, “Continuous detection of the self-initiated walking pre-movement state from EEG correlates without session-to-session recalibration,” Journal of Neural Engineering, vol. 12, p. 036007, apr 2015.
  * [54] J. Pereira, A. I. Sburlea, and G. R. Müller-Putz, “Eeg patterns of self-paced movement imaginations towards externally-cued and internally-selected targets,” Scientific Reports, vol. 8, p. 13394, Sep 2018.
  * [55] P. Good and P. Good, Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses.  Springer series in statistics, Springer, 2000.
  * [56] M. Ojala and G. C. Garriga, “Permutation tests for studying classifier performance,” J. Mach. Learn. Res., vol. 11, p. 1833–1863, Aug. 2010\.
  * [57] M. Fatourechi, R. Ward, and G. Birch, “Evaluating the performance of a self-paced bci with a new movement and using a more engaging environment,” Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Conference, vol. 2008, pp. 650–3, 02 2008.
  * [58] M. Ruchsow, B. Herrnberger, P. Beschoner, G. Grön, M. Spitzer, and M. Kiefer, “Error processing in major depressive disorder: Evidence from event-related potentials,” Journal of Psychiatric Research, vol. 40, no. 1, pp. 37 – 46, 2006.
  * [59] D. M. Olvet, D. N. Klein, and G. Hajcak, “Depression symptom severity and error-related brain activity,” Psychiatry Research, vol. 179, no. 1, pp. 30 – 37, 2010.
  * [60] Z. Khazaeipour, S.-M. Taheri-Otaghsara, and M. Naghdi, “Depression following spinal cord injury: Its relationship to demographic and socioeconomic indicators,” Topics in spinal cord injury rehabilitation, vol. 21, no. 2, pp. 149–155, 2015.  26364284[pmid].
  * [61] C. Migliorini, B. Tonge, and G. Taleporos, “Spinal cord injury and mental health,” Australian & New Zealand Journal of Psychiatry, vol. 42, no. 4, pp. 309–314, 2008.  PMID: 18330773.
  * [62] M. W. M. Post and C. M. C. van Leeuwen, “Psychosocial issues in spinal cord injury: a review,” Spinal Cord, vol. 50, no. 5, pp. 382–389, 2012.
  * [63] C. M. C. van Leeuwen, L. H. V. van der Woude, and M. W. M. Post, “Validity of the mental health subscale of the sf-36 in persons with spinal cord injury,” Spinal Cord, vol. 50, no. 9, pp. 707–710, 2012.
  * [64] M. Falkenstein, J. Hoormann, and J. Hohnsbein, “Changes of error-related ERPs with age,” Experimental Brain Research, vol. 138, no. 2, pp. 258–262, 2001.
  * [65] S. Hoffmann and M. Falkenstein, “Aging and error processing: age related increase in the variability of the error-negativity is not accompanied by increase in response variability,” PloS one, vol. 6, pp. e17482–e17482, Feb 2011.
  * [66] S. Nieuwenhuis, K. R. Ridderinkhof, D. Talsma, M. G. H. Coles, C. B. Holroyd, A. Kok, and M. W. van der Molen, “A computational account of altered error processing in older age: Dopamine and the error-related negativity,” Cognitive, Affective, & Behavioral Neuroscience, vol. 2, no. 1, pp. 19–36, 2002\.


[a]A. Sinha

# Search for enhanced TeV gamma-ray emission from Giant Molecular Clouds using
H.E.S.S.

V. Baghmanyan  G. Peron  Y. Gallant  S. Casanova  M. Holler  A. Mitchell

###### Abstract

Cosmic Ray (CR) interactions with the dense gas inside Giant Molecular Clouds
(GMCs) produce neutral pions, which in turn decay into gamma rays. Thus, the
gamma ray emission from GMCs is a direct tracer of the cosmic ray density and
the matter density inside the clouds. Detection of enhanced TeV emission from
GMCs, i.e., an emission significantly larger than what is expected from the
average Galactic cosmic rays illuminating the cloud, can imply a variation in
the local cosmic ray density, due to, for example, the presence of a recent
accelerator in proximity to the cloud.

Such gamma-ray observations can be crucial in probing the cosmic ray
distribution across our Galaxy, but are complicated to perform with present
generation Imaging Atmospheric Cherenkov Telescopes (IACTs). These studies
require differentiating between the strong cosmic-ray induced background, the
large scale diffuse emission, and the emission from the clouds, which is
difficult to the small field of view of present generation IACTs.

In this contribution, we use H.E.S.S. data collected over 16 years to search
for TeV emission from GMCs in the inner molecular galacto-centric ring of our
Galaxy. We implement a 3D FoV likelihood technique, and simultaneously model
the hadronic background, the galactic diffuse emission and the emission
expected from known VHE sources to probe for excess TeV gamma ray emission
from GMCs.

##  1 Introduction

The paradigm of origin and propagation of cosmic rays (CRs) is based on direct
observations [1] of CRs in the vicinity of the solar system [2]. However, in
specific regions of the Galaxy, the level of CRs differs from the spectrum
measured in the vicinity of Earth. Fermi-LAT [3, 4] observations show an up to
$\sim$4 times higher density in the region around 4 kpc from the Galactic
center, where the spectrum is also slightly harder. This challenges our
current understanding on the propagation of CRs in the Galaxy. While a radial
dependence of the propagation can give rise to such effects [5], these effects
can also be generated by the higher density of accelerators in the inner
galactic regions [6].

Giant Molecular Clouds (GMCs) provide us with unique conditions for testing
the density of CR in the Galactic disc. CRs interaction with the dense gas
produces neutral pions, which in turn decay into gamma rays. Probing GMCs at
different spatial points across our galaxy makes it possible to trace the
distribution of galactic CRs point by point [7]. It is interesting to
understand if the enhancement and hardening seen at different galacto-centric
rings are localized in regions coincident with clouds or if they are
characteristic of a larger region; it is only by directly probing GMCs that
such information can be obtained. Such studies have been performed extensively
with Fermi-LAT [6, 4, 8] and enhanced cosmic ray density have been found near
many clouds, especially in the inner galacto-centric rings.

Probing GMCs at Very High Energies (VHE) will be crucial to understand if the
hardening seen at GeV energies continues in the TeV regime as well,
constraining whether it is a local or global behaviour and if it is related to
accelerators or to propagation effects.

Studies of TeV emission from GMCs are, however, relatively rare. Enhanced
emission has been reported from the galactic ridge by the H.E.S.S.
collaboration [9] spatially coincident with a complex of GMCs, while studies
of high latitude clouds by HAWC have yielded upper limits consistent with
expectations of a homogeneous “sea”of CRs [10]. Detection of such clouds with
Imaging Atmospehric Cherenkov Telescopes (IACT) is challenging because it
requires making a separation between the large residual hadronic background,
the large scale diffuse emission, and the emission from the cloud, which is
complicated due to the limited field of view of IACTs.

In this work, we use a full 3D likelihood method to separate the diffuse
emission from the hadronic background, showing that this allows us to directly
probe GMCs at TeV energies with a method similar to Fermi-LAT. As a proof of
concept, we use this to investigate the region of cloud 877 from Rice et al.’s
catalog [11], located at 5.5 kpc from the Galactic Center. This cloud has been
detected by Fermi-LAT to have a higher and harder spectrum than locally
expected [6]. We show that this deviation continues at TeV energies, and thus,
report on the detection of significant emission from the direction of this
cloud at VHE.

##  2 H.E.S.S observations

The High Energy Stereoscopic System (H.E.S.S) is an array of five imaging
atmospheric Cherenkov telescopes (IACTs) located in the Khomas Highland of
Namibia, 1800 m above sea level. Results presented here contain data collected
between 2004 - 2019, using the four 12m diameter telescopes, CT1-4, which have
a field of view of $5\degree$. Observations used here have been carried out
either as dedicated pointings on multiple individual objects or during the
Galactic plane scan for the HGPS campaign. Calibration, reconstruction and
$\gamma$-hadron separation were performed using the ImPACT chain [12] within
the H.E.S.S. analysis package, and the counts and instrument response
functions (IRFs) exported to a DL3 format as specified in the Gamma Astro Data
Format111https://gamma-astro-data-formats.readthedocs.io/en/v0.2/

##  3 3D Field of View likelihood technique

A crucial element in any IACT analysis is the rejection of the residual
hadronic background from $\gamma$-like cosmic ray induced showers. Standard
background estimation techniques, e.g. the ring background or the reflected-
region [13] rely on a measurement of the background in supposedly source-free
regions within the observed field of view. However, such approaches are poorly
suited for our present case. They exclude faint emission at the level of the
diffuse gas, and therefore are not suitable to trace molecular clouds, for
which we expect a factor $\lesssim$ 10% difference from the diffuse emission
at TeV energies. A distinction between the of the hadronic component and the
true diffuse emission is essential.

A three dimensional (3D) likelihood analysis provides us with a much more
sensitive technique. While such techniques have been routinely used for high
energy $\gamma$-ray data processing, its implementation in the field of IACT
analysis is relatively recent [14, 15]. A spectro-morphological template model
is constructed for the cosmic ray-induced background from archival runs with
mostly empty fields of view. These runs are grouped according to zenith angle
(and other observation conditions like optical efficiency, and an acceptance
model is derived for each bin. Then, given an observation corresponding to
certain parameters, a model background map is created which can then be
directly fit as a separate component to the observed data. For details of
background modelling, see [15].

The observed data is then described by a combination of many 3D (energy and 2
spatial dimensions) models, one for each expected emission component plus the
constructed background model. The models are fitted to the data via a
likelihood formalism and the significance of specific components determined by
means of likelihood ratio tests.

Such a likelihood analysis technique is implemented within the CTA Science
tools, gammapy [16] . In this contribution, we use
gammapy0.18222https://docs.gammapy.org/ [17] to simultaneously model the
residual hadronic background, the large scale diffuse emission and the
emission from GMCs, and show that, for the first time, we achieve the
detection of emission from the direction of a molecular cloud at TeV energies
illuminated by background CRs.

###  3.1 Analysis set-up

The DL3 data are reduced following the standard gammapy analysis procedure. An
energy range of $0.5-20$ TeV is chosen for the analysis, with events falling
beyond an offset of 2.0°from the camera centre, or with energies below the
peak of the background spectrum ignored.

The background models, as constructed in the previous section, suffer from
systematic fluctuations, and before stacking multiple runs, corrections to the
model are required for each individual observation by fitting it outside
exclusion region. For analysis of molecular clouds at low Galactic latitudes,
it is crucial to make an optimal choice for the background exclusion region -
exclude the maximal possible gamma-ray emission, while retaining enough off
counts to do a do run by run background fitting. Hence, we choose an exclusion
mask based on the Planck dust maps [18] which trace the interstellar gas. We
excluded the dense gas, i.e., all pixels with a value of density above
$2\times 10^{22}$ cm-2. Also, all known TeV sources are masked (Figure 1). A
spectral correction to the background, (normalisation and tilt), is then fit
for each run, and final stacked counts, background and IRFs constructed.

Figure 1: Exclusion mask used for the field of view normalisation: (left) A
mask on the known HGPS sources, (middle) Mask constructed from dust column
density for the diffuse emission, (right) The full mask, a product of the
previous two.

We use a technique similar to what is followed in [4, 8, 6, 10] to construct
spatial templates for the pion decay emission. We use the dust 353 GHz opacity
map which is a tracer of both molecular and atomic gas, and realize a
rectangular cutout (Figure 2) around the ellipse that represents the cloud
extension, with $\sigma_{minor}=0.22$ °, $\sigma_{major}=0.48$ °, and
position_angle = 166 °as given in the clouds catalog [11]. To probe the
emission from the cloud, we use the cut-out on the cloud as computed above.
The extracted spectrum refers to the dust template region which includes the
cloud. Even though the extraction region formally includes the entire gas on
the line of sight, in the case of cloud 877, the dominant contribution arises
from a region in the velocity range coincident with the cloud which accounts
for $\sim$60% of the gas column. As demonstrated in [19], observations of the
entire column are a viable alternative in these cases, and can be successfully
compared to the theoretical prediction that we have for the flux of galactic
cosmic rays.

The rest of the dust map, with the cloud cut out, is used as a template for
the large scale diffuse emission. Not modelling this correctly can lead to an
over-prediction of the hadronic background, and in turn, an under-prediction
of the flux from the direction of the cloud. Using a cut-out on the cloud
allows us to analyze the diffuse emission as separate component and to extract
the spectrum from there. The use of dust allows us to reduce the uncertainties
related to the $X_{CO}$ conversion factor, the HI spin temperature, and the
untraced (a.k.a. dark) gas. However, the downside is that a kinematic
separation of the cloud is not possible, and the entire gas column in the
given direction is taken as the cloud template. The maps are normalised such
that they integrate to unity on the target geometry. The spectral model is
assumed to be a simple power law,
$F_{\gamma}(E)=K(\frac{E}{E_{0}})^{-\alpha}$. We keep the known sources masked
during the likelihood fit.

Figure 2: The templates used for modelling the diffuse emission and the cloud
in the left and right panel, respectively. The elliptical size as quoted in
[11] is plotted in red.

##  4 Fermi-LAT analysis

We performed an analysis of the region of the cloud 877 similar to the one in
[6], using the same spatial template applied in the H.E.S.S. analysis, namely
a rectangular cutout based on dust templates. With respect to the published
results [6], based on 9 years of observations, we used here a collection of
$\sim$ 12 years of data accumulated from August 4th 2008 to January 8th 2020.
The standard quality cuts have been applied (DATA_QUAL==1 && LAT_CONFIG==1 and
$z_{max}=90^{\circ}$). We used SOURCE Pass8 data converted both at the front
and at the back of the detector. Our model included all the sources of the
4FGL catalog. As done in [6], we proceeded by fitting the diffuse sources
first, with the parameters of the point sources fixed to the cataloged value.
After a first fit, we proceeded with the optimization of the sources in the
catalog. Then we performed a likelihood fit, freeing all the parameters of the
sources within $3\degree$ from the center and leaving the normalization of the
brightest sources (TS>1000) free. Finally, we investigated the residuals and
included in the model any residual spots with TS >25\. To investigate the
influence of these new spots, we calculated the difference between the SED
obtained for the cloud before and after adding the new sources. The difference
was found to be lower than 2% in each energy bin.

##  5 Results

Significant emission ($\Delta(TS)=101$) is detected from the direction of
cloud 877. In Figure 3, we show the significance map obtained before and after
modelling the cloud, overlaid with $5\sigma$ significance contour lines and
dust column densities. The H.E.S.S. significance contours well coincide with
the gas contours, suggesting that the emission is truly correlated with the
cloud.

Figure 3: Top: Significance map of the cloud on a full RoI (left) and a cutout
on the cloud region (right) with 4-sigma significance contour in green and
9e22 cm$ˆ{-2}$ dust column density contour in blue); Bottom: Significance maps
after modelling the emission from the cloud

The fitted spectrum ($\alpha=2.73\pm 0.16$) is well consistent with the
extrapolation of the Fermi-LAT spectrum ($\alpha=2.45\pm 0.01$). In Figure 4,
we show the results of a joint fit ($\alpha=2.48\pm 0.01$) between the
H.E.S.S. and Fermi-LAT spectral points. With the opportune scaling, the
expected flux of a cloud can be factorized in the emissivity $\phi(E)$ per H
atom, and a factor $A\equiv\frac{M_{5}}{d_{kpc}^{2}}$, where $M_{5}$=M/105 M⊙
and $d_{kpc}$ is the distance of the cloud in kpc, which accounts for the
column density in the cloud. We followed [6] to compute the $\gamma$-ray
emissivity, expected from a cloud illuminated by the local CR spectrum. We
used a nuclear enhancement factor of 1.8 and a fit of the proton spectrum
which interpolates the data of AMS-02 (up to $\sim$ 1 TeV) and DAMPE (up to
$\sim$ 100 TeV) [1]. For cloud 877, $M_{5}=20$ and $d_{kpc}=3.4$, implying
$A=1.8$. Here the mass has been calculated for the cutout region used for the
analysis, approximating the entire column to be located at the cloud distance,
which is taken from [11].

Interestingly, the TeV observation confirms the enhancement and hardening with
respect to the local flux already reported for this region at GeV energies
with Fermi-LAT. The observed GeV-TeV emission is enhanced with respect to the
local emissivity as derived from AMS-DAMPE spectral points, and significantly
harder. With a joint fit of the Fermi-LAT and H.E.S.S. flux points using a
pion decay model [20] as implemented in the python package naima [21], we
estimate the differential cosmic-ray proton number density inside the cloud to
be $\sim 1.5\times 10^{-17}$/GeV/cm3 (at 1 TeV), with an index of $2.58\pm
0.01$. This corresponds to an energy density of cosmic-ray protons with
energies in the range 3–30 TeV (roughly corresponding to the gamma-ray energy
range of the H.E.S.S. flux points) of $1.0\times 10^{-2}$ eV/cm3, which is
$\sim 5-6$ times the local density.

Figure 4: Joint fit between Fermi-LAT and H.E.S.S. spectral points using Left:
An analytic power law spectral model. The gamma-ray spectrum expected from the
local CR measurements is plotted in violet, the shaded region representing
uncertainties in the cloud column density. Right A pion decay model as
implemented in Naima.

##  6 Conclusions

We show that a full 3D likelihood analysis setup provides us with a powerful
method of directly studying emission from GMCs. This technique is validated
applied on Cloud 877 which is known to show excess emission in gamma rays, and
we report on the first detection excess emission coincident with a passive GMC
at VHE energies. This opens up exciting possibilities to probe the cosmic ray
distribution across our Galaxy, and a detailed study with more clouds is in
progress.

## Acknowledgements

The support of the Namibian authorities and of the University of Namibia in
facilitating the construction and operation of H.E.S.S. is gratefully
acknowledged, as is the support by the German Ministry for Education and
Research (BMBF), the Max Planck Society, the German Research Foundation (DFG),
the Helmholtz Association, the Alexander von Humboldt Foundation, the French
Ministry of Higher Education, Research and Innovation, the Centre National de
la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à
l’énergie atomique et aux énergies alternatives (CEA), the U.K. Science and
Technology Facilities Council (STFC), the Knut and Alice Wallenberg
Foundation, the National Science Centre, Poland grant no. 2016/22/M/ST9/00382,
the South African Department of Science and Technology and National Research
Foundation, the University of Namibia, the National Commission on Research,
Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of
Education, Science and Research and the Austrian Science Fund (FWF), the
Australian Research Council (ARC), the Japan Society for the Promotion of
Science and by the University of Amsterdam.

We appreciate the excellent work of the technical support staff in Berlin,
Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tûbingen and in Namibia in the
construction and operation of the equipment. This work benefitted from
services provided by the H.E.S.S. Virtual Organisation, supported by the
national resource providers of the EGI Federation.

## References

  * [1] P. Lipari and S. Vernetto, _The shape of the cosmic ray proton spectrum_ , _Astroparticle Physics_ 120 (2020) 102441.
  * [2] S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando et al., _The origin of Galactic cosmic rays: Challenges to the standard paradigm_ , _International Journal of Modern Physics D_ 28 (2019) 1930022 [1903.11584].
  * [3] F. Acero, M. Ackermann, M. Ajello, A. Albert, L. Baldini, J. Ballet et al., _Development of the Model of Galactic Interstellar Emission for Standard Point-source Analysis of Fermi Large Area Telescope Data_ , _The Astrophysical Journal Supplement Series_ 223 (2016) 26 [1602.07246].
  * [4] R.-z. Yang, E. de Oña Wilhelmi and F. Aharonian, _Probing cosmic rays in nearby giant molecular clouds with the fermi large area telescope_ , _Astronomy & Astrophysics_ 566 (2014) A142.
  * [5] M. Pothast, D. Gaggero, E. Storm and C. Weniger, _On the progressive hardening of the cosmic-ray proton spectrum in the inner galaxy_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 045.
  * [6] F. Aharonian, G. Peron, R. Yang, S. Casanova and R. Zanin, _Probing the sea of galactic cosmic rays with fermi-lat_ , _Phys. Rev. D_ 101 (2020) 083018.
  * [7] S. Casanova, F.A. Aharonian, Y. Fukui, S. Gabici, D.I. Jones, A. Kawamura et al., _Molecular Clouds as Cosmic-Ray Barometers_ , _pasj_ 62 (2010) 769 [0904.2887].
  * [8] V. Baghmanyan, G. Peron, S. Casanova, F. Aharonian and R. Zanin, _Evidence of Cosmic-Ray Excess from Local Giant Molecular Clouds_ , _The Astrophysical Journal Letters_ 901 (2020) L4 [2009.08893].
  * [9] F. Aharonian, A.G. Akhperjanian, A.R. Bazer-Bachi, M. Beilicke, W. Benbow, D. Berge et al., _Discovery of very-high-energy $\gamma$-rays from the Galactic Centre ridge_, _Nature_ 439 (2006) 695 [astro-ph/0603021].
  * [10] A. Albert, R. Alfaro, C. Alvarez, J. Camacho, J. Arteaga-Velázquez, K. Arunbabu et al., _Probing the sea of cosmic rays by measuring gamma-ray emission from passive giant molecular clouds with hawc_ , _arXiv preprint arXiv:2101.08748_ (2021) .
  * [11] T.S. Rice, A.A. Goodman, E.A. Bergin, C. Beaumont and T.M. Dame, _A Uniform Catalog of Molecular Clouds in the Milky Way_ , _The Astrophysical Journal_ 822 (2016) 52 [1602.02791].
  * [12] R.D. Parsons and J.A. Hinton, _A Monte Carlo Template based analysis for Air-Cherenkov Arrays_ , _Astropart. Phys._ 56 (2014) 26 [1403.2993].
  * [13] D. Berge, S. Funk and J. Hinton, _Background modelling in very-high-energy $\gamma$-ray astronomy_, _Astronomy and Astrophysics_ 466 (2007) 1219 [astro-ph/0610959].
  * [14] H.E.S.S. collaboration, _Evidence of 100 TeV $\gamma$-ray emission from HESS J1702-420: A new PeVatron candidate_, 2106.06405.
  * [15] L. Mohrmann, A. Specovius, D. Tiziani, S. Funk, D. Malyshev, K. Nakashima et al., _Validation of open-source science tools and background model construction in $\gamma$-ray astronomy_, _Astron. Astrophys._ 632 (2019) A72 [1910.08088].
  * [16] C. Deil, R. Zanin, J. Lefaucheur, C. Boisson, B. Khelifi, R. Terrier et al., _Gammapy - A prototype for the CTA science tools_ , in _35th International Cosmic Ray Conference (ICRC2017)_ , vol. 301 of _International Cosmic Ray Conference_ , p. 766, Jan., 2017 [1709.01751].
  * [17] C. Deil, A. Donath, R. Terrier, J.E. Ruiz, J. King, QRemy et al., _gammapy/gammapy: v.0.18.2_ , Nov., 2020.  10.5281/zenodo.4701500.
  * [18] C. Planck, P. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont et al., _Planck early results. xix. all-sky temperature and dust optical depth from planck and iras. constraints on the" dark gas" in our galaxy_ , _Astronomy and Astrophysics_ 536 (2011) 19.
  * [19] G. Peron, F. Aharonian, S. Casanova, R. Yang and R. Zanin, _Probing the cosmic-ray density in the inner galaxy_ , _The Astrophysical Journal Letters_ 907 (2021) L11.
  * [20] E. Kafexhiu, F. Aharonian, A.M. Taylor and G.S. Vila, _Parametrization of gamma-ray production cross sections for p p interactions in a broad proton energy range from the kinematic threshold to PeV energies_ , _PhysRevD_ 90 (2014) 123014 [1406.7369].
  * [21] V. Zabalza, _naima: a python package for inference of relativistic particle energy distributions from observed nonthermal spectra_ , _Proc. of International Cosmic Ray Conference 2015_ (2015) 922 [1509.03319].

## Full Authors List: H.E.S.S. Collaboration

H. Abdalla1, F. Aharonian2,3,4, F. Ait Benkhali3, E.O. Angüner5, C. Arcaro6,
C. Armand7, T. Armstrong8, H. Ashkar9, M. Backes1,6, V. Baghmanyan10, V.
Barbosa Martins11, A. Barnacka12, M. Barnard6, R. Batzofin13, Y. Becherini14,
D. Berge11, K. Bernlöhr3, B. Bi15, M. Böttcher6, C. Boisson16, J. Bolmont17,
M. de Bony de Lavergne7, M. Breuhaus3, R. Brose2, F. Brun9, T. Bulik18, T.
Bylund14, F. Cangemi17, S. Caroff17, S. Casanova10, J. Catalano19, P.
Chambery20, T. Chand6, A. Chen13, G. Cotter8, M. Curyło18, H. Dalgleish1, J.
Damascene Mbarubucyeye11, I.D. Davids1, J. Davies8, J. Devin20, A. Djannati-
Ataï21, A. Dmytriiev16, A. Donath3, V. Doroshenko15, L. Dreyer6, L. Du
Plessis6, C. Duffy22, K. Egberts23, S. Einecke24, J.-P. Ernenwein5, S.
Fegan25, K. Feijen24, A. Fiasson7, G. Fichet de Clairfontaine16, G.
Fontaine25, F. Lott1, M. Füßling11, S. Funk19, S. Gabici21, Y.A. Gallant26, G.
Giavitto11, L. Giunti21,9, D. Glawion19, J.F. Glicenstein9, M.-H. Grondin20,
S. Hattingh6, M. Haupt11, G. Hermann3, J.A. Hinton3, W. Hofmann3, C.
Hoischen23, T. L. Holch11, M. Holler27, D. Horns28, Zhiqiu Huang3, D. Huber27,
M. Hörbe8, M. Jamrozy12, F. Jankowsky29, V. Joshi19, I. Jung-Richardt19, E.
Kasai1, K. Katarzyński30, U. Katz19, D. Khangulyan31, B. Khélifi21, S.
Klepser11, W. Kluźniak32, Nu. Komin13, R. Konno11, K. Kosack9, D. Kostunin11,
M. Kreter6, G. Kukec Mezek14, A. Kundu6, G. Lamanna7, S. Le Stum5, A.
Lemière21, M. Lemoine-Goumard20, J.-P. Lenain17, F. Leuschner15, C. Levy17, T.
Lohse33, A. Luashvili16, I. Lypova29, J. Mackey2, J. Majumdar11, D.
Malyshev15, D. Malyshev19, V. Marandon3, P. Marchegiani13, A. Marcowith26, A.
Mares20, G. Martí-Devesa27, R. Marx29, G. Maurin7, P.J. Meintjes34, M.
Meyer19, A. Mitchell3, R. Moderski32, L. Mohrmann19, A. Montanari9, C.
Moore22, P. Morris8, E. Moulin9, J. Muller25, T. Murach11, K. Nakashima19, M.
de Naurois25, A. Nayerhoda10, H. Ndiyavala6, J. Niemiec10, A. Priyana Noel12,
P. O’Brien22, L. Oberholzer6, S. Ohm11, L. Olivera-Nieto3, E. de Ona
Wilhelmi11, M. Ostrowski12, S. Panny27, M. Panter3, R.D. Parsons33, G. Peron3,
S. Pita21, V. Poireau7, D.A. Prokhorov35, H. Prokoph11, G. Pühlhofer15, M.
Punch21,14, A. Quirrenbach29, P. Reichherzer9, A. Reimer27, O. Reimer27, Q.
Remy3, M. Renaud26, B. Reville3, F. Rieger3, C. Romoli3, G. Rowell24, B.
Rudak32, H. Rueda Ricarte9, E. Ruiz-Velasco3, V. Sahakian36, S. Sailer3, H.
Salzmann15, D.A. Sanchez7, A. Santangelo15, M. Sasaki19, J. Schäfer19, H.M.
Schutte6, U. Schwanke33, F. Schüssler9, M. Senniappan14, A.S. Seyffert6,
J.N.S. Shapopi1, K. Shiningayamwe1, R. Simoni35, A. Sinha26, H. Sol16, H.
Spackman8, A. Specovius19, S. Spencer8, M. Spir-Jacob21, Ł. Stawarz12, R.
Steenkamp1, C. Stegmann23,11, S. Steinmassl3, C. Steppa23, L. Sun35, T.
Takahashi31, T. Tanaka31, T. Tavernier9, A.M. Taylor11, R. Terrier21, J. H.E.
Thiersen6, C. Thorpe-Morgan15, M. Tluczykont28, L. Tomankova19, M. Tsirou3, N.
Tsuji31, R. Tuffs3, Y. Uchiyama31, D.J. van der Walt6, C. van Eldik19, C. van
Rensburg1, B. van Soelen34, G. Vasileiadis26, J. Veh19, C. Venter6, P.
Vincent17, J. Vink35, H.J. Völk3, S.J. Wagner29, J. Watson8, F. Werner3, R.
White3, A. Wierzcholska10, Yu Wun Wong19, H. Yassin6, A. Yusafzai19, M.
Zacharias16, R. Zanin3, D. Zargaryan2,4, A.A. Zdziarski32, A. Zech16, S.J.
Zhu11, A. Zmija19, S. Zouari21 and N. Żywucka6.

1University of Namibia, Department of Physics, Private Bag 13301, Windhoek
10005, Namibia
2Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2,
Ireland
3Max-Planck-Institut für Kernphysik, P.O. Box 103980, D 69029 Heidelberg,
Germany
4High Energy Astrophysics Laboratory, RAU, 123 Hovsep Emin St Yerevan 0051,
Armenia
5Aix Marseille Université, CNRS/IN2P3, CPPM, Marseille, France
6Centre for Space Research, North-West University, Potchefstroom 2520, South
Africa
7Laboratoire d’Annecy de Physique des Particules, Univ. Grenoble Alpes, Univ.
Savoie Mont Blanc, CNRS, LAPP, 74000 Annecy, France
8University of Oxford, Department of Physics, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK
9IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
10Instytut Fizyki Ja̧drowej PAN, ul. Radzikowskiego 152, 31-342 Kraków, Poland
11DESY, D-15738 Zeuthen, Germany
12Obserwatorium Astronomiczne, Uniwersytet Jagielloński, ul. Orla 171, 30-244
Kraków, Poland
13School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue,
Braamfontein, Johannesburg, 2050 South Africa
14Department of Physics and Electrical Engineering, Linnaeus University, 351
95 Växjö, Sweden
15Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, D
72076 Tübingen, Germany
16Laboratoire Univers et Théories, Observatoire de Paris, Université PSL,
CNRS, Université de Paris, 92190 Meudon, France
17Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité,
CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4
Place Jussieu, F-75252 Paris, France
18Astronomical Observatory, The University of Warsaw, Al. Ujazdowskie 4,
00-478 Warsaw, Poland
19Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for
Astroparticle Physics, Erwin-Rommel-Str. 1, D 91058 Erlangen, Germany
20Université Bordeaux, CNRS/IN2P3, Centre d’Études Nucléaires de Bordeaux
Gradignan, 33175 Gradignan, France
21Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France
22Department of Physics and Astronomy, The University of Leicester, University
Road, Leicester, LE1 7RH, United Kingdom
23Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-
Strasse 24/25, D 14476 Potsdam, Germany
24School of Physical Sciences, University of Adelaide, Adelaide 5005,
Australia
25Laboratoire Leprince-Ringuet, École Polytechnique, CNRS, Institut
Polytechnique de Paris, F-91128 Palaiseau, France
26Laboratoire Univers et Particules de Montpellier, Université Montpellier,
CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France
27Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität
Innsbruck, A-6020 Innsbruck, Austria
28Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149,
D 22761 Hamburg, Germany
29Landessternwarte, Universität Heidelberg, Königstuhl, D 69117 Heidelberg,
Germany
30Institute of Astronomy, Faculty of Physics, Astronomy and Informatics,
Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland
31Department of Physics, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima-
ku, Tokyo 171-8501, Japan
32Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warsaw, Poland
33Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, D 12489
Berlin, Germany
34Department of Physics, University of the Free State, PO Box 339,
Bloemfontein 9300, South Africa
35GRAPPA, Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Science Park 904, 1098 XH Amsterdam, The Netherlands
36Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan,
Armenia