In the following lines, we will briefly present some specific ECG methods proposed in the literature over the past years, which contain results similar to the methods we presented in this paper, except for the fact of using patient-specific dictionaries or involving updating the dictionary when there are changes in the ECG signal. In general, there is a big inconvenience in the situation of using such a system in practice, because it involves

**Citation:** Fira, M.; Costin, H.-N.; Goras, , L. A Study on Dictionary Selection in Compressive Sensing for ECG Signals Compression and Classification. *Biosensors* **2022**, *12*, 146. https://doi.org/10.3390/ bios12030146

Received: 28 December 2021 Accepted: 24 February 2022 Published: 27 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

resubmitting the dictionary and necessary calculations in real time to see if the dictionary is good or needs to be updated. All these calculations imply additional hardware needs, which can make the method less practical in real-time acquisition situations. On the other hand, our approach is based on the use of non-modified patient-specific dictionaries or pathology-specific dictionaries; these are established once and updating can be done less frequently than in other techniques and does not require real-time decisions.

In one paper [33], the presented method uses an over complete wavelet dictionary, a dictionary that is later reduced due to a training phase. In addition, it is proposed to align the beats according to the position of the R-peak. This alignment aims to exploit the different scaling characteristics of ECG waves in the wavelet dictionary optimization process. Three different methods are tested for dictionary optimization. It should be mentioned that this optimized dictionary is specific to the patient and for its construction, the first 5 minutes of registration are taken. For acquisition, the authors use a matrix optimized for the ECG signal to be acquired through CS. The use of an optimized compression matrix leads to improved results, but has the disadvantage that once this matrix is changed it must be sent together with the compressed ECG signal. That means both the compressed ECG signal and the compression matrix must be sent to restore the ECG signal.

Another approach is presented in [34], where the quasi-periodic character of the ECG signal is used to detect similarities between ECG pulses and to transmit segments that show dissimilarities normally, without compression. This approach is proposed because abnormal frames, which could be signs of heart disease, are not similar to normal frames. Thus, only the ECG segments considered normal are transmitted by CS, the rest being transmitted normally. Once it is determined, whether the heartbeat is acquired normally or by CS, a quantization step follows and then a Huffman compression. These two steps lead to improved compression results. A critical point in the method is the correct detection of normal vs. abnormal beats, because this automated detection is debatable in the light of the fact that normality or abnormality is determined by a cardiologist and the accuracy of the acquisition should not be influenced by this decision.

In paper [35], the authors also used CS associated to dictionaries built specifically for the ECG signal, thus using the dictionary learning technique to construct a better sparsifying basis to improve the compression ratio. Moreover, the authors consider the change of ECG signal characteristics and propose a physiological variation detection technique and a low-complexity dictionary refreshing algorithm to update the dictionary from time to time when the current dictionary is no longer suitable for the patient.

Many papers in the CS field focus on optimizing the measurement matrix, i.e., the matrix is used in the acquisition stage or on optimizing the necessary calculations in this stage by arranging this matrix in a way that allows easy hardware implementation of the necessary calculations. In practical implementations, the simple random or Bernoulli matrix may have the inconvenience of the required number of operations. Thus, in paper [36], the authors propose an optimized algorithm for collecting the compressed ECG signal, based on the proposed optimization of a deterministic binary block diagonal matrix. The blocks, which make up the diagonal of the matrix, are identical and contain m = N/M elements each, where M and N, respectively, represent the number of rows and columns.

In paper [37], a new method of compressive sampling of ECG signals is presented, which is based on the idea of building the compression matrix adapted to the frame of the ECG signal to be compressed. Thus, a circulating matrix is proposed, containing zeros and ones, obtained by quantizing (with 1-bit resolution) the size of the ECG signal. The detection matrix adapted in this way guarantees that the significant portions of the waveform of the compressed ECG signal are in fact contained in the compressed version. In this way, a more precise reconstruction is guaranteed in relation to the methods already available in the literature. For the reconstruction stage, the acquisition matrix is then used in combination with a modified wavelet dictionary, which also allows the reconstruction of the signal deviation for each processed frame. The big disadvantage of the method is that

whenever the acquisition matrix has to be updated, it has to be sent to the receivers and for reconstruction we have to know each frame with which matrix was collected.

In this paper, we propose a detailed comparative study of two different approaches regarding the possibility of compressed sensing specific ECG signals. This study considers several acquisition techniques/projection matrices used in the acquisition stage and several dictionaries used in the ECG signal reconstruction stage. We will also analyze the effect of preprocessing on the results.

Broadly speaking, we analyze and discuss two CS approaches dedicated to ECG signals, namely:

For both approaches from above, we will analyze several projection matrices, namely, matrices with random independent and identically distributed (i.i.d.) elements taken from the Gaussian or Bernoulli distribution and project matrices optimized for the particular dictionary used in the reconstruction. To optimize the projection matrix, the method presented in [7] will be used.

Furthermore, we will pay special attention to the way the dictionary is built. We will also present the advantages and disadvantages of each and the choice of the method that depends on the available hardware and software resources.

The paper is organized as follows: Section 2 is dedicated to the types of sampling vectors, projection matrices and dictionary construction methods. Section 3 presents the CS methods dedicated to ECG signals. Section 4 shows the results obtained. In Section 5 the results from the previous section are compared and in Section 6 conclusions are drawn.

#### **2. Compressed Sensed Overview**

Traditionally, signals are acquired according to the sampling theorem [8] that states that an f0-bandlimited signal can be recovered from its samples if the sampling frequency is at least 2 f0, i.e., twice the highest frequency of the signal spectrum. Thus, in a time window W, an f0-bandlimited analog signal can be represented by N = 2 f0 W samples equally spaced at T = 1/2 f0, i.e., as a vector belonging to the space RN. Such a signal can be alternatively defined by using any complete set of orthogonal functions in RN. In fact, sampling is nothing else than taking projections (scalar products) on the elements of the canonical basis. In the general case the signal can be reconstructed from its projections on N orthogonal (or only linear independent) elements in RN the canonical basis being the most frequently used. However, in practice, there are cases in which a signal can be reconstructed from fewer samples or projections on an appropriate set of signals, compared to the number prescribed by the sampling theorem. This is possible since the samples contain unnecessary information, and thus, these signals can be compressed and recovered using projections and previous known information. An example would be the class of sparse or rare signals [9–11] that allow a representation based on a small number of elements/atoms in RN. In signal processing literature, the name "k-sparse" denotes signals that can be reconstructed by means of k of elements of RN, the most significant situation being that in which *k<<N*. A discrete signal or vector *x* ∈ *RN* is k-sparse if there exists a base Ψ = {Ψ*i*, *i* = 1, . . . , *N*} in RN so that most of the elements *α* = {*αi*, *i* = 1, . . . , *N*} of its representation in that basis, *x* = Ψ*α*, are zero. Alternatively, they can be approximately zero, so that the signal can be

represented accurate enough with the k's largest terms *αi* from its expansion with respect to that basis. The CS concept is based on theory that a k-sparse signal, i.e., a signal that can be compressed into a base (or, more general, dictionary) Ψ can be recuperated with very good quality from a number *m* of the order of scale *m* = *O*(*k* log(*N*/*k*)) of non-adaptive linear projections on a set of vectors Φ, which are not comprehensible with the first, i.e., their elements cannot be used for a compressed representation of any Ψ*i*, *i* = 1, ... , *N*. Therefore, for obtaining the measurement signal instead of measuring the *N* components of the signal in the canonic base, a number of m (*k* < *m* << *N*) linear projections on the elements of the matrix Φ*N*∗*m* are acquired:

$$y = \Phi \mathbf{x} = \Phi \Psi \mathbf{x} = \Theta \mathbf{x} \tag{1}$$

where the measurement noise was not taken into account. If we use as a projection matrix (noted with Φ) a matrix with dimensions *mxN*, with *m<N*, then it means that we will make a number of m measurements, each measurement of size *N*. That is, the vectors on which *x* is projected represent the rows of the projection matrix.

The main idea regarding (1) is that, because *m* < *N*, the rebuilding of the original signal cannot be realized, but only under the compressibility hypothesis. It has been shown that if Φ and Ψ satisfy certain conditions, the original vector *α* can be obtained as the unique result to the optimization problem:

$$
\mathfrak{A} = \operatorname\*{argmin}\_{} ||a||\_{l\_0} \operatorname\*{subject\, to} \, y = \Phi \Psi \mathfrak{a}, \tag{2}
$$

where *l0* is the (pseudo)norm consisting of the number of nonzero entries of *α*.

The reconstructed signal has the form:

$$\mathfrak{X} = \Psi \mathfrak{X} \tag{3}$$

corresponding to the sparsest representation of *y* in terms of the dictionary ΦΨ. To circumvent the problems of combinatorial nature and noise effect in the case of almost sparse signals, two directions evolved:

$$
\mathfrak{A} = \arg\min\_{\mathfrak{u}} \|\mathfrak{a}\|\_{l1} \text{ subject to } \mathfrak{y} = \Phi \Psi \mathfrak{a} \tag{4}
$$

Let us stress the fact that although pure sparse signals (built of exactly *k*<<*N* atoms from a specified dictionary) are difficult to find, conventional results are valid for signals that are "almost sparse" (which can be built of *k*<<*N* non-negligible atoms) with respect to dictionaries that can be overcomplete (contain more atoms than their intrinsic dimension), as in the case of some classes of biomedical signals. Taking into consideration this fact, it has been found useful to adapt the theory of CS to the field of processing ECG and electroencephalographic (EEG) signals [2–4] as well as for applications [5] such as compression, transmission, reconstruction of ECG signals, ECG filtering and monitoring [6,27,30–32].

For a better understanding of the algorithm, in the following we present a pseudocode summary.

INPUTS: ORIGINAL SIGNAL = *x*

Acquisition Stage:

Step 1: Compute random measurements

*y* = Φ*x*, where Φ is a MxN matrix of random independent and identically distributed (i.i.d.) entries.

Reconstruction Stage:

Step 2: Compute α coefficients using L1 minimization

*α*ˆ = arg*α*min *α l*1 *subject to y* = ΦΨ*α*

Step 3: Reconstruct original signal *x x* = αΨ OUTPUTS: RECONSTRUCTED SIGNAL = *x-*

#### **3. Sample Vectors, Projection Matrices and Dictionaries**

Here, we briefly show several ways of segmentation, a couple of projection matrices, as well as some several ways of building various types of dictionaries specific to the ECG signal. Depending on the chosen CS method, the way of building the dictionary which is used to reconstruct the ECG signal is different.

#### *3.1. Sample Vectors*

First of all, let us mention that we will refer to ECG signals with a sampling frequency according to the Nyquist–Shannon sampling theorem of 360 Hz and 300 (or 301 for case with R-wave centered) samples/vector, respectively. Each vector is projected on a number of random vectors with identic size and the obtained values are utilized for recovering through a dictionary.

In the simplest way, the first 300 samples of the ECG signal set up the first vector; then, the succeeding 300 samples form the second vector, etc. The place of the R-wave can be anywhere in a vector or it may be missing sometimes, which is, obviously, not desirable.

In order to take advantage of the cyclicity of the ECG signal and of the changes produced on the ECG signal in case of some diseases, we proposed some modified acquisition techniques that requires preprocessing [13–16]. Thus, samples of the ECG signal are stored in a buffer zone and a series of preprocessing can be performed on these stored signals. The R-waves can be detected, and based on them, the ECG signal can be segmented into cardiac patterns. A cardiac pattern is delimited by the halves of the RR intervals of two adjacent intervals and re-sampled by interpolation so that the pattern has a fixed number of 301 samples. The above segmentation and preprocessing technique contain simple calculations and, as will be presented, notably increases results for the compression and reconstruction processes.

Starting from the method described above, an improvement of the cardiac model can be obtained by centering the R-wave on sample 151. Thus, a resampling to the left of the R-wave will be performed and another resampling to the right of the R-wave and the final cardiac model will have 301 samples with the R-wave centered. This alignment of the R-wave can be a reversible process, provided that the reduction/stretch ratio from left to right is known. To make a much clearer picture of the re-sampling and alignment effect, we provide in the following examples of unfocused (misaligned) heartbeats and the same cycles prepared to be aligned. These segments constitute atoms in the dictionary or preprocessed sample vectors.

Figure 1 shows examples of cardiac models with and without a centered R-wave.

In conclusion, the sampling vectors and the atoms of the dictionary can be: (i) unprocessed or pre-processed through segmentation and resampling or (ii) segmented and resampled with a centered R-wave.

#### *3.2. Projection Matrices*

A key element in the CS method is the projection matrix for the acquisition of the ECG signal. The reconstruction quality of the ECG signal is considerably decided by the kind of the matrix used in the compression stage [7,9,10,13].

Moreover, the number of random vectors (and respectively the number of calculated scalar products) considered is based on the tolerated tradeoff between the compression ratio and the reconstruction error: thus, the compression ratio is directly related to the reconstruction error.

In Section 3, we will analyze and determine which is a stop ratio and we will determine in the case of our ECG signals how many projections we need for a good ECG compression.

