# Permutation time irreversibility in sleep electroencephalograms: Dependence
on sleep stage and the effect of equal values

Wenpo Yao  State Key Laboratory of Organic Electronics and Information
Displays, Institute of Advanced Materials, School of Chemistry and Life
Sciences, School of Geographic and Biologic Information, Nanjing University of
Posts and Telecommunications, Nanjing 210023, China  Key Laboratory of
Computational Neuroscience and Brain-Inspired Intelligence(Fudan University),
Ministry of Education

###### Abstract

Time irreversibility (TIR) refers to the manifestation of nonequilibrium brain
activity influenced by various physiological conditions; however, the
influence of sleep on electroencephalogram (EEG) TIR has not been sufficiently
investigated. In this paper, a comprehensive study on permutation TIR (pTIR)
of EEG data under different sleep stages is conducted. Two basic ordinal
patterns (i.e., the original and amplitude permutations) are distinguished to
simplify sleep EEGs, and then the influences of equal values and forbidden
permutation on pTIR are elucidated. To detect pTIR of brain electric signals,
5 groups of EEGs in the awake, stages I, II, III, and rapid eye movement (REM)
stages are collected from the public Polysomnographic Database in PhysioNet.
Test results suggested that the pTIR of sleep EEGs significantly decreases as
the sleep stage increases (p$<$0.001), with the awake and REM EEGs,
demonstrating greater differences than others. Comparative analysis and
numerical simulations support the importance of equal values. Distribution of
equal states, a simple quantification of amplitude fluctuations, significantly
increases with the sleep stage (p$<$0.001). If these equalities are ignored,
incorrect probabilistic differences may arise in the forward-backward and
symmetric permutations of TIR, leading to contradictory results; moreover, the
ascending and descending orders for symmetric permutations also lead different
outcomes in sleep EEGs. Overall, pTIR in sleep EEGs contributes to our
understanding of quantitative TIR and classification of sleep EEGs.

time irreversibility; amplitude permutation; distribution of equal states;
sleep classification; EEG

††preprint: APS/123-QED

##  I Introduction

The human brain, which contains large numbers of neurons and interacts with
other physiological organs Andrea  _et al._ (2018); Bashan  _et al._ (2012),
is a highly complex system. Brain activity exhibits evidently nonlinear,
nonequilibrium properties and is influenced by various internal and external
factors. Sleep is a vital physiological activity for living beings; various
changes occur in the physiological systems during sleep, with brain electric
activities demonstrating pronounced changes. Therefore, the complex
characteristics of sleep electroencephalograms (EEGs) can help in analyzing
the sleep mechanisms Zhao  _et al._ (2019); Phan and Mikkelsen (2022). Various
nonlinear methods, such as fractal approaches and detrended fluctuation
analysis, have been used to explore the dynamics of EEGs during different
sleep stages. Entropy measures Ma  _et al._ (2018); Xiong  _et al._ (2017) are
widely employed to detect the complexity of sleep brain signals, such as the
Shannon entropy Hou  _et al._ (2021); Bandt (2017), sample entropy Liang  _et
al._ (2021); Wang  _et al._ (2021), and multiscale entropy Miskovic  _et al._
(2019). Among these complex characteristics, the loss of time reversibility is
a manifestation of nonequilibrium brain electric activity Fang  _et al._
(2019); Lynn  _et al._ (2021) and has been applied to several neural
pathological conditions, including epilepsy Yao  _et al._ (2020a, 2021a);
Martinez  _et al._ (2018), Alzheimer’s disease Martin Gonzalo  _et al._
(2019), and alcoholism Zanin (2021).

Time irreversibility (TIR) Weiss (1975), also defined as temporal asymmetry
(TAS) Kelly (1979), is an important characteristic of nonequilibrium EEGs. To
quantity TIR, the probabilistic differences between the forward-backward
processes or symmetry vectors must be measured Weiss (1975); Kelly (1979);
Ramsey and Rothman (1995); both are nontrivial. In real-world signal
processing, TIR is generally quantified by coarse-graining the time series.
Probabilistic differences between up and down were introduced by Costa Costa
_et al._ (2005, 2008), Porta Porta  _et al._ (2008), and Ehlers et al. Guzik
_et al._ (2006); Ehlers  _et al._ (1998) as simplified temporal asymmetries.
Lacasa et al. Lacasa  _et al._ (2012); Lacasa and Flanagan (2015); Flanagan
and Lacasa (2016) estimated the irreversibility considering the
distinguishability between the in-out distributions of the visibility graph
Lacasa  _et al._ (2008). The in-out difference was then employed to detect
irreversible characteristics in effective interactions and particular
networked connectivity Donges  _et al._ (2013); Zou  _et al._ (2019). Given
nonequilibrium statistical mechanics in neuronal spike trains, the network
inference and couplings of asymmetric models were investigated via
asynchronous update Zeng  _et al._ (2013, 2011). Symbolic TIR based on a way
of coarse-graining or reduction of description is widely adopted in
quantitative TIR owing to its computational effectiveness and simplified
statistical analysis Daw  _et al._ (2003); Cammarota and Rogora (2007). Among
these coarse-graining methods, ordinal patterns convey structural dynamics and
do not impose further model assumptions Bandt and Pompe (2002); Bandt (2016,
2020); Yao  _et al._ (2022); thus, they are popular in the quantification of
TIR Yao  _et al._ (2020a, 2021a); Martinez  _et al._ (2018); Martin Gonzalo
_et al._ (2019); Zanin (2021). However, several challenges are observed in
permutation TIR (pTIR). It should be noted that there are two basic ordinal
patterns, i.e., original permutation and amplitude permutation; they
differently reflect the structural dynamics of any series Yao  _et al._
(2022). Amplitude permutation directly reflects the vector temporal structure,
while the application of original permutation might result in conceptual
errors in pTIR Yao  _et al._ (2021a, 2019a, 2023a). Forbidden permutation
refers to the missing ordinal pattern in the simplified transformation of a
temporal structure. Amigo et al. Amigo  _et al._ (2006, 2007, 2008)
extensively analyzed the features of these forbidden ordinal patterns and
proposed several methods for chaotic determination in real-world time series
analysis. The distributions of forbidden permutation convey important system
information, such as chaotic dynamics, correlation, and nonlinearity Zanin
(2008); Carpi  _et al._ (2010); Kulp  _et al._ (2017); Cuesta Frau (2020).
Forbidden permutation is an adverse factor in pTIR because it generates
individual permutations and make division-based parameters (e.g., the
Kullback-Leibler distance) unsuitable to calculate the probabilistic
difference between forbidden and individual permutations Yao  _et al._ (2020a,
2021a, 2019a, 2023a). Equal value is another implicated factor because it
significantly affects the construction of ordinal patterns and permutation
analysis. Ma et al. Bian  _et al._ (2012) focused on the indexes of equal
values in permutation and modified them into the same symbol (rank) for a more
accurate characterization of the system structure. Zunino et al. Zunino  _et
al._ (2017) reported erroneous conclusions on permutation entropy (PEn) in the
event of equalities in time series; David et al. Cuesta Frau  _et al._ (2018)
further identified the weakness of PEn, i.e., the possible ambiguities
introduced by equal values in the subsequences. However, the influence of
equal values on pTIR has not been sufficiently analyzed. Equal values are
generally assumed to be rare if a process has a continuous distribution Bandt
and Pompe (2002). This is partly true; however, equal values are observed in
physiological signals such as heartbeats Yao  _et al._ (2019b, 2020b) and raw
EEGs Yao  _et al._ (2021b, 2023b), and they significantly impact permutation
analysis. Moreover, equal values might generate self-symmetry vectors
containing important physical implication (i.e., time reversibility or
temporal symmetry) in TIR and produce contradictory findings in real-world
series analysis Yao  _et al._ (2019b, 2020b). Furthermore, the different
treatment of equal values might lead to inconsistencies between forward-
backward and symmetric permutations, thus yielding differences in the
quantification of TIR and TAS of time series. Therefore, the effects on the
pTIR during signal processing, particular the TIR and TAS in real-world series
analysis, should be comprehensively studied. Owing to these unfavorable
factors in the quantification of TIR, the time irreversible characteristics of
sleep EEGs have not been given the deserved attention.

To address this problem, a comprehensive analysis of pTIR and its application
in sleep EEG classification is performed. Accordingly, two basic ordinal
patterns are compared to simplify the time series, and several crucial factors
in pTIR are clarified. Subsequently, the distributions of equal values and
individual permutations that affected the pTIR of sleep EEGs are detected. For
comparison, nonequal permutations of TIR are employed, and PEn is applied to
evaluate the complexity of sleep EEGs. Furthermore, several issues are
discussed such as equal values in sleep EEGs and the relationship between
entropy and TIR. The contributions of this research are as follows: (1) it
elucidates the key factors of pTIR in physiological signal processing,
especially the influence of equal values, (2) it explores the time
irreversible characteristics of EEGs under different sleep conditions, which
is helpful for sleep classification.

##  II Methods

###  II.1 Time irreversibility

A process is defined as time reversible if it is invariant under the reversal
of the timescale; otherwise, it is time irreversible. Given below are two
statistical definitions of time reversibility:

Definition 1. According to Weiss Weiss (1975), a stationary process, $X(t)$,
is time reversible if $\\{X(t_{1}),X(t_{2}),\ldots,X(t_{m})\\}$ and
$\\{X(-t_{1}),X(-t_{2}),\ldots,X(-t_{m})\\}$ have the same joint probability
distributions for every $t_{1},t_{2},\ldots,t_{m}$ and $m$; otherwise, $X(t)$
is time irreversible.

Definition 2. Based on the works of Kell Kelly (1979), if $X(t)$ is time
reversible, $\\{X(t_{1}),X(t_{2}),\ldots,X(t_{m})\\}$ and
$\\{X(-t_{1}+n),X(-t_{2}+n),\ldots,X(-t_{m}+n)\\}$ have the same probability
distributions for every $t_{1},t_{2},\ldots,t_{m}$ and $n$. Particularly,
under $n=t_{1}+t_{m}$, symmetric $\\{X(t_{1}),X(t_{2}),\ldots,X(t_{m})\\}$ and
$\\{X(t_{m}),\ldots,X(t_{2}),X(t_{1})\\}$ have the same joint probabilities.
Moreover, the symmetric form of a vector is the same as its counterpart in the
time-reversal series. Therefore, TIR is also defined as TAS.

To quantify the TIR or TAS of a process, the forward-backward probabilistic
difference and symmetric vectors’ probabilistic divergence should be
equivalent Ramsey and Rothman (1995); however, their applications in real-
world processes are different Yao  _et al._ (2020a, 2021a). The forward-
backward approach for TIR is operationally convenient and more reliable in
application. Meanwhile, it should have the entire process to obtain the
reversed one, which is not feasible if the process is large or uninterrupted;
therefore, TIR based on forward-backward differences does not satisfy real-
time requirement. By contrast, TAS based on the probabilistic difference of
symmetric vectors demonstrates real-time performance; therefore, it has higher
applicability in physiological and environmental condition monitoring. Note
that TAS based on symmetric vectors has a limitation, i.e., the vector should
be faithfully associated with its alternative; otherwise, conceptual
misleading may occur Yao  _et al._ (2021a, 2019a, 2023a).

###  II.2 Original and amplitude permutations

TIR can be measured by alternatively calculating the joint probabilistic
differences of simplified processes instead of the raw process. Among these
simplified measures, pTIR is particularly popular. The ordinal pattern comes
naturally from the series and does not pose further model assumptions Bandt
and Pompe (2002); Yao  _et al._ (2022); hence, it plays an important role in
quantitative TIR.

According to the representation of a vector structure, there are two basic
ordinal patterns, i.e., the original permutation (OrP) and amplitude
permutation (AmP) Yao  _et al._ (2022, 2023a). The OrP consists of the indexes
of reorganized values in the original series, while the AmP comprises the
positions of the original values in the reordered series. Given series
$X(i)=\\{x(i_{1}),\ldots,x(i_{i}),\ldots,x(i_{m})\\}$, it is reordered in the
ascending or descending order to
$X(j)=\\{x(j_{1}),\ldots,x(j_{j}),\ldots,x(j_{m})\\}$ e.g., increased as
$x(j_{1})<\cdots<x(j_{j})<\cdots<x(j_{m})$. Then, the OrP and AmP are
generated according to the organized indexes of the original and reordered
series as Eq. (3).

