# An Adaptive Task-Related Component Analysis Method for SSVEP recognition

Vangelis P. Oikonomou

###### Abstract

Steady-state visual evoked potential (SSVEP) recognition methods are equipped
with learning from the subject’s calibration data, and they can achieve extra
high performance in the SSVEP-based brain–computer interfaces (BCIs), however
their performance deteriorate drastically if the calibration trials are
insufficient. This study develops a new method to learn from limited
calibration data and it proposes and evaluates a novel adaptive data-driven
spatial filtering approach for enhancing SSVEPs detection. The spatial filter
learned from each stimulus utilizes temporal information from the
corresponding EEG trials. To introduce the temporal information into the
overall procedure, an multitask learning approach, based on the bayesian
framework, is adopted. The performance of the proposed method was evaluated
into two publicly available benchmark datasets, and the results demonstrated
that our method outperform competing methods by a significant margin.

###### keywords:

Steady State Visual Evoked Potentials , EEG , Task related Component Analysis
, Multi-task learning , Spatial Filtering , Brain Computer Interfaces

††journal:

[inst1]organization=Information Technologies Institute, Centre for Research
and Technology Hellas,addressline=Thermi-Thessaloniki, city=Thessaloniki,
postcode=57001, country=Greece

##  1 Introduction

A Brain Computer Interface (BCI) system is a device that it translates human
brain activity into artificially generated control signals, providing us with
an alternative communication medium, other than physical communication, which
can be used to help people with motor disabilities [1], to augment
communication abilities of healthy individuals,[2], for entertainment [2],
and, for neuromarketing purposes[3]. Brain activity can be measured with
various specialized devices such as MRI scanners and electroencephalograms
(EEG)-based devices, where, EEG devices is widely used since the required
equipment is simple and inexpensive. EEG-based BCI systems utilize various
brain responses such as motor imagery and visual responses, from which, the
use of Steady State Visual Evoked Potentials (SSVEPs) have attracted the
attention of many researchers due to their lower training requirements for the
end-user and higher information transfer rates [4, 5]. When an individual is
looking into a visual stimulus, which is flashing as a fixed frequency, then a
brain response is revealed in occipital and occipital - parietal areas of
individual’s brain which is called SSVEP response[6]. A SSVEP response
contains sinusoidal components which are related to the fundamental frequency
of the visual stimulus as well as its harmonics. The overarching goal of a
SSVEP BCI system is to detect the different frequency components corresponding
to the visual stimuli and translate them into commands, by using an EEG-based
pattern recognition algorithms.

The recognition of SSVEP responses involves the use of Machine Learning (ML)
algorithms. Linear classifiers such as Support Vector Machines (SVMs) and the
Linear Discriminant Analysis (LDA) have been used to detect SSVEPs [7]. In
addition, in [8] the use of Multivariate Linear Regression (MLR) was proposed
to learn discriminative features for improving SSVEP classification, while, in
[9] kernel - based extensions of MLR were proposed using SSVEP-related kernels
as an integral part of the Sparse Bayesian Learning (SBL) framework.
Furthermore, Deep Learning (DL) approaches using Convolutional Neural Networks
(CNN), based on time frequency analysis, are used to discriminate SSVEP
responses [10, 11].

However, SSVEP responses present specific frequency and spatial
characteristics, hence methods utilizing these characteristics have been
proposed. More specifically, methods based on Power Spectrum Density Analysis
(PSDA) were widely used for frequency detection [12], where the target
frequency is assigned to the frequency corresponding to maximum value of PSD.
However, the efficient calculation of PSD requires a relatively large time
window, and it is sensitive to noise [13, 14]. To avoid the above shortcomings
of PSDA, spatial filtering approaches have been proposed. In [14] the Minimum
Energy Combination (MEC) method has been proposed while in [13] the Canonical
Correlation Analysis (CCA) method was introduced. Both methods use sinusoids
waves as reference templates and they solve an optimization problem, based on
multi-channel SSVEP data, in order to obtain optimal spatial filters. Finally,
extensions of CCA have been proposed in [15, 16, 17, 18, 19, 20], extracting
the subject-specific and task-related information from the individual
calibration data and reduce the effect of spontaneous background EEG
activities.

From spatial filtering methods, the task-related component analysis
(TRCA)-based method [18] shows great potential since it has achieved superior
performance among various spatial filtering methods. The core idea of TRCA is
to acquire the spatial filters by strengthening the task-related SSVEP
components and suppressing the noise. TRCA-based methods are followed by the
target detection step where the similarity between the filtered test signal
and the filtered template is calculated via the correlation coefficient. All
spatial filtering - based methods are based on the basic (generalized)
eigenvalue problem [21, 22]. However, difference between various approaches
can be observed and these difference are reflected to the way that the
matrices, involved in the eigenvalue problem, are constructed [22]. In [23]
Correlated Component Analysis (CORCA)assumes that the task related component
is shared among the subjects by adopting a transfer learning procedure in the
construction of covariance matrices. While, in [24] a task discriminant
component analysis was applied which involves the construction of within and
between SSVEP targets covariance matrices.

One defect of TRCA is that it can only deal with limited noise components. For
other kinds of noise, such as locally occurring noises that have the same
profile, the TRCA-based method is powerless[25, 26]. Additionally, its
performance deteriorate drastically if the calibration trials are
insufficient. To deal with more general noises and the number of trials, we
introduce a novel adaptive time-domain filter resulting in more reliable
similarity measurement. By introducing the temporally-based filter into the
objective function of the TRCA-based method we construct a time filter that
acts together with the spatial filter to suppress more general noises.
Furthermore, the filter adapts to the statistical properties of SSVEP trials.

The rest of this paper is organized as follows. In section 2, first, we
provide a short description of CCA and TRCA concentrating on the mathematical
formulation of them, and then, we describe our approach for SSVEP recognition.
Section 3, we describe the SSVEP datasets that are used in our study, and then
we provide details about our experiments and the performance of our method.
Also, a comparison with competing methods is provided. Finally, a short
discussion and some concluding remarks are are provided in section 4.

##  2 Materials and Methods

###  2.1 Problem Description

When a SSVEP experiment is taken place, the subject is seated in front of a
screen where visual stimuli are flashing in different frequencies, During the
experiment raw EEG data are collected in order to calibrate the overall
system. The segmentation of raw EEG data (using event triggers), results into
a set of trials for each visual stimulus (or class). Using these EEG trials
the experimenter can calibrated the BCI system (for example, by training the
classifier). Let us assume that the SSVEP dataset is a collection of multi-
channel EEG trials
$\\{\mathbf{X}_{1}^{(s)},\mathbf{X}_{2}^{(s)},\cdots,\mathbf{X}_{M}^{(s)}\\}_{s=1}^{N_{s}}$
for each participant, where $M$ is the number of trials of a SSVEP target,
$(s)$ is the index of the SSVEP target, $N_{s}$ is the number of SSVEP targets
(or classes). Each $\mathbf{X}_{m}^{(s)},m=1,\cdots,M,s=1,\cdots,N_{s}$ is a
matrix of $N_{ch}\times N_{t}$, where $N_{ch}$ is the number of channels and
$N_{t}$ the number of samples. Additionally, we assume that the multi-channel
EEG signals are centralized since in practise the EEG trials are bandpass
filtered or detrended.

###  2.2 Canonical Correlation Analysis (CCA)

Spatial filtering attempts to maximize the SNR between the raw EEG data and
the spatial filtered version of them. In typical cases, such as bipolar
combination or laplacian filtering, the spatial filters are determined
manually. However, this approach does not take into account any prior
knowledge about SSVEPs or any subject-specific information. One of the first
approach that take into consideration the structure of SSVEPs was based on
Canonical Correlation Analysis (CCA)[13]. The CCA is a multivariate
statistical method attempting to discover underlying correlations between two
sets of data[13, 27]. These two sets of data is assumed to be only a different
view (or representation) of the same original (hidden) data. More
specifically, CCA finds a linear projection for each set such that these two
set are maximally correlated in the hidden (dimensionality-reduced) space.

In the SSVEP problem these two views is the test EEG trial
$\mathbf{X}_{m}^{(s)}$ and the reference templates for $s$-th stimulus
$\mathbf{Y}_{(f_{s})}$, where

$\mathbf{Y}_{(f_{s})}=\begin{bmatrix}\sin(2\pi\cdot 1\cdot f_{s}t)\\\
\cos(2\pi\cdot 1\cdot f_{s}t)\\\ \vdots\\\ \sin(2\pi\cdot N_{h}\cdot
f_{s}t)\\\ \cos(2\pi\cdot N_{h}\cdot f_{s}t)\end{bmatrix}^{\top}$

$\mathbf{Y}_{(f_{s})}\in\mathbb{R}^{N_{t}\times 2N_{h}}$, $f_{s}$ is the
frequency of $s$-th stimulus

Typically, CCA methods maximize the linear correlation between the projections
$\mathbf{w}_{s}^{T}\mathbf{X}_{m}^{(s)}$ and
$\mathbf{v}_{s}^{T}\mathbf{Y}_{f_{s}}$, where
$\mathbf{w}_{s}\in\mathbb{R}^{N_{ch}}$ and
$\mathbf{v}_{s}\in\mathbb{R}^{N_{t}}$. At the end, we solve the following
optimization problem:

$\displaystyle\max\rho_{s}=\max_{\mathbf{w}_{s},\mathbf{v}_{s}}\frac{\mathbf{w}_{s}^{\top}\mathbf{X}_{m}^{(s)}\mathbf{Y}_{f_{s}}^{\top}\mathbf{v}_{s}}{\sqrt{\mathbf{w}_{s}^{\top}\mathbf{X}_{m}^{(s)}(\mathbf{X}_{m}^{(s)})^{\top}\mathbf{w}_{s}\mathbf{v}_{s}^{\top}\mathbf{Y}_{f_{s}}\mathbf{Y}_{f_{s}}^{\top}\mathbf{v}_{s}}}$
(1)

Since $\rho_{s}$ is invariant to the scaling of $\mathbf{w}_{s}$ and
$\mathbf{v}_{s}$ the above optimization problem can be also formulated as the
following generalized eigenvalue problem:

$\displaystyle\mathbf{X}_{m}^{(s)}\mathbf{Y}_{f_{s}}^{\top}(\mathbf{Y}_{f_{s}}\mathbf{Y}_{f_{s}}^{\top})^{-1}\mathbf{Y}_{f_{s}}(\mathbf{X}_{m}^{(s)})^{\top}\mathbf{w}_{s}=\lambda_{s}\mathbf{X}_{m}^{(s)}(\mathbf{X}_{m}^{(s)})^{\top}\mathbf{w}_{s}$
(2)

where $\lambda_{s}$ is the eigenvalue corresponding to the eigenvector
$\mathbf{w}_{s}$.

In order to find the stimulus of test EEG trial $\mathbf{X}_{m}^{(s)}$, that
the subject intends to select, we find $\rho_{s}$ for all available stimuli
and the stimulus-target, $c$, is then identified by finding the index of the
maximum feature among $N_{s}$ features: $c=\arg\max_{s}\\{\rho_{s}\\}$. It
must be observed here that there is no need for training (or calibration)
since the templates $\mathbf{Y}_{f_{s}}$ are artificially generated.

###  2.3 Task Related Component Analysis (TRCA)

Task-related component analysis (TRCA) enhances reproducibility of SSVEPs
across multiple trials and the intuition of TRCA is to maximize the
reproducibility of SSVEP target-related components after spatial filtering.
More specifically, the TRCA method find the spatial filters $\mathbf{w}_{s}$
by solving a generalized eigenvalue problem which is described by the
following equation:

$\displaystyle\max_{\mathbf{w}_{s}}\frac{\mathbf{w}_{s}^{\top}AA^{\top}\mathbf{w}_{s}}{\mathbf{w}_{s}^{\top}BB^{\top}\mathbf{w}_{s}}$
(3)

where $A=\frac{1}{M}\sum_{m=1}^{M}\mathbf{X}_{m}^{(s)}$, and $B$ is a
concatenated matrix contains all trials of $s$-th stimulus,
$B=[\mathbf{X}_{1}^{(s)},\mathbf{X}_{2}^{(s)},\cdots\mathbf{X}_{M}^{(s)}]$.

In order to find the target of the test trial, $\mathbf{X}_{test}$ we apply
the following discriminant function:

$c=\arg\max_{s}\\{corr(\mathbf{w}_{s}^{\top}\mathbf{X}_{test},\mathbf{w}_{s}^{\top}A)\\}$
(4)

where $corr(\cdot,\cdot)$ denotes the Pearson’s correlation coefficient.

###  2.4 Adaptive Task-Related Component Analysis (adTRCA)

In our work we propose a new generalized eigenvalue problem for SSVEP
detection which is described by the following equation:

