# Adaptive Gated Graph Convolutional Network for Explainable Diagnosis of
Alzheimer’s Disease using EEG Data

Dominik Klepl  Centre for Computational Science and Mathematical Modelling,
Coventry University, Coventry CV1 2JH, UK  Fei He  Corresponding authors: Fei
He<EMAIL_ADDRESS>and Min Wu<EMAIL_ADDRESS>Centre for
Computational Science and Mathematical Modelling, Coventry University,
Coventry CV1 2JH, UK  Min Wu  Corresponding authors: Fei He
<EMAIL_ADDRESS>and Min Wu<EMAIL_ADDRESS>Institute for
Infocomm Research (I${}^{\text{2}}$R), Agency for Science, Technology and
Research (A*STAR), Singapore)  Daniel J. Blackburn  Department of
Neuroscience, University of Sheffield, Sheffield, S10 2HQ, UK  Ptolemaios G.
Sarrigiannis  Department of Neurophysiology, Royal Devon and Exeter NHS
Foundation Trust, Exeter, EX2 5DW, UK

###### Abstract

Graph neural network (GNN) models are increasingly being used for the
classification of electroencephalography (EEG) data. However, GNN-based
diagnosis of neurological disorders, such as Alzheimer’s disease (AD), remains
a relatively unexplored area of research. Previous studies have relied on
functional connectivity methods to infer brain graph structures and used
simple GNN architectures for the diagnosis of AD. In this work, we propose a
novel adaptive gated graph convolutional network (AGGCN) that can provide
explainable predictions. AGGCN adaptively learns graph structures by combining
convolution-based node feature enhancement with a well-known correlation-based
measure of functional connectivity. Furthermore, the gated graph convolution
can dynamically weigh the contribution of various spatial scales. The proposed
model achieves high accuracy in both eyes-closed and eyes-open conditions,
indicating the stability of learned representations. Finally, we demonstrate
that the proposed AGGCN model generates consistent explanations of its
predictions that might be relevant for further study of AD-related alterations
of brain networks.

Keywords: Alzheimer’s disease, graph neural network, classification, EEG

##  I Introduction

The brain is a complex, densely connected system that operates across multiple
spatial and temporal scales. Neurological diseases, such as Alzheimer’s
disease (AD), can alter the connectivity of the brain and thus disrupt brain
function delbeuck2003disrupted ; de2012disrupted ; pievani2011functional ;
konig2005decreased . AD is the most common cause of dementia and affects
millions of patients worldwide. Currently, the diagnosis of AD is typically
made using a combination of cognitive and neurological assessments, as well as
neuroimaging techniques, such as positron emission tomography (PET) or
magnetic resonance imaging (MRI), which can be time-consuming and expensive.
The development of rapid, economical, and explainable diagnosis methods is
becoming increasingly important.

Electroencephalography (EEG) is an economical and non-invasive neuroimaging
method that records the sum of electrical potentials generated by various
brain areas. EEG is extensively used in the research of AD-related alterations
in brain function and functional connectivity. Although EEG is not yet widely
used in clinical settings, numerous studies have demonstrated the high
effectiveness of an EEG-based diagnosis of AD klepl2022GNN ;
klepl2022bispectrum ; klepl2021energylandscape ; shan2022STGCN-AD ;
nobukawa2020classificationFC ; fouladi2022contwaveletCNN .

AD causes disruption of synaptic connections across multiple scales
jeong2004eeg ; pievani2011functional ; klepl2021bispectrum  and can thus be
viewed as a network disorder delbeuck2003disrupted . The synaptic
disconnection can be observed in EEG signals as alterations of synchronisation
and functional connectivity (FC) pievani2011functional ;
babiloni2016synchronisationreview . Furthermore, the slowing of EEG signals is
a reliable characteristic of AD jeong2004eeg ; ghorbanian2015slowing ,
observed as a shift of spectral power towards low-frequency components. Graph-
theoretic studies of AD also report reduced complexity, disruption of small-
world properties, decreased integration, and increased segregation
supekar2008smallworld ; kabbara2018integrationsegregation ;
vecchio2021coherence-smallworld ; stam2007correlation ; klepl2021bispectrum ;
klepl2021energylandscape . However, one of the challenges in EEG-based
predictive models is the efficient utilisation of the information collected
over multiple electrodes, since there is information to be gained both at the
electrode level, e.g. frequency spectrum, and the cross-electrode level, e.g.
FC.

Machine learning-based approaches often require domain knowledge and rely on
manual feature extraction. For example, Oltu et al. oltu2021PSD-coherence
calculate power spectrum density (PSD) and coherence across multiple EEG
electrodes, and then use descriptive statistics, such as sum and variance, as
input features. Other feature-based methods use FC
nobukawa2020classificationFC ; song2018FC ; yu2019fuzzyFC . These methods
first reconstruct the brain graph via measures of FC, such as phase lagging
index nobukawa2020classificationFC , generalised composite multiscale entropy
vector song2018FC , or phase synchronisation index yu2019fuzzyFC . The
features can then be extracted via statistics song2018FC  or graph-theoretic
measures nobukawa2020classificationFC ; yu2019fuzzyFC .

In contrast, deep learning methods can extract features automatically from the
input. However, utilising the information from multiple electrodes with
classical deep learning methods is challenging. To overcome this issue,
several studies have transformed EEG signals into images to make use of
convolutional neural networks (CNN) ieracitano2019CNN ; deepthi2020CNN ;
bi2019spectralimage ; huggins2021deep ; fouladi2022contwaveletCNN , which are
efficient in image classification. For instance, Ieracitano et al.
ieracitano2019CNN  compute the PSD across channels and compose them to form a
channel by PSD image. Bi et al. bi2019spectralimage  use spectral topology
images and leverage the colour channels of an image to represent three
frequency bands. Finally, Huggins et al. huggins2021deep  create tiled images
where each tile contains the continuous wavelet transform of an EEG electrode.
Although these methods utilise multiple channels, the cross-electrode
information is still omitted. To address this limitation, a CNN trained on FC-
based adjacency matrices has been proposed alves20FC-CNN . However, CNN is not
well suited for such input since the adjacency matrix is irregular and non-
euclidean.

Graph neural network (GNN) is an extension of CNN to process graph-structured
inputs. Multiple studies propose GNN-based architectures to process EEG.
However, GNN methods for EEG-based diagnosis of AD are limited klepl2022GNN ;
shan2022STGCN-AD . GNN-EEG implementations often include several steps: (1)
input construction, i.e. graph structure and node features; (2) GNN encoder to
learn node embeddings; and (3) aggregation of node embeddings to a graph
embedding, which can be used in the final classification step.

There are various approaches to realise the graph construction in step (1).
Node features are commonly defined as EEG time-series signal shan2022STGCN-AD
; li2021multidomain ; asadzadeh2022bayes ; li2022bimodal , or a statistical
summary of the signal in the time domain raeisi2022neonatal-epilepsy ;
tang2021selfsupervised , the frequency domain klepl2022GNN ;
sun2022transformer , or the differential entropy li2021multidomain ;
yin2021LSTMfusion ; chen2021bayesianGNN ; jia2020graphsleepnet ; zeng2022siam
; sun2022transformer . Based on network neuroscience literature, many
approaches define the brain graph using FC measures klepl2022GNN ;
shan2022STGCN-AD ; liu2022MST-GNN ; chang2021MMN-schizophrenia ;
raeisi2022neonatal-epilepsy ; tang2021selfsupervised ; li2021multidomain ;
li2022bimodal . The graph structure can also be based on the distance between
EEG electrodes tang2021selfsupervised ; yin2021LSTMfusion ;
chen2021bayesianGNN . However, such an approach largely ignores brain
connectivity information. Alternatively, the brain graph can be automatically
learned by the model, either as a learnable mask shared across samples
li2019GNN-motor-movement ; li2021multidomain ; sun2022transformer  or by
pairwise node feature distance minimisation regularised by an additional graph
loss function zhong2020GNN-regularised-emotion ; jia2020graphsleepnet ;
zeng2022siam . While such approaches are flexible and should converge to an
optimal graph structure with respect to a given learning task, the learned
brain graph might not be representative of the underlying brain connectivity,
thus limiting the interpretability of such learned brain graphs. In this work,
we propose an adaptive graph learning mechanism based on node feature
enhancement via CNN and subsequent graph construction. This is achieved by
using a standard FC measure (Pearson’s correlation) and sparsified via
k-nearest neighbour (KNN) edge selection. Thus, it combines the strength of
the FC-based and automated graph learning methods.

The design of GNN encoders in step (2) for EEG applications has been mainly
limited to simple architectures, such as the Chebyshev graph convolution
(ChebConv) chang2021MMN-schizophrenia ; asadzadeh2022bayes ;
tang2021selfsupervised ; yin2021LSTMfusion ; chen2021bayesianGNN ;
jia2020graphsleepnet ; li2022bimodal , and simple graph convolution (GCN)
klepl2022GNN ; li2019GNN-motor-movement ; zhdanov2022schizophrenia ;
li2021multidomain ; sun2022transformer ; zhong2020GNN-regularised-emotion .
However, we hypothesise that such node embedding updating mechanisms are not
optimal for EEG tasks. These graph convolutions update node embeddings by
summing the initial embedding and the aggregated messages from the
neighbouring nodes. Such updating implies that information from different
scales contributes equally to the final node embeddings, hence graph
embeddings as well. While brain disruptions caused by AD occur across multiple
spatial scales, their predictive power is likely different. Therefore, a
gating mechanism is crucial for filtering and weighting the information
collected across different scales. We propose to adopt the gated graph
convolution li2016GRGCN  to address this issue.

Finally, we implement the aggregation of node embeddings in step (3) by
adopting the adaptive structure-aware pooling (ASAP) node pooling mechanism
ranjan2020asap  to first learn the most important clusters of nodes, which are
in turn concatenated to form the graph embedding. This is in contrast to the
previous approaches that do not use any node pooling and form graph embeddings
via simple element-wise readout layers klepl2022GNN ; xu2022dagam ;
raeisi2022neonatal-epilepsy ; li2021multidomain ; zhdanov2022schizophrenia ;
liu2022MST-GNN ; zhong2020GNN-regularised-emotion  or concatenating all nodes
of the graph shan2022STGCN-AD ; chang2021MMN-schizophrenia . Other node
pooling approaches were tested for EEG applications demir2021ERP-comparison ;
xu2022dagam . In contrast to ASAP pooling, these approaches pool the graph by
selecting a specified number of nodes without considering their local context
within the graph. Therefore, important information might be lost due to such
node pooling.

In this paper, we propose a novel GNN model for explainable AD classification,
which can adaptively enhance node features and dynamically construct brain
graph structures as shown in Fig. 1. The learned brain graphs can then be used
for the interpretation of predictions. Moreover, a clustering-based node
pooling mechanism is adopted to coarsen the brain graph, thus localising the
brain regions that contribute to the predictions. Finally, we carry out
extensive ablation and parameter sensitivity experiments to elucidate the
importance of the individual blocks within the proposed model architecture.

Figure 1: Architecture of the proposed adaptive gated graph convolutional
network. Node features are defined as power spectral density from 1 to 45 Hz.
The node features are then used as input to the graph learning module (green
box), where they are enhanced by a 1D convolutional neural network. The brain
graph structure is then constructed as a correlation graph and made sparse by
a k-nearest-neighbour edge selection. The enhanced node features and the
learned graph structure are then passed to an encoder (purple box) consisting
of a gated graph convolutional layer (repeated for R iterations) and the ASAP
node pooling mechanism. The node pooling coarsens the graph. The node features
of the coarsened graph are then flattened and passed to a multilayer-
perceptron classifier, which outputs the predictions.

##  II Data

EEG recordings were collected from 20 AD patients and 20 healthy control
participants (HC) younger than 70 years. A detailed description of the
experimental design and confirmation of the diagnosis is provided in
blackburn2018synchronisation . All the AD participants were recruited from the
Sheffield Teaching Hospital memory clinic. AD participants were diagnosed
between one month and two years before data collection. All of them were in
the mild to moderate stage of the disease at the time of recording, with an
average Mini Mental State Examination (MMSE) score of 20.1 (sd = 4). High-
resolution structural magnetic resonance imaging (MRI) scans of all patients
were acquired to eliminate alternative causes of dementia. Age and gender-
matched HC participants with normal neuropsychological tests and structural
MRI scans were recruited.