$\displaystyle\left\\{\begin{array}[]{lr}\mbox{OrP}:(x_{1,j},\ldots,x_{i-1,j},x_{i,j},x_{i+1,j},\ldots,x_{m,j})\\\
\mbox{AmP}:(x_{i,1},\ldots,x_{i,j-1},x_{i,j},x_{i,j+1},\ldots,x_{i,m})\end{array}\right.$
(3)

In the construction of OrP, $i$ increases from 1 to $m$ and $j$ is the
location of reorganized values in the original series $X(i)$, i.e.,
$OrP_{j}=(j_{1},j_{2},\ldots,j_{i-1},j_{i},j_{i+1},\cdots,j_{m-1},j_{m})$. In
the generation of AmP, $j$ increases from 1 to $m$ and $i$ is the position of
original values in the reordered series $X(j)$, i.e.,
$AmP_{i}=(i_{1},i_{2},\ldots,i_{j-1},i_{j},i_{j+1},\cdots,i_{m-1},i_{m})$.
Taking a series with five values $X(i)=\\{5,1,7,3,9\\}$ as an example, the
reordered series can be represented as $X(j)=\\{1,3,5,7,9\\}$ in the ascending
order. The indexes of reorganized $X(j)$ values in the original $X(i)$ is
OrP=(2,4,1,3,5); the positions of the original $X(i)$ values in the
reorganized $X(j)$ is AmP=(3,1,4,2,5).

Equal values are not rare in real-world signals and are usually generated
owing to limitations in signal collection Yao  _et al._ (2021b, 2023b),
especially the quantization error in analog-to-digital conversion (ADC). Equal
values have an important role in the construction of ordinal patterns and
permutation analysis Bian  _et al._ (2012); Zunino  _et al._ (2017); Cuesta
Frau  _et al._ (2018); Yao  _et al._ (2019b, 2020b); therefore, their indexes
in ordinal patterns should be improved accordingly. If there are equal values
in series $X(i)$, they can be organized in neighboring orders according to
their order of occurrence; for example,
$\cdots<x(i_{1},j_{1})=x(i_{2},j_{2})<\cdots<x(i_{3},j_{3})=x(i_{4},j_{4})=x(i_{5},j_{5})<\cdots$.
Then, the indexes of equal values can be rewritten to be the same in each
group, such as to the smallest indexes as
$\cdots<x(i_{1},j_{1})=x(i_{1},j_{1})<\cdots<x(i_{3},j_{3})=x(i_{3},j_{3})=x(i_{3},j_{3})<\cdots$
Bian  _et al._ (2012) or the largest ones as
$\cdots<x(i_{2},j_{2})=x(i_{2},j_{2})<\cdots<x(i_{5},j_{5})=x(i_{5},j_{5})=x(i_{5},j_{5})<\cdots$,
and modify the OrP and AmP accordingly. Further taking a series with five
values $X(i)=\\{5,1,9,1,7\\}$ as an example, the second ‘1’ could be treated
as ‘2’ according to its order of occurrence; then, the OrP and AmP become
(2,4,1,5,3) and (3,1,5,2,4) in the ascending order, respectively. In equal-
value ordinal patterns, the indexes of equal values are modified to be the
smallest; consequently, the OrP and AmP are improved to be (2,2,1,5,3) and
(3,1,5,1,4), respectively. Without equal values, if the length of the series
is $m$, there exist $m$! ordinal patterns, and there are more motifs if the
equality is considered. Equal-value permutation is necessary for the
comprehensive reflection of the series’ temporal structure Yao  _et al._
(2022, 2023a). Figure 1 illustrates the structures of triple-value series and
their OrPs and AmPs.

Figure 1: OrP and AmP of triple-value series. AmPs are bold and underlined.
Indexes of red equal values in the ascending order are modified to be the
smallest ones in their corresponding groups. a) OrPs and AmPs of all-up, all-
down, and all-equal vectors are always the same. b) Symmetric AmPs of time-
symmetric ($y$-axis) series are connected by dashed black arrows; symmetric
OrPs of $x$-axis symmetric series are connected by solid red arrows. c) AmPs
of time self-symmetry vectors are symmetric.

Figure 1 displays the comprehensive temporal structures of triple-value series
where equal values have an important role. Furthermore, OrP and AmP are two
basic ordinal patterns that differently convey the temporal structure of the
series Yao  _et al._ (2022). AmP directly reflects the temporal structure of
the series because its elements directly correspond to the amplitude of the
original sequence element. Therefore, AmPs of symmetric vectors (e.g.,
$y$-axis or time symmetric) are all symmetric in Fig. 1, and AmPs of the three
temporally self-symmetric vectors are also self-symmetric. By contrast, OrP
indirectly reflects the original series structure owing to targets on the
reordered series, and OrPs of amplitude-symmetric (i.e., $x$-axis symmetry)
series are symmetric. The differences between OrP and AmP differently affect
permutation analysis. When ordinal patterns are used as an alternative label
or symbol of the series, there are no differences between the two basic
permutations, such as in PEn Bandt (2017); Bandt and Pompe (2002). Otherwise,
if ordinal patterns are constructed as a direct replacement of the series
structure, OrP and AmP should be dealt with carefully to avoid possible
errors, such as in pTIR Yao  _et al._ (2021a, 2019a, 2023a). In this study,
AmP is used as a direct alternative for the temporal structure of the series
to avoid possible errors in pTIR.

###  II.3 Permutation time irreversibility

Permutation TIR is the probabilistic difference between forward-backward or
symmetry permutations as an alternative to the original series. Owing to the
advantages of ordinal patterns and their application in quantitative TIR, pTIR
is widely used in time series analysis. In real world signal analysis, the
application of pTIR encounters some challenges, as listed below:

(1) OrP is not a direct reflection of the temporal structure of vectors, as
shown in Fig. 1, and its effect on TAS is such that if the probabilistic
difference between symmetric vectors is calculated, symmetric OrPs should not
be used. The AmPs of time-symmetric vectors are symmetric, and symmetric AmPs
could be employed as alternatives for evaluating TIR Yao  _et al._ (2022).
Meanwhile, the OrPs of amplitude-symmetric vectors are symmetric, and
symmetric OrPs should be employed as alternatives for measuring amplitude
irreversibility Yao  _et al._ (2021a, 2023a). It should be noted that when
measuring the joint probabilistic difference of the forward-backward process
for TIR, there is no difference between the OrP and AmP. Overall, AmP directly
reflects the temporal structure of vectors and is recommended for pTIR.

(2) The consideration of equal values is necessary to construct comprehensive
and reliable vector structures; more importantly, equal values can generate
self-symmetric vectors. Self-symmetric vectors, e.g., the vectors and their
AmPs in boxes in Fig. 1a) and c), have a special physical implication, i.e.
time reversibility or temporal symmetry Yao  _et al._ (2022, 2023a).
Furthermore, equal values widely exist in physiological data (either directly
collected EEGs or indirect heartbeats derived from electrocardiography) owing
to the nonlinear and irreversible quantization process of signal collection
Yao  _et al._ (2019b, 2021b, 2023b). Therefore, equal values should not be
broken by adding small random perturbations or ranked according to their order
of appearance, and their indexes in ordinal patterns should be modified
accordingly.

(3) Forbidden permutations are a special type of forbidden symbol, i.e.,
symbols that do not exist for a process. Forbidden permutations are closely
related to system characteristics Zanin (2008); Carpi  _et al._ (2010); Kulp
_et al._ (2017); Cuesta Frau (2020); Amigo  _et al._ (2006, 2007, 2008), but
they negatively impact quantitative TIR. Among the pairs of permutations for
quantitative TIR, if there exists a forbidden permutation, its corresponding
permutation is an individual permutation. The mathematical difference between
the probabilities of forbidden and individual permutations is zero or infinite
if using division-based parameters, such as the Kullback-Leibler distance.
Therefore, such division-based parameters are not suitable for quantitative
TIR Yao  _et al._ (2020a, 2021a, 2019a, 2023a). It is recommended to use
subtraction-based parameters, such as $Y_{s}$ in Eq. (4), in quantitative TIR
to calculate the probabilistic differences, where $p_{i}$ and $p_{j}$
probability of corresponding permutations and $p_{i}$ should not be less than
$p_{j}$.

$\displaystyle Y_{s}\langle
p_{i},p_{j}\rangle=p_{i}\frac{p_{i}-p_{j}}{p_{i}+p_{j}},$ (4)

(4) TIR and TAS are statistically consistent; theoretically, the probabilistic
differences between symmetric permutations and those between forward-backward
permutations should be the same. However, the existence of equal values and
the traditional treatment that ranks them according to their order of
emergence make TIR and TAS different. Figure 2 illustrates a comparative
construction of nonequal and equal-value AmPs considering equal values.

Figure 2: AmPs of triple-value series considering equalities. Equal values are
represented in red and the crossed greens are their alternatives according to
their order of occurrence in the ascending order. In equal-value AmPs (bold
and underlined), indexes of red equal values are modified to be the smallest
in their corresponding groups.

As evident in Fig. 2, nonequal AmP represents a temporal structure of vector
without equal values. Equal values are transformed into ‘false up’ in nonequal
ordinal patterns and their symmetric permutation is ‘real down’, but their
corresponding permutations in the backward series are ‘false up’ and ‘real
up’. Note that ‘real down’ in the forward series is ‘real up’ in the backward
series; moreover, the probability distributions of ‘up’ in the forward and
backward series are both incorrect. Furthermore, if there are equal values in
the process and nonequal permutation is applied, neither TIR nor TAS can be
correctly quantified; in addition, the AmPs of symmetric and forward-backward
vectors are not the same, leading to inconsistent and even contradictory TIR
and TAS in the quantitative nonequilibrium of signal processing, which will be
confirmed in this study.

##  III Results

EEG is a typically complex signal and subject to different physiological
activities. In this section, EEG data under different sleep conditions are
collected from the publicly available PhysioNet Goldberger  _et al._ (2000) to
test the sleep stages on pTIR. The probabilistic differences of AmPs are
calculated using the subtraction-based $Y_{s}$ to quantify TIR; TIR and TAS
based on equal-value AmP are denoted as pTIR and pTAS, while those based on
nonequal AmP are represented as noeTIR and noeTAS, respectively.

###  III.1 Sleep EEGs

The MIT-BIH Polysomnographic Database is a collection of sleep physiologic
data. This database contains information related to the sleep physiological
signals of 16 male subjects (age ranging from 32 to 56, mean age 43; weight
ranging from 89 to 152 kg, mean weight 119 kg). The database contains over 80
h of four-, six-, and seven-channel polysomnographic recordings, with a
standard expert annotation for sleep stages after every 30 s according to the
criteria of Rechtschaffen and Kales Rechtschaffen and Kales (1968). EEG
signals are recorded from the C4-A1, O2-A1, or C3-O1 channel at a sampling
rate of 250 Hz and 12-bit quantization. Referring to the annotation files,
sleep EEGs in five stages, i.e., awake and sleep stages 1 (SI), 2 (SII), 3
(SIII), and rapid eye movement (REM) are extracted; each stage contained 45
sets of EEG signals with a duration of 60 s (15000 points) after visual
inspection for artifacts. More detailed information can be found in Ref
Ichimaru and Moody (1999); Goldberger  _et al._ (2000). The nonparametric
Mann–Whitney U test is performed to test the statistical differences in pTIR
between each of the two stages of sleep EEGs, and Kruskal–Wallis analysis of
variance is applied to measure those in pTIR of sleep EEGs in five stages.

Equal values play an important role in the construction of ordinal patterns
and might significantly change the probability distribution of permutations;
moreover, they convey important physical implication, i.e., time reversibility
Yao  _et al._ (2020a, 2021a, 2022, 2023a). The distribution of equal states
(DES) Yao  _et al._ (2021b) is measured by Eq. (5), where $L$ denotes the
length of series, $N(s(t)=s(t+\tau))$ represents the number of neighboring
equal states with delay $\tau$. The states of DES indicate EEG values in this
report. DES of five groups of sleep EEGs are shown in Fig. 3.

$\displaystyle\mbox{DES}=\frac{N(s(t)=s(t+\tau))}{L-\tau}$ (5)