**Figure 1.** Examples of cardiac patterns obtained by centered or non-centered R-wave: (**a**) Cardiac patterns with a centered R-wave; (**b**) Cardiac patterns without a centered R-wave.

In the following, we use and discuss three types of projection matrices.

In this paper, we examine the consequence of these three types of projection matrices on various dictionaries.

#### *3.3. Dictionaries*

Using standard Discrete Cosinus/Sinus Transform (DCT/DST), Wavelet or other typical dictionaries is not always the best choice if we are referring to ECG signal reconstruction errors [15]. Thus, we will analyze the use of dictionaries dedicated to ECG signals, dictionaries that can be specific to the patient, specific to the pathology or universal. The way dictionaries are built is closely related to the segmentation methods of the ECG signals presented above. Thus, concerning the preprocessing stage, we used dictionaries with three types of atoms: (1) Unprocessed (patient-specific only) and processed atoms; (2) Segmented atoms; (3) Segmented plus R-wave centered. The last two types contain either patient-specific beats, or normal beats and/or seven types of pathological beats.

#### 3.3.1. Patient-Specific Dictionaries

In order to build patient-specific dictionaries, we used the first minutes of each patient's record and then the rest of the ECG signal was used for testing. Thus, the atoms represent ECG segments of size 300, successive segments of vectors, without any processing. In our studies, such dictionaries were constructed (only) from the first few minutes of the

patient's records (patient-specific dictionary), the atoms being further used for CS with various projection matrices.

In order to maintain uniformity in the size of the dictionaries, we chose to build patient-specific dictionaries of 700 atoms, each atom having a size of 300. A size of 300 for atoms was determined considering the sampling frequency (360 Hz) and the average beat frequency heart rate (~70 beats/min for normal patients). In this way, the dictionary is actually a matrix with a size of 300 × 700. We highlight that the atoms of the dictionary were aleatory sequences of the ECG recording, and therefore, the R-wave can appear anywhere in the 300 samples or even be missing (not a happy case).

We note that besides the simplicity of the ECG signal segmentation method, another advantage is the capture of the specificity of the patient's ECG particularities in the moment the recording has started.

An improved version of the method is to preprocess the ECG segments to build the dictionary. Thus, segmentation can be performed by detecting heartbeats (i.e., R-waves) and then the R-wave centers. Therefore, patient-specific dictionaries can be constructed without or with preprocessing for R-wave alignment. However, in all cases, the first portion of an ECG recording is used to construct the dictionary, while the rest of the signal (the unused part in the dictionary) was used in the testing techniques.

The next two types of dictionaries contain only atoms obtained through segmentation, normalized to 301 elements with or without a centered R-wave.

#### 3.3.2. Universal Mega-Dictionaries

The mega-dictionary used consists of 1472 atoms (i.e., 184 beats from each of the 8 classes discussed, 7 pathological and the normal beat class). Depending on the preprocessing tested, the atoms of the dictionary may or may not have a centered R-wave.

#### 3.3.3. Pathology-Specific Dictionaries

When the reconstruction stage considers the pathologic class that the cardiac beat belongs to, a particular or specific dictionary has been constructed for each pathological class. Because the ECG recordings include heartbeats from several pathological classes, we tested the variant in which, for each pathological class, we made a specific dictionary. Thus, analyzing 7 pathological classes and the normal class, we built 8 dictionaries, each with 700 atoms specific to each class. Atoms may or may not have a centered R-wave. Thus, we note that the number of atoms in each of the dictionaries is higher than the number of atoms related to a certain pathology contained in the mega-dictionary.

#### **4. Proposed Methods for Dictionary-Based ECG Compression**

In the Introduction, we talked about the presentation of two totally different methods of CS specific to ECG signals, but both have in common the need of building specific dictionaries. However, the use of ECG signal characteristics and how to build dictionaries differ remarkably.

Thus, the PSCCS method is based on ECG signal specific features of each patient, while CPCS on the cyclical patterns of the heartbeat.

In the next subsection, we present two methods for CS of ECG signals with some dissimilarity associated to the projection matrices.

#### *4.1. Patient-Specific Classical Compressed Sensing—PSCCS*

A first variant of compressed acquisition of the ECG signal is presented in Figure 2. It can be implemented even on hardware system and involves the compressed collection of the ECG signal using the CS technique and a patient-specific dictionary together with the Basis Pursuit technique [14].

**Figure 2.** Principle of the PSCCS method.

In this method, the compression of the ECG signal involves the classic use of the CS technique, without any additional signal processing. The advantage of the method is that it speculates on the specific features of the patient. Another advantage is the reduced complexity equal to that of the traditional CS algorithm. The particularity of this procedure is the need for a classic 6-minute ECG acquisition to build the dictionary. In order to obtain improved results, the dictionary can be upgraded in case of long recordings or in case the patient has undergone changes on the ECG signal from one recording to another.

#### *4.2. Cardiac Patterns Compressed Sensing—CPCS*

Below, we present a different approach from the classic CS, which involves a preprocessing stage used both for segmentation of the ECG signal for compressed acquisition and for building useful dictionaries in the signal reconstruction stage.

Figure 3 shows the block diagram of the method. As we can see, at the level of the reconstruction stage there are two approaches, namely, a way of reconstruction using a mega-dictionary or another variant in which dictionaries specific to pathologies are used. The first two operations are common to both approaches and are colored in yellow in the block diagram.

**Figure 3.** Block diagram of the CSCP method, using the mega-dictionary and/or a pathologyspecific dictionary.

The upper branch of the block scheme, colored in green, is for the version with the universal mega-dictionary and the lower part of the figure, colored in blue, is for the version with dictionaries specific to pathologies.

In the case of reconstruction with dictionaries specific to pathologies, it is necessary to know the pathological class to which each cardiac pattern belongs. Therefore, it is necessary to classify the heartbeats. One option is to use a KNN classifier or any other classifier trained with various compressed beats [15,17]. Another option for classifying the heartbeats is a first reconstruction with the mega-dictionary on the upper branch of Figure 3 and the analysis of alpha coefficients corresponding to the mega-dictionary, i.e., the pathological class associated with the heartbeats is the same as the class in which the atom in the mega-dictionary with the highest coefficient belongs at reconstruction with the BP algorithm. Once the pathological class is established, the final reconstruction will be performed with the dictionary specific to that pathology [16].

For the classification of the ECG pattern and the establishment of the dictionary with which the signal will be reconstructed, the KNN classifier trained with the compressed version of the heartbeat from the universal mega-dictionary can be used.

Thus, a first step is to establish the class of the pattern. For this, we will use the KNN classifier based on the highest coefficient corresponding to the mega-dictionary, shown in light blue in Figure 3. Once the membership class is established, the Basis Pursuit algorithm together with the calculation of α coefficients necessary for the reconstruction of the ECG pattern are used. In addition, the almost insignificant distortions due to the centering of the R-wave can be improved by means of the knowledge about the original location of the R-wave.

#### *4.3. Acceptance of the Compression Methods*

To evaluate the compression and reconstruction performances, we assess the distortion between the original and the reconstructed signals by standard *PRD* and *PRDN* measures. Most ECG compression algorithms in the literature evaluate the errors using the percentage root-mean-square difference (*PRD*) measure and its normalized version, *PRDN*, defined as:

$$PRD\% = 100 \left\lceil \frac{\sum\_{n=1}^{N} \left( \mathbf{x}(n) - \tilde{\mathbf{x}}(n) \right)^2}{\sum\_{n=1}^{N} \mathbf{x}^2(n)} \right\rceil$$

and:

$$PRDN\% = 100 \left\lceil \frac{\overline{\sum\_{n=1}^{N} \left(\mathbf{x}(n) - \widetilde{\mathbf{x}}(n)\right)^2}}{\sum\_{n=1}^{N} \left(\mathbf{x}(n) - \overline{\mathbf{x}}\right)^2} \right\rceil$$

where *x*(*n*) and *x* (*n*) are the samples of the original and the reconstructed signals, respectively, *x* is the mean value of the original signal and *N* is the length of the window over which the *PRD* is calculated.

For the evaluation of the compression, we used the compression rate (*CR*) defined as the ratio between the number of bits needed to represent the original and the compressed signal:

$$CR = \frac{b\_{orig}}{b\_{comp}}$$

where *borig* and *bcomp* represent the number of the bits required for the original and compressed signals, respectively.

We also used an alternative measure defined in [19], the Quality Score (*QS*), which is the ratio between the *CR* and the *PRD*:

$$QS = \frac{CR}{PRD}.$$

In addition to the quantitative measure related to the reconstruction of ECG signals, we also used a qualitative evaluation of the signals by classifying them. For classification, we used the KNN classifier. Thus, in the CPCS method version with a pathology-specific dictionary, in order to estimate the signal classification ratio in one of the eight possible classes, we used a KNN classifier to evaluate the reconstructed models. We mention that the KNN was trained with the models from the mega-dictionary used in the classification block (models that were not subjected to compression with the known class for each atom) and tested with the models reconstructed with class-specific dictionaries.

In addition to the KNN classifier, a neural network was used to test the reconstructed signals. The neural network was a multilayer perceptron (MLP) with 10 neurons in the hidden layer with backpropagation gradient descent for training.

However, the final verdict on the fidelity and clinical acceptability of the reconstructed signal should be validated by visual inspection by the cardiologist.

#### **5. Experimental Results**

In this study, we used 24 ECG recordings from the MIT-BIH Arrhythmia database acquired at a sampling frequency of 360 Hz, with 11 bits/sample [18]. Besides the ECG signals, the database also includes annotation files containing the index of the R-wave and the class to which each ECG pattern belongs.

In the CPCS method, we used the annotation databases in the preprocessing step (segmentation of cardiac cycles and forming of dictionaries) and in the reconstructed signal validation phase (KNN classifier-training stage).

The PSCCS technique used only the ECG signals from the MIT-BIH database, without requiring additional knowledge (ECG annotated files).

#### *5.1. Results for the Patient-Specific Classical Compressed Sensing (PSCCS) Method*

To test the PSCCS procedure, we used several compression ratios, namely, 4:1, 10:1 and 15:1. We also used several types of projection matrices (Bernoulli, Gaussian distribution random and dictionary specifics). The data used are 24 records from the MIT-BIH Database. In Table 1, we present the average results for 24 ECG records.

**Table 1.** Average results for 24 ECG records processed with the PSCCS method.

In addition to the average results reported for the MIT-BIH database, a number of authors reported the results for record no. 117 (in Table 2), which is why we will report these results as well.

**Table 2.** Results for the 117 records processed with the PSCCS method.

In Figure 4a, we present a part of the registration no. 117 in the initial version and its version reconstructed following the compression of 4:1, 10:1 and 15:1 for the application of a Bernoulli type projection matrix. It is observed that for CR = 15:1, especially in the noisy region (sample from 2000 to 2200), there are some visible reconstruction differences due to this noise. There are no significant differences in the rest of the signal.

In Figure 4b, we also present from the recording 117 an original ECG signal segment and its variant reconstructed subject to a CR = 10:1 (for random projection matrix with Gaussian distribution). The segment shown is the segment with the highest noise in the entire recording. In this way, we wanted to highlight the robustness of the method to noise and artifacts due to the patient's movement and breathing.

The results obtained on 14 ECG signals, for a compression ratio of 15:1, for centered and non-centered R-wave are shown in Table 3. We used the KNN and MLP algorithm for the evaluation by classification.

**Table 3.** Average results for 14 ECG records with the PSCCS method.

The KNN and MLP classifiers were trained with normal and abnormal heart beats evenly distributed on both classes. The beats used to train the classifier were extracted from the dictionary constructed for the compressed acquisition. In this case, the classification was on two classes, normal or abnormal, and it did not follow the seven pathological classes.

The advantage of the KNN classifier is the simplicity of the calculations, this classifier assuming only the calculation of some Euclidean distances. In the case of MLP networks, the calculations are more complex, but the results are better compared to the KNN classifier.

## *5.2. Results for the Cardiac Patterns Compressed Sensing (CPCS) Method*

#### 5.2.1. Universal Mega-Dictionary

For the construction of a mega-dictionary, from all the 24 ECG recordings, we randomly chose 184 patterns from the 8 cardiac classes, thus obtaining a dictionary with 1472 patterns with the size 1472 × 301.

The testing was performed on 200 patterns from each class, chosen at random from the 24 records, with the mention that special attention was paid to random choice, namely, the models used to build the dictionary could no longer be used for testing.

Table 4 shows the average results obtained on all 24 records, with R-wave alignment and centering and without R-wave centering, for all the projection matrices presented.

**Table 4.** Average results for 24 ECG records processed with the CSCP method with the mega dictionary.

5.2.2. Pathology-Specific Dictionaries

Each of the eight pathology-specific dictionaries is made up of 700 atoms that actually represent patterns with or without a centered R-waves. Dictionaries are matrices of size 700 × 301.

For testing, we used a number of 2000 cardiac patterns chosen at random from the 24 records with the mention that the patterns used for testing are different from those used for training (see Table 5).

