\section{Introduction}
\label{intro}

In the review \cite{anderson_surgical_2013}, the authors have identified 7 publications between 1991 and 2008 reporting adverse events (AEs). Over all these studies, the median number of patients having undergone one or several AEs was 14.4\%, and a third (37.9\%) of those AEs were considered preventable. An AE is defined by the World Health Organization (WHO) as ``an injury related to medical management, in contrast to complications of disease. Medical management includes all aspects of care, including diagnosis and treatment, failure to diagnose or treat, and the systems and equipment used to deliver care'' \cite{world_health_organization_who_2005}. In the surgical field, we can distinguish between different events categorized as postoperative adverse events (pAEs) for AEs occurring following surgery, and intraoperative adverse events (iAEs), when AEs occur during surgery.

Hospitals use risk management to prevent AEs. This consist in identifying and characterizing AEs along with their severity, with the aim to propose strategies designed to reduce the likelihood they will occur again. The identification consists in determining when an AE occurred and which anatomic structure was affected. The characterization consists in determining the AE's severity. However, since both steps are performed manually, this is a costly and time-consuming process, prone to both subjectivity and mistakes. In this paper, we have analyzed the relevance of surgical process models (SPMs) to help identification of iAEs.

A SPM is ``a simplified pattern of a surgical process that reflects a predefined subset of interest of the surgical process in a formal or semi-formal representation'' \cite{jannin_modeling_2001}. A SPM describes a surgical procedure at different granularity levels: phases, steps, and activities \cite{lalys_surgical_2013}. A surgical procedure is divided into successive phases corresponding to the procedure's main periods. A phase is composed of one or several steps. A step is a sequence of activities deployed to achieve a surgical objective. An activity is a physical action performed by the surgeon. Each activity is deconstructed into different components, including the action verb, anatomic structure concerned by the action, and surgical instrument employed to perform this action.

Surgical process modeling has been used in various applications, such as surgical skills evaluation \cite{riffaud_recording_2010,forestier_classification_2012}, operating room management optimization \cite{sandberg_deliberate_2005, padoy_-line_2008}, or robotic assistance \cite{ko_surgery_2007,nomm_recognition_2008}. However, SPMs have rarely been applied for surgical quality assessment. A method was presented in \cite{bouarfa_workflow_2012} for detecting modifications from the standard process called deviations. The authors employed surgical tool information  to create a standard surgical process and draw correlations between this standard surgical process and a specific surgery using the Needleman-Wunsch global alignment algorithm. One limitation of this method is that the reasons for the deviations are not identified.

Deviation detection has been studied in other domains, such as bank \cite{jadhav_anomaly_2013} or software security \cite{tan_hidden_2008}. The principle of deviation detection relies on constructing a standard process and detecting deviations using a comparison between this standard process and a new one. To the best of our knowledge, these authors did not distinguish different types of deviations either.

To overcome this limit, we propose for the surgical field, as illustrated in table~\ref{tab:deviations}, three types of surgical deviations based on observations from participating surgeons. It is important to note that the notion of deviation is independent of the occurrence of AE, this notion only reflects the modification of the standard surgical process model. An AE could occur as a result of any type of deviation or even if no deviation is visible.

\begin{table*}[!ht]

	\centering
	\caption{Surgical deviation types and definitions. }
    \begin{tabular}{>{\centering\arraybackslash}m{.20\linewidth}|>{\centering\arraybackslash}m{.70\linewidth}}
         \shortstack{Surgical \\deviation type}   &  Definition  \\\hline
         Context deviations & Deviations due to patient's particularities as anatomic specificities, patient's pathology, and co-morbidity; this category also considers all deviations due to the surgical context, as operating room disruptions. \\\hline
         Expert deviations &  Deviations due to the surgeon who performs the surgery; this category includes deviations due to surgical expert knowledge, and surgeons’ habits or preferences.\\\hline
         Event deviations & Deviations from the usual surgical process to correct or limit the impact of iAEs.
    \end{tabular}
    \label{tab:deviations}
\end{table*}

To help the identification of iAES, this work aims to detect surgical deviations from a standard surgical process and classify them according to the above categories. For this purpose, we propose using an extension to non-linear temporal scaling \cite{forestier_non-linear_2014}, called Multi-Dimensional Non-Linear Temporal Scaling (MD-NLTS), with the aim to detect deviations and a hidden semi-Markov model (HsMM) designed to classify them.

\section{Material and methods}
\label{Method}
This section presents our offline method to detect and classify deviations in rectopexy surgery for skill analysis, as summarized in Fig~\ref{fig:method}.

\begin{figure*}[h]
 \centering
 \includegraphics[width=\textwidth]{method2.png}
 \caption{{\bf Classification of surgical deviation types based on four modules.}
Module A enables the creation of individual Surgical Process Models (iSPMs) based on surgical video annotations. Each color in the iSPMs represents one type of activity. Module B provides two outputs, a standard surgical process $\overline{S}$, and the sequences $S^*_{1..n}$ temporally aligned to $\overline{S}$. Module C compares one aligned sequence to the standard surgical process $\overline{S}$ for each instant to detect a potential deviation. The deviations are highlighted by dots in deviation sequences. Module D classifies each type of surgical deviation, i.e., blue diagonal crosshatch for a context deviation and red horizontal crosshatch for an event deviation.}
 \label{fig:method}
\end{figure*}

Our method is composed of four modules: A) the creation of individual surgical process models (iSPMs) based on clinical data; B) the creation of a standard surgical process; C) the detection of surgical deviations; D) the classification of deviation types. Each of these modules is described in the following subsections.

\subsection{Creation of individual surgical process models based on clinical data (Module A)}

The objective of this module is to describe surgical procedures with individual surgical process models (iSPMs) based on observations of surgical videos.

\subsubsection{Data}
\label{Data}
The dataset used in this paper consists of 11 endoscopic videos of laparoscopic rectopexies. A rectopexy is a digestive surgery that consists of correcting the anal prolapse by fixing the rectum to the sacrum through meshes. The operations were performed by a single expert surgeon at the Grenoble University Hospital, France, involving 11 women who had not undergone a hysterectomy during previous hospitalization.