Figure 3: Exemplary EEGs in awake, SI–SIII, and REM states and the
distribution of equal states (mean$\pm$standard error) of sleep EEGs. The
state in DES is the direct amplitude value during sleep EEGs. # indicates
p$<$0.0001 across all stages using the Kruskal–Wallis test and * suggests
p$<$0.01 between each two stages of sleep EEGs using Mann–Whitney U test.

The results in Fig. 3 are consistent with those obtained in our previous
report on the distribution of equal values in sleep EEGs Yao  _et al._
(2023b). Raw sleep EEGs have numerous neighboring equal values. Although awake
EEGs have minimum equal values, their DES are 12%, and as the sleep stages
increase, the DES increase significantly; hence, REM EEGs have almost 34%
equal values. Under acceptable data recording resolution, awake EEGs have
larger amplitude fluctuations that produce fewer neighboring equal values, and
as the depth of sleep increases, amplitude fluctuations of the brain’s
electrical activity decrease, thus significantly decreasing the DES. Equal
values are generated owing to the limitation of the ADC, i.e., the zero-
amplitude fluctuation Yao  _et al._ (2021b, 2023b). Test results suggest that
the DES of EEGs under the five sleep conditions are statistically
significantly different (p$<$0.01 for Mann–Whitney U test; p$<$0.0001 for
Kruskal–Wallis test), and the DES of REM EEGs are particularly different
between others (p$<$0.00001). Therefore, equal values not only have important
effects on pTIR, but their distribution also serves as a simple parameter for
time-domain feature extraction and should not be ignored.

The existence of forbidden permutations makes division-based parameters
unsuitable for pTIR, and their distribution has been widely proved to be
closely associated with systematic information Zanin (2008); Carpi  _et al._
(2010); Kulp  _et al._ (2017); Cuesta Frau (2020); Amigo  _et al._ (2006,
2007, 2008). In sleep EEGs, all permutations have their corresponding forms
when $m$=2 and 3, while the distribution of individual permutations is rare
when $m$=4. When $m$ is 5 and larger, there exist forbidden as well as
individual permutations. Further, forbidden permutations contain false
forbidden permutations that do not exist because the data length is short, and
they decay with the sequence length Amigo  _et al._ (2007). Considering false
forbidden permutations, if the selected EEG data length is short, the pTIR of
classified sleep EEGs might not be reliable. The distribution of individual
permutations (DIPs) is expressed in Eq. (6), where $N(\pi_{I})$ and $N(\pi)$
represent the amount of individual permutation $\pi_{I}$ and existing
permutation $\pi$, respectively. Taking $m$=5 as an example, the effect of
data length on the DIPs of five groups of sleep EEGs is illustrated in Fig. 4.

$\displaystyle\mbox{DIP}=N(\pi_{I})/N(\pi)$ (6)

Figure 4: Distribution of individual permutations (mean$\pm$standard error) of
wake, SI–SIII, and REM EEGs.

Fig. 4 shows that sleep EEG signals contain individual permutations whose
probability distributions are affected by the signal length. With the increase
in data length from 1 to 10 s (2500 points), the DIP values of sleep EEGs
decrease stepwise and tend to converge when the data length was larger than 20
s (5000 points). Therefore, when the data length is short, there exist false
forbidden as well as false individual permutations, and they decay when the
data length increases. The DIP is also related to sleep conditions when the
EEG length is more than 20 s. SIII and REM EEGs have larger DIPs when awake
and SI EEGs exhibit smaller DIPs. The existence of forbidden and individual
permutations suggests that subtraction-based index is a necessity in
quantitative TIR. Moreover, the length of EEG signals cannot be too low; they
should be greater than 30 s in this study, otherwise these false individual
permutations will lead to unreliable pTIR for the nonequilibrium analysis of
sleep EEGs.

Overall, there are equal values and individual permutations in sleep EEGs that
affect pTIR analysis. Hence, the two basic ordinal patterns must be
distinguished considering the necessity of equal-value permutation, the
subtraction-based probabilistic difference and requirement of data length also
should be paid attentions in pTIR of sleep EEGs.

###  III.2 pTIR in sleep EEGs

Given the DIPs in sleep EEGs, dimension is set to 2, 3, and 4 in an
enumerative manner. The increase in delay is equivalent to the reduction in
signal sampling frequency, i.e., 250/$\tau$ Hz Yao  _et al._ (2021b, 2023b).
To satisfy the Nyquist sampling rate (i.e., more than twice the signal band),
the delay is set to 1–4, thus maintaining sufficient information about sleep
conditions in EEG data. The pTIR of sleep EEGs in awake, SI–SIII, and REM
stages are shown in Fig. 5.

Figure 5: pTIR (mean$\pm$standard error) of wake, SI–SIII, and REM EEGs. #
indicates p$<$0.0001 across all stages using Kruskal–Wallis test and *
suggests p$<$0.01 between pTIR of sleep EEGs and others using Mann–Whitney U
test.

Contrary to the DES in Fig. 3, the pTIR of sleep EEGs exhibit decreasing
trends as the sleep depth increase, as shown in Fig. 5. As the subjects come
into sleep from wakefulness, the pTIR of EEG signals significantly decrease;
moreover, as the sleep depth increase, the pTIR consistently decrease.
Particularly, when subjects enter the REM state, the pTIR of EEGs show a
considerable decrease. The choice of dimension and delay do not affect the
trend of EEGs’ pTIR with the increasing sleep depth, but it affect the
statistical discriminations. According to statistical results, when $m$=3 and
$\tau$=1, the pTIR of EEGs in the five sleep stages exhibit optimal
classification (p$<$0.0001 for Kruskal–Wallis test). As shown in Fig. 5, more
significant reductions are observed in the pTIR of EEGs from the awake to
sleep and SIII to REM stages when $m$=2 and 3 and $\tau$=1. Therefore, time
irreversible features of brain electric activity decrease as the sleep stages
advance, and more significant differences are observed in nonequilibrium
features during the awake–sleep transformation and REM state. For comparison,
the probabilistic differences between symmetric AmPs, i.e., pTAS, are also
calculated. The pTAS and pTIR yielded the same results, suggesting that equal-
value AmPs reliably characterize the temporal structure of vectors in sleep
EEGs and are not affected by the reverse process in backward time series.

It should be noted that equal values significantly affect the pTIR, and sleep
EEGs generally contain numerous equal values Yao  _et al._ (2023b). Given the
amplitude fluctuations measured by DES in Fig. 3, the noeTIR and noeTAS of
EEGs under the five sleep stages are calculated and shown in Fig. 6.

Figure 6: Contradictory noeTIR and noeTAS (mean$\pm$standard error) in wake,
SI–SIII, and REM EEGs. # indicates p$<$0.0001 across all stages using
Kruskal–Wallis test and * suggests p$<$0.01 between each of the two stages of
sleep EEGs using Mann–Whitney U test.

As evident in Fig. 6, noeTIR and noeTAS of the five groups of sleep EEGs
considerably differed and even showed completely contradictory results based
on the nonequal AmP. The comparison shows that noeTIR results exhibited a
consistent classification of sleep stages with the pTIR, while noeTAS had
contradictory results. Although pTIR and noeTIR demonstrated similar
decreasing trends with the increase in sleep depth, they were not the same in
these EEGs. Statistically, noeTIR of sleep EEGs was not significantly
different, while noeTAS with $m$=3 and $\tau$=1 effectively differed in the
five groups of EEGs (p$<$0.01 for Mann–Whitney U test; p$<$0.0001 for
Kruskal–Wallis test). Figure 2 shows that when the signal contained equal
elements and was reorganized in ascending order, the probabilistic differences
of permutations between the forward and backward series were closer to the
pTIR, while those of symmetric permutations had a greater deviation. Note that
neither noeTIR nor noeTAS yielded correct results. Taking $m$=2 as an example,
only up, down, and equal forms are observed, among which, equality was
transformed into up in both the forward and backward series. The pTIR, noeTIR,
pTAS, and noeTAS were calculated using Eq. (11), where $up_{p}$, $down_{p}$,
and $equal_{p}$ represent the probability of up, down, and equal values,
respectively, and $up_{p}/down_{p}$ denotes the probabilistic difference
between $up_{p}$ and $down_{p}$ using $Y_{s}$.

$\displaystyle\left\\{\begin{array}[]{lr}\mbox{pTIR}:0.5*(up_{p}/down_{p}+down_{p}/up_{p})\\\
\mbox{noeTIR}:0.5*(\frac{up_{p}+equal_{p}}{down_{p}+equal_{p}}+down_{p}/up_{p})\\\
\mbox{pTAS}:up_{p}/down_{p}\\\
\mbox{noeTAS}:(up_{p}+equal_{p})/down_{p}\end{array}\right.$ (11)

In Eq. (11), pTIR is the same as pTAS, but different from noeTIR and noeTAS.
Due to equal being mistaken as up, both noeTIR and noeTAS wrongly measure the
time reversibility of sleep EEGs. Evidently, if there is no equality, pTIR,
noeTIR, pTAS, and noeTAS are the same.

To compare with the pTIR, PEn Bandt and Pompe (2002); Yao  _et al._ (2020b),
i.e., the Shannon entropy of permutation probability $p(\pi)$, given in Eq.
(12), is further employed in the analysis of sleep EEGs. It is a widely
applied statistical parameter used to measure the nonlinear characteristics of
complex systems Bandt and Pompe (2002). The selection of OrP and AmP does not
affect the calculation of PEn because PEn is the mean information contained in
the existing ordinal patterns Yao  _et al._ (2022), while equal values affect
PEn in signal analysis Bian  _et al._ (2012); Zunino  _et al._ (2017); Cuesta
Frau  _et al._ (2018). PEn based on equal-value AmP with delay from 1 to 4 and
dimension ranging from 2 to 4 of sleep EEGs is displayed in Fig. 7.

$\displaystyle\mbox{PEn}=-\sum p(\pi)lnp(\pi)$ (12)

Figure 7: Permutation entropy (mean$\pm$standard error) of wake, SI-SIII, and
REM EEGs based on equal-value AmP.

As evident from Figs. 7 and 5, PEn is observed to be the opposite of pTIR in
the five groups of sleep EEGs. When $m$=2, PEn of EEGs increases as the
subjects fell asleep and sleep depth increased, suggesting that EEGs have
larger complexity with the increase of sleep stages. As $m$ increases, no
consistent trend is observed in the PEn complexity of EEGs with the increasing
sleep stages. Statistically, the PEn of five groups of sleep EEGs is not
significantly different. Therefore, pTIR is more reliable for sleep stage
classification owing to the quantification of EEG nonequilibrium features.

The contradictory results of pTIR and PEn can be explained by their
statistical concepts. Shannon entropy and TIR compute the difference in
probability distributions, but they focus on different sets of probabilistic
differences. PEn calculates the average amount of information contained in the
distribution of permutations, i.e., the probabilistic differences in all
existing permutations, whereas pTIR measures the difference between forward
and backward permutation series or that between symmetric permutations of a
series. Therefore, when the permutation probability difference is smaller, the
entropy complexity of signals will be higher, while the TIR will be smaller.
This is the fundamental reason for the opposite outcomes of pTIR and PEn in
sleep EEGs, which is consistent with the contradictory results obtained in our
previous report on heartbeats Yao  _et al._ (2020b). Moreover, because they
characterize complexity and nonequilibrium features, the results exhibit the
diversity of features from different perspectives of complex physiological
signals and enable us to explore complex systems more comprehensively.

The comparative analysis of sleep EEGs demonstrated that the TIR, TAS, and
Shannon entropy have a complex relationship, as shown in Fig. 8. TIR and TAS
based equal-value AmP share same results while they yielded contradictory
outcomes based on nonequal AmP. The trend of noeTIR is consistent with that of
pTIR, while noeTAS exhibited different results in the five groups of sleep
EEGs. PEn also presented contradictory findings in comparison to pTIR. To
determine the difference between these associated measures and the influence
of equal values, probability distributions of AmPs are further analyzed
numerically. Taking $m$=2 and $\tau$=1 as an example, the probabilities of
AmPs with and without equal values of five groups of sleep EEGs are shown in
Fig. 8.

Figure 8: Relationship of permutation analysis in sleep EEGs and probability
distributions (mean$\pm$standard error) of AmPs when dimension is 2 and delay
is 1. a) Consistent results of pTIR and pTAS are linked by black solid arrows,
while contradictory results of noeTIR and noeTAS are connected by crossed red
dashed arrows. b) Probability distributions of down ‘2,1’ in the two types of
AmPs are the same, while those of bold ‘1,2’ in nonequal AmPs are the sum of
those of up ‘1,2’ and equal ‘1,1’ in equal-value AmPs.