$\displaystyle\max_{\mathbf{w}_{s}}\frac{\mathbf{w}_{s}^{\top}ACA^{\top}\mathbf{w}_{s}}{\mathbf{w}_{s}^{\top}BDB^{\top}\mathbf{w}_{s}}$
(5)

where $C$ and $D$ are ”filtering” matrix that acts on the time dimension of
trials. The matrices $C$ and $D$ can be defined using various approaches and
their goal is to remove noise in time domain.

In our study we make some critical assumptions about the generation model of
SSVEP responses, which affect the data analysis procedure. More specifically,
SSVEP responses contains strong sinusoids components [13], hence the SSVEP
signal in each channel is modeled as a linear combination of sinusoids
described the following matrix:

$\boldsymbol{\Phi}=[\mathbf{Y}_{(f_{1})}\mathbf{Y}_{(f_{2})}\cdots\mathbf{Y}_{(f_{N_{s}})}]\in\mathbb{R}^{N_{t}\times(2N_{s}N_{h})}.$

Additionally, SSVEP responses belonging to the same visual stimulus share
common components. From the above we can observe that the generation of SSVEP
responses can be modeled as multiple regression tasks that share common
information.

EEG trials from the $s$-th stimulus are collected in matrix
$\mathbb{B}=[{\mathbf{X}_{1}^{(s)}}^{\top},{\mathbf{X}_{2}^{(s)}}^{\top},\cdots{\mathbf{X}_{M}^{(s)}}^{\top}]^{\top}$,
$\mathbb{B}\in\mathbb{R}^{N_{t}\times(N_{ch}N_{s})}$, where each column of
$\mathbb{B}$ contains the data from one channel or each column of $\mathbb{B}$
contains the data from one task. Hence we have
$\boldsymbol{y}_{i}\in\mathbb{R}^{N_{t}\times 1},i=1,\cdots,N_{ch}N_{s}$ tasks
(the i-th column of $\mathbb{B}$). Each learning task can be described by the
following linear regression model:

$\boldsymbol{y}_{i}=\boldsymbol{\Phi}\boldsymbol{w}_{i}+\boldsymbol{e}_{i}$
(6)

where $\boldsymbol{w}_{i}$ $2N_{s}N_{h}\times 1$ vector of weights (or
parameters), and, $\boldsymbol{e}_{i}$ $N_{t}\times 1$ vector of noise coming
from a zero mean Gaussian random variable with unknown precision (inverse
variance) $a_{0}$. We can observe that each of the mapping yields a
corresponding regression task, and performing multiple such learning tasks has
been referred to as multitask learning [28], which aims at sharing information
effectively among multiple related tasks. In a more abstract view of our
problem we can see that each learning task is a linear regression problem, and
sinusoids components from one regression task affect the fitting procedure of
another regression task.

The likelihood function for parameters $\boldsymbol{w}_{i}$ and $a_{0}$ is
given by:

$p(\boldsymbol{y}_{i}|\boldsymbol{w}_{i},a_{0})=(2\pi
a_{0})^{-\frac{N_{t}}{2}}\exp\Big{\\{}-\frac{a_{0}}{2}\|\boldsymbol{y}_{i}-\boldsymbol{\Phi}\boldsymbol{w}_{i}\|^{2}_{2}\Big{\\}}$
(7)

The parameters of a regression task, $\boldsymbol{w}_{i}$, are assumed to be
drawn from a product of zero-mean Gaussian distributions that are shared by
all tasks. Letting $w_{i,j}$ be the $j$-th parameters for $i$-th task then we
have:

$p(\boldsymbol{w}_{i}|\boldsymbol{a})=\prod_{j=1}^{2N_{s}N_{h}}{\cal
N}(w_{i,j}|0,a^{-1}_{i})$ (8)

where the hyperparameters
$\boldsymbol{a}=\\{a_{j}\\}_{j=1,2,\cdots,2N_{s}N_{h}}$ are shared among
$N_{ch}N_{s}$ regression tasks, hence, data from all regression tasks
contribute to learning these hyperparameters. To promote sparsity over
parameters, we place Gamma priors over hyperparameters $\boldsymbol{a}$[29,
28]. Also, the same type of prior is placed over noise precision $a_{0}$.

$\displaystyle p(a_{0}|\alpha,\beta)$ $\displaystyle=Ga(a_{0}|\alpha,\beta)$
$\displaystyle=\frac{\beta^{\alpha}}{\Gamma(\alpha)}a_{0}^{\alpha-1}\exp\Big{\\{}-\beta
a_{0}\Big{\\}}$ (9)

$p(\boldsymbol{a}|c,d)=\prod_{j=1}^{2N_{s}N_{h}}Ga(a_{j}|c,d)$ (10)

In addition, we can observed here, that noise properties are shared among
different tasks (i.e. the noise vectors in Eq. (6) are drawn from the same
Gaussian distribution). Finally, it must be noted that we have an hierarchical
model, and these types of models are natural to be ”dealt” within the bayesian
framework.

Given hyperparameters $\boldsymbol{a}$ and noise precision $a_{0}$, we can
apply Bayes theorem to find the posterior distribution over
$\boldsymbol{w}_{i}$, which is a Gaussian distribution:

$\displaystyle p(\boldsymbol{w}_{i}|\boldsymbol{y}_{i},\boldsymbol{a},a_{0})$
$\displaystyle=\frac{p(\boldsymbol{y}_{i}|\boldsymbol{w}_{i},a_{0})p(\boldsymbol{w}_{i})|\boldsymbol{a})}{\int
p(\boldsymbol{y}_{i}|\boldsymbol{w}_{i},a_{0})p(\boldsymbol{w}_{i})|\boldsymbol{a})d\boldsymbol{w}_{i}}$
$\displaystyle={\cal
N}\Big{(}\boldsymbol{w}_{i}\Big{|}\boldsymbol{\mu}_{i},\boldsymbol{\Sigma}_{i}\Big{)}$
(11)

where

$\displaystyle\boldsymbol{\mu}_{i}$
$\displaystyle=a_{0}\boldsymbol{\Sigma}_{i}\boldsymbol{\Phi}^{T}\boldsymbol{y}_{i}$
(12) $\displaystyle\boldsymbol{\Sigma}_{i}$
$\displaystyle=\Big{(}a_{0}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}+\boldsymbol{A}\Big{)}^{-1}$
(13)

and $\boldsymbol{A}=diag(a_{1},a_{2},\cdots,a_{M})$.

In order to find hyperparameters $\boldsymbol{a}$ and promote sparsity in
parameters, the type-II Maximum Likelihood procedure is adopted [29, 30],
where the objective is to maximize the marginal likelihood (or its logarithm).
Also, a similar procedure is followed for the noise precision. The marginal
likelihood ${\cal L}(\boldsymbol{a},a_{0})$ is given by:

$\displaystyle{\cal L}(\boldsymbol{a},a_{0})=\sum_{i=1}^{L}\log\int
p(\boldsymbol{y}_{i}|\boldsymbol{w}_{i},a_{0})p(\boldsymbol{w}_{i}|\boldsymbol{a})d\boldsymbol{w}_{i}$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{L}\Big{(}N_{i}\log(2\pi)+\log|\boldsymbol{C}_{i}|+\boldsymbol{y}_{i}^{T}\boldsymbol{C}_{i}^{-1}\boldsymbol{y}_{i}\Big{)}$
(14)

where
$\boldsymbol{C}_{i}=a_{0}^{-1}\boldsymbol{I}+\boldsymbol{\Phi}\boldsymbol{A}\boldsymbol{\Phi}^{T}$

Differentiating ${\cal L}(\boldsymbol{a},a_{0})$ with respect to
$\boldsymbol{a}$ and $a_{0}$ and setting the results into zero [29, 30, 28]
(after some algebraic manipulations) we obtain:

$\displaystyle a_{j}^{(new)}$
$\displaystyle=\frac{(N_{ch}N_{s})-a_{j}\sum_{i=1}^{(N_{ch}N_{s})}\Sigma_{i,(j,j)}}{\sum_{i=1}^{(N_{ch}N_{s})}\mu_{i,j}},j=1,2,\cdots,2N_{s}N_{h}$
(15) $\displaystyle a_{0}^{(new)}$
$\displaystyle=\frac{\sum_{i=1}{N_{ch}N_{s}}\Big{(}N_{t}-2N_{s}N_{h}+\sum_{j=1}^{2N_{s}N_{h}}a_{j}\Sigma_{i,(j,j)}\Big{)}}{\sum_{i=1}^{N_{ch}N_{s}}\|\boldsymbol{y}_{i}-\boldsymbol{\Phi}_{i}\boldsymbol{\mu}_{i}\|^{2}_{2}}$
(16)

where $\mu_{i,j}$ the $j$-th element of $\boldsymbol{\mu}_{i}$ and
$\Sigma_{i,(j,j)}$ the $j$-th diagonal element of covariance matrix
$\boldsymbol{\Sigma}_{i}$. The above analysis suggests an iterative algorithm
that iterates between Eqs. (12), (13), (15) and (16), until a convergence
criterion is satisfied. Also, the same algorithm can be derived by adopting
the EM framework and treating parameters $\boldsymbol{w}_{i}$ as hidden
variables[29]. Finally, based on the above bayesian formulation, we can derive
a fast version of the above algorithm. The fast version provides an elegant
treatment of feature vectors by constructing adaptively the matrix
$\boldsymbol{\Phi}$ through three basic operators: addition, deletion and re-
estimation. More information on this subject can be found in [29, 28].

Now, SSVEP components in each task can be represented as:

$\boldsymbol{\hat{y}}_{i}=\boldsymbol{\Phi}\boldsymbol{\mu}_{i},i=1,\cdots,N_{ch}N_{s}$

rearranging filtered EEG signals, $\boldsymbol{\hat{y}}_{i}$, each filtered
EEG trial is represented as:
$\mathbf{X^{f}}_{m}^{(s)}=\mathbf{X}_{m}^{(s)}a_{0}\boldsymbol{\Phi}\boldsymbol{\Sigma}_{i}\boldsymbol{\Phi}^{\top}$.
Due to filtered trials we find the spatial filters $\mathbf{w}_{s}$ by solving
the following generalized eigenvalue problem:

$\displaystyle\max_{\mathbf{w}_{s}}\frac{\mathbf{w}_{s}^{\top}A_{f}A_{f}^{\top}\mathbf{w}_{s}}{\mathbf{w}_{s}^{\top}B_{f}B_{f}^{\top}\mathbf{w}_{s}}$
(17)

where $A_{f}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{X^{f}}_{m}^{(s)}$, and $B_{f}$
is a concatenated matrix contains all trials of $s$-th stimulus,
$B_{f}=[\mathbf{X^{f}}_{1}^{(s)},\mathbf{X^{f}}_{2}^{(s)},\cdots\mathbf{X^{f}}_{M}^{(s)}]$.
The above generalized eigenvalue problem can be connected by that of Eq. 5.
After some algebraic manipulations Eq. 17 can be written as:

$\displaystyle\max_{\mathbf{w}_{s}}\frac{\mathbf{w}_{s}^{\top}ACA^{\top}\mathbf{w}_{s}}{\mathbf{w}_{s}^{\top}BDB^{\top}\mathbf{w}_{s}}$
(18)

where
$C=(a_{0}\boldsymbol{\Phi}\boldsymbol{\Sigma}_{i}\boldsymbol{\Phi}^{\top})(a_{0}\boldsymbol{\Phi}\boldsymbol{\Sigma}_{i}\boldsymbol{\Phi}^{\top})^{\top}$
and $D=\begin{bmatrix}C&\cdots&0\\\ \vdots&\ddots&\vdots\\\ 0&\cdots&C\\\
\end{bmatrix}.$ We can observe an interesting connection between the proposed
method and the TRCA. When the $C=\mathbf{I}$, where $\mathbf{I}$ is the
unitary matrix, the proposed approach degraded to the TRCA method. We see that
the TRCA method is a limiting case of the proposed method. In addition, we can
observe that matrices $C$ and $D$ acts on the time dimension of the EEG
trials, hence time samples are differently treated according to the time
dimension rather than equally weighted. Also, we can observe that filters,
represented by matrix $C$, are adapted to the statistical properties of EEG
trials. Finally, after finding the spatial filters, to find the target of the
test trial, $\mathbf{X}^{f}_{test}$ we apply the following discriminant
function:

$c=\arg\max_{s}\\{corr(\mathbf{w}_{s}^{\top}\mathbf{X}^{f}_{test},\mathbf{w}_{s}^{\top}A_{f})\\}$
(19)

###  2.5 Ensemble case