EEG data was acquired using an XLTEK 128-channel headbox, Ag/AgCL electrodes
with a sampling frequency of 2 kHz using a modified 10-10 overlapping a 10-20
international electrode placement system with a referential montage with a
linked earlobe reference. The recordings lasted 30 minutes, during which the
participants were instructed to rest and not think about anything specific.
Within each recording were two-minute-long epochs during which the
participants had their eyes closed, alternating with an equal duration of
eyes-open epochs.

All the recordings were reviewed by an experienced neurophysiologist on the
XLTEK review station with time-locked video recordings (Optima Medical LTD).
For each participant, three 12-second-long artefact-free epochs were isolated.
Finally, the following 23 bipolar channels were created: F8–F4, F7–F3, F4–C4,
F3–C3, F4–FZ, FZ–CZ, F3–FZ, T4–C4, T3–C3, C4–CZ, C3–CZ, CZ–PZ, C4–P4, C3–P3,
T4–T6, T3–T5, P4–PZ, P3–PZ, T6–O2, T5–O1, P4–O2, P3–O1 and O1–O2
blackburn2018synchronisation .

As a neurophysiologist confirmed the EEG signal to be artefact-free, we did
not perform any further cleaning of the signals. The signals are filtered
using a band-pass Butterworth filter to a range of 0.5 Hz and 45 Hz and down-
sampled to 200 Hz. Finally, 1-second long windows with 50% overlap are created
to increase the sample size.

##  III Methods

The proposed adaptive gated graph convolutional network (AGGCN) model consists
of three blocks, namely, a graph learning module, a GNN encoder and a
classifier. The graph learning module receives a node feature matrix as input,
enhances it using a 1D-CNN and learns the brain graph structure. The GNN
encoder then uses the output of the graph learning module as input, i.e. a
featured, weighted, undirected graph. The encoder generates a graph embedding
used by the classifier to output the predicted probabilities.

###  III.1 Node feature and graph learning

The node features are defined as power spectral density computed from
1-second-long EEG signals with 1 Hz increments from 1 to 45 Hz. Hence, the
input is a node feature matrix $X\in\mathbb{R}^{N\times D_{in}},D_{in}=45$.

The input is then passed to a convolutional neural network (CNN) with batch
normalisation, $L_{CNN}$ 1D convolutional layers and a maximum pooling with
kernel size 2 and step size 2. The output is flattened and fed to a fully
connected layer with hidden size $h_{CNN}$ and batch normalisation. This
neural network outputs a matrix of enhanced node features
$X^{\prime}\in\mathbb{R}^{N\times D_{h_{CNN}}}$.

A graph structure is then inferred from the enhanced node features by
computing the absolute value of Pearson’s correlation for each pair of nodes.
Thus, a unique graph structure is learned for each input sample and is defined
by an adjacency matrix $A\in\mathbb{R}^{N\times N}$ with $N=23$ being the
number of EEG channels. In order to produce sparse graphs, the k-nearest-
neighbours algorithm is utilised. This means that the $k$ strongest edges are
preserved for each node.

This proposed graph learning module has multiple hyperparameters that control
its architecture. Namely, these are the number of convolutional layers
$L_{CNN}$, the kernel size (which is equal to the step size), the number of
filters, the hidden size $h_{CNN}$, the dropout rate $drop_{CNN}$ and the
$k_{KNN}$ parameter that controls the graph sparsity.

###  III.2 Graph neural network encoder and classifier

A graph convolution extends the classical convolution from the Euclidean
domain to the graph domain. The input graph is given by $G=(N,A,X^{\prime})$
where $N$ is the set of nodes, $A$ is the learned graph, and $X^{\prime}$ is
the enhanced node feature matrix. A simple graph convolution is defined by the
message-passing mechanism wherein the node embedding of node $i$ is learned by
aggregating information from its 1-hop neighbourhood, i.e. nodes connected
with an edge, as follows:

$x^{l+1}_{i}=x^{l}_{i}+\Theta\sum_{j\in N(i)}e_{ij}x^{l}_{j},$ (1)

where $x^{l}_{i}$ are the node features of node $i$ at the $l$th layer,
$x^{0}_{i}$ is the $i$th row of the input node feature matrix $F^{\prime}$,
and $\Theta$ is a learnable linear transformation. $N(i)$ and $e_{ij}$ are the
neighbourhood of node $i$ and the edge weight connecting nodes $i$ and $j$
given by the adjacency matrix $A$, respectively. Stacking $L$ graph
convolutional layers then means aggregating information iteratively from 1-hop
to $L$-hop neighbourhoods, thus gradually going from local to global
information about the graph.

Note that the aggregated message is added to the initial node embedding
$x^{l}_{i}$. Thus, the entire information collected from each L-hop
neighbourhood is always fully integrated into the node embedding. However,
information might be distributed unequally across spatial scales in brain
graphs. The gated graph convolution (GGCN) li2016GRGCN  addresses this problem
by introducing a mechanism to decide what information should be retained at
each scale selectively:

$\displaystyle m_{i}^{(r+1)}$ $\displaystyle=\sum_{j\in
N(i)}e_{ji}\cdot\Theta^{r+1}\cdot x_{j}^{(r)},$ (2) $\displaystyle
x_{i}^{(r+1)}$ $\displaystyle=\textrm{GRU}(m_{i}^{(r+1)},x_{i}^{(r)}),$ (3)