Figure 8 shows that the distribution of false up (bold ‘1,2’) in nonequal AmPs
is a combination of real up (denoted by ‘1,2’) and equal (denoted by ‘1,1’) in
equal-value AmPs because equal values are replaced by up patterns, and that of
equal values is the DES of sleep EEGs in Fig. 3. The distributions of downs
(denoted by ‘2,1’) are the same in the two types of ordinal patterns. In
nonequal AmP, the probability distribution of up increases with the sleep
depth, while that of down presents the opposite trend. In equal-value AmP, the
probability distributions of up and down are significantly close, and both
decreased as the sleep depth increased. Equal-value AmPs further confirm our
conclusion, as shown in Fig. 3, that fluctuations in the EEG amplitude
decreased as the sleep depth increased. Owing to the decrease in amplitude
fluctuations in sleep EEGs, the probabilities of both ups and downs in sleep
EEGs decrease, while those of equal values (i.e., DES) significantly increase
with the sleep depth.

The probabilistic difference in AmPs and different ways to use them directly
influenced the different outcomes of pTIR, pTAS, and PEn in sleep EEGs.
According to Eq. (11), numerical simulations suggest that the probabilistic
difference between false up (‘1,2’) and down (‘2,1’) in nonequal AmPs, i.e.,
(‘1,2’+‘1,1’)/‘2,1’, increased with the sleep depth; hence, noeTAS in Fig. 6
also increased with the sleep depth. In the backward EEG signal, equal values
were transformed into up ‘1,2’. Therefore, noeTIR is measured using the sum of
forward-backward probabilistic differences between up and down, i.e.,
0.5*[(‘1,2’+‘1,1’)/(‘2,1’+‘1,1’)]+‘2,1’/‘1,2’, exhibiting a decreasing trend
in sleep EEGs, as shown in Fig. 6. However, neither noeTIR nor noeTAS truly
quantified the nonequilibrium characteristic of sleep EEGs. Meanwhile, in
equal-value AmPs, pTIR represented the probabilistic difference between up and
down because equality indicated time reversibility and temporal symmetry,
while PEn measured the probabilistic difference between up, down, and
equality. Numerical results indicated that the probabilistic difference in up
and down decreased as the sleep depth increased, while the average amount of
information contained in up, down, and equality increased. Therefore, the pTIR
in Fig. 5, i.e., up-down probabilistic differences, decreased with the
increasing sleep depth, while the PEn of EEGs in Fig. 7 increased with the
sleep depth. Hence, the probability distribution of sleep EEG permutation and
different usage of probabilistic difference of pTIR, pTAS, and PEn are
responsible for the conflicting results in sleep EEGs.

In ordinal patterns and permutation analysis, elements are reordered in the
ascending order; however, it remains unclear how these methods perform of a
series reordered in the descending order. In ordinal patterns, if equal values
are organized according to their order of occurrence, they will be treated as
down in nonequal permutations in the descending order. If indexes of equal
values are modified to be the same in their corresponding groups, the AmP and
OrP have the same relationship as that shown in Fig. 1; moreover, the results
for TIR, TAS, and PEn show no difference in the ascending order. Otherwise, if
the indexes of equal values are not modified, these permutation methods might
still generate incorrect results as the ascending order. Particularly, noeTAS
yields different outcomes in the ascending and descending orders. The results
of noeTAS ($m$=2 and $\tau$=1) for sleep EEGs are shown in Fig. 9.

Figure 9: Different noeTAS (mean$\pm$standard error) in wake, SI-SIII, and REM
EEG data in ascending and descending ordinal patterns.

In the ascending order, noeTAS is the probabilistic difference between
‘$up_{p}+equal_{p}$’ and ‘$down_{p}$’, while in the descending order, noeTAS
is measured by the probabilistic difference between ‘$down_{p}+equal_{p}$’ and
‘$up_{p}$’. The two noeTAS share consistent results, i.e., their value
increases as the sleep stage progresses; meanwhile, they are both incorrect.

Ordinal pattern is a simplified alternative to the vector; their construction
is a coarse-grained procedure to signals. Weak noises and artifacts are
eliminated during the generation of permutations; therefore, permutation
analysis exhibits evident noise insensitivity Bandt and Pompe (2002); Bandt
(2016, 2020). To further test the pTIR, Gaussian noises are added and some
frequency components are removed from the sleep EEG data of the five groups.
These procedures through software (e.g., MATLAB in this work) eliminate equal
values and the information conveyed by DES Yao  _et al._ (2023b) from sleep
EEG data. Therefore, no difference was observed between TIR and TAS based the
ordinal pattern, regardless of the ascending or descending order. Taking $m$=2
and $\tau$=1 as an example, when Gaussian noises with signal noise ratio (SNR)
from 0-0.01 dB are added, the results are not significantly affected that the
pTIR values of EEG data decrease as the sleep stage progresses. Alternatively,
when a 0.3–35 Hz bandpass filter is implemented to the sleep EEG data, pTIR
also shows a decreased trend with the increase of sleep stages, confirming the
robustness of permutation analysis Bandt and Pompe (2002); Bandt (2016, 2020).

According to the research on pTIR for sleep classification, pTIR effectively
quantifies nonequilibrium characteristics in EEGs, yielding more reliable
results than the entropy measure. Owing to the important influence of equal
values on ordinal patterns, equal-value permutation is necessary in the pTIR
analysis of sleep EEGs. The TIR, TAS, and PEn based on equal-value and
nonequal permutations yielded different or even opposite results in the 5
groups of sleep EEG signals; however, this helped us in gaining a deeper
understanding of the differences between analytical methods, important role of
equal values, and multifaceted characteristics of complex physiological
signals.

##  IV Discussions

In the research on pTIR in sleep EEGs, several issues should be further
discussed.

Equal values in time series and their effect on permutation analysis should be
given sufficient attention. Equal values are generally ignored considering the
limitation resolution of ADC, particularly coarse-grain quantization. In
traditional signal processing theory, an arbitrary time series with a weak
stationarity exhibits a continuous distribution; therefore, equal values are
rare and can be broken numerically by adding small random perturbations Bandt
and Pompe (2002). The traditional treatment to equality is based on the
assumptions that equal values have no significant effect on signal processing
and they do not contain information about systems. These assumptions are not
correct. Equal values significantly impact permutation analysis. Moreover,
they are necessary for the comprehensive construction of ordinal patterns Yao
_et al._ (2022). As shown in Fig. 1, irrespective of the OrP or AmP, the
structural information of vectors can be fully displayed only when equal
values are considered. The distribution of equal values in some signals
conveys important information and has significant effects on the construction
and probability of ordinal patterns, thus yielding different or even
contradictory results, such as the results presented in Figs. 5 and 6.
According to our previous report, PEn and pTIR exhibited contradictory results
for PhysioNet heartbeats data with and without equal values Yao  _et al._
(2020b). In time reversibility, vectors with equal values might be self-
symmetric, such as those in Figs. 1 and 2, and have definitive physical
implication, i.e., time reversibility and temporal symmetry Yao  _et al._
(2022, 2023a). Under acceptable ADC resolution, the DES are an effective
parameter for quantifying signal amplitude fluctuations in extreme forms,
i.e., zero fluctuation. The advantages of DES have been observed in the
characterization of sleep and epileptic EEGs Yao  _et al._ (2021b, 2023b);
they also exhibit reliable performance in heartrate data (derived from
electrocardiography) to reflect the decrease in heart rate variability with
age and heart failure Yao  _et al._ (2019b). If the ADC has a high resolution
or the data undergo preprocessing (e.g., software filtering), equal states are
rarely observed and amplitude fluctuation information cannon be detected; in
this case, a low-pass threshold can be established to filter the differential
states for measuring the amplitude fluctuation Yao  _et al._ (2023b).
Otherwise, the equalities can be increased by further coarse graining the
digital signal, such as the partition symbolic transformation that resembles
the ADC. It should be noted that the results of DES were consistent with those
of pTIR for sleep EEGs. The amplitude fluctuation of signals is the most
direct representation of system information, which is influenced by various
factors such as frequency composition and dynamic characteristics. According
to Figs. 3 and 8, the DES of EEG signals increased with the sleep depth,
indicating that the amplitude fluctuation decreased. In the previous analysis
of heartrate signals, amplitude fluctuations in heartrate decreased with age
and heart failure, consistent with the theory of loss of complexity Yao  _et
al._ (2019b). Our results further confirmed the positive correlation between
amplitude fluctuations and time irreversible features, and whether there exist
other factors closely related to amplitude fluctuations require further
research. Overall, the distribution of equal values conveyed important
information about amplitude fluctuation and significantly affected signal
processing; moreover, it is required to construct reliable ordinal patterns
and has explicit physical meaning in TIR; hence, equal values cannon be
ignored.

Time reversibility and temporal symmetry are equivalent in statistical
definitions, even though they may generate different results in real-world
quantification. In pTIR and pTAS, the construction of ordinal patterns,
particularly the treatment of equal values, plays an important role. In
traditional ordinal patterns, equal values are generally ranked according to
their order of emergence. Considering double values, there are three kinds of
permutations, namely up, down, and equal. In the traditional ordinal, equal
values are neglected and treated as up. As shown in Figs. 2 and 8, the
probabilistic difference between ’up+equal’ and down is calculated in noeTAS,
while that between ‘up+equal’ and ‘down+equal’ is measures in noeTIR.
Numerical calculations demonstrated that the difference in forward-backward
permutations was closer to the pTIR than that in symmetric permutations, but
both were wrong. The noeTIR and noeTAS were both incorrect and not the same.
Taking an extreme example, for a series of all-equal values, the pTAS is 1,
indicating the difference between all up and zero down, while the noeTIR is 0
for the forward and backward series (it is the same as the pTIR of 1). If more
values are considered, the situation will become more complex.

Next, the relationship between PEn and pTIR requires more discussion. Test
results indicated that pTIR enables the more effective classification of sleep
stages than PEn, consistent with the results of our previous report on
epilepsy EEGs Yao  _et al._ (2020a). TIR and Shannon entropy are both
statistical parameters for measuring probabilistic differences; they are both
widely employed in complex process analysis. Shannon entropy quantifies the
static complexity and unpredictability considering the amount of information,
i.e., the mean logarithmic calculation of all probabilities of permutations.
TIR measures nonequilibrium features considering the sum of probabilistic
differences between vectors in forward-backward series or pairs of symmetric
vectors. Mathematically, if two permutations have larger probabilistic
differences, they convey less information while being more nonequilibrium,
thus leading to smaller PEn and bigger pTIR. In the special case where all
permutations of symmetric vectors have the same probability distribution, PEn
reaches a maximum value while pTIR is 0. In another special case where all
permutations are single permutations, pTIR is the maximum 1 while PEn varies
with the difference among probability distributions. This is the fundamental
reason for different or even contradictory results in the same process.
Similar results have also been reported in heartrate analysis Yao  _et al._
(2020b); irrespective of nonequal or equal-value permutation, PEn and pTIR
yielded contradictory results. Such discrepant results of pTIR and PEn
inspired us to gain a more comprehensive and profound understanding of the
characteristics of complex systems from different perspectives.