According to the previous discriminant rule, described by Eq. 19, we can
observe that to calculate the similarity of the test trial to the stimulus $s$
first we apply spatial filter $\mathbf{w}_{s}$. However, since we have enough
calibration trials we can obtained $N_{s}$ spatial filters for each stimulus
[18, 23]. Hence, we can extend our method using an ensemble approach, similar
to [18, 23], where all spatial filters for stimulus $s$ are concatenated to
create an ensemble spatial filter $\mathbf{W}_{s}\in\mathbb{R}^{N_{ch}\times
N_{s}}$. Now, in order to find the target of the test trial,
$\mathbf{X}^{f}_{test}$ we apply the following discriminant function:

$c=\arg\max_{s}\\{corr(\mathbf{W}_{s}^{\top}\mathbf{X}^{f}_{test},\mathbf{W}_{s}^{\top}A_{f})\\},$
(20)

where, now, the function $corr(\cdot)$ depicts correlation between matrices.

##  3 Results

In our study two, publicly available, SSVEP datasets are used to evaluate our
method. We call these datasets, the Speller dataset and the EPOC dataset, and
we provide a description of them in the following paragraphs.

Speller dataset [31]: The Speller dataset [31] contains 40-target visual
stimuli were presented on a 23.6-in LCD monitor. Thirty five healthy subjects
(with normal or corrected-to-normal vision) were recruited for the SSVEP
experiment. Furthrmore, eight of them had experience with using a SSVEP-based
BCI speller. EEG signals were recorded with 64 channels according to an
extended 10–20 system. In our experiments we have used electrodes from the
occipital and parietal-occipital areas ($P_{z}$, $PO_{5}$, $PO_{3}$, $PO_{z}$,
$PO_{4}$, $PO_{6}$, $O_{1}$, $O_{z}$, and $O_{2}$) . During the experiment,
each subject completes 6 blocks, where in each block the subject was looking
at one of the visual stimuli, indicated by the stimulus program, in a randomg
order for 5s and he/she complete 40 trials corresponding to all 40 targets.
Using event triggers, EEG trials were extracted and down-sampled to 250Hz.
Futhermore, there trials have been band-pass filterd to 7-90Hz with an
infinite impulse response (IIR) filter using the filtfilt function from
MATLAB. Additionally, a delay of 140ms in the visual system was considered in
our experiments[31].

EPOC dataset [7]: EEG trials, from 11 subjects executing a SSVEP-based
experimental protocol, were acquired. The frequencies of visual stimulation
were: 6.66Hz, 7.50Hz, 8.57Hz, 10.00Hz and 12.00Hz. EEG trials were recorded
using the Emotiv Epoc device, with 14 wireless channels and a sampling rate of
128 Hz. Each subject was asked to gaze at one of the visual stimuli, indicated
by the stimulus program, in a random order, and complete 20 trials for each of
the five targets. The EEG data have been band-pass filtered from 5Hz to 45Hz.
More information about this dataset can be found in [7] and
https://physionet.org/content/mssvepdb/1.0.0/. The visual latency is an
important factor in a SSVEP BCI system. An accurate estimation of the visual
latency ensures that the extracted data epochs only contain SSVEP responses to
the stimulation. In the Speller dataset the visual latency is defined using
the classification accuracy[31], while, in the EPOC dataset an alignment
procedure is used to extract data epochs containing the SSVEP responses.

###  3.1 Performance metrics

To calculate the performance metrics we adopt a leave-one-out cross-validation
scheme where we use $B-1$ blocks for training and 1 block for testing, similar
to [17, 8, 31, 9]. The performance of the examined algorithms was calculated
by using the following metrics:

Classification Accuracy: Classification Accuracy is the ratio between the
number of correctly classified trials to the total number of trials.

Information Transfer Rate (ITR): When a BCI system is designed, we are
interested, besides the accuracy, to the number of classes offered, as well
as, to the time required for the classification process. These two factors are
particularly important in the design of BCI systems since they determine the
amount of information that can be transmitted. Hence, we use an additional
metric called Information Transfer Rate (ITR). ITR is used to quantify the
rate of information transmitted by the BCI system and is measured by the
following equation:

$ITR=\Bigg{(}\log_{2}K+P\cdot\log_{2}P+(1-P)\log_{2}\Big{[}\frac{1-P}{K-1}\Big{]}\Bigg{)}\cdot\Bigg{(}\frac{60}{T}\Bigg{)}$
(21)

with $K$ being the number of classes, $P$ the classification accuracy, and $T$
the time in seconds required for the classification process to complete.

We compare the proposed spatial filtering method with two spatial filtering
approaches: the CCA[13] and the TRCA[18], and with two ML methods: the MLR
approach [8] and the Graph-based Sparse Representations Classification (MLR-
SRC) [32]. We have calculate the above metrics with respect the time window of
the test EEG trial, the number of EEG channels and the number of training
trials.

###  3.2 Performance Comparison versus the Time Window Length

In the first series of experiments we follow a typical analysis of SSVEP
datasets by examine the performance of methods with respect to the time
window. More specifically, the performance of methods are evaluated for
variable time from 0.5sec to 4sec with step of 0.5sec. In addition for the
Speller dataset we use the 9 channels from occipital and parietal - occipital
areas (Pz, PO5, PO3, POz, PO4, PO6, O1, Oz, O2), while, for the EPOC dataset
we use all available 14 channels, covering the entire brain. In Fig. 1 we
provide the obtained results for all comparative methods for the two datasets
using the basic channel configuration. We can observed the adTRCA and TRCA
methods provides better results from all other methods in both datasets. Also,
we can observed that between the adTRCA and TRCA, the adTRCA method has
marginally better detection accuracy in the Speller dataset, and significant
better detection accuracy in the EPOC dataset. Similar conclusions can be
drawn with respect to the ITR (see Tables 1 and 2). We can observed that for
most TWs, the adTRCA methods provides the best ITR among all methods. However,
we must point out, that for the Speller dataset the best ITR is achieved at
TW=0.5s from the adTRCA method, while, for the EPOC dataset, the best ITR
value is achieved at TWs=0.5s from the TRCA method.

|
---|---
(a) | (b)
Figure 1: Average Classification over all subjects (a) for the Speller dataset and (b) the EPOC dataset 14 using the basic configuration with respect to the EEG channels. In both cases, the time window ranges from 0.5s to 4s (0.5s interval). * indicates statistically significant difference between the TRCA and adTRCA methods, using paired sample t-test for Speller dataset and Wilcoxon signed rank test for the EPOC dataset ($p<0.05$). Table 1: ITR on Speller dataset - 9 channels TW | CCA | MLR | MLR-SRC | TRCA | adTRCA
---|---|---|---|---|---
0.5 | 27.3581 | 141.5439 | 98.0103 | 270.9573 | 271.3097
1 | 101.8150 | 137.3785 | 123.1899 | 260.0963 | 260.3216
1.5 | 122.8637 | 118.8840 | 116.4590 | 189.4191 | 190.8013
2 | 116.7986 | 98.5886 | 98.1738 | 149.0613 | 148.8653
2.5 | 101.4344 | 85.9023 | 85.3361 | 122.2336 | 122.1442
3 | 90.6986 | 79.0954 | 77.8773 | 103.1242 | 103.3462
3.5 | 81.1131 | 71.3472 | 71.3823 | 89.2204 | 89.0756
4 | 72.5963 | 65.1814 | 65.0779 | 78.2199 | 78.3041
Table 2: ITR on EPOC dataset - 14 channels TW | CCA | MLR | MLR-SRC | TRCA | adTRCA
---|---|---|---|---|---
0.5 | 13.7550 | 61.4288 | 81.6831 | 113.1879 | 109.7278
1 | 16.4944 | 45.1581 | 49.0298 | 67.8476 | 70.7059
1.5 | 13.3129 | 35.7231 | 37.9279 | 52.7534 | 55.7158
2 | 13.6286 | 31.7665 | 31.8755 | 44.0878 | 48.9809
2.5 | 12.4521 | 29.7089 | 28.1211 | 38.7167 | 42.2460
3 | 12.0830 | 25.1981 | 26.2000 | 35.0314 | 38.0493
3.5 | 12.1995 | 23.8428 | 24.4427 | 31.7552 | 34.8263
4 | 12.6687 | 23.6196 | 23.6422 | 29.8733 | 32.1765

###  3.3 Performance Comparison using the minimal number of EEG channels

In the second series of experiments we have used the minimal number of EEG
channels to evaluate the performance of methods. These EEG channels covering
the occipital area depending of the EEG device that had been used in each
dataset. More specifically, for the Speller dataset O1, O2 and Oz EEG channels
are used, while, for the EPOC dataset, we use O1 and O2 channels. This study
could correspond to cases where we are not able to use high density EEG
devices such as in BCI applications outside a controlled environment.

In Fig. 2 we provide the obtained results for all comparative methods for the
two datasets using the minimal channels’ configuration. We can observed the
adTRCA and TRCA methods provides better results from all other methods in both
datasets. Also, we can observed that between the adTRCA and TRCA, the adTRCA
method has better significant better detection accuracy in both datasets. In
Tables 3 and 4, we provide the results with respect to the ITR measure for all
methods. We can observed that for all TWs, besides ones, the adTRCA methods
provides the best ITR among all methods. However, while for the Speller
dataset the best ITR is achieved at TW=0.5s from the adTRCA method, for the
EPOC dataset, the best ITR value is achieved at the same TW from the TRCA
method.

|
---|---
(a) | (b)
Figure 2: Average Classification over all subjects (a) for the Speller dataset (b) for the EPOC dataset using the EEG channels covering the occipital areas. In both cases, the time window ranges from 0.5s to 4s (0.5s interval). * indicates statistically significant difference between the TRCA and adTRCA methods, using paired sample t-test for Speller dataset and Wilcoxon signed rank test for the EPOC dataset ($p<0.05$). Table 3: ITR on Speller dataset - 3 channels TW | CCA | MLR | MLR-SRC | TRCA | adTRCA
---|---|---|---|---|---
0.5 | 26.6513 | 82.8903 | 39.9318 | 162.8981 | 163.3228
1 | 71.6644 | 93.4390 | 74.2354 | 189.5521 | 190.6936
1.5 | 89.8807 | 85.1997 | 82.2758 | 148.7973 | 150.0271
2 | 93.5801 | 76.2840 | 73.1798 | 121.8152 | 123.6819
2.5 | 85.4184 | 69.7438 | 67.2738 | 102.4240 | 103.7595
3 | 80.0041 | 67.1743 | 66.4061 | 90.6681 | 91.6374
3.5 | 73.1382 | 62.4800 | 62.2848 | 79.6201 | 81.0671
4 | 67.1804 | 58.2207 | 57.9068 | 71.3779 | 72.9298
Table 4: ITR on EPOC dataset - 2 channels TW | CCA | MLR | MLR-SRC | TRCA | adTRCA
---|---|---|---|---|---
0.5 | 20.4758 | 17.7383 | 41.2960 | 55.2841 | 54.7309
1 | 20.8587 | 13.2992 | 29.1648 | 34.8970 | 38.5391
1.5 | 20.2541 | 16.2247 | 26.1393 | 29.4451 | 31.4028
2 | 18.8570 | 16.9955 | 21.9714 | 27.6870 | 31.5392
2.5 | 16.7189 | 17.8696 | 20.5254 | 23.4503 | 27.1407
3 | 15.7776 | 16.6921 | 18.1873 | 20.4231 | 25.5093
3.5 | 14.5984 | 15.0986 | 16.8087 | 19.7346 | 23.7536
4 | 14.6751 | 15.3353 | 17.0919 | 19.9582 | 23.4404

###  3.4 Performance Comparison versus the Number of training trials

In the last series of experiments, we investigate how the performance of TRCA
and adTRCA are affected by the number of training trials when the time window
is 1sec. The obtain results are provided in Tables 5 and 6. The proposed
method clearly has better accuracy from TRCA. Additionally, the proposed
scheme achieves the best performance among them and the margin provided by it
is more distinct when a smaller number of training blocks are utilized.
Especially, in the case Speller dataset, when only three training blocks are
utilized, the proposed scheme achieves a classification accuracy of 63%, while
the TRCA method achieve an accuracy of 61%. Furthermore, in the case of EPOC
dataset, when three training blocks are used the proposed scheme achieves a
classification accuracy of 33%, while the TRCA method achieve an accuracy of
22% (slighlty above the random guess).

Table 5: Classification accuracy (%) on Speller dataset with respect to the number of training trials - 9 channels Num | TRCA | adTRCA
---|---|---
3 | 61.5714 | 63.0238
4 | 82.9429 | 83.3000
5 | 82.9429 | 83.3000
6 | 86.0476 | 86.1429
Table 6: Classification accuracy (%) on EPOC dataset with respect to the number of training trials - 2 channels Num | TRCA | adTRCA
---|---|---
3 | 22.4242 | 33.3333
7 | 36.8831 | 39.4805
10 | 37.6364 | 39.8182
15 | 42.7879 | 42.7879
20 | 54.5455 | 55.2727

