ORIGINAL RESEARCH

published: 30 October 2018

doi: 10.3389/fphar.2018.01188

Non-linear Entropy Analysis in EEG

to Predict Treatment Response to
Repetitive Transcranial Magnetic
Stimulation in Depression
Reza Shalbaf 1 , Colleen Brenner 2 , Christopher Pang 1 , Daniel M. Blumberger 3,4 ,
Jonathan Downar 4,5 , Zafiris J. Daskalakis 3,4 , Joseph Tham 6 , Raymond W. Lam 6 ,
Faranak Farzan 7 and Fidel Vila-Rodriguez 1*
1
Non-Invasive Neurostimulation Therapies (NINET) Laboratory, Department of Psychiatry, University of British Columbia,
Vancouver, BC, Canada, 2 Department of Psychology, Loma Linda University, Loma Linda, CA, United States, 3 Temerty
Centre for Therapeutic Brain Intervention and Campbell Family Research Institute, Centre for Addiction and Mental Health,
Toronto, ON, Canada, 4 Department of Psychiatry, University of Toronto, Toronto, ON, Canada, 5 MRI-Guided rTMS Clinic
and Krembil Research Institute, University Health Network, Toronto, ON, Canada, 6 Department of Psychiatry, University
of British Columbia, Vancouver, BC, Canada, 7 School of Mechatronic Systems Engineering, Simon Fraser University, Surrey,
BC, Canada

Background: Biomarkers that predict clinical outcomes in depression are essential for
Edited by: increasing the precision of treatments and clinical outcomes. The electroencephalogram
Hector J. Caruncho, (EEG) is a non-invasive neurophysiological test that has promise as a biomarker sensitive
University of Victoria, Canada
to treatment effects. The aim of our study was to investigate a novel non-linear index of
Reviewed by:
resting state EEG activity as a predictor of clinical outcome, and compare its predictive
Jose Manuel Olivares,
Hospital Alvaro Cunqueiro, Spain capacity to traditional frequency-based indices.
Gopalkumar Rakesh,
Duke University, United States Methods: EEG was recorded from 62 patients with treatment resistant depression
Alejandro A. García Caballero, (TRD) and 25 healthy comparison (HC) subjects. TRD patients were treated with
Universidade de Santiago
de Compostela, Spain
excitatory repetitive transcranial magnetic stimulation (rTMS) to the dorsolateral
*Correspondence:
prefrontal cortex (DLPFC) for 4 to 6 weeks. EEG signals were first decomposed using
Fidel Vila-Rodriguez the empirical mode decomposition (EMD) method into band-limited intrinsic mode
[email protected]
functions (IMFs). Subsequently, Permutation Entropy (PE) was computed from the
Specialty section:
obtained second IMF to yield an index named PEIMF2. Receiver Operator Characteristic
This article was submitted to (ROC) curve analysis and ANOVA test were used to evaluate the efficiency of this index
Neuropharmacology,
(PEIMF2) and were compared to frequency-band based methods.
a section of the journal
Frontiers in Pharmacology Results: Responders (RP) to rTMS exhibited an increase in the PEIMF2 index compared
Received: 22 June 2018 to non-responders (NR) at F3, FCz and FC3 sites (p < 0.01). The area under the curve
Accepted: 28 September 2018
Published: 30 October 2018 (AUC) for ROC analysis was 0.8 for PEIMF2 index for the FC3 electrode. The PEIMF2
Citation: index was superior to ordinary frequency band measures.
Shalbaf R, Brenner C, Pang C,
Blumberger DM, Downar J,
Conclusion: Our data show that the PEIMF2 index, yields superior outcome prediction
Daskalakis ZJ, Tham J, Lam RW, performance compared to traditional frequency band indices. Our findings warrant
Farzan F and Vila-Rodriguez F (2018)
further investigation of EEG-based biomarkers in depression; specifically entropy indices
Non-linear Entropy Analysis in EEG
to Predict Treatment Response applied in band-limited EEG components. Registration in ClinicalTrials.Gov; identifiers
to Repetitive Transcranial Magnetic NCT02800226 and NCT01887782.
Stimulation in Depression.
Front. Pharmacol. 9:1188. Keywords: EEG, rTMS, major depressive disorder, permutation entropy, empirical mode decomposition,
doi: 10.3389/fphar.2018.01188 biomarker

Frontiers in Pharmacology | www.frontiersin.org 1 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

INTRODUCTION Complexity and non-linear behavior are characteristics of