Forbidden permutation is an important influencing factor that is generally
overlooked in the pTIR analysis of real-world signals. The existence of
forbidden permutations might generate individual permutations, which have no
symmetric form in forward series and de not exist in backward series;
moreover, their probabilistic difference is zero or infinite considering
division-based parameters, which is not appropriate in quantitative TIR.
Associated with forbidden permutations, individual permutations also convey
information about sleep EEGs. In Fig. 4, the DIPs generally increased with the
sleep stage when the data length was larger than 20 s. This may be because as
the sleep depth increased, amplitude fluctuations, type of temporal structure,
and number of ordinal patterns decreased, yielding more forbidden as well as
individual permutations. Similar associations have also been reported in
epileptic EEG analysis Yao  _et al._ (2020a, 2021b). Seizure ictal EEGs
exhibit abnormally large amplitude fluctuations and TIR owing to the
development of synchronous neuronal firings, while seizure-free postictal
brain activity exhibits rather smooth amplitude fluctuations as well as
smaller TIR Yao  _et al._ (2021b). Consistently, brain electric signals under
ictal and postictal states demonstrate larger and smaller DIP values,
respectively Yao  _et al._ (2020a). In a series of related studies, Amigo et
al. found that the existence of forbidden patterns is a feature of chaotic
dynamics and can be used to distinguish random from pseudorandom orbit
generation Amigo  _et al._ (2006); subsequently, they identified false
forbidden patterns Amigo  _et al._ (2007) and detected determinism in noisy
time series based on the properties of topological PEn Amigo  _et al._ (2008).
The decay rate of forbidden permutations Carpi  _et al._ (2010) in stochastic
processes has been reported to be associated with their correlation
structures; moreover, nonlinearity in time series can be possibly detected by
the number of forbidden permutations Kulp  _et al._ (2017); Cuesta Frau
(2020). Given the systematic information conveyed by forbidden permutations
and association of DIPs with sleep stages, individual permutations might also
be potentially used for feature detection in complex systems.

##  V Conclusions

In this study, pTIR in sleep EEG particular the dependence on sleep stage and
the effect of equal EEG values are analyzed. Main findings are summarized
below:

Permutation TIR is an important measure for the quantification of
nonequilibrium EEGs. When symmetric vector differences are used for real-time
requirements, AmP is more suitable as a direct alternative to vector.
Moreover, equal-value ordinal patterns are required because they construct
comprehensive vector structures and self-symmetric vectors convey an important
physical implication, i.e., time reversibility.

Brain electrical activity demonstrates nonequilibrium features that are
influenced by sleep conditions. EEGs during wakefulness exhibit higher pTIR;
during sleep and with the advancement of sleep stages, the pTIR values of EEGs
significantly decrease. These findings suggest that when people fall asleep,
their brain electric activity contains less nonequilibrium characteristics.
The effective classification of sleep EEGs suggested that pTIR could serve as
an aid for the expert manual annotation of sleep stages.

When constructing ordinal patterns, if equal values are ordered according to
their order of occurrence, the values of noeTIR and noeTAS may be
inconsistent, and even if noeTIR is closer to the actual result, they would
both be wrong. Moreover, noeTAS performs differently over series reordered in
the ascending and descending orders. Therefore, the consideration of equal
values is necessary to construct reliable permutations, and it is important to
modify the indexes of equal values to same forms in permutation TIR.

Equal values and individual permutations affect the pTIR, but both DES and DIP
contain important information about sleep EEGs. The DES characterize amplitude
fluctuations in sleep EEGs such that as the sleep stages advance, DES values
significantly increase as amplitude fluctuations decrease. The DIP may be
related to structures that require further investigation.

Comparative analysis of pTIR and PEn as well as numerical simulations of
permutation probability distributions verified the different and even
contradictory results of time reversibility and entropy complexity, thus
providing us with valuable insights on statistical measures and enabling us to
explore complex physiological signals more comprehensively.

##  VI Acknowledgment

The project is supported by the Natural Science Foundation of Jiangsu Province
(Grant No.BK20220383), Natural Science Research of Jiangsu Higher Education
Institutions of China (Grant No.22KJB110003), Natural Science Foundation of
Nanjing University of Posts and Telecommunications (Grant Nos. NY221142,
NY222172), Shanghai Municipal Science and Technology, China Major Project
(Grant No. 2018SHZDZX01), Key Laboratory of Computational Neuroscience and
Brain-Inspired Intelligence (LCNBI) and ZJLab.

## References

  * Andrea  _et al._ (2018) A. K. Andrea, B. Misic, and O. Sporns, Nature Reviews Neuroscience 19, 17 (2018).
  * Bashan  _et al._ (2012) A. Bashan, R. P. Bartsch, J. W. Kantelhardt, S. Havlin, and P. C. Ivanov, Nature Communications 3, 702 (2012).
  * Zhao  _et al._ (2019) D. Zhao, Y. Wang, Q. Wang, and X. Wang, Computer methods and programs in biomedicine 175, 53 (2019).
  * Phan and Mikkelsen (2022) H. Phan and K. Mikkelsen, Physiological Measurement 43, 04TR01 (2022).
  * Ma  _et al._ (2018) Y. Ma, W. Shi, C.-K. Peng, and A. C. Yang, Sleep medicine reviews 37, 85 (2018).
  * Xiong  _et al._ (2017) W. Xiong, L. Faes, and P. C. Ivanov, Physical Review E 95, 062114 (2017).
  * Hou  _et al._ (2021) F. Hou, L. Zhang, B. Qin, G. Gaggioni, X. Liu, and G. Vandewalle, Sleep 44 (2021).
  * Bandt (2017) C. Bandt, Entropy 19, 197 (2017).
  * Liang  _et al._ (2021) X. Liang, J. Xiong, Z. Cao, X. Wang, J. Li, and C. Liu, Physiological Measurement 42, 044001 (2021).
  * Wang  _et al._ (2021) Y.-H. Wang, I.-Y. Chen, H. Chiueh, and S.-F. Liang, IEEE Transactions on Instrumentation and Measurement 70, 1 (2021).
  * Miskovic  _et al._ (2019) V. Miskovic, K. J. MacDonald, L. J. Rhodes, and K. A. Cote, Human brain mapping 40, 538 (2019).
  * Fang  _et al._ (2019) X. Fang, K. Kruse, T. Lu, and J. Wang, Reviews of Modern Physics 91, 045004 (2019).
  * Lynn  _et al._ (2021) C. W. Lynn, E. J. Cornblath, L. Papadopoulos, M. A. Bertolero, and D. S. Bassett, Proceedings of the National Academy of Sciences 118 (2021).
  * Yao  _et al._ (2020a) W. P. Yao, J. Dai, M. Perc, J. Wang, D. Yao, and D. Guo, Nonlinear Dynamics 100, 907 (2020a).
  * Yao  _et al._ (2021a) W. P. Yao, J. Wang, M. Perc, W. L. Yao, J. Dai, D. Guo, and D. Yao, Communications in Nonlinear Science and Numerical Simulation 96, 105688 (2021a).
  * Martinez  _et al._ (2018) J. H. Martinez, J. L. Herrera-Diestra, and M. Chavez, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 123111 (2018).
  * Martin Gonzalo  _et al._ (2019) J. A. Martin Gonzalo, I. Pulido Valdeolivas, Y. Wang, T. Wang, G. Chiclana Actis, M. d. C. Algarra Lucas, I. Palmí Cortés, J. Fernández Travieso, M. D. Torrecillas Narváez, A. A. Miralles Martinez, E. Rausel, D. Gómez Andrés, and M. Zanin, Entropy 21, 868 (2019).
  * Zanin (2021) M. Zanin, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 103118 (2021).
  * Weiss (1975) G. Weiss, Journal of Applied Probability 12, 831 (1975).
  * Kelly (1979) F. P. Kelly,  _Reversibility and stochastic networks_ (Cambridge University Press, Chichester, 1979).
  * Ramsey and Rothman (1995) J. B. Ramsey and P. Rothman, Journal of Money Credit & Banking 28, 1 (1995).
  * Costa  _et al._ (2005) M. D. Costa, A. L. Goldberger, and C. K. Peng, Physical Review Letters 95, 198102 (2005).
  * Costa  _et al._ (2008) M. D. Costa, C. K. Peng, and A. L. Goldberger, Cardiovascular Engineering 8, 88 (2008).
  * Porta  _et al._ (2008) A. Porta, K. R. Casali, A. G. Casali, T. Gnecchi-Ruscone, E. Tobaldini, N. Montano, S. Lange, D. Geue, D. Cysarz, and P. Van Leeuwen, American Journal of Physiology Regulatory Integrative & Comparative Physiology 295, 550 (2008).
  * Guzik  _et al._ (2006) P. Guzik, J. Piskorski, T. Krauze, A. Wykretowicz, and H. Wysocki, Biomedizinische Technik Biomedical Engineering 51, 272 (2006).
  * Ehlers  _et al._ (1998) C. L. Ehlers, J. Havstad, D. Prichard, and J. Theiler, Journal of Neuroscience the Official Journal of the Society for Neuroscience 18, 7474 (1998).
  * Lacasa  _et al._ (2012) L. Lacasa, A. Nunez, E. Roldan, J. M. R. Parrondo, and B. Luque, European Physical Journal B 85, 217 (2012).
  * Lacasa and Flanagan (2015) L. Lacasa and R. Flanagan, Physical Review E 92, 022817 (2015).
  * Flanagan and Lacasa (2016) R. Flanagan and L. Lacasa, Physics Letters A 380, 1689 (2016).
  * Lacasa  _et al._ (2008) L. Lacasa, B. Luque, F. Ballesteros, L. Jordi, and J. C. Nuno, Proceedings of the National Academy of Sciences of the United States of America 105, 4972 (2008).
  * Donges  _et al._ (2013) J. F. Donges, R. V. Donner, and J. Kurths, Europhysics Letters 102, 381 (2013).
  * Zou  _et al._ (2019) Y. Zou, R. V. Donner, N. Marwan, J. F. Donges, and J. Kurths, Physics Reports 787, 1 (2019).
  * Zeng  _et al._ (2013) H. Zeng, M. Alava, E. Aurell, J. Hertz, and Y. Roudi, Physical Review Letters 110, 210601 (2013).
  * Zeng  _et al._ (2011) H. Zeng, E. Aurell, M. Alava, and H. Mahmoudi, Physical Review E 83, 041135 (2011).
  * Daw  _et al._ (2003) C. S. Daw, C. E. A. Finney, and E. R. Tracy, Review of Scientific Instruments 74, 915 (2003).
  * Cammarota and Rogora (2007) C. Cammarota and E. Rogora, Chaos, Solitons & Fractals 32, 1649 (2007).
  * Bandt and Pompe (2002) C. Bandt and B. Pompe, Physical Review Letters 88, 174102 (2002).
  * Bandt (2016) C. Bandt, “Permutation entropy and order patterns in long time series,” in  _Time Series Analysis and Forecasting_ (Springer, 2016) pp. 61–73.
  * Bandt (2020) C. Bandt, Statistical Papers 61, 1565 (2020).
  * Yao  _et al._ (2022) W. Yao, W. Yao, and J. Wang, Physics Letters A 430, 127977 (2022).
  * Yao  _et al._ (2019a) W. Yao, W. Yao, J. Wang, and J. Dai, Physics Letters A 383, 738 (2019a).
  * Yao  _et al._ (2023a) W. Yao, W. Yao, R. Xu, and J. Wang, Communications in Nonlinear Science and Numerical Simulation 117, 106925 (2023a).
  * Amigo  _et al._ (2006) J. M. Amigo, L. Kocarev, and J. Szczepanski, Physics Letters A 355, 27 (2006).
  * Amigo  _et al._ (2007) J. M. Amigo, S. Zambrano, and M. A. Sanjuán, Europhysics Letters 79, 50001 (2007).
  * Amigo  _et al._ (2008) J. M. Amigo, S. Zambrano, and M. A. Sanjuán, Europhysics Letters 83, 60005 (2008).
  * Zanin (2008) M. Zanin, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 013119 (2008).
  * Carpi  _et al._ (2010) L. C. Carpi, P. M. Saco, and O. Rosso, Physica A: Statistical Mechanics and its Applications 389, 2020 (2010).
  * Kulp  _et al._ (2017) C. W. Kulp, L. Zunino, T. Osborne, and B. Zawadzki, Physical Review E 96, 022218 (2017).
  * Cuesta Frau (2020) D. Cuesta Frau, Entropy 22, 494 (2020).
  * Bian  _et al._ (2012) C. Bian, C. Qin, Q. D. Ma, and Q. Shen, Physical Review E 85, 021906 (2012).
  * Zunino  _et al._ (2017) L. Zunino, F. Olivares, F. Scholkmann, and O. A. Rosso, Physics Letters A 381, 1883 (2017).
  * Cuesta Frau  _et al._ (2018) D. Cuesta Frau, M. Varela Entrecanales, A. Molina Pico, and B. Vargas, Complexity 2018, 1324696 (2018).
  * Yao  _et al._ (2019b) W. Yao, W. Yao, and J. Wang, Physics Letters A 383, 1764 (2019b).
  * Yao  _et al._ (2020b) W. Yao, W. Yao, D. Yao, D. Guo, and J. Wang, Applied Physics Letters 116, 014101 (2020b).
  * Yao  _et al._ (2021b) W. Yao, W. Yao, Y. Ju, Y. Xia, D. Guo, and D. Yao, Biomedical Signal Processing and Control 69, 102738 (2021b).
  * Yao  _et al._ (2023b) W. Yao, W. Yao, and J. Wang, Physiological Measurement 44, 095004 (2023b).
  * Goldberger  _et al._ (2000) A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C. Peng, and H. E. Stanley, Circulation 101, e215 (2000).
  * Rechtschaffen and Kales (1968) A. Rechtschaffen and A. Kales,  _A manual of standardized terminology, techniques and scoring system for sleep stages of human subjects_ (U.S. Government Printing Office, National Institute of Health Publication, Washington, D.C.,, 1968).
  * Ichimaru and Moody (1999) Y. Ichimaru and G. B. Moody, Psychiatry and Clinical Neurosciences 53, 175 (1999).