###  3.5 Ensemble case - Experiments

In the previous subsection we have presented experiments and we performed
comparisons using the basic version of our method. In the current subsection,
we provide the obtained results of the ensemble version of our method (see
2.5), and also, we perform a comparison with ensemble TRCA [18]. More
specifically, in Fig. 3 we provide the obtained results for the ensemble TRCA
(EnsembleTRCA) and ensemble Adaptive TRCA (Ensemble_adTRCA) methods. The
comparison between the two methods is performed with respect to the number of
channels and the datasets. For the Speller dataset we can observe that in the
case of 9 channels the two methods present similar performance. However, in
the case of 3 channels the Ensemble_adTRCA method provides significantly
better performance than the EnsembleTRCA method. Additionally, when we use the
EPOC dataset, the Ensemble_adTRCA method provides significantly better
performance, either using all 14 channels of the EPOC device or the two
channels covering the occipital lobe.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 3: Average Classification over all subjects by using for the Speller
dataset with (a) 9 channels and (b) 3 channels and for the EPOC dataset with
(c) 14 channels and (d) 2 channels, respectively. In both cases, the time
window ranges from 0.5s to 4s (0.5s interval). * indicates statistically
significant difference between the two methods using paired sample t-test for
Speller dataset and Wilcoxon signed rank test for the EPOC dataset ($p<0.05$).

##  4 Discussion and Conclusions

Enhancing the performance of SSVEP recognition is a significant issue for BCI
applications. In this study, we develop a multitask learning scheme to
strengthen the TRCA method. The idea behind the proposed learning scheme is to
develop a adaptive time - domain filter which can be used to a more general
eigenvalue problem than the corresponding problem of TRCA method. The proposed
method is able to deal with more general noises and with reduced number of
trials as the experiments have shown. However, this increase in performance
from our method has an increase also in the computation time of the overall
procedure since an iterative method is used to find the adaptive time domain
filters.

A significant part of our study is the use of two SSVEP datasets to evaluate
our method. Most of SSVEP studies use the Speller dataset to evaluate the
proposed methods. However, this dataset was created into a controlled
environment with high cost EEG equipment which make difficult to replicate the
study for real BCI applications with low cost equipment and in very noisy
environment. Hence the aforementioned methods that are evaluated on the
Speller dataset tend to underestimate the noise part of SSVEP EEG trials. An
effect which can be observe by comparing the performance of adTRCA and TRCA on
both datasets. We can see that adTRCA method provides much better performance
than the TRCA in the case of EPOC dataset. In the Speller dataset the adTRCA
method has around 1% better accuracy than the TRCA, while, in the EPOC dataset
this difference is increase to 5%. Additionally, the ensemble version of our
method presents better performance than the ensemble version of TRCA in both
datasets, especially when we have a limited number of channels.

The spatial filters and the SSVEP templates play important roles in the target
recognition methods. When the spatial filters and the SSVEP templates can not
be accurately computed, e.g. in the case of small calibration data or noisy
EEG recordings, the resulting recognition performance will be dramatically
decreased. Hence, to this challenge the key is how to estimate reliable
spatial filters. In this study, we present a novel spatial filtering approach
to recognize SSVEP signals. Our method use the multi-task idea to construct
adaptive time - domain filters resulting into a generalized eigenvalue problem
from where the final spatial filters are obtained. Extensive experiments,
using two SSVEP datasets, have been shown the usefulness of our method. The
proposed method significantly outperformed the TRCA, the CCA, the MLR and the
SRC methods in terms of classification accuracy and ITR. Finally, future
extensions of our approach could include transfer learning approaches
utilizing the data from all subjects to construct the recognition model.

## Acknowledgments

This work was part of project NeuroMkt that had been co-financed by the
European Regional Development Fund of the European Union and Greek National
Funds through the Operational Program Competitiveness, Entrepreneurship and
Innovation, under the call RESEARCH CREATE INNOVATE (Project code T2EDK-03661)

## Data Availability Statement

The datasets that have been used in this study are available in the Internet.
The Speller dataset can be found in
http://bci.med.tsinghua.edu.cn/download.html. The EPOC dataset can be found in
https://physionet.org/content/mssvepdb/1.0.0/.

## References

  * [1] J. R. Wolpaw, N. Birbaumer, D. J. McFarland, G. Pfurtscheller, T. M. Vaughan, Brain computer interfaces for communication and control, Clinical Neurophysiology 113 (6) (2002) 767 – 791.
  * [2] L. F. Nicolas-Alonso, J. Gomez-Gil, Brain computer interfaces, a review, Sensors 12 (2) (2012) 1211–1279.  doi:10.3390/s120201211.
URL https://www.mdpi.com/1424-8220/12/2/1211

  * [3] F. P. Kalaganis, K. Georgiadis, V. P. Oikonomou, N. A. Laskaris, S. Nikolopoulos, I. Kompatsiaris, Unlocking the subconscious consumer bias: A survey on the past, present, and future of hybrid eeg schemes in neuromarketing, Frontiers in Neuroergonomics 2 (2021).  doi:10.3389/fnrgo.2021.672982.
URL https://www.frontiersin.org/article/10.3389/fnrgo.2021.672982

  * [4] G. Bin, X. Gao, Y. Wang, B. Hong, S. Gao, VEP-based brain-computer interfaces: time, frequency, and code modulations , IEEE Computational Intelligence Magazine 4 (4) (2009) 22–26.
  * [5] R. Zerafa, T. Camilleri, O. Falzon, K. P. Camilleri, To train or not to train? a survey on training of feature extraction methods for SSVEP-based BCIs, Journal of Neural Engineering 15 (5) (2018) 051001.  doi:10.1088/1741-2552/aaca6e.
URL https://doi.org/10.1088/1741-2552/aaca6e

  * [6] S. Gao, Y. Wang, X. Gao, B. Hong, Visual and Auditory Brain Computer Interfaces, IEEE Transactions on Biomedical Engineering 61 (5) (2014) 1436–1447.
  * [7] V. P. Oikonomou, G. Liaros, K. Georgiadis, E. Chatzilari, K. Adam, S. Nikolopoulos, I. Kompatsiaris, Comparative evaluation of state-of-the-art algorithms for SSVEP-based BCIs, arXiv:1602.00904 (February 2016).
  * [8] H. Wang, Y. Zhang, N. R. Waytowich, D. J. Krusienski, G. Zhou, J. Jin, X. Wang, A. Cichocki, Discriminative feature extraction via multivariate linear regression for ssvep-based bci, IEEE Transactions on Neural Systems and Rehabilitation Engineering 24 (5) (2016) 532–541.  doi:10.1109/TNSRE.2016.2519350.
  * [9] V. P. Oikonomou, S. Nikolopoulos, I. Kompatsiaris, A bayesian multiple kernel learning algorithm for ssvep bci detection, IEEE Journal of Biomedical and Health Informatics 23 (5) (2019) 1990–2001.
  * [10] H. Cecotti, A time-frequency convolutional neural network for the offline classification of steady-state visual evoked potential responses, Pattern Recognition Letters 32 (8) (2011) 1145 – 1153.
  * [11] N.-S. Kwak, K. Muller, S.-W. S.-W.Lee, A convolutional neural network for steady state visual evoked potential classification under ambulatory environment, PLOS ONE 12 (2) (2017) 1–20.  doi:10.1371/journal.pone.0172578.
  * [12] H. Wang, T. Li, Z. Huang, Remote control of an electrical car with SSVEP-Based BCI, in: 2010 IEEE International Conference on Information Theory and Information Security, 2010, pp. 837–840.  doi:10.1109/ICITIS.2010.5689710.
  * [13] Z. Lin, C. Zhang, W. Wu, X. Gao, Frequency Recognition Based on Canonical Correlation Analysis for SSVEP-Based BCIs, IEEE Transactions on Biomedical Engineering 53 (12) (2006) 2610–2614.
  * [14] O. Friman, I. Volosyak, A. Graser, Multiple Channel Detection of Steady-State Visual Evoked Potentials for Brain-Computer Interfaces, IEEE Transactions on Biomedical Engineering 54 (4) (2007) 742–750.
  * [15] Y. Zhang, G. Zhou, J. Jin, M. Wang, X. Wang, A. Cichocki, L1-Regularized Multiway Canonical Correlation Analysis for SSVEP-Based BCI, IEEE Transactions on Neural Systems and Rehabilitation Engineering 21 (6) (2013) 887–896.
  * [16] G. Bin, X. Gao, Y. Wang, Y. Li, B. Hong, S. Gao, A high-speed BCI based on code modulation VEP, Journal of Neural Engineering 8 (2) (2011) 025015.
URL http://stacks.iop.org/1741-2552/8/i=2/a=025015

  * [17] M. Nakanishi, Y. Wang, Y. Wang, T. Jung, A Comparison Study of Canonical Correlation Analysis Based Methods for Detecting Steady-State Visual Evoked Potentials, PLoS ONE (2015) e0140703doi:10.1371/journal.pone.0140703.
  * [18] M. Nakanishi, Y. Wang, X. Chen, Y.-T. Wang, X. Gao, T.-P. Jung, Enhancing detection of ssveps for a high-speed brain speller using task-related component analysis, IEEE Transactions on Biomedical Engineering 65 (1) (2018) 104–112.  doi:10.1109/TBME.2017.2694818.
  * [19] C. Tong, H. Wang, C. Yang, X. Ni, Group ensemble learning enhances the accuracy and convenience of ssvep-based bcis via exploiting inter-subject information, Biomedical Signal Processing and Control 68 (2021) 102797.
  * [20] X. Yuan, Q. Sun, L. Zhang, H. Wang, Enhancing detection of ssvep-based bcis via a novel cca-based method, Biomedical Signal Processing and Control 74 (2022) 103482\.
  * [21] V. P. Oikonomou, S. Nikolopoulos, I. Kompatsiaris, Machine-learning techniques for eeg data, Signal Processing to Drive Human-Computer Interaction: EEG and eye-controlled interfaces (2020) 145.
  * [22] C. M. Wong, B. Wang, Z. Wang, K. Lao, A. Rosa, F. Wan, Spatial filtering in ssvep-based bcis: Unified framework and new improvements, IEEE Transactions on Biomedical Engineering 67 (11) (2020) 3057–3072.  doi:10.1109/TBME.2020.2975552.
  * [23] Y. Zhang, D. Guo, F. Li, E. Yin, Y. Zhang, P. Li, Q. Zhao, T. Tanaka, D. Yao, P. Xu, Correlated component analysis for enhancing the performance of ssvep-based brain-computer interface, IEEE Transactions on Neural Systems and Rehabilitation Engineering 26 (5) (2018) 948–956.
  * [24] B. Liu, X. Chen, N. Shi, Y. Wang, S. Gao, X. Gao, Improving the performance of individually calibrated ssvep-bci by task- discriminant component analysis, IEEE Transactions on Neural Systems and Rehabilitation Engineering 29 (2021) 1998–2007.  doi:10.1109/TNSRE.2021.3114340.
  * [25] C. M. Wong, F. Wan, B. Wang, Z. Wang, W. Nan, K. F. Lao, P. U. Mak, M. I. Vai, A. Rosa, Learning across multi-stimulus enhances target recognition methods in SSVEP-based BCIs, Journal of Neural Engineering 17 (1) (2020) 016026.  doi:10.1088/1741-2552/ab2373.
URL https://doi.org/10.1088/1741-2552/ab2373

  * [26] J. Jin, Z. Wang, R. Xu, C. Liu, X. Wang, A. Cichocki, Robust similarity measurement based on a novel time filter for ssveps detection, IEEE Transactions on Neural Networks and Learning Systems (2021) 1–10doi:10.1109/TNNLS.2021.3118468.
  * [27] L. Sun, S. Ji, J. Ye, Multi-Label Dimensionality Reduction, CRC Press, Taylor and Francis Group, 2014.
  * [28] S. Ji, D. Dunson, L. Carin, Multitask compressive sensing, IEEE Transactions on Signal Processing 57 (1) (2009) 92–106.  doi:10.1109/TSP.2008.2005866.
  * [29] M. E. Tipping, Sparse Bayesian Learning and the Relevance Vector Machine, Journal of Mach. Learn. Research 1 (2001) 211–244.
  * [30] D. J. MacKay, Bayesian interpolation, Neural Computation 4 (1992) 415–447.
  * [31] Y. Wang, X. Chen, X. Gao, S. Gao, A benchmark dataset for SSVEP-based brain - computer interfaces, IEEE Transactions on Neural Systems and Rehabilitation Engineering 25 (10) (2017) 1746–1752.  doi:10.1109/TNSRE.2016.2627556.
  * [32] V. P. Oikonomou, S. Nikolopoulos, I. Kompatsiaris, Sparse graph-based representations of ssvep responses under the variational bayesian framework, in: 2021 IEEE 21st International Conference on Bioinformatics and Bioengineering (BIBE), 2021, pp. 1–6.  doi:10.1109/BIBE52308.2021.9635427.