typical brain functioning, (Elbert et al., 1994) and therefore the
Major depressive disorder (MDD) is a global public health application of non-linear dynamics analyses to the EEG signal
concern (World Health Organization [WHO], 2016) as the may prove to be a better measure of neural activity (Hosseini
disease is the leading cause of disability (Zarate et al., 2013) and et al., 2010). Recently, a non-linear parameter called permutation
affects approximately 5.4% of the population worldwide (Ferrari entropy (PE) (Bandt and Pompe, 2002; Cao et al., 2004) has
et al., 2013). Furthermore, it is estimated that 40% of those living been developed to dissect the complexity of EEG signals by
with the disorder do not respond to first line treatments such as deciphering the local order structure of a dynamical time series
pharmacological or psychosocial treatments and have treatment (Cao et al., 2004). In addition, PE properly tracks the dynamics
resistant depression (TRD) (Berlim et al., 2008). of brain activity (Shalbaf et al., 2013), is conceptually simple,
Repetitive transcranial magnetic stimulation (rTMS) is a safe computationally efficient, and robust against artifacts (Shalbaf
and effective treatment for TRD with 50–55% response and et al., 2015a).
30–35% remission rates (Galletly et al., 2011; McDonald et al., However, PE can be underestimated if the signals are
2011), and rTMS is considered a first-line treatment option superimposed with local or global trends. Research has suggested
for TRD (Milev et al., 2016). rTMS induces an electric field that properly removing the trends in biological signals with
in the brain strong enough to depolarize neurons and trigger a decomposition approach may improve the performance of
action potentials, and the treatment is delivered non-invasively non-linear signal analysis (Lo et al., 2009; Tsai et al., 2012).
by applying a coil in contact with the scalp. Neuroimaging studies Different variations of decompositions are suitable due to their
have shown a degree of hypoactivity in the left dorsolateral ability to derive dynamical features from the signals with
cortex (L-DLPFC) in MDD (Hallett, 2007), and therefore rTMS an enhanced resolution (Kevric and Subasi, 2017). One such
protocols that increase cortical excitability have been applied to method, empirical mode decomposition (EMD), was developed
the L-DLPFC for the treatment of MDD. for analyzing non-stationary data (Huang et al., 1998). EMD
The prescription of rTMS, similar to antidepressant can decompose a complicated signal without a basis function,
medication prescription, is currently based on clinical such as sine or wavelet functions, into several intrinsic mode
assessment and a process of trial and error. Identification functions (IMFs) that are embedded in the original signal. The
of effective biomarkers that can inform clinical decisions is decomposition procedure is adaptive, data-driven and highly
lacking, and this absence may contribute to higher health-care efficient (Kevric and Subasi, 2017). Therefore, entropy index
costs (Silverstein et al., 2015). Developing reliable biomarkers applied in band-limited EEG component extracted with EMD
may have profound implications for clinical practice as it would method may optimally quantify non-linear neuronal oscillations.
shift the prescription process to a more precise and personalized The purpose of this study is to examine rsEEG features
approach that would further improve clinical outcomes and as predictors of treatment response in TRD patients receiving
efficiency during treatment initiation (Collins et al., 2011). excitatory rTMS to the L-DLPFC. We hypothesized that rsEEG
The search for biomarkers of response has expanded decomposition components will hold different energies for
to molecular, neurophysiological and neuroimaging methods different patients and that these would differentiate responders
(Silverstein et al., 2015). The resting-state electroencephalogram (RP) from non-responders (NR). Furthermore, we hypothesized
(rsEEG) has merited particular interest due to its ease of that non-linear methods would be better suited and more efficient
use, cost-effectiveness and non-invasive nature which are predictors of rTMS treatment response compared to traditional
optimal characteristics for its implementation in clinical settings linear frequency-band power metrics.
(Baskaran et al., 2012; Shalbaf et al., 2015b; Wade and Iosifescu,
2016).
Several frequency-based rsEEG measures have been proposed MATERIALS AND METHODS
as predictors of response in TRD in the context of rTMS.
Examples include theta (4–7 Hz) activity in the subgenual Participants and Experimental
zone of the anterior cingulate cortex (Narushima et al., 2010), Procedures
anterior alpha (8–12 Hz) peak frequency (Arns et al., 2012), The neurophysiology dataset was part of two randomized, single-
prefrontal cordance (combination of absolute and relative EEG blinded trials in which patients with TRD were assigned to
power at different bands), (Bares et al., 2015; Erguzel et al., receive either intermittent theta burst stimulation (iTBS) or
2015) and Lempel-Ziv analysis on the alpha band (Arns et al., high frequency left (HFL) rTMS protocols to the left DLPFC
2014). However, these frequency-based methods are susceptible (parameters discussed below). Patients received a 4–6 week
to artifacts and are more suitable for the analysis of stationary course of rTMS.
signals. Furthermore, these frequency-based measures require Participants were part of two separate clinical trials
a Fourier transform, and this transform precludes precise with identical inclusion and exclusion criteria registered in
estimation of temporal patterns in EEG (Huang et al., 1998). ClinicalTrials.Gov, identifier NCT02800226 and NCT01887782
(Blumberger et al., 2018). All participants provided informed
Abbreviations: DLPFC, dorsolateral prefrontal cortex; EMD, empirical mode
consent and both experimental protocols were approved by both
decomposition; HDRS, hamilton depression rating scale; MDD, major depressive
disorder; ROC, receiver operating characteristic; rTMS, repetitive transcranial the UBC Clinical Research Ethics Board as well as the Vancouver
magnetic stimulation; TRD, treatment-resistant depression. Coastal Health Research Institute.

Frontiers in Pharmacology | www.frontiersin.org 2 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

Subjects were 62 TRD patients and 25 healthy comparison Pre-treatment EEG Acquisition
(HC) subjects. Demographic characteristics were collected at rsEEG was collected using Brain Products EEG systems (Brain
baseline (Table 1). All participants completed rsEEG at baseline, Products, Gilching, Germany) at two UBC sites part of the
prior to receiving treatment. Six patients did not complete Canadian Biomarker Integration Network in Depression (Lam
the full 4 weeks of rTMS treatment and were excluded from et al., 2016). The process of acquiring data with different
analysis. Data from an additional 5 patients were removed from systems has been carefully considered and addressed (Farzan
analysis due to random noise after quality control analysis. et al., 2017a). Continuous rsEEG was recorded using 31 (site
Thus, a dataset with a total of 76 participants was used in A) or 64 (site B) recording sites determined using the 10–20
this study, including 51 TRD patients (26 randomized to iTBS system of electrode placement, an EasyCap electrode cap, and
and 25 to HFL) and 25 HC. The inclusion and exclusion sintered Ag-AgCl electrodes. rsEEG data were recorded using
criteria of patients and HC are outlined in Supplementary a QuickAmp amplifier (Brain Products, Gilching, Germany;
Material. 1000 Hz A/D rate; 0.10 Hz high pass, 499 Hz low pass; common
average reference; impedances ≤ 10 k). rsEEG was obtain
Stimulation Technique and Parameters within 7 days of treatment initiation in all participants (mean
Stimulation techniques have been previously described (Ge 3.7 days).
et al., 2017). Briefly, a MagPro X100 stimulator with a Cool- Participants were given the same resting state instructions,
B70 fluid-cooled coil was used to deliver rTMS for all patients “Please close your eyes for 3 min while we collect your brain
(Magventure, Farum, Denmark). Resting motor threshold was activity at rest. Let your mind wander and try not to fall asleep.”
determined by visual inspection of right interpolicis brevis All rsEEGs were conducted in a sound-attenuated room with
muscle contraction with the aid of the TMS Motor Threshold reduced lighting to limit distraction and noise. Two sets of bipolar
Assessment Tool (Dobek et al., 2016). All treatments were electrodes were placed around the participant’s eyes for collecting
delivered at 120% resting motor threshold (Blumberger et al., Electrooculogram (EOG) to track eye movements for artifact
2018). Following randomization, TRD patients received either rejection.
HFL stimulation or iTBS over the left DLPFC, using a
Neuronavigation system (Visor 2.0, ANT Neuro, Enschede, EEG Preprocessing
Netherlands) and the target location specified by reverse Two levels of pre-processing steps were implemented in order
coregistration from a stereotaxic coordinate on the standard to standardize the EEG data collected from both sites. These
Montreal neurological institute (MNI-152) template brain steps are done in MATLAB (The Mathworks, Inc., Natick, MA,
[x − 38, y + 44, z + 26] identified as optimal based on United States) via the open-source EEGLAB toolbox (Delorme
functional connectivity and clinical outcome (Dobek et al., and Makeig, 2004).
2016). The aim of the first pre-processing step was to minimize
raw data heterogeneity across two sites and prepare the data
Clinical Measures for integration. First, data from Site B were reduced because
Primary clinical outcome was measured using the 17-item site B had more recording electrodes, and thus some electrodes
Hamilton Depression Rating Scale (HDRS). For each patient, were removed from analyses to match the same number of
HDRS scores were collected at baseline and at the end of the electrodes as Site A. The electrode locations in both caps
rTMS course. Interviewers were blinded to patient treatment were identical as per manufacturer description (Brain Products,
allocation. Responders (RP) were defined as those having a Gilching, Germany in both site A and site B). Second, data from
50% or greater reduction in HDRS scores between baseline Site B were re-referenced to common average reference such that
and end of treatment. Out of the 51 patients included in the data from the two sites possess the equivalent electrode reference.
analysis, there were 31 responders and 20 non-responders to Then, length of recording for all participants modified to have the
rTMS treatment. same length. Also, separate frequency analysis and statistical tests