**Table 5.** Average results for 24 ECG records for CSCP method with a specific dictionary and classification based on the largest coefficient of the sparsest decomposition for the mega-dictionary.

In this variant, with dictionaries specific to the pathological class, in the reconstruction stage, it is necessary to identify the class to which the pattern belongs. The reconstruction results are strongly influenced by the correctness of establishing the pathological class to which the model belongs. Thus, for patterns classification, a KNN type classifier will be used or it will be made based on the highest alpha coefficient. Once the pathological class is established, the Basis Pursuit algorithm, the dictionary specific to that pathology and the projection matrix will be used for reconstruction.

Thus, using the classification of patterns based on the highest alpha coefficient in the mega-dictionary version, a pattern classification rate of 88.75% is obtained [16]. Using the KNN classifier with training on 1472 compressed cardiac patterns (uniformly distributed in the eight classes), a classification rate of 93.77% is obtained [15].

In Figure 5, we present examples of reconstructed cardiac beats for every pathology class.

Qualitative estimation of reconstructed signals based on classification. In addition to the quantitative measures of the distortions between the original and reconstructed ECG signals, for a further verification of the quality of the proposed compression scheme, we performed a classification of reconstituted models with the KNN algorithm. The classifier was trained with the atoms from the mega-dictionary. A first check of the method is to test the performance of the KNN classifier, and for this, we initially tested the original models (i.e., the uncompressed models that we used to test the compression scheme). For these patterns, we obtained a classification rate of 93.75%. The results presented below are obtained on the reconstructed patterns [28].

In addition to KNN, an MLP classifier was also tested. This second classification aims to strengthen the correctness of the idea of testing the reconstructed patterns from a qualitative point of view. This test is based on a classifier and is needed to compare the results obtained with these two different classifiers. Thus, there is a slight and almost insignificant improvement of the classification rate in the case of MLP compared to KNN. However, in practical implementations, the MLP classifier should be chosen according to the available hardware resources. Table 6 shows obtained results for dictionaries with a centered R-wave.

**Figure 5.** Original and reconstructed signals with pathology-specific dictionaries.

**Table 6.** Results summary for dictionaries with a centered R-wave.

It is known that in a classification process, especially when it applies to several classes, special attention must be assigned to the confusion matrix, to see if the classification is uniform on all classes or only certain classes are detected. For this we have exemplified in Table 7 a confusion matrix for the classification variant with a mega-dictionary. It can be seen that the classification rate is evenly distributed over all eight classes.

**Table 7.** Confusion matrix for KNN classification of the reconstructed patterns with a mega-dictionary.

5.2.3. Patient-Specific Dictionaries

The patient-specific dictionaries were constructed from the patient's first 700 heartbeats, and preprocessed as previously described (i.e., with or without R-wave alignment). Thus, the dictionary is made up of 700 atoms, each of size 300, i.e., it is actually a matrix of size 301 × 700. This method has the advantage of speculating quasi-periodicity and the particular characteristics of the ECG signal of a particular patient. Table 8 shows average results for 24 ECG recordings for the CSCP method and it can be seen that the best results are obtained if we refer to QS.

**Table 8.** Average results for 24 ECG Records for the CSCP method with a patient-specific dictionary built from the first 700 cardiac cycles.

Because our results are generally obtained by mediating the results obtained by processing 24 records from MIT-BIH Arrhythmia database, we present in Figure 6 the histograms of PRD and PRDN, respectively, for the method of CS with patient-specific dictionaries with a centered R-wave and projection matrix by type Gaussian distribution Random \* Dict †. For this case, PRD\_average = 0.51 and PRDN\_average = 9 (see Table 8).

**Figure 6.** Histogram of PRD and PRDN for 24 ECG records for the CSCP method with a patientspecific dictionary with projection matrix by type of Gaussian distribution Random \* Dict †.

#### **6. Discussions**

In Table 9, we resume the results previously presented for the two analyzed methods, for a CR = 15:1 with all investigated projection matrices and with all discussed reconstruction and preprocessing dictionaries. We marked in bold the best results obtained on QS (Quality Score) for each method.

It can be seen that the best QS result is obtained for dictionary specific to the patient in which the R-wave is centered and a projection matrix is optimized to the dictionary. In addition, it has been also found that in all cases optimization of dictionaries improves the results. Moreover, it has also been observed that preprocessing improves the results, namely, for PSCCS (i.e., without preprocessing) for CR = 15:1 the best QS equals 15.46, i.e., almost half of the value obtained with CPCS with a patient-specific centered R-wave dictionary when QS = 29.13.

It should be noted that any preprocessing means hardware resources and choosing a method with preprocessing means additional hardware resources. However, we must mention that the detection of the QRS complex and the R-wave is a problem that can be implemented in real time in the Matlab® environment, an example of implementation being even available in Help Matlab® [29].

Table 10 shows the average results on the 24 records and for the 117 records obtained by other authors.

**Table 10.** Average values for 24 records and 117 record for other compression algorithms.

We note that Mamaghanian in [22] presents a classical CS compression method followed by Huffman coding, the final CR being higher due to the additional Huffman compression. For a more accurate comparison, we must compare our results with those obtained by Mamaghanian before Huffman compression. Additionally, the same author uses in [22] the compression ratio defined as:

$$CR = \frac{b\_{orig} - b\_{comp}}{b\_{orig}} \* 100\%$$

which is not the same as ours and gives a very different gamut of values compared with ours.

The results we obtained with the proposed method are compared in Table 11 with the results of other compression methods in the literature.

NOTE: The results reported in [26] marked with \* in Table 11 were obtained using a combined ECG compression method consisting of a preprocessing stage with quad level vector (QLV) for the extraction of the ECG skeleton achieving an 8.4:1 compression and a coding block (consisting of delta and Huffman Coding). The results referenced in Table 3 are the final one, improved by the Huffman coding stage.

#### **7. Conclusions**

The results presented in this paper reveal several interesting aspects, as follows.

It has been revealed that the first stage of the CS method, i.e., the signal acquisition part, based on the projection matrices, has only a relatively small influence on the decompression or classification results.

On the other hand, for the second stage, namely signal reconstruction, the dictionary used for reconstruction of the compressed sensed ECG signals has an essential role in obtaining good results. Therefore, depending on the application targeted with the used CS technique, namely, Holter monitoring or recorded ECG signal classification, a dictionary that leads to optimal final results can be selected.

Thus, in a Holter monitoring application, where the ECG signal is recorded for 24 h from the same patient, one can choose the Patient-Specific Classical Compressed Sensing (PSCCS) method. By analyzing the first minutes of the recording, a dictionary specific to the patient will be built, and then it will be used to reconstruct the ECG segments of interest to the specialist.

Otherwise, if the CS-based application aims at classifying heartbeats for ECG monitoring or abnormality identification applications, the Cardiac Patterns Compressed Sensing (CPCS) method will be chosen, where each pathological heart beat class will be associated to a specific dictionary.

The above discussed methods are primarily based on waveform segmentation (cardiac beats) with no preprocessing. Yet, depending on the available hardware resources and the time constraints in which the application should run, the results can be significantly improved by centering the R-wave using ECG preprocessing i.e., segmented cardiac patterns with a centered R-wave.

This choice is related to the idea that any ECG signal preprocessing leads to higher hardware requirements and slowdowns in the acquisition and reconstruction processes over time. However, these aspects can be easily dealt with, aiming at better results. However, we must mention that the detection of the QRS complex and the R-wave is a problem that can be implemented in real time in the Matlab® environment, with an example of implementation even being available in Help Matlab®.