* [72] M. Schauer and R. Heinrich “Quasi-steady-state approximation in the mathematical modeling of biochemical reaction networks” In _Mathematical Biosciences_ 65.2, 1983, pp. 155–170 DOI: https://doi.org/10.1016/0025-5564(83)90058-5
  * [73] LA Segel “On the validity of the steady state assumption of enzyme kinetics” In _Bull Math Biology_ 50, 1988, pp. 579–593 DOI: https://doi.org/10.1007/BF02460092
  * [74] Lee A. Segel and Marshall Slemrod “The Quasi-Steady-State Assumption: A Case Study in Perturbation” In _SIAM Review_ 31.3, 1989, pp. 446–477
  * [75] Abraham Seidenberg “Some basic theorems in differential algebra (characteristic p arbitrary)” In _Transactions of the American Mathematical Society_ 73.1, 1952, pp. 174–190
  * [76] Yoav D Shaul and Rony Seger “The MEK/ERK cascade: from signaling specificity to diverse functions” In _Biochimica et Biophysica Acta (BBA)-Molecular Cell Research_ 1773.8 Elsevier, 2007, pp. 1213–1226
  * [77] Peng Shi, El-K Boukas and Ramesh K Agarwal “Kalman filtering for continuous-time uncertain systems with Markovian jumping parameters” In _IEEE Transactions on Automatic Control_ 44.8 IEEE, 1999, pp. 1592–1597
  * [78] Ralph C Smith “Uncertainty quantification: theory, implementation, and applications” Siam, 2013
  * [79] Eduardo D Sontag “For differential equations with r parameters, 2r+ 1 experiments are enough for identification” In _Journal of Nonlinear Science_ 12.6 Springer, 2002, pp. 553–583
  * [80] Mark A. Sweeney “Conditions for Solvability in Chemical Reaction Networks at Quasi-Steady-State”, 2017 arXiv:1712.05533 [math.DS]
  * [81] Koichi Takahashi, Sorin Tănase-Nicola and Pieter Rein Ten Wolde “Spatio-temporal correlations can drastically change the response of a MAPK pathway” In _Proceedings of the National Academy of Sciences_ 107.6 National Acad Sciences, 2010, pp. 2473–2478
  * [82] Clinton A. Taylor et al. “Functional divergence caused by mutations in an energetic hotspot in ERK2” In _Proceedings of the National Academy of Sciences_ 116.31 National Academy of Sciences, 2019, pp. 15514–15523 DOI: 10.1073/pnas.1905015116
  * [83] Clinton A. Taylor et al. “Functional divergence caused by mutations in an energetic hotspot in ERK2” In _Proceedings of the National Academy of Sciences_ 116.31 National Academy of Sciences, 2019, pp. 15514–15523 DOI: 10.1073/pnas.1905015116
  * [84] Matthew Thomson and Jeremy Gunawardena “The rational parameterisation theorem for multisite post-translational modification systems” In _Journal of theoretical biology_ 261.4 Elsevier, 2009, pp. 626–636
  * [85] Matthew Thomson and Jeremy Gunawardena “Unlimited multistability in multisite phosphorylation systems” In _Nature_ 460.7252 Nature Publishing Group, 2009, pp. 274–277
  * [86] Sandor Vajda, Herschel Rabitz, Eric Walter and Yves Lecourtier “Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic models” In _Chemical Engineering Communications_ 83.1, 1989, pp. 191–219 URL: http://dx.doi.org/10.1080/00986448908940662
  * [87] Alejandro F. Villaverde, Neil D. Evans, Michael J. Chappell and Julio R. Banga “Input-Dependent Structural Identifiability of Nonlinear Systems” In _IEEE Control Systems Letters_ 3.2, 2019, pp. 272–277 DOI: 10.1109/LCSYS.2018.2868608
  * [88] Alejandro F. Villaverde, Neil D. Evans, Michael J. Chappell and Julio R. Banga “Sufficiently Exciting Inputs for Structurally Identifiable Systems Biology Models’Åé’ÅéThis project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 686282 (“CANPATHPRO”).” 7th Conference on Foundation of Systems Biology in Engineering FOSBE 2018 In _IFAC-PapersOnLine_ 51.19, 2018, pp. 16–19 DOI: https://doi.org/10.1016/j.ifacol.2018.09.015
  * [89] Sean T. Vittadello and Michael P. H. Stumpf “Model comparison via simplicial complexes and persistent homology”, 2020 arXiv:2012.13039 [math.AT]
  * [90] Sean T. Vittadello and Michael P. H. Stumpf “Model comparison via simplicial complexes and persistent homology” In _Royal Society Open Science_ 8.10 The Royal Society, 2021 DOI: 10.1098/rsos.211361
  * [91] Eberhard Voit “A first course in systems biology” Garland Science, 2017
  * [92] Liming Wang and Eduardo D Sontag “On the number of steady states in a multiple futile cycle” In _Journal of mathematical biology_ 57.1 Springer, 2008, pp. 29–52
  * [93] Shu Wang, Jia-Ren Lin, Eduardo D Sontag and Peter K Sorger “Inferring reaction network structure from single-cell, multiplex data, using toric systems theory” In _PLoS computational biology_ 15.12 Public Library of Science San Francisco, CA USA, 2019, pp. e1007311
  * [94] E Yeung et al. “Inference of Multisite Phosphorylation Rate Constants and Their Modulation by Pathogenic Mutations” In _Current Biology_ https://doi.org/10.1016/j.cub.2019.12.052, 2019

## Appendix

###  A.1. QSSA Model Reduction via the Classical Singular Perturbation Theory
Approach

Recall the Full ERK model (Equations (1a) to (1f)):

$\displaystyle\frac{\mathrm{d}S_{0}}{\mathrm{d}t}$
$\displaystyle=-k_{f_{1}}E\cdot S_{0}+k_{r_{1}}C_{1},$
$\displaystyle\frac{\mathrm{d}C_{1}}{\mathrm{d}t}$
$\displaystyle=k_{f_{1}}E\cdot S_{0}-(k_{r_{1}}+k_{c_{1}})C_{1},$
$\displaystyle\frac{\mathrm{d}C_{2}}{\mathrm{d}t}$
$\displaystyle=k_{c_{1}}C_{1}-(k_{r_{2}}+k_{c_{2}})C_{2}+k_{f_{2}}E\cdot
S_{1},$ $\displaystyle\frac{\mathrm{d}S_{1}}{\mathrm{d}t}$
$\displaystyle=-k_{f_{2}}E\cdot S_{1}+k_{r_{2}}C_{2},$
$\displaystyle\frac{\mathrm{d}S_{2}}{\mathrm{d}t}$
$\displaystyle=k_{c_{2}}C_{2},$ $\displaystyle\frac{\mathrm{d}E}{\mathrm{d}t}$
$\displaystyle=-k_{f_{1}}E\cdot S_{0}+k_{r_{1}}C_{1}-k_{f_{2}}E\cdot
S_{1}+(k_{r_{2}}+k_{c_{2}})C_{2},$

with initial conditions $S_{0}=S_{tot}$, $E=E_{tot}$, and
$S_{1}=S_{2}=C_{1}=C_{2}=0$ at time $t=0$. Note that the ODE for $S_{2}$
decouples from the system. It is straightforward to show that the system
satisfies the following conservation laws:

$\displaystyle S_{tot}$ $\displaystyle=S_{0}+S_{1}+S_{2}+C_{1}+C_{2},$
$\displaystyle E_{tot}$ $\displaystyle=E+C_{1}+C_{2}.$

We then non-dimensionalise by substituting $\tilde{S}_{i}:=S_{i}/S_{tot}$ for
$i=0,1,2$; $\tilde{C}_{j}:=C_{j}/E_{tot}$ for $j=1,2$; and
$\tilde{E}=1-\tilde{C}_{1}-\tilde{C}_{2}$ which is analogous to Equation (7)
and removes one equation from the ODE system. We also observe that $S_{2}$
decouples, leaving us with a system of four differential equations. Then

$\displaystyle\frac{\mathrm{d}\tilde{S}_{0}}{\mathrm{d}t}$
$\displaystyle=-k_{f_{1}}E_{tot}\tilde{E}\cdot\tilde{S}_{0}+k_{r_{1}}\frac{E_{tot}}{S_{tot}}\tilde{C}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{C}_{1}}{\mathrm{d}t}$
$\displaystyle=k_{f_{1}}S_{tot}\tilde{E}\cdot\tilde{S}_{0}-(k_{r_{1}}+k_{c_{1}})\tilde{C}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{C}_{2}}{\mathrm{d}t}$
$\displaystyle=k_{c_{1}}\tilde{C}_{1}-(k_{r_{2}}+k_{c_{2}})\tilde{C}_{2}+k_{f_{2}}S_{tot}\tilde{E}\cdot\tilde{S}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{S}_{1}}{\mathrm{d}t}$
$\displaystyle=-k_{f_{2}}E_{tot}\tilde{E}\cdot\tilde{S}_{1}+k_{r_{2}}\frac{E_{tot}}{S_{tot}}\tilde{C}_{2}.$

Subsequently, substituting $T:=(E_{tot}k_{f_{1}})\cdot t$ for $t$:

$\displaystyle\frac{\mathrm{d}\tilde{S}_{0}}{\mathrm{d}T}$
$\displaystyle=-\tilde{E}\cdot\tilde{S}_{0}+k_{r_{1}}\frac{1}{S_{tot}k_{f_{1}}}\tilde{C}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{C}_{1}}{\mathrm{d}T}$
$\displaystyle=\frac{S_{tot}}{E_{tot}}\tilde{E}\cdot\tilde{S}_{0}-\frac{1}{E_{tot}k_{f_{1}}}(k_{r_{1}}+k_{c_{1}})\tilde{C}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{C}_{2}}{\mathrm{d}T}$
$\displaystyle=\frac{k_{c_{1}}}{E_{tot}k_{f_{1}}}\tilde{C}_{1}-\frac{k_{r_{2}}+k_{c_{2}}}{E_{tot}k_{f_{1}}}\tilde{C}_{2}+\frac{k_{f_{2}}}{k_{f_{1}}}\frac{S_{tot}}{E_{tot}}\tilde{E}\cdot\tilde{S}_{1},$
$\displaystyle\frac{\mathrm{d}\tilde{S}_{1}}{\mathrm{d}T}$
$\displaystyle=-\frac{k_{f_{2}}}{k_{f_{1}}}\frac{S_{tot}}{E_{tot}}\tilde{E}\cdot\tilde{S}_{1}+\frac{k_{r_{2}}}{S_{tot}k_{f_{1}}}\tilde{C}_{2}.$

Then

$\displaystyle\varepsilon\frac{\mathrm{d}\tilde{C}_{1}}{\mathrm{d}T}$
$\displaystyle=\tilde{E}\cdot\tilde{S}_{0}-\frac{1}{S_{tot}k_{f_{1}}}(k_{r_{1}}+k_{c_{1}})\tilde{C}_{1},$
$\displaystyle\varepsilon\frac{\mathrm{d}\tilde{C}_{2}}{\mathrm{d}T}$
$\displaystyle=\frac{k_{c_{1}}}{S_{tot}k_{f_{1}}}\tilde{C}_{1}-\frac{k_{r_{2}}+k_{c_{2}}}{S_{tot}k_{f_{1}}}\tilde{C}_{2}+\frac{k_{f_{2}}}{k_{f_{1}}}\tilde{E}\cdot\tilde{S}_{1},$

where $\varepsilon:=E_{tot}/S_{tot}$. As $E_{tot}\ll S_{tot}$, we have
$\varepsilon\approx 0$ and assume $\varepsilon\to 0$. This gives

$\displaystyle\tilde{C}_{1}$
$\displaystyle=\frac{k_{f_{1}}S_{tot}K_{2}\tilde{S}_{0}}{((k_{f_{1}}\tilde{S}_{0}+k_{f_{2}}\tilde{S}_{1})k_{c_{1}}+\tilde{S}_{0}k_{f_{1}}K_{2}+\tilde{S}_{1}k_{r_{1}}k_{f_{2}})S_{tot}+K_{1}K_{2}},$
$\displaystyle\tilde{C}_{2}$
$\displaystyle=\frac{((k_{f_{1}}\tilde{S}_{0}+k_{f_{2}}\tilde{S}_{1})k_{c_{1}}+k_{r_{1}}k_{f_{2}}\tilde{S}_{1})S_{tot}}{(S_{tot}(k_{f_{1}}\tilde{S}_{0}+k_{f_{2}}\tilde{S}_{1})+k_{c_{2}}+k_{r_{2}})k_{c_{1}}+(k_{r_{1}}k_{f_{2}}\tilde{S}_{1}+k_{f_{1}}K_{2}\tilde{S}_{0})S_{tot}+k_{r_{1}}K_{2}},$