TABLE 1 | Demographics and clinical characteristics of HC and TRD patients by responder and non-responder groups.

TRD patients (n = 51)

Responders (n = 31) Non-responders (n = 20) Healthy volunteers (n = 25) p

Sex (F/M) 19/12 12/8 18/7 nsa

Age (SD) 43.4 (11.6) 42.5 (14.1) 39.8 (13.1) nsb
Years of education (SD) 15.2 (2.5) 15.0 (1.9) 16.0 (1.9) nsb
Handedness (R/L/A) 26/5/0 15/3/2 21/2/2 nsa
HDRS (SD) 22.5 (4.3) 2249 (4.0) − nsb
Treatment (HFL/iTBS) 13/18 12/8 − nsa

M, male; F, female; SD, standard deviation; L, left; R, right; A, ambidextrous; HDRS, 17-item Hamilton Depression Rating Scale. a Chi-square test. b Two sample t-test to
compare healthy controls to patients, since in the original analysis, a one-way ANOVA was used to compare age and education across the 3 groups.

Frontiers in Pharmacology | www.frontiersin.org 3 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

from the same two healthy volunteers were done to show data non-stationary series can be transformed to an almost stationary
were equivalent in quality across the two sites. ordinal series. The smallest and the largest values of PE are zero
The second pre-processing step was to implement EEG artifact and one, with zero reflecting a highly regular time series and
removal, since these artifacts interfere with the identification of one reflecting equal probability of all permutations. The detailed
true neurophysiological signal. First, the sampling rate of all data algorithm and parameter sets have been previously published
was decreased from 1000 to 256 samples per second to reduce (Shalbaf et al., 2013).
white noise. Second, data segments contaminated with large-
amplitude or random noise sources that cannot be extracted Permutation Entropy Intrinsic Mode Functions
through filtering were removed with GUI workflow. Third, the (PEIMF2)
high and low band pass filters were set to at 0.5 and 55 Hz The EEG signals were first decomposed by applying the
respectively, to remove low and high frequency noise. The notch EMD method into symmetric and band-limited IMFs which
filter was also set at 60 Hz to remove industrial noise. Finally, are arranged from high to low frequency components. EMD
blind source separation techniques via independent component decomposition was computed over segments of 8 s with an
analysis (ICA) (Makeig et al., 1996) were used to extract eye overlap of 6 s in order to consistently track the transient changes
movements and blinks, muscle activity, and cardiac signals in in the EEG recording. Then, PE was computed from each of the
order to separate neural activity from these sources of noise. 8-s IMF2 segments. The average PE of all segments in a recording
was considered the PEIMF2 index for each participant.
Linear EEG Analyses
Time-frequency analyses using short time windowed Fast Fourier Statistical Analysis
Transform (FFT) applied to resting EEG epochs were computed Differences in neurophysiological variables between RP, NR, and
using MATLAB and EEGLAB software (Delorme and Makeig, HC groups were examined using 1-way analysis of variance
2004). Some conventional frequency band measures such as Delta (ANOVA). The normality of the data was investigated before
(1–4 Hz), Theta (4–8 Hz), Alpha (8–12 Hz), Beta (12–24 Hz) performing analyses, and a p-value of 0.01 was set as the
and Gamma (30–50 Hz) relative powers were extracted from the criteria for statistical significance for greater stringency than
rsEEG signals for both HC and TRD groups. 1024 points discrete conventional levels.
FFT with a 100% Hanning window was computed over segments A Receiver Operator Characteristic (ROC) curve was plotted
of 8 s with an overlap of 6 s and the average of all segments in as a two-dimensional depiction of the classifier’s performance in
a recording was considered the relative powers index for each predicting treatment outcomes using the proposed biomarkers.
participant. The two axes of this graph represent tradeoffs between errors
(false positives) and successes (true positives) that a classifier
Non-linear EEG Analyses makes between two classes. To corroborate the results of this
Non-linear features were extracted from the eyes-closed resting analysis, the area under the ROC curve, abbreviated as AUC was
state EEG signals of both HC and TRD groups. Details of the calculated.
algorithm will be described in following sections.

Empirical Mode Decomposition (EMD) RESULTS