where $m_{i}$ are the aggregated messages, $\sum$ is the aggregation function,
$\Theta^{r}$ is a learnable matrix for iteration $r$, which maps the node
features from shape $[1,D_{h_{CNN}}]$ to $[1,D_{h_{GNN}}$, and $GRU$ is the
gated recurrent unit cho2014GRU . Briefly, a GRU is a recurrent neural network
layer with an update, reset, and new gates that allow the network to
recursively update or forget information about the input. The node embeddings
are learned recursively up to $R$ iterations with a shared $GRU$ gate, which
is equivalent to stacking $R$ GCN layers.

The node embeddings are then passed through an activation function and a batch
normalisation layer. Finally, the node embeddings are passed to the node
pooling module. The hyperparameters of the proposed encoder are the number of
iterations $R$, the hidden size $h_{GNN}$, the activation function, the
aggregation function and the dropout rate $drop_{GNN}$ applied after the
encoder.

####  III.2.1 Node pooling

After learning the node embeddings, the model learns a coarsened graph using
the ASAP pooling mechanism ranjan2020asap . This pooling first learns $N$
clusters, each centred at one node, also named ego-graphs. The membership of a
node $j$ in the ego-cluster centred at node $i$ is given by the $S_{ij}$
matrix. Note that this is a soft-cluster assignment matrix; thus, each node
can belong to multiple clusters with varying membership strengths. The
clusters are learned as follows:

$\displaystyle S_{ij}$ $\displaystyle=a_{ij},$ (4) $\displaystyle a_{ij}$
$\displaystyle=\mathrm{softmax}\left(\theta^{\mathrm{T}}\sigma\left(\Theta
x^{m}_{i}\|x_{j}\right)\right),$ (5) $\displaystyle x^{m}_{i}$
$\displaystyle=\max_{j\in N(i)}x_{j},$ (6)

where $a_{ij}$ is the attention score and the membership strength, $\theta$
and $\Theta$ are learnable vector and matrix, respectively. $\sigma$ is the
LeakyReLU activation function, and $x^{m}_{i}$ is the master query
representing the initial cluster embedding. The attention scores are also
subject to a dropout probability $drop_{pool}$. The final cluster embedding is
then calculated as an attention-weighted sum, which is additionally weighted
by the cluster score $\phi_{i}$:

$x^{c}_{i}=\phi_{i}\sum_{j\in N(i)}a_{ij}x_{j},$ (7)

where the cluster score $\phi_{i}$ is computed by the local extremum graph
convolution ranjan2020asap :

$\phi_{i}=\Theta_{1}\cdot x_{i}+\sum_{j\in
N(i)}e_{ji}\cdot\left(\Theta_{2}x_{i}-\Theta_{3}x_{j}\right),$ (8)

which is designed to measure the relative importance of each cluster.

The cluster embedding $x^{c}_{i}$ is then used to select the top $k$ scoring
clusters, which will be included in the coarsened graph:

$\displaystyle\bar{i}$ $\displaystyle=Top_{k}(X^{c}),k\in[1,2,...N],$
$\displaystyle\bar{S}$ $\displaystyle=S(:,\bar{i})$ (9) $\displaystyle A^{p}$
$\displaystyle=\bar{S}^{\mathrm{T}}\cdot A\cdot\bar{S},$ $\displaystyle
X^{\mathrm{p}}$ $\displaystyle=X^{c}(:,\bar{i})$ (10)

where $Top_{k}$ is a function that returns the indices of clusters $\bar{i}$.
$\bar{S}$ and $X^{\mathrm{p}}$ is the pruned soft-cluster assignment matrix
and the pruned cluster embedding matrix, respectively, and $A^{p}$ is the
adjacency matrix of the coarsened graph.

The graph pooling module has the following hyperparameters: the size of the
pooled graph $k_{pool}$, the dropout rate $drop_{pool}$ and the negative slope
of the LeakyReLU activation.

####  III.2.2 Multilayer perceptron classifier

The cluster embedding matrix $X^{p}$ of the coarsened graph returned by the
node pooling module is flattened and fed to a multilayer perceptron (MLP)
classifier. Specifically, a $L_{MLP}$-layer MLP with hidden size $h_{MLP}$ is
utilised with a block of batch normalisation, activation function, and dropout
layers utilised between the fully connected layers. The final layer outputs a
two-dimensional vector of log probabilities for each class.

The classifier has the following hyperparameters: the number of layers
$L_{MLP}$, hidden size $h_{MLP}$, activation function and dropout rate
$drop_{MLP}$.

Table 1: Performance of the proposed AGGCN in eyes closed (EC), eyes open (EO) and combined (EC+EO) conditions. Condition | Accuracy | AUC | Sensitivity | Specificity
---|---|---|---|---
EC | 90.83 $\pm$ 1.18 | 94.93 $\pm$ 0.93 | 93.37 $\pm$ 1.86 | 88.24 $\pm$ 1.92
EO | 88.77 $\pm$ 1.67 | 93.75 $\pm$ 1.27 | 90.41 $\pm$ 3.21 | 87.06 $\pm$ 2.88
EC+EO | 89.61 $\pm$ 1.54 | 94.05 $\pm$ 1.03 | 89.65 $\pm$ 2.79 | 89.57 $\pm$ 2.22

###  III.3 Model implementation and evaluation

The proposed AGGCN model was implemented using PyTorch 1.10 pytorch , and
PyTorch Geometric 2.0.2 pytorch-geometric  and trained on a laptop with Intel
i7 CPU, 16 GB RAM and an NVIDIA RTX 2070 GPU.

The model is trained by minimising the binary cross-entropy loss. The model
performance is evaluated using repeated (100 times) 10-fold stratified cross-
validation and trained on the dataset collected during the eyes-closed
condition. Since there are multiple samples from the same participant, keeping
all samples within the same fold is crucial to prevent information leakage. A
stochastic gradient descent (SGD) optimiser and an exponential learning rate
scheduler are used to train the model with a batch size of 90 for 100 epochs.
Additionally, zero-mean Gaussian noise with standard deviation $\sigma$ is
added to the input during training with probability $p_{noise}$ to improve the
generalisability of the model. Eventually, the best model was identified using
the average cross-validated accuracy measured on the test folds. The selected
model was then tested on the dataset obtained during the eyes-open condition
as well as on the combined dataset from both conditions.

Note that the hyperparameters of the proposed model are optimised using
Bayesian optimisation. Ten warm-up random iterations were used to initialise
the optimisation, followed by 100 optimisation iterations. Moreover, we carry
out parameter-sensitivity experiments to verify the influence of a few key
hyperparameters of the proposed model architecture. Specifically, these are
the number of iterations of the GGCN encoder, the size of the pooled graphs,
the sparsity of the learned graph and the choice of aggregation function of
the GGCN encoder. The hyperparameters of the model are reported in our
supplementary materials.

##  IV Results and Discussion

In this section, we demonstrate the experimental results of our AGGCN. As
illustrated in Table 1, our AGGCN has shown robustness across all the
conditions, indicating its potential for real-world applications. Note that
the highest accuracy was reached during the EC condition. This is likely
because with eyes closed, the ocular artefacts are minimised; thus, the
underlying dynamics are easier to detect. However, the performance remains
high even in EO and EC+EO conditions suggesting that the proposed model can
detect underlying patterns in all the conditions. In addition, the
hyperparameter values of the optimised model are reported in Table S1 in
Supplementary Materials.

###  IV.1 Comparison with the baselines

The proposed model was compared to three baseline models that were proposed in
the literature across the three conditions. The first baseline is the best-
performing model from our previous work klepl2022GNN . It is a GNN with two
spatial graph convolutional layers, maximum readout and brain graph defined
using the amplitude-envelope-correlation (AEC-GNN). The second baseline model
is the spatio-temporal GNN (STGCN) that uses temporal convolutions and
ChebConv layers and defines the brain graphs using wavelet coherence
shan2022STGCN-AD . Finally, we use a naive baseline model where the input node
feature matrix is flattened and classified using an MLP klepl2022GNN  without
using graph-domain information.

Table 2 shows the AUC values of various methods across different conditions.
Note that all four methods were evaluated under the same setting (e.g. the
same 1-second EEG window). We can observe that our proposed AGGCN outperforms
the baselines significantly. Moreover, STGCN was evaluated using non-
stratified cross-validation in their original paper. It is well expected that
its performance drops significantly when evaluated using stratified cross-
validation in our experiments.

Table 2: The AUC values of the baseline models and the proposed method across conditions. The best-performing model is highlighted in bold. Model | EC | EO | EC+EO
---|---|---|---
AEC-GNN klepl2022GNN  | 87.15 $\pm$ 2.45 | 81.26 $\pm$ 2.58 | 83.16 $\pm$ 1.84
MLP klepl2022GNN  | 89.09 $\pm$ 0.69 | 86.43 $\pm$ 1.03 | 89.21 $\pm$ 0.94
STGCN shan2022STGCN-AD  | 69.50 $\pm$ 3.74 | 63.82 $\pm$ 2.53 | 67.16 $\pm$ 2.68
Proposed | 94.93 $\pm$ 0.93 | 93.75 $\pm$ 1.27 | 94.05 $\pm$ 1.03

###  IV.2 Model ablation study

We perform ablation experiments to determine the contribution of each module
of the proposed model. The following seven ablated variants of the proposed
model were tested in our experiments.

  * •

A: no node pooling;

  * •

B: graph learning replaced with a fully connected graph;

  * •

C: GGCN replaced with a $R$th-order ChebConv;

  * •

D: variants A and B combined;

  * •

E: variants A and C combined;

  * •

F: variants B and C combined;

  * •

G: variants A, B and C combined.

Figure 2: Accuracy of model variants. The asterisks report the p-value of a
nonparametric Mann-Whitney U test measuring the difference between AGGCN and
the ablated variants.

The ablation results in Fig. 2 reveal that each of the proposed modules
contributes significantly to the high performance of the proposed
architecture. For variant A, we can observe that the contribution of the node
pooling module is significant, albeit relatively small. However, this module
reduces the number of parameters of the model and helps to produce explainable
predictions (Fig. 4 and Fig. 5 in the next subsection). Without the node
pooling, the final MLP classifier would have $N\times h_{GNN}\times h_{MLP}$
parameters ($N=23$), but node pooling reduces it to $k_{pool}\times
h_{GNN}\times h_{MLP}$ ($k_{pool}=5$). For variant B, it is not surprising
that its performance decreases significantly, as the graph learning module is
replaced with a fully connected graph and variant B cannot leverage any graph-
domain information.

Next, we demonstrate that the GGCN encoder improves performance significantly
compared to a ChebConv encoder according to variant C. A ChebConv layer is
similar to a GGCN in its iterative nature, i.e. ChebConv iteratively updates
node embeddings by approximating the eigendecomposition of graph Laplacian.
However, ChebConv does not have any gating mechanism, which means that
information across scales contributes to the final embedding equally. Since
all of the major modules of the proposed are shown to contribute to the final
performance significantly, it is unsurprising that the rest of the ablated
models with more than one of these modules perform significantly worse as well
(Variants D-G in Fig. 2).

The parameter sensitivity experiments also confirm the optimal values of
crucial hyperparameters of the proposed model (Supplementary Materials, Figs
S1-S4). It is worth noting that the proposed architecture allows training
relatively deep models (using up to ten GGCN iterations) with only a minor
performance decrease (Fig. S1).

###  IV.3 Explainability of AGGCN

The proposed model generates plausible and consistent explanations for its
predictions. Specifically, the graph learning module learns a clear difference
between the AD and HC cases, as shown in Fig. 3. The learned brain graphs show
that AD cases have increased connectivity globally, while HC graphs seem more
sparse with few strongly connected clusters. A well-defined cluster is present
in both groups within the parietal, temporal and occipital regions.

Figure 3: Average adjacency matrix of learned graphs of AD and HC cases.
Figure 4: Average adjacency matrix of pooled graphs of AD and HC cases.

The node pooling mechanism then highlights the importance of this cluster, as
most of the nodes within this cluster are included in the coarsened graph in
Fig. 4. The coarsened graph also captures the pattern of increased
connectivity of AD cases compared to HC cases. Note that the nodes of the
pooled graphs are, in fact, cluster representations, i.e. attention weighted
sum of node embeddings as shown in Fig. 5. These clusters are primarily
located in the parietal and temporal channels and a few central and frontal
channels. Surprisingly, the cluster assignments are quite similar in AD and HC
in Fig. 5 and are similar to the learned graph structure shown in Fig. 3.

Finally, the role of the gating mechanism is further elucidated by analysing
the amount of information gathered at each scale, i.e. iteration of GGCN (Fig.
6). We measure this by computing the average Euclidean distance between the
initial and updated node embedding at each iteration, i.e. $x_{i}^{(r)}$ and
$m_{i}^{(r+1)}$ in Eq. 3. For instance, a small distance means a small amount
of information was gathered from that scale, i.e. iteration. Local information
contributes highly to the node embeddings of the HC cases, and then the degree
of contributions linearly decreases with increasing graph scale. The opposite
pattern is observed for AD cases, where the later iterations influence the
node embeddings. This highlights the distributed nature of the neural
disruptions caused by AD as the later iterations gather more global and
distributed information.

Figure 5: Average attention scores of preserved clusters in the pooled graphs
of AD and HC cases. Figure 6: The average distance between initial node
embedding and updated node embeddings shows the amount of information retained
in each iteration of GGCN, i.e. going from local to global information. The
asterisks denote the p-value of statistical tests comparing the average
distance between AD and HC cases.

###  IV.4 Limitations and future work

Although our approach achieves competitive performance, we identify a few
drawbacks. First, the relatively small size of our dataset imposes a limit on
fitting complex models. We address this issue by segmenting the EEG signals
into short windows. The short window length means that the model might not be
able to represent information from low-frequency components of the signal.

Next, we do not explore alternative node feature representations beyond PSD in
this study. PSD is merely a linear frequency-domain representation of the
signal. Including time-domain and nonlinear information in the node features
might improve the expressiveness of the model. Similarly, the proposed graph
learning mechanism is limited to linear patterns of FC, because (1) it is
inferred from the node features and (2) it is expressed as Pearson’s
correlation coefficient. Future work should explore other forms of FC that
might be integrated into the graph learning mechanism; and study ways to
include more complex frequency-dependent coupling information.

##  V Conclusion

This work proposes a novel graph learning model that performs highly in the AD
diagnosis task. Additionally, we show that the model produces robust and
clinically relevant explanations for its predictions via the novel graph
structure learning module and the node pooling mechanism. Finally, we
highlight the importance of utilising the gating mechanism within a message-
passing encoder. This allows the model to accurately represent the multiscale
distributed network disruptions displayed in the AD cases.

## Acknowledgement

The EEG data was funded by a grant from Alzheimer’s Research UK (ARUK-
PPG20114B-25). The views expressed are those of the author(s) and not
necessarily those of the NHS, the NIHR or the Department of Health. The work
was supported by A*STAR, AI, Analytics and Informatics (AI3) Horizontal
Technology Programme Office (HTPO) seed grant C211118015.

## References

  * [1] Xavier Delbeuck, Martial Van der Linden, and Fabienne Collette.  Alzheimer’disease as a disconnection syndrome?  Neuropsychology review, 13(2):79–92, 2003.
  * [2] Willem de Haan, Wiesje M van der Flier, T Koene, Lieke L Smits, Philip Scheltens, and Cornelis J Stam.  Disrupted modular brain dynamics reflect cognitive dysfunction in Alzheimer’s disease.  Neuroimage, 59(4):3085–3093, 2012.
  * [3] Michela Pievani, Willem de Haan, Tao Wu, William W Seeley, and Giovanni B Frisoni.  Functional network disruption in the degenerative dementias.  The Lancet Neurology, 10(9):829–843, 2011.
  * [4] Thomas König, L Prichep, Thomas Dierks, Daniela Hubl, LO Wahlund, ER John, and V Jelic.  Decreased eeg synchronization in alzheimer’s disease and mild cognitive impairment.  Neurobiology of aging, 26(2):165–171, 2005.
  * [5] Dominik Klepl, Fei He, Min Wu, Daniel J Blackburn, and Ptolemaios Sarrigiannis.  Eeg-based graph neural network classification of alzheimer’s disease: An empirical evaluation of functional connectivity methods.  IEEE Transactions on Neural Systems and Rehabilitation Engineering, 30:2651–2660, 2022.
  * [6] Dominik Klepl, Fei He, Wu Min, Daniel Blackburn, and Ptolemaios Sarrigiannis.  Bispectrum-based cross-frequency functional connectivity: Classification of alzheimer’s disease.  In 2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 305–308. IEEE, 2022\.
  * [7] Dominik Klepl, Fei He, Min Wu, Matteo De Marco, Daniel J Blackburn, and Ptolemaios G Sarrigiannis.  Characterising alzheimer’s disease with eeg-based energy landscape analysis.  IEEE Journal of Biomedical and Health Informatics, 26(3):992–1000, 2021.
  * [8] Xiaocai Shan, Jun Cao, Shoudong Huo, Liangyu Chen, Ptolemaios Georgios Sarrigiannis, and Yifan Zhao.  Spatial–temporal graph convolutional network for alzheimer classification based on brain functional connectivity imaging of electroencephalogram.  Human Brain Mapping, 2022.
  * [9] Sou Nobukawa, Teruya Yamanishi, Shinya Kasakawa, Haruhiko Nishimura, Mitsuru Kikuchi, and Tetsuya Takahashi.  Classification methods based on complexity and synchronization of electroencephalography signals in alzheimer’s disease.  Frontiers in psychiatry, 11:255, 2020.
  * [10] Saman Fouladi, Ali A Safaei, Nadia Mammone, Foad Ghaderi, and MJ Ebadi.  Efficient deep neural networks for classification of alzheimer’s disease and mild cognitive impairment from scalp eeg recordings.  Cognitive Computation, pages 1–22, 2022.
  * [11] Jaeseung Jeong.  Eeg dynamics in patients with alzheimer’s disease.  Clinical neurophysiology, 115(7):1490–1505, 2004.
  * [12] Dominik Klepl, Fei He, Min Wu, Daniel J Blackburn, and Ptolemaios G Sarrigiannis.  Cross-frequency multilayer network analysis with bispectrum-based functional connectivity: A study of alzheimer’s disease.  bioRxiv, pages 2021–08, 2021.
  * [13] Claudio Babiloni, Roberta Lizio, Nicola Marzano, Paolo Capotosto, Andrea Soricelli, Antonio Ivano Triggiani, Susanna Cordone, Loreto Gesualdo, and Claudio Del Percio.  Brain neural synchronization and functional coupling in alzheimer’s disease as revealed by resting state eeg rhythms.  International Journal of Psychophysiology, 103:88–102, 2016.
  * [14] Parham Ghorbanian, David M Devilbiss, Terry Hess, Allan Bernstein, Adam J Simon, and Hashem Ashrafiuon.  Exploration of eeg features of alzheimer’s disease using continuous wavelet transform.  Medical & biological engineering & computing, 53(9):843–855, 2015\.
  * [15] Kaustubh Supekar, Vinod Menon, Daniel Rubin, Mark Musen, and Michael D Greicius.  Network analysis of intrinsic functional brain connectivity in Alzheimer’s disease.  PLoS Comput Biol, 4(6):e1000100, 2008.
  * [16] A Kabbara, H Eid, W El Falou, M Khalil, F Wendling, and M Hassan.  Reduced integration and improved segregation of functional brain networks in Alzheimer’s disease.  Journal of neural engineering, 15(2):026023, 2018.
  * [17] Fabrizio Vecchio, Francesca Miraglia, Francesca Alú, Alessandro Orticoni, Elda Judica, Maria Cotelli, and Paolo Maria Rossini.  Contribution of graph theory applied to EEG data analysis for Alzheimer’s disease versus vascular dementia diagnosis.  Journal of Alzheimer’s Disease, 82(2):871–879, 2021.
  * [18] Cornelis J Stam, BF Jones, G Nolte, M Breakspear, and Ph Scheltens.  Small-world networks and functional connectivity in Alzheimer’s disease.  Cerebral cortex, 17(1):92–99, 2007.
  * [19] Burcu Oltu, Mehmet Feyzi Akşahin, and Seda Kibaroğlu.  A novel electroencephalography based approach for alzheimer’s disease and mild cognitive impairment detection.  Biomedical Signal Processing and Control, 63:102223, 2021.
  * [20] Zhenxi Song, Bin Deng, Jiang Wang, and Ruofan Wang.  Biomarkers for alzheimer’s disease defined by a novel brain functional network measure.  IEEE Transactions on Biomedical Engineering, 66(1):41–49, 2018\.
  * [21] Haitao Yu, Xinyu Lei, Zhenxi Song, Chen Liu, and Jiang Wang.  Supervised network-based fuzzy learning of eeg signals for alzheimer’s disease identification.  IEEE Transactions on Fuzzy Systems, 28(1):60–71, 2019.
  * [22] Cosimo Ieracitano, Nadia Mammone, Alessia Bramanti, Amir Hussain, and Francesco C Morabito.  A convolutional neural network approach for classification of dementia stages based on 2d-spectral representation of eeg recordings.  Neurocomputing, 323:96–107, 2019.
  * [23] L Dharshana Deepthi, D Shanthi, and M Buvana.  An intelligent alzheimer’s disease prediction using convolutional neural network (cnn).  International Journal of Advanced Research in Engineering and Technology (IJARET), 11(4):12–22, 2020.
  * [24] Xiaojun Bi and Haibo Wang.  Early alzheimer’s disease diagnosis based on eeg spectral images using deep learning.  Neural Networks, 114:119–135, 2019.
  * [25] Cameron J Huggins, Javier Escudero, Mario A Parra, Brian Scally, Renato Anghinah, Amanda Vitória Lacerda De Araújo, Luis F Basile, and Daniel Abasolo.  Deep learning of resting-state electroencephalogram signals for three-class classification of alzheimer’s disease, mild cognitive impairment and healthy ageing.  Journal of Neural Engineering, 18(4):046087, 2021.
  * [26] Caroline L Alves, Aruane M Pineda, Kirstin Roster, Christiane Thielemann, and Francisco A Rodrigues.  Eeg functional connectivity and deep learning for automatic diagnosis of brain disorders: Alzheimer’s disease and schizophrenia.  Journal of Physics: Complexity, 3(2):025001, 2022.
  * [27] Rui Li, Yiting Wang, and Bao-Liang Lu.  A multi-domain adaptive graph convolutional network for eeg-based emotion recognition.  In Proceedings of the 29th ACM International Conference on Multimedia, pages 5565–5573, 2021.
  * [28] Shiva Asadzadeh, Tohid Yousefi Rezaii, Soosan Beheshti, and Saeed Meshgini.  Accurate emotion recognition using bayesian model based eeg sources as dynamic graph convolutional neural network nodes.  Scientific Reports, 12(1):1–14, 2022.
  * [29] Lu-Bao LI, Tian CHEN, Fu-Ji REN, and Bei-Bei LUO.  Bimodal emotion recognition method based on graph neural networks and attention.  Journal of Computer Applications, page 0, 2022.
  * [30] Khadijeh Raeisi, Mohammad Khazaei, Pierpaolo Croce, Gabriella Tamburro, Silvia Comani, and Filippo Zappasodi.  A graph convolutional neural network for the automated detection of seizures in the neonatal eeg.  Computer Methods and Programs in Biomedicine, page 106950, 2022\.
  * [31] Siyi Tang, Jared Dunnmon, Khaled Kamal Saab, Xuan Zhang, Qianying Huang, Florian Dubost, Daniel Rubin, and Christopher Lee-Messer.  Self-supervised graph neural networks for improved electroencephalographic seizure analysis.  In International Conference on Learning Representations, 2021.
  * [32] Mingyi Sun, Weigang Cui, Shuyue Yu, Hongbin Han, Bin Hu, and Yang Li.  A dual-branch dynamic graph convolution based adaptive transformer feature fusion network for eeg emotion recognition.  IEEE Transactions on Affective Computing, 2022.
  * [33] Yongqiang Yin, Xiangwei Zheng, Bin Hu, Yuang Zhang, and Xinchun Cui.  Eeg emotion recognition using fusion model of graph convolutional neural networks and lstm.  Applied Soft Computing, 100:106954, 2021.
  * [34] Jianhui Chen, Hui Qian, and Xiaoliang Gong.  Bayesian graph neural networks for eeg-based emotion recognition.  In Clinical Image-Based Procedures, Distributed and Collaborative Learning, Artificial Intelligence for Combating COVID-19 and Secure and Privacy-Preserving Machine Learning, pages 24–33. Springer, 2021\.
  * [35] Ziyu Jia, Youfang Lin, Jing Wang, Ronghao Zhou, Xiaojun Ning, Yuanlai He, and Yaoshuai Zhao.  Graphsleepnet: Adaptive spatial-temporal graph convolutional networks for sleep stage classification.  In IJCAI, pages 1324–1330, 2020.
  * [36] Hong Zeng, Qi Wu, Yanping Jin, Haohao Zheng, Mingming Li, Yue Zhao, Hua Hu, and Wanzeng Kong.  Siam-gcan: a siamese graph convolutional attention network for eeg emotion recognition.  IEEE Transactions on Instrumentation and Measurement, 2022.
  * [37] Hanjie Liu, Jinren Zhang, Qingshan Liu, and Jinde Cao.  Minimum spanning tree based graph neural network for emotion classification using eeg.  Neural Networks, 145:308–318, 2022.
  * [38] Qi Chang, Cancheng Li, Qing Tian, Qijing Bo, Jicong Zhang, Yanbing Xiong, and Chuanyue Wang.  Classification of first-episode schizophrenia, chronic schizophrenia and healthy control based on brain network of mismatch negativity by graph neural network.  IEEE Transactions on Neural Systems and Rehabilitation Engineering, 29:1784–1794, 2021.
  * [39] Xiaoyu Li, Buyue Qian, Jishang Wei, An Li, Xuan Liu, and Qinghua Zheng.  Classify eeg and reveal latent graph structure with spatio-temporal graph convolutional neural network.  In 2019 IEEE International Conference on Data Mining (ICDM), pages 389–398. IEEE, 2019.
  * [40] Peixiang Zhong, Di Wang, and Chunyan Miao.  Eeg-based emotion recognition using regularized graph neural networks.  IEEE Transactions on Affective Computing, 2020.
  * [41] Maksim Zhdanov, Saskia Steinmann, and Nico Hoffmann.  Investigating brain connectivity with graph neural networks and gnnexplainer.  arXiv preprint arXiv:2206.01930, 2022.
  * [42] Yujia Li, Richard Zemel, Marc Brockschmidt, and Daniel Tarlow.  Gated graph sequence neural networks.  In Proceedings of ICLR’16, 2016.
  * [43] Ekagra Ranjan, Soumya Sanyal, and Partha Talukdar.  Asap: Adaptive structure aware pooling for learning hierarchical graph representations.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 5470–5477, 2020.
  * [44] Tao Xu, Wang Dang, Jiabao Wang, and Yun Zhou.  Dagam: A domain adversarial graph attention model for subject independent eeg-based emotion recognition.  arXiv preprint arXiv:2202.12948, 2022.
  * [45] Andac Demir, Toshiaki Koike-Akino, Ye Wang, Masaki Haruna, and Deniz Erdogmus.  Eeg-gnn: Graph neural networks for classification of electroencephalogram (eeg) signals.  In 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pages 1061–1067. IEEE, 2021\.
  * [46] Daniel J. Blackburn, Yifan Zhao, Matteo De Marco, Simon M. Bell, Fei He, Hua-Liang Wei, Sarah Lawrence, Zoe C. Unwin, Michelle Blyth, and Jenna Angel.  A pilot study investigating a novel non-linear measure of eyes open versus eyes closed eeg synchronization in people with Alzheimer’s disease and healthy controls.  Brain sciences, 8(7):134, 2018.
  * [47] Kyunghyun Cho, Bart Van Merriënboer, Dzmitry Bahdanau, and Yoshua Bengio.  On the properties of neural machine translation: Encoder-decoder approaches.  arXiv preprint arXiv:1409.1259, 2014.
  * [48] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [49] Matthias Fey and Jan E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.

##  Appendix A Hyperparameters of proposed model

The optmised and allowed values of the various hyperparameters of the proposed
adaptive gated graph convolutional network (AGGCN) are reported in Tables 3
and 4, respectively.

$L_{CNN}$ | kernel size | CNN filters | $h_{CNN}$ | $drop_{CNN}$ | $k_{KNN}$
---|---|---|---|---|---
2 | 2 | 45 | 403 | 0 | 11
$R$ | $h_{GNN}$ | activation | aggregation | $drop_{GNN}$ |
4 | 403 | ReLU | sum | 0.04 |
$k_{pool}$ | $drop_{pool}$ | negative slope | $L_{MLP}$ | $h_{MLP}$ | $drop_{MLP}$
5 | 0.353 | 0 | 2 | 323 | 0.084
learning rate | momentum | weight decay | $\gamma$ | $\sigma$ | $p_{noise}$
0.037 | 0.192 | 0 | 0.95 | 0.346 | 0.6

Table 3: Hyper-parameter values of the optimised model Hyperparameter | Values
---|---
$L_{CNN}$ | $[1,\dots,4]$
kernel size | $[2,\dots,4]$
CNN filters | $[16,\dots,100]$
$h_{CNN}$ | $[16,\dots,1024]$
$drop_{CNN}$ | $[0,0.9]$
$k_{KNN}$ | $[1,\dots,23]$
$R$ | $[1,\dots,10]$
$h_{GNN}$ | $[16,\dots,1024]$
activation | $ReLU,Tanh,ELU,LeakyReLU$
aggregation | $add,mean,max$
$drop_{GNN}$ | $[0,0.9]$
$k_{pool}$ | $[1,\dots,23]$
$drop_{pool}$ | $[0,0.9]$
negative slope | $[0,0.5]$
$L_{MLP}$ | $[1,\dots,5]$
$h_{MLP}$ | $[16,\dots,2048]$
$drop_{MLP}$ | $[0,0.9]$
learning rate | $[0.001,0.1]$
momentum | $[0,0.9]$
weight decay | $[0,0.1]$
$\gamma$ | $[0.8,0.95]$
$\sigma$ | $[0,0.5]$
$p_{noise}$ | $[0,0.6]$
Table 4: Hyper-parameter value ranges allowed during optimisation

##  Appendix B Parameter sensitivity parameters

Multiple parameter sensitivity experiments were performed to test the
influence of the selected crucial hyperparameters of AGGCN. The results of
these experiments are reported in Figures 7, 8, 9 and 10 for the number of
GGCN iterations, K-nearest neighbour edges kept in the sparse learned graph
structure, size of the coarsened (pooled) graph and aggregation function,
respectively.

Figure 7: Sensitivity of the proposed model to the number of iterations of the
GGCN encoder. The error bars show the standard deviation of accuracies
measured across 30 repetitions. The optimal value showed in blue. Figure 8:
Sensitivity of the proposed model to the k-nearest-neighbour edges kept in the
learned graph structure. The error bars show the standard deviation of
accuracies measured across 30 repetitions. The optimal value showed in blue.
Figure 9: Sensitivity of the proposed model to the size of the pooled graph.
The error bars show the standard deviation of accuracies measured across 30
repetitions. The optimal value showed in blue. Figure 10: Sensitivity of the
proposed model to the choice of the aggregation function. The error bars show
the standard deviation of accuracies measured across 30 repetitions. The
optimal value showed in blue.


# Mathematical modelling of zika virus in Brazil

Saúl E. Buitrago Boret<EMAIL_ADDRESS>René Escalante<EMAIL_ADDRESS>Minaya Villasana<EMAIL_ADDRESS>Dpto. Cómputo Científico y Estadística,
Universidad Simón Bolívar, Caracas, Venezuela

###### Abstract

In this paper we study some deterministic mathematical models that seek to
explain the expansion of zika virus, as a viral epidemic, using published data
for Brazil. SIR type models are proposed and validated using the epidemic data
found, considering several aspects in the spread of the disease. Finally, we
confirmed that the crucial epidemic parameter such as $R_{0}$ is consistent
with those previously reported in the literature for other areas. We also
explored variations of the parameters within Brazil for different federal
entities. We concluded that a parsimonious model that includes both human and
vector populations best describe the epidemic parameters.

###### keywords:

Mathematical modelling, zika, SIR, epidemiological modelling

††journal: Journal of LaTeX Templates

##  1 Introduction

The zika virus (ZIKV) was first isolated in 1947 from a sentinel rhesus monkey
in the Zika forest in Uganda (see [9]) and was classified by sequence analysis
into two genotypes, African and Asian (see [15]). In April 2007, a large
epidemic of Asian genotype ZIKV was reported in Yap Island and Guam,
Micronesia. Between 2013-2014 the Asian genotype caused epidemics reported in
several Pacific Islands, including French Polynesia, New Caledonia, Cook
Islands, Tahiti, and Easter Island. (see [14] and [1]).

In a general review published in 2014 by Ioos et al. [18] it was reported that
the ZIKV infection caused two major epidemics in Pacific previously naive
territories, in less than a decade. This emergent arbovirosis transmitted by
mosquitoes of the Aedesgenus has a high potential for spreading in countries
where the vector is present.

In March 2015, the first endogenous transmission of Zika virus in Brazil was
reported. Beginning 2015, many patients came to health services in the state
of Rio Grande do Norte presenting dengue-like symptoms. However, physicists
and initial laboratory tests ruled out the presence of both dengue virus and
chikungunya. Some samples of these patients were transferred to the Instituto
Oswaldo Cruz (world leader in research on tropical diseases) where they
performed very sensitive screening DNA virus, the results were consistent with
those expected for the presence of Zika. The issue becomes very serious and
worrisome from May 2015, when the first case of a fetus with microcephaly
whose mother had suffered a Zika infection was detected. It may have been a
fluke, but cases of microcephaly have grown exponentially in the geographical
areas where the Zika virus is present (see [16, 3, 28]). In fact, Zika
intrauterine infection has been reported by detecting the virus in the
amniotic fluid of some pregnant patients (see [24, 30, 8]). A recently
published modelling study conducted in French Polynesia, where a documented
outbreak of Zika virus infection occurred in 2013, retrospectively estimated
the risk to be 95 cases of microcephaly per 10000 women infected during the
first trimester (see [19]). Moreover, a rare complication of Zika virus
infection can lead to neurological and autoimmune disorders such as Syndrome
Gillain-Barre. In fact there are certain press releases that summarize the
increase of this syndrome in patients who have overcome the disease (see
[35]).

In the following months ZIKV was reported in several countries of South and
Central America, and the Caribbean. Fig. 1 (Source: Washington Post; Credit:
Laris Karklis) shows graphically the propagation of zika virus from 1947 to
2015.

Recent studies in mice suggest that the virus could attack the adult brain as
well (see [2] and [25]).

Figure 1: Propagation of zika virus from Zika forest in Uganda (1947) to
Brazil (2015). Source: Washington Post; Credit: Laris Karklis.

ZIKV infection, together with denge and chikungunya, are one of the leading
causes of illness in the tropics and subtropics, where it inflicts substantial
health, economic and social burdens. Humans are infected with zika virus by
the bite of an infective female mosquito Aedes aegypti, the principal vector
of zika. Once a person gets bitten by an infective mosquito, the virus
undergoes an incubation period of about 3 to 12 days, after which the person
enters the acute phase of infection. The acute phase can be as short as 2 days
and as long as 7 days. If other female Aedes aegypti mosquitoes bites the ill
person during this acute phase, those mosquitoes may become infected and
subsequently begin the transmission cycle anew. Fig. 2 (Source: CDC, PLOS,
Reuters; Credits: David Foster, Laurie Garrett, Doug Halsey and Gabriela
Meltzer) shows graphically how zika virus enters the human population.

Figure 2: How the zika virus enters the human population. Source: CDC, PLOS,
Reuters; Credits: David Foster, Laurie Garrett, Doug Halsey and Gabriela
Meltzer.

Many people infected with zika virus will not have symptoms or will only have
mild symptoms. The most common symptoms of zika are fever, rash, joint pain,
conjuntivitis (red eyes), muscle pain and headache. Zika is usually mild with
symptoms lasting for several days to a week. People usually do not get sick
enough to go to the hospital, and they very rarely die of zika. For this
reason, many people might not realize they have been infected. Symptoms of
zika are similar to other viruses spread through mosquito bites, like dengue
and chikungunya.

There is scientific consensus that zika virus is a cause of microcephaly and
Guillain-Barré syndrome (see [35]). Links to other neurological complications
are also being investigated. Sexual transmission of zika virus is also
possible. Other modes of transmission such as blood transfusion are being
investigated.

The foundations of the entire approach to epidemiology, based on compartmental
models, were laid by public health physicians such as Sir R.A. Ross [29], W.H.
Hamer, A.G. McKendrick, and W.O. Kermack [20, 21, 22] between 1900 and 1935,
along with important contributions from a statistical perspective by J.
Brownlee [6, 12].

Mathematical models have been extensively used to study the dynamics of
infectious diseases at population level. Most continuous time models are in
the form of ordinary differential equations (ODEs). Such ODE models assume
that the population is well mixed, and the transmission is instantaneous (see
[4, 5]).

Mathematical modelling is typically the only way to examine the possible
impact of different release and control scenarios. Questions that can be
addressed are, for instance, what fraction of the population should be
quarantined and/or vaccinated? How fast can control measures to be
implemented?, etc.

The basic reproductive ratio (see [10, 11] and [17]), $R_{0}$, is defined as
the expected number of secondary infections arising from a single individual
during his or her entire infectious period, in a population of susceptibles.
$R_{0}$ often serves as a threshold parameter that predicts whether an
infection will spread. Determinig $R_{0}$ is vital to understand and
characterize the dynamics of the disease. However this crucial parameter is
model dependent. We study and calculate $R_{0}$ using different modelling
perspectives that can allow us to draw better conclusions on its validity and
range. Also, this parameter is investigated at different granularity levels:
country wide and state wide.

This paper, based on the preliminary work by Buitrago et al. [7], is organized
as follows: section 2 outlines the mathematical models used, section 3 briefly
describe details around the basic reproductive ratio $R_{0}$ for the models
developed in section 2, section 4 provides information about the data used and
in the last section, section 5, a discussion of the methodologies and their
application to the data sets are summerized.

##  2 Mathematical Modelling

We formulate our descriptions as compartmental models, with the population
under study being divided into compartments and with assumptions about the
nature and time rate of transfer from one compartment to another.

In formulating models in terms of the derivatives of the sizes of each
compartment we are also assuming that the number of members in a compartment
is a differentiable function of time. This assumption is plausible once a
disease outbreak has become established but is not valid at the beginning of a
disease outbreak when there are only a few infectives.

In this work we describe models for epidemics, acting on a sufficiently rapid
time scale that demographic effects, such as births, natural deaths, maintain
a constant level of the overoll population, and migration may be ignored.

All the models considered in this work satisfy the following assumptions:

  * 1.

There is homogeneous mixing, which means that individuals of the population
make contact at random and do not mix mostly in a smaller subgroup.

  * 2.

The disease is novel, so no vaccination is available and or applied.

  * 3.

Any recovered person has permanent immunity or least considered as such within
the time-frame of the disease.

  * 4.

The population size is constant for the models.

###  2.1 SIR model

Consider a population in which a small number of its members suffer from an
infectious disease that can be transmitted to other members of the same
population. The objective we are pursuing now is to determine what proportion
of the total population will be infected and for how long, using a
mathematical model that incorporates into their structure the transmission
mechanisms that we consider important.

In order to model such an epidemic we divide the population being studied into
three classes labeled $S$, $I$, and $R$.

Let $S(t)$ denote the number of individuals who are susceptible to the
disease, which can acquire the infection through contacts with infectious,
that is, who are not (yet) infected at time $t$. $I(t)$ denotes the number of
infected individuals, assumed infectious and able to spread the disease by
contact with susceptibles. $R(t)$ denotes the number of individuals who have
been infected and then removed from the possibility of being infected again or
of spreading infection (see Fig. 3).

Figure 3: Structure of the SIR model.

Removal is carried out through isolation from the rest of the population,
through immunization against infection, through recovery from the disease with
full immunity against reinfection, or through death caused by the disease.
These characterizations of removed members are different from an
epidemiological perspective but are often equivalent from a modelling point of
view that takes into account only the state of an individual with respect to
the disease.

We will use the terminology SIR to describe a disease that confers immunity
against reinfection, to indicate that the passage of individuals is from the
susceptible class $S$ to the infective class $I$ to the removed class $R$. The
mathemaical model is:

$\displaystyle\frac{dS}{dt}(t)$ $\displaystyle=$ $\displaystyle\mu N-\beta
S(t)I(t)/N-\mu S(t),$ (1) $\displaystyle\frac{dI}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta S(t)I(t)/N-\gamma I(t)-\mu I(t),$ (2)
$\displaystyle\frac{dR}{dt}(t)$ $\displaystyle=$ $\displaystyle\gamma I(t)-\mu
R(t),$ (3) $\displaystyle\frac{dC_{i}}{dt}$ $\displaystyle=$ $\displaystyle
p\beta S(t)I(t),$ (4)

with the initial conditions in $t=t_{0}$

$\displaystyle S(t_{0})$ $\displaystyle=$ $\displaystyle S_{0}>0,$
$\displaystyle I(t_{0})$ $\displaystyle=$ $\displaystyle I_{0}>0,$
$\displaystyle R(t_{0})$ $\displaystyle=$ $\displaystyle R_{0}>0,$

where $S_{0}$, $I_{0}$ and $R_{0}$ are, respectively, the initial number of
suceptible, infected and recovered people with $\beta$, $\gamma$ and $\mu$
positive constants. As is often the case, not all infectives are symptomatic,
specially in the Zika virus, and thus not all cases are reported as such
making the determination of the real number of infectives a difficult task. In
equation 4, $C_{i}$ accounts for the cumulative infectives and is a smooth
monotone function that is used for identification purposes. The parameter $p$
is a proportion of the infectives that are reported. $\beta$ is the
transmission rate from mosquitoes to humans, $\gamma$ is the per capita rate
of recovery in humans such that $1/\gamma$ is the mean infectious period for
humans, $\mu$ is the per capita rate of mortality in humans such that $1/\mu$
is the life expectancy of humans, and $N$ is the human population size.

###  2.2 SIR/SI model

This model is an extension of the SIR model, and has been used before in the
study of the dynamics of dengue in Thailand by Pandey et al. [27].

The population is divided into three classes for humans and two clases for
mosquitoes or vectors that transmit the disease. $S_{H}$ represents the number
of susceptible, $I_{H}$ the number of infectious, and $R_{H}$ the number
recovered individuals in the human sub-population. Similarly, $S_{v}$
represents the proportion of mosquitoes currently susceptible, and $I_{v}$ the
proportion of infectious mosquitoes infectious (see Fig. 4).

Figure 4: Structure of the SIR/SI model.

Mosquitoes are assumed to remain infectious for life. $\beta_{v}$ is the
transmission rate from humans to mosquitoes, $\beta_{H}$ is the transmission
rate from mosquitoes to humans, while $1/\gamma_{H}$ and $1/\mu_{H}$ are the
mean infectious period and the mean lifespan of humans, $1/\mu_{v}$ is the
mean lifespan of mosquitoes. $N_{H}$ stands for the human population size and
$N_{v}$ is the mosquito population size. The mathematical model is written as
follows:

$\displaystyle\frac{dS_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mu_{H}N_{H}-\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\mu_{H}S_{H}(t),$
(5) $\displaystyle\frac{dI_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\gamma_{H}I_{H}(t)-\mu_{H}I_{H}(t),$
(6) $\displaystyle\frac{dR_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\gamma_{H}I_{H}(t)-\mu_{H}R_{H}(t)$ (7)
$\displaystyle\frac{dI_{v}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{v}S_{v}(t)I_{H}(t)/N_{H}-\mu_{v}I_{v}(t),$ (8)
$\displaystyle\frac{dS_{v}}{dt}(t)$ $\displaystyle=$
$\displaystyle-\beta_{v}S_{v}(t)I_{H}(t)/N_{H}+\mu_{v}I_{v}(t),$ (9)
$\displaystyle\frac{dC_{i}}{dt}$ $\displaystyle=$ $\displaystyle
p\beta_{H}S_{H}(t)I_{v}(t)/N_{v},$ (10)

with the initial conditions in $t=t_{0}$

$\displaystyle I_{H}(t_{0})$ $\displaystyle=$ $\displaystyle I_{H0}>0,$
$\displaystyle R_{H}(t_{0})$ $\displaystyle=$ $\displaystyle R_{H0}>0,$
$\displaystyle I_{v}(t_{0})$ $\displaystyle=$ $\displaystyle I_{v0}>0,$

where $I_{H0}$, $R_{H0}$ and $I_{v0}$ are respectively, the initial number of
infected people, the initial number of recovered people, and the initial
number of infectious mosquitoes, respectively. As in the previous model, we
use the cummulative number of infectives for identification purposes (equation
10).

Given the fact that the human population remains constant, one can express
$R_{H}$ in terms of the variables $S_{H}$ and $I_{H}$, i.e.
$R_{H}=N_{H}-S_{H}-I_{H}$, therefore equation 7 can be discarded and we can
reduce the dimensionality of the system. A similar argument is true for the
case of $S_{v}$ in equation 9, i.e. $S_{v}=N_{v}-I_{v}$.

###  2.3 SEIR/SEI model

This model in based upon the work of Kucharski et al. (see [23]). This model
incorporates a new compartment, exposed, for the human and the mosquito
subpopulation which represents the number of individuals (and mosquitoes) that
are incubating the virus, $E_{H}$ (and $E_{v}$ in the case of the vector
population), i.e. where individuals (mosquitoes) are infected but are not able
yet to transmit the virus. The inclusion of such a compartment into the model
is due to the fact that it is known that vector and human populations incubate
the virus for a number of days. Fig. 5 depicts the structure of this new
model.

Figure 5: Structure of the SEIR/SEI model.

All parameters have the same connotation as in the SIR/SI model, and here
$1/\kappa_{H}$ and $1/\kappa_{v}$ are the mean latent periods for humans and
mosquitoes population respectively. The governing equations of the model
follow:

$\displaystyle\frac{dS_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mu_{H}N_{H}-\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\mu_{H}S_{H}(t),$
(11) $\displaystyle\frac{dE_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\kappa_{H}E_{H}(t)-\mu_{H}E_{H}(t),$
(12) $\displaystyle\frac{dI_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\kappa_{H}E_{H}(t)-\gamma_{H}I_{H}(t)-\mu_{H}I_{H}(t),$ (13)
$\displaystyle R_{H}$ $\displaystyle=$ $\displaystyle
N_{H}-S_{H}-I_{H}-E_{H},$ (14) $\displaystyle\frac{dE_{v}}{dt}(t)$
$\displaystyle=$
$\displaystyle\beta_{v}(S_{v}+E_{v})I_{H}(t)/N_{H}-\kappa_{v}E_{v}(t)-\mu_{v}E_{v}(t),$
(15) $\displaystyle\frac{dI_{v}}{dt}(t)$ $\displaystyle=$
$\displaystyle\kappa_{v}E_{v}(t)-\mu_{v}I_{v}(t),$ (16) $\displaystyle S_{v}$
$\displaystyle=$ $\displaystyle N_{v}-I_{v}-E_{v},$ (17)
$\displaystyle\frac{dC_{i}}{dt}$ $\displaystyle=$ $\displaystyle
p\beta_{H}S_{H}(t)I_{v}(t)/N_{v},$ (18)

with the initial conditions in $t=t_{0}$

$\displaystyle I_{H}(t_{0})$ $\displaystyle=$ $\displaystyle I_{H0}>0,$
$\displaystyle R_{H}(t_{0})$ $\displaystyle=$ $\displaystyle R_{H0}>0,$
$\displaystyle I_{v}(t_{0})$ $\displaystyle=$ $\displaystyle I_{v0}>0,$
$\displaystyle E_{H}(t_{0})$ $\displaystyle=$ $\displaystyle E_{H0}>0,$
$\displaystyle E_{v}(t_{0})$ $\displaystyle=$ $\displaystyle E_{v0}>0,$

where $I_{H0}$, $R_{H0}$, $I_{v0}$, $E_{H0}$ and $E_{v0}$ are respectively,
the initial number of infected people, the initial number of recovered people,
the initial number of infectious mosquitoes, the initial number of people
incubating the virus, and the initial number of mosquitoes incubating the
virus, respectively. In equation 18, $C_{i}$ accounts for the cumulative
infectives and is a smooth monotone function that is used for identification
purposes.

##  3 The basic reproductive ratio $R_{0}$

The basic reproductive ratio (see [10, 11] and [17]), $R_{0}$, is defined as
the expected number of secondary infections arising from a single individual
during his or her entire infectious period, in a population of susceptibles.
$R_{0}$ often serves as a threshold parameter that predicts whether an
infection will spread. Determinig $R_{0}$ is vital to understand and
characterize the dynamics of the disease. However this crucial parameter is
model dependent.

$R_{0}$ is the dominant eigenvalue of the so call “next generation matrix”. It
is shown that, if $R_{0}<1$, then the disease free equilibrium is locally
asymptotically stable; whereas if $R_{0}>1$, then it is unstable (see [32,
33]).

The basic reproduction number for the SIR model is known to be calculated as

$R_{0}=\frac{\beta}{\mu+\gamma}.$

The basic reproduction number for the SIR/SI model is known to be calculated
(see [27]) as:

$R_{0}=\frac{\beta_{H}\beta_{v}}{\mu_{v}(\mu_{H}+\gamma_{H})}.$

The basic reproduction number for the SEIR/SEI model is known to be calculated
(the dominant eigenvalue of the next generation matrix, see [10, 11]) as:

$R_{0}=\frac{\beta_{H}\beta_{v}\kappa_{v}}{\mu_{v}\gamma_{v}(\kappa_{v}+\mu_{v})}.$

We study and calculate $R_{0}$ using different modelling perspectives that can
allow us to draw better conclusions on its validity and range. Also, this
parameter is investigated at different granularity levels: country wide and
state wide.

##  4 Results

The used data used to validate the models, available on the internet, was
published by Faria et al. [14, 13] while gathering information for their
research. The data is partitioned by municipalities in Brazil, thus one may
have different levels of granularity in the visualization of the data. By
aggregation, of cases in different federal entities one can determine those
entities that have the most reported cases. Table 1 shows the total number of
cases throughout 2015 as reported in the epidemiological data provided.

In table 1 only those entities that report more than 100 cases in a year are
displayed. One can see that the federal entity that reported the most number
of cases is Bahía, followed by far by Alagoas, Ceará, and Rio Grande du Norte.
All of these entities are located in the northeastern part of the country. The
northeartern region is characterized by high temperatures (annual averages
between 20 and 28 oC (68.0 and 82.4 oF), maxima of around 40 oC (104 oF)).
During the months of June and July temperatures vary between 12 and 16 oC
(53.6 and 60.8 oF) in the coastal regions, where most cases are reported.

Table 1: Reported cases in some entities of Brazil Entities | Reported cases
---|---
Bahía (BA) | 27290
Alagoas (AL) | 497
Ceará (CE) | 416
Rio Grande do Norte (RN) | 264
Pará (PA) | 155
Espíritu Santo (ES) | 125
Total Brazil | 29639

The incidence data (left panels, images (a) to (d)) and the cummulative data
used for identification purposes (right panels, images (e) to (h)) for the
agreggated country data and most reported federal entities are shown in Fig.
6.

Relevant demographic information taken from the web page of the Instituto
Brasileiro de Geografia e Estatística (IBGE) (2013) are given in table 2,
where proportion refers to the proportion of reported cases to total
population.

Table 2: Relevant demographic data Parameter | Brazil | Bahía | Alagoas | Ceará
---|---|---|---|---
Population (MM) | 207 | 15.15 | 8.86 | 3.32
Life expectancy (years) | 75 | 71.9 | 69.2 | 72.4
Reported cases | 29639 | 27290 | 416 | 497
Proportion | 1.43e-4 | 1.80e-3 | 4.69e-5 | 1.491e-4

The life expectancy was used in all parameter estimations, thus reducing the
number of free parameters to be identified.

A brief description on how to solve the system of ODEs associated to any of
the three models proposed (SIR, SIR/SI and SEIR/SEI) is given in the appendix.
The numerical models, which allow us to calculate the cumulative number of
cases for each model, were implemented in language M of MATLAB. These
functions will be used in the identification process of the unknown
parameters, being 3, 6 and 10 the number of parameters for the SIR, SIR/SI and
SEIR/SEI models respectively (see table 3).

It is important to note that the identification was performed over
nondimensionalized systems.

The problem to be solved for the estimation of the unknown parameters is the
following: Given a function $f(t)$ which represents the cumulative weekly
number of cases for a period of 49 weeks during 2015, find the parameters
$x_{i}$, $i=1,\dots,n$ such that the answer $\tilde{f}(t)$ given by the model
and the cumulative incidence $f(t)$ are closed, that is

$g(x_{min})=\mathop{\rm global\;min}_{x\in\Omega}g(x),$

where $g:\Omega\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined by

$g(t)=\int_{0}^{T}(f(t)-\tilde{f}(t))^{2}dt$

the mean squared error (MSE), $\Omega=\prod_{i=1}^{n}[a_{i},b_{i}]$ and $T$ is
the time corresponding to the data.

It is important to point out that it is possible to find more than one set of
parameters which satisfy the minimality condition required and related to the
size of $\Omega$. Also relevant is that every time the objective function
needs to be evaluated, one run of the model has to be carried out.

The identification was carried through an exhaustive search procedure within
the range of the different parameters involved in the model, not being the
best available methodology because of the high number of evaluations of the
objective function.

We considered feasible ranges for the parameters of the model whenever
possible. Thus, according to epidemiological data, previously reported results
[23, 31] and some sensibility carried out, the range considered for the
parameters are in table 3.

Table 3: Ranges for the estimation of the models unknown parameters. Model | SIR | SIR/SI | SEIR/SEI | reference for
---|---|---|---|---
Parameter |  |  |  | the reciprocal
$\beta_{H}$ (day-1) | [0.01,0.6] | [0.01,0.6] | [0.1,1.2] | -
$\gamma_{H}$ (day-1) | [0.1,0.4] | [0.1,0.5] | [0.1,0.5] | 2-9 days
$\beta_{v}$ (day-1) | - | [0.01,0.4] | [0.01,0.3] | -
$\mu_{v}$ (day-1) | - | [0.03,0.3] | [0.03,0.3] | 4-30 days
$\kappa_{H}$ (day-1) | - | - | [0.083,0.34] | 3-12 days
$\kappa_{v}$ (day-1) | - | - | [0.20,0.50] | 2-5 days
$I_{v0}$ | - | [1.e-9,5.e-5] | [8.e-7,5.e-5] | -
$E_{0}$ | - | - | [8.e-7,9.e-6] | -
$E_{v0}$ | - | - | [8.e-7,9.e-6] | -
$p$ | [1.e-5,6.e-3] | [1.e-5,6.e-3] | [4.e-5,6.e-3] | -

We estimated the total of three unknown parameters for the SIR model (see
table 4), six parameters for the SIR/SI model (see table 5), and ten
parameters for the SEIR/SI model (see table 6), using the cumulative incidence
suspicious data. Cumulative incidence is generally smoother than the original
incidence data and thus easier to fit.

The numerical solutions of the models were performed using Matlab and compared
to the cummulative data to obtain a mean squared error (MSE). Those parameters
that produced the smallest errors are reported as the best identification for
the data. In all cases the basic reproduction number is reported as well as
the MSE.

Table 4: Estimated parameters for the SIR model. Parameter | Brazil | Bahía | Alagoas | Ceará
---|---|---|---|---
$\beta_{H}$ | 0.4340 | 0.1954 | 0.26767 | 0.46124
$\gamma_{H}$ | 0.3816 | 0.1493 | 0.21576 | 0.40667
$p$ | 5.7931e-4 | 4.0204e-3 | 4.0912e-4 | 2.0506e-4
$R_{0}$ | 1.137 | 1.274 | 1.2382 | 1.1336
Error | 2.2406e-3 | 9.4618e-4 | 4.1092e-3 | 2.4544e-3
Table 5: Estimated parameters for the SIR/SI model. Parameter | Brazil | Bahía | Alagoas | Ceará
---|---|---|---|---
$\beta_{H}$ | 0.49368 | 0.48947 | 0.31889 | 0.54211
$\gamma_{H}$ | 0.21263 | 0.11 | 0.36333 | 0.31053
$\beta_{v}$ | 0.038947 | 0.037474 | 0.36000 | 0.15474
$\mu_{v}$ | 0.0333 | 0.05 | 0.22 | 0.15667
$I_{v0}$ | 3.4526e-5 | 1e-7 | 3e-8 | 1e-9
$p$ | 1.5263e-4 | 1.6632e-3 | 2.7778e-4 | 6.8737e-5
$R_{0}$ | 2.715 | 3.3338 | 1.436 | 1.724
Error | 7.0798e-4 | 6.0564e-3 | 2.9883e-3 | 7.6298e-4

Fig. 6 overlays the cummulative reported data points with the simulated
cummulative cases with the best set of parameters for each entity (right
panels) and the 3 models investigated. Although incidence data was not used in
the identification purposes, left panel on the left show an overlay of the
incidence data with the simulated incidence.

(a) (e)

(b) (f)

(c) (g)

(d) (h)

Figure 6: Incidence (images (a) to (d)) and cumulative (images (e) to (h)) data and model results for Brazil ((a) and (e)), Bahía ((b) and (f)), Alagoas ((c) and (g)) and Ceará ((d) and (h)). Points marked with circles in (a) to (d) correspond to incidence data and in (e) to (h) correspond to cumulative data. Dashed, solid, and dashdotted lines represent SIR, SIR/SI and SEIR/SEI models respectively. Table 6: Estimated parameters for the SEIR/SEI model. Parameter | Brazil | Bahía | Alagoas | Ceará
---|---|---|---|---
$\beta_{H}$ | 0.86372 | 0.84 | 0.8 | 1.1
$\gamma_{H}$ | 0.28333 | 0.19667 | 0.11 | 0.305
$\beta_{v}$ | 0.045 | 0.04 | 0.040556 | 0.06
$\mu_{v}$ | 0.03333 | 0.03333 | 0.03333 | 0.03333
$\kappa_{H}$ | 0.083 | 0.083 | 0.083 | 0.083
$\kappa_{v}$ | 0.425 | 0.2 | 0.2 | 0.5
$I_{v0}$ | 5.e-5 | 5.e-5 | 8.e-7 | 1.e-8
$E_{0}$ | 9.e-6 | 9.e-6 | 8.e-7 | 6.e-7
$E_{v0}$ | 4.9e-6 | 6.95e-6 | 8e-7 | 6.e-7
$p$ | 0.00014 | 0.00173 | 0.00014 | 0.000046667
$R_{0}$ | 3.8165 | 4.3936 | 7.5854 | 6.0867
Error | 0.00081171 | 0.00089141 | 0.0015467 | 0.00052088

##  5 Discussion

In this work we identified three models for the Zika virus in Brasil. The
models were identified using the aggregated data for the country and then for
3 federal entities, namely Bahia, Alagoas and Ceará, which were the areas with
the greatest number of cases reported. One aspect worth mentioning is that the
mayority of the cases came from Bahía. One can see from the results in tables
5, and 6 that the model parameters are similar for the aggregated data
(Brazil) and the federal region with most numbered cases. This is so for those
models that include the vector within the models, namely the SIR/SI and
SEIR/SEI. However, this is not maintained in the more basic SIR model (see
Table 4).

For the SIR model, the results for Ceará and the aggregated Brazil cases are
similar according to the estimated parameters, while Bahía and Alagoas seemed
tuned into one another. This result is in agreement with the proportions
calculated in Table 2 for Brazil and Ceará, but they do not agree when
comparing Bahía and Alagoas since there is a difference in two orders of
magnitude. But the models adjust best to similar contact and infectious rates.
It is noteworthy however that the $R_{0}$ is within similar ranges (1.1 -
1.3), see [27] with $R_{0}=1.10$ for SIR model, but far from the ranges
reported in other Zika outbreaks, see [26] ($R_{0}=4.3-5.8$ for the Yap Island
epidemic and $R_{0}=1.8-2.2$ for the French Polynesia epidemic), and [34]
($R_{0}=1.57$ for the SIR model, $R_{0}=1.65$ using a SEIR model, and
$R_{0}=1.66$ for a more complicated model, all applied to influenza A). One
reason could be that the $R_{0}$ is model dependent and in this case this very
simple model does not include important infectious factor such as the vector
population. Therefore, we regard this model to be non-informative and too
simplistic for the situation being modelled, and it was considered for
comparison purposes and as a parsimonious model. When the vector population is
included the ranges of the values for $R_{0}$ expand and the values are
bigger. For the SIR/SI model the range is (1.4 - 3.3) and for the SEIR/SEI
model the range is (3.8 - 7.5), this could be due to model dependency of
$R_{0}$. Although the ranges for the $R_{0}$ are disjunct, these are
consistent with other reported values in the literature, see [31] with
$R_{0}=4.4$ with $95\%$ CI $[3.0,6.2]$ and a one standard deviation
uncertainty of $0.9$ for the outbreak of ZIKV that began in 2015 in
Barranquilla, Colombia using the SEIR/SEI model, [27] with $R_{0}=1.57$ using
a vector-host model, and [23] with $R_{0}$ ranged for from $2.6$ (95Marquises
to $4.8$ (95% CI: 3.2 8.4) in Moorea for French Polynesia ZIKV outbreak using
the SEIR/SEI model.

Even though the models that included the vector population have more free
parameters to adjust and identify, allowing for greater degrees of freedom,
the errors in the adjustments are similar for all three models used.
Furthermore, all three models follow the data very well as one can see on the
(e) to (h) panels in Fig. 6. Even though the incidence data was not used
during the identification process, panels (a) to (d) from Figure Fig. 6
overlays the data with the induced incidence curve with the model parameters
found. In the case of Brazil, the SIR model was able to best capture the rise
in the incidence around the peak of infectiousness (day 200), but for the
subregions the models that incorporated the vector population best captured
such rise.

One point to mention is that the Brazil data aggregates information from
varying regions with varying forms and times of outbreak, which makes the
interpretation of the results difficult.

## Acknowledgements

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

## References

  * [1] B. Berkowitz, L. Karklis, S. Tan, D. Lu, P. Clark, What you need to know about the zika virus, The Washington Post, https://www.washingtonpost.com/graphics/health/zika-virus/ (2016).
  * [2] A. Bernardo, Sospechan que el zika podría afectar al cerebro de adultos, Hipertextual, https://hipertextual.com/2016/08/zika-cerebro-adultos (8 2016).
  * [3] P. Brasil, J. Pereira, C. Raja Gabaglia, L. Damasceno, M. Wakimoto, R. Ribeiro Nogueira, P. Carvalho de Sequeira, A. Machado Siqueira, L. Abreu de Carvalho, D. Cotrim da Cunha, G. Calvet, E. Neves, M. Moreira, A. Rodrigues-Bai/ ao, P. Nassar de Carvalho, C. Janzen, S. Valderramos, J. Cherry, A. Bispo de Filippis, K. Nielsen Saines, Zika virus infection in pregnant women in rio de janeiro - preliminary report, The New England Journal of Medicine, doi: 10.1056/NEJMoa1602412 (3 2016).
  * [4] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, 2nd ed., Springer, 2012, ISBN 978-1-4614-1685-2.
  * [5] F. Brauer, C. Castillo-Chavez, E. De La Pava-Salgado, K. Barley, C. Castillo-Garsow, D. Chowell, B. Espinoza, P. Gonz lez Parra, C. Hern ndez Su rez, V. Moreno, Modelos de la propagación de enfermedades infecciosas, Universidad Aut noma de Occidente, Cali, Colombia, 2015, ISBN 978-958-8713-65-6.
  * [6] J. Brownlee, Certain considerations regarding the epidemiology of phthisis pulmonalis, Public Health 29 (1916) 130–145, doi: 10.1016/S0033-3506(15)80528-9.
  * [7] S. Buitrago, R. Escalante, M. Villazana, Modelización matemática da la expansión del virus zika en la región Latinomericana, in: Simulación y Aplicaciones Recientes para Ciencia y Tecnología, Ed: Y. G nzales, E. D vila, V. Duarte, M.V. Candal, O. Pelliccioni, J. Darias, M. Cerrolaza, 2016, ISBN 978-980-7161-05-3.
  * [8] G. Calvet, R. Aguiar, A. Melo, S. Sampaio, I. de Filippis, A. Fabri, E. Araujo, P. de Sequeira, M. de Mendon a, L. de Oliveira, D. Tschoeke, C. Schrago, F. Thompson, P. Brasil, F. dos Santos, R. Nogueira, A. Tanuri, A. de Filippis, Detection and sequencing of zika virus from amniotic fluid of fetuses with microcephaly in brazil: a case study, The Lancet Infectious Diseases 16 (6) (2016) 653–660, doi: 10.1016/S1473-3099(16)00095-5.
  * [9] G. Dick, S. Kitchen, A. Haddow, Zika virus (i). isolations and serological specificity, Trans. R. Soc. Trop. Med. Hyg. 46 (1952) 509–520, doi: 10.1016/0035-9203(52)90042-4.
  * [10] O. Diekmann, J. Heesterbeek, J. Matz, On the definition and computation of the basic reproduction ratio ${R}_{0}$ in models for infectious diseases in heterogeneous populations, Journal of Mathematical Biology 28 (4) (1990) 365–382, doi: 10.1007/BF00178324.
  * [11] O. Diekmann, J. Heesterbeek, M. Roberts, The construction of next-generation matrices for compartmental epidemic models, Journal of the Royal Society Interface 7 (47) (2010) 873–885, doi: 10.1098/rsif.2009.0386.
  * [12] V. Farewell, T. Johnson, Commentary: Dr john brownlee ma, md, dsc, dph (cantab), frfps, fss, frmets (1868 1927), public health officer, geneticist, epidemiologist and medical statistician, International Journal of Epidemiology 42 (2013) 935–943, doi: 10.1093/ije/dyt067.
  * [13] N. Faria, R. Azevedo, M. Kraemer, et al., Data from: Zika virus in the americas: early epidemiological and genetic findings, Dryad Digital Repository, doi: 10.5061/dryad.6kn23 (2016).
  * [14] N. Faria, R. Azevedo, M. Kraemer, et al., Zika virus in the americas: early epidemiological and genetic findings, Science 352 (6283) (2016) 345–349, doi: 10.1126/science.aaf5036.
  * [15] O. Faye, O. Faye, D. Diallo, M. Diallo, M. Weidmann, A. Sall, Quantitative real-time pcr detection of zika virus and evaluation with field-caught mosquitoes, Virol. J. 10 (311) (2013) 1–8, doi: 10.1186/1743-422X-10-311.
  * [16] V. Greenwood, New insights into how zika harms the brain, Quanta Magazine, https://www.quantamagazine.org/20160707-how-zika-interferes-with-brain-development/ (2016).
  * [17] J. Heffernan, R. Smith, L. Wahl, Perspectives on the basic reproductive ratio, Journal of The Royal Society Interface, doi: 10.1098/rsif.2005.0042 (2005).
  * [18] S. Ioos, H. Mallet, I. Leparc-Goffart, V. Gauthier, T. Cardoso, M. Herida, Current zika virus epidemiology and recent epidemics, Med. Mal. Infect. 44 (7) (2014) 302–307, doi: 10.1016/j.medmal.2014.04.008.
  * [19] M. Johansson, L. Mier y Teran Romero, J. Reefhuis, S. Gilboa, S. Hills, Zika and the risk of microcephaly, New England Journal of Medicine 375 (2016) 1–4, doi: 10.1056/NEJMp1605367.
  * [20] W. Kermack, A. McKendrick, A contribution to the mathematical theory of epidemics, Proc. Royal Soc. London 115 (1927) 700–721, doi: 10.1098/rspa.1927.0118.
  * [21] W. Kermack, A. McKendrick, A contribution to the mathematical theory of epidemics, part ii, Proc. Royal Soc. London 138 (1932) 55–83, http://www.math.utah.edu/ bkohler/Journalclub/kermack1932.pdf.
  * [22] W. Kermack, A. McKendrick, A contribution to the mathematical theory of epidemics, part iii, Proc. Royal Soc. London 141 (1933) 94–112, http://www.uvm.edu/pdodds/files/papers/others/1933/kermack1933.pdf.
  * [23] A. Kucharski, S. Funk, R. Eggo, H. Mallet, W. Edmunds, E. Nilles, Transmission dynamics of zika virus in island populations: A modelling analysis of the 2013- 14 french polynesia outbreak, PLoS Negl Trop Dis 10 (5) (2016) 1–15, doi: 10.1371/journal.pntd.0004726.
  * [24] C. Li, D. Xu, Q. Ye, S. Hong, Y. Jiang, X. Liu, N. Zhang, L. Shi, C.-F. Qin, Z. Xu, Zika virus disrupts neural progenitor development and leads to microcephaly in mice, Cell Stem Cell 19 (1) (2016) 120–126, doi: 10.1016/j.stem.2016.04.017.
  * [25] H. Li, L. Saucedo Cuevas, J. Regla Nava, A. Terskikh, S. Shresta, J. Gleeson, Zika virus infects neural progenitors in the adult mouse brain and alters proliferation, Cell Stem Cell 19 (2016) 1–6, doi: 10.1016/j.stem.2016.08.005.
  * [26] H. Nishiura, R. Kinoshita, K. Mizumoto, Y. Y., K. Nah, Transmission potential of zika virus infection in the south pacific, International Journal of Infectious Diseases 45 (2016) 95–97, doi: 10.1016/j.ijid.2016.02.017.
  * [27] A. Pandey, A. Mubayi, J. Medlock, Comparing vector-host and sir models for dengue transmission, Mathematical Biosciences 246 (2) (2013) 252–259, doi: 10.1016/j.mbs.2013.10.007.
  * [28] H. Retallack, E. Di-Lullo, C. Arias, K. Knopp, C. Sandoval Espinosa, M. Laurie, Y. Zhou, M. Gormley, W. Mancia Leon, R. Krencik, E. Ullian, J. Spatazza, A. Pollen, K. Ona, T. Nowakowski, J. DeRisi, S. Fisher, A. Kriegstein, Zika virus in the human placenta and developing brain: Cell tropism and drug inhibition, bioRxiv, The preprint server for biology, doi: 10.1101/058883 (2016).
  * [29] R. Ross, The Prevention of Malaria, 2nd ed., John Murray, London, 1911.
  * [30] H. Tang, C. Hammack, S. Ogden, Z. Wen, X. Qian, Y. Li, B. Yao, J. Shin, F. Zhang, E. Lee, K. Christian, R. Didier, P. Jin, H. Song, G. Ming, Zika virus infects human cortical neural progenitors and attenuates their growth, Cell Stem Cell 18 (2016) 587–590, doi: 10.1016/j.stem.2016.02.016.
  * [31] S. Towers, F. Brauer, C. Castillo-Chavez, A. Falconar, A. Mubayi, C. Romero-Vivas, Estimation of the reproduction number of the 2015 zika virus outbreak in barranquilla, colombia, and a first estimate of the relative role of sexual transmission, ResearchGate, https://www.researchgate.net/publication/303821853 (2016).
  * [32] P. van den Driessche, W. J., Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmision, Mathematical Biosciences 180 (1-2) (2002) 29–48, doi: 10.1016/S0025-5564(02)00108-6.
  * [33] P. van den Driessche, W. J., Further notes on the basic reproduction number, in: Mathematical Epidemiology, vol. 1945 of Lecture Notes in Mathematics, chap. 6, Springer, Berlin Heidelberg, 2008, pp. 159–178, doi: 10.1007/978-3-540-78911-6_6.
  * [34] J. Wallinga, M. Lipsitch, How generation intervals shape the relationship between growth rates and reproductive numbers, Proceeding of the Royal Society B 274 (2007) 599–604, doi: 10.1098/rspb.2006.3754.
  * [35] WHO, Zika strategic response plan, revised for july 2016 december 2017, Tech. rep., World Health Organization, http://reliefweb.int/sites/reliefweb.int/files/resources/WHO-ZIKV-SRF-16.3-eng.pdf (june 2016).

## Appendix - Discretizing the models

In order to exemplify the discretization methodology, the SIR/SI model
(equations 5 to 10) will be used. Given the fact that the human population
remains constant, one can express $R_{H}$ in terms of the variables $S_{H}$
and $I_{H}$, i.e. $R_{H}=N_{H}-S_{H}-I_{H}$, therefore, equation 7 for this
model can be discarded and we can reduce the dimensionality of the system. A
similar argument is true for the case of $S_{v}$ in equation 9, i.e.
$S_{v}=N_{v}-I_{v}$. Finally, the following system of ordinary differential
equations with initial conditions has to be solved.

$\displaystyle\frac{dS_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mu_{H}N_{H}-\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\mu_{H}S_{H}(t),$
(19) $\displaystyle\frac{dI_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{H}S_{H}(t)I_{v}(t)/N_{v}-\gamma_{H}I_{H}(t)-\mu_{H}I_{H}(t),$
(20) $\displaystyle R_{H}$ $\displaystyle=$ $\displaystyle N_{H}-S_{H}-I_{H}$
(21) $\displaystyle\frac{dI_{v}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{v}S_{v}(t)I_{H}(t)/N_{H}-\mu_{v}I_{v}(t),$ (22)
$\displaystyle S_{v}$ $\displaystyle=$ $\displaystyle N_{v}-I_{v},$ (23)
$\displaystyle\frac{dC}{dt}$ $\displaystyle=$ $\displaystyle
p\beta_{H}S_{H}(t)I_{v}(t)/N_{v},$ (24)

with the initial conditions in $t=t_{0}$

$\displaystyle I_{H}(t_{0})$ $\displaystyle=$ $\displaystyle I_{H0}>0,$
$\displaystyle R_{H}(t_{0})$ $\displaystyle=$ $\displaystyle R_{H0}>0,$
$\displaystyle I_{v}(t_{0})$ $\displaystyle=$ $\displaystyle I_{v0}>0.$

The following step is the normalization of this system, i.e. with

$\displaystyle s_{H}$ $\displaystyle=S_{H}/N_{H},$ $\displaystyle i_{h}$
$\displaystyle=I_{H}/N_{H},$ $\displaystyle r_{H}$
$\displaystyle=R_{H}/N_{H},$ $\displaystyle i_{v}$
$\displaystyle=I_{v}/N_{v},$ $\displaystyle s_{v}$
$\displaystyle=S_{v}/N_{v},$ $\displaystyle c$ $\displaystyle=C/N_{H},$

the following equivalent system of ODEs arises

$\displaystyle\frac{ds_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\mu_{H}-\beta_{H}s_{H}(t)i_{v}(t)-\mu_{H}s_{H}(t),$ (25)
$\displaystyle\frac{di_{H}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{H}s_{H}(t)i_{v}(t)-\gamma_{H}i_{H}(t)-\mu_{H}i_{H}(t),$
(26) $\displaystyle r_{H}$ $\displaystyle=$ $\displaystyle 1-s_{H}-i_{H}$ (27)
$\displaystyle\frac{di_{v}}{dt}(t)$ $\displaystyle=$
$\displaystyle\beta_{v}s_{v}(t)i_{H}(t)-\mu_{v}i_{v}(t),$ (28) $\displaystyle
s_{v}$ $\displaystyle=$ $\displaystyle 1-i_{v},$ (29)
$\displaystyle\frac{dc}{dt}$ $\displaystyle=$ $\displaystyle
p\beta_{H}s_{H}(t)i_{v}(t).$ (30)

The basic idea of any approximation method is to replace the original problem
by another problem that is easier to solve and whose solution is, is in some
sense, close to the solution of the original problem.

Given $M\in\mathbb{N}$, let $\\{t_{i}\\}_{1\leq i\leq M}$ be an uniform
subdivision of the time domain, with mesh length $h=\Delta t=t_{i+1}-t_{i}$
along the direction $t$.

The following is the forward finite difference for the first order operator

$\left(\frac{df}{dt}\right)_{i+1}=\frac{f(t_{i}+h)-f(t_{i})}{h}\ \ {\rm for}\
1\leq i\leq M-1$

where $f$ represents the functions $s_{H}$, $i_{H}$, $r_{H}$, $i_{v}$, $s_{v}$
and $c$.

Finally the following set of linear equations arises upon susbtitution of the
forward finite difference operator

$\displaystyle s_{H}(t_{i}+h)$ $\displaystyle=$ $\displaystyle
s_{H}(t_{i})+\mu_{H}\,h-\beta_{H}s_{H}(t_{i})i_{v}(t_{i})\,h-\mu_{H}s_{H}(t_{i})\,h,$
(31) $\displaystyle i_{H}(t_{i}+h)$ $\displaystyle=$ $\displaystyle
i_{H}(t_{i})+\beta_{H}s_{H}(t_{i})i_{v}(t_{i})\,h-\gamma_{H}i_{H}(t_{i})\,h-\mu_{H}i_{H}(t_{i})\,h,$
(32) $\displaystyle r_{H}(t_{i}+h)$ $\displaystyle=$ $\displaystyle
1-s_{H}(t_{i}+h)-i_{H}(t_{i}+h)$ (33) $\displaystyle i_{v}(t_{i}+h)$
$\displaystyle=$ $\displaystyle
i_{v}(t_{i})+\beta_{v}s_{v}(t_{i})i_{H}(t_{i})\,h-\mu_{v}i_{v}(t_{i})\,h,$
(34) $\displaystyle s_{v}(t_{i}+h)$ $\displaystyle=$ $\displaystyle
1-i_{v}(t_{i}+h),$ (35) $\displaystyle c(t_{i}+h)$ $\displaystyle=$
$\displaystyle c(t_{i})+p\beta_{H}s_{H}(t_{i})i_{v}(t_{i})\,h,$ (36)

for $\ 1\leq i\leq M-1$.

The numerical model (equations 31 to 36) was implemented using the M languages
of MATLAB.