where $K_{i}:=k_{r_{i}}+k_{c_{i}}$.

Finally, by substituting the above expressions for $\tilde{C}_{1}$ and
$\tilde{C}_{2}$:

$\displaystyle\frac{\mathrm{d}\tilde{S}_{0}}{\mathrm{d}T}$
$\displaystyle=\frac{-k_{c_{1}}K_{2}\tilde{S}_{0}}{S_{tot}(k_{c_{1}}+k_{c_{2})}(k_{f_{1}}\tilde{S}_{0}+k_{f_{2}}\tilde{S}_{1})+K_{1}K_{2}},$
$\displaystyle\frac{\mathrm{d}\tilde{S}_{1}}{\mathrm{d}T}$
$\displaystyle=\frac{-(k_{f_{2}}/k_{f_{1}})k_{c_{2}}K_{1}\tilde{S}_{1}+k_{r_{2}}k_{c_{1}}\tilde{S}_{0}}{S_{tot}(k_{c_{1}}+k_{c_{2})}(k_{f_{1}}\tilde{S}_{0}+k_{f_{2}}\tilde{S}_{1})+K_{1}K_{2}}.$

Then, by substituting $S_{i}$ for $\tilde{S}_{i}$, and dividing both numerator
and denominator by $K_{1}K_{2}$

$\displaystyle\frac{\mathrm{d}S_{0}}{\mathrm{d}t}$
$\displaystyle=\frac{-\kappa_{1}S_{0}}{\gamma_{1}S_{0}+\gamma_{2}S_{1}+1},$
$\displaystyle\frac{\mathrm{d}S_{1}}{\mathrm{d}t}$
$\displaystyle=\frac{-\kappa_{2}S_{1}+(1-\pi)\kappa_{1}S_{0}}{\gamma_{1}S_{0}+\gamma_{2}S_{1}+1},$

where $\kappa_{i}$, $\pi$, and
$\gamma_{i}:=(k_{c_{1}}+k_{c_{2}})k_{f_{i}}/(K_{1}K_{2})$ are as defined
previously in this paper.

This system reduces to the linear case if the denominator is approximately
constant. In particular, this is the case if
$0<\gamma_{1}S_{0}+\gamma_{2}S_{1}\ll 1$. This is the case, for our setting,
which can be seen by writing

$\gamma_{1}=\frac{k_{c_{1}}+k_{c_{2}}}{k_{c_{2}}+k_{r_{1}}}\frac{1}{k_{M_{1}}}$

(similarly for $\theta_{2}$), as we expect the first factor to be
approximately 1 [8] or smaller and $k_{M_{1}}$ (similarly $k_{M_{2}}$), the
Michaelis-Menten constant, to be approximately $25\mu M$.

###  A.2. QSSA Model Reduction via the Algebraic Approach

We use the algebraic approach to derive two reduced models.

####  A.2.1. Deriving the Rational ERK Model

We substitute $E$ using Equation (7) into the Full ERK model (Equations (1)).
Using the algebraic reduction, we follow the steps of Appendix A.1, except
instead of taking a limit $\varepsilon\to 0$, simply set
$\varepsilon(\mathrm{d}\tilde{C}_{1}/\mathrm{d}T)=\varepsilon(\mathrm{d}\tilde{C}_{2}/\mathrm{d}T)=0$,
which leads to the Rational ERK model. In other words, we treat the assumption
that $\varepsilon(\mathrm{d}\tilde{C}_{1}/\mathrm{d}T)$ and
$\varepsilon(\mathrm{d}\tilde{C}_{2}/\mathrm{d}T)$ converge to zero as an
algebraic fact and without regard of any analytic implications the assumption
$\varepsilon\to 0$ could have on the right-hand-side of these differential
equations. Note that we do not need to non-dimensionalise as we did in
Appendix A.1. We justify the accuracy of such an approximation Appendix A.3.

####  A.2.2. Deriving the Linear ERK Model

We now substitute $E$ using Equation (8),

$E=E_{tot}-S_{tot}+S_{0}+S_{1}+S_{2}=:E^{\prime},$

into the Full ERK model (Equations (1)). Setting
$\mathrm{d}C_{1}/\mathrm{d}t=0$ and $\mathrm{d}C_{2}/\mathrm{d}t=0$, we get

$\displaystyle C_{1}$
$\displaystyle=\frac{k_{f_{1}}}{k_{r_{1}}+k_{c_{1}}}E^{\prime}\cdot S_{0},$
$\displaystyle C_{2}$
$\displaystyle=\frac{1}{k_{r_{2}}+k_{c_{2}}}\left(k_{f_{2}}E^{\prime}\cdot
S_{1}+k_{c_{1}}C_{1}\right)=\frac{E^{\prime}}{k_{r_{2}}+k_{c_{2}}}\left(k_{f_{2}}S_{1}+\frac{k_{f_{1}}k_{c_{1}}}{k_{r_{1}}+k_{c_{1}}}S_{0}\right),$

which can be substituted into $\mathrm{d}S_{0}/\mathrm{d}t$,
$\mathrm{d}S_{1}/\mathrm{d}t$ and $\mathrm{d}S_{2}/\mathrm{d}t$:

$\displaystyle\frac{\mathrm{d}S_{0}}{\mathrm{d}t}$
$\displaystyle=-\frac{k_{f_{1}}k_{c_{1}}}{k_{r_{1}}+k_{c_{1}}}E^{\prime}\cdot
S_{0},$ $\displaystyle\frac{\mathrm{d}S_{1}}{\mathrm{d}t}$
$\displaystyle=-\frac{k_{f_{2}}k_{c_{2}}}{k_{r_{2}}+k_{c_{2}}}E^{\prime}\cdot
S_{1}+\frac{k_{r_{2}}}{k_{r_{2}}+k_{c_{2}}}\frac{k_{f_{1}}k_{c_{1}}}{k_{r_{1}}+k_{c_{1}}}E^{\prime}\cdot
S_{0},$ $\displaystyle\frac{\mathrm{d}S_{2}}{\mathrm{d}t}$
$\displaystyle=\frac{k_{f_{2}}k_{c_{2}}}{k_{r_{2}}+k_{c_{2}}}E^{\prime}\cdot
S_{1}+\frac{k_{c_{2}}}{k_{r_{2}}+k_{c_{2}}}\frac{k_{f_{1}}k_{c_{1}}}{k_{r_{1}}+k_{c_{1}}}E^{\prime}\cdot
S_{0}.$

We rewrite the above equations in terms of the parameter combinations
$\kappa_{1}$, $\kappa_{2}$, and $\pi$ (see Equation (11)) then gives the
Linear ERK model (Equations (10a)–(10c)).

Observe that the Linear ERK model obeys the conservation law
$S_{0}+S_{1}+S_{2}=\tilde{S}_{tot}$. Here, $\tilde{S}_{tot}$ is a constant
close to $S_{tot}$. Hence, $E^{\prime}=E_{tot}-(S_{tot}-\tilde{S}_{tot})$ is a
constant close to $E_{tot}$. The implication that $E^{\prime}=E_{tot}$ at all
times $t$ is physically infeasible. Whether $\tilde{S}_{tot}=S_{tot}$ depends
on whether one updates the initial conditions for $S_{j}$ through an inner
solution, which is typically assumed in singular-perturbation-theory
approaches. The algebraic approach does not require computing an inner
solution.

It is important to reiterate at this point that this reduction is not the
result of a singular perturbation analysis. Indeed, if we had followed a
singular perturbation analysis, similar to the one in Appendix A.1, but with
substitution (8), we could not factor out $E_{tot}/S_{tot}$ in the non-
dimensionalised differential equations for $C_{1}$ and $C_{2}$ as nicely as we
were able to in the previous subsection. In this instance, we would have to
leave factors of $\varepsilon^{-1}$ in the algebraic equations, which would be
ambiguous when taking a limit $\varepsilon\to 0$.

###  A.3. Accuracy of Algebraic QSSA

We start by providing the full statement of Proposition 2 of [29] (i.e., the
full version of Proposition 4 of this manuscript):

Let
$K^{*}\subset\operatorname{\mathbb{R}}_{+}^{n}\times\operatorname{\mathbb{R}}_{+}^{m}$
satisfy the following:

  * •

There exists $(\hat{y},\hat{\theta})$ in the interior of $K^{*}$ such that
$f^{[2]}(\hat{y},\hat{\theta})=0$.

  * •

$D_{2}f^{[2]}(x,\theta)$ is invertible for all $(x,\theta)\in K^{*}$.

  * •

There exist $y_{0}\in\operatorname{\mathbb{R}}^{n}$ and $r>0$ with the
following property:
Whenever $(x,\theta)\in K^{*}$ for some $x\in\operatorname{\mathbb{R}}^{n}$
and some $\theta\in\operatorname{\mathbb{R}}^{m}_{+}$ then
$\overline{B_{r}(y_{0})}\times\\{\theta\\}\subseteq K^{*}$.

  * •

There exists an $R>0$ such that $\|f(x,\theta)\|\leq R$ and
$\|f_{\mathrm{red}}(x,\theta)\|\leq R$ for all $(x,\theta)\in K^{*}$.

  * •

There exists an $L>0$ such that $\|Df(x,\theta)\|\leq L$ and
$\|Df_{\mathrm{red}}(x,\theta)\|\leq L$ for all $(x,\theta)\in K^{*}$.

These conditions imply that every $V_{\theta^{*}}$, with ${\theta^{*}}$ near
$\hat{\theta}$, is a submanifold. Note that every $(y,\theta^{*})$ with $y$ in
the interior of $\operatorname{\mathbb{R}}^{n}_{+}$ is contained in some
$K^{*}$ that satisfies the last three of the above conditions.

######  Proposition 18.

Assume that the above conditions are satisfied for $K^{*}$.

(a) Let $\theta^{*}$ be given such that $V_{\theta^{*}}\times\\{\theta^{*}\\}$
has non-empty intersection with $\mathrm{int}\,K^{*}$, let $(y,\theta^{*})$ be
a point in this intersection and
$V^{\prime}_{\theta^{*}}\subset\operatorname{\mathbb{R}}^{n}$ be some open
neighborhood of $y$ such that $(V_{\theta^{*}}\cap
V^{\prime}_{\theta^{*}})\times\\{\theta^{*}\\}\subset K^{*}$. Moreover let
$t^{*}>0$ such that the solution of (4) with initial value $y$ exists and
remains in $V^{\prime}_{\theta^{*}}$ for all $t\in[0,t^{*}]$. Then there
exists a compact neighbourhood $A_{\theta^{*}}\subseteq
V^{\prime}_{\theta^{*}}$ of $y$ with the following properties: (i) For every
$z\in A_{\theta^{*}}$ the solution of (4) with initial value $z$ exists and
remains in $V^{\prime}_{\theta^{*}}$ for all $t\in[0,t^{*}]$. (ii) For every
$\varepsilon^{\prime}>0$ there is a $\delta_{1}>0$ such that the solution of
(6) with initial value $z\in A_{\theta^{*}}\cap V_{\theta^{*}}$ exists and
remains in $V_{\theta^{*}}$ for $t\in[0,t^{*}]$ whenever
$\|f-f_{\mathrm{red}}\|<\delta_{1}$ on $V^{\prime}_{\theta^{*}}$. (iii) For
every $\varepsilon^{\prime}>0$ there is a $\delta\in(0,\delta_{1}]$ such that
the difference of the solutions of (4) resp. of (6) with initial value $z\in
A_{\theta^{*}}\cap V_{\theta^{*}}$ has norm less than $\varepsilon^{\prime}$
for all $t\in[0,t^{*}]$ whenever $\|f-f_{\mathrm{red}}\|<\delta$ on
$V^{\prime}_{\theta^{*}}$.

(b) Let $y\in V_{\theta^{*}}$ and let $\rho_{0}>0$ such that

$\overline{B_{\rho_{0}/2L}(y)}\times\\{{\theta^{*}}\\}\subseteq K^{*}.$

Let $\rho\leq\rho_{0}$ such that
$\|f(y,{\theta^{*}})-f_{\mathrm{red}}(y,{\theta^{*}})\|\geq 2\rho$. Then for
$t^{\prime}:=\rho/(2LR)$ the solutions of (4) resp. of (6) with initial value
$y$ exist and remain in $\overline{B_{\rho_{0}/2L}(y)}$ for $0\leq t\leq
t^{\prime}$, and their difference has norm at least $\rho^{2}/(2LR)$ at
$t=t^{\prime}$.