EMD is a method of signal analysis (Huang et al., 1998) that has
recently been applied to biological signals (Shalbaf et al., 2012). Demographic and Clinical
Using EMD, any complex signal can be decomposed into a small
number of intrinsic mode functions (IMFs) through a sifting
Characteristics
Table 1 summarizes the demographic and clinical characteristics
process. The IMFs should fulfill two requirements:
of all participants. The responders, non-responders, and healthy
(1) The IMF is symmetric with respect to the local mean. comparison groups had similar age, sex, years of education, and
(2) Each IMF has the same number of zero crossings and handedness.
extrema, or they differ at most by one. There were no differences between RP and NR as well as
between MDD and HC in age and sex. Also Baseline HDRS scores
The detailed algorithm has been previously described and the
of responders were similar to that of non-responders.
number of IMFs was established to be six (Shalbaf et al., 2012).
An EEG epoch of 8 s from one participant in FC3 electrode is
plotted in Figure 1A. The EMD of this EEG epoch is composed
Response Prediction Based on PEIMF2
of six IMFs which are shown in Figure 1B. These IMFs are almost Index
orthogonal components. The calculated PEIMF2 indices for the RP, NR and HC groups are
plotted onto scalp topographic maps in Figure 2 with scales to the
Permutation Entropy (PE) right of the maps. Greater PEIMF2 index are observed in RP and
Permutation entropy (PE) is a new non-linear parameter that HC groups compared to NR, especially at left frontal sites.
quantifies the amount of regularity in EEG data (Bandt and One-way ANOVAs were calculated to investigate whether
Pompe, 2002; Cao et al., 2004). This feature converts a given there were significant differences between RP and NR groups on
EEG series into a sequence of ordinal patterns, meaning that a PEIMF2 index for each electrode site (Table 2). The ANOVA

Frontiers in Pharmacology | www.frontiersin.org 4 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

FIGURE 1 | A segment of EEG signal from one participant in FC3 site (A) [X(t)] and EMD of the same segment (B, Imf 1 to Imf 6).

FIGURE 2 | Scalp topographical maps of PEIMF2 index (resting state, eyes closed). From left to right: topographies of the (A) RP, (B) the NR, and (C) HC. PEIMF2
index in RP and HC groups is higher than NR group especially at left frontal electrodes.