11institutetext: Stanford University, Stanford CA, USA 22institutetext:
Philips Research North America, Cambridge MA, USA 33institutetext: Global
Health Laboratories, Bellevue, WA, USA 44institutetext: Oregon Health &
Science University, Portland, OR, USA 55institutetext: Columbia University
Medical Center, New York, NY, USA 66institutetext: New York University, New
York, NY, USA 77institutetext: Brooke Army Medical Center, San Antonio, TX,
USA 88institutetext: Tripler Army Medical Center, Honolulu, HI, USA
99institutetext: MedStar Washington Hospital Center, Washington, DC, USA
1010institutetext: University of Kentucky, Lexington, KY, USA
1111institutetext: Warren Alpert Medical School of Brown University,
Providence, RI, USA 1212institutetext: University of Vermont Larner College of
Medicine, Burlington, VT, USA

# Weakly Semi-Supervised Detection in Lung Ultrasound Videos

Jiahong Ouyang Work completed during internship at Philips Research North
America.11 Li Chen Corresponding author. Email<EMAIL_ADDRESS>Gary
Y. Li 22 Naveen Balaraju 22 Shubham Patil 22 Courosh Mehanian 33 Sourabh
Kulhare 33 Rachel Millin 33 Kenton W. Gregory 44 Cynthia R. Gregory 44 Meihua
Zhu 44 David O. Kessler 55 Laurie Malia 55 Almaz Dessie 55 Joni Rabiner 55 Di
Coneybeare 55 Bo Shopsin 66 Andrew Hersh 77 Cristian Madar 88 Jeffrey Shupp 99
Laura S. Johnson 99 Jacob Avila 1010 Kristin Dwyer 1111 Peter Weimersheimer
1212 Balasundar Raju 22 Jochen Kruecker 22 Alvin Chen 22

###### Abstract

Frame-by-frame annotation of bounding boxes by clinical experts is often
required to train fully supervised object detection models on medical video
data. We propose a method for improving object detection in medical videos
through weak supervision from video-level labels. More concretely, we
aggregate individual detection predictions into video-level predictions and
extend a teacher-student training strategy to provide additional supervision
via a video-level loss. We also introduce improvements to the underlying
teacher-student framework, including methods to improve the quality of pseudo-
labels based on weak supervision and adaptive schemes to optimize knowledge
transfer between the student and teacher networks. We apply this approach to
the clinically important task of detecting lung consolidations (seen in
respiratory infections such as COVID-19 pneumonia) in medical ultrasound
videos. Experiments reveal that our framework improves detection accuracy and
robustness compared to baseline semi-supervised models, and improves
efficiency in data and annotation usage.

###### Keywords:

Weakly Supervised Learning Semi-Supervised Learning Object Detection Medical
Ultrasound.

##  1 Introduction

Despite the remarkable performance of deep learning networks for object
detection and other computer vision tasks [10, 19, 2], most models rely on
large-scale annotated training examples, which are often unavailable or
burdensome to generate in the medical imaging domain. This is especially true
for video-based imaging modalities such as medical ultrasound, where frame-by-
frame annotation of bounding boxes or other localization labels is extremely
time-consuming and costly, and even more so if annotations must be done by
clinical experts.

Reducing the annotation burden of training object detectors on medical images
has been a focus of much recent work. Semi-supervised and weakly supervised
approaches have been proposed to address the annotation challenge, where
unlabeled or inexactly/inaccurately labeled data are used to supplement
training, often in combination with a small amount of fully labeled data [13,
9, 11, 15, 21, 1]. Examples of weak supervision for object detection include
point annotations [22, 4, 17, 14, 3] and image-level class labels [17, 14,
12], both of which are applied on individual image frames. However, even these
methods of weak supervision may not be practical in the video domain, where
hundreds or thousands of image frames requiring interpretation may be
collected in a single clinical exam.

In this work, we propose a weakly semi-supervised framework for training
object detection models based on video-level supervision, where only a single
label is provided for each video. Video-level labels represent a significantly
weaker form of supervision than instance- or image-level labels but can be
generated much more efficiently. Our approach extends teacher-student models
adopted for semi-supervised object detection [11] to the weakly semi-
supervised video-based detection task. Our main contributions are as follows:

  1. 1.

We introduce a simple mechanism during teacher-student training which
aggregates individual detections from the teacher (pseudo-labels) into video-
level confidence predictions. This allows video-level weak supervision using
any standard classification loss.

  2. 2.

We improve the reliability of pseudo-labels generated during the mutual
learning stage by introducing techniques to re-weigh pseudo-labels based on
video-level weak supervision.

  3. 3.

We investigate the learning dynamics between the teacher and student, and
propose several improvements to the underlying teacher-student mechanism to
increase training stability. These include a method to better initialize
models in the "burn-in" stage and a set of adaptive updating schemes to
optimize knowledge transfer bidirectionally during mutual learning.

We demonstrate the effectiveness of our approach on the task of detecting lung
consolidations in medical ultrasound videos, which is an important step in
aiding diagnosis and management of patients with respiratory infections such
as bacterial and viral pneumonia (including COVID-19 infection). Computer-
aided detection of lung consolidation in medical ultrasound is a uniquely
challenging problem where the appearance of pathology varies dramatically
across disease types, patient populations, and training levels of the
personnel acquiring images. Experimental results on a large, multi-center,
clinical ultrasound dataset for lung consolidation demonstrate that the
proposed framework leads to improved detection performance, robustness, and
training stability compared to existing semi-supervised methods.

##  2 Related Work

Semi-Supervised Object Detection: Semi-supervised object detection aims to
utilize large amounts of unlabeled data together with a small set of labeled
data. These efforts generally fall into two categories: (1) consistency
regularization, which regularizes the prediction of the detector for images
undergoing different augmentations [6, 7], and (2) pseudo-labeling, where a
teacher model is trained on labeled data to generate pseudo-labels for
unlabeled data, and a student model is then trained on both the labeled and
pseudo-labeled data [11, 23, 20, 16, 18]. Unbiased Teacher (UBT) [11] is one
of the state-of-the-art methods in this category. Our work is inspired by the
framework of UBT, but extends the method to the weakly semi-supervised
scenario, where we leverage frame-level pseudo-labels and video-level weak
supervision.

Weakly Semi-supervised Object Detection: Weakly semi-supervised object
detection is usually based on instance-level weak supervision, e.g., a point
on the object [22, 4, 17, 14, 3, 8], or image-level supervision, e.g., the
class of the image [17, 14, 12]. Video-level labels are a significantly weaker
supervisory signal than instance- or frame-level labels, since the only
information provided is that the object class exists somewhere in at least one
frame of the video, but in which specific frame(s) and at what location(s) is
unknown.

##  3 Method

The basic intuition behind our approach is to simultaneously learn from both
frame-level and video-level labels to improve object detection performance. We
denote a fully labeled set $\mathcal{D}_{f}$ and a weakly labeled set
$\mathcal{D}_{w}$ as our training data. The fully labeled set
$\mathcal{D}_{f}=\\{x^{f}_{i},y^{f}_{i}\\}_{i=1}^{N_{f}}$ comprises a set of
$N_{f}$ frames $x^{f}_{i}$ and their paired frame-level annotations
$y^{f}_{i}$ (i.e., the coordinates of all bounding boxes present in each
frame). The weakly labeled set
$\mathcal{D}_{w}=\\{x^{w}_{j,1:T_{j}},z^{w}_{j}\\}_{j=1}^{N_{w}}$ consists of
$N_{w}$ videos $x^{w}_{j,1:T_{j}}$ and their video-level class labels
$z^{w}_{j}$ (indicating whether a video contains at least one instance of the
object class, which is annotated on the whole video of $T_{j}$ frames). Due to
its high accuracy, computational efficiency, and capacity for real-time
inference, YOLO-v5 [5] is adopted as the backbone detector. Non-maximum
suppression (NMS) is applied on the outputs of YOLO-v5 as the final output of
the detector to remove duplicate prediction. Here, we use $c$ to denote the
confidence vector (a part of YOLO-v5’s output) of all predicted boxes from a
frame $x$.

###  3.1 Teacher-Student Training

We adopt the UBT framework from [11] for the weakly semi-supervised detection
task. UBT uses two training stages: a burn-in stage for model initialization
and a mutual learning stage for teacher-student training (Fig. 1).

Figure 1: Overview of the proposed method.

Burn-In Stage: The burn-in stage aims to initialize the detector in a
supervised manner using the fully labeled data $\mathcal{D}_{f}$. The original
burn-in method in UBT adopts model weights after a fixed number of early
training epochs. However, we observed that detection performance can vary
dramatically during training. Here, we improve the training stability during
burn-in by applying hierarchical exponential moving average (EMA) updates at
each iteration and epoch (yellow block in Fig. 1). Specifically, the
iteration-based model $\theta_{I}$ is introduced to transfer knowledge from an
initial detection model $\theta$ after each iteration (that is, the weights of
$\theta_{I}$ are updated per batch). The epoch-based model $\theta_{E}$ is
added to transfer information from $\theta_{I}$ after each epoch (weights
updated after all batches). This is given by:

$\cases{\theta}_{I}\leftarrow\alpha_{i}\theta_{I}+(1-\alpha_{i})\theta,&\textrm{for
each iteration}\\\
\theta_{E}\leftarrow\alpha_{e}\theta_{E}+(1-\alpha_{e})\theta_{I},\textrm{for
each epoch}\vspace{-5pt}$ (1)

where $\alpha_{i}$ and $\alpha_{e}$ are the iteration- and epoch-based EMA
keep rates, respectively. The EMA keep rates define a trade-off between the
rate of knowledge being transferred from the preceding model versus the
stability of the succeeding model. Major oscillations in model performance are
seen while training $\theta_{I}$ alone, even with a carefully selected
$\alpha_{i}$ [5]. In contrast, the addition of $\theta_{E}$ serves to
stabilize the training and results in a better-initialized detector after
burn-in.

Mutual Learning Stage: The mutual learning stage combines the fully labeled
data $\mathcal{D}_{f}$ and weakly labeled data $\mathcal{D}_{w}$ for teacher-
student training. Both the student and teacher models are initialized from the
last checkpoint of $\theta_{E}$ trained in the burn-in stage. During mutual
learning, the student $\theta_{S}$ is optimized via backpropagation using a
combination of full, weak, and semi-supervised losses (see Sections 3.3 and
3.4), and the teacher $\theta_{T}$ is updated via a gradual EMA transfer of
weights from the student. Analogous to the model updates in the burn-in stage,
the student is updated with an iteration-based EMA during mutual learning,
while the teacher is updated with an epoch-based EMA. An additional mutual
refinement scheme is introduced to adaptively adjust the EMA keep rate of the
teacher, as well as to conditionally allow transfer of weights back to the
student (see Section 3.5). As mutual learning progresses, the accuracy and
stability of pseudo-labels produced by the teacher are continuously improved,
which in turn improves knowledge distillation to and from the student. At the
end of the training, only the teacher $\theta_{T}$ is kept for evaluation and
deployment.

###  3.2 Weakly Semi-Supervised Learning

Frame-Level Full Supervision: For both burn-in and mutual learning, fully
labeled data $\mathcal{D}_{f}$ are used in training, with supervision by a
detection loss $\mathcal{L}_{f\\_sup}$:

$\mathcal{L}_{f\\_sup}=\sum_{i=1}^{N_{f}}\lambda_{coord}\mathcal{L}_{coord}(\mathcal{T}_{s}(x^{f}_{i}),y^{f}_{i})+\lambda_{conf}\mathcal{L}_{conf}(\mathcal{T}_{s}(x^{f}_{i}),y^{f}_{i}),\vspace{-5pt}$
(2)

where $\mathcal{L}_{coord}$ is the bounding box coordinate error, and
$\mathcal{L}_{conf}$ is the binary cross-entropy loss between predicted box
confidences and corresponding box labels. $\lambda_{coord}$ and
$\lambda_{conf}$ balance the two losses. $\mathcal{T}_{s}$ denotes the data
augmentation.

Frame-Level Semi-Supervision: In the mutual learning stage, the weakly labeled
data $\mathcal{D}_{w}$ are added to allow frame-level semi-supervision based
on pseudo-labels from the teacher $\theta_{T}$. Specifically, we generate sub-
clips $\tilde{x}^{w}_{j,1:N_{fpv}}$ by uniformly sampling $N_{fpv}$ frames
from each video $x^{w}_{j,1:T_{j}}$ in the weak dataset. These sub-clips are
fed through the teacher with reduced augmentation $\mathcal{T}_{r}$ to obtain
box predictions $\hat{y}^{T}_{j,1:N_{fpv}}$ with confidence scores
$c^{T}_{j,1:N_{fpv}}$ from each frame. Only predicted boxes with confidence
above a threshold $\beta$ are kept as pseudo-labels. We use a second detection
loss, similar to Eq. 2, to train the student $\theta_{S}$ against the
teacher’s pseudo-labels:

$\mathcal{L}_{f\\_semi}=\sum_{j=1}^{N_{w}}\sum_{t=1}^{N_{fpv}}\left[\lambda_{coord}\mathcal{L}_{coord}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})+\lambda_{conf}\mathcal{L}_{conf}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})\right]\vspace{-5pt}$
(3)

Video-Level Weak Supervision: Finally, we utilize the confidence of boxes
predicted by the student $c^{S}_{j,1:N_{fpv}}$ to obtain a frame-level
confidence by computing the maximum confidence among all detected boxes in the
frame, i.e., $max(c^{S}_{j,t})$. A final video-level prediction
$\tilde{c}^{S}_{j}$ is computed as the averaged frame-level confidence over
the sub-clip, i.e.,
$\tilde{c}^{S}_{j}=\frac{1}{N_{fpv}}\sum_{t=1}^{N_{fpv}}max(c^{S}_{j,t})$.

Here, we apply a video-level binary cross-entropy classification loss to
supervise the student $\theta_{S}$ against the video-level labels $z^{w}_{j}$
from the weak data $\mathcal{D}_{w}$:

$\mathcal{L}_{v\\_weak}=-\sum_{j=1}^{N_{w}}\left[z^{w}_{j}log(\tilde{c}^{S}_{j})+(1-z^{w}_{j})log(1-\tilde{c}^{S}_{j})\right]\vspace{-5pt}$
(4)

Combined Loss: The final loss function for training the student model combines
the fully supervised detection loss $\mathcal{L}_{f\\_sup}$, the frame-level
semi-supervised detection loss $\mathcal{L}_{f\\_semi}$, and the video-level
weakly supervised loss $\mathcal{L}_{v\\_weak}$:

$\mathcal{L}=\lambda_{f\\_sup}\mathcal{L}_{f\\_sup}+\lambda_{f\\_semi}\mathcal{L}_{f\\_semi}+\lambda_{v\\_weak}\mathcal{L}_{v\\_weak}\vspace{-5pt}$
(5)

with $\lambda_{f\\_sup}$, $\lambda_{f\\_semi}$, and $\lambda_{v\\_weak}$ to
balance the three loss components.

###  3.3 Weighted Pseudo-labels

A critical factor in the effectiveness of mutual learning is the quantity and
quality of pseudo-labels from the teacher $\theta_{T}$. We propose two pseudo-
label re-weighting techniques to increase the number of high-quality pseudo-
labels during training.

Weakly Supervised Pseudo-label Filtering: The first approach utilizes the weak
video-level label $z^{w}_{j}$ to filter false pseudo-labels. For a negative
video ($z^{w}_{j}=0$), we can simply remove all pseudo-labels from every
frame, which could be considered as re-weighting the pseudo-label as 0. For a
positive video ($z^{w}_{j}=1$), if no pseudo-label confidence exceeds $\beta$,
we keep the pseudo-label with the highest confidence if it exceeds a lower
threshold $\beta_{l}$, where $\beta_{l}<\beta$, which could be considered as
re-weighting the pseudo-label from 0 to its confidence.

Soft Pseudo-labels: The second approach assigns a weight to each pseudo-label
based on its prediction confidence. That is, we re-weigh the loss component
for each pseudo-label $\hat{y}^{T}_{j,t}$ by the square of its confidence
$(c^{T}_{j,t})^{2}$ to create "soft pseudo-labels" ($\hat{y}^{T}_{j,t}$,
$c^{T}_{j,t}$). The semi-supervised detection loss is reformulated as:

$\aligned\mathcal{L}_{coord}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})&=\sum_{k}^{N_{j,t}}(c_{j,t,k}^{T})^{2}\mathcal{L}_{coord}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})_{k}\\\
\mathcal{L}_{conf}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})=\sum_{k}^{N_{j,t}}(c_{j,t,k}^{T})^{2}\mathcal{L}_{conf}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{T}_{j,t})_{k}\vspace{-5pt}$
(6)

where $N_{j,t}$ is the number of pseudo-labels in frame $t$ of video $j$.
$c^{T}_{j,t,k}$ denotes the confidence of the $k$-th pseudo-label in a given
frame.
$\mathcal{L}_{coord}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{w}_{j,t})_{k}$ and
$\mathcal{L}_{conf}(\mathcal{T}_{s}(x^{w}_{j,t}),\hat{y}^{w}_{j,t})_{k}$ are
the bounding box coordinate and confidence loss components for the $k$-th
pseudo-label.

###  3.4 Bidirectional, Adaptive, Teacher-Student Mutual Refinement

The learning dynamics between the teacher and student also play a significant
role in determining training stability and model robustness. First, the
teacher should be updated at a sufficient rate such that it can catch up to
the student before the student overfits, i.e., the EMA keep rate $\alpha_{e}$
cannot be too large. At the same time, the teacher should be characterized by
gradual changes in the training curve as opposed to rapid oscillations. Thus
the teacher also cannot be made to update too quickly, i.e., $\alpha_{e}$
cannot be too small. Finally, the student’s training curve should not have
sudden drops, for example, due to a bad training batch. Here, we introduce two
additional techniques to dynamically balance the rate and direction of
knowledge transfer during mutual learning:

Adaptive EMA (Student $\rightarrow$ Teacher): When using a fixed EMA keep rate
$\alpha_{e}$, there is a trade-off between training stability and rate of
knowledge transfer from the student. Instead, we propose an adaptive EMA keep
rate that is conditioned on the relative performance of the teacher and
student after each training epoch. We use a sigmoid-shaped function for
$\alpha_{e}$, given by:

$\alpha_{e}=\alpha_{e,min}+(\alpha_{e,max}-\alpha_{e,min})\cdot\frac{1}{1+e^{-\tau_{0}(m_{T}-m_{S})-\tau_{1}}}\vspace{-5pt}$
(7)

where $\alpha_{e,min}$, $\alpha_{e,max}$, $\tau_{0}$ and $\tau_{1}$ are hyper-
parameters defining the function shape. $m_{T}$ and $m_{S}$ denote teacher and
student performance on a validation set according to some evaluation metric.
The adaptive scheme allows $\alpha_{e}$ to be dynamically adjusted, that is,
$\alpha_{e}$ is decreased (higher rate of knowledge transfer) as the student
outperforms the teacher, and increased (lower rate of knowledge transfer) as
the student underperforms compared to the teacher.

Inverse Adaptive EMA (Teacher $\rightarrow$ Student): To avoid sudden drops in
performance by the student during mutual learning, we further introduce a
mechanism which allows knowledge transfer in the reversed direction, i.e.,
from the teacher to the student. We design a similar sigmoidal function for
the inverse EMA keep rate $\alpha_{inv}$, given by:

$\alpha_{inv}=\cases{1},&m_{T}\leq m_{S}\\\
\alpha_{inv,min}+(2-2\alpha_{inv,min})\cdot\frac{1}{1+e^{-\tau_{2}(m_{S}-m_{T})}},m_{T}>m_{S}\vspace{-5pt}$
(8)

where $\alpha_{inv,min}$ and $\tau_{2}$ are hyper-parameters of the function.
Here, knowledge transfer to the student is increased (lower $\alpha_{inv}$)
when the teacher outperforms the student, and decreased (higher
$\alpha_{inv}$) when the teacher underperforms.

##  4 Experiments

###  4.1 Experimental Settings

Data. An extensive retrospective, multi-center clinical dataset of 7,998 lung
ultrasound videos were used in this work. The data were acquired from 420
patients with suspicion of lung consolidation or other related pathology
(e.g., pneumonia, pleural effusion) from 8 U.S. clinical sites between 2017
and 2020. The videos were each at least 3 seconds in length and contained at
least 60 frames. 385 (fully labeled training set), 337 (validation set) and
599 (test set) videos were annotated for lung consolidation regions using
bounding boxes. All data were partitioned at the subject level. The remaining
6,677 videos were annotated only for the presence or absence of lung
consolidation at the video level (weakly labeled set). Annotation was carried
out by a multi-center team of expert physicians with medical training in lung
ultrasound. Each video was annotated by two experts and adjudicated by a third
expert when a disagreement between the first two annotators occurred.

Implementation Details. We used the PyTorch Ultralytics implementation of the
YOLO-v5 object detector [5] with default training settings (Adam optimizer
with learning rate of 0.001). The weights for the confidence and coordinate
losses were set to $\lambda_{conf}=1.0$, $\lambda_{coord}=0.05$. The weights
for the frame-level fully supervised, frame-level semi-supervised, and video-
level weakly supervised losses were set to
$\lambda_{f\\_sup}=\lambda_{f\\_semi}=1$, and $\lambda_{v\\_weak}=0.05$. For
training from pseudo-labels without re-weighting, the confidence threshold was
set to $\beta=0.5$. Otherwise, hyperparameters were set to
$\beta=\beta_{l}=0.1$ when using weighted pseudo-labels. To train with a fixed
EMA keep rate, we used $\alpha_{e}=0.95$. Otherwise, when applying
bidirectional adaptive EMA for mutual teacher-student refinement, the
hyperparameters were set to $\alpha_{e,min}=0.75$, $\alpha_{e,max}=0.99$,
$\alpha_{inv,min}=0.85$, $\tau_{0}=180$, $\tau_{1}=3$, and $\tau_{2}=180$.

Table 1: Validation and test mAP for fully, semi-, and weakly semi-supervised models. Mean $\pm$ standard deviation based on five repeated experiments. All methods were significantly superior to YOLO (first row) and significantly inferior to the proposed method (last row)(paired two-way t-test, p-value $<$ 0.05) Category |  Method |  Validation mAP |  Test mAP
---|---|---|---
Fully supervised |  YOLO [5] |  0.435 $\pm$ 0.012 |  0.412 $\pm$ 0.013
YOLO+HE |  0.452 $\pm$ 0.015 |  0.440 $\pm$ 0.017
Semi-supervised |  YOLO+HE+Unlabeled |  0.468 $\pm$ 0.005 |  0.447 $\pm$ 0.003
Weakly Semi-supervised |  YOLO+HE+Weak |  0.505 $\pm$ 0.004 |  0.476 $\pm$ 0.002
YOLO+HE+Weak+Pseudo |  0.508 $\pm$ 0.004 |  0.479 $\pm$ 0.004
YOLO+HE+Weak+TSMR |  0.515 $\pm$ 0.005 |  0.480 $\pm$ 0.004
YOLO+HE+Weak+Pseudo+TSMR |  0.519 $\pm$ 0.003 |  0.484 $\pm$ 0.003