Returning to our two model reductions of interest, we use the notation
Sections 2 and 3 and label the two possible substitutions for the variable $E$
as (7) and (8), as before. Moreover, we write $\dot{S_{0}}$ and $\dot{S_{1}}$
as shorthand for the algebraic expressions for the time derivatives of $S_{0}$
and $S_{1}$ respectively. Recall that $K_{i}:=k_{c_{i}}+k_{r_{i}}$.
Throughout, we will view our set of polynomials $f$, which govern the ODE
system, as a smooth function from $\operatorname{\mathbb{R}}^{5}$ to
$\operatorname{\mathbb{R}}$. Any norm in this subsection will refer to the
$\|\,\cdot\,\|_{\infty}$-norm on the respective space.

Define a domain $K^{*}$ such that $k_{f_{i}}\ll K_{i}$ for $i=1,2$ and such
that it satisfies all hypotheses in the list at the start of this section
(i.e., from Proposition 2 of [29]).

Note that for both substitutions, $E$ is bounded above by $E_{tot}\approx
0.65\mu E$. Hence, the norms on the slow variables are going to be
approximately the same for both substitutions.

For the case (7), we get that $\|f_{\mathrm{red}}\|$ (showing only the
components that have changed compared to $f$) is given by

$\begin{Vmatrix}k_{f_{1}}E(S_{0}-\det^{-1}\dot{S_{0}}(k_{f_{2}}S_{1}+K_{2})+\dot{S_{1}}k_{f_{2}}S_{0})-K_{1}C_{1}\\\
k_{f_{2}}E(S_{1}-\det^{-1}(k_{f_{1}}(\dot{S_{1}}S_{0}+\dot{S_{0}}S_{1}))+\det^{-1}E(\dot{S_{0}}k_{f_{1}}k_{c_{2}}-\dot{S_{1}}k_{f_{2}})+k_{c_{1}}C_{1}-K_{2}C_{2}\end{Vmatrix}$

on variables $C_{1}$ and $C_{2}$. Here, the quantity $\det$ is given by

$\det=k_{f_{1}}(k_{c_{1}}+K_{2})S_{0}+K_{1}(k_{f_{2}}S_{1}+K_{2})=\mathcal{O}(K_{1}K_{2}).$

Similarly, for (8), we get

$\begin{Vmatrix}k_{f_{1}}\left(ES_{0}-\frac{E+S_{0}}{K_{1}}\dot{S_{0}}\right)-K_{1}C_{1}+\frac{k_{c_{1}}k_{f_{2}}}{K_{1}K_{2}}(E+S_{1})\dot{S_{1}}\\\
k_{f_{2}}\left(ES_{1}-\frac{E+S_{1}}{K_{2}}\dot{S_{1}}\right)+k_{c_{1}}C_{1}-K_{2}C_{2}\end{Vmatrix}.$

In both cases, applying the triangle inequality shows that the upper entry is
of order $|K_{1}C_{1}|$ while the lower entry is of order
$|k_{c_{1}}C_{1}-K_{2}C_{2}|$. It follows that the quantity $R$, as in the
assumptions of Proposition 18, is of a similar order for both substitutions.

For $\|Df_{\mathrm{red}}\|$, we get

$\begin{Vmatrix}k_{f_{1}}E\left(1+\frac{k_{f_{1}}(\det+\dot{S_{0}}(k_{c_{1}}+K_{2})}{\det^{2}}(k_{f_{2}}S_{1}+K_{2})+k_{f_{2}}\dot{S_{1}}\right)\\\
-k_{f_{2}}^{2}S_{0}\\\
-k_{f_{1}}\left(S_{0}-\frac{\dot{S_{0}}(k_{f_{2}}S_{1}+K_{2})}{\det}+\dot{S_{1}}k_{f_{2}}S_{0}-E\frac{k_{r_{1}}(k_{f_{2}}S_{1}+K_{2})}{\det}\right)-K_{1}\\\
0\end{Vmatrix}^{T}$

for the row of the Jacobian corresponding to the entry of $f_{\mathrm{red}}$
giving the rate of change of $C_{1}$ and

$\begin{Vmatrix}Ek_{f_{1}}\frac{k_{f_{2}}((k_{f_{1}}S_{1}-\dot{S_{1}})+(\dot{S_{1}}S_{0}+\dot{S_{0}}S_{1})k_{f_{1}}(k_{c_{1}}+K_{2}))+k_{f_{1}}k_{c_{2}}\det+(\hat{S_{1}}k_{f_{2}}-k_{f_{1}}k_{c_{2}}\hat{S_{0}})(k_{c_{1}}+K_{2})}{\det^{2}}\\\
Ek_{f_{2}}\frac{\det^{2}+k_{f_{1}}k_{f_{2}}S_{0}\det+k_{f_{1}}k_{f_{2}}(\dot{S_{1}}S_{0}+\dot{S_{0}}S_{1})K_{1}+k_{f_{2}}\det-(\dot{S_{0}}k_{f_{1}}k_{c_{2}}-\dot{S_{1}}k_{f_{2}})K_{1}}{\det^{2}}\\\
-\frac{k_{f_{2}}k_{r_{1}}ES_{1}}{\det}+\frac{k_{f_{1}}k_{c_{2}}k_{r_{1}}E}{\det}+k_{c_{1}}\\\
-\frac{k_{f_{2}}k_{r_{2}}ES_{0}}{\det}+E\frac{k_{f_{2}}k_{r_{2}}}{\det}+K_{2}\end{Vmatrix}^{T}$

for the row of the Jacobian corresponding to the entry of $f_{\mathrm{red}}$
giving the rate of change of $C_{2}$ for the substitution (7). Here, the four
columns (or rows in the given transposed presentation) correspond to partial
derivatives with respect to $S_{0}$, $S_{1}$, $C_{1}$ and $C_{2}$ (in that
order). Note that the rows of $Df_{\mathrm{red}}$ corresponding to the rates
of change of $S_{0}$, $S_{1}$, and $S_{2}$ are identical for both
substitutions. Thus, we omit their calculations.

Similarly, for (8)

$\begin{Vmatrix}k_{f_{1}}\left(S_{0}+E-\frac{2\dot{S_{0}}-k_{f_{1}}(E+S_{0})}{K_{1}}\right)+\frac{k_{c_{1}}k_{f_{2}}}{K_{1}K_{2}}\dot{S_{1}}&k_{f_{1}}\left(S_{1}-\frac{\dot{S_{1}}}{K_{2}}\right)\\\
k_{f_{1}}\left(S_{0}-\frac{\dot{S_{0}}}{K_{1}}\right)-2k_{f_{2}}\frac{k_{c_{1}}k_{f_{2}}}{K_{1}K_{2}}&k_{f_{1}}\left(E+S_{1}-\frac{2\dot{S_{1}}-k_{f_{1}}(E+S_{1})}{K_{2}}\right)\\\
-k_{f_{1}}\frac{\dot{S_{0}}+k_{r_{1}}(E+S_{0})}{K_{1}}-K_{1}&-k_{f_{2}}\frac{\dot{S_{1}}+k_{r_{2}}(E+S_{1})}{K_{2}}+k_{c_{1}}\\\
0&K_{2}\end{Vmatrix}^{T}.$

We note that the dominant terms are exactly $K_{1}$, $k_{c_{1}}$, and $K_{2}$
for both substitutions. We conclude that the quantity $L$, as in the
assumptions of Proposition 18, is of similar a order for both substitutions.

For the quantity $\|f-f_{\mathrm{red}}\|$, for substitution (7), we get

$\begin{Vmatrix}k_{f_{1}}E(\dot{S_{1}}k_{f_{2}}S_{0}-\det^{-1}\dot{S_{0}}(k_{f_{2}}S_{1}+K_{2}))\\\
\det^{-1}k_{f_{1}}k_{c_{2}}E\dot{S_{0}}-\det^{-1}k_{f_{1}}k_{f_{2}}E(\dot{S_{0}}S_{1}+\dot{S_{1}}S_{0})-\det^{-1}k_{f_{2}}K_{1}E\dot{S_{1}}\end{Vmatrix}.$

Similarly, for (8) we get

$\begin{Vmatrix}\frac{k_{c_{1}}k_{f_{1}}}{K_{1}K_{2}}(E+S_{1})\dot{S_{1}}-k_{f_{1}}\frac{E+S_{0}}{K_{1}}\dot{S_{0}}\\\
-k_{f_{2}}\frac{E+S_{1}}{K_{2}}\dot{S_{1}}\end{Vmatrix}.$

Both of the above norms are $\mathcal{O}(1)$ (as $\dot{S_{0}}$ and
$\dot{S_{1}}$ contain terms $K_{1}$ and $K_{2}$ respectively).

Hence, all quantities use to bound accuracy in Proposition 18 (above, by part
(a), or below, by part (b)) are of a similar order. We conclude that the
accuracy of both reductions is approximately the same when assuming $0\leq
k_{f_{i}}\ll K_{i}$. This follows from the measurements $k_{M_{i}}\approx
25\,\mu M$ [8]. Moreover, we have seen previously (in the main text) that
$K_{i}>0$ for $i=1,2$ is a sufficient condition for a parameter value $\theta$
to be a QSS-parameter value in the sense of [29].

###  A.4. Deriving Equality of ideals

In this section of the appendix, we derive that
$I(\mathcal{V}_{0})=I_{\Sigma}$ for the ERK models with given initial
condition. In Saccomani et al. [71], $\mathcal{V}_{0}$ is defined as all
trajectories of a given ODE system with initial condition $x_{0}$. That is, if
$x(t)$ satisfies $S(0)=x_{0}$ and $\dot{S}=f(S,\theta)$, then
$\mathcal{V}_{0}=\\{x(t)\,|\,t\geq
0\\}\subset\operatorname{\mathbb{C}}(\theta)^{n}$. Then $I(\mathcal{V}_{0})$
is defined to be the ideal of all polynomials in
$\operatorname{\mathbb{C}}(\theta)[S]$ vanishing on $\mathcal{V}_{0}$.

######  Proposition 19.

Assume we are given a complex-valued ODE model, including an initial
condition. Let $\theta$ denote the set of its parameters and $S$ the set of
its variables. Consider the affine space $\mathcal{A}$ given by its variables
and their derivative terms (of all orders), viewed as an affine space over the
fraction field $\mathbb{C}(\theta)$. Define $\mathcal{V}_{0}$ to be the
manifold in $\mathcal{A}$ given by the ODE trajectories under the given
initial condition.

Then for the Full, Rational and Linear ODE model considered in this
manuscript, we have that $I(\mathcal{V}_{0})$, the differential ideal of all
polynomials in the differential ring $\mathbb{C}(\theta)[S]$ vanishing on
$\mathcal{V}_{0}$ (c.f. [71]), is exactly $I_{\Sigma}$.

###### Proof.

By definition, $I(\mathcal{V}_{0})\supset I_{\Sigma}$. Hence, we need only
show the inclusion $I_{\Sigma}\subset I(\mathcal{V}_{0})$. by way of
contradiction, suppose there are elements of $I(\mathcal{V}_{0})$ which are
not elements of $I_{\Sigma}$. Without loss of generality, we may assume that
such polynomials contain no derivatives of variables. Indeed, starting with
any polynomial in $I(\mathcal{V}_{0})$, we can use the differential equations
in $I_{\Sigma}$ to cancel out all terms containing derivatives of variables
and so obtain an element of $I(\mathcal{V}_{0})$ without such terms.

We need to take some extra care for the Rational ERK model: strictly speaking,
the polynomial $A(\gamma_{1}S_{0}+\gamma_{2}S_{1}+1)-1$ is in
$I(\mathcal{V}_{0})$ but not in $I_{\Sigma}$. However, we know this relation a
priori and it is merely a result of us converting a rational ODE model into a
polynomial one. In particular, it is independent of an initial condition.
Hence, in the context of this section, we will assume that the Rational ERK
model contains an output variable
$y_{3}:=A(\gamma_{1}S_{0}+\gamma_{2}S_{1}+1)$ (in addition to $y_{1}=S_{0}$
and $y_{2}=S_{1}$). Whilst $y_{3}$ adds polynomials to $I_{\Sigma}$, we will
omit it from any polynomials we study henceforth as it is equivalent to the
constant 1 (given the definition of $A$).

For the Rational ERK model, assume there exists a non-zero
$p\in\mathbb{C}(\theta)[S_{0},S_{1},A]$ such that $p(S_{0},S_{1},A)=0$ for all
$t\geq 0$ given our initial condition. By using
$A(\gamma_{1}S_{0}+\gamma_{2}S_{1}+1)=1$, we can turn $p$ into a rational
function in variables $S_{0},S_{1}$ which vanishes on the ODE trajectories.
Hence, there must exist a polynomial $q\in\mathbb{C}(\theta)[S_{0},S_{1}]$
such that $q(S_{0},S_{1})=0$ for all $t\geq 0$ given our initial condition.
Any properties we will use in the following argument are satisfied by both the
Rational and Linear ERK models unless stated otherwise. As we need to disprove
the existence of such polynomial $q$ when studying the Linear ERK model the
following conclusion will hold for both the Linear and Rational ERK models.

Without loss of generality, assume that $q$ is a polynomial of the smallest
non-zero degree satisfying $q(S_{0},S_{1})=0$ on the ODE trajectories and that
$q$ is irreducible. If it were reducible, simply replace $q$ with an
irreducible factor of $q$ such that $q(S_{0},S_{1})=0$ for all $0\leq t<t_{1}$
for some $t_{1}>0$. The following argument will still hold.

The Linear and Rational ERK models have thus far been presented as an ODE in
$\operatorname{\mathbb{R}}^{2}$ in the basis $S_{0}$, $S_{1}$. At generic
parameter values (i.e.,if $\kappa_{1}\neq\kappa_{2}$), we may change our basis
to $X_{0}=S_{1}$,
$X_{1}=\kappa_{1}(1-\pi)/(\kappa_{1}-\kappa_{2})S_{0}+S_{1}$. We then find
$\dot{X_{0}}=-\kappa_{1}X_{0}$ and $\dot{X_{1}}=-\kappa_{2}X_{1}$ for the
Linear ERK model and $\dot{X_{0}}=-\kappa_{1}AX_{0}$ and
$\dot{X_{1}}=-\kappa_{2}AX_{1}$ for the Rational ERK model.

As this change of basis is an invertible affine transformation, we may assume
that $q$ is a polynomial in variables $X_{0}$ and $X_{1}$. We note that this
diagonalisation of the ODE models would not have been possible if we had not
removed $S_{2}$ from our systems, as $S_{2}$ decouples.

We will write

$q(X_{0},X_{1})=\sum_{i,j\geq 0}a_{ij}X_{0}^{i}X_{1}^{j}=0.$

Taking the derivative of this equation with respect to $t$ gives

$\sum_{i,j\geq
0}a_{ij}(i\dot{X_{0}}X_{0}^{i-1}X_{1}^{j}+j\dot{X_{1}}X_{0}^{i}X_{1}^{j-1})=0,$

which, after substitution for the derivative variables, yields

$q^{\prime}(X_{0},X_{1}):=\sum_{i,j\geq
0}a_{ij}(i\kappa_{1}+j\kappa_{2})X_{0}^{i}X_{1}^{j}=0$

(for the Rational ERK model, we need to divide by $A$ to obtain such
$q^{\prime}$).

Note that, given our initial condition, both $X_{0}$ and $X_{1}$ vary, hence
the intersection of $V(q)$ and $V(q^{\prime})$ must contain a smooth point, as
both contain the ODE trajectory. This would imply that $q^{\prime}$ is a
constant multiple of $q$, which is not true at generic parameter values.
Hence, no such $q$ can exist.

For the Full ERK model, assume that there is a non-zero polynomial
$p\in\mathbb{C}(\theta)[S_{0},S_{1},C_{1},C_{2}]$ such that
$p(S_{0},S_{1},C_{1},C_{2})=0$ for all $t\geq 0$ given our initial condition.
Let $p$ have degree $n$. Again, we may assume without loss of generality that
$n$ is minimal. We know that there are no conserved quantities in the Full ERK
model (after removing $E$ and $S_{2}$) and thus $n\geq 2$. Let

$p(S_{0},S_{1},C_{1},C_{2})=\sum_{i,j,k,l\geq
0}a_{ijkl}S_{0}^{i}S_{1}^{j}C_{1}^{k}C_{2}^{l}=0.$

By taking the derivative with respect to $t$, we get

$\sum_{i,j,k,l\geq
0}a_{ijkl}(i\dot{S_{0}}S_{0}^{i-1}S_{1}^{j}C_{1}^{k}C_{2}^{l}+j\dot{S_{1}}S_{0}^{i}S_{1}^{j-1}C_{1}^{k}C_{2}^{l}+k\dot{C_{1}}S_{0}^{i}S_{1}^{j}C_{1}^{k-1}C_{2}^{l}+l\dot{C_{2}}S_{0}^{i}S_{1}^{j}C_{1}^{k}C_{2}^{l-1})=0.$

For $S_{0}$, this rearranges to

$\dot{S_{0}}=-\frac{\sum_{i,j,k,l\geq
0}a_{ijkl}(j\dot{S_{1}}S_{0}^{i}S_{1}^{j-1}C_{1}^{k}C_{2}^{l}+k\dot{C_{1}}S_{0}^{i}S_{1}^{j}C_{1}^{k-1}C_{2}^{l}+l\dot{C_{2}}S_{0}^{i}S_{1}^{j}C_{1}^{k}C_{2}^{l-1})}{\sum_{i,j,k,l\geq
0}a_{ijkl}iS_{0}^{i-1}S_{1}^{j}C_{1}^{k}C_{2}^{l}},$

(assuming there are $S_{0}$ terms in $p$) and we can derive similar
expressions for the other variables. Note that the denominator is of a smaller
degree than $p$, hence it will be non-zero at almost all points of our ODE
trajectory. In addition, for our given initial condition, all of our variables
vary and hence the numerator must be a non-zero polynomial and must not vanish
on the ODE trajectory. Denote $i^{*}$ the highest power of $S_{0}$ in $p$.
Then $a_{i^{*}jkl}\neq 0$ implies $k=0$, as otherwise the highest power of
$S_{0}$ in the numerator is $i^{*}+1$ and $i^{*}-1$ in the denominator,
implying that we can find an $S_{0}^{2}$-term in $\dot{S_{0}}$, a
contradiction. We can derive similar statements for $S_{1}$, $C_{1}$, and
$C_{2}$. Then, in the above equation, the highest power of $C_{1}$ in the
numerator is $k^{*}+1$, while the highest power in the denominator is
$k^{*}-1$ (at generic parameter values), implying that $\dot{S_{0}}$ is
quadratic in $C_{1}$, a contradiction. If $p$ does not contain terms in
$S_{0}$, we can apply a similar argument to $S_{1}$, $C_{1}$ or $C_{2}$.

In conclusion, no such $p$ can exist. ∎

###  A.5. The Linear ERK model is globally structurally identifiable

######  Proposition 20.

For any choice of three time points $0<t_{1}<t_{2}<t_{3}$ the model prediction
map $\phi_{t_{1},t_{2},t_{3}}$ is injective and so the Linear ERK model is
globally structurally identifiable.

###### Proof.

For a parameter $(\kappa_{1},\kappa_{2},\pi)$, we denote the analytic solution
in the Linear ERK model for species $i$ at time $t$ by ${S}_{i}(t)$ as in
Subsection 4.4. Then the model prediction map $\phi_{t_{1},t_{2},t_{3}}$ is
given by

$(\kappa_{1},\kappa_{2},\pi)\mapsto({S}_{0}(t_{1}),{S}_{1}(t_{1}),{S}_{2}(t_{1}),{S}_{0}(t_{2}),{S}_{1}(t_{2}),{S}_{2}(t_{2}),{S}_{0}(t_{3}),{S}_{1}(t_{3}),{S}_{2}(t_{3}))$

Suppose that $(\kappa_{1},\kappa_{2},\pi)$ and
$(\kappa_{1}^{\prime},\kappa_{2}^{\prime},\pi^{\prime})$ are two parameters
such that

$\phi_{t_{1},t_{2},t_{3}}(\kappa_{1},\kappa_{2},\pi)=\phi_{t_{1},t_{2},t_{3}}(\kappa_{1}^{\prime},\kappa_{2}^{\prime},\pi^{\prime}).$

Looking at the first component we find that
$5e^{-\kappa_{1}t_{1}}=5e^{-\kappa_{1}^{\prime}t_{1}}$, and since $t_{1}\neq
0$ it follows that $\kappa_{1}=\kappa_{1}^{\prime}$.

There are three cases to consider. The first case is if
$\kappa_{2}=\kappa_{1}=\kappa_{2}^{\prime}$. In this case, looking at the
second component, we get that
$5\kappa_{1}(1-\pi)t_{1}e^{-\kappa_{1}t_{1}}=5\kappa_{1}(1-\pi^{\prime})t_{1}e^{-\kappa_{1}t_{1}}$,
and so $\pi=\pi^{\prime}$, since $t_{1}\neq 0$.

Without loss of generality, the second case is if
$\kappa_{2}=\kappa_{1}\neq\kappa_{2}^{\prime}$. Suppose for a contradiction
that
$(\kappa_{1},\kappa_{1},\pi)\neq(\kappa_{1},\kappa_{2}^{\prime},\pi^{\prime})$.
Let ${S}_{1}(t)$ be the analytic solution of the Linear ERK model for species
1 at time $t$ with parameter $(\kappa_{1},\kappa_{1},\pi)$ and
${S}^{\prime}_{1}(t)$ be the analytic solution at time $t$ with parameter
$(\kappa_{1},\kappa^{\prime}_{2},\pi^{\prime})$. For the moment, consider
$t>0$ to be a variable. Then ${S}_{1}(t)={S^{\prime}}_{1}(t)$ is equivalent to

$(1-\pi)te^{-\kappa_{1}t}=\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa_{2}^{\prime}}\left(e^{-\kappa^{\prime}_{2}t}-e^{-\kappa_{1}t}\right).$
(16)

Dividing both sides by $e^{-\kappa_{1}t}$ and rearranging the above gives

$(1-\pi)t-\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa_{2}^{\prime}}e^{-(\kappa_{2}^{\prime}-\kappa_{1})t}=-\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa_{2}^{\prime}}$