showed that there were significant differences in the PEIMF2 (0.028), HC = 0.533 (0.027) (mean (standard deviation))], FCz
index between the two groups in FC4 (p = 0.009), FCz (p = 0.001), (RP = 0.511 (0.021), NR = 0.496 (0.029), HC = 0.508 (0.024)),
F3 (p = 0.002), F4 (p = 0.003), Fz (p = 0.003), CP3 (p = 0.0049), and F3 (RP = 0.534 (0.037), NR = 0.501 (0.022), HC = 0.525
P3 (p = 0.005) and FC3 (p < 0.001) electrodes. Also, as (0.028)) electrode sites. These results suggest that the PEIMF2
shown in Figure 3, the largest differences between RP and index may be able to differentiate between RP and NR groups,
NR were observed at FC3 [RP = 0.543 (0.033), NR = 0.521 particularly at the frontal regions. There were no significant

Frontiers in Pharmacology | www.frontiersin.org 5 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

differences between RP and HC participants at these three have been consistently associated with healthy states where the
electrode sites (p-value > 0.01), but there was a marginally nervous system is able to respond and adapt to dynamic changes.
significant difference between NR and HC (p-value << 0.01), Conversely, lower entropy values (more regular, less information)
where aHC exhibited a higher PEIMF2 index as compared to NR. are associated with pathological states and loss of the prime
Results were unchanged when taking in consideration treatment ability of the nervous system to respond to changes (Kevric and
group as a covariate (i.e., HFL vs. iTBS stimulation). Subasi, 2017). A plausible hypothesis to explain the absence of
In this study, ROC curve analyses were used to explore difference between RP and HC would be that RP still have a
the optimum component of EMD to undergo PE calculation system that is capable of such changes and this would make these
in order to best differentiate between RP and NR patients. patients amenable to respond to rTMS. Lower entropy levels in
PE was computed for IMF1 to IMF4 to extract indices called NR participants may indicate a reduction in typical intra-cortical
PEIMF1 to PEIMF4 respectively, for FC3 electrode and multiple information flow, a more regular, less complex EEG in the left
surrounding electrodes since this area best differentiated RP and frontal region and fewer chances for the system to change in
NR. The AUC value for PEIMF2 index is 0.8, compared with response to rTMS. Therefore, lower levels of entropy (i.e., more
0.71 for PEIMF1, 0.76 for PEIMF3, and 0.74 for PEIMF4 in FC3 regularity) may reflect a less preserved brain function that is
electrode and similar result gained from other electrodes. The not amenable to the effect of rTMS. This raises the question
AUC value of the ROC analysis classifying RP and NR was the as to whether these patients would be amenable to a different
greatest for PEIMF2 during resting state EEG, suggesting that PE type of stimulation (e.g., inhibitory rTMS), a different anatomical
calculated on the second IMF yields the best results for prediction target (e.g., dorsomedial prefrontal cortex), or a higher dose
of treatment response with moderate accuracy. of excitatory rTMS (e.g., more pulses per session or per day,
accelerated protocols). Furthermore, our findings are convergent
Response Prediction Based on with those recently reported by Jaworska et al. (2018) who
Frequency Band Measures found that increased diffuse multi-scale entropy was predictive of
Some conventional frequency band measures such as Delta, treatment response to antidepressants and Farzan et al. (2017b)
Theta, Alpha, Beta and Gamma relative powers were extracted who showed that non-linear complexity measures was superior
from EEG signal of all electrodes. The efficiency of these indices to power in explaining the therapeutic efficacy of seizure therapy.
was evaluated via AUC values on the best electrode for the The electrode sites that significantly distinguish RP from NR
classification of RP or NR. The AUC values of Delta at Oz, Theta are located in left frontal (F3 and FC3), right frontal (F4, FC4),
at O1, Alpha at Oz, Beta at CPz and Gamma at O2 are 0.67, left parietal (CP3, P3) and central (FCz, Fz) sites (Table 2).
0.64, 0.63, 0.60, and 0.68 respectively (Figure 4). The result show A plausible explanation of our results is that the antidepressant
that there was a considerable difference between predictive value effects of rTMS to the DLPFC are not restricted only to the local
of PEIMF2 index (AUC = 0.8) and traditional frequency band effect on the L-DLPFC, but rather that target circuits that underlie
measures on best electrode site. complex brain functions (e.g., affect regulation) including the
frontal gyrus, anterior cingulate cortex, amygdala, and insula
(Sheline et al., 2010; Smart et al., 2015). The current result
DISCUSSION converge with previous findings of brain areas related to MDD
(Silverstein et al., 2015; Milev et al., 2016; Ge et al., 2017) and
The current study evaluated predictors of treatment response further extends it by adding new evidence that an entropy index
to rTMS administered to the left DLPFC in patients with TRD indicating the complexity in a network could potentially serve as
based on rsEEG signal. Our data shows that a non-linear entropy a predictor of clinical response to rTMS.
index, PEIMF2, yields superior outcome prediction performance There are conflicting reports regarding the predictive capacity
compared to traditional frequency band indices. Our data of frequency band metrics in the EEG such as alpha and
indicate that TRD patients who responded to rTMS had higher theta activity. Some studies showed there was no correlation
entropy compared to NR, with the most prominent differences in alpha activity and treatment response (Price et al., 2008;
appearing in prefrontal areas (FCz, F3, and FC3 electrodes). Widge et al., 2013) whereas another showed that there was a
Our findings extend previous the investigation of EEG-based negative correlation between the two (Micoulaud-Franchi et al.,
biomarkers in depression, and position entropy indices applied 2012). Moreover, one group reported that increased slow theta
in band-limited EEG components extracted with EMD method activity in the subgenual zone of the anterior cingulate cortex
as a potential predictor for clinical use. was correlated with positive response to rTMS (Hallett, 2007)
Entropy is becoming a valuable tool for the analysis of EEG while another group reported that theta rhythm increase in
activity and has received much attention in recent years in the frontal cortex is associated with non-response to rTMS
the study of brain disorders (Mizuno et al., 2010). Our data treatment (Silverstein et al., 2015). Our own data would be
consistently show that NR patients have significantly lower convergent with the idea of a moderate predictive ability of
entropy values in the prefrontal areas, and particularly the frequency band (linear) metrics (Figure 4). One possible reason
DLPFC region, compared to RP and HC subjects (Figure 3). for the inferior predictive ability of linear metrics may be
Entropy indicates the complexity in a system, (Erguzel et al., related to the non-linear nature of neural processes, as threshold
2015) and is also associated with the amount of “information” and saturation phenomena control the dynamical behavior of
the signal carries. In the nervous system, higher levels of entropy individual neurons (Narushima et al., 2010). In contrast, our

Frontiers in Pharmacology | www.frontiersin.org 6 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

TABLE 2 | Results of the ANOVA investigating differences in PEIMF2 index between RP and NR groups at all electrode sites (α = 0.01).

Electrode site P-value Electrode site P-value Electrode site P-value

FPz 0.107 FT7 0.108 TP8 0.024

FP2 0.123 FT8 0.481 Pz 0.275
FP1 0.120 Cz 0.022 P3 0.005
Fz 0.003 C3 0.01 P4 0.040
F3 0.002 C4 0.038 P7 0.014
F4 0.003 T7 0.524 P8 0.014
F7 0.053 T8 0.889 Oz 0.043
F8 0.054 CPz 0.251 O1 0.076
FCz 0.001 CP3 0.004 O2 0.012
FC3∗ <0.001 CP4 0.084
FC4 0.009 TP7 0.566

FC3 asterisk denotes the most significant site for differentiating between RP and NR. Bold terms denote sites that significantly differentiate RP from NR.

FIGURE 3 | PEIMF2 index as a function of electrode sites for RP and NR with comparison to HC participant groups. Error bars represent ± 1 standard error. RP
differ than NR to rTMS treatment especially at FC3.

data supports the utility of non-linear metrics in predicting between RP and NR groups are observed in the left frontal
treatment outcome. Our measure quantified higher order non- electrodes (Figure 3), which is consistent with previous findings
linear complexity, which is not obtained using traditional EEG of brain areas related to MDD. Considering this information,
spectral-band analyses such as alpha or theta band power. we would speculate that a plausible mechanism mediating the
PEIMF2 measure is based on non-linear dynamics and has response to rTMS may be neuroplastic changes on relevant
been found to indirectly index neuroplasticity (Hayley et al., circuits involved in affect regulation and other symptomatic
2005). PEIMF2 may represent the excitatory and inhibitory domains.
balance of the related networks in MDD which would also EMD adaptively and locally decomposes non-stationary EEG
be associated with neuroplasticity. The scalp topography for signals into a sum of IMFs that represent amplitude- and
the PEIMF2 values show that the most prominent differences frequency- modulated components specific to the energy levels of

Frontiers in Pharmacology | www.frontiersin.org 7 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

FIGURE 4 | ROC curve analysis and AUC value of Delta, Theta, Alpha, Beta, Gamma relative powers to discriminate between RP an NR on best electrodes. The low
AUC value of frequency band measures indicates weak prediction accuracy of these linear approaches.

individual patients (Shalbaf et al., 2012). EMD has shown better 2007). The former could lead to uniform or deviated IMFs, and
properties over other methods of EEG decomposition such as the the latter may result in problems overshooting or distorting the
short-time Fourier transform (STFT), independent component beginning and ending of signals. Second, we believe that a multi-
analysis (ICA) (Makeig et al., 1996) and wavelet transform modal approach of response prediction that holistically integrates
(WT) (Zikov et al., 2006). For instance, the STFT excludes a variety of sources of data including clinical, neuroimaging,
EEG features with a short duration or narrow frequency band. and neurophysiological measures may be most reliable because
ICA is hampered by the intrinsically non-stationary nature and many clinical factors may affect the central nervous system
the non-linear couplings involved in neural signal generation. including baseline physiological and neurological differences,
Moreover, WT forces the decomposition of the signal into a thus decreasing the predictive power of related EEG measures.
pre-defined set of basic functions, therefore temporal patterns of Finally, we acknowledge this work provides a preliminary proof-
EEG signals cannot be obtained precisely. Conversely, EMD is a of-principle evidence and the real consistency will only be
decomposition technique that is completely data-driven and thus determined with future larger trials or replication of the two
utilizes empirical knowledge of oscillations intrinsic to the given methods of analysis.
time series (Baskaran et al., 2012). Furthermore, unlike wavelet To conclude, this study addresses a new method to decompose
analysis, EMD does not depend on a fixed set of basic functions; the neuronal oscillations with EMD to obtain a series of IMFs.
instead it searches for IMF embedded within the data. Therefore, The results show that the permutation entropy measure applied
these pitfalls of the other methods, or the strength of the data- to the second IMF yields the optimal result in predicting
drive EMD method may make it a better biomarker of treatment treatment response. This seems to indicate that second IMF
response. oscillation plays an important role in discriminating between RP
Some limitations of our work should be considered in order to and NR in the context of rTMS treatment. The method described
better interpret our results. First, the decomposition procedure of herein certainly merits further research in larger samples to
EMD requires the arbitrary choice of the stopping criteria for the replicate and improve its predictive power. Prospective studies
sifting and the spline-fitting scheme (Sweeney-Reed and Nasuto, applying our predictor to decide treatment interventions may be

Frontiers in Pharmacology | www.frontiersin.org 8 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

used in the future, perhaps facilitating the precise prescription FUNDING

of rTMS as a first-line treatment when several favorable
neurophysiological predictive factors are present, thus avoiding The study was partly funded by the Philanthropic
unsuccessful pharmacotherapy trials and expediting recovery. donation to the Non-Invasive Neurostimulation Therapies
laboratory at UBC and the Canadian Institutes of
Health Research. MagVenture provided in-kind equipment
ETHICS STATEMENT support; however, MagVenture had no role in the study
design, data analysis, interpretation, or preparation of
This study was carried out in accordance with the this manuscript and none of the investigators receive any
recommendations and approvals of UBC Clinical Research Ethics financial compensation or have any financial interests in
Board as well as the Vancouver Coastal Health Research Institute MagVenture.
with written informed consent from all subjects in accordance
with the Declaration of Helsinki. Trials were registered in
ClinicalTrials.Gov, identifier NCT02800226 and NCT01887782. ACKNOWLEDGMENTS
The authors thank the clinical research staff and the patient
AUTHOR CONTRIBUTIONS participants.

RS, CB, and FV-R conceived and designed the study. DB, ZD,
and RL provided input on the study design. RS, CB, and FV-R SUPPLEMENTARY MATERIAL
developed the plan for statistical analyses. RS analyzed the data.
All authors contributed to the interpretation of data. RS and The Supplementary Material for this article can be found
FV-R, drafted the manuscript. All authors made revisions to the online at: https://www.frontiersin.org/articles/10.3389/fphar.
manuscript. 2018.01188/full#supplementary-material

REFERENCES Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for
analysis of single-trial EEG dynamics including independent component
Arns, M., Cerquera, A., Gutiérrez, R. M., Hasselman, F., and Freund, J. A. (2014). analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.
Non-linear EEG analyses predict non-response to rTMS treatment in major 10.009
depressive disorder. Clin. Neurophysiol. 125, 1392–1399. doi: 10.1016/j.clinph. Dobek, C. E., Dipinto, A. A., Kuan, A. J., Lam, R. W., Blumberger, D. M.,
2013.11.022 Downar, J., et al. (2016). Sinus tachycardia induced by methocarbamol and
Arns, M., Drinkenburg, W. H., Fitzgerald, P. B., and Kenemans, J. L. (2012). repetitive transcranial magnetic stimulation (rTMS). Brain Stimul. 9, 156–158.
Neurophysiological predictors of non-response to rTMS in depression. Brain doi: 10.1016/j.brs.2015.10.012
Stimul. 5, 569–576. doi: 10.1016/j.brs.2011.12.003 Elbert, T., Ray, W. J., Kowalik, Z. J., Skinner, J. E., Graf, K. E., and
Bandt, C., and Pompe, B. (2002). Permutation entropy: a natural complexity Birbaumer, N. (1994). Chaos and physiology: deterministic chaos in
measure for time series. Phys. Rev. Lett. 88:174102. doi: 10.1103/PhysRevLett. excitable cell assemblies. Physiol. Rev. 74, 1–47. doi: 10.1152/physrev.1994.
88.174102 74.1.1
Bares, M., Brunovsky, M., Novak, T., Kopecek, M., Stopkova, P., Sos, P., et al. Erguzel, T. T., Ozekes, S., Gultekin, S., Tarhan, N., Hizli Sayar, G., and Bayram, A.
(2015). QEEG theta cordance in the prediction of treatment outcome to (2015). Neural network based response prediction of rTMS in major depressive
prefrontal repetitive transcranial magnetic stimulation or venlafaxine ER in disorder using QEEG cordance. Psychiatry Investig. 12, 61–65. doi: 10.4306/pi.
patients with major depressive disorder. Clin. EEG Neurosci. 46, 73–80. doi: 2015.12.1.61
10.1177/1550059413520442 Farzan, F., Atluri, S., Frehlich, M., Dhami, P., Kleffner, K., Price, R., et al.
Baskaran, A., Milev, R., and McIntyre, R. S. (2012). The neurobiology of (2017a). Standardization of electroencephalography for multi-site, multi-
the EEG biomarker as a predictor of treatment response in depression. platform and multi-investigator studies: insights from the Canadian biomarker
Neuropharmacology 63, 507–513. doi: 10.1016/j.neuropharm.2012. integration network in depression. Sci. Rep. 7:7473. doi: 10.1038/s41598-017-07
04.021 613-x
Berlim, M., Fleck, M., and Turecki, G. (2008). Current trends in the assessment and Farzan, F., Atluri, S., Mei, Y., Moreno, S., Levinson, A. J., Blumberger, D. M.,
somatic treatment of resistant/refractory major depression: an overview. Ann. et al. (2017b). Brain temporal complexity in explaining the therapeutic and
Med. 40, 149–159. doi: 10.1080/07853890701769728 cognitive effects of seizure therapy. Brain 140, 1011–1025. doi: 10.1093/brain/
Blumberger, D. M., Vila-Rodriguez, F., and Thorpe, K. E. (2018). Effectiveness awx030
of theta burst versus high-frequency repetitive transcranial magnetic Ferrari, A., Charlson, F., Norman, R., and Patten, S. (2013). Burden of depressive
stimulation in patients with depression (THREE-D): a randomised non- disorders by country, sex, age, and year: findings from the global burden of
inferiority trial. Lancet 391, 1683–1692. doi: 10.1016/S0140-6736(18) disease study 2010. PLoS Med. 10:e1001547. doi: 10.1371/journal.pmed.100
30295-2 1547
Cao, Y., Tung, W., Gao, J. B., Protopopescu, V. A., and Hively, Galletly, C., Gill, S., Clarke, P., Burton, C., and Fitzgerald, P. B. (2011).
L. M. (2004). Detecting dynamical changes in time series using the A randomized trial comparing repetitive transcranial magnetic stimulation
permutation entropy. Phys. Rev. E 70:46217. doi: 10.1103/PhysRevE.70.04 given 3 days/week and 5 days/week for the treatment of major depression:
6217 is efficacy related to the duration of treatment or the number of treatments?
Collins, P. Y., Patel, V., Joestl, S. S., March, D., Insel, T. R., and Daar, A. S. Psychol. Med. 13, 1–8.
(2011). Scientific advisory board and the executive committee of the grand Ge, R., Blumberger, D. M., Downar, J., Daskalakis, Z. J., Dipinto, A. A., Tham,
challenges on global mental health. Nature 475, 27–30. doi: 10.1038/47 J. C. W., et al. (2017). Abnormal functional connectivity within resting-
5027a state networks is related to rTMS-based therapy effects of treatment resistant

Frontiers in Pharmacology | www.frontiersin.org 9 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

depression: a pilot study. J. Affect. Disord. 218, 75–81. doi: 10.1016/j.jad.2017. Shalbaf, R., Behnam, H., and Jelveh Moghadam, H. (2015a). Monitoring
04.060 depth of anesthesia using combination of EEG measure and hemodynamic
Hallett, M. (2007). Transcranial magnetic stimulation: a primer. Neuron 55, 187– variables. Cogn. Neurodyn. 9, 41–51. doi: 10.1007/s11571-014-
199. doi: 10.1016/j.neuron.2007.06.026 9295-z
Hayley, S., Poulter, M. O., Merali, Z., and Anisman, H. (2005). The Shalbaf, R., Behnam, H., Sleigh, J. W., Steyn-Ross, D. A., and SteynRoss,
pathogenesis of clinical depression: stressor- and cytokine-induced alterations M. L. (2015b). Frontal-temporal synchronization of EEG signals quantified
of neuroplasticity. Neuroscience 135, 659–678. doi: 10.1016/j.neuroscience. by order patterns cross recurrence analysis during propofol anesthesia. IEEE
2005.03.051 Trans. Neural Syst. Rehabil. Eng. 23, 468–474. doi: 10.1109/TNSRE.2014.235
Hosseini, P. T., Shalbaf, R., and Nasrabadi, A. M. (2010). Extracting a seizure 0537
intensity index from one-channel EEG signal using bispectral and detrended Shalbaf, R., Behnam, H., Sleigh, J. W., Steyn-Ross, A., and Voss, L. J. (2013).
fluctuation analysis. J. Biomed. Sci. Eng. 3, 253–261. doi: 10.4236/jbise.2010. Monitoring the depth of anesthesia using entropy features and an artificial
33034 neural network. J. Neurosci. Methods 218, 17–24. doi: 10.1016/j.jneumeth.2013.
Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., et al. (1998). 03.008
The empirical mode decomposition and the Hilbert spectrum for nonlinear Shalbaf, R., Behnam, H., Sleigh, J. W., and Voss, L. J. (2012). Using
and non-stationary time series analysis. Proc. R Soc. London A Math. Phys. the Hilbert–Huang transform to measure the electroencephalographic
Eng. Sci. 454, 903–995. Available at: http://rspa.royalsocietypublishing.org/ effect of propofol. Physiol. Meas. 33, 271–285. doi: 10.1088/0967-3334/33/
content/454/1971/903.short?casa_token=xk10C2yVljYAAAAA:UlYn6Ov8s83 2/271
FZ4j8FVKWnNwZBqc0QD9zWiHWiv1n-s8NbdfbIZBb878bjpipOZDaFSYux Sheline, Y. I., Price, J. L., Yan, Z., and Mintun, M. A. (2010). Resting-state functional
kn9KUNjTEINAg MRI in depression unmasks increased connectivity between networks via the
Jaworska, N., Wang, H., Smith, D. M., Blier, P., Knott, V., and Protzner, A. B. dorsal nexus. Proc. Natl. Acad. Sci. U.S.A. 107, 11020–11025. doi: 10.1073/pnas.
(2018). Pre-treatment EEG signal variability is associated with treatment 1000446107
success in depression. Neuroimage Clin. 17, 368–377. doi: 10.1016/j.nicl.2017. Silverstein, W. K., Noda, Y., Barr, M. S., Vila-Rodriguez, F., Rajji, T. K.,
10.035 Fitzgerald, P. B., et al. (2015). Neurobiological predictors of response to
Kevric, J., and Subasi, A. (2017). Comparison of signal decomposition dorsolateral prefrontal cortex repetitive transcranial magnetic stimulation in
methods in classification of EEG signals for motor-imagery BCI system. depression: a systematic review. Depress. Anxiety 32, 871–891. doi: 10.1002/da.
Biomed. Signal. Process. Control 31, 398–406. doi: 10.1016/j.bspc.2016. 22424
09.007 Smart, O. L., Tiruvadi, V. R., and Mayberg, H. S. (2015). Multimodal approaches
Lam, R. W., Milev, R., Rotzinger, S., Andreazza, A. C., Blier, P., Brenner, C., et al. to define network oscillations in depression. Biol. Psychiatry 77, 1061–1070.
(2016). Discovering biomarkers for antidepressant response: protocol from doi: 10.1016/j.biopsych.2015.01.002
the Canadian biomarker integration network in depression (CAN-BIND) and Sweeney-Reed, C. M., and Nasuto, S. J. (2007). A novel approach to the detection
clinical characteristics of the first patient cohort. BMC Psychiatry 16:105. doi: of synchronisation in EEG based on empirical mode decomposition. J. Comput.
10.1186/s12888-016-0785-x Neurosci. 23, 79–111. doi: 10.1007/s10827-007-0020-3
Lo, M.-T., Novak, V., Peng, C.-K., Liu, Y., and Hu, K. (2009). Nonlinear phase Tsai, P. H., Lin, C., Tsao, J., Lin, P. F., Wang, P. C., Huang, N. E., et al.
interaction between nonstationary signals: a comparison study of methods (2012). Empirical mode decomposition based detrended sample entropy in
based on Hilbert-Huang and Fourier transforms. Phys. Rev. E 79:61924. doi: electroencephalography for Alzheimer’s disease. J. Neurosci. Methods 210, 230–
10.1103/PhysRevE.79.061924 237. doi: 10.1016/j.jneumeth.2012.07.002
Makeig, S., Bell, A., and Jung, T. (1996). Independent component analysis of Wade, E., and Iosifescu, D. (2016). Using electroencephalography for treatment
electroencephalographic data. Adv. Neural 8:8. guidance in major depressive disorder. Biol. Psychiatry Cogn. Neurosci. 1,
McDonald, W. M., Durkalski, V., Ball, E. R., Holtzheimer, P. E., Pavlicova, M., 411–422.
Lisanby, S. H., et al. (2011). Improving the antidepressant efficacy of Widge, A. S., Avery, D. H., and Zarkowski, P. (2013). Baseline and treatment-
transcranial magnetic stimulation: maximizing the number of stimulations emergent EEG biomarkers of antidepressant medication response do not
and treatment location in treatmentresistant depression. Depress. Anxiety 28, predict response to repetitive transcranial magnetic stimulation. Brain Stimul.
973–980. doi: 10.1002/da.20885 6, 929–931. doi: 10.1016/j.brs.2013.05.001
Micoulaud-Franchi, J.-A., Richieri, R., Cermolacce, M., Loundou, A., World Health Organization [WHO] (2016). World Health Statistics: Monitoring
Lancon, C., and Vion-Dury, J. (2012). Parieto-temporal alpha EEG Health for the SDGs Sustainable Development Goals. Geneva: WHO.
band power at baseline as a predictor of antidepressant treatment Zarate, C. A., Mathews, D. C., and Furey, M. L. (2013). Human biomarkers of rapid
response with repetitive transcranial magnetic stimulation: a preliminary antidepressant effects. Biol. Psychiatry 73, 1142–1155. doi: 10.1016/j.biopsych.
study. J. Affect Disord. 137, 156–160. doi: 10.1016/j.jad.2011. 2012.11.031
12.030 Zikov, T., Bibian, S., and Dumont, G. (2006). Quantifying cortical activity during
Milev, R. V., Giacobbe, P., Kennedy, S. H., Blumberger, D. M., Daskalakis, general anesthesia using wavelet analysis. IEEE Trans. 53, 617–632.
Z. J., Downar, J., et al. (2016). Canadian Network for Mood and Anxiety
Treatments (CANMAT) 2016 Clinical guidelines for the management Conflict of Interest Statement: DB reports research grants from the Canadian
of adults with major depressive disorder: section 4. neurostimulation Institutes of Health Research (CIHR), US National Institutes of Health, Weston
treatments. Can. J. Psychiatry 61, 561–575. doi: 10.1177/070674371666 Brain Institute, Brain Canada, the Temerty Family Foundation (through the
0033 Centre for Addiction and Mental Health Foundation and the Campbell Research
Mizuno, T., Takahashi, T., Cho, R., and Kikuchi, M. (2010). Assessment Institute), and Brainsway; reports receiving in-kind equipment support for
of EEG dynamical complexity in Alzheimer’s disease using multiscale investigator-initiated studies (including this study) MagVenture; is the site
entropy. Clin. Neurophysiol. 121, 1438–1446. doi: 10.1016/J.CLINPH.2010. principal investigator for three sponsor-initiated studies for Brainsway; and has
03.025 been on an advisory board for Janssen Pharmaceutical. ZD reports research grants
Narushima, K., McCormick, L. M., Yamada, T., Thatcher, R. W., and Robinson, and equipment in-kind support for an investigator-initiated study from Brainsway
R. G. (2010). Subgenual Cingulate theta activity predicts treatment response and Magventure. JD reports research grants from CIHR, the National Institute
of repetitive transcranial magnetic stimulation in participants with vascular for Mental Health, Brain Canada, the Canadian Biomarker Integration Network
depression. J. Neuropsychiatry Clin. Neurosci. 22, 75–84. doi: 10.1176/jnp.2010. in Depression, the Ontario Brain Institute, the Klarman Family Foundation, the
22.1.75 Arrell Family Foundation, and the Edgestone Foundation; reports travel stipends
Price, G. W., Lee, J. W., Garvey, C., and Gibson, N. (2008). Appraisal of from Lundbeck and ANT Neuro; reports in-kind equipment support for this
sessional EEG features as a correlate of clinical changes in an rTMS investigator-initiated trial from MagVenture; and is an advisor for and is an advisor
treatment of depression. Clin. EEG Neurosci. 39, 131–138. doi: 10.1177/ for BrainCheck. RL reports research grants or consulting or speaking honoraria
155005940803900307 from Akili Interactive, Asia-Pacific Economic Cooperation, Allergan, AstraZeneca,

Frontiers in Pharmacology | www.frontiersin.org 10 October 2018 | Volume 9 | Article 1188

Shalbaf et al. Entropy for Depression Treatment-Response Prediction

Bristol-Myers Squibb, Canadian Depression Research and Intervention Network, Copyright © 2018 Shalbaf, Brenner, Pang, Blumberger, Downar, Daskalakis, Tham,
Canadian Network for Mood and Anxiety Treatments, Johnson and Johnson, Lam, Farzan and Vila-Rodriguez. This is an open-access article distributed under
Lundbeck, Lundbeck Institute, MagVenture, Pfizer, St Jude Medical, Otsuka, and the terms of the Creative Commons Attribution License (CC BY). The use,
Takeda. FV-R reports research grants from CIHR, Brain Canada, Michael Smith distribution or reproduction in other forums is permitted, provided the original
Foundation for Health Research, and Vancouver Coastal Health Research Institute; author(s) and the copyright owner(s) are credited and that the original publication
reports receiving in-kind equipment support for this investigator-initiated trial in this journal is cited, in accordance with accepted academic practice. No
from MagVenture; and has been on an advisory board for Janssen. CB, FF, RS, JT, use, distribution or reproduction is permitted which does not comply with these
and CP declare no competing interests. terms.

Frontiers in Pharmacology | www.frontiersin.org 11 October 2018 | Volume 9 | Article 1188