Experiments. We first trained a baseline YOLO-v5 detector on the fully
supervised data $\mathcal{D}_{f}$, denoted as YOLO. We compared the baseline
to a fully supervised training experiment with hierarchical (iteration and
epoch-based) EMA training during burn-in, as proposed in Section 3.2; this is
denoted by +HE. All subsequent experiments involving teacher-student mutual
learning were initialized from the same YOLO+HE model checkpoint. We
implemented the semi-supervised approach from Unbiased Teacher [11] by adding
all videos from the weakly supervised dataset $\mathcal{D}_{w}$, but without
providing video-level labels (i.e., treating these as unlabeled data); this is
denoted as +Unlabeled. Note, to demonstrate the effectiveness of the proposed
method, +Unlabeled a derived version of the Unbiased Teacher [11], using its
way of utilizing unlabeled data while keeping the same setting of the model
and train strategy as the rest of competing methods. We then introduced our
proposed method of weak semi-supervision, described in Section 3.3, by
including the video-level labels from $\mathcal{D}_{w}$, denoted as +Weak.
Finally, experiments using our methods for pseudo-label re-weighting (Section
3.4) and bidirectional, adaptive, teacher-student mutual refinement (Section
3.5) are denoted as +Pseudo and +TSMR respectively. We compared mean Average
Precision (mAP) on the validation and test sets described above, with each
experiment repeated five times to assess repeatability. Experimental results
are summarized in Table 1.

###  4.2 Results & Discussion

Figure 2: Learning curves of teacher models $\theta_{T}$ on validation set
during teacher-student mutual learning. Solid lines show mean validation mAP
across five repeated experiments. Ranges indicate 95% confidence intervals.

Contribution of Hierarchical EMA Training: The baseline, fully-supervised
detector (YOLO) achieved validation and test mAP of 0.435 and 0.412,
respectively. These improved to 0.452 and 0.440, respectively, with the
inclusion of the hierarchical EMA training strategy during burn-in (YOLO+HE).

Contribution of Semi- and Weak (Video-level) Supervision: The addition of
unlabeled data and semi-supervision based on YOLO+HE+Unlabeled improved
validation mAP from 0.452 to 0.468 and test mAP from 0.440 to 0.447, which was
a statistically significant increase (p-value $<$ 0.05). Furthermore, the
standard deviation of mAP values over repeated experiments decreased (0.015 to
0.005 for validation, 0.017 to 0.003 for test), suggesting that the semi-
supervised model is more stable and repeatable across runs. This was also
reflected in the tighter 95% confidence intervals for validation mAP learning
curves across repeated runs (Fig. 2). Model performance again increased with
the introduction of weak supervision of video-level labels (YOLO+HE+Weak)
(validation mAP 0.505, test mAP 0.476, p-value $<$ 0.05), with corresponding
decreases in mAP standard deviation (0.004 in validation, 0.002 in test) and
95% confidence intervals over repeat runs (Fig. 2). To further investigate the
contribution of video-level supervision, we trained models with all fully
labeled data $\mathcal{D}_{f}$ but utilized a proportion of video labels
$z^{w}_{j}$ with the remainder of $\mathcal{D}_{w}$ treated as unlabeled. mAP
improved consistently with increased video-level supervision (Fig. 3).

Figure 3: Contribution of video-level supervision during teacher-student
mutual learning. Models were trained with varying proportions of video labels
relative to unlabeled data.

Contribution of Weighted Pseudo-labels: The use of weighted pseudo-labels
(YOLO+HE+Weak+Pseudo) further improved validation mAP from 0.505 to 0.508 and
test mAP from 0.476 to 0.479. The re-weighting mechanism eliminated all false
positive pseudo-labels in negative videos and increased the number of pseudo-
labels to better match the overall number of true labels. In comparison, fixed
pseudo-label thresholds resulted in worse detection performance (test mAP
0.477, 0.474, 0.476, and 0.474 for thresholds of 0.1, 0.3, 0.5, and 0.7).

Contribution of Bidirectional Teacher-Student Mutual Refinement: Teacher-
student mutual refinement (YOLO+HE+Weak+TSMR) further boosted validation mAP
from 0.508 to 0.515 and test mAP from 0.479 to 0.480. Model repeatability was
improved, as shown in Fig. 2, where variation between experiments was greatly
reduced (narrow 95% confidence intervals) and model convergence occurred more
quickly. Ablation experiments confirmed that a fixed EMA keep rate
$\alpha_{e}$ was unable to achieve comparable detection performance compared
to the proposed bidirectional adaptive EMA updates ($\alpha_{e}=0.9$, lower
fixed rate: 0.506 and 0.471; $\alpha_{e}=0.95$, moderate fixed rate: 0.505 and
0.476; and $\alpha_{e}=0.99$, high fixed rate: 0.489 and 0.470, for validation
mAP and test mAP respectively).

Final Results for Proposed Method: Finally, best-performing models
incorporating all proposed components achieved validation mAP of 0.519 and
test mAP of 0.484, which were statistically significant improvements to both
fully supervised (0.452 and 0.440) and semi-supervised (0.468 and 0.447)
baselines (p-value $<$ 0.05). Furthermore, an extra experiment suggested that
comparable test mAP of baseline YOLO detector (0.412) could be achieved by the
proposed method using merely one-third of the labeled data (0.395). The
reference speed was 2.9ms per frame on a NVIDIA GeForce RTX 3090 GPU, enabling
real-time detection. Examples of lung consolidation detection in ultrasound
are seen in Fig. 4, where the proposed method demonstrates successful
detection of challenging pathology not identified (or falsely identified) by
the baseline YOLO detector.

##  5 Conclusion

This is the first study to introduce a weakly semi-supervised framework for
object detection on medical video data. Our method extends a teacher-student
training strategy to provide weak supervision via a video-level loss. We also
introduce improvements to the underlying teacher-student mutual learning
mechanism, including methods to improve the quality of pseudo-labels and
optimize knowledge transfer between the student and teacher. Empirical results
on a lung ultrasound pathology detection task demonstrate that the framework
leads to improved detection accuracy and robustness compared to existing
baseline models, while also being more efficient in data and annotation usage.
One limitation of the method is the need to empirically select
hyperparameters, which could be resolved using adaptive hyperparameter tuning
techniques as part of future work. Moreover, we considered the proposed
components as orthogonal research directions as other pseudo-label refinement
techniques in the state-of-the-art semi-supervised detection methods, which
could be further combined in achieving better performance. Lastly, the
proposed improvements to the teacher-student mechanism could potentially be
adapted for other semi-supervised and weakly supervised learning tasks,
including classification and segmentation.

Figure 4: Lung consolidation detection with baseline YOLO (top) and proposed
method (YOLO+HE+Weak+Pseudo+TSMR, bottom). White boxes show expert
annotations. Orange boxes show model detections. mAPs calculated for each
video are also shown.

####  5.0.1 Acknowledgements

We would like to acknowledge the contributions from the following people for
their efforts in data curation and annotations: Zohreh Laverriere, Xinliang
Zheng (Lia), Annie Cao, Katelyn Hostetler, Yuan Zhang, Amber Halse, James
Jones, Jack Lazar, Devjani Das, Tom Kennedy, Lorraine Ng, Penelope Lema, Nick
Avitabile.

## References

  * [1] Bakalo, R., Goldberger, J., Ben-Ari, R.: Weakly and semi supervised detection in medical imaging via deep dual branch net. Neurocomputing 421, 15–25 (2021)
  * [2] Bassiouny, R., Mohamed, A., Umapathy, K., Khan, N.: An Interpretable Object Detection-Based Model For The Diagnosis Of Neonatal Lung Diseases Using Ultrasound Images. In: 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 3029–3034. IEEE (nov 2021). https://doi.org/10.1109/EMBC46164.2021.9630169, https://ieeexplore.ieee.org/document/9630169/
  * [3] Chai, Z., Lin, H., Luo, L., Heng, P.A., Chen, H.: Orf-net: Deep omni-supervised rib fracture detection from chest ct. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 238–248. Springer (2022)
  * [4] Ge, Y., Zhou, Q., Wang, X., Shen, C., Wang, Z., Li, H.: Point-teaching: Weakly semi-supervised object detection with point annotations. arXiv preprint arXiv:2206.00274 (2022)
  * [5] Glenn Jocher, e.a.: ultralytics/yolov5: v6.0 - YOLOv5n ’Nano’ models, Roboflow integration, TensorFlow export, OpenCV DNN support (Oct 2021). https://doi.org/10.5281/zenodo.5563715
  * [6] Jeong, J., Lee, S., Kim, J., Kwak, N.: Consistency-based semi-supervised learning for object detection. Advances in neural information processing systems 32 (2019)
  * [7] Jeong, J., Verma, V., Hyun, M., Kannala, J., Kwak, N.: Interpolation-based semi-supervised learning for object detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 11602–11611 (2021)
  * [8] Ji, H., Liu, H., Li, Y., Xie, J., He, N., Huang, Y., Wei, D., Chen, X., Shen, L., Zheng, Y.: Point beyond class: A benchmark for weakly semi-supervised abnormality localization in chest x-rays. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 249–260. Springer (2022)
  * [9] Jiao, R., Zhang, Y., Ding, L., Cai, R., Zhang, J.: Learning with Limited Annotations : A Survey on Deep Semi-Supervised Learning for Medical Image Segmentation. arXiv pp. 1–19 (2022)
  * [10] Kulhare, S., Zheng, X., Mehanian, C., Gregory, C., Zhu, M., Gregory, K., Xie, H., McAndrew Jones, J., Wilson, B.: Ultrasound-based detection of lung abnormalities using single shot detection convolutional neural networks. In: Simulation, image processing, and ultrasound systems for assisted diagnosis and navigation, pp. 65–73. Springer (2018)
  * [11] Liu, Y.C., Ma, C.Y., He, Z., Kuo, C.W., Chen, K., Zhang, P., Wu, B., Kira, Z., Vajda, P.: Unbiased teacher for semi-supervised object detection. arXiv preprint arXiv:2102.09480 (2021)
  * [12] Meethal, A., Pedersoli, M., Zhu, Z., Romero, F.P., Granger, E.: Semi-weakly supervised object detection by sampling pseudo ground-truth boxes. arXiv preprint arXiv:2204.00147 (2022)
  * [13] Peng, J., Wang, Y.: Medical Image Segmentation with Limited Supervision : A Review of Deep Network Models. arXiv pp. 1–24 (2021)
  * [14] Ren, Z., Yu, Z., Yang, X., Liu, M.Y., Schwing, A.G., Kautz, J.: Ufo2: A unified framework towards omni-supervised object detection. In: European Conference on Computer Vision. pp. 288–313. Springer (2020)
  * [15] Shao, F., Chen, L., Shao, J., Ji, W., Xiao, S., Ye, L., Zhuang, Y., Xiao, J.: Deep learning for weakly-supervised object detection and localization: A survey. Neurocomputing (2022)
  * [16] Tang, Y., Chen, W., Luo, Y., Zhang, Y.: Humble teachers teach better students for semi-supervised object detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 3132–3141 (2021)
  * [17] Wang, P., Cai, Z., Yang, H., Swaminathan, G., Vasconcelos, N., Schiele, B., Soatto, S.: Omni-detr: Omni-supervised object detection with transformers. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 9367–9376 (2022)
  * [18] Wang, Z., Li, Y., Guo, Y., Fang, L., Wang, S.: Data-uncertainty guided multi-phase learning for semi-supervised object detection. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4568–4577 (2021)
  * [19] Xing, W., Li, G., He, C., Huang, Q., Cui, X., Li, Q., Li, W., Chen, J., Ta, D.: Automatic detection of a-line in lung ultrasound images using deep learning and image processing. Medical Physics (2022)
  * [20] Xu, M., Zhang, Z., Hu, H., Wang, J., Wang, L., Wei, F., Bai, X., Liu, Z.: End-to-end semi-supervised object detection with soft teacher. In: Proceedings of the IEEE/CVF International Conference on Computer Vision. pp. 3060–3069 (2021)
  * [21] Zhang, D., Zeng, W., Guo, G., Fang, C., Cheng, L., Han, J.: Weakly supervised semantic segmentation via alternative self-dual teaching. arXiv preprint arXiv:2112.09459 (2021)
  * [22] Zhang, S., Yu, Z., Liu, L., Wang, X., Zhou, A., Chen, K.: Group r-cnn for weakly semi-supervised object detection with points. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 9417–9426 (2022)
  * [23] Zhou, Q., Yu, C., Wang, Z., Qian, Q., Li, H.: Instant-teaching: An end-to-end semi-supervised object detection framework. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 4081–4090 (2021)