**Author Contributions:** Conceptualization, M.F., H.-N.C. and L.G.; methodology, M.F.; software, M.F.; validation, M.F., H.-N.C. and L.G.; formal analysis, M.F.; investigation, M.F.; resources, M.F.; data curation, M.F.; writing—original draft preparation, M.F., H.-N.C. and L.G.; writing—review and editing, M.F., H.-N.C. and L.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in [physionet] at [10.1109/51.932724 and 10.1161/01.cir.101.23.e215], reference number [18]. The webpage of the MIT-BIH Arrhythmia Database is "https://www.physionet.org/content/mitdb/1.0.0/" (accessed on 27 December 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Machine Learning Based Lens-Free Shadow Imaging Technique for Field-Portable Cytometry**

**Rajkumar Vaghashiya 1,†, Sanghoon Shin 2,†, Varun Chauhan 1, Kaushal Kapadiya 1, Smit Sanghavi 1, Sungkyu Seo 2,\* and Mohendra Roy 3,\***

**Abstract:** The lens-free shadow imaging technique (LSIT) is a well-established technique for the characterization of microparticles and biological cells. Due to its simplicity and cost-effectiveness, various low-cost solutions have been developed, such as automatic analysis of complete blood count (CBC), cell viability, 2D cell morphology, 3D cell tomography, etc. The developed auto characterization algorithm so far for this custom-developed LSIT cytometer was based on the handcrafted features of the cell diffraction patterns from the LSIT cytometer, that were determined from our empirical findings on thousands of samples of individual cell types, which limit the system in terms of induction of a new cell type for auto classification or characterization. Further, its performance suffers from poor image (cell diffraction pattern) signatures due to their small signal or background noise. In this work, we address these issues by leveraging the artificial intelligence-powered auto signal enhancing scheme such as denoising autoencoder and adaptive cell characterization technique based on the transfer of learning in deep neural networks. The performance of our proposed method shows an increase in accuracy >98% along with the signal enhancement of >5 dB for most of the cell types, such as red blood cell (RBC) and white blood cell (WBC). Furthermore, the model is adaptive to learn new type of samples within a few learning iterations and able to successfully classify the newly introduced sample along with the existing other sample types.

**Keywords:** artificial intelligence; lens-free shadow imaging technique; cell-line analysis; cell signal enhancement; deep learning

#### **1. Introduction**

The lens-free shadow imaging technique (LSIT) is a well-established technique for the characterization of microparticles and biological cells [1]. This technique is widely popular for its simple imaging structure and cost-effectiveness. It comprises a lens-less detector, such as a complementary metal-oxide semiconductor (CMOS) image sensor, a semi-coherent light source, such as light-emitting diode (LED), and a disposable cell chip (C-Chip). The absence of a lens or other optical arrangements allows it to fit into a very small space, thereby reducing the size of the overall system (as described in Figure 1a in the LSIT platform (Cellytics) built within a dimension of 100 × 120 × 80 mm3). Since this arrangement consists of a few components, most of which are easily available at a low price, it therefore reduces the overall cost of the system [2]. This simple and cost-effective

**Citation:** Vaghashiya, R.; Shin, S.; Chauhan, V.; Kapadiya, K.; Sanghavi, S.; Seo, S.; Roy, M. Machine Learning Based Lens-Free Shadow Imaging Technique for Field-Portable Cytometry. *Biosensors* **2022**, *12*, 144. https://doi.org/10.3390/ bios12030144

Received: 30 January 2022 Accepted: 25 February 2022 Published: 27 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

nature facilitates the feasibility of the LSIT for the applications in the fields of point-of-care systems or telemedicine systems [3–5].

**Figure 1.** Schematics of the LIST setup and the proposed neural network architecture for the auto characterization of LSIT micrographs. (**a**) LSIT platform (Cellytics) (**b**) schematic of the principle of diffraction, i.e., shadow, pattern generation of a micro-object, (**c**) schematic of the LSIT imaging setup showing the simplicity of the setup, (**d**) schematic of the dataset creation process by automatic cropping individual cell diffraction pattern from the whole LSIT micrograph, and (**e**) the schematic of the proposed denoising and classification architecture. Here, the denoising autoencoder enhances the signal of the individual cells which is then fed to the CNN module for classification.

Recent advancements in machine learning, especially deep learning, have facilitated many applications concerning medical diagnostics [6–12], and have been widely adopted in the field of microscopy [13–15]. In particular, deep learning has been incorporated with the LSIT [14], where it is has been used to enhance the resolution of the LSIT micrographs [16] and enabled polarization-based holographic microscopy [17].

In our previous work, we have successfully developed the LSIT imaging system for the complete blood count using an analytical model based on handcrafted features [3] that can automatically segment out the individual cells from a whole frame LSIT micrograph and subsequently analyze them based on the handcrafted parameters. However, the performance of the system is dependent on the uniform illumination as well as the strong signatures of the microparticle samples. Since the diffraction signature of a microparticle depends on the size as well as the signal-to-noise ratio of the particle, therefore any background noise can affect the overall performance of the auto characterization system. Further, the handcrafted approach of finding the features for every additional cell line is time-consuming and prone to subjective errors. To address these limitations, in this work, we have developed an artificial intelligence (AI) powered signal enhancement scheme for the LSIT micrographs that can enhance the signal quality (signal to noise ratio (SNR)) for various cell lines in a heterogeneous cell sample. For this, we employed the autoencoderbased denoising scheme [18]. Further, we have developed an auto characterization method based on a convolutional neural network [19,20] (CNN) architecture to classify the various

cell lines from the LSIT micrograph. Here, we have first introduced the transfer of learning scheme in a neural network, which can leverage the feasibility to introduce new cell types to the algorithm and thus learn their characteristics within a few iterations. Thus, the LSIT platform saves time and computation resources required to learn to classify the additional cell types along with the existing ones.

In this article, we have described the detailed methods adopted for the designing as well as optimization of various parameters to design a suitable model with better accuracy. These optimized models are simple and light-weight, and require a smaller number of samples for effectively learning the cell signatures. The details are as given in the following sections.

#### **2. Methods**

#### *2.1. LSIT Imaging Setup*

The schematic of our proposed setup (Figure 1a) is as shown in Figure 1. When light from the coherent or semi-coherent source passes through a micro-object, it produces characteristic diffraction, i.e., shadow, pattern of the object [21,22] as shown in Figure 1b. These diffraction patterns are prominent just beneath the sample, typically a few hundred micrometers away from the sample plane, from where they are captured using a highdensity image sensor such as CCD or CMOS [21] (Figure 1c). As these signatures are significant enough to be captured by the bared image sensor, it does not require any kind of lens arrangement [23]. In our proposed setup, we used a pinhole conjugated semicoherent LED light source with a peak wavelength of 470 ± 5 nm (HT-P318FCHU-ZZZZ, Harvatek, Hsinchu, Taiwan). The diffraction patterns were captured using a 5-megapixel CMOS image sensor (EO-5012M, Edmund Optics, Barrington, NJ, USA), and a customdeveloped C-Chip (Infino, Seoul, Korea) was used to hold the cell samples [2,5]. All of these components can fit in a compact dimension of 100 mm × 120 mm × 80 mm. Due to the absence of a lens-based setup, the field-of-view of this system is about 20 times that of a conventional optical microscope at 100×. This high-throughput nature provides an extra advantage to characterize several thousand cells within a single digital frame.

#### *2.2. Preparation of Various Cell Lines*

In this work we used various cell lines, starting from red blood cell (RBC), white blood cell (WBC), cancer cell lines HepG2 (human liver cell-line) and MCF7 (human breast cancer cell-line), and polystyrene microbeads of 10 μm and 20 μm. The preparations of these cell lines are as follows [2,3]. The use of human whole blood in the experiment was approved by the Institutional Review Board (Approval No. # 2021AN0040 of Korea University Anam Hospital (Seoul, Korea).

**RBC:** The RBC samples were prepared from the whole blood samples that were collected from the Korea University Anam Hospital under IRB approval. The samples were diluted about 16,000 times by using RPMI solution (Thermo Scientific, Waltham, MA, USA) [2,3].

**WBC:** First, Ficoll solution (Ficoll-Paque™ Plus, GE Healthcare, Chicago, IL, USA) was used to isolate mononuclear cells from the whole blood. The samples of peripheral blood mononuclear cells (PBMCs) obtained using the Ficoll solution are mixtures of lymphocytes and monocytes. To separate these two cell types, the MACS (Magnetic-activated cell sorting) device and antibodies (Miltenyi Biotec, Bergisch Gladbach, Germany) were utilized. The helper-T cells in the lymphocytes were separated using the CD4 antibody (#130-090-877), and the cytotoxic-T cells with the CD8 antibody (# 130-090-878). Finally, 10 μL of this solution was then loaded into the unruled C-Chip cell counting chamber [2,3].

**HepG2:** The HepG2 cell lines were prepared from the American Type Culture Collection (ATCC HB-8065) and incubated in a high-glucose medium (DMEM, Merck, Darmstadt, Germany) with 10% heat-inactivated fetal bovine serum, 0.1% gentamycin, and a 1 penicillin/streptomycin solution under 95% relative humidity and 5% CO2 at 370 ◦C. The developed cells were then trypsinized and separated from 24 well pate and incubated from 2–5 min at 370 ◦C. These cells were then diluted with DMEM solution [2,3].

**MCF7:** The MCF7 cell samples were prepared from the American Type Culture Collection (ATCC HTB-22). The cells were preserved in a solution of DMEM containing 1% penicillin/streptomycin solution, 0.1% gentamycin, and 10% calf serum at 95% relative humidity and 5% CO2 at 370 ◦C. These cells were then trypsinized and separated from the 24 well pate. These separated cells were then incubated for 2–5 min at 370 ◦C. The cells were then washed with DMEM solution. 10 μL of this solution was then loaded in the C-Chip [2,3].

**Polystyrene microbead:** The 10 μm and 20 μm bead samples were prepared by diluting the respective polystyrene microbeads (Thermo Scientific, Waltham, MA, USA) with de-ionized water [2,3].

#### *2.3. Dataset Creation*

A whole frame LSIT image (of cell diffraction patterns) contains an average of ~500 diffraction patterns of microparticles. Deep learning-based architectures utilize the features of each class, and typically require a minimum of a few hundred diffraction patterns of each cell type for optimal learning. Therefore, we cropped individual diffraction patterns of each cell type (that were verified using a traditional microscope) with a window of 66 × 66 pixels as shown in Figure 1d. This window size included the complete sample signature along with a minimal background that would provide complete cell-line information during the auto-feature selection process in learning algorithms. We further augment this base sample set by rotating the individual diffraction patterns with an increasing angle of 10 degrees clockwise. Finally, a dataset of 1980 samples for each of the six cell lines and microparticles was created, totaling 11,880 samples for all of the classes under study. The typical architecture of a CNN is illustrated in Figure 1e. As many learning algorithms are black-box models, it is difficult to ascertain the optimal cell-signal size that covers the majority of the information and minimal background. Naturally, a smaller cell-size would need lesser computation and have lesser noise. Hence, we created the dataset for 60 × 60, 56 × 56, 50 × 50, 46 × 46, 40 × 40, and 36 × 36 cell sizes as input sets, with each set further divided into training and test folds. As the data augmentation used sample rotation, the splitting of the dataset into the train, validation, and test folds needs to be carried out while keeping a check on data leakage. Augmented samples distributed across the train and test sets may bias the model and may give a wrong estimate of its performance as the test data may not be of entirely "unseen" samples. Accounting for this, the 1980 samples of each class were carefully split into 1490 training samples, 166 validation samples, and 324 testing samples.

Though the cell-lines may seem visually similar, there are significant differences in the statistical distributions of the pixel illumination intensity in the cell diffraction pattern as revealed in our exploratory data analysis. The 2D contour plots (in Figure 2) show the observed variances. Hence, it is possible for intelligent algorithms to automatically identify and utilize the descriptive features in signal enhancement as well as classification.

#### *2.4. Denoising Modality*

For denoising of the LSIT micrographs, we adopted the concept of autoencoder [24]. An autoencoder is an unsupervised scheme that scuffles to recreate the input at its output. It consists of an input layer (x), an output layer (r), and a hidden layer (h). The hidden layer h termed as a code layer stands for the input in a revised dimension. The whole network structure can be labelled into two parts. The first part is an encoder, which tries to code the input as h = f(x), and the second part is a decoder which tries to recreate the input from the reduced code layer as r = g(h), where r is the recreated assortment of input x. Basically, it tries to attain r = g(f(x)). However, this is not a linear transformation since the model is enforced to learn the significant features of the input to encode it into the code layer.

In this work, we specifically used the denoising version of the autoencoder. Traditionally, the autoencoders try to reduce the loss as L(x, g(f(x))). However, the denoising autoencoder attempts to reduce the cost as L(x, g(f(x'))) where x' is the noisy form of the input x. We tried two different methods to design the denoising architectures, namely, extreme learning machine (ELM) and convolutional neural network (CNN).

#### 2.4.1. ELM

This is a single hidden layer fully connected architecture [25]. In this method, the input weights are initiated randomly and kept intact. Only the output weights take part in the learning process through a straightforward learning method [25–27]. For N arbitrary input samples xi ∈ Rn and their counterpart targets ti ∈ Rm, the ELM achieves this mapping using the following relation as shown in Equation (1).

$$H\beta = T\tag{1}$$

Here, *H* is the hidden layer output matrix, *β* is the output weight matrix (i.e., between the hidden layer and the output layer) and *T* is the target matrix or matrix of desired output [26]. From Equation (1), we can obtain the *β* using Moore–Penrose pseudoinverse [25] as shown in Equation (2).

$$\mathcal{J} = \left(H^T H\right)^{-1} H^T T \tag{2}$$

In the extended sequential learning form of ELM, the *β* can update sequentially. This provides an added advantage of updating the learning whenever a new type of sample is available, thus providing the flexibility of transfer of learning. The *β* update mechanism [28,29] is as shown in Equation (3).

$$\beta\_n = \beta\_{n-1} + P\_n^{-1} \left( \ T\_n - H\_n \beta\_{n-1} \right) H\_n^T \tag{3}$$

Here

$$P\_n = P\_{n-1} + H\_n H\_n^T \tag{4}$$

For *n* = 1,

$$P\_{n-1} = P\_0 = \left(\frac{1}{\mathcal{C}} + H\_0 H\_0^T\right) \tag{5}$$

Here *H*0 is the hidden layer output with the first sample or first batch of samples [30].

#### 2.4.2. CNN

Convolutional neural networks (CNN) [6,19,20] are a type of neural network widely used in the analysis of spatial data such as image classification and object segmentation. In this network, two-dimensional kernels are used to extract the spatial features from the input patterns, using a convolution operation between the kernel and the input. The typical architecture of a CNN is as shown in Figure 1e. Here, the kernel is shared spatially by the input or by the feature map. The feature at the location (*i*, *j*) in the kth feature map of the lth layer can be evaluated as shown in Equation (6).

$$Z\_{i,j,k}^{l} = \left(W\_K^l\right)^T X\_{i,j}^l + b\_k^l \tag{6}$$

Here, *Wl K* and *bl k* are the weights and the bias vector of the kth filter in the lth layer. Here the weight layer is shared spatially which reduces the complexity. *Xl i*,*j* is the value of the input at location (*i*, *j*) of the lth layer. The nonlinearity in this network can be obtained by introducing the activation function, denoted here as *g*(.). The activated output can be represented as shown in Equation (7).

$$a\_{i,j,k}^l = \lg\left(Z\_{i,j,k}^l\right) \tag{7}$$

Additionally, there are pooling layers that introduce shift-invariance by reducing the resolution of the activated feature maps. Each pooling layer connects the feature map to the preceding convolutional layer. The expression for pooling is as shown in Equation (8).

$$y\_{i,j,k}^{l} = P\left(a\_{n,m,k}^{l}\right), \forall (m,n) \in \mathcal{R}\_{ij} \tag{8}$$

Here *P*(.) is a pooling operation for the local neighborhood *Rij* around the location (*i*, *j*). In this work, we used CNN for both denoising as well as classification. The details of their architectures and their impacts are discussed in the Section 3.

#### **3. Results and Discussion**

#### *3.1. Performance of Denoising Algorithms*

For efficient and adaptive denoising, we analyzed various autoencoder schemes, starting with the fully-connected autoencoder. In our first iteration, we experimented with the fully connected network having three hidden layers with 512, 256, and 512 neurons, respectively. The input layer is the 1D vectorized array of the input cell diffraction pattern, e.g., of 66 × 66 pixels. The input to the model is the noisy version of the input cell diffraction pattern and the expected target output is the original cell diffraction pattern. The noisy cell diffraction patterns were created using a Gaussian distribution with variance ranging from 100 to 600 with zero mean (refer to the supplementary section for detail). Further, we experimented with an increased network size having five hidden layers with 256, 128, 64, 128, and 256 neurons, respectively. In all of these networks, rectified linear unit (ReLU) was used as the activation function while mean squared error (MSE) [31–34] was used to calculate the loss. The Adam optimizer [35,36] was found to deliver better convergence and hence used to perfect the weight and biases. The denoising performance was quantified in terms of the improvement in *SNR*, measured in dB, denoted here by *SNRimp,* as given by Equation (9) [37].

$$SNR\_{imp} = SNR\_{out} - SNR\_{in} \tag{9}$$

where *SNRout* = 10*log*10 ∑*N n*=1 *x*2 *i* ∑*N n*=1(*x*ˆ*ι*−*xi* ) 2 , and *SNRin* = 10*log*10 ∑*N n*=1 *x*2 *i* ∑*N n*=1(*x ι*−*xi* ) 2 . Here *xi* is the value of sampling point *i* in the original LSIT signal, *x ι* is the value of sampling point *i* in the noisy LSIT, and *x*ˆ*ι* is the value of sampling point *i* in the denoised version of the same cell diffraction pattern. *N* is the total number of sample points in that LSIT image (cell diffraction patterns).

The fully connected network for both the above configuration shows no significant improvement in *SNRimp* after reaching saturation at around −10.08 dB. For further improvement, we experimented with CNN architecture using various models with a different number of convolution layers and distinct kernel sizes. The configuration of the model which accomplished the best outcomes is 3 × 3, 3 × 3, 5 × 5, 5 × 5, 7 × 7, 7 × 7, 1 × 1 with 32 filters in each layer except the last layer. The last layer consists of a single pixel filter (1 × 1 filter) that is used to condense the output across all the 32 filters. Here, the input and output size are the same. Padding was used to maintain the original size after the output of each convolutional layer. The Adam optimizer was used to optimize the network to reduce the mean squared error loss. The CNN results show a better reconstruction as shown in Figure 2.

The CNN network has been optimized for various parameters. First, the optimization of the network for various design parameters, such as varying the convolution layers and the kernel sizes, was carried out. The results in Figure 3a show that the architecture with kernel sizes 3 × 3, 3 × 3, 5 × 5, 5 × 5, 7 × 7, 7 × 7, 1 × 1 has a better performance in terms of *SNRimp*. The performance of the optimized network for various noise parameters, as shown in Figure 3b, indicates the network performs better reconstruction with increasing noise variance in the image (cell diffraction pattern). An increase in the variance results in a noisier image (cell diffraction pattern), which warrants a detailed reconstruction to reverse it to the original form, and hence larger the value of *SNRimp*. Therefore, a higher improvement in *SNRimp* implies the network has learned the optimal representational features for the cell types which enables it to perform a better qualitative reconstruction. Figure 3c compares the reconstruction performance of the model on different sizes of the input image (cell diffraction pattern). Due to the black-box nature of deep learning methods, we had to create datasets with multiple cell-signature dimensions, such that the smallest size just covered the central signature of the cell and increased the window size till it covered a significant background portion as well. The models were evaluated across

varying cell sizes to determine the optimal signal to background ratio, the spatial extent up to which the models covered the features, and to study its effects on the model performance. This analysis is critical in understanding the model explainability and interpretability since having a size larger than the optimum increases the inclusion of background artifacts that affect denoising as well as overpower the cell signal while having a smaller one could exclude the important deterministic features of the cell signature. The convergence in the training phase of the network is as shown in Figure 3d. The results depict that the loss across the first epoch, with a high variation in the initial phase, gets smoother towards the end of the first iteration. The advantage of this system is that it generalizes well for all of the types of cell lines using the same model.

**Figure 2.** Reconstructed results from the optimized CNN. The top row is the original LSIT image (cell diffraction pattern) of a single RBC, WBC, MCF7, HepG2, 10 μm beam, and 20 μm bead. The second row is the noisy version (with variance 100) of the corresponding original images (cell diffraction pattern). The third row is the denoised version of the corresponding original images (cell diffraction pattern) from the noisy image (cell diffraction pattern). The fourth row is the 2D intensity contour plot of the original image (cell diffraction pattern) to show the unique signature of each of these cell lines.

Further, we tried the ELM architecture which is well known for its fast convergence [25]. The results in Figure 3e–h show the performance of the ELM architecture with varying number of neurons in the hidden layer. As it can be concluded from Figure 3e, the model with 2000 neurons provides better performance in terms of *SNRimp*. Further, the optimized model has been used to test the performance across various noise levels as shown in Figure 3f. It is observed that the model maintains the *SNRimp* value on increasing the noise in the input image (cell diffraction pattern), i.e., the image (cell diffraction pattern) quality of the output relative to the input remains the same. The results of the model performance across different sizes of the input image (cell diffraction pattern), as shown in Figure 3g, indicate that the 40 × 40 is having a higher value of *SNR*. However, the variation is of 2 as compared to the variation for the size 50 × 50, which is of 1.5, representing the lowest compared to all of the other sizes. Since CNN shows a substantial performance with lower variance for the 50 × 50 input size, therefore we fixed it as optimal for all of the further studies and comparisons. Figure 3h shows the loss across the first epoch for ELM which is remarkably high initially but converges faster, as compared to CNN, after training with only a few thousand samples. This faster convergence may help save time and resources during incremental training phases for newer cell types. The performances of these optimized models have been compared with the traditional denoising methods as shown in Table 1 (The comparative visual reconstructions are provided in the supplementary document). It is concluded from the data that CNN shows better performance

compared to the other modalities (The details of the traditional methods are provided in the supplementary document). Therefore, we prefer to use CNN for denoising.

**Figure 3.** Results from the CNN and ELM autoencoder. (**a**) The performance of the CNN autoencoder for improved *SNR* (average of all of the classes) across various layers and kernel sizes. (**b**) The performance of the CNN autoencoder across varying noise levels. Here, the variance of the Gaussian noise ranges from 100 to 600, and is evaluated on the optimal network architecture, i.e., 3,3,5,5,7,7,1. (**c**) The performance of the CNN autoencoder across various input sizes (cropping size). Here, the sizes vary from 66 × 66 to 36 × 36. (**d**) The convergence of the optimal CNN network with the number of samples for the first epoch. (**e**) The number of hidden layer neurons in ELM autoencoder vs. improvement in *SNR*, (**f**) variance vs. improved SNR for ELM autoencoder, (**g**) Input size vs. improved SNR for ELM autoencoder, and (**h**) The convergence of ELM autoencoder within the first epoch.

**Table 1.** Comparison of improved *SNR* with respect to the variance for various denoising modalities (here we keep input and output size as 50 × 50).

#### *3.2. Performance of Classification Algorithm*

Since the diffraction patterns of cells and microparticles in a LSIT micrograph depend upon their physical and optical properties, therefore, the diffraction patterns carry the unique signatures of each of the cell types as shown in the 2D contour plot in Figure 2. These unique signatures can be utilized for the classification of these cell types. Since our previous inference concludes that CNN works better for denoising, therefore we experimented with the same modality for the classification as well. In this work, in order to determine the optimal architecture of CNN for cell-line recognition, we first proceeded to find the optimal depth of the network by studying the classification performance of the model on increasing the depth, by adding convolutional and pooling layers, as well as by varying number of kernels and kernel size, till we reached performance saturation. We have experimented with and evaluated various shallow and deep CNN models to classify cell lines. The details of the model architecture are as described in Figure 4.

**Figure 4.** The CNN Architectures with varying depth. The models in the solid box are for the optimization of the model with varying depth. The models inside the dotted box are for the optimization of the parameters. Here, the orange line is the input layer and the rectangle represents the 2D convolution layer. The kernel size and number of kernels are indicated inside the round brackets, the dropout rate in square brackets, and the number of neurons in fully connected layers are placed directly within the rectangle. The aqua blue, green, and purple lines are the max-pool layers. The red line is the output layer and uses a SoftMax activation function. The final optimized architecture is as represented by the last 3D figure.

The Deep Model starts with a convolutional (Conv2D) layer having 512 kernels of size 3 × 3, followed by a max-pooling layer of the same kernel size. The output from this is further convoluted with 128 kernels of 3 × 3 size with a dropout rate of 0.5, and then a max pool with 2 × 2 kernel. This output goes to a Conv2D layer with 64 of 3 × 3 sized kernels and a dropout rate of 0.2. We further reduce the dimension using a 2 × 2 max pool kernel. The output of this layer further convolves with 32 of 3 × 3 kernels, then a dropout of 0.2. The output dimension from this convoluted layer is further reduced by using the max pool with a 3 × 3 kernel. This again convolves with 16 of 3 × 3 kernels, and a max pool layer with 3 × 3 kernel. The output of which is again convoluted with 8 of 3 × 3 kernels, followed by a 3 × 3 max pool. This output is then vectorized and input to a fully connected (FC) layer having 256 nodes, and then to another FC with 128 nodes and having a dropout of 0.2. The final layer is a SoftMax function, with six output nodes. The model architectures used for studying the impact of the network depth and breadth on performance are well described in Figure 4. Once the approximate optimal depth and breadth had been determined, we proceeded to fine-tune the hyper-parameters such as the number of kernels, kernel size, and dropouts, to reach the best performance of the models across varying cell sizes (i.e., input dimension of cells). In all of the models, Adam [38] provided better convergence as compared to other optimizers and has been used as the model optimizer, with categorical cross-entropy [39,40] as a loss estimator. The result for depth and breadth optimization indicates that the average accuracy of the intermediate model is >0.85 (including all of the classes), whereas it is <0.85 for shallow and deep networks. Therefore, we proceed with this intermediate model as an optimum model for further study. Considering the intermediate model has the optimum breadth and depth, the further optimization of the parameters and the results are as shown in Figure 5a.

**Figure 5.** Results for the optimization of the CNN model. (**a**) Performance of the fine-tuned intermediate architecture. (**b**) The confusion matrix showing classification accuracy on the test dataset across all cell lines using the fine-tuned and optimized model. (**c**) The receiver operating characteristic (ROC) curve for each of the cell lines for the optimized model.

From Figure 5a, it is inferred that Model 3 shows better classification performance, on the validation fold, of all of the models. The results depict that there is consistency in performance for the input sizes 40 × 40 to 66 × 66, with a very small variance in the accuracy. The performance of this optimized model is further evaluated over the test dataset containing 324 samples of each cell type. The per-class performance of this model is shown in the confusion matrix of Figure 5b. The results depict that the model can classify

RBC, WBC, 10 μm, and 20 μm bead with over 99% accuracy. However, the comparatively poor performance of about 90% for the cancer cells, HepG2 and MCF7, can be attributed to the non-homogeneity in their signature characteristics as well as the lack of sufficient original samples which further complicates the issue. This is well depicted in our previous work [4] (see Figure 2 of the reference). From the receiver operating characteristic (ROC) curve for all of the cell lines shown in Figure 5c, the area under the curve (AUC) for all of the cell lines is >0.99, except MCF7 (~0.95) and HepG2 (~0.96). From these results, it can be inferred that the classifier is working well, especially for RBC, WBC, 10 μm, and 20 μm beads. The visualization of the internal activation maps, as shown in the Supplementary Information, implies that the network is learning core descriptive features in the diffraction signatures rather than using some random features.

The performance evaluation of the proposed Al model with various matrices such as true positive (TP), true negative (TN), false positive (FP), false negative (FN), accuracy, recall, specificity, sensitivity, F1 score, positive predictive value (PPV) and negative predictive value (NPV) is summarized in Table 2.

**Table 2.** Performance evaluation scores of the proposed AI model using various metrics for different cell types.

Additionally, we also investigate the transfer of learning to gauge the ability of the trained network to adapt to newer cell types (Figure 6). For this scenario, the CNN was initially trained with all of the cell lines except RBC. From the epoch vs. accuracy graph in Figure 6a, the transfer training achieved higher accuracy with the same number of epochs compared to the initial training. This is also validated by the epoch vs. loss graph in Figure 6b. From these results, it can be inferred that the network can be effectively used to adapt to newer cell lines with very less amount of training. From the per-class test accuracy shown in Figure 6c, it is observed that the model misclassified all of the RBC samples as WBC. In the transfer of the learning phase, the initially trained network is frozen except for the last layer, which is modified to accommodate the newer class and kept trainable. The network is then re-trained with a mix of RBC samples. The per-class test accuracy of the re-trained model is shown in the confusion matrix of Figure 6d, where it is inferred that the re-trained network can classify RBC correctly with substantial accuracy.

The comparison between the proposed AI method and the manual method for counting various cell types from a heterogeneous sample is shown in Figure 7. This comparison shows the robustness of the model.

**Figure 6.** The results from the transfer of learning study. (**a**) Epoch vs. accuracy graph for the initial epochs in the initial training phase without RBC and then the transfer of training with RBC. The blue and saffron colors represent the initial training accuracy and validation curves. The green and red lines represent the transfer learning accuracy and validation curves, (**b**) epoch vs. loss, (**c**) the confusion matrix from the pre-trained model, and (**d**) the confusion matrix with the re-trained model.

**Figure 7.** Comparison between the proposed AI method and the manual method for counting. (**a**) A whole frame LSIT micrograph, (**b**) magnified region-of-interest (ROI) of the whole frame LSIT micrograph, (**c**) automatically classified cell diffraction patterns of the ROI, and (**d**) comparison of the cell counts between those two different modalities.

#### **4. Conclusions**

In conclusion, we have explored the advantages of using neural networks in the characterization of LSIT micrographs. Here, we have perfected neural networks that can automatically improve the signal quality and classify the cell types. We find that this neural network can classify the RBC and WBC with great accuracy (i.e., over 98%), and the cancer cell with an accuracy of about 90%. This network is also flexible for adapting to newer cell lines by retraining the trained network with very few samples. Together with this algorithm, the lightweight and cost-effective LSIT setup can be utilized as a point of care system for the diagnosis of pathological disorders in the resource-limited setup of our world. In our future work, we aim to combine the denoising and classification modalities due to the significant overlap in their operation. This will remove the dual training times as well as minimize computation costs. Also, we aim to work on improving their performance and deploying it in real scenarios.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/bios12030144/s1, 1. Convolutional Neural Network workflow; 2. Gaussian Noise; 3. Traditional denoising methods; 4. Comparison of the denoised outputs from various modalities; 5. SNR of individual samples of various cell types; 6. Grad-CAM and Saliency maps of the individual cell-lines.

**Author Contributions:** Conceptualization, S.S. (Sungkyu Seo) and M.R.; methodology, R.V., S.S. (Sanghoon Shin), V.C., K.K. and S.S. (Smit Sanghavi); validation, R.V., V.C., K.K., S.S. (Smit Sanghavi) and M.R.; formal analysis, R.V., K.K., S.S. (Sungkyu Seo) and M.R.; investigation, R.V., V.C., K.K. and S.S. (Smit Sanghavi); resources, S.S. (Sungkyu Seo) and M.R.; data curation, R.V., S.S. (Sanghoon Shin), V.C., K.K. and S.S. (Smit Sanghavi); writing—original draft preparation, R.V., V.C., K.K., S.S. (Smit Sanghavi) and M.R.; writing—review and editing, R.V., S.S. (Sanghoon Shin), S.S. (Sungkyu Seo) and M.R.; visualization, R.V., S.S. (Sanghoon Shin), V.C., S.S. (Sungkyu Seo) and M.R.; supervision, M.R.; project administration, S.S. (Sungkyu Seo) and M.R.; funding acquisition, S.S. (Sungkyu Seo) and M.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** M.R. acknowledges the seed grant No. ORSP/R&D/PDPU/2019/MR/RO051 of PDPU (the AI software development part) and the core research grant No. CRG/2020/000869 of the Science and Engineering Research Board (SERB), India. S.S. acknowledges the support of the Basic Science Research Program (Grant#: 2014R1A6A1030732, Grant#: 2020R1A2C1012109) through the National Research Foundation (NRF) of Korea, the Korea Medical Device Development Fund (Project#: 202012E04) granted by the Korean government (the Ministry of Science and ICT, the Ministry of Trade, Industry and Energy, the Ministry of Health & Welfare, the Ministry of Food and Drug Safety), and the project titled 'Development of Management Technology for HNS Accident' and 'Development of Technology for Impact Assessment and Management of HNS discharged from Marine Industrial Facilities', funded by the Ministry of Oceans and Fisheries of Korea.

**Institutional Review Board Statement:** All human blood samples were approved by the Institutional Review Board of Anam Hospital, Korea University (# 2021AN0040).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

## *Article* **Decoding Vagus-Nerve Activity with Carbon Nanotube Sensors in Freely Moving Rodents**

**Joseph T. Marmerstein, Grant A. McCallum and Dominique M. Durand \***

Neural Engineering Center, Biomedical Engineering, Case Western Reserve University, Cleveland, OH 44106, USA<EMAIL_ADDRESS>(J.T.M<EMAIL_ADDRESS>(G.A.M.) **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The vagus nerve is the largest autonomic nerve and a major target of stimulation therapies for a wide variety of chronic diseases. However, chronic recording from the vagus nerve has been limited, leading to significant gaps in our understanding of vagus nerve function and therapeutic mechanisms. In this study, we use a carbon nanotube yarn (CNTY) biosensor to chronically record from the vagus nerves of freely moving rats for over 40 continuous hours. Vagal activity was analyzed using a variety of techniques, such as spike sorting, spike-firing rates, and interspike intervals. Many spike-cluster-firing rates were found to correlate with food intake, and the neural-firing rates were used to classify eating and other behaviors. To our knowledge, this is the first chronic recording and decoding of activity in the vagus nerve of freely moving animals enabled by the axon-like properties of the CNTY biosensor in both size and flexibility and provides an important step forward in our ability to understand spontaneous vagus-nerve function.

**Keywords:** vagus nerve; intraneural; decoding; intrafascicular; recording; carbon nanotube

#### **1. Introduction**

The vagus nerve innervates nearly every internal organ, providing sensory input to the brain and parasympathetic-control inputs to the viscera. Therefore, abnormal vagus-nerve activity has been linked to many chronic diseases, such as epilepsy, diabetes, hypertension, and cancer [1–5]. Vagus-nerve stimulation has been used to treat a wide variety of diseases [6], most successfully implemented for the treatment of epilepsy [7], even while the mechanisms are not well understood and direct recordings of vagal activity associated with disease are not available [8]. The majority of vagal afferent fibers come from the gut [9,10], and abnormal vagal activity has been clearly implicated in eating and metabolic disorders [11–15]. In this study, we analyze the first chronic recordings of vagal spikes and the correlation of signals to several behaviors in healthy rats.

The chronic recording of vagal signals has been limited, partially due to the difficulty in chronically recording high-quality signals in small autonomic nerves. Extraneural cuff electrodes have proven to be very effective peripheral nerve interfaces, allowing for selective stimulation [16] and some selectivity in recording [17–20]. However, the insulating perineurium layer between the active nerve fibers and the recording electrodes results in a low signal-to-noise ratio (SNR) or requires desheathing of the nerve. Thus, intraneural or intrafascicular electrodes placed much closer to active fibers may be necessary for certain recording applications, providing a higher SNR and higher selectivity for multifascicular nerves. However, intraneural electrodes are more invasive and thus, have issues related to long-term stability [21–25]. In particular, small autonomic nerves necessitate smaller and more flexible neural sensors. We have previously shown that carbon nanotube yarn (CNTY) electrodes have favorable properties for nerve interfacing: specifically, their small size, high flexibility, and low impedance. Thus, they provide a stable, high-SNR interface for chronic recording in small autonomic nerves in rats, with high-quality signals continuing up to four months after implantation [26]. Furthermore, we have developed techniques

**Citation:** Marmerstein, J.T.; McCallum, G.A.; Durand, D.M. Decoding Vagus-Nerve Activity with Carbon Nanotube Sensors in Freely Moving Rodents. *Biosensors* **2022**, *12*, 114. https://doi.org/10.3390/ bios12020114

Received: 8 January 2022 Accepted: 5 February 2022 Published: 11 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

for recording activity in the cervical vagus nerves of rats without anesthesia, allowing for the first chronic recordings of truly spontaneous vagal activity. This technology has been successfully applied to make the first direct measurements of vagal tone in freely moving animals [27].

Due to the prevalence of gastric afferents in the vagus nerve, we expect signaling from the gut to be the dominant activity present in vagal recordings. Previous studies have shown that vagal afferents are sensitive to mechanical stimulation of the gut [28] and to gastric hormones which regulate food intake and gastric motility [29–31]. There are also reflex pathways which modulate efferent vagal activity in response to gastric distention and contractions [32]. However, such signals have not been reported from chronic, unanesthetized animals, and this CNTY-electrode biosensor demonstrates the ability to decode vagal activity related to various animal behaviors, such as eating.

In this study, we continuously record spontaneous vagal-spiking activity from awake, freely moving rats for >48 h up to two weeks after implantation. To our knowledge, this is the first time this has been successfully demonstrated. The neural-recording data was synchronized with continuous video recording of the subjects. Spike sorting is used to separate semi-distinct spike clusters, which are then correlated to animal behavior identified from the video recordings. Interspike interval distributions are also found to change in response to food intake, presenting another neural feature that can be used to decode spontaneous vagal activity. We report several spike clusters that show tuning to animal eating, and the firing dynamics of multiple decoded spike clusters can be used to classify eating compared to drinking, grooming, and resting behaviors.

#### **2. Materials and Methods**

#### *2.1. CNTY Electrode Manufacture*

CNT yarns were manufactured at Case Western Reserve University, as described previously [26]. CNTYs were then connected to 35NLT®-DFT® wire (Fort Wayne Metals, Fort Wayne, IN, USA) with silver conductive epoxy (H20E, EPO-TEK), creating a CNTY-DFT® junction. Dacron mesh and silicone elastomer (MED-4211/MED-4011, NuSil Silicone Technology, Carpinteria, CA, USA) were added to seal the junction, confirmed by measuring the impedance of the junction at 1kHz in a saline bath. The free end of the CNTY was tied to the end of an 11-0 nylon suture (S&T 5V33) using a fisherman's knot, as shown in Figure 1A. The entire CNTY was coated with parylene-C (5 μm thick vapor deposition coating, SMART Microsystems, Elyria, OH, USA) on a custom rack which masks the suture needle from coating. Then, a small section (~200 μm long) of parylene-C was removed approximately 500 μm behind the CNTY-suture knot using a laser spot welder (KelanC Laser, set to 1A current, 0.3 ms pulse width, and 300 μm diameter), as shown in Figure 1B. Figure 1C shows the CNTY-suture knot outside of the nerve after implantation. Electrode viability was confirmed by measuring the impedance of the recording site before and after using the laser.

**Figure 1.** Electrode implantation, histology, and recording methods. (**A**) Diagram of CNTY electrode mated with an 11-0 nylon suture with a fisherman's knot. (**B**) Section of CNTY electrode deinsulated by laser. (**C**) Vagus nerve with two implanted CNTY electrodes. CNTY-suture knots are shown with arrows. (**D**/**E**) Diagram showing the setup for continuous recording of vagal activity and video for behavior identification. Signals travel from the implants to the headcap connector mounted on the animal's skull, where they are digitized and amplified by the custom amplifier board shown. These signals are then routed through a commutator, which can rotate and allows the animal to move freely without twisting or pulling on the cable. From the commutator, the signals are sent to an Intan USB interface board, which is powered by an external DC-power source and finally sends the signals to a computer, where they are saved and can be viewed in real time. A video camera is manually synced to the vagal recordings. (**F**) Fluorescent images showing collagen + cellular encapsulation of CNTY electrodes implanted in the vagus nerve for seven days. (**G**) Toluidine blue-stained nerve section showing encapsulation of a CNTY electrode implanted for two weeks.

#### *2.2. Surgery*

All surgical and experimental procedures were done with the approval and oversight of the Case Western Reserve University Institutional Animal Care and Use Committee to ensure compliance with all federal, state, and local animal welfare laws and regulations. Electrodes were implanted in male Sprague Dawley rats between 7–12 weeks of age.

To expose the left cervical vagus nerve, a midline incision was made along the neck. The muscles and salivary glands were separated and held in place, revealing the carotid sheath which contains the carotid artery and vagus nerve. The vagus nerve was carefully separated from the carotid artery using blunt dissection and held in slight tension using a glass hook. CNTY electrodes were implanted by sewing the suture through the nerve for ~2 mm, then pulling the suture until the CNTY-suture knot was pulled through. Then, the electrode was pulled back so that the knot sat against the epineurium, ensuring the recording site remained inside the nerve, as shown in Figure 1C. Two electrodes were implanted with ~2 mm separation; the extra suture and needles were cut off after implantation, and the nerve, electrodes, and junctions were covered with ~1 mL of fibrin glue (Tisseel, Baxter International Inc., Deerfield, IL, USA) to help secure the area for recovery. Next, the DFT wires were tunneled from the neck to the back of the skull and soldered to a 5-pin Omnectics connector (Omnetics Connector Corporation MCP-5-SS). The skin on top of the skull was opened, and the connector was fixed on top of the skull with dental cement. The amplifier ground was connected to a screw placed in the skull, which also helps keep the headcap in place. Electrodes were implanted for chronic recording in two animals, and animals were given one week for recovery before recording.

#### *2.3. Recording*

Recordings were carried out continuously in awake, behaving animals for 56 and 40 h (Rat 1 and 2, respectively). A custom-built PCB with an Intan RHD2216 recording chip was attached to the headcap connector, which was secured to the animal with a 3D-printed locking mechanism and attached to a PlasticsOne® (Roanoke, VA, USA) commutator, allowing the rat to move around the cage without tangling or pulling on the connector cable [27]. Input signals were routed to eight amplifier channels, using 8-channel hardware averaging to decrease amplifier noise. Output from the amplifier board was run through the commutator into an Intan RHD USB Interface board (Intan part #C3100), which is powered by an external battery supplying 5V DC power. Signals are then routed to a computer where they are saved for offline analysis and can be viewed in real time.

Neural recordings were sampled at 20 kHz with a 5 kHz low-pass filter. Recordings were started around 10 AM (approximately four hours after the start of the light cycle). During ENG recording, a video camera was used for simultaneous video recording. The camera was equipped with an infrared light and infrared sensor, allowing for filming even during the dark cycle. The camera was connected to the recording computer and manually synced to the recording. A diagram of the recording setup can be seen in Figure 1D,E.

#### *2.4. Signal Processing*

ENG data were imported into MATLAB, where they were further processed. ENG was band-pass filtered from 500–5000 Hz to minimize interference from EMG, ECG, or other possible sources. The filter bandwidth was kept relatively wide to minimize distortion of spike waveforms. Spikes were detected and sorted into clusters using the UltraMega-Sort2000 software in MATLAB, using a threshold of eight times the RMS of baseline. Spike waveforms (3 ms long) were transformed into the principal component space, and principal components accounting for 95% of the total waveform variance were used for spike clustering. Spike clustering was done using k-means clustering of spike waveform principal components, with a maximum of k = 256 clusters. Using the UMS2000 software, clusters were further analyzed for better separation and exclusion of artifacts. First, outliers were removed if they had a z-score greater than 500 on the χ2 distribution of distance to the cluster center. Clusters were removed from analysis if the spike waveform contained a second, larger threshold crossing (i.e., removing of spikes which were detected twice due to threshold-crossing of the spike tail). Clusters were also removed if spike width was less than 0.3 ms and amplitude was greater than 1mV (presumed recording artifacts) or if spike width was greater than 2 ms. Spike waveform values were used to calculate

spike amplitudes (difference between the maximum and minimum voltage values) and the spike RMS. Spike-cluster-firing timings were also used to calculate cluster-firing rates and interspike intervals (ISI). Average spike amplitudes over time are shown in Supplementary Figure S1, and spike RMS was used to calculate average SNR, shown in Supplementary Figure S2. Animal behaviors (eating, drinking, grooming, and resting) were identified via video recording. The overall data processing and analysis workflow is diagrammed in Figure 2.

**Figure 2.** Diagram of data processing and analysis workflow. Vagal ENG and video are recorded simultaneously from freely moving rats. Spike sorting is used to decode spike metrics, which are analyzed with respect to animal behaviors identified from the video.

#### *2.5. Histology*

Toluidine blue staining: the image shown in Figure 1G was obtained from an implanted nerve which was fixed, sectioned, and stained with toluidine blue. Two weeks after implantation, animals were perfused with 1.25% glutaraldehyde, 1% formalin, and 0.1 M phosphate buffer. This fixative solution is approximately 640 mOsM/kg. Animals were injected with 0.2–0.5 mL of 1% procaine at 37 ◦C through the left ventricle. Followed by 200 mL of the fixative solution perfused at 37 ◦C using a variable speed peristaltic pump. After completing the perfusion process, the vagus nerve was dissected at the implant location. The complete nerve section was transferred into a postfixative solution (1% osmium tetroxide in 100-mM phosphate buffer) for two hours at room temperature before being transferred to 4 ◦C. Following postfixation, the nerve tissue was dissected in 1-mm-long pieces and embedded in an epoxy resin. Sections (0.7 m) were cut from the epoxy blocks using a diamond knife (DiATOM) microtome. Toluidine blue (1% toluidine blue and 2% borate) was used to stain the nerve axons.

Fluorescent staining: the image shown in Figure 1F was obtained from an implanted nerve which was optically cleared using the CLARITY protocol [33,34]. Seven days after implantation, the vagus nerve was extracted and immediately placed into hydrogel monomer solution. The sample was passively cleared and stained with a collagen antibody, as described by our group previously [26]. DAPI staining was done by placing the sample in VectaShield with DAPI (Vector Laboratories) on a glass-bottom petri dish (Ted Pella, Inc., Redding, CA, USA). Samples were imaged on a Leica SP8 gSTED Super-Resolution Confocal microscope (Leica Microsystems, Wetzlar, Germany).

#### *2.6. Statistical Methods*

Where relevant, results are reported as mean ± standard deviation. Average spike waveforms in Figure 2 are shown with shaded areas representing the 95% confidence interval. Overall spike-firing rate, median spike amplitudes, and average spike SNR over time were fitted with a linear regression to determine if the slope was different from zero, with slopes and *p*-values shown in Figure 2 and Supplemental Figures S1 and S2. Spike clusters were grouped based on their response to eating, and firing rate changes before, during, or after eating for each group were compared to baseline group firing rates using a one-sample *t*-test, with a significance level of 0.01 and a Bonferroni correction (α = 5.6 × 10−4), as shown in Table 1. ISI distributions of the before, during, and after

eating periods were compared to noneating periods using a two-sample Kolmogorov– Smirnov test, with a significance level of 0.01 and a Bonferroni correction for the number of tested distributions (α = 2.2 × 10−5, Supplementary Table S1). All tests performed were two tailed.

**Table 1.** Firing rates of cluster groups relative to eating. Sorted clusters are separated into five cluster groups based on their response to eating. Table shows the number of clusters of each group recorded in both animals and the behavior of those cluster groups before, during, and after eating: up arrow/green color means an increased firing rate, dash/yellow color means no change in firing rate, and down arrow/red color means a decreased firing rate for the cluster group.

#### **3. Results**

#### *3.1. CNTY Electrodes Record Stable Spikes from Freely Moving Animals*

We have previously shown that CNTY electrodes can record spikes from the glossopharyngeal and vagus nerves in anesthetized rats and can be used to measure vagal tone in freely moving animals [26,27]. Here, we demonstrate a novel continuous chronic-recording setup (shown in Figure 1D,E) to record unanesthetized spiking activity which can be sorted into semi-distinct clusters. A total of four electrodes were implanted, two each in the left cervical vagus nerves of two rats, with an average impedance of 11.7 ± 6.5 kΩ at the time of implantation (measured at 1 kHz). Further measurements of CNTY electrode impedances for long-term implants have been published previously [26,27]. Figure 3A shows an example of filtered ENG with several recorded spikes, and Figure 3B–E show several example spike clusters from two animals. A total of 132 spike clusters were identified (56 in Rat 1, and 76 in Rat 2). Clusters are referred to as RatNumber.ClusterNumber (e.g., Cluster 1.21 is Cluster 21 from Rat 1). Average peak-to-peak amplitude of recorded spikes was 152 ± 97 μV for Rat 1 and 180 ± 162 μV for Rat 2. Spike SNR, defined as the average RMS of the spike waveforms compared to the RMS of the baseline, was 7.0 ± 4.9 for Rat 1, and 9.1 ± 5.3 for Rat 2. This is significantly larger than published SNR for acute recording with either the TIME or the LIFE electrodes [33,34]. Furthermore, median spike amplitude for all recorded spikes was stable over the recording time for Rat 2 and slightly increased over time for Rat 1, as shown in Supplementary Figure S1. Overall spike-firing rates were also consistent over the recording periods for both animals: Figure 3F,G show the average firing rates for each hour of recording, with least-squares regression lines showing no significant change in firing rate over time. Similarly, average spike SNR was stable over the recording time for both animals, as shown in Supplemental Figure S2. Thus, we are able to continuously record vagal spikes which have stable amplitude, SNR, and firing rates over time.

**Figure 3.** Spontaneous spikes recorded in freely moving animals. (**A**) Filtered ENG showing example recording spikes. (**B**–**E**) Example clusters sorted from recorded spikes in two animals. (**F**,**G**) Spike firing rate over recording time for two animals. Neither animal had a significant change in firing rate over time.

#### *3.2. Spike Clusters' Activity Is Correlated with Eating*

Identifying the function of spontaneous spikes in freely moving animals is important to understanding how vagal fibers modulate their activity during normal animal behavior. Given the high ratio of gastric afferents in the vagus, most vagal spiking is involved with gastric signaling.

After animal-eating times were identified from video recordings, they were compared to the firing rates of individual spike clusters. In both animals, several clusters show a significant increase in firing rate that occurs <25 min before eating. Some clusters also had increased or decreased firing that occurred during eating, while others had increased firing that occurred <10 min after eating. Figure 4A,B show raster plots for one such spike cluster from each animal, with each row representing one eating event (shown by the shaded grey area). Figure 4C,D show the average firing rate of these clusters relative to the eating events, along with the overall average firing rate for each cluster. Cluster 1.36 (Figure 4A,C) had higher-than-average activity in the 25 min before eating, and higher-than-average activity in the 10 min following eating, with no change occurring during food consumption. Similarly, the firing rate of Cluster 2.1 is increased before and during eating, and unchanged after eating.

**Figure 4.** Example spiking activity related to eating. (**A**) Raster plot for Cluster 1.36. Grey-shaded boxes represent eating events, with dots representing spikes. (**B**) Raster plot for Cluster 2.1. (**C**) Firing rate of Cluster 1.36 relative to eating, averaged for all eating events. Red line represents the overall average firing rate of Cluster 1.36. (**D**) Firing rate of Cluster 2.1 relative to eating.

Many clusters exhibited a mix of behaviors, showing firing rates before, during, or after eating that were significantly different from baseline activity (*p* < 0.01 with Bonferroni correction). To analyze cluster behavior related to eating, clusters were sorted into groups based on their firing rate response before eating (from 25 min before, until the start of eating), during eating, and after eating (end of the eating event, until 10 min after eating). These data are summarized in Table 1 for both rats, which show how the cluster-firing rates changed for each group and the number of clusters from each animal which make up each group. The table shows the direction of change and associated p-value for the changes in firing rate of each group in the different eating-related periods (sum of the spiking activity in all clusters within a group compared to the baseline firing activity for the clusters in that group). Only 3 of the 132 recorded clusters did not showing any significant tuning to eating behavior. While specific spiking correlations are unique to each subject, they are consistent within each animal, and Figure 3 and Table 1 show that we can identify spike clusters that exhibit firing rate changes before, during, and after eating in both subjects.

#### *3.3. Spike Cluster Interspike Intervals Show Changes in Bursting Related to Eating*

Spikes are often observed exhibiting bursting behavior, where fibers tend to fire at specific frequencies. Bursting behavior can be seen in Figure 4A,B, where spikes appear in clumps. To quantify bursting, spike cluster interspike intervals (ISIs) were calculated for noneating, pre-eating, during eating, and posteating time windows. Eating-related distributions were compared to noneating distributions using a two-sample Kolmogorov– Smirnov test and were plotted in a histogram. Figure 5 shows ISI distributions for noneating, pre-eating, during eating, and posteating periods for one example cluster (Cluster 1.8, which is part of Group II and has increased activity before and after eating). In Figure 5A, we can see that the peak ISI of this cluster during noneating times is around 21 ms or a 48 Hz firing rate. However, in the 25 min before eating, this distribution shifts to the left, peaking instead at 7 ms or 143 Hz, signifying an increase in the bursting firing rate before eating. In the 10 min following eating, the bursting rate returns to the noneating value, though the ISI peak is more pronounced, meaning that bursting is a more prevalent spike behavior after eating. After eating, we also observe a secondary peak around 47 ms (21 Hz). During eating, the ISI distribution is not significantly different from noneating; thus, the bursting activity of Cluster 1.8 is changed before and after, but not during, eating behavior. In total,

10 clusters in Rat 1 and 18 clusters in Rat 2 demonstrated changes in ISI distribution related to eating.

**Figure 5.** Interspike interval histograms for Cluster 1.8. (**A**) ISI histogram for noneating periods, which has a peak around 21 ms. (**B**) ISI histogram for pre-eating periods, which has a peak around 7ms, and a significantly different ISI distribution compared to noneating periods. (**C**) ISI histogram for eating periods, which has a peak around 23 ms and is not significantly different from noneating periods. (**D**) ISI histogram for post-eating periods, which has a peak around 21 ms and a secondary peak around 47 ms, and a significantly different ISI distribution compared to noneating periods.

These data are summarized in Supplemental Table S1, which shows *p*-values comparing noneating and eating-related ISI distributions for any cluster which showed a significant change. The 18 clusters in Rat 2 only showed a change in ISI distribution during eating, with no changes either before or after. The 10 clusters in Rat 1 each showed changes before eating, while some also had a significantly different ISI distribution during or after eating as well. Figure 5 and Supplementary Table S1 show that some of the spike clusters which are tuned to eating are observed to change bursting activity related to eating, though not all the clusters which show changes in overall activity have altered ISI/bursting behavior.

#### *3.4. Spike-Cluster-Firing Rates Can Be Used to Classify Eating Compared to Other Behavior*

In addition to showing that individual spike clusters are correlated with food intake, we also examined whether spike-firing rates are sufficient to classify the times during which the animal is eating, compared to other behaviors, such as drinking, grooming, and resting. A multinomial logistic regression model was constructed, with behaviors and spike-cluster-firing rates averaged over 30 s. The model uses firing rates from each of the recorded clusters, as well as firing rates during peak delayed or preceding correlations with eating. The models were trained on the first 2/3 of recording data and tested on the final 1/3 of recording data. Figure 6 shows the confusion matrices for both animals, which show the performance of the model for classifying behavior with a probability threshold of π = 0.5 for classification. Percentages on the y-axis show the amount of time spent doing each behavior as a percentage of total recording time. In Rat 1, the model was able to classify eating most accurately, with a 73.1% accuracy. In Rat 2, the model performed best at classifying resting, with a 93.8% accuracy. Additionally, we can see by plotting the receiver operating characteristic (ROC) curves and the associated areas under the curve (AUC) in Supplementary Figure S3 that both models performed better than random chance for

almost all behaviors (the only exception being classifying other activity in Rat 1). Overall, these results show that the firing rates of spontaneous vagal spikes sorted into clusters are sufficient to classify eating behavior in freely moving animals.

**Figure 6.** Confusion matrix for classifying animal behavior based on spike firing rates. Blue-colored cells show rates of correct classification, and orange-colored cells show rates of incorrect classification, such that each row sums to 100%. *Y*-axis labels show the percentage of recording time spent doing each behavior. (**A**) Confusion matrix for the classification of behavior in Rat 1. (**B**) Confusion matrix for the classification of behavior in Rat 2.

#### **4.  Flexural rigidity, measured with an atomic-force microscope, shows that the CNTY electrodes are >10 times more flexible than PtIr electrodes of the same diameter (3.3 ± 1.5 × 10−12 N× m2 for the CNTY compared to 2.0 ± 0.57 × 10−10 N× m2 for the PtIr) [35]. Partly due to its small size (10 μm diameter) and flexibility, this axon-like biosensor has demonstrated stable, low impedance with chronic implantation, stable high SNR, and minimal evidence of chronic inflammation or nerve damage [26]. Furthermore, we have shown that CNTYs can be used for chronic recording in small autonomic nerves, such as the vagus and glossopharyngeal nerves, and for stimulation in larger somatic nerves and fascicles, such as the rat sciatic nerve [26]. Compared to previous intrafascicular interfaces, CNTYs provide higher SNR and improved stability and should be further investigated as a component of neural sensor devices.

Vagus-nerve-stimulation therapy is a rapidly growing field, with a wide variety of companies and studies investigating its use for treatment of a wide variety of diseases, including epilepsy, obesity, and heart failure [7,36–38]. However, many VNS studies have reported ambiguous results, pointing to the need for an improved understanding of natural vagal function, VNS mechanisms, and closed loop control of stimulation. Neural interfaces which allow for stable, high-SNR recordings are necessary for high-fidelity closed-loop control, and chronic recording in animal models may be used to better understand vagal function and response to therapy. In this study, we utilize the CNTY neural interface to show that eating-related spikes can be decoded from continuous chronic recordings in the vagus nerve, providing the first demonstration that spontaneous, physiologically specific signals can be recorded. We also show that spike-firing rate and interspike interval distributions show differing responses to physiological changes, which may be used as neural features for long-term recording and closed-loop systems, though further histological analysis may also be necessary to show the impact of chronic implantation on nerve health.

The vagus nerve contains afferent and efferent fibers that sense and control nearly every internal organ, playing a vital role in homeostasis, reflex pathways, and responses to physiological changes. While individual vagal-spiking activity has been recorded from isolated fibers and from acute intraneural recordings, to our knowledge, this is the first time spikes have been recorded in the vagus nerve in a chronic model. Combined with previous studies on recording average vagal RMS [27], we have demonstrated that various spontaneous and physiologically relevant signals can be recorded from the rat vagus nerve from freely moving animals using CNTY electrodes. Spike-firing rates stayed consistent for up to 56 recording hours (7–10 days after implantation), suggesting that the electrode interface is relatively stable during that time. Though spike clusters detected in the peripheral nervous system likely represent multiunit activity, Figure 2 shows that CNTY-recorded spikes can be sorted into clusters, allowing for more specialized decoding of relevant signals. Individual spike clusters recorded in freely moving rats show changes in firing rate before, during, and after eating; clusters were further sorted into different groups based on their firing behavior, as shown in Figure 3 and Table 1. Almost all recorded spike clusters showed increased activity up to 25 min before the start of eating. Firing rate also increased during and after eating for several cluster groups. Additionally, Figure 4 and Supplemental Table S1 show that some spike clusters exhibit changes in their ISI distribution during different eating phases. ISI distribution and peak ISI values are an important metric for describing spike-bursting behavior, and bursting frequencies may change independently of overall spiking activity in the cluster (due to clusters containing recordings from multiple individual axons). Thus, analysis of overall cluster activity and cluster ISI distributions are important metrics of vagal activity. Finally, Figure 5 and Supplementary Figure S3 shows that spike-cluster-firing rates can be used to accurately classify eating, drinking, grooming, and resting behavior, an important proof of concept for a closed-loop VNS system. Future studies could analyze how overall vagal-spiking activity and bursting rates respond to models of chronic disease, physiological stimuli (such as fasting), or to VNS treatment.

There are several important questions regarding the behavior of these fibers and how they might respond to physiological changes. One possibility is that these fibers may be related to the secretion or sensing of ghrelin and cholecystokinin (CCK) in the gut, peptide hormones that are known to regulate hunger and satiety via the vagus nerve and can even be found in small concentrations in the brain [30,39]. Ghrelin increases food intake and weight gain in rats and suppresses vagal activity of some gastric afferent fibers [30]. CCK, on the other hand, suppresses appetite while stimulating gastric afferent discharge [40,41]. However, the full picture relating the secretion of these hormones and their effects to vagus-nerve activity is not known. Future studies could investigate how administration of these peptides alters vagal activity in freely moving animals or how changes in vagalspiking behavior correlate with changes in the concentration of gastric hormones in the gastrointestinal tract, in the blood, or in the CNS. Furthermore, changes in diet, such as high-fat or high-carbohydrate models of diet-induced obesity, have been shown to alter vagal satiety signaling and may have an effect on hunger signaling as well [14,42]. Utilizing the CNTY chronic-recording model described here would allow for insights into the effects of diet on vagal gastric signaling, and a more detailed measurement of food intake would allow for analysis of how differences in vagal spiking relate to differences in dietary behavior.

One important application for chronic recording of vagal activity is measuring the acute and chronic response of vagal signaling to various therapy approaches. High frequency VNS (vBLoc®) is thought to suppress afferent hunger signaling to reduce overeating in obese patients [43]. vBLoc® stimulation is performed with a fixed on/off cycle; when stimulation is on, vagal activity is suppressed, which would decrease hunger signaling (e.g., the activity observed in cluster groups 1–5), and satiety signaling (e.g., cluster groups 1–3). However, studies have shown that vBLoc® reduces hunger and *increases* satiety, though it is not clear how significant this effect is compared to a placebo [44]. More thorough investigation of the effects of vBLoc® on hunger, satiety, and vagal activity are necessary to better understand therapy performance. The abiliti® gastric stimulation system, on the other hand, is closed-loop, producing satiety during preprogrammed periods and in response to eating, thus reducing food intake and helping create a more stable meal schedule [45]. Patients with this system report improved self-control while eating, decreased binge eating, and reduced sensitivity to hunger [46], likely as a result of increased afferent activity during and after eating. Combining animal models of these stimulation systems with chronic spike recording using CNTY electrodes would allow for direct investigation of the mechanisms of these therapies [43]. Furthermore, chronic recording could be paired with other VNS paradigms, such as VNS for epilepsy or depression, to investigate possible off-target side effects.

Overall, our results show that it is possible to record chronic signals in the rat vagus nerve continuously, opening the door for studies that were previously not possible. Furthermore, the ability to detect and decode spontaneous spiking activity from chronic vagal recordings could allow for a more detailed analysis of vagus-nerve response to changes in diet, therapy, and behavior. This technology could be used to develop closed-loop VNS for metabolic disorders, which adapt stimulation based on recorded vagal activity. Chronic recording in animal models can also be used to further study the vagal pathways that control food intake and how they respond to VNS and other treatments.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/bios12020114/s1, Supplementary Figure S1. Median spike amplitude over time. (A). Median spike amplitude for spikes recorded in Rat 1 slightly increase over the recording time, with a statistically significant slope of 0.39μV per hour. (B). Median spike amplitude of spikes recorded in Rat 2 did not change over the recording time. Supplementary Figure S2. Average spike SNR (spike RMS divided by noise RMS) over time. (A). SNR of spikes recorded in Rat 1 did not change over the recording time. (B). SNR of spikes recorded in Rat 2 did not change over the recording time. Supplementary Table S1. Differences in ISI distributions for before, during, and after eating periods, compared to non-eating periods, for all clusters which had at least one group with a significant change. Cluster groups are shown for each cluster (see Table 1), and non-significant p-values are not shown. Supplementary Figure S3. Receiver operating characteristic (ROC) curves and area-under-the-curve (AUC) values to assess performance of a multinomial logistic regression model to classify animal behaviors based on spike cluster firing rates. Dotted lines show the expected ROC curve for a random classifier. (A). ROC curve for classifying drinking in Rat 1, with AUC = 0.86. (B). ROC curve for classifying drinking in Rat 2, with AUC = 0.62. (C). ROC curve for classifying eating in Rat 1, with AUC = 0.94. (D). ROC curve for classifying eating in Rat 2, with AUC = 0.82. (E). ROC curve for classifying grooming in Rat 1, with AUC = 0.88. (F). ROC curve for classifying grooming in Rat 2, with AUC = 0.90. (G). ROC curve for classifying resting in Rat 1, with AUC = 0.82. (H). ROC curve for classifying resting in Rat 2, with AUC = 0.86. I: ROC curve for classifying other activity in Rat 1, with AUC = 0.47. (J). ROC curve for classifying other activity in Rat 2.

**Author Contributions:** Conceptualization, J.T.M., G.A.M. and D.M.D.; methodology, J.T.M., G.A.M. and D.M.D.; software, J.T.M.; validation, J.T.M. and G.A.M.; formal analysis, J.T.M.; investigation, J.T.M.; resources, D.M.D.; data curation, J.T.M. and G.A.M.; writing—original draft preparation, J.T.M.; writing—review and editing, J.T.M., G.A.M. and D.M.D.; visualization, J.T.M.; supervision, G.A.M. and D.M.D.; project administration, D.M.D.; funding acquisition, D.M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Congressionally Directed Medical Research Programs (Autonomic nervous system activity and the implications on breast cancer metastasis, W81XWH-18-1- 0581), and the National Institutes of Health, (Nerve reshaping for improved selectivity, 5R01NS032845- 22). Confocal microscopy was performed at the Light Microscopy Imaging Facility at Case Western Reserve University, made available through the Office of Research Infrastructure (NIH-ORIP) Shared Instrumentation Grant (S10OD016164).

**Institutional Review Board Statement:** All surgical and experimental procedures were done with the approval and oversight of the Case Western Reserve University Institutional Animal Care and Use Committee to ensure compliance with all federal, state, and local animal welfare laws and regulations.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request. Some custom-analysis code was created to analyze spiking activity and is available from the corresponding author upon reasonable request.

**Acknowledgments:** Assistance for tissue fixation, sectioning, staining, and imaging for Figure 1G was provided by Hisashi Fujioka of the Case Western Reserve University Cryo-Electron Microscopy Core. Special thanks to William Marcus, who helped monitor animals during surgery and recovery and assisted with recording.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Detection of Liver Dysfunction Using a Wearable Electronic Nose System Based on Semiconductor Metal Oxide Sensors**

**Andreas Voss 1,2,\*, Rico Schroeder 1,3, Steffen Schulz 1, Jens Haueisen 2, Stefanie Vogler 4, Paul Horn 4, Andreas Stallmach 4 and Philipp Reuken 4**

**Abstract:** The purpose of this exploratory study was to determine whether liver dysfunction can be generally classified using a wearable electronic nose based on semiconductor metal oxide (MOx) gas sensors, and whether the extent of this dysfunction can be quantified. MOx gas sensors are attractive because of their simplicity, high sensitivity, low cost, and stability. A total of 30 participants were enrolled, 10 of them being healthy controls, 10 with compensated cirrhosis, and 10 with decompensated cirrhosis. We used three sensor modules with a total of nine different MOx layers to detect reducible, easily oxidizable, and highly oxidizable gases. The complex data analysis in the time and non-linear dynamics domains is based on the extraction of 10 features from the sensor time series of the extracted breathing gas measurement cycles. The sensitivity, specificity, and accuracy for distinguishing compensated and decompensated cirrhosis patients from healthy controls was 1.00. Patients with compensated and decompensated cirrhosis could be separated with a sensitivity of 0.90 (correctly classified decompensated cirrhosis), a specificity of 1.00 (correctly classified compensated cirrhosis), and an accuracy of 0.95. Our wearable, non-invasive system provides a promising tool to detect liver dysfunctions on a functional basis. Therefore, it could provide valuable support in preoperative examinations or for initial diagnosis by the general practitioner, as it provides noninvasive, rapid, and cost-effective analysis results.

**Keywords:** electronic nose; liver dysfunction; cirrhosis; semiconductor metal oxide gas sensor

#### **1. Introduction**

Metabolic disorders are sometimes connected with typical odors which can be measured on breath, sweat, or other excreta from humans. Examples are ammonia odor, which is related to renal diseases, and acetone odor, which is related to diabetes.

The beginnings of the use of electronic noses (e-noses) date back to pioneering work by a few research groups, such as Hartman, Wilkens, Dodd, and Moncrieff [1–4]. Here, the foundation was laid for specific odors to be detectable and, thus, evaluable with suitable electronics and analysis technology. The concept of sampling breath for health monitoring was initially conceived in the 20th century. In 1952, Henderson [5] reported on the increased acetone content of breath samples from young diabetics, promoting an interest in the content of breath [6].

In recent decades, improvements in materials, sensors, electronics, and signal processing technologies have led to a rapid increase in the development and application of e-noses [7–9]. E-noses are used, among other things, to analyze, detect, discriminate, classify, and monitor gas components or odors in many fields of science and industry, and are

**Citation:** Voss, A.; Schroeder, R.; Schulz, S.; Haueisen, J.; Vogler, S.; Horn, P.; Stallmach, A.; Reuken, P. Detection of Liver Dysfunction Using a Wearable Electronic Nose System Based on Semiconductor Metal Oxide Sensors. *Biosensors* **2022**, *12*, 70. https://doi.org/10.3390/ bios12020070

Received: 29 December 2021 Accepted: 24 January 2022 Published: 26 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of interest for numerous applications. For example, e-noses are used in the food and beverage industry to monitor processing and determine the quality of the final product [10,11], in pharmaceutical science for formulation development and quality assurance [8], and for air quality monitoring [12]. In addition, e-noses are also used in agriculture, water management, medicine, security systems, and many other fields [13].