This study was approved by an ethics committee and declared to the French authorities (CCTIRS\footnote{Comit\'e Consultatif sur le Traitement de l'Information en mati\`ere de Recherche dans le domaine de la Sant\'e.} and CNIL\footnote{Commission Nationale de l'Informatique et des Libert\'es.}). All patients operated by the participating surgeon between January 2015 and December 2017 were included as long as they met the declared inclusion criteria and signed a written informed consent authorizing data collection and data utilization for this study.

Since the surgeon performed the surgeries in a limited period, we considered that his knowledge, habits, and preferences did not vary enough to introduce expert deviations. We considered that the dataset thus contained only context and event deviations. However, some rare variations of performance of an expert may occur due to personal reasons. This was not considered here, but this could be checked in a prospective study.

\subsubsection{Creation of individual surgical process models}
\label{iSPM}
For the creation of iSPMs, we have focused on the two following phases: dissection and resection. The objectives of these phases were to respectively access both fixation points (rectum and sacrum) and remove the Pouch of Douglas. According to the participating surgeons, these phases are the most difficult to perform and most likely to cause AEs. In another hand, these phases are the most standardized according to patients’ particularities. One element able to modify the process is the presence or not of a uterus. In our dataset, no patient underwent a previous hysterectomy.

To create iSPMs, first, a Cognitive Task Analysis (CTA) \cite{clark_cognitive_2008} was conducted by a bio-medical engineer familiarized with this methodology and involved two expert surgeons including the one who performed the surgical procedures. The objective of the CTA was to capture and understand the expert knowledge to allow the annotation of the activities and iAEs. The latter were manually recorded by the bio-medical engineer thanks to the ``Surgery Workflow Toolbox [annotate]" software \cite{garraud_ontology-based_2014}. The annotation of the 11 surgeries represent a total of 671 activities (e.g. cut the rectum with a monopolar hook), and 16 iAEs. All of these iAEs were bleeding events, some on which lasting only a few seconds, though others over two minutes. Each of them was validated by an additional surgeon who did not perform the surgeries.

At this point, the iSPM is a label sequence composed of a succession of activities defined by three components (action verb, surgical instrument, and anatomic target) \cite{lalys_surgical_2013} and characterized by a duration. Thus, each iSPM is a continuous sequence of activities characterized by a duration (see Module A in Fig~\ref{fig:method}).

\subsection{Creation of a standard surgical process by multi-dimensional non-linear temporal scaling (Module B)}\label{subsec:creationRef}
The objective of this module is to create a standard surgical process that represents the most typical sequence of activities performed by the surgeon, to be used for detecting deviations in the third module. This second module is composed of three steps: 1) sampling the iSPMs; 2) aligning the sampled iSPMs to get the same length for all iSPMs; 3) creating the standard surgical process itself.

\subsubsection{Sampling the individual surgical process models}
As explained previously, iSPMs are continuous sequences. However, to perform the following steps, we need discrete sequences. We thus sampled the iSPMs to achieve this goal. The impact of the sampling rate is analyzed in the validation section.

\subsubsection{Alignment of iSPMs by multi-dimensional non-linear temporal scaling}
\label{MD_NLTS}
To create a standard surgical process according to activity sequences rather than their durations, we needed to temporally align the iSPMs. To this end, we have proposed a new approach called Multi-Dimensional Non-Linear Temporal Scaling (MD-NLTS), inspired by the Non-Linear Temporal Scaling (NLTS) proposed in \cite{forestier_non-linear_2014}.

NLTS is a multiple alignment method developed for one-dimensional surgical processes alignment. It is derived from dynamic time warping (DTW) and involves three steps:
\begin{enumerate}[label=\alph*)]
	\item An average sequences of the set of sequences is computed by DTW Barycenter Averaging (DBA) \cite{petitjean_global_2011};
	\item The average sequence is independently aligned to each sequence of the set, with the aim of defining which elements of the sequence correspond to each element $l$ of the average sequence. Thus, for each element $l$ we have the corresponding set of elements of all sequences and $widths[l]$ and the maximum number of elements in the set sequence corresponding to each element $l$ of the average sequence;
	\item The alignments are finally ``unpacked'': All sequences are warped to include the same number of elements, defined by $widths[l]$, in a way that avoids information loss.
\end{enumerate}

NTLS was created to overcome one DTW limitation. To perform multiple alignments with DTW, one sequence must be chosen as the reference, with the other sequences aligned to this reference. The alignment is thus dependent on the chosen reference. On the contrary, NLTS enables alignment between three or more sequences by computing an average sequence using DBA \cite{petitjean_global_2011}. NLTS realizes a local alignment by focusing on regions with string similarity rather than on all sequences' durations. Moreover, with this alignment, there is no loss of information, given that the sequences are extended during step c of NLTS, so that even an item with few samples will be retained anyway.

Despite these advantages, NLTS only enables the alignment of one-dimensional surgical processes. Thus, when we seek to align activity sequences in which activities are composed of three components (action verb, surgical instrument, and anatomic target), NLTS considers these three components as one dimension. If we have three activities, defined as follows: A\_1: \textless verb\_1, instrument\_1, target\_1\textgreater, A\_2: \textless verb\_1, instrument\_2, target\_1\textgreater{} and A\_3: \textless verb\_3, instrument\_3, target\_3\textgreater, NLTS will consider A\_2 and A\_3 equally different from A\_1, even though the instrument only differs between A\_1 and A\_2.

In \cite{shokoohi-yekta_non-trivial_2015}, two approaches were proposed to achieve multi-dimensional warping: either dependent warping or independent warping. However, these approaches were only applied to classic DTW. Thus, to take into account the benefits of multi-dimensional warping and NLTS, we propose a dependent multi-dimensional warping applicable to NLTS. We chose to develop a dependent warping approach, because the three components are strongly linked within the activity.

We adapted the NLTS cost matrix, used in step a) of NLTS, to develop MD-NLTS. In NLTS, the cost matrix between sequence $Q$ and sequence $C$ is defined as:

\begin{equation}
    d(q_{i},c_{j})=
    \begin{cases}
      0, & \text{if}\ q_{i}=c_{j} \\
      1, & \text{if}\ q_{i} \neq c_{j}
    \end{cases}
\end{equation}

where $q_{i}$ is the label of sequence $Q$ at $t=i$, and $c_{j}$ the label of sequence $C$ at $t=j$.

In MD-NLTS, each sequence is composed of $M$ dimensions. We define element $D$ of the cost matrix as the sum of the distance of each dimension:

\begin{equation}
    D(q_{i},c_{j})=\sum_{m=1}^{M} d(q_{i,m},c_{j,m})
    \label{equ:MD-NLTS}
\end{equation}
with,

\begin{equation}
    d(q_{i,m},c_{j,m})=
    \begin{cases}
      0, & \text{if}\ q_{i,m}=c_{j,m} \\
      1, & \text{if}\ q_{i,m} \neq c_{j,m}
    \end{cases}
\end{equation}

For the three activities, A\_1, A\_2, and A\_3,  previously defined, MD-NLTS will consider A\_1 more similar to A\_2 ($D(_{A\_1,A\_2})=1$) than A\_3 ($D(_{A\_1,A\_3})=3$). This difference will impact step a) by influencing the average sequence created by DBA. The other MD-NLTS steps are similar to NLTS. Following the alignment, all aligned sequences $S_{1..n}^*$ exhibit the same length (see Module B in Figure~\ref{fig:method}).

\subsubsection{Computation of the standard surgical process} \label{creationSSP}
The standard surgical process $\overline{S}$ is created by computing the more frequent activity in all aligned sequences $S_{1..n}^*$, at each instant. Let’s assume that we have the following activities at instant t: $s_1^* [t]= A\_1$, $s_2^* [t]= A\_1$ and $s_3^* [t]= A\_2$. Activity $A\_1$ is more frequent than $A\_2$, so $\overline{S}[t] = A\_1$.

Although MD-NLTS computes an average sequence in its first step, we did not select this as the standard surgical process given that this average sequence does not have the same length as the aligned sequences, rendering it impossible to detect deviations by comparing the activities. The coherence of this proposed standard surgical process was validated with the surgeons (see Section \ref{validationStdSurgPro}).

\subsection{Detection of surgical deviations (Module C)}
The objective of this third module is to detect deviations by comparing an aligned sequence $s$ of $S_{1..n}^*$ to the standard surgical process $\overline{S}$. To compare these sequences, we compute the distance $D(\overline{S}_{t},S^*_{s,t})$ between these two surgical processes at each time-step $t$. Contrary to the computation of $D$ in Eq~(\ref{equ:MD-NLTS}), the distance between the two sequences is computed for the same time-step since the surgical sequences are aligned. Similarly to step \ref{MD_NLTS}, the distance is multidimensional, i.e. the three components of the sequence are taken into account. Deviations are detected at each instant $t$ when $D(\overline{S}_{t},S^*_{s,t})>0$.

\subsection{Classification of deviation types by a hidden semi-Markov model (Module D)}
\label{HsMM}
The objective of a hidden semi-Markov model (HsMM), also called explicit-duration HMM \cite{yu_efficient_2003}, is to explain a non-observable sequence (i.e., the hidden state sequence) using an observable sequence (i.e., the observation sequence) in which it is theoretically conceivable to stay in the same state for an infinite duration. An HsMM is characterized by $\lambda = (\pi, A, B, P)$, where:
\begin{enumerate}
 \item $\pi$ is the initialization matrix containing probabilities to start the sequence at each state;
 \item $A$ is the transition matrix between hidden states containing the probabilities of changing states between two instants. Self-transitions are impossible ($A_{i,i}=0$, for each $i\in[0, nb\_hidden\_state]$);
 \item $B$ is the emission matrix containing the probabilities of producing given observations knowing that we are in a specific state;
 \item $P$ is the state duration matrix defining the probability to stay in a specific state for each possible duration.
\end{enumerate}

In a first step, the HsMM is trained to define model $\lambda$ using observation sequences and true hidden sequences. These sequences, thus, need to be defined for each aligned iSPM. The observation sequences are defined as the concatenation of:
\begin{itemize}
	\item the three components of each aligned sequence's activities;
	\item the distance used to detect deviations between this sequence and the standard surgical process.
\end{itemize}
We define three different hidden states: ``no deviation'' (compared to the standard surgical process), ``context deviation,'' and ``event deviation.'' These latter two are defined according to the definition given in Table~\ref{tab:deviations}.  In our case, we have not ``expert deviation'' due to the fact that all surgeries were performed by a single surgeon. Figure~\ref{fig:TrueStatesRepresentation} presents the creation of the true hidden sequences. If the distance between a specific surgery $s$ and the standard surgical process at instant \textit{t} is not null ($D(\overline{S}_{t},S^*_{s,t})>0$ in Figure~\ref{distToMean}) and an intraoperative adverse event (in Figure~\ref{event}) occurs, the true hidden state is defined as ``event deviation'' (in Figure~\ref{trueStates3states}). If the distance between a specific surgery and the standard surgical process at instant \textit{t} is not null ($D(\overline{S}_{t},S^*_{s,t})>0$), but there is no iAE, the true hidden state is defined as ``context deviation.'' Any other case corresponds to a ``no deviation'' true hidden state.

\begin{figure}[h]
  \subfloat[$D(\overline{S},S_s^*)$ : Distance between a specific surgery and the standard surgical process $\overline{S}$. Blue=0, yellow=1, grey=2, red=3.]{
    \includegraphics[width=\linewidth,height=0.2cm]{DistToMean2.png}
    \label{distToMean}
  }\newline
  \subfloat[Event in a specific surgery. Blue= no event, yellow = intraoperative adverse event.]{
    \includegraphics[width=\linewidth,height=0.2cm]{event2.png}
   \label{event}
  }\newline
  \subfloat[True hidden state. Blue=no deviation, yellow= context deviation, grey= event-deviation.]{
    \includegraphics[width=\linewidth,height=0.2cm]{trueStates3states2.png}
    \label{trueStates3states}
  }
  \caption{{\bf Representation of hidden state creation for sequence $s$.}
	If $D(\overline{S}_{t},S^*_{s,t})=0$ (in a), the true hidden state is ``no deviation'' (in c). If $D(\overline{S}_{t},S^*_{s,t})>0$ (in a), and no event occurs (in b), the true hidden state is ``context deviation'' (in c). If $D(\overline{S}_{t},S^*_{s,t})>0$ (in a) and an intraoperative adverse event occurs (in b), the true hidden state is ``event-deviation'' (in c).}
  \label{fig:TrueStatesRepresentation}
\end{figure}

The training step of our HsMM was performed using the forward-backward algorithm developed by Yu and Kobayashi \cite{yu_efficient_2003} and the enabled creation of detection model $\lambda$.

\section{Results}
\label{Results}
Our method was validated using a leave-one-out cross-validation. The HsMM training was performed on all patients, except one. The remaining operation was sampled and aligned with the standard surgical process $\overline{S}$ in order to create an aligned test surgery $S^*_{test}$. We have computed the distance between $\overline{S}$ and $S^*_{test}$. With this distance and $S^*_{test}$, we have computed the true hidden state sequence and observation sequence, as explained in Section \ref{HsMM}. Deviations were detected by feeding the observation sequence to the trained detection model $\lambda$. Each model was evaluated by comparing the detected deviation sequence with the true hidden state sequence.

We investigated the impact of the sampling rate (Section \ref{subsec:creationRef}) on the results. To this end, we varied the sampling rates between 2 to 12 samples-per-second in 1-second steps. We also studied the results at 12.5 samples-per-second, given that this sampling rate corresponds to half of the video frequency (25Hz).

The distribution of observations between each hidden state is very heterogeneous  (Table~\ref{tab:heterogeneity}): 68\% of them belonging to the ``no deviation'' state, 26\% to the ``context deviation'' state, and only 6\% to the ``event deviation'' state. Moreover, this distribution is also very heterogeneous  between surgeries especially for ``event deviation'' state, with a standard deviation of 6.80\% and a large range (minimum of 0\% and a maximum superior of 18\%). Due to these heterogeneities, we could not be satisfied with accuracy only, as performance metrics. Given the small amount of ``event deviation'' occurrences, we can reach an accuracy of 94\% if all observations belonging to the ``no deviation'' and ``context deviation'' states are correctly classified, even if none of the ``event deviation'' states were detected. We, thus, used recall and precision to accurately estimate our model's ability to classify event deviation types. All results are given with a confidence interval (CI) of 95\%.

\begin{table*}[!ht]

	\centering
	\caption{{\bf Distribution of observations for each hidden states.} }

    \begin{tabular}{|c|c|c|c|c|c|}\hline
    	&Mean (\%)&STD(\%)&Median(\%)&Min(\%)&Max(\%)\\\hline
        No deviation&	68.41&	4.22&	69.52&	60.05&	74.91\\\hline
		Context deviation &	25.86&	5.73&	26.91&	12.19&	32.74\\\hline
		Event deviation&	5.73&	6.80&	2.83&	0.00&	18.13\\\hline
    \end{tabular}
  \label{tab:heterogeneity}
\end{table*}

Kendall's Tau, a non-parametric test, was performed to examine a possible statistical correlation between the sampling rate and accuracy, precision and recall for each hidden state. We chose a level of 0.05 to consider the correlation statistically significant. With seven statistical tests (one for accuracy, three for recall, and three for precision), we were in the context of multiple comparisons. To counteract the problem of false-positive results in multiple comparisons, we employed the Bonferroni correction method \cite{bonferroni_teoria_1936,dunn_multiple_1961}. Therefore, the statistical significance level was set at 0.0071 (0.05/7).

We additionally investigated the model errors more closely in order to understand the reasons underlying failed classification of deviation types, with the aim to classify event deviations (Section \ref{error_identification}).

\subsection{Validation of the standard surgical process}\label{validationStdSurgPro}
The deviation detection method proposed is based on the learning phase of our HsMM. It is thus dependent on the hidden state sequences defined, due to the distance between a specific surgery and the standard surgical process (Figure \ref{distToMean}). It is essential to validate the consistency of the standard surgical process (section \ref{creationSSP}) in terms of surgical workflow, i.e. the sequence of computed activities that could happen in a real sequence. If the standard surgical process is not consistent, such as closure of the body occurring before the first skin incision, the deviation detection will not be correct.
To carry out this validation, we computed multiple standard surgical processes. One was computed with all iSPMs available, and the others by removing one or more iSPMs before computation. We have randomly shown to two surgeons graphical representations of real surgeries (iSPMs) and graphical representations of standard surgical processes. Several examples are available as supplementary material. The surgeons were asked if they deemed that, according to their surgical expertise, each representation was consistent in terms of surgical workflow, and whether each representation corresponded to either real surgery or computed standard surgical process. They considered all representations consistent, being unable to distinguish standard surgical processes from real surgeries.
With this validation, based on experts’ opinions, we assessed the computation of standard surgical processes produces realistic sequences in terms of surgical workflow.

\subsection{Deviation classification results}
\label{observations_Results}
Figure~\ref{fig:NDvsSDvsED} and Table~\ref{tab:NDvsSDvsED} present the results of the deviation classification for the three metrics (accuracy, recall, and precision) at different sampling rates. A table providing all results is available as supplementary material.

\begin{table*}[!ht]

	\centering
	\caption{{\bf Results of the classification of surgical deviations.} }

    \begin{tabular}{|c|c|c|c|c||c|c|}\hline
    	\multicolumn{2}{|c|}{Samples/sec}& 2 & 8 & 12 & tau & p-value\\\hline
    	\multicolumn{2}{|c|}{Accuracy(\%)}& 75.69 - 83.65 & 84.63 - 93.43 & 82.31 - 91.23 & 0.2727 &0.1248\\\hline
    	\multirow{3}{*}{\rotatebox{90}{Recall}}& ND(\%)&94.67 - 99.29 & 96.40 - 99.34& 95.67 - 99.21&0.24.24 &0.1554\\\cline{2-7}
    	&CD(\%)&42.39 - 64.77& 53.92 - 89.66& 46.87 - 79.33&0.0606 &0.4203\\\cline{2-7}
    	&ED(\%)&11.32 - 29.72& 33.09 - 82.53& 42.38 - 88.32&0.3636 &0.0580\\\hline
    	\multirow{3}{*}{\rotatebox{90}{Precision}}&ND(\%)&80.41 - 88.43&96.19 - 99.19& 98.83 - 99.56&0.6667 & 0.0009*\\\cline{2-7}
    	&CD(\%)&61.04 - 86.68 & 78.21 - 94.33& 71.97 - 100&0.4848 &0.0155\\\cline{2-7}
    	&ED(\%)&10.38 - 46.46& 9.73 - 57.29& 13.33 - 37.05&-0.0909&0.6808\\\hline

    \end{tabular}

    \begin{flushleft}  ND: ``no deviation.'' CD: ``context deviation.'' ED: ``event deviation.'' The star (*) represents a significant relationship between the sampling rate and results.
		\end{flushleft}

  \label{tab:NDvsSDvsED}
\end{table*}

\begin{figure}[h]

\centering
 \subfloat[Accuracy]{%
    \includegraphics[height=3.5cm]{NDvsSDvsED-accuracy.png}%
    \label{NDvsSDvsED-accuracy}

 }
 \subfloat[Recall]{
    \includegraphics[height=3.5cm]{NDvsSDvsED-PC.png}
    \label{NDvsSDvsED-PC}
 }

 \subfloat[Precision]{
    \includegraphics[height=3.5cm]{NDvsSDvsED-PV.png}
    \label{NDvsSDvsED-PV}
 }
 \caption{{\bf Graph of the classification of surgical deviations.}
 Dotted lines represent the trend curve of the 95\% confidence intervals. In Figures~(b) and (c), the blue solid line corresponds to the ``no deviation'' state, the red dotted-solid line to the ``context deviation'' state, and the green dashed line to the ``event deviation'' state.
 }
 \label{fig:NDvsSDvsED}

\end{figure}
\FloatBarrier

The accuracy ranged between 80 \% and 90 \% with a one-half the width of the CI inferior to 6 \% for two samples-per-second or more. We could observe a non-significant (p-value=0.1248) upward trend in accuracy along with the sampling rate (Figure~\ref{NDvsSDvsED-accuracy}).

\modif{The recall for the ``no deviation'' state ranged between 96 \% and 99 \% with a CI range inferior to 5 \% for all sampling rates}. On the other hand, the recall for the ``context deviation'' state fluctuated between a mean of 50 \% and 72 \% without any specific trend, however the CI range (red dotted line) trend to disperse with the sampling rate. For the ``event deviation'' state, the sampling rate appeared to impact the recall (green dashed line in Figure~\ref{NDvsSDvsED-PC}). For two samples-per-second, the recall as a 95\% CI of [11.32 - 29.72], whereas it was [42.38 - 88.32] for 12 samples-per-second. However, this trend proved to be not significant (p-value=0.0580) . Results are very different between sequences as we could shown with the CI.

The precision for the ``no deviation'' state significantly (p-value = 0.0009) increased from [80.41 - 88.43] to [98.83 - 99.21]  with the sampling rate and a half CI range inferior to 6 \%. The precision for the ``context deviation'' state fluctuated between a mean of 73 \% and 98 \% with a half CI range inferior to 15 \%, except for 12.5 samples-per-second where the half CI range is 21.16 \%, the trend being statistically not significant (p-value=0.0155). \modif{For the ``event deviation'' state, the mean precision was less than 40 \% for all sampling rates}. As for the recall, the CI demonstrate a high variability between sequences for the precision of ``event deviation'' state, with for example result between [9.73 - 57.29] for 8 samples per second.

\subsection{Analysis of the model's errors in classifying event deviations}
\label{error_identification}
To understand why the mean precision for ``event deviation'' state was less than to 40 \%, we analyzed the classification errors. We noticed that over 98 \% of the observations falsely classified as ``event deviation'' were actually representative of  ``context deviation.'' These observations were then classified into four categories:

\begin{enumerate}
 \item Rarely wrongly classified: less than 1\% of observations belonging to this observation type (same action verb, surgical instrument, anatomic target, and distance $D$) were wrongly classified by the model.
 \item Untrained: this observation type was not present in the training dataset.
 \item Correctly trained: in the training dataset, this observation type was characteristic of the ``event deviation'' state.
 \item Other: observations that did not belong to the previous three categories.
\end{enumerate}

The distribution of the observations within of these four categories has been provided in Figure~\ref{fig:NDvsSDvsED-FalseED}. Category 1 observations (rarely wrongly classified) were negligible for all sampling rates, accounting for less than 1 \% of errors. Category 2 (untrained) and 3 (correctly trained) observations accounted for less than 20 \% of falsely-classified observations, with a downward trend in the sampling rate for category 2, and upward trend for category 3 (Figure~\ref{fig:NDvsSDvsED-FalseED}). Category 4 (other) represented more than 50 \% of the falsely-classified observations, though we observed a downward trend with increasing sampling rates.

\begin{figure}[h]
\centering
     \includegraphics[height=3.5cm]{NDvsSDvsED-FalseED.png}
     \caption{{\bf Graphical distribution of the four categories of falsely classified observations for different sampling rates.} The green line corresponds to the ``rarely wrongly classified'' category, blue line to the ``untrained'' category, yellow line to the ``correctly trained'' category, and red line to the ``other'' category.}
     \label{fig:NDvsSDvsED-FalseED}
\end{figure}

\section{Discussion}
Our validation study's results have clearly shown the high accuracy of our approach in detecting deviations, for all sampling rates, during the surgical process. However, the distinction between the deviation types does not prove to be effective, and results have high variability between sequences. Our model classifies a number of "context deviation" states as "event deviation" states. However, from a patient's safety point of view, it proves crucial that deviations due to intraoperative events are not missed, even if false detections do occur.

We have deeply analyzed the misclassification between the two deviation types (Figure~\ref{fig:NDvsSDvsED-FalseED}). Based on this analysis, four observation categories were extracted, with each category interpreted as follows:
\begin{itemize}
    \item Category 1 errors (rarely wrongly detected observations) are caused by the time taken by  the model  to perform a state transition;
    \item Category 2 errors (untrained observations) are due to a lack of data: because the model did not encounter the observation in the training phase, it provides an arbitrary result. A larger dataset might reduce the number of observations pertaining to this category;
	\item Category 3 errors (correctly trained observations) could be caused by activities that would have led to an ``event deviation,'' as observed in the training dataset, but were corrected by the surgeon before they occurred. A study of surgical behavior is warranted to confirm this hypothesis;
	\item a large number of observations were classified as Category 4 (other). Our interpretation for this error type is that surgical activities, such as defined today (action verb, surgical instrument, and anatomic target), may not capture enough information concerning the surgical scene. To solve this issue, a more refined description of activities appears necessary. For instance, by combining the laparoscopic images with registered complementary per-operative imaging modalities, such as ultrasounds \cite{lanchon_augmented_2016,billings_system_2012} or fluorescence imaging \cite{voros_devices_2013}, it may be easier to identify vasculature information and express the distance between the surgical instrument and this underlying vasculature as a label (too close, close or safe). This information could be treated similarly to any activity component in our approach.
\end{itemize}

The CI range demonstrates a high variability of the results between sequences, especially for ``event deviation'' results. This could be explained by the high variability on the distribution of each type of deviation, as shown in Table~\ref{tab:heterogeneity}, and the limited size of our database. Indeed, when we test the sequence with the more important distribution of event deviations, we have a higher probability to encounter observations never encounter during the training phase.

Our study presents several limitations. First, the entire study has been based on manual annotations performed by one observer only. Neumuth et al. \cite{neumuth_validation_2009} studied the reliability of the annotation process and concluded that ``granularity was reconstructed correctly by 90\%, content by 91\%, and the mean temporal accuracy was 1.8 s." This temporal variability is reinforced by Huaulm\'e et al. \cite{huaulme_automatic_2019} which demonstrated a temporal inter-variability which could have a relative standard deviation superior to 18\% for activities when multiple observers are implied. According to the literature, we could thus consider having similar results, i.e. if we implied several observers, the identified activities will be very similar but we will introduce temporal variability. As explained in section \ref{MD_NLTS}, this works was based on activity sequences rather than their durations. Consequently, this variability would not be relevant information for our method.

Second, although we have validated the consistency of the standard surgical process $\overline{S}$ (Section \ref{creationSSP}), in terms of surgical workflow, its creation can be a source of false deviation classification. Remember, we have determined the activity at each instant \textit{t} by selecting the most frequent activity in the aligned sequences. However, if the aligned sequences are too heterogeneous at one instant, this most frequent activity might not be present in the majority of aligned sequences (e.g. it might only represent a small percentage of the activities at time t, even though it is the most frequent activity). For further developments of this approach, it will be paramount to consider the probability of the chosen activity or allow for alternative surgical paths within the standard surgical process.

The third limitation concerns the choice of only considering the dominant hand of the surgeon for classifying deviations. It could be of interest to add the information provided by the non-dominant hand to further investigate its influence on deviation detection. To this end, the activities of the second hand must be annotated, while our method has to be improved by modifying our Multi-Dimensional Non-Linear Temporal Scaling method and using a coupled Hidden semi-Markov Model with two observation sequences (one per hand), and one hidden state sequence.

Finally, our dataset comprises of surgeries performed by a single surgeon, while we have not considered differing habits or expertise levels among the surgeons. We have, therefore, removed one level of complexity. Furthermore, the dataset includes bleeding events only. We cannot predict the performance of our method for other types of iAEs. However, our annotation methodology, which only relies on activity annotations, would be identical whatever the type of iAE, as long as the iAE start and end times can be identified by the expert surgeons involved in the annotation process. Moreover, we did not take into account at this stage organizational/context factors into our annotation process. Indeed, operating room disruptions due to e.g. the composition of the surgical team could also result in deviations. They are currently considered as ``operating room disruptions" and belong to context deviations. Future works would examine the robustness of our approach by including more types of iAEs, multi-surgeon data to include the expert deviation classification and study the influence of the surgical team composition on deviation occurrences. Collect data from multiple surgeons and other surgical team members will allow using more complex approaches than a simple ``leave-one-out" one, e.g. a ``leave-one-user-out" or one where the couple surgeon/assistant is excluded from the training.

\section{Conclusion}
Surgical deviation classification is challenging and should enable us to understand the hidden processes underlying their occurrence. We have, herein, proposed the first offline method for automatically classifying deviations based on their type (event deviation or context deviation). The method, namely Multi-Dimensional Non-Linear Temporal Scaling followed by a Hidden semi-Markov Model, has provided interesting initial results, whereas its precision still needs to be improved.

The detection of event deviations is an important preliminary step towards the identification of iAEs. Indeed, event deviations are a marker of the occurrence of one or multiples iAEs. This could help determine the exact moment when iAEs occur. Moreover, the objective of an event deviation is to ``correct or limit the impacts of iAEs'' (Table~\ref{tab:deviations}), so by studying the anatomical structure concerned by event deviations, it will be possible to determine which one is impacted. Of course, to make a complete identification of iAEs, further work will be necessary.

To propose routine surgery applications of our method, two further improvements are required. The first is to develop an on-line multi-dimensional alignment method. Recently, Forestier \textit{et al.} \cite{forestier_optimal_2015} proposed a method designed to create an online one-dimensional alignment. The second aspect pertains to creating a reliable and automatic online activity recognition method \cite{katic_ontology-based_2013,despinoy_unsupervised_2015,dergachyova_automatic_2016}. With these two developments available, a real-time implementation of our method will be rendered possible.

\section*{Conflict of interest statement}
The authors declare that they have no conflict of interest.

\section*{Acknowledgements}
\noindent
This work was partially supported by French state funds managed by the ANR within the Investissements d'Avenir programme (Labex CAMI) under reference ANR-11-LABX-0004.

\noindent
Authors thanks the IRT b\textless \textgreater com for the provision of the software ``Surgery Workflow Toolbox [annotate]'' , used for this work.

\bibliographystyle{elsarticle-num}


\section{Conclusions}
\label{sec:conclusion}
In this paper, we present our four-class sleep stage classification trained by two BLSTM layers. We use our low-cost wearable multi-sensor system, SensEcho, to acquire subjects' RR intervals and breathing signals. Three novel features are designed during the feature extraction to detect the sudden variation on RR intervals. Our robust prediction accuracy is 80.25\% on a large public dataset (417 subjects), and 80.75\% on our 32 enrolled subjects, respectively. These results are better than the previous works which either used small data sets with the potential over-fitting issues, or used the conventional machine learning methods on large data sets.
\section{Discussions, Potentials And Future Work}
\label{sec:discussion}

As we have shown that our approach compared to the latest contributions has an outstanding sleep stage classification accuracy on a large dataset, there are many potential directions in our future works to further pursue a better performance. Authors in \cite{zhao2017learning} leverage CNN to create effective features which are difficult interpreted and extracted manually. We are currently improving the BLSTM model and designing an appropriate CNN architecture to integrate with our current feature sets. Meanwhile, alongside recruiting healthy subjects, we are enrolling actual patients using our wearable SensEcho who have heart, breathing, or sleep issues/diseases, in order to build our comprehensive sleep datasets for the intra-cohort and inter-cohort sleep researches. As aforementioned, a 3-axis accelerometer is embedded in SensEcho which can infer the body movement behaviors in previous works \cite{muzet2016assessing,hong2018noncontact}. Our future works will also involve the use of accelerometer data.

\section{Evaluations}
\label{sec:eval}

\begin{figure*}[t]
	\centering
	\subfloat[Experiment set-up]{%
		\includegraphics[width=0.18\linewidth]{zhangyuan.pdf}%
		\label{fig:zhangyuan}
	}
	\hfill
	\subfloat[Four-class confusion matrix on the validation set of the SHHS dataset]{%
		\includegraphics[width=0.185\linewidth]{confmat-SHHS-crop.pdf}%
		\label{fig:confmat-SHHS}
	}
	\hfill
	\subfloat[Four-class confusion matrix on the test set (our own dataset)]{%
		\includegraphics[width=0.185\linewidth]{confmat-own-crop.pdf}%
		\label{fig:confmat-ours}
	}
	\hfill
	\subfloat[CDF of the prediction accuracies in different datasets]{%
		\includegraphics[width=0.18\linewidth]{cdf-crop.pdf}%
		\label{fig:cdf}
	}
	\hfill
	\subfloat[Visualization of sleep epochs using t-SNE]{%
		\includegraphics[width=0.185\linewidth]{tsne-crop.pdf}%
		\label{fig:tsne}
	}
		\vspace{-5pt}
	\caption{\small{Experiment set-up and visualization of results}}
	\label{fig:results}
		\vspace{-15pt}
\end{figure*}

\subsection{Datasets}
\subsubsection{Dataset from the public database}
To mitigate the training over-fitting issue using a small dataset, we leverage the large public sleep PSG database \emph{Sleep Heart Health Study} (SHHS) \cite{quan1997sleep} to build our training model. SHHS consists of the PSG monitoring during sleep of more than 6,000 subjects (52.4\% females) in the U.S. across the age range of $\sim$40+ years, including the records of six-class sleep stage classification for each subject manually determined by clinical specialists (wake, REM, S1, S2, S3, S4). It is worth noting that in the SHHS database numerous subjects have various sleep-related diseases, such as breathing disorders or irregularities, insomnia, etc. Those subjects would introduce serious bias to the training model. To minimize these impacts, we filter the data based on the following steps.
\begin{enumerate}
	\item We exploit the metric Apnea Hypopnea Index (AHI) score to indicate subjects' apnea levels, and further help to select the appropriate cohort for the training model. AHI score is calculated by the summation of manually identified breathing irregularity events. It could be commonly divided into four levels \cite{nandakumar2015contactless}: \emph{No apnea}: AHI $<$ 5; \emph{Mild}: 5 $\leq$ AHI $<$ 15; \emph{Medium}: 15 $<$ AHI $\leq$ 30; and \emph{Severe}: AHI $>$ 30. The subjects who have no apnea (AHI $<$ 5) are chosen as candidates to build our data set from SHHS.

	\item Among the subject candidates, we then select the subjects who have at least 5\% of S3 and S4 stages, and at least 15\% REM stage, during their sleep periods, since sleep with such quality is regarded as \emph{regular sleep} \cite{fonseca2017cardiorespiratory}

	\item We combine the S1 and S2 stages as the light sleep stage, and combine the S3 and S4 stages as the deep sleep stage \cite{berry2012aasm}. With the existing stages of wake and REM, the six-class sleep stages are converted into four-class ones.

	\item As a result, 417 subjects are selected to construct the dataset from the SHHS database. We randomly take 70\% of them for the training set, while the remaining 30\% contribute to the validation set. The subjects constituting the training set never appear in the validation set.
\end{enumerate}

\begin{figure}[t]
	\begin{center}

		\includegraphics[width=\linewidth]{best-average-worst-crop.pdf}
				\vspace{-13pt}
		\caption{\small{Examples of best, average, worst classifications in the validation set of SHHS (referring to the validation set in Fig.~\ref{fig:cdf})}}
		\label{fig:three-cases}
				\vspace{-20pt}
	\end{center}
\end{figure}

\subsubsection{Our own dataset}
Aside from using the public database, we are also building our own dataset for sleep stage classification. Till now, 32 subjects (22 -- 45 ages and 15 females) with no apnea (AHI $<$ 5) are enrolled to wear SensEcho during the PSG tests at Chinese PLA General Hospital (shown in Fig.~\ref{fig:zhangyuan}).
Their four-class sleep stages are identified by three clinical specialists' consent. Since the size of the dataset is now still limited, it is served as the test set only. In other words, we leverage the trained model from the SHHS training set described above to make the prediction on our own dataset. Note that this study is reviewed and approved by the Ethics Committee of Chinese PLA General Hospital (IRB number: S2018-095-01)

\subsubsection{Consistency of the datasets}
As aforementioned in Sec.~\ref{sec:feature-extraction}, we align our data with the sleep stage determination based on the standard PSG records using the 30-second epoch. Such processing approach guarantees that the model trained from SHHS dataset can be applied to predict the test set (i.e., our own dataset).

\subsection{Sleep Stage Classification}

Table \ref{tab:result} lists the recent literature which contributes to the four-class sleep stage classification using RR interval/heart rate and breathing effort, where some works also involve the information of body movements. As we see, the sizes of datasets used in their works are very diverse. Compared to the previous works using large dataset \cite{tataraidze2016sleep,aggarwal2018sleep}, the BLSTM method in our work significantly outperform them. The reason is that BLSTM is able to learn the hidden relationship in time-series correlated signals. Meanwhile, aside from overcoming the over-fitting problem for the model training, we keep the similar accuracy level compared to the works which use small data sets to build the learning model \cite{zhao2017learning,hong2018noncontact}.
Note that we let \Large{$\kappa$} \normalsize denote \emph{Cohen's kappa coefficient}, which is a more robust metric than the accuracy percentage as it considers the occurrence by chance \cite{cohen1960coefficient}.

Fig.~\ref{fig:confmat-SHHS} presents that the confusion matrix of four sleep stage classes on the validation set of the SHHS dataset. We observe that except deep sleep, all other three stages achieve great classification accuracies. However, plenty of deep sleep epochs are classified into light sleep epochs. This observation aligns with the examples of these two stages in Figs.~\ref{fig:light} and \ref{fig:deep}, where the patterns of RR intervals and breathing signals are difficult to distinguish. We then apply the prediction model to our test set to obtain another confusion matrix result (shown in Fig.~\ref{fig:confmat-ours}). While the high prediction accuracies of wake and REM stages are maintained, it is interesting that the accuracy of deep sleep increases much but the accuracy of light sleep declines. It is worth recalling the different age ranges of the training and validation sets and the test set. As observed in \cite{peters2014age}, the light sleep of young subjects is less than it of senior subjects, while the deep sleep of young subjects is more. Furthermore, the ethnic differences (White in SHHS vs. Asian in our dataset) and gender distributions also lead to dissimilar sleep patterns \cite{rao2009ethnic}. Due to these reasons, we believe that there is some over-fitting issue during the light sleep stage training but that the younger people's deep sleep patterns are clearer to be identified and predicted by the prediction model.
Fig.~\ref{fig:cdf} shows that the Cumulative Distribution Function (CDF) of the accuracies of the three sets in terms of subjects. Even though there is some generational bias in the training model, the overall prediction accuracy on the test set still maintains the same level as it on the validation set. The smaller size of test set might result in the slightly better prediction result.

\begin{figure}[t]
	\centering
	\includegraphics[width=\linewidth]{features-crop.pdf}%
	\caption{\small{Top 20 significant features ranked in XGBoost. (\textit{RMSSD}: root mean square of the successive differences}}
	\label{fig:feature-importance}
	\vspace{-10pt}
\end{figure}

We then select the best case, the worst case, and one average case from the validation set to view the prediction details. Fig.~\ref{fig:three-cases} shows that the major impact on the prediction accuracy is the transition between the light sleep stage and the deep sleep stage, which is not well recognized, especially in the case that such transitions occur many times during the sleep period (i.e. the worst case). However, the transition between the light sleep stage and the wake stage is well predicted across these three cases. Considering the distinct difference in RR intervals and breathing signals of these two stages in Figs.~\ref{fig:wake} and \ref{fig:light}, our three novel features are designed for capturing such a sudden transition even though it lasts for only one epoch. We will evaluate them in the next section.

\subsection{Knowledge of Features}

In this section, we investigate how the extracted features contribute to the learning process. Fig.~\ref{fig:tsne} visualizes the sleep epochs of one subject using the t-SNE \cite{maaten2008visualizing}. We notice that the system is able to learning much coherent information across the four classes. The interweaving between the deep sleep epochs and the light sleep epochs are more than all other two stage pairs. Since it is still difficult to interpret the significances of specific features, we leverage the recent common used tree model, XGBoost \cite{chen2016xgboost}, to sort all the features' relative importance according to their \emph{weights}. The weight refers to the times a feature is used to split the data across all the trees. Fig.~\ref{fig:feature-importance} lists the top 20 important features which contribute to the classification most significantly. It is worth noting that our three designed features $f_1$, $f_2$, and $f_3$ are ranked at No. 4, 8, and 19, respectively. We also observe that several features extracted from CPC are top significant, indicating the indispensable synergistic effects of cardiorespiratory signals on different sleep stages.
Note that the prediction accuracy using XGBoost is 75.6\% on the validation set.

\section*{Acknowledgments}
This work was supported by NSF of China (61471398), Beijing Municipal Science and Technology (Z181100001918023), CERNET Innovation Project (NGII20160701), Translational Medicine Project of Chinese PLA General Hospital (2016TM-041) and Big Data Research \& Development Project of Chinese PLA General Hospital (2018MBD-09).

\bibliographystyle{IEEEtran}

\section{Introduction}
\label{sec:intro}

Sleep quality is essential to human being's health. To analyze it, sleep is usually classified into multiple stages in which physiological signals have different patterns, related to various physiological functions.
Although there are many criteria to classify the sleep stages, the most common four-class sleep stages criterion is our focus, which includes \emph{Wake}, \emph{Light sleep}, \emph{Deep sleep}, and \emph{Rapid Eye Movement (REM)} \cite{berry2012aasm}. For example, deep sleep helps body metabolism and wound healing, while dreaming occurs during the REM period most. The accurate classification of sleep stages is necessary to properly analyze one's sleep architecture, further for the sleep disorder diagnose and recovery.

Polysomnography (PSG) is a standard method for sleep quality analysis, which has been used for sleep stage classification in past decades, as it provides the electrical activities of various body parts. Electroencephalography (EEG) with the assistance of Electrooculography (EOG) and Electromyography (EMG) during the PSG test is used to determine the different sleep stages.
However, the entire procedure of sleep stage classification using PSG is expensive and labor-intensive, as it requires the dedicated equipment and laboratory as well as the clinical specialists' expertise. Rather than relying on EEG signals from PSG, recent researchers focus on using cardiorespiratory signals because they can also indicate the different sleep stages and can be captured in plenty of the low-cost simplified devices. This approach is easily adapted into the various scenarios of mobile and Internet of Things (IoT) healthcare, sleep health monitoring in smart home \cite{yang2017vital} and other accommodations \cite{dor2012experiences}.
As the deep learning methods are rapidly developing in recent years, some latest researches have demonstrated the effectiveness that applying Recurrent Neural Network (RNN) to capture the implicit patterns of the time-series cardiorespiratory signals  \cite{zhang2017sleep,zhao2017learning}.

In this paper, we present our ongoing work that uses our designed wearable multi-sensor system to capture subjects' ECG and breathing signals, and that applies a Bi-directional RNN architecture with Long Short-term Memory (BLSTM) \cite{graves2005framewise} to classify the four sleep stages. An overview of our contributions is provided below:

\begin{enumerate}
	\item We propose and demonstrate three novel features which are effective to detect the sudden variation of RR intervals\footnote{The term of \textit{RR interval} will be explained in Sec~\ref{sec:rr-intervals}.}. A total of 152 features are extracted for the learning procedure.
	\item Compared to the recent works, we leverage a large public dataset and our own dataset to achieve the outstanding sleep stage prediction in terms of both accuracy (80.25\% on the public dataset, and 80.75\% on our own dataset) and sample size (total 449 subjects). These results outperform the previous works on the four-class sleep stage classification, which either used small data sets having the potential over-fitting issues \cite{fonseca2015sleep,muzet2016assessing,zhao2017learning,hong2018noncontact}, or used the conventional machine learning methods on large data sets \cite{tataraidze2016sleep,aggarwal2018sleep}.
	\item We design our own low-power and low-cost wearable multi-sensor vest, \emph{SensEcho}, with the support of the secure cloud system. We validate that it is robust to acquire the cardiorespiratory signals compared to the standard PSG test.
\end{enumerate}

This paper is organized as follows. Related works are discussed in Section~\ref{sec:related-works}. We describe our hardware and methods in Section~\ref{sec:methods}. The evaluations are presented in Section~\ref{sec:eval}. Section~\ref{sec:discussion} proposes further potentials and future works. We finally conclude our work in Section~\ref{sec:conclusion}.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.6\linewidth]{main-board-smaller.pdf}
	\caption{\small{Hardware of SensEcho}}
	\label{fig:main-board-smaller}
	\vspace{-15pt}
\end{figure}
\section{Hardware and Methods}
\label{sec:methods}

\begin{figure*}[t]
	\centering
	\hspace{.1pt}
	\subfloat[Wake]{%
		\includegraphics[width=0.24\linewidth]{wake-crop.pdf}%
		\label{fig:wake}
	}
	\hfill
	\subfloat[Light sleep]{%
		\includegraphics[width=0.24\linewidth]{light-crop.pdf}%
		\label{fig:light}
	}
	\hfill
	\subfloat[Deep sleep]{%
		\includegraphics[width=0.24\linewidth]{deep-crop.pdf}%
		\label{fig:deep}
	}
	\hfill
	\subfloat[REM]{%
		\includegraphics[width=0.24\linewidth]{rem-crop.pdf}%
		\label{fig:rem}
	}
		\vspace{-5pt}
	\caption{\small{Examples of RR intervals and breathing signals in different sleep stages}}
	\label{fig:types}
	\vspace{-18pt}
\end{figure*}

\subsection{Hardware}
Our wearable multi-sensor system, SensEcho \cite{zhang2018automated,zhang2018breathing,lan2019tachy}, exploits three comfortable electrode patches to achieve the single-lead ECG signal monitoring at a 200 Hz sampling rate. Two sensing wires are embedded at the locations of chest and abdomen for the two corresponding breathing behaviors monitoring at a sampling rate of 25 Hz. In addition, an ultra-low power, 3-axis digital accelerometer component ADXL345 at a sampling rate of 25 Hz is also embedded in the vest. The main board of SensEcho is shown in Fig.~\ref{fig:main-board-smaller}. Fig.~\ref{fig:signal-overview} presents an example fragment to show the signal quality of our system compared to the PSG records. A wrist oximeter wirelessly communicates with SensEcho via Bluetooth, whose sampling rate is 1 Hz. SensEcho provides the local and secure cloud data storage options
The 2GB-capacity local storage option enables us to backtrack the raw data when cloud storage is neither stable nor available. The embedded battery power can support continuous monitoring for five consecutive days.

\subsection{Data Pre-processing}
\label{sec:rr-intervals}
\textbf{ECG signals:} Rather than tackling the raw ECG signals, we extract \emph{RR intervals} from ECG signals.
RR interval is the time duration between two adjacent R-peaks of ECG signals
It can be used to calculate heart rate as well as Heart Rate Variability (HRV).

\textbf{Breathing signals:} As breathing signals always have the offset from the zero line of amplitude, we apply wavelet decomposition technologies to remove the offset effects. A Butterworth low-pass filter with the frequency of 1 Hz is next applied to remove signal noises at high frequency.

Fig.~\ref{fig:types} presents that the RR intervals and breathing signals on a subject along with the four-class sleep stage determination, respectively. It is obvious that the wake and the REM stages have distinguishing patterns on both RR intervals and breathing signals, while the light and deep sleep stages look similar to each other but minute differences on the RR intervals can be still found. Note that a specialist can easily distinguish these two sleep stages relying on the information from EEG, EOG, and EMG.

\subsection{Feature Extraction}
\label{sec:feature-extraction}

In this section, we elaborate our process of feature extraction. As a sleep stage is determined based on every 30-second \emph{epoch} in PSG tests \cite{quan1997sleep}, we define our sliding processing window unit as the same size. This way enables our time-series data to align the sleep stage labels in PSG records. Our feature extraction is processed on either one sliding window unit (one epoch), or a larger window consisting of several consecutive epochs. The moving step size is defined as one epoch. Note that when multiple consecutive epochs (a larger sliding window) are needed, the amount of epochs is an odd number such that the calculated feature would be exactly associated with the middle epoch.
For example, extracting a feature from 3 consecutive epochs means that a feature is calculated inside a 90-second window and assigned to the second epoch.

\begin{table*}[t]
	\caption{\small{Comparison of four-class sleep stage classification (\textit{BE}: breathing effort; \textit{BLD}: Bayesian linear discriminant; \textit{BM}: body movement; \textit{BT}: boosted trees; \textit{CNN}: convolutional Neural Network; \textit{CRF}: conditional random field; \textit{GAN}: generative adversarial network; \textit{GBM}: gradient boosting machine; \textit{HB}: heart beat; \textit{RandF}: random forest; \textit{RR}: RR interval; \textit{SKNN}: subspace KNN; \textit{SVM}: support vector machine)}}
	\vspace{-5pt}
	\begin{tabular}{cccccccc}
		\hline
		\textbf{Authors} & \textbf{Data type} & \textbf{Dataset source} & \textbf{Data acquisition} & \textbf{Number of subjects} & \textbf{Method} & \textbf{Accuracy} & \Large{$\kappa$} \\
		\hline
		\cite{fonseca2015sleep} & RR+BE & SIESTA\cite{klosh2001siesta}+Own & PSG & 15+27 & BLD & 69.0\% & 0.49 \\
		\cite{tataraidze2016sleep} & RR+BE & SHHS & PSG & 625 & GBM+RandF & 71.4\% & 0.57 \\
		\cite{aggarwal2018sleep} & BE & MESA\cite{bild2002multi} & Ventilator & 400 & Deep CRF & 74.1\% & 0.57 \\
		\cite{muzet2016assessing} & RR+BM & Own & PSG & 48 & SVM & 74.9\% & n/a \\
		\cite{zhao2017learning} & HB+BE & Own & RF signals & 25 & CNN+RNN+GAN & 79.8\% & 0.70 \\
		\cite{hong2018noncontact} & HB+BE+BM & Own & RF signals & 13 & BT, SKNN, etc. & 81.0\% & n/a \\
		\multirow{4}{*}{\textbf{Ours}} & \multirow{4}{*}{\textbf{RR+BE}} & \multirow{4}{*}{\textbf{SHHS+Own}} & & \multirow{4}{*}{\textbf{417+32}} & \multirow{4}{*}{\textbf{BLSTM}} & \multicolumn{2}{c}{\textbf{Validation set from SHHS}}\\
		& & & \textbf{PSG (SHHS)} & & & \textbf{80.25\%} & \textbf{0.71} \\
		\cline{7-8}
		& & & \textbf{Wearable (own)} & & & \multicolumn{2}{c}{\textbf{Test set from own}} \\
		& & & & & & \textbf{80.75\%} & \textbf{0.69} \\
		\hline
	\end{tabular}
	\label{tab:result}
	\vspace{-15pt}
\end{table*}

\subsubsection{Features from RR intervals}
We extract numerous features from the time domain and the frequency domain on RR intervals.

\textbf{Time domain features:} We extract 10 features commonly used for the HRV analysis \cite{redmond2007sleep}, and 34 conventional statistical features on RR intervals \cite{fonseca2015sleep}, such as mean, quantiles, range, etc. We also extract 5 non-linear features including sample entropy, zero crossing analysis \cite{xiao2013sleep}. Nevertheless, the sudden changes of RR intervals are not well captured by using these features only. To address this issue, we design three novel features as follows:
\begin{equation}
f_1 = \overline{R^{mid}_n} - \overline{R_n}
\label{eq:f1}
\end{equation}
\vspace{-10pt}
\begin{equation}
f_2 = \overline{R^{mid}_n} - \widetilde{R_n}
\label{eq:f2}
\end{equation}
	\vspace{-10pt}
\begin{equation}
f_3 = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(\overline{R^i_n} - \overline{R_n})^2}
\label{eq:f3}
\end{equation}
where $R_n$ represents raw RR intervals in the consecutive $n$ epochs; $R^i_n$ denotes the $i$th epoch of $R_n$; $R^{mid}_n$ refers to the middle epoch of $R_n$; $\overline{R_n}$ represents the average value of $R_n$; $\widetilde{R_n}$ denotes the median value of $R_n$. Such three features investigate the impact of sudden variation of RR intervals in one epoch over the longtime series. In practice, we set $n=119$ in $f_1$ to approximate one-hour-long time series, and set $n=9$ in $f_2$ and $f_3$ to consider 9 consecutive epochs (4.5 minutes).

\textbf{Frequency domain features:} In order to obtain the frequency domain information, we first perform the linear interpolation on RR intervals, and then use Fast Fourier Transformation (FFT) on the given epoch(s). According to \cite{xiao2013sleep}, 21 features of the frequency domain are extracted, such as mean, spectrum power, entropy, etc.

\subsubsection{Features from breathing signals}
Similar to the feature extraction in RR intervals, we extract 25 statistical features from breathing signals according to \cite{redmond2007sleep,fonseca2015sleep}.
For example, in the time domain, mean and standard deviance of breathing peak sequence, kurtosis, skewness, etc. are extracted; features in the frequency domain include dominant peak, energy, etc.

\subsubsection{Features from cardiopulmonary coupling}
Cardiopulmonary coupling (CPC) is a technology to analyze the connection between HRV and breathing volume variability. CPC index \cite{thomas2005electrocardiogram} is a set of magnitudes in the frequency domain to present CPC degrees at a given time spot. We calculate CPC indexes in $R_n$ and obtain the summation of CPC indexes at the high-frequency (HF) band (0.1--0.4Hz), at the low-frequency (LF) band (0.01--0.1Hz), and at the very-low-frequency (VLF) band (0--0.01Hz), respectively. The ratios of these summations at the specific frequency bands over the CPC index summation across the entire frequency domain are then calculated.

We finally extract 152 features from RR intervals, breathing signals, and cardiopulmonary coupling effects. At last, the Z-score normalization is performed on each feature such that its mean value is 0 and its variance value is 1. The feature vectors are fed into our learning model architecture. We then leverage two BLSTM layers, each of which has 16 units, and one 4-unit output layer corresponding to the four sleep stage classes. The advantage of BLSTM over conventional LSTM is that the BLSTM architecture is able to learn the target's both previous and future information.

\section{Related Works}
\label{sec:related-works}

\subsection{EEG for Sleep Stage Classification}
Since specialists rely on EEG, EOG, EMG signals to determine sleep stages, automated classifying sleep stages based on those signals are well studied to reduce the input of labor and cost.
Authors in \cite{pan2012transition} presented a transition-constrained discrete hidden Markov models for automatic sleep stage classification using the EEG, EOG and EMG signals of 20 subjects.
A deep learning model, DeepSleepNet, was proposed in \cite{supratak2017deepsleepnet} for automated grading sleep stages using the single-channel EEG. DeepSleepNet consists of Convolutional Neural Networks (CNN) to extract time-invariant features and BLSTM to learn stage-to-stage transition rules from EEG signals. Again, automated classification relying on EEG signals still requires PSG tests.

\subsection{Cardiorespiratory Signals for Sleep Stage Classification}
\begin{figure}[t]
	\includegraphics[width=\linewidth]{PSG-sensecho.pdf}
	\caption{\small{Overview of signals from PSG and SensEcho}}
	\vspace{-5pt}
	\label{fig:signal-overview}
\end{figure}

Since it is difficult to compare the works when the number of sleep stages is not consistent, we mainly review the previous works only on the four sleep stage classification as our focus.

Large datasets \cite{tataraidze2016sleep,aggarwal2018sleep} containing the cardiorespiratory recordings were leveraged to predict the sleep stages. Due to the limitation of conventional machine learning algorithms, the prediction accuracies were moderate. Authors in \cite{fonseca2015sleep} extracted a total of 142 features from electrocardiogram and thoracic respiratory effort of 25 subjects, and applied Bayesian linear discriminant for the four-class sleep stage classification on a 30-second epoch basis. Similarly, the traditional Support Vector Machine algorithm was applied to the cardiorespiratory features of 48 subjects \cite{muzet2016assessing}. Recently, the contactless sensing methods such as using Radio Frequency (RF) signals to acquire the heart beats and breathing efforts have gained much attention. Authors in \cite{zhao2017learning} proposed a new learning framework to involve various deep learning methods on the RF signals of 25 subjects to classify sleep stages into four classes with the high prediction accuracy. The BLSTM-based five-class sleep stage classification was proposed in \cite{zhang2017sleep}, using heart rate and wrist actigraphy acquired by a wristband. They achieved the higher accuracies when the sizes of datasets were much smaller. Compared to the works mentioned above, our results accomplish the significant improvement on the prediction accuracy, and the robust performance on a large dataset. We will have a detailed comparison in Table~\ref{tab:result} in Section \ref{sec:eval}.

\section{System Overview}
\label{sec:system-overview}

\todo{system flow diagram; }


\section{Introduction}

Strong attractive interactions between the $\eta$ meson and nucleons mean that there is a chance to form $\eta$ meson bound states in nuclei \cite{Haider:1986sa}.
If discovered in experiments, these mesic nuclei would be a new state of matter bound just by the strong interaction without electromagnetic Coulomb effects playing a role. Strong interaction bound states are formed in a different way as compared to exotic atoms which involve binding of electrically charged mesons with nuclei.
For the latter, negatively charged pions or kaons could replace an electron in an outer orbital in a standard atom and get bound in the atom due to the Coulomb interaction.~The charged meson in such an excited state quickly undergoes transitions to the lower states until it is close enough to the nucleus and is either absorbed by the nucleus or lost in a nuclear reaction.
For strong interactions,
in contrast to the pion, the neutral $\eta$ meson is special due to the strong attractive nature of this meson-nucleon interaction~\cite{Haider:1986sa}.~An off-shell $\eta$ meson produced in nuclear reactions such as the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} 2\gamma$ and $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} 6\gamma$ below the $\eta$ production threshold may form a bound state with the nucleus within which it is produced. Thus the absence of the electromagnetic interaction and the attractive nature of the $\eta$-nucleon interaction, makes the case of the neutral $\eta$ meson different from that of the pion or the kaon and opens the possibility for an exotic nucleus made up of the meson and nucleons.
Early experiments with low statistics
using photon~\cite{Pheron:2012aj,Baskov:2012yd},
pion~\cite{Chrien:1988gn}, proton~\cite{Budzanowski:2008fr}
or deuteron~\cite{Afanasiev:2011zza,Moskal:2010ee,Adlarson2013,Adlarson:2016dme}
beams gave hints for possible $\eta$ mesic bound states but no clear signal \cite{Metag:2017yuh,Kelkar:2013lwa}.

Here we present a new high statistics search for \mbox{$^3$He-$\eta$} bound states with data from the WASA-at-COSY experiment. We focus on the two main neutral decay channels of the $\eta$ meson:
$\eta \to 2\gamma$ with branching ratio
39.41$\pm$0.20\% and $\eta \to 3\pi^0 \to 6\gamma$
with branching ratio
$31.54 \pm 0.22$\%~\cite{PDG}.
These processes constitute more than 70\% of the $\eta$ decays.
The choice of neutral decay channels
minimizes final state interactions
involving charged particles.
Concurrent measurement of the two channels increases the statistics and enables one to control systematic uncertainties in photons detection.
The two-photon decay was previously suggested
in \cite{Bass:2005hn}
as a clean probe of the $\eta$ in nuclear media.

Considering the $\eta$-nucleus interaction, bound states can be formed by the attractive interaction with finite level width corresponding to the finite lifetime of the state due to the absorptive interaction with the nucleus. The momentum distribution of the bound $\eta$ meson determines the sum of the momenta of the emitted photons.
Nuclear absorption and the additional $\eta$ decay (disappearance) processes, reduces significantly the in-medium branching ratio of 2$\gamma$ and 6$\gamma$ decay channels~\cite{Skurzok_Hirenzaki2019}.

$\eta$ meson interactions with nucleons and nuclei are a topic
of great experimental and
theoretical interest.
For recent reviews see
\cite{Metag:2017yuh,Kelkar:2013lwa,Bass:2018xmz,Krusche:2014ava,Wilkin:2016mfn}.
Possible $\eta$-nucleus binding energies are related
to the $\eta$-nucleon optical potential and to the
value of $\eta$-nucleon scattering length $a_{\eta N}$
\cite{Ericson:1988gk}.
Phenomenological
estimates for the real part
of $a_{\eta N}$ are typically
between 0.2 and 1 fm
depending on the model assumptions.
$\eta$ bound states in helium require a large $\eta$-nucleon
scattering length
with real part greater than about 0.7--1.1~fm
\cite{Barnea:2017epo,Barnea:2017oyk,Fix:2017ani}.
Recent calculations in the framework of
optical potential~\cite{Xie:2016zhs},
multi-body calculations~\cite{Barnea:2017oyk},
and pionless effective field theory~\cite{Barnea:2017epo}
suggest a possible $^3$He-$\eta$ bound state.

Modifications of meson properties
are expected in medium.~In studies of the transparency of nuclei to propagating mesons produced in photoproduction experiments one finds strong $\eta$ absorption in nuclei~\cite{Nanova:2012vw}.
For the $\eta'$ one finds weaker interaction with the nucleus. An effective mass shift for the $\eta'$ in medium has been
observed by the CBELSA/TAPS Collaboration~\cite{Nanova:2013fxl}.
The $\eta'$-nucleus optical potential
$V_{\rm opt} = V_{\rm real} + iW$
deduced from these photoproduction experiments with a
carbon target is
$V_{\rm real} (\rho_0)
= m^* - m
= -37 \pm 10 \pm 10 \ {\rm MeV}
$
and
$
W(\rho_0) = -10 \pm 2.5 \ {\rm MeV}$
at nuclear matter density $\rho_0$.
This mass shift is very close to the prediction of the
Quark Meson Coupling mode (QMC)
with mixing angle -20 degrees \cite{Bass:2005hn,Bass:2013nya},
which also predicts a potential depth about -100 MeV
for the $\eta$ at $\rho_0$.
The $\eta'$ results are also
consistent with scattering length estimates from COSY-11~\cite{Czerwinski:2014yot} and Bonn~\cite{Anisovich:2018yoo}. Experimental search for $\eta'$ - nucleus bound states has also been performed with results reported in Ref. \cite{Tanaka}.

Hints for possible $\eta$ helium bound states are inferred from the observation of strong interaction in the $\eta$ helium system. One finds a sharp rise in the cross section
at threshold for $\eta$ production in both
photoproduction from $^3$He
\cite{Pheron:2012aj,Pfeiffer:2003zd}
and in the proton-deuteron reaction
$dp\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} \eta$
\cite{Adlarson:2018rgs}.
These observations may hint at a reduced $\eta$ effective mass in the nuclear medium.

Previous bound state searches at COSY have been focused on the reaction
$dd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} N \pi$~\cite{Adlarson2013,Adlarson:2016dme}. Studies of the excitation function around the threshold for $dd\rightarrow$ $^{4}\hspace{-0.03cm}\mbox{He} \eta$
did not reveal a structure
that could be interpreted as a narrow mesic nucleus.
Upper limits for the total cross sections for bound state production and decay in the processes
$dd\rightarrow$ ($^{4}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} n \pi^{0}$
and
\mbox{$dd\rightarrow$ ($^{4}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} p \pi^{-}$}
were deduced to be about 5~nb and 10~nb for
the $n \pi^0$ and $p \pi^-$ channels respectively \cite{Adlarson:2016dme}.
The bound state production
cross sections
for $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound}$ \cite{Wilkin:2014mla}
are expected to be more than 20 times larger than for
$dd\rightarrow$ ($^{4}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound}$~\cite{Wycech:2014wua}.

In May 2014 the experiment searching for $\eta$ mesic $^3$He nuclei was performed at the COSY accelerator~\cite{COSY_description1,COSY_description5} in J\"ulich, Germany.
The measurements were carried out using the WASA-at-COSY detector~\cite{WASA_description1,WASA_description2,WASA_description3,WASA_description4,3Heeta_cross_section}.
The mesic nuclei are supposed to be formed in proton-deuteron collisions.
A ramped proton beam with beam momentum varying in the range from $1.426$ to $1.635$~GeV/c corresponding to $^3$He$\eta$ excess energy range from $-70$ to $30$~MeV and a pellet deuterium target~\cite{WASA_pellett_target} were used.
The $^3$He-$\eta$ bound state was searched for in the $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $ ^{3}\hspace{-0.03cm}\mbox{He}2\gamma$  and $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ decay channels. These channels that manifest the direct decay of $\eta$ bound in $^3$He nucleus have been investigated for the first time. The existence of the bound $^3$He-$\eta$ state would manifest itself as a maximum or interference pattern in the excitation function for both of the studied reactions below the $pd\rightarrow ^3$He$\eta$ reaction threshold.

For the normalization of the excitation functions, the integrated luminosity was determined as a function of the excess energy.~The analysis is presented in the next section. Further on, the data selection and efficiency determination is described. The data analysis is followed by the interpretation of the achieved excitation functions in view of the possible signal from the $\eta$-mesic $^3$He.

\section{Luminosity determination}
Luminosity was determined based on the $pd\rightarrow ^3$He$\eta$ and $pd\rightarrow ppn_{spectator}$ reactions.
The $pd\rightarrow ^3$He$\eta$ reaction analysis allows one to estimate the integrated luminosity for $^3$He$\eta$ excess energy $Q_{^3\hspace{-0.05cm}He\eta}$ above zero. The $^3$He particles were registered in the forward detector~\cite{WASA_description1} and identified using the $\Delta E-E$ method based on energy losses in scintillator layers (see~Fig.~\ref{fig_pid}).

\begin{figure}[h!]
\centering
	\includegraphics[width=13.0cm,height=6.0cm]{3He_ident.png}
	\caption{2-D histograms of energies deposited in the first layer of Forward Trigger Hodoscope (FTH1) and the first layer of Forward Range Hodoscope (FRH1) for all events with signal in Forward Proportional Chamber (FPC) (left panel) and events that were identified as $^{3}\hspace{-0.03cm}\mbox{He}$ (right panel).}\label{fig_pid}
\end{figure}

The count of events originating from this reaction was obtained based on the $^3$He missing mass spectra for each excess energy interval separately. An example spectrum is shown in Fig.~\ref{he3_mm_spectra_fit}.
The reconstruction efficiency was calculated using Monte Carlo simulations taking into account the experimental data on cross sections and angular distributions~\cite{3Heeta_cross_section,3Heeta_cross_section_older1,3Heeta_cross_section_older2,3Heeta_cross_section_older3}.

\begin{figure}[h!]
	\includegraphics[width=10.0cm,height=6.0cm]{He3eta-fit.png}
	\caption{\small
		$^3$He missing mass spectrum obtained from data for the excess energy range of $Q_{^3\hspace{-0.05cm}He\eta}~\in~[20.0;22.5]$~MeV. The part of the spectrum that is considered to be background is shown with green color and is fitted with a polynomial of fourth power (orange).
	}
	\label{he3_mm_spectra_fit}
\end{figure}

The $pd\rightarrow ppn_{spectator}$ reaction analysis allows one to determine the integrated luminosity for the whole beam momentum range.~As far as the target overlapping by the beam is changing during the acceleration cycle, the integrated luminosity value can change depending on the beam momentum.~The registration efficiency for the $pd\rightarrow ppn_{spectator}$ reaction was obtained with dedicated Monte Carlo simulations described  in Refs.~\cite{OleksandrPhD-arXiv,czyzyk}.
The distribution of relative proton-neutron motion inside the target deuteron was calculated based on the parametrisation of the Paris potential~\cite{PARIS_model}. Data on the proton-proton elastic scattering cross section and
the angular distribution
~\cite{PP_CS_link} were used for simulating the quasi-elastic scattering in the framework of the spectator model.~The calculated cross section was multiplied by the factor $0.96$
to take into account the shading effect~\cite{shading_effect}. It is worth noting that above the $\eta$ production threshold, the two estimates of luminosity are in agreement (based on the $pd\rightarrow ppn_{spectator}$ and $pd\rightarrow ^3$He$\eta$ reactions~\cite{OleksandrPhD-arXiv}). The total integrated luminosity was determined to be $2446\pm3(stat.)\pm66(syst.)\pm4$(norm.)~nb$^{-1}$
where the statistical, systematic and normalisation errors are indicated, respectively~\cite{OleksandrPhD-arXiv}. This is the largest statistics ever obtained for these experimental conditions.

\section{The analysis of $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ and $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reactions}

As a first step, in order to establish the optimal selection criteria, Monte Carlo simulations for the $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ and $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reactions were performed in the framework
of the spectator model with the assumption of an isotropic distribution of bound $\eta$ meson decay products
in its rest frame. The momentum of the $\eta$ meson was simulated using the recent
model~\cite{Skurzok_Hirenzaki2019} in which the $^3$He-$\eta$ relative momentum distribution was calculated by solving the Klein-Gordon equation assuming the potential of $\eta$-nucleus interaction based on Hiyama's density distribution in $^3$He~\cite{Hiyama1,Hiyama2,Hiyama3}.

For the $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ reaction analysis, the events containing a $^3$He track in the forward detector and at least two photons in the central detector were selected.
If there were more than two photons, the pair with the invariant mass closest to the $\eta$ mass corrected by $Q_{^3\hspace{-0.05cm}He\eta}$ value was chosen.
Then the restrictions on $^3$He missing mass, $\gamma$-$\gamma$ missing mass, and $\gamma$-$\gamma$ invariant mass were applied using selection ranges based on the simulated distributions~\cite{OleksandrPhD-arXiv}.
The excitation function obtained for the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ reaction is shown in the left panel of Fig.~\ref{two_gamma_events}.
\begin{figure}[h!]
	\includegraphics[width=6.0cm,height=6.0cm]{He3gg-events-final.png}
	\includegraphics[width=6.0cm,height=6.0cm]{He36g-events-final.png}
\caption{The dependence of determined events count on $Q_{^3\hspace{-0.05cm}He\eta}$ for $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ reaction (left panel) and $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reaction (right panel). The error bars include both statistical and systematic uncertainties.}
	\label{two_gamma_events}
\end{figure}

The signal from the bound state is expected for excess energies around or below zero. The increase of events above 10~MeV is due to the $pd\rightarrow^3$He$\eta$ reaction. It starts at 10~MeV because of a hole for the COSY beam in the geometrical acceptance of the WASA-at-COSY detector (see Fig.~\ref{effic}).

\begin{figure}[h!]
	\includegraphics[width=6.0cm,height=6.0cm]{eff_2gamma.png}
	\includegraphics[width=6.0cm,height=6.0cm]{eff_6gamma.png}
	\caption{The efficiency for different reactions when applying selection criteria defined for the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} 2\gamma$ (left) and $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He} 6\gamma$ (right) reaction analysis.}\label{effic}
\end{figure}

For the $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reaction analysis, the events containing a $^3$He track in the forward detector and at least six photons in the central detector were selected.
For each combination forming three pairs, to identify the $\eta\rightarrow3\pi^0\rightarrow6\gamma$ decay, the following quantity is calculated:
\begin{equation}
	D= \sum_{i=1}^{3} (m_{\gamma_{(2i-1)}\gamma_{2i}} - m_{\pi^0})^2
	\label{6gamma_3pi0_identification}
\end{equation}
where $m_{\gamma_{(2i-1)}\gamma_{2i}}$ is the $\gamma$ pair invariant mass and
$m_{\pi^0}$ is $\pi^0$ mass. The combination of six photons
that minimises $D$ was choosen.
Then analogous to the $2\gamma$ case, the selection conditions on the $^3$He missing mass, $6\gamma$ invariant mass, and $6\gamma$ missing mass were applied based on the simulated distributions~\cite{OleksandrPhD-arXiv}. The excitation function obtained for the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reaction is shown in the right panel of Fig.~\ref{two_gamma_events}.

The excitation curves have been normalised using the integrated luminosity values calculated based on the $pd\rightarrow ppn_{spectator}$ reaction and the efficiency determined based on Monte Carlo simulations.
The results for both studied reactions are shown in Fig.~\ref{excitation_fit_linear}.

\begin{figure}[h!]
		\includegraphics[width=10.0cm,height=6.0cm]{UpperLimit-LinearFit1.png}
		\includegraphics[width=10.0cm,height=6.0cm]{UpperLimit-LinearFit2.png}
	\caption{
		Excitation curves determined for the
		$pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ (upper panel) and $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ (lower panel) reactions.
		Superimposed lines indicate result of the fit of the line.
		The points above the $\eta$ production threshold are excluded from the analysis.
	}
	\label{excitation_fit_linear}
\end{figure}

\section{The upper limit for the $\eta$ mesic \mbox{$^3$He} production cross section}

The excitation curves obtained in the analysis (Fig.~\ref{excitation_fit_linear}) did not reveal any resonance-like structures and the fit with linear functions results in $\chi^2$ value $<$ 1 when normalized to the number of degrees of freedom.
This indicates that no strong signal from the bound $^3$He-$\eta$ state is observed.

Further on, for the quantitative estimates of the upper limits for the bound state production, a fit to the excitation curves with a linear function (for background) plus  a Breit-Wigner function (for the signal) was performed. The fit was done for different combinations of the assumed $\eta$-mesic $^3$He binding energies $B_s$ and widths $\Gamma$. The value of $\Gamma$ was tested in the range from 1.25~MeV to 38.75~MeV (with the step of 2.5~MeV) and  $B_s$ in the range from 1.25~MeV to 63.75~MeV (with the step of 2.5~MeV).
\begin{figure}[h!]
		\includegraphics[width=10.0cm,height=6.0cm]{UpperLimit--137-287Fit1.png}
		\includegraphics[width=10.0cm,height=6.0cm]{UpperLimit--137-287Fit2.png}
	\caption{
		Exemplary result of the simultaneous fit of functions~\ref{excitation_curves_fit_function1} and \ref{excitation_curves_fit_function2} to the experimental data for the assumed $B_s$ and $\Gamma$ values as indicated above the figures. Superimposed
		black line shows the full fit result, and the green line shows the background function only.
	}
	\label{excitation_fit_BW}
\end{figure}

For a given $B_s$ and $\Gamma$ pair, the following functions were fit simultaneously for the two studied reaction channels:
\begin{equation}
	\rho^{fit}_{^3\hspace{-0.05cm}\mathrm{He}2\gamma}(Q_{^3\hspace{-0.05cm}\mathrm{He}\eta}) = P_{\eta\rightarrow2\gamma} \cdot \sigma \cdot \sigma_b(Q_{^3\hspace{-0.05cm}\mathrm{He}\eta}) + p_1 Q_{^3\hspace{-0.05cm}\mathrm{He}\eta} + p_2,
	\label{excitation_curves_fit_function1}
\end{equation}
\begin{equation}
\rho^{fit}_{^3\hspace{-0.05cm}\mathrm{He}6\gamma}(Q_{^3\hspace{-0.05cm}\mathrm{He}\eta}) = P_{\eta\rightarrow6\gamma} \cdot \sigma \cdot \sigma_b(Q_{^3\hspace{-0.05cm}\mathrm{He}\eta}) + p_3 Q_{^3\hspace{-0.05cm}\mathrm{He}\eta} + p_4.
\label{excitation_curves_fit_function2}
\end{equation}
Here $\sigma$, $p_1$, $p_2$, $p_3$, and $p_4$ are the free fit parameters, $P_{\eta\rightarrow2\gamma}$ and $P_{\eta\rightarrow6\gamma}$ are the branching ratios for the $\eta\rightarrow 2\gamma$
and $\eta\rightarrow 6\gamma$ decays.
Assuming that the ratio of branching ratios for the $\eta\rightarrow2\gamma$ and $\eta\rightarrow3\pi^0$ decay channels for the bound $\eta$ meson remain the same as in vacuum, the vacuum branching ratio values of
	$P_{\eta\rightarrow2\gamma}$~=~0.3941 and
	$P_{\eta\rightarrow3\pi^0\to 6\gamma}$~=~0.3268
were used for performing the fit~\cite{PDG}. The function $\sigma_b(Q_{^3\hspace{-0.05cm}He\eta})$ in the fit formulae represents a Breit-Wigner shape which for a given values of $B_s$ and $\Gamma$ reads:
\begin{equation}
	\sigma_b(Q_{^3He\eta},B_s,\Gamma)
	= \sigma~\frac{\Gamma^2/4}{(Q_{^3\hspace{-0.05cm}He\eta}-B_s)^2+\Gamma^2/4}.
	\label{Breight_Wigner_formulae}
\end{equation}
Example results of the fit are shown in Fig.~\ref{excitation_fit_BW}. The figure shows results for the $B_s$ and $\Gamma$ values (indicated above the plots) for which the fitted values of $\sigma$ differ from zero with the largest statistical significance.
Fig.~\ref{upper_limit_1d} indicates the results of the fit as a function of the $B_s$ for the most promising value of $\Gamma~=~28.75$~MeV.

The upper limit of the total cross section was determined based on the fit parameter uncertainty $\Delta\sigma^{stat}$:

\begin{equation}
\sigma^{CL=90\%}_{upper}(B_{s},\Gamma)=\sigma+k\Delta\sigma^{stat},
\end{equation}

\noindent where $k$ is the statistical factor equal to 1.64 corresponding to 90\% confidence level as given in PDG~\cite{PDG}).
Fig.~\ref{upper_limit_1d} shows the systematic
limits (blue lines) in addition to the statistical uncertainties (green lines).
Systematic errors were estimated by changing the parameters of all cuts applied in the data analysis, and changing the values of assumed potential parameters for the $^3$He-$\eta$ interaction that determines the Fermi momentum distribution for relative motion in the bound state. The highest contribution to the systematic error is connected with the background fit function. The uncertainty due to the fit of quadratic or linear function estimated as $\sigma_{quad}-\sigma_{lin}$ varies from about $2$ to $5~nb$.

\begin{figure}[h!]
		\includegraphics[width=10.0cm,height=6.0cm]{UpperLimit-Gconst287-Parameter.png}
	\caption{
		Upper limits for the bound state production cross section via $pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}$($\eta~ \mbox{decays})$
		as function of binding energy for fixed width $\Gamma$=28.75~MeV.
		The values of the Breit-Wigner amplitude $\sigma$ are shown with statistical uncertainties. The range of possible bound state production cross section obtained based on statistical uncertainty corresponding to $90\%$ confidence level is shown by blue lines.
		The range of possible bound state production cross section including  systematic uncertainty is shown by green lines.
	}
	\label{upper_limit_1d}
\end{figure}

In the obtained excitation functions one can see a slight signal from the possible bound state for $\Gamma>20$~MeV and $B_{s}\in[0;15]$~MeV corresponding to the optical potential parameters $-100<V_{0}<-70$ MeV and $|W_{0}|>20$ MeV in the model described in \cite{Skurzok_Hirenzaki2019}.
The result is also consistent
with the QMC prediction of a
potential depth about -100 MeV at nuclear matter density \cite{Bass:2005hn} and with the models in Refs.~\cite{Barnea:2017epo,Barnea:2017oyk,Xie:2016zhs,kelkar:prl2007}. The allowed $V_{0}$-$W_{0}$ area is however different to those deduced from the $\eta$-$^4$He
system~\cite{Skurzok_PLB2018}
using the optical model of \mbox{Ikeno et al. \cite{Ikeno:2017xyb}} where most of the model parameter space was excluded allowing values of the real and imaginary parts of the potential only between zero and about -60~MeV and -7~MeV respectively.
However, the observed signal is within the range of the
systematic uncertainty. Hence one cannot make definite conclusions whether $\eta$-mesic $^3$He exists with the decay mechanism studied here.

\section{Conclusions}
The analysis of the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ and $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reactions has been performed in order to search for the existence of an $\eta$-mesic $^3$He state.
The analysis of the obtained excitation functions for the $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}2\gamma$ and $pd\rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}6\gamma$ reactions shows slight indication of
the signal from the bound state for $\Gamma>20$~MeV and $B_{s}\in[0;15]$~MeV.
However, the observed indication is within the range of the systematic error which does not allow one to make a definite conclusion on a possible bound state formation.

The upper limit for the cross section of the bound state production varies between $2$ and $15$~nb depending on the bound state parameters.
It is however important to stress that the determined upper limit concerns the production of the $(^3\mbox{He}$-$\eta)_{bound}$ state and its subsequent disintegration via decay of the $\eta$ meson. The branching ratio for the latter in the nuclear medium remains to be estimated theoretically.

This is the first result obtained for the direct decay of bound $\eta$ meson.
The upper limit is much lower than the limit of $70~nb$ for \mbox{$pd\rightarrow$ ($^{3}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}\pi^{0}$}
reaction obtained by the \mbox{COSY-11} Collaboration~\cite{light_cosy11_prev}
and is comparable with the upper limits obtained for the
$dd\rightarrow$ ($^{4}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}n\pi^{0}$ and $dd\rightarrow$ ($^{4}\hspace{-0.03cm}\mbox{He}$-$\eta)_{bound} \rightarrow$ $^{3}\hspace{-0.03cm}\mbox{He}p\pi^{-}$reactions~\cite{Adlarson:2016dme}.
The much improved constraint will help tuning theoretical modelling of the $\eta$-nucleon and $\eta$-nucleus
interactions.

\section*{Acknowledgements}
We acknowledge the support from the Polish National Science Center through grant No. 2016/23/B/ST2/00784, and from the Foundation for Polish Science through the MPD and TEAM POIR.04.04.00-00-4204/17 programmes. Theoretical parts of this work was partly supported by the Faculty of Science, Universidad de los Andes, Colombia, through project number P18.160322.001-17,
and by JSPS KAKENHI Grant Numbers JP16K05355 (S.H.) in Japan.

\section*{References}