Taking the derivative with respect to $t$ yields

$(1-\pi)=(1-\pi^{\prime})e^{-(\kappa_{2}^{\prime}-\kappa_{1})t}.$

Rearranging and taking a $\log$ then gives

$t=\frac{1}{\kappa_{1}-\kappa_{2}^{\prime}}\log\left(\frac{1-\pi}{1-\pi^{\prime}}\right).$

As the derivative of Equation (16) has exactly one solution in $t$, Equation
(16) has at most two solutions in $t$ by Rolle’s theorem. Therefore, the
second, fifth and eighth components of $\phi_{t_{1},t_{2},t_{3}}$ cannot take
the same value. This is a contradiction, and so we should have
$(\kappa_{1},\kappa_{1},\pi)=(\kappa_{1},\kappa_{2}^{\prime},\pi^{\prime})$,
meaning that this case is simply not possible.

We now consider the third and final case,
$\kappa_{2}\neq\kappa_{1}\neq\kappa_{2}^{\prime}$. Suppose for a contradiction
that
$(\kappa_{1},\kappa_{2},\pi)\neq(\kappa_{1},\kappa_{2}^{\prime},\pi^{\prime})$.
Let ${S}_{1}(t)$ be the analytic solution of the Linear ERK model for species
1 at time $t$ with parameter $(\kappa_{1},\kappa_{2},\pi)$ and
${S^{\prime}}_{1}(t)$ be the analytic solution at time $t$ with parameter
$(\kappa_{1},\kappa^{\prime}_{2},\pi^{\prime})$. For the moment, consider
$t>0$ to be a variable. Then ${S}_{1}(t)={S}^{\prime}_{1}(t)$ is equivalent to

$\frac{1-\pi}{\kappa_{1}-\kappa_{2}}\left(e^{-\kappa_{2}t}-e^{-\kappa_{1}t}\right)=\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa_{2}^{\prime}}\left(e^{-\kappa^{\prime}_{2}t}-e^{-\kappa_{1}t}\right).$
(17)

Dividing both sides by $e^{-\kappa_{1}t}$ and rearranging the above gives

$\frac{1-\pi}{\kappa_{1}-\kappa_{2}}e^{-(\kappa_{2}-\kappa_{1})t}-\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa^{\prime}_{2}}e^{-(\kappa^{\prime}_{2}-\kappa_{1})t}=\frac{1-\pi}{\kappa_{1}-\kappa_{2}}-\frac{1-\pi^{\prime}}{\kappa_{1}-\kappa^{\prime}_{2}}.$

Taking the derivative with respect to $t$ yields

$(1-\pi)e^{-(\kappa_{2}-\kappa_{1})t}-(1-\pi^{\prime})e^{-(\kappa^{\prime}_{2}-\kappa_{1})t}=0.$

Taking a $\log$ and rearranging gives

$t=\frac{1}{\kappa_{1}-\kappa_{2}^{\prime}}\log\left(\frac{1-\pi}{1-\pi^{\prime}}\right).$

As the derivative of Equation (17) has exactly one solution in $t$, Equation
(17) has at most two solutions in $t$ by Rolle’s theorem. Therefore, the
second, fifth and eighth components of $\phi_{t_{1},t_{2},t_{3}}$ cannot take
the same value. This is a contradiction, and so we must have
$(\kappa_{1},\kappa_{2},\pi)=(\kappa_{1},\kappa_{2}^{\prime},\pi^{\prime})$.

∎