# **Biomedical Sensing and Imaging**

Edited by Wuliang Yin and Yuedong Xie Printed Edition of the Special Issue Published in *Biosensors*

www.mdpi.com/journal/biosensors

## **Biomedical Sensing and Imaging**

## **Biomedical Sensing and Imaging**

Editors

**Wuliang Yin Yuedong Xie**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Wuliang Yin University of Manchester UK

Yuedong Xie Beihang University China

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Biosensors* (ISSN 2079-6374) (available at: https://www.mdpi.com/journal/biosensors/special issues/Biomed Sens Imaging).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-2454-2 (Hbk) ISBN 978-3-0365-2455-9 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

**Wuliang Yin** was appointed as an MT Sponsored Lecturer at The University of Manchester, U.K., in 2012, and promoted to a Senior Lecturer in 2016 and Reader in 2020. He has authored one book and more than 300 articles, and was granted more than ten patents in the area of electromagnetic sensing and imaging. Mr. Yin was a recipient of the National Scientific and Technical Progress Award by the Chinese Government in 1999, the Excellent Paper Award–Agilent Measurement Forum in 2000, the Williams Award from The Institute of Materials, Minerals and Mining in 2014 and 2015, the 2016 Best Graduate Student Paper from the IEEE I2MTC as a supervisor, the 2018 IEEE Instrumentation and Measurement Society Graduate Fellowship as a Supervisor, the 2020 Best Application Award from the IEEE Instrumentation and Measurement Society, and the IEEE gold medal as the most productive reviewer in the U.K. for IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT in 2020.

**Yuedong Xie** received a Ph.D. degree in Electrical and Electronic Engineering from the University of Manchester, Manchester, U.K, in 2016. After that, he worked as a Research Associate at the University of Manchester and a Research Fellow at the University of Warwick for 3 years. At present, he works as an Associate Professor in Beihang University in China. His research interests include electromagnetic NDT, ultrasonic phased array, electromagnetic acoustic transducers, nano-devices and energy harvesting. He was a recipient of the Best Student Paper Award (First Place) and the Student Travel Grant Award from the IEEE International Instrumentation and Measurement Technology Conference in 2016. He is the author of more than 30 peer-reviewed publications, and he has also given several invited oral presentations in the last five years.

## **Preface to "Biomedical Sensing and Imaging"**

Biomedical sensing and imaging technology are of great significance for the early detection, rapid diagnosis and precise treatment of diseases. In the last three decades, considerable research efforts have been devoted to biomedical sensors and imaging systems. Biomedical sensing and imaging systems can be roughly divided into optical imaging technology, electrical imaging technology, acoustic imaging technology, radiographic imaging technology, etc. As each technology has its own unique advantages, various new types of sensing and monitoring equipment have been invented, and improving the performance of the existing equipment is essential. In recent years, biomedical sensing and imaging have been applied in many aspects of biomedical research, such as biomolecular testing, health monitoring, and disease diagnosis and treatment. Biosensors convert biomedical signals to electrical signals for further acquisition and processing using downstream devices and algorithms, which are also crucial for rapid and accurate biomedical diagnosis. At present, biosensor design emphasizes portability, networked detection, integration and intelligentization. Moreover, research methods are being developed regarding multi-modality and multi-function fusion, which has become a trend in the development of biological imaging systems. Some multi-modal imaging equipment has been successfully used in clinical applications, which can provide more accurate and high-resolution anatomical, functional and biochemical information. With the development of computer technology, physics, molecular biology, materials chemistry and other disciplines, new software and hardware devices have been developed, opening up a new world for the development and application of biological imaging technology. This book, which is a Special Issue of Biomedical Sensing and Imaging, aims to collect recent advances in the above topics regarding biosensing and imaging, especially in wearable/smart biosensors and advanced imaging algorithms.

> **Wuliang Yin, Yuedong Xie** *Editors*

### *Review* **Biomedical Applications of Electromagnetic Detection: A Brief Review**

**Pu Huang 1, Lijun Xu <sup>2</sup> and Yuedong Xie 1,2,\***


**Abstract:** This paper presents a review on the biomedical applications of electromagnetic detection in recent years. First of all, the thermal, non-thermal, and cumulative thermal effects of electromagnetic field on organism and their biological mechanisms are introduced. According to the electromagnetic biological theory, the main parameters affecting electromagnetic biological effects are frequency and intensity. This review subsequently makes a brief review about the related biomedical application of electromagnetic detection and biosensors using frequency as a clue, such as health monitoring, food preservation, and disease treatment. In addition, electromagnetic detection in combination with machine learning (ML) technology has been used in clinical diagnosis because of its powerful feature extraction capabilities. Therefore, the relevant research involving the application of ML technology to electromagnetic medical images are summarized. Finally, the future development to electromagnetic detection for biomedical applications are presented.

**Keywords:** electromagnetic detection and biosensors; electromagnetic biological theory; biomedical application; frequency; machine learning

**Citation:** Huang, P.; Xu, L.; Xie, Y. Biomedical Applications of Electromagnetic Detection: A Brief Review. *Biosensors* **2021**, *11*, 225. https://doi.org/10.3390/bios11070225

Received: 28 May 2021 Accepted: 3 July 2021 Published: 7 July 2021

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

**Copyright:** © 2021 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/).

#### **1. Introduction**

Recently, the biological detection technology has developed rapidly. It can be used to detect the relevant parameters that characterize the performance of the organism due to their unique identification ability to the organism [1]. When the organism is stimulated externally, the organism converts these stimulation signals into electrical signals that can be received and processed. This attribute allows it to obtain nutrients or stay away from danger. Human beings use the sensitivity of biometrics to observe and understand the living environment. In other words, the core of biological detection is converting observable things into measurable physical quantities as the result of biological simulation. In fact, biological detection technology has made greater breakthroughs in bio-electric signal conversion because of the continuous integration and penetration of biotechnology, chemistry, biology, physics, electronics, and information technology [2–4]. According to the detection mechanism, biological detection technology can be divided into electromagnetic biological detection technology, electrochemical biological detection technology, optical biological detection technology, semiconductor biological detection technology, thermal biological detection technology, piezoelectric crystal biological detection technology, etc.

Compared with other detection technology, electromagnetic biological detection technology has received attention from an increasing number of researchers due to the characteristics of high sensitivity and high detection speed [5]. The basic principle of electromagnetic detection technology is electromagnetic biological effect. Specifically, organisms can produce biological tissue effects when they are affected by the external electromagnetic field. People can get relevant information based on the relevant reactions of organisms. Since the inextricable connection between biological electromagnetic effects and life systems, it can be used to explore the positive effects of electromagnetic fields to promote the development of humanity.

More important, the biological action mechanism and biological effects of electromagnetic fields are closely related to the excitation frequency, intensity, and duration of electromagnetic fields [6]. Different electromagnetic fields have different biological targets, scope of action, action time, and energy and information exchange methods [7]. Therefore, electromagnetic detection has a wide range of applications. For example, electromagnetic biosensors have played an important role in increasing the germination rate of seeds [8], improving the adaptability of organisms to the environment [9], controlling environmental pollution [10], etc. At the same time, bioelectromagnetic detection is also applied in body health monitoring, especially the emergence of nuclear magnetic resonance technology, which has pushed bioelectromagnetic detection to a new peak. In addition, treatment methods based on electromagnetic biological effects have emerged one after another, such as magnetized water and magnetic acupuncture. The study of bioelectromagnetic detection is of great significance to promote the level of social medical care and improve the quality of human life.

This article mainly introduces the biomedical applications of electromagnetic detection in recent years. Specifically, this review first introduces the electromagnetic biological effects and their mechanism of action. In fact, the excitation frequency and intensity are the main parameters affecting electromagnetic biological effects. Taking the excitation frequency of the electromagnetic field as a clue, the application of bioelectromagnetic detection, including health monitoring and disease treatment, is discussed. Furthermore, electromagnetic detection in combination with ML technology has been used in clinical diagnosis since the ML can extract the deep features of electromagnetic signals to obtain deep biological information. The relevant research involving the application of ML technology to electromagnetic medical images is summarized. Finally, the author presents a conclusion on the reviewed biosensing techniques, and provides an opinion on the future development trend of electromagnetic detection.

#### **2. Electromagnetic Biological Effect and Mechanism**

When the electromagnetic field acts on an organism, various substances in the organism will react to the alternating electric and magnetic fields in the way of induction, which in turn results in the biological effect of the electromagnetic field. The biological effects of electromagnetic fields can be divided into thermal effects, non-thermal effects, and cumulative effects [11]. The biological thermal effect of electromagnetic wave refers to the biological effect caused by the Joule heat converted from the electromagnetic wave penetrating the biological system [12]. Compared with thermal effects, the non-thermal effects are changed in tissues or systems, that is, there is no direct relationship with heat after the organism absorbs electromagnetic radiation energy [13]. In addition, if an organism damaged by thermal and non-thermal effects is not repaired in time, the degree of damage will accumulate when it is exposed to electromagnetic radiation again. The above process is called the cumulative effect.

A series of complex bioelectrical activities are regarded as an important part of the whole organism. The shape, structure, and function of different biological levels will change when a certain intensity of electromagnetic waves penetrate the biological system. In fact, the biological action mechanism is closely related to the frequency, intensity, and duration of electromagnetic fields [14]. The scope of action of biological objects, action time, energy and information exchange methods, and exchange values are all different due to the different characteristics of electromagnetic waves, which brings different applications. Table 1 summarizes different electromagnetic effects and their applications. In order to analyze and understand the electromagnetic biological effect, it is necessary to summarize and sort out the hypothesis of the electromagnetic biological effect mechanism.


**Table 1.** Electromagnetic biological effects and related applications.

#### *2.1. Biological Thermal Effects*

Thermal effect refers to the effect produced by electromagnetic radiation heating biological tissues or systems. The strength of the thermal effect is monotonously and positively related to the amount of electromagnetic power absorbed by the biological system. There are three ways to generate the thermal effect, namely, friction heat generation, dielectric loss heat generation, and conduction current heat generation [25].

Friction heat generation is caused by the redistribution of electric charges of water molecules, amino acids, DNA, and other biological molecules under the action of alternating electromagnetic fields [26]. It leads to changes in the structure and function of biomolecules. If the frequency of the alternating electromagnetic field is high, it will have two effects. On the one hand, various polar molecules in the biological system will undergo rapid and periodic changes. In the process of change, these polar molecules collide and rub against the surrounding molecules violently, which generates a large amount of heat energy. On the other hand, various ions in the biological system undergo periodic motion near their equilibrium position due to the action of alternating electromagnetic fields. During the movement, these ions collide with other molecules to generate heat energy.

There are positive and negative charges in the molecule. Under the action of the electric field, these positive and negative charges are subjected to electric field forces in different directions, resulting in corresponding displacements [27,28]. The coupling dipole moment and polarization direction of polar molecules changes due to charge displacement, which may lead to changes in the molecular structure. In fact, it takes time to change the polarization of the dielectric. More specifically, if the electric field changes faster than the polarization direction, the polarization of the dielectric will not keep up with the electric field change, and the electric field energy will be lost.

The last electromagnetic heating effect is conduction current heat generation. There are charged particles that can migrate freely in organisms, which makes the organisms generate an electric current under the action of an electric field; the conduction current passing through the tissue can generate heat. The generation of thermal effects generally requires a higher intensity electromagnetic radiation. At present, the research on the thermal effect mechanism is relatively clear.

#### *2.2. Biological Non-Thermal Effects*

There are weak electromagnetic fields in the organs and tissues of organisms. These electromagnetic fields are stable and orderly. The non-thermal effect of electromagnetic field will interfere with the inherent weak electromagnetic field inside of the organism. Once the inherent electromagnetic field in a balanced state is disturbed by the external electromagnetic field, the tissues and organs in the organism are damaged. Compared with the electromagnetic thermal effect, the non-thermal effect has four main characteristics. The first one is coherence, that is, the biological effect can only be produced when the electromagnetic wave parameters and the relevant parameters of the biological system meet a certain deterministic relationship. The second feature is the window characteristics, which

mainly consists of frequency windows and intensity windows. The former refers to the effect that the target in the biological system only produces certain discrete electromagnetic waves with a very narrow frequency range. The latter refers to the biological effect that the target in the biological system only produces on certain discrete electromagnetic waves with extremely narrow field strength. Another feature is synergy. The weak electromagnetic field and the metabolism of the biological system can stimulate extremely strong biological effects. The last feature is non-linearity. The intensity of the effect will not increase significantly when the intensity of the electromagnetic wave is greater than a certain threshold. Non-thermal effects can occur under extremely low frequency electromagnetic fields and low-intensity radio frequency radiation. In recent years, the non-thermal effects of electromagnetic fields on biological systems have become a research hotspot due to the extensive application of power equipment and radio frequency radiation.

#### 2.2.1. Coherent Oscillation Theory

Regarding the mechanism of the non-thermal effects of electromagnetic fields, Frohlich first proposed the coherent oscillation theory [29]. The theory believes that biological molecules with dipole characteristics have a fixed vibration frequency due to the interaction of the Coulomb gravitational force. Therefore, the polarized wave is formed, and the coherent effect is produced. When the external electromagnetic field energy is introduced into this fixed vibration, it will be stored in the oscillation mode in a purely mechanical manner. The biomolecules are strongly excited and stay away from the thermal equilibrium state. If the external electromagnetic field close to the coherent oscillation frequency of the biological system acts on the biological system, it will produce two kinds of "frequency window" effects. These two kinds of effects are destructive coherence and constructive coherence, which leads to obvious biological reactions.

#### 2.2.2. Ion Transmembrane Cyclotron Resonance Theory

The theory of ion transmembrane cyclotron resonance is another theory referring the non-thermal effects. The theory considers that charged ions in a static magnetic field move in a circular motion under the action of the Lorenz force [30,31]. In the electromagnetic radiation, the alternating magnetic field parallel to the static magnetic field will increase the angular velocity and orbit radius of the charged particles, which in turn produce cyclotron resonance. When the cyclotron resonance frequency is consistent with the ELF frequency of the intracellular Ca2<sup>+</sup>, it will induce the receptors on the cell surface and the ions in the transmembrane channel to move circularly or spirally. This phenomenon will affect many enzymes regulated by calmodulin and induce various physiological and biochemical changes. At the same time, the electromagnetic field changes the dipole moment of the ion channel, and the cyclotron motion of the ion in the channel will eventually interfere with the process of ion permeability.

#### 2.2.3. Free Radical Mechanism Theory

The currently commonly recognized mechanism is the free radical pair mechanism [32]. According to the quantum theory, the chemical bond formed between two reactive free radicals is requested to be in a single state. Since electrons have a magnetic moment, the local magnetic field generated by other electrons or atomic nuclei in the molecule may reverse the spin direction of an electron. The radical pair changes from a single state to a triplet state. More important, the external magnetic field will maintain the magnetic moment of the electrons and reduce the possibility of spin reversal, thereby affecting the reaction rate. Although the experiment shows that the reaction exists in biological systems, its biological significance is not yet fully defined [33].

#### *2.3. Cumulative Effects*

It takes a certain time for the magnetic field to change the metabolism, function, and structure of the organism. This phenomenon is called hysteresis. The main reason for the hysteresis is the existence of various dipoles, ion genes, and inorganic ions in the biological macromolecule [34]. They are restricted by steric hindrance and potential energy. In addition, only a few stable equilibrium positions are available for the particles in space. They can jump from one equilibrium position to another under the effect of external electric fields. There is also an electric dipole moment between each other because ionic groups and complex ions exist in the form of different chemical bonds. When the electric dipole moment rotates towards the outward electric field, a certain potential barrier must be overcome, so there must be a time course. The electric dipole moment must overcome a certain barrier to rotate towards the external electric field, which must take a certain amount of time. From the energy point of view, a higher the product between the field strength and the action time increases the biological effect, which reflects the cumulative process.

#### **3. Application of Electromagnetic Detection Technology at Different Frequencies**

The utilization of electromagnetic biological effects involves many aspects, including disease treatment [35], food quality inspection [36], food sterilization [37], plant cultivation [38], and so on. Since the main parameters affecting the electromagnetic biological effects are frequency and intensity, electromagnetic waves at different frequencies produce different biological effects. According to the interaction mechanisms between biological tissues and electromagnetic fields at different frequencies, the radiated electromagnetic fields are usually divided into four frequency bands: electrostatic field (0 Hz), extremely low frequency (ELF, 0–300 Hz), intermediate frequency (300 Hz–10 MHz), and radio frequency (RF, 10 MHz–300 GHz). The frequency band of the electromagnetic field is shown in Figure 1. Therefore, this review discusses the related application of electromagnetic biological detection by using their working frequency band as a clue.

**Figure 1.** The frequency band of the electromagnetic field.

#### *3.1. High-Voltage Electrostatic Field (0 Hz)*

In recent years, high-voltage electrostatic fields have been widely applied in agricultural production and food processing. This technology has an effect on water molecules, enzymes, and other substances inside crop seeds to optimize their original properties and achieve good conditions for seed growth. When the seeds are treated with high-voltage electrostatic fields, ion mist, *O*3, and nitrogen oxides will be produced. Nitrogen oxides on the surface of the seed compound with hydration to form nitric acid and nitrous acid, which can stimulate embryos in the dormant state of the seed. Moreover, the *O*<sup>3</sup> is a strong bactericide, which can greatly reduce the incidence of plant diseases. Therefore, the application of electrostatic field helps to obtain high-yield, high-quality, and non-toxic agricultural products [39]. As for the aspect of food processing, the high-voltage electrostatic field technology can not only be used for the fresh-keeping of fruits and vegetables, but also in food sterilization and thawing to meet human dietary needs.

#### 3.1.1. Seed Germination and Growth

The earliest and most extensive field of electrostatic biological effects research is the use of high-voltage electrostatic fields to treat plant seeds. In 1962, Krueger et al. first discovered that the air environment ionized by positive or negative ions can significantly boost the germination rate of oat seeds [40]. In addition, the growth rate and dry weight of oat seeds after germination are also significantly improved. The experiment equipment is shown in Figure 2 [41]. Sidaway et al. placed lettuce seeds between energized aluminum electrodes for 24 h [42]. The results showed that both positive and negative electrostatic fields from 180 V/cm to 360 V/cm can change the germination rate of lettuce seeds. Figure 3 is the result of high-voltage electrostatic field treatment of safflower seeds, which further proves the effectiveness of the technology [43]. More important, they found that the positive electric field has an inhibitory effect on plant growth, while the negative electric field has a promoting effect. Huang R. et al. treated cucumber seeds with swelling treatment and electric field treatment, respectively, and then carried out related experiments about germination and membrane permeability [44]. The results showed that these two treatment methods can improve the germination rate of cucumber seeds. In subsequent studies, it was further discovered that the combination of a high-voltage electrostatic field and magnetic field technology has a more significant effect on increasing the germination rate of plants. Based on this technology, Li et al. applied a high-voltage electrostatic field with an intensity of 2.25–2.5 kv/cm to tomato seeds, while applying an alternating magnetic field over 8 h [45]. It was found that the germination rate of tomato seeds was increased by 1.1 to 2.8 times compared to the control group. Shao applied Helmholtz coils with the intensity of 2 A to induce the magnetic field. The magnetic field with discharge plasma can improve the germination rate and activity of old spinach seeds. In fact, the germination rate and activity of seeds depend on the flux of the excitation magnetic field [46]. Xu et al. applied corona discharge produced by a multi-needle-plate HVEF to treat naked oat seeds. According to the Fourier Transform infrared spectroscopy (FTIR) of the seed coat, highvoltage electrostatic fields treatment can form a new absorption peak at 1740 cm<sup>−</sup>1, which is closely related to the hydrophilicity of the seed; that is to say, the high-voltage electrostatic fields treatment can improve the hydrophilicity of seeds [47].

**Figure 3.** Seed germination processed by a high-voltage electrostatic field. Reprinted with permission from ref. [43]. Copyright 2006 Elsevier.

A high-voltage electrostatic field not only improves the germination rate of seeds, but also improves the physiological indicators of seedlings. Kiatgamjorn et al. found that the germination rate of seeds irradiated with electric fields of 10 kv/m and 25 kv/m are improved, and the seedling height and root length of seedling are significantly higher than seedlings of the control group [48]. In addition, the effect of high-voltage electrostatic field on plant seeds is also manifested in improving the germination index, vigor index, and the activity of various enzymes. Isobe et al. used hydrogen spectrum technology to find that high-voltage electrostatic field treatment can enhance the combined water content of morning glory seeds [49]. In fact, different plants have different sensitivities to a high-voltage electrostatic field environment, and the unsuitable electric field environment may have a negative effect on seed germination. Although Isobe et al. applied a highvoltage electrostatic field of suitable strength to increase the germination rate of morning glory seeds, the germination rate of morning glory seeds is significantly reduced when the electric field intensity is 500 kv/m. A suitable high-voltage electrostatic field can improve seed vigor and make seedlings grow vigorously, thereby improving the ability of seed growth to absorb nutrients and transform and utilize energy substances.

#### 3.1.2. Food Thawing

As shown in Figure 4, the high-voltage electrostatic field thawing technology is a new type of non-thermal processing method. In the small area around the needle electrode, the electric field is used to accelerate the generation of ions. The generated momentum is transferred from the air ions to the neutral air molecules, and corona wind drives a large amount of fluid to the surface of the material, which in turn leads to an increase in the heat transfer coefficient [50,51]. Foods placed between the electrodes are processed by the applied electromagnetic field to achieve thaw. Therefore, the high-voltage electrostatic field thawing has the characteristics of high efficiency and energy saving. Compared with air thawing, high-voltage electrostatic field thawing can shorten the thawing time by half without affecting the thawing quality. Hsieh et al. used the electrostatic field to thaw tilapia meat [52]. They found that the high-voltage electrostatic field can effectively reduce the cooking loss of chicken by 1% to 3% and increase the water holding capacity by 1% to 2%. He et al. discovered the current under negative voltage is greater than the current under positive voltage for the same applied voltage. The corona wind speed changes within a certain range under a given applied voltage and electrode distance, and it is more stable under negative voltage than under positive voltage. Moreover, the current increases with the increase of the applied voltage, and since the corona wind speed increases with the square root of the current, the thawing rate increases under a high voltage and large current [53]. The rapid thawing characteristics of the high-voltage electrostatic field have also attracted the attention of scholars. Dalvi-Isfahan et al. investigated the voltages and freezing temperatures of thawing meat in a high-voltage electrostatic field [54]. The

results demonstrated that a slight increase in the hardness value could be seen with an increased voltage or decreased freezing temperature. However, only changing the intensity of the voltage plays a significant role in the equivalent diameter, hardness, and drip loss. Li et al. used the negative electrons generated by the high-voltage electrostatic field to thaw the samples [55]. Defrosting under the −12 kv high-voltage electrostatic field treatment significantly reduced the water loss and microbial contamination of the fish pieces. In addition, it can enhance AMP-deaminase activity, reduce ACP activity, and delay IMP degradation. Therefore, the method can not only effectively maintain the freshness of the meat, but it also has a faster thawing speed. Rahbari et al. discovered that the fibrils of chicken breast stayed relatively intact without any gaps under 2.25 kv/cm. However, the largest deformation and loss of collagen fibrillation network appear under 3 kv/cm. The full results are shown in Figure 5 under different intensities voltages [56].

**Figure 4.** The schematic diagram of high-voltage electrostatic field thawing technology.

**Figure 5.** Sections of fresh chicken breast. (**a**) Thawed chicken breast by air method. (**b**–**d**) Thawed chicken breast by a high-voltage electrostatic field: (**b**) 1.5 kv/cm (**c**) 2.5 kv/cm (**d**) 3 kv/cm. Adapted with permission from ref. [56]. Copyright 2020 Wiley.

As shown in Table 2, the high-voltage electrostatic field has a significant effect on thawing meat products according to the related indicators. Some scholars have discussed the quality of food thawing on the structure and size parameters of the electrode plate. Amiri et al. determined that increasing the number of electrode needles can improve the oxidation degree of beef by fixing the electrode plate spacing [57]. Mousakhani-Ganjeh et al. showed

that the oxidation degree of tuna will increase as the distance between electrode plates decreases [58]. Jia et al. found that there is no significant difference between high-voltage electrostatic field thawing and air thawing [59]. Consequently, the oxidation problem caused by the high-voltage electrostatic field can be solved by selecting the appropriate voltage intensity, plate spacing, or the number of electrode plates.

**Table 2.** Comparison of the thawing meat products under high-voltage electrostatic field.


#### 3.1.3. Food Preservation

The mechanism of electrostatic field preservation is quite different from that of plant growth. The electrostatic field plays a positive role in the germination of seeds and plants to a certain extent, which strengthens the metabolic process [60]. However, preservation is to delay the aging process of fruits and vegetables by inhibiting the respiratory metabolism of fruits and vegetables. There is still a big controversy in the preservation mechanism. Some researchers believe that the electric field affects the metabolic activity of biological cell membranes by changing the transmembrane potential. Others hold the view that the biological electric field inside the fruits and vegetables affects the electron transporter in the respiratory system.

Zhao et al. used a high-voltage electrostatic field to treat green mature tomatoes for 2 h at 20 ◦C. The high electric field apparatus contains a power generator, treatment chamber, two stainless steel plates, one high-voltage controller, one voltmeter, and one amperometer. It can output −30–0 kv and 0–30 kv maximum, respectively. The results shown that the respiration rate of vegetables after electric field treatment decreased, which prolonged its storage time [61]. Kao et al. treated fresh-cut broccoli with 50–400 kv/m high-voltage electric field at 4 ◦C. The results show that the high-voltage electric field can prolong the storage time of fresh-cut broccoli according to the analysis of color and vitamins [62]. Takaki et al. introduced vegetables and fruits with a high-voltage electrostatic field. Compared with the control group, the storage time of the juice after the treatment was prolonged [63]. As for meat products, Wen-Ching Koa et al. used 100 kv/m high-voltage electrostatic treatment on tilapia. The equipment included an electrostatic generator, whose output voltage was 50 kv. This technology can delay the water and salt soluble to prolong the fresh-keeping period of fish by ~2 days [64]. Moreover, the 300 kv/m high-voltage electrostatic can inhibit actomyosin Ca2+-ATPase activity. If a 600 kvrt/m high-voltage electrostatic is applied, a longer extension time can be obtained.

#### *3.2. Extremely Low Frequency Electromagnetic Field (0–300 Hz)*

The extremely low frequency electromagnetic field (ELF-EMF) is the most common electromagnetic field because the power frequency of most countries in the world is in this frequency band. ELF-EMF has a certain effect on humans, including immune function, cell differentiation, cell proliferation, and transmission of intracellular signal substances. It can achieve cell differentiation, tissue repair, differentiation of stem cells, and healing of injured cells. Specific ELF-EMF can stimulate the production of local bone factors, and significantly improve the proliferation and differentiation indicators of bone density. In addition, ELF-EMF can also inhibit tumor growth and induce cell apoptosis [65,66]. More important, the human nervous system is the most sensitive issue to ELF-EMF. A large number of experiments have shown that the effect of ELF-EMF on the central nervous system is double-sided. On the one hand, the ELF-EMF may affect the normal brain

development. It has an important impact on the growth and differentiation of neurons and the direction of neurites. On the other hand, electromagnetic fields may cause oxidative stress in the brain tissue, which leads to the activation of astrocytes. Based on the above research, this review focuses on the application of extremely low frequency electromagnetic sensors in cancer treatment, fracture healing, peripheral nerve repair, and transcranial magnetic stimulation.

#### 3.2.1. Cancer Treatment

ELF-EMF can induce apoptosis of cancer cells to some extent, thus providing a new way for the treatment of cancer. Studies have shown that ELF-EMF have significant therapeutic effects on a variety of tumors. They can induce growth inhibition and tumor killing effects on lung cancer, liver cancer, breast cancer, prostate cancer, and melanoma cells [67–69]. The main anti-tumor mechanisms of electromagnetic fields can be roughly divided into the following four categories:


The rotating electromagnetic sensor has a significant effect on tumor suppression. In 2011, Wang et al. found that the magnetic field with an intensity of 0.4 T and excitation frequency of 7.5 Hz for 4 days can inhibit the proliferation of a variety of human tumor cells, including gastric cancer BGC-823, gastric cancer MKN-28 cells, and lung cancer A549 [70]. However, this technology does not inhibit poorly differentiated gastric cancer MKN-45 cells and lung adenocarcinoma SPC-A1 cells. It cannot inhibit poorly differentiated gastric cancer MKN-45 cells and lung adenocarcinoma SPC-A1 cells. Ren et al. found that the rotating magnetic field with an intensity of 0.4 T and excitation frequency of 7.5 Hz for 6 days can induce the senescence of lung cancer cell A549 and inhibits its iron metabolism [71]. If the magnetic field intensity and frequency are changed, for example, the rotating magnetic field with an intensity of 0.4 T and excitation frequency of 7.5 Hz for 15 min does not have a significant effect on the proliferation of HepG2 liver cancer cells. The research shows that the intensity of the magnetic field, the processing time and the cell type are all key factors for the effect of magnetic field on cells.

Zhang et al. designed a gyromagnetic sensor composed of two groups of 30 neodymium iron boron permanent magnets which are anti-parallelly stacked. The structure of sensor presented in Figure 6a [72]. The magnet array rotates on the vertical axis at a frequency of 8–10 Hz. The MF intensity at the top of the processing table is 0.6 T, and the MF intensity at 5 cm above the surface is 0.4 T. Moreover, the horizontal spatial change in the animal cage region is 0.6–0.32 T. The rotation introduces an approximately sinusoidal 0.21 T average fluctuation in the MF intensity with a frequency of 8–10 Hz. As shown in Figure 6b, Du et al. investigated the gyromagnetic sensor with a cylindrical permanent magnet rotating around the long axis [73]. Some scholars studied the gyromagnetic sensor composed of NdFeB Halbach array magnets. As shown in Figure 6c,d, a uniform 1 T magnetic field is generated radially in the center part of the device, and a magnet is installed on the motor to control its rotation [74]. Under a rotating magnetic field, the internalized MPs can effectively destroy the cell membrane through mechanical force and promote programmed cell death in vitro and in vivo. Magnetic field therapy successfully reduced the size of brain tumors and prolonged survival rates. The development of the above-mentioned sensors can significantly improve the uniformity of the magnetic field and accelerate the induction of apoptosis of some cancer cells to a certain extent.

**Figure 6.** The structure of the gyromagnetic sensor. (**a**) Two anti-parallel stacks of 30 neodymium-iron-boron permanent magnets. Reprinted with permission from ref. [72]. Copyright 2005 Wiley-Liss. (**b**) Cylindrical permanent magnet rotating around the long axis. (**c**) NdFeB Halbach array magnet. (**d**) The changes of the magnetic field corresponding to (**c**).

Sisken et al. used pulsed electromagnetic fields with an intensity of 0.3 mT and a frequency of 2 Hz to treat rat sciatic nerve defects, which proves that electromagnetic fields can stimulate the regeneration of the nervous system [75]. Experimental results show that the application of pulsed electromagnetic fields can accelerate the formation of myelin and connective tissue [76]. Zhong et al. measured the proliferation and differentiation of bone marrow mesenchymal stem cells under condition of the electromagnetic field with an intensity of 0.5 mT and frequency of 50 Hz [77]. The results confirmed that prolonging the stimulation time can promote the cell proliferation and induce the cell differentiation. After transplantation into animals, the original cells of bone marrow mesenchymal stem cells stimulated by the magnetic field performed obvious biological characteristics. Diniz et al. used a Helmholtz copper coil to analyze the effect of pulsed magnetic fields on the proliferation of osteoblasts, which also confirmed the effectiveness of electromagnetic stimulation. The diameter and number of turns are set to 12 cm and 250, respectively. The excitation frequency of the coil is set to 15 Hz [78]. Liu et al. designed a set of ELF-EMF biosensors consisting of a driving unit and a pair of Helmholtz coils. The diameter of the coil is 14 cm, and the distance between the coils is 13 cm. The induced peak voltage and the change intensity of the magnetic field are ±0.7 mV and 0.4 T/s, respectively [79]. They found that the pulsed magnetic field generated by the electromagnetic sensor can protect the synthesis of chicken sternum cartilage proteoglycan in the embryo.

#### 3.2.2. Transcranial Magnetic Stimulation (TMS)

TMS is based on the principle of electromagnetic induction, which applies a rapidly changing current to the coil to generate a strong and short-lived magnetic field, as shown in Figure 7 [80]. When the voltage of the induced electric field reaches a certain strength, the resting membrane potential of the cell will change. Thereby, it has an excitatory or inhibitory effect on the nervous tissue to affect the physiological and biochemical activities of the brain, and finally achieves the therapeutic effect on neurological and mental diseases.

**Figure 7.** Magnetic field distribution stimulated by TMS. (**a**) Front view of the induced electric field distribution on the scalp surface. (**b**) Sectional view of the x-y plane.

In the application of cranial nerve stimulation, the magnetic field generates an induced electric field in the brain tissue to form an intracranial current, thereby changing the stimulation location and related neuronal activity. In 1985, Anthony Baker et al. were the first to study TMS. Nowadays, TMS has been widely used in clinical experiments and research in the fields of neurology, neuroscience, psychiatry, rehabilitation, psychocognitive science, and sports medicine [81]. It has a good therapeutic effect in depression, obsessivecompulsive disorder, schizophrenia, and peripheral nerve rehabilitation [82].

Jiang et al. used TMS to treat 40 patients with chronic primary insomnia. The results show that TMS treatment significantly improved the sleep cycle of stage III sleep and rapid eye movement sleep [83]. Levitt et al. used 10 Hz TMS to treat patients with refractory depression [84]. Magnetic resonance spectroscopy is used to assess the changes in GABA levels of the left dorsolateral prefrontal cortex (DLPFC). The results indicate that TMS treatment is associated with elevated GABA levels, and subjects who receive GABA agonist treatment have less response to TMS. Leblhuber et al. studied the changes in the utilization of neurotransmitter precursor amino acids in elderly patients with depression after stimulation of the prefrontal cortex, and found that TMS has a significant effect on the treatment of depression [85].

The performance of the TMS system is mainly determined by the focus, stimulation depth, and intensity of the induced electric field. In fact, the structure of the TMS coil is directly related to the electric field distribution in the brain tissue. Among them, the depth and focus of the induced electric field distribution are the two most critical indicators for the optimal design of the TMS coil [86]. In order to locate the optimal stimulation point of the sensor to the intracranial cerebral cortex, the structure of the sensor needs to be precisely designed. It is helpful to make continuous breakthroughs in the performance of TMS by simulating the performance of the coil in terms of stimulus intensity, stimulus depth and focus.

(1) Circular coil

The structure of the circular coil is simple, and the induced electric field on the surface of the brain is distributed in a "circular shape." The stimulation performance of the circular coil is poor, especially in focus behavior. In order to improve the focus, one side of the coil can be folded up along the center of the circle. The novel structure of coil is called the deformed circular coil. The simulation results show that the magnetic field induced by the deformed circular coil has a large attenuation. The change in the magnetic field leads to a change in the electric field distribution, which in turn results in an improved focusing. Figure 8a shows the structure of a circular coil and a deformed circular coil.

**Figure 8.** The structure of the TMS coil. (**a**) Circular coil. (**b**) Figure-of-eight coil. (**c**) H-shaped coil.

#### (2) Figure-of-Eight coil

Figure-of-eight coils are widely applied in commercial applications. Two coils with opposite currents cause the peak of the induced electric field to appear at the tangent position of the coil [87,88]. Figure 8b shows the figure-of-eight coil. Compared with the circular coil, the focus of the figure-of-eight coil is greatly improved. However, it can only be applied to the stimulation of the superficial areas of the brain. In order to further improve the performance of the figure-of-eight coil, a butterfly coil and a quadruple butterfly coil are proposed. The structure of the butterfly coil and quadruple butterfly coil are improved based on the figure-of-eight coil. By enhancing the magnetic field at the tangent of the coil, the focusing performance is improved and the target point has greater stimulation intensity. In addition, the quadruple butterfly coil can be used in combination with the shielding plate, which can further reduce the focal area of the brain surface.

#### (3) H-shaped coil

The H-shaped coil does not cause high electric field strength in non-target areas because it has a small electric field attenuation rate and surface charge accumulation. As shown in Figure 8c, it can effectively stimulate the deep neuron area, so it belongs to the deep stimulation coil. In addition, the penetration depth of each area can be flexibly controlled by adjusting the coil input or changing the distance between the coil and the brain [89,90].

(4) Halo coil

A large diameter circular coil placed in the middle of the head is called the Halo coil. The special position of the Halo coil enables it to stimulate deeper tissues in the brain; thus, the Halo coil belongs to the deep stimulation coil category. Figure 9 and Table 3 show the structure and model characteristics of the Halo coil, respectively. Halo coils can also be combined with other coils, such as the Halo-circular assembly coil (HCA coil), Halo coil working with two circular coils (HTC coil), Halo and Fo8 coils (HFA coil), double-cone coil with the Halo coil (HAD coil), and the triple-halo coil (THC coil) [90–92]. Although different coils have different structures and characteristics, they all realize the superposition and cancellation of the magnetic field by changing the direction of the input current, which increases or reduces the intensity of the induced electric field.

**Figure 9.** The structure of the Halo coil and related combinations. (**a**) Halo coil. (**b**) HCA coil. (**c**) HTC coil. (**d**) HFA coil. (**e**) HDA coil. (**f**) THC coil.


**Table 3.** Characteristics of the Halo coil model.

Moreover, there are linear coil arrays, cap-formed coil arrays, multi-layer coil arrays, etc. [93,94]. The structures of the coil arrays are shown in Figure 10. The linear coil array has a small number of input control units, and the stimulation mode is flexible and changeable by changing the size and direction of the input current. The cap-formed coil array has a greater stimulation depth and weakens the impact on non-target areas due to the increased coupling with the head surface. The multi-layer coil array increases the stimulation intensity to a certain extent. As the number of coils increases, the number of channels of the magnetic field generator should also increase accordingly. In other words, the overall structure and working conditions are more complicated, and hence, the multi-layer coil array has poor applicability.

**Figure 10.** The structure of the TMS coil arrays. (**a**) Circular coil array. (**b**) Linear coil array. (**c**) Capformed coil array. (**d**) Multi-layer coil array.

#### *3.3. Intermediate Frequency Electromagnetic Field (300 Hz–10 MHz)*

Intermediate frequency electromagnetic detection are widely applied in areas such as food quality and spatial location monitoring. In food quality monitoring, the organization, composition, structure and state of meat products, fruits, and vegetables are closely related to their electromagnetic properties [95]. For example, the impedance value has a certain change rule in the process of meat decay. In addition, the impedance value of the electromagnetic sensor is also closely related to its position, so the electromagnetic sensor can detect changes in the relative position in space. At present, this technology is mainly used for the monitoring of human joints, limbs, and eye movements. The applications of electromagnetic sensors in food quality and spatial location monitoring are introduced as follows.

#### 3.3.1. Food Quality Inspection

The decay of meat products produces many complex biochemical and physical and chemical processes, which will affect its organization and structural characteristics [96]. These phenomena also change the electrical parameters. Biological tissue can be considered as a nonmagnetic material, and its electrical properties are the result of the joint action of its tissue and structure. Among them, the most important component that affects the dielectric properties of the tissue is the ion. As for the structure of the dipole moment, the main factors are water tissue molecules, proteins, and lipids. The movement of electric charges causes conduction effects, and the polarization of dipoles produces dielectric relaxation. The conductivity of most tissues is low when electromagnetic waves are excited at low frequencies, which mainly depends on the volume fraction of extracellular fluid. As the excitation frequency increases, the conductivity of intracellular and extracellular ions rises sharply due to the dielectric relaxation of water. The whole process can be represented by Figure 11.

**Figure 11.** The dispersions of biological samples.

The above-mentioned increase in conductivity is called dispersion. In biological systems, there are mainly three kinds of dispersions (*α*, *β*, *γ*). According to the phenomenon of dispersion, it is divided into three frequency bands. The frequency range of α dispersion is from millihertz to hundreds of hertz, and it is mainly related to the phenomenon of tissue polarization. The frequency range of *β* dispersion caused by the Maxwell-Wagner effect is from hundreds of hertz to megahertz. It is a method of measuring the integrity of cell membranes during meat aging because of reduced insulation properties. The frequency range of the *γ*-dispersion is in gigahertz, which is mainly caused by the permanent dipole relaxation of free water in muscle tissue. This characteristic also further proves the reliability of the electrical method for detecting meat products. In fact, electrical detection can be divided into contact measurement methods and non-contact measurement methods. The response speed of the electrode contact measurement method is fast and accurate. However, the contact method may cause unnecessary measurement errors, such as the

chemical reaction and food tissue damage. This problem does not exist in non-contact measurements. As a non-contact measurement method, the electromagnetic detection method has a wide range of applications in food quality detection [97].

Swatland et al. reported the changes in rheological properties, electrical properties, and optical properties during the ossification of pig and cattle carcasses [98]. This process can be detected by its low reactance and low capacitance, and it is found that these three physical detection characteristics have similar variation curves. Ibba et al. studied the influence of frequency on the bio-impedance of fruits using electrical impedance spectroscopy. They pointed out that the bio-impedance can be used to measure fruit quality. The results demonstrated the potential of this technique for monitoring the ripening of fruit, obtaining a good correlation of the impedance evolution with low frequency points [99]. Sun et al. combined chemometrics and impedance characteristics to create a multi-mass index for Atlantic salmon. Moreover, they applied a partial least squares method to create a quantitative prediction model of the bioimpedance spectroscopy and total volatile basic nitrogen value. The correlation coefficients of the training set and test set are 0.9447 and 0.9387, respectively, which means that the impedance characteristics measurement method has the advantages of convenience, economy, and high precision. Castro-Giraldez et al. noted that the certain frequencies (500 Hz, 300 kHz) are effective at 24 h after slaughtering to detect the dark firm dry (DFD) meats during the post-mortem period. In fact, there may be systematic errors in the measurement [100]. Kim positioned the monopolar injection needle to the dermis, subcutaneous tissue, or muscle layer of pork tissue, and then measured electrical impedance in the frequency range of 10 Hz to 10 kHz based on an impedance analyzer [101]. This method has also been proven to be used for pork quality monitoring. Tang et al. introduced a non-contact induction measurement system for bioimpedance spectroscopy [102–104]. The electrical conductivity of the food is measured by an excitationretraction electromagnetic sensor in the frequency sweeping mode, and the freshness of the food is detected according to the change of the electrical conductivity signal. Apart from that, Tang et al. used the fast finite element method to establish the cell model, as shown in Figure 12. The solution of the numerical model can obtain the vortex field distribution inside the cell and further proves the accuracy of the measurement results.

**Figure 12.** Three-dimension finite element model. (**a**) Measurement set-up. (**b**) Conductivity spectroscopy of fresh potato and defrosted potato. (**c**) Simulation model. Adapted with permission from ref. [102]. Copyright 2019 Elsevier.

#### 3.3.2. Spatial Position Measurement

Electromagnetic sensors can be used for spatial position measurement. When a small receiving sensor moves in space, the electromagnetic positioning system can accurately calculate its position and orientation, thereby providing dynamic and real-time measurement of the position and orientation angle [105,106]. In biomedical research, it is an ideal choice for measuring range of motion and limb rotation due to the advantages of being fast and precise. The electromagnetic tracking system generally consists of a magnetic field transmitter and receiver. A magnetic field generator is placed near the target, and the magnetic field covers a certain range. The receiver detects the strength and phase of

the magnetic field and sends the signal to the control unit. The control unit can calculate multiple degrees of freedom for tracking the target. The advantage of the electromagnetic tracking system is that the detection result is not restricted by the line of sight. Nothing can block the tracking of the electromagnetic tracking system except for electrical conductors or permeable magnets. The disadvantage is that it is easily interfered by metal and the working range is small.

Recently, a three-dimensional electromagnetic tracking system was launched, as shown in Figure 13. It can measure the range of motion of the lumbar spine, cervical spine, and shoulders of humans. The "Flock of Birds" electromagnetic tracking device is used for three-dimensional motion analysis, and the patient's arm is extended in three directions [107]. The slope of the regression line between the rotation of the scapula joint and the rotation of the scapula is calculated, which reflects the rhythm of the scapula. According to the measurement result, the movement information of the shoulder joint can be obtained in time.

**Figure 13.** Three-dimensional electromagnetic tracking system. (**a**) Schematic diagram of detection device. (**b**) Movement measurement signal. Adapted with permission from ref. [107]. Copyright 2018 Elsevier.

In accordance with the idea, Xie et al. proposed a new method of biomechanical motion monitoring based on electromagnetic sensing technology [108]. The proposed method has the advantages of non-invasiveness, easy implementation, low cost, and high efficiency. A theoretical model was established to model the sensor response, and the wearable sensor system and corresponding measurement signal shown in Figure 14 was designed to monitor various biomechanical movements, including blinking frequency, finger/wrist bending level, and frequency. The whole system is measured by the mutual impedance of the coil. According to the change frequency of the signal per unit time, it can reflect the body's movement state in time. Therefore, this technology can be used to detect early signs of these abnormal biomechanical behaviors, so that relevant treatment can be carried out in time. Robison invented an electromagnetic method of eye movement recording, which has the characteristic of low noise; however, the method is limited by blinking and fluctuations in pupil size [109]. Remmel designed a cheaper version. The system uses three alternating magnetic fields (48, 60, and 80 kHz) in the X, Y, and Z directions to record eye movement information [110]. The subject wears a contact lens with a small loop of wire, and his eyes are located close to the center of the test field. The three alternating magnetic fields introduce three voltages into the coil through Faraday's law of induction. As the eye rotation, more magnetic field passes through the coil, which generates induces larger voltage. Three analog (direct current) voltages of the vector components in the X, Y, and Z directions are obtained by amplifying and demodulating the coil signal. According to the measurement signal, the angle information of the eye can be inverted. Later, Remmel extended this method to arm recording [111]. The two small coils fixed at the elbow and two coils fixed at the wrist are applied to measure the movement of the arm. For each coil in the system, the voltage is demodulated to give three DC voltages to obtain the vector direction of the coil axis. At the same time, the inhomogeneity of the magnetic field is corrected. This method can describe five angles of arm movement, and the noise is only 2 arcsec rms.

**Figure 14.** The biomechanical motion monitoring based on electromagnetic biosensors. (**a**) Wearable EM sensor. (**b**) Finger bending. (**c**) Wrist bending. Adapted with permission from ref. [108]. Copyright 2020 IEEE.

#### *3.4. Radio Frequency Electromagnetic Field (10 MHz–300 GHz)*

Radio frequency electromagnetic fields (RF-EMFs) generated by communication base stations, mobile phones, wireless phones, and wireless routers have become an indispensable part of modern life. The non-thermal effects of radio frequency electromagnetic fields have a serious impact on daily lives [112]. Although the non-thermal effects of radio frequency electromagnetic fields have a serious impact on people's daily life, electromagnetic fields can react with biological tissues for disease detection. Therefore, this review will introduce some applications and developments of high-frequency electromagnetic detection in disease detection.

#### 3.4.1. Disease Detection Based on MRI

Nuclear magnetic resonance (NMR) technology is mainly based on the interaction between magnetic nuclei and magnetic fields. The nuclear magnetic moment of the lowenergy state absorbs the energy provided by the alternating field strength, and then

continuously transforms to the high-energy state, thereby generating nuclear magnetic resonance signals [113]. This transition phenomenon caused by resonance absorption is called the nuclear magnetic resonance phenomenon.

MRI is based on the nuclear magnetic resonance signal emitted by the hydrogen nuclei in the human body. The strength of this signal depends on the position of the hydrogen nucleus in the molecular structure and the surrounding environment, so it can reflect a large amount of media information. The strength of this signal depends on the position of the hydrogen nucleus in the molecular structure and the surrounding environment, so it can reflect a large amount of media information. In the late 1980s, the phased array coil concept was introduced into the field of magnetic resonance imaging, which revolutionized the reception of MR signals [114]. The receiving array can use multiple small surface RF coils to obtain a higher signal-to-noise ratio image in a wider range by appropriately combining the signals from a single element. Since the efficiency of parallel imaging and the achievable acceleration factor depend on the number of coils in the receiving array, the number of receiving channels in the MR system has been increasing over time [115]. In fact, most research is based on phased array coil technology to reduce the coil size and increase the number of channels, which can improve the imaging quality [116]. Compared with traditional mechanical coils, wearable technology can maximize the fill factor of the magnetic resonance receiving coil to significantly improve the imaging quality [117].

The Institute of Biomedical Engineering of the University of Zurich in Switzerland has made various attempts to implement wearable technology [118,119]. As shown in Figure 15, the fabricated coil adopts an adjustable mechanical structure to adapt to the size change of the measured part. This kind of coil has a complex structure and large volume, which determines that it is difficult to have a wide range of applications. Subsequently, the team used thin wires to achieve the softness and scalability in a single direction by braiding a wire harness. In in vivo knee imaging, it was demonstrated that a stretchable array of eight receiver coils allows adaptation to different anatomical dimensions and different knee flexion angles.

**Figure 15.** The adjustable mechanical structure for MRI. (**a**) Wrist coil and corresponding detection results. Adapted with permission from ref. [118]. Copyright 2009 Wiley-Liss. (**b**) Knee joint coil, consisting of the reference coil made from circuit and stretchable coils made from 5 mm, 7 mm, and 9 mm copper braids, and detection results. Adapted with permission from ref. [119]. Copyright 2012 Wiley-Liss.

In order to explore the relationship between the number of phased array coils and the ultimate performance, the GE global research team created a 128-channel "flexible" coil for torso imaging. As shown in Figure 16a, this coil can change the shape of the coil within a certain range. However, the huge geometric size and limited deformation space are not suitable for infants and weak patients [117]. The signal-to-noise ratio (SNR) in the proposed system is 1.03 in the center of the elliptical loading model and 1.7 in the outer region on average. Furthermore, the maximum g factor of 4 × 4 acceleration is reduced to 2.0, and the 5 × 5 acceleration is reduced to 3.3. These characteristics significantly reduce the residual aliasing artifacts in human imaging.

As shown in Figure 16b, Jia et al. deployed coils on light and soft materials, but the coil array cannot be stretched. The feature of extremely flexible makes it suitable for inspections in many scenes [120]. The experiments show that the proposed MRI system can achieve a high (SNR) of 336 in vitro, and the phantom images are evenly distributed. As for in vivo experiments, it has a higher SNR of 169 in the region of interest, and the longitudinal coverage exceeds 48.5 cm. However, the bulky system line is still the biggest obstacle affecting its application. As shown in Figure 16c, Corea et al. used inkjet printing and screen-printing techniques to make coils [121]. At the same time, the relationship between the printing coil inductance and the "distributed capacitance" in different printing materials was analyzed. The results shown that the relative dielectric constant increases linearly with the concentration of barium titanate in the ink. In addition, the fabricated coil can be attached to the tested tissue, which is beneficial to improve SNR. The cost of this technology is relatively high, and the coils made by printed materials are less flexible. Therefore, materials that can meet the printing technology and high flexibility are still the bottleneck of the realization of this magnetic resonance coil.

**Figure 16.** *Conts*.

**Figure 16.** The wearable MRI sensor. (**a**) The 128-channel torso array coils and corresponding imaging results. Adapted with permission from ref. [117]. Copyright 2008 Wiley-Liss. (**b**) Lightness and softness coils and corresponding imaging results. Adapted with permission from ref. [120]. Copyright 2015 John Wiley & Sons. (**c**) Coil printing process flow. (**d**) Schematic of a printed coil. (**e**) Photograph of a printed coil.

#### 3.4.2. Disease Detection Based on Microwave

The main task of microwave detection is to simulate the scattering and diffraction of the electromagnetic field passing through the medium [122]. This technology measures the electromagnetic field scattering information of the object medium to reconstruct the dielectric properties and geometric characteristics of the object. Since the quality of the microwave image is greatly affected by the dielectric properties of the tested target, the existing microwave imaging methods are mainly used for breast cancer detection. Only a few studies are focused on the detection of strokes.

In breast cancer detection, Wang introduced a fast microwave sensor for breast cancer detection through near-field imaging methods [123]. The method of detecting cancer cells based on microwave imaging is highly correlated with the contrast of dielectric properties between healthy tissue and malignant tissue.

Bassi et al. investigated the confocal microwave imaging (CMI) for breast tumor detection. The results show that the 2D CMI system can detect small tumors (diameter 2 mm), and the 3D CMI system can identify tumors of medium size (diameter greater than 6 mm). Moreover, they applied delayed multiplication and summation methods in CMI to reduce artifacts and noise and enhance the image [124]. Yang et al. first proposed the detection method of the Microwave Image for breast cancer detection based on the principle of ground penetrating radar [125]. It uses the transmitting antenna to emit Ultra-Wideband (UWB) microwaves to irradiate the breast, and then the receiving antenna to receive the scattered signal. According to the scattered signal, the strong scattering area in the breast is recovered to determine the position of the tumor. Rappaport et al. used the time reversal algorithm of Finite-Difference Time-Domain (FDTD) to realize the spatial reconstruction of the dielectric constant distribution inside the breast. They deduced the iterative formula of FDTD to obtain the initial distribution of the electromagnetic field in space [126,127]. The method is an ill-conditioned solution problem, for which huge challenges exist in robustness and computational cost.

As for brain stroke detection, researchers at the University of Queensland in Australia have designed an ultra-wideband microwave system for head imaging [128]. According to the construction of the human brain model, the confocal algorithm is used to realize the detection of stroke. Ireland et al. used the Born iterative method to detect hemorrhagic stroke in brain tissue [129]. Furthermore, a high-contrast brain conductivity distribution map is reconstructed through microwave imaging technology. The internal information can be detected according to the distribution of brain conductivity. The reconstruction results of conductivity can be used to detect strokes.

#### **4. Application of Machine Learning Technology in Electromagnetic Medical Images**

Artificial intelligence is a technological science which can be used to simulate and expand the theories and methods of human intelligence. Machine learning (ML) is a method of realizing artificial intelligence, which applies algorithms to analyze data, make decisions, and predict specific events. Deep learning (DL) is a technology that realizes machine learning. It is a data-driven method of automatically learning high-level features hidden in images, which can greatly reduce the interference of subjective factors in feature selection [130]. In addition, the model applies a non-linear layer structure, which can be used to build a more complex model. Shallow neural networks are prone to under-fitting, while DL can improve the learning ability by increasing the depth of the network to solve more complex problems.

According to the electromagnetic biological effect, the signal detected by the electromagnetic biosensor can reflect the information of the biological body. However, manual analysis and interpretation are still required. In fact, machine learning can be used to learn the characteristics of biological signals, instead of manual operation. In particular, machine learning algorithms have been widely used in medical image diagnosis. Therefore, this review summarizes the application of machine learning in microwave breast cancer images, MRI prostate cancer images, and MRI brain tumor segmentation images.

#### *4.1. Application of Machine Learning in Microwave Breast Cancer Images*

In recent years, researchers have begun to use ML methods to directly process the signals obtained by the microwave breast cancer detection system, and thus, achieve the purpose of breast cancer screening [131]. The method of machine learning is mainly divided into two research directions. The first one is to assume that the tumor has been confirmed to exist in the breast. According to the trained model, an evaluation of the physical characteristics is needed, such as the location and size of the tumor. The other is to use the model to determine whether the tumor exists.

Under the assumption that the tumor exists, Davis et al. used a combination of Local discriminant bases (LDB) and Principal Component Analysis (PCA) to estimate the morphological characteristics of the tumor, such as the shape and size of the tumor [132]. The detection results are shown in Figure 17. Under the condition of a 10 dB signal-to-noise ratio, the method can reach an accuracy of 97% for breast cancer detection, and the shape accuracy is 70%. Conceicao et al. took the shape of the tumor as a feature [133]. The features combined with support vector machine (SVM) are used to determine whether it is a malignant tumor. They first constructed a breast prosthesis and tumors of different shapes, then placed tumors of different shapes on the breast prosthesis for measurement to obtain training samples. After that, the PCA was used to reduce the dimensionality and feature extraction of the acquired data. Finally, they analyzed the data by an SVM classifier to establish the model.

On the assumption that the tumor cannot be determined to exist, Alshehri et al. first used Discrete Cosine Transform (DCT) to extract the characteristics of the received signal of the microwave antenna [134]. Then, an artificial neural network was used as a model for judging whether the tumor was presented. They verified the effectiveness of the algorithm in uniform and uneven breast implants. The results show that the tumor existence, size, and location detection rates for both cases are highly satisfactory, which are approximately: (i) 100%, 95.8%, and 94.3% for homogeneous cases and (ii) 100%, 93.4%, and 93.1% for heterogeneous cases, respectively. Byrne et al. used PCA to extract features of the received signal on the digital breast model. They selected the SVM as classifier on the signal of each channel [135]. The final decision result was voted by the classifiers of all channels according to the majority principle. The investigated method can detect tumors as small

as 4 mm in diameter with an accuracy over 71% in a dielectrically heterogeneous breast. In 2014, Santorelli et al. used SVM linear discriminant analysis to classify breast cancer tumors, and achieved an accuracy of 76.71% for identifying signals from the tumorous phantoms [136]. The above studies have confirmed that ML technology can efficiently assist doctors in disease diagnosis. For one thing, it can reduce the burden on doctors; for another, it can reduce the subjectivity of diagnosis, thereby making the diagnosis more objective and accurate.

**Figure 17.** The classification results of the morphological characteristics of a breast tumor. (**a**) Shape. (**b**) Size.

#### *4.2. Application of Machine Learning in MRI Prostate Cancer Images*

Prostate cancer (PCa) is the cancer with the highest incidence among men. It ranks as the second-leading cause of cancer death in men. PCa has an insidious onset, but the disease progresses rapidly. The tumor cells are prone to local infiltration or distant metastasis. Therefore, early diagnosis and timely treatment have more positive meanings for PCa patients [137]. MRI is a mature and non-invasive imaging technique. Relying on the good display ability of prostate tissue structure and surrounding tissues, it occupies an important position in the diagnosis of PCa. Machine learning technology is also used in prostate cancer MRI image monitoring, mainly including prostate segmentation and cancer detection.

The segmentation of the prostate is usually a necessary step for clinical setting and further image analysis. In order to save time and increase reproducibility, it is very meaningful to realize automated prostate contour and lesion segmentation. Perfect automatic segmentation can get a fully automatic post-processing pipeline. Yan et al. proposed a model based on deep propagation neuron network by extracting the data from different levels of complexity, which can obtain more credible segmentation results of glands and their boundaries as shown in Figure 18 [138]. In addition, Alkadi et al. applied a system based on unimodal deep learning. The system can automatically segment prostate and PCa lesions [139]. Their algorithm has high segmentation accuracy and can obtain an area under the curve (AUC) of 0.995.

In terms of prostate cancer detection, ML-based computer-aided diagnosis software has important applications in the PCa field. It can improve diagnosis accuracy and repeatability. McGarry et al. used ML to detect epithelial cells and luminal density [140]. The data of 10 patients remained stable, which means the method can be used to diagnose high-grade PCa. In addition, volume of index lesion ROI can evaluate histogram parameters obtained from T2WI, DWI, and DCE images. The combination of ROI and SVM can improve the PI-RAD Sv2 score.

**Figure 18.** Segmentation results of glands and their boundaries. Red: ground truth; yellow: segmentation contour by the proposed algorithm. Adapted with permission from ref. [138]. Copyright 2018 Elsevier.

#### *4.3. Application of Machine Learning in MRI Brain Tumor Segmentation Image*

As a typical non-invasive imaging technology, magnetic resonance imaging (MRI) can produce high-quality brain images without damage and skull artifacts. MRI can provide more comprehensive information for the diagnosis and treatment of brain tumors [141]. Doctors can quantitatively analyze brain tumors to measure the maximum diameter with the help of multi-modal brain imaging to segment tumors and determine the volume and number of brain lesions. MRI brain tumor segmentation has also become an important component in the field of medical imaging because of the needs of clinical application and scientific research [142]. Researchers have proposed many methods for the segmentation of brain tumors in MRI images, but the existing methods still cannot always achieve satisfactory results due to the complexity of MRI brain images. The deep neural network model has been successfully applied to many computer vision tasks, which has received extensive attention from the scholars [143]. In specific implementation, the deep learning method segmented a large number of small-scale image blocks from the original MRI brain image. Furthermore, the image blocks are used to train a classification network. During segmentation, the corresponding image blocks are classified pixel by pixel in the context of a sliding window, and then post-processing is used to complete the segmentation of the brain tumor. This method obtains a large number of re-labeled small image blocks from the limited labeled brain image samples. Therefore, it can solve the contradiction between the few brain tumor-labeled images and the large demand for deep neural network training samples to a certain extent. The basic framework of the model for MRI brain tumor segmentation can be described as in Figure 19. The whole network architecture consisted by convolution, activation, pooling, fully connected layers, and post-processing, including conditional random fields (CRF), long short-term memory (LSTM), etc. In view of the powerful ability of automatic extraction of high discriminative features, it is quickly applied to MRI brain tumor segmentation.

The structure of a single network is used for brain tumor segmentation. The input data are processed by convolution, pooling, and nonlinear layers to obtain the feature map. Subsequently, the fully connected layer and classification layer are used to predict the label of the corresponding category based on the extracted feature. Single-network brain tumor segmentation methods can be divided into 2D and 3D due to convolution dimensions. Table 4 lists the performance of some typical single-network MRI brain tumor segmentation methods on the BraTS dataset.

**Figure 19.** The basic method of MRI brain tumor segmentation based on CNN.

The 2D convolutional neural network has the advantages of occupying less resources and fast training for brain tumor segmentation. Some researchers have carried out 2D network brain tumor segmentation research. For example, Zikic et al. proposed a 2D brain tumor segmentation network by using the Alex-Net network structure [144]. Pereira et al. used a similar VGG-Net network structure to construct an automatic segmentation method for exploring the advantages of small 3 × 3 convolution kernels [35]. In addition, they explored that small convolution kernels can design deeper network architectures, which can effectively avoid overfitting by reducing the number of parameters. Dvorak and Menze proposed a local structure prediction method to improve the local segmentation ability of MRI brain images [145]. Randhawa et al. improved boundary classification of MRI brain images [146]. In order to design a more effective brain tumor segmentation model, Ben Naceur et al. adopted the idea of multi-network ensemble learning to achieve accurate brain tumor segmentation [147]. The segmentation results are shown in Figure 20. Cui et al. focused on dividing the image block into 32 × 32 or 13 × 13-pixel scales according to the different label, so as to achieve better segmentation accuracy [148]. This method improves the accuracy of brain tumor segmentation and lays the foundation for a series of later methods.

**Figure 20.** Brain tumor segmentations by multi-network. (**a**) Flair. (**b**) T1. (**c**) T1c. (**d**) T2. (**e**) Segmentation results. Adapted with permission from ref. [147]. Copyright 2018 Elsevier.


**Table 4.** Evaluation results of MRI brain tumor segmentation based on single CNN network.

Compared with the 2D network model, the 3D network can capture more characteristic information. The related research results are shown in Table 5. Kamnitsas et al. and Casamitjana et al. used a combination of the 3D convolution kernel and dual-pass model to improve the ability of the network to capture 3D features in MRI brain images. Among them, Kamnitsas et al. constructed a dense CNN segmentation architecture based on 3D local blocks [149]. As shown in Figure 21, the method achieved a high segmentation accuracy for brain tumor image. The proposed network combines multi-scale information by using different input sizes for each path. Casamitjana et al. employed the same input size for the two paths so that different receiving areas can focus on different information [150]. For the purpose of further improving performance, advanced network, such as Resnet and DenseNet, are gradually being introduced. Meanwhile, 3D convolution and dual-path modes are gradually being adopted in combination. Kamnitsas et al., (2017) combined dualpath and residual connection blocks to construct a 3D dual-path residual segmentation network for brain tumor segmentation [151]. It can obtain better segmentation performance while using less training data and fewer filters. Chen et al. proposed a new structure that can hierarchically segment necrotic and non-enhanced tumors and edema around tumors [152]. Therefore, 3D-CNN can accurately classify each pixel in MRI brain images.


**Table 5.** Evaluation results of MRI brain tumor segmentation based on the multi-CNN network.

The generative adversarial network GAN was proposed by Google Brain researcher Goodfellow and others. It can be considered as the most representative unsupervised generative method in recent years [155]. As shown in Figure 22, the generative confrontation network consists of generator G, discriminator D, and real data X. The generator is used to generate approximately realistic objects, so that it is difficult to distinguish true from false. The discriminator is mainly used to identify the difference between the real object and the false object produced by the generator. In the training process, the generator and the discriminator overlap each other, and continuously optimize and produce the best results in the confrontation. GAN has the inherent advantages of finding the data distribution of the model by paying attention to the potential probability density of the data, which overcomes the shortcomings of relying on large amounts of data. As the most creative unsupervised generative network, the GAN has gained wide attention in many fields such as image classification and segmentation. In recent years, it has also been effectively migrated to medical image analysis tasks, including MRI brain tumor segmentation. Table 6 is the performance of the MRI brain tumor segmentation method based on the GAN.

**Figure 21.** Brain tumor segmentations by 3D CNN. (**a**) Flair. (**b**) T2. (**c**) Manual. (**d**) Deep Medic. (**e**) Deep Medic + CRF.

**Figure 22.** The MRI brain tumor segmentation based on GAN.

At the end of 2016, the semantic segmentation model based on adversarial networks is proposed. The method treats the segmentation network as a generator to train the network [156]. Li et al. tried to use GAN for MRI brain image segmentation at the MICCAI conference [157]. They trained the generator network and the discriminator network together to reduce the difference between the synthetic label and the real label. Furthermore, the high-order loss terms are used to enhance the spatial continuity of the segmentation results. Rezaei et al. proposed a new end-to-end training conditional generation confrontation network CGAN based on the GAN network [158]. The segmentation network and the discriminator network are changed to U-Net and Markovian GAN, respectively. The two networks are simultaneously trained through backpropagation to improve the segmentation accuracy. In 2019, Rezaei et al. proposed a new confrontation network called voxel-GAN [159]. It is mainly used to alleviate the problem of data imbalance in the brain tumors segmentation, since the proportion of tumors or unhealthy areas is small. In addition, there are methods that also introduce a 3D conditional generation confrontation network to alleviate the problem of brain tumor data imbalance based on the weighted adversarial loss generated by network training.


**Table 6.** Evaluation results of MRI brain tumor segmentation based on GAN network.

#### **5. Future Research Directions**

Currently, the stability and reproducibility of biological detection are still poor due to the shortcomings of instability and variability of biologically active units. Nowadays, biological detection will inevitably be greatly developed with the rapid development of biology, informatics, materials science, and microelectronics [160–164]. It is foreseeable that electromagnetic detection will have the following characteristics in the future:

1. Integration and intelligentization

In the future, electromagnetic detection will realize information storage and memory, logical judgment, two-way communication, self-check, self-calibration, self-compensation, numerical processing, and other functions. The sensor measurement unit can become a standard unit, which can increase the portability and interchangeability of the detection system. In addition, electromagnetic detection can also be closely integrated with computers. Machine learning technology automatically analyzes the data collected by the sensor and provides analysis results more quickly and accurately, thus forming a full automation of analysis and detection.

2. Multi-sensor fusion technology

Electromagnetic detection technology can be combined with other detection technology, such as capacitive detection technology, ultrasonic detection technology, etc. The fusion data processing methods is applied to improve measurement performance. It can integrate electrical, magnetic, acoustic, thermal, force, and other multi-physical field measurements. Based on the redundant data between multiple sensors, the stability and robustness of the system can be enhanced, and the performance and stability of the system can be improved.

3. Portability

At present, flexible sensing technology has made considerable progress based on the development of soft materials, flexible electronics, and nanotechnology. The combination of electromagnetic sensors and flexible sensing technology is conducive to online and realtime monitoring. People can use portable electromagnetic sensors for daily autonomous measurement, making it possible for people to diagnose diseases and monitor food quality in daily life.

#### 4. Mechanism Research

At present, there is no reasonable explanation for the mechanism of some electromagnetic biological effects. Mechanism research should focus on reflecting biological issues from the perspective of physics. The establishment of standardized experimental models is particularly important for mechanism research. At the same time, new theories and

methods learned from the fields of mathematics, physics, chemistry, and life sciences are used to achieve a major breakthrough in the study of biological electromagnetic mechanism. Based on the electromagnetic biological effect mechanism, the application range of electromagnetic detection can be further expanded. More and more parameter information can be obtained according to the measured signal inversion. More important, the study of mechanism also contributes to the improvement of sensor sensitivity and stability.

5. Networked detection

With the rapid development of network communication technology, networked detection systems that link distributed sensor nodes and fusion centers through communication networks have gradually received widespread attention. The networked sensor system improves the flexibility and scalability of system design due to its advantages, such as less wiring, low cost, easy expansion, and maintenance. The biosensor integrates a wireless transmitter chip, which can automatically transmit the collected data to the server and display it on the mobile terminal in real time. People can check biological information anytime and anywhere to facilitate daily life. Therefore, sensor networking is also an important development direction with the application of the Internet.

#### **6. Conclusions**

Electromagnetic detection technology has developed rapidly during the 40-year development process, which has attracted an increasing number of scholars. Since bioelectromagnetic detection has different mechanisms at different frequency ranges, they are widely used in disease diagnosis and treatment, food monitoring, environmental testing, etc. In particularly, the development of MRI technology has pushed electromagnetic detection to its peak. Therefore, this review summarizes the electromagnetic biological effects and their mechanism of action. Considering the excitation frequency and intensity significantly affects the electromagnetic biological effects; this review takes the excitation frequency of the electromagnetic field as a clue to introduce the application of bioelectromagnetic detection. Moreover, electromagnetic detection in combination with ML technology has been applied in clinical diagnosis, since ML can extract the deep information of data. The research involving the application of ML technology to electromagnetic medical images is also discussed. Finally, we present an opinion on the future development trend of electromagnetic detection and make a conclusion on the reviewed biosensing techniques.

**Author Contributions:** P.H. prepared and wrote the manuscript; L.X. edited the manuscript; Y.X. reviewed the manuscript and prepared the figure/table configuration and literature research. All authors have read and agreed to the published version of the manuscript.

**Funding:** Yuedong Xie received the National Natural Science Foundation of the China Youth Fund Project (61901022) as well as Fundamental Scientific Research Business Fees of the Central Universities-Top-notch Young Talents Program (YWF-19-BJ-J-358).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **A Novel E**ffi**cient FEM Thin Shell Model for Bio-Impedance Analysis**

#### **Jiawei Tang 1, Mingyang Lu 1, Yuedong Xie 2,3 and Wuliang Yin 1,\***


Received: 4 May 2020; Accepted: 15 June 2020; Published: 17 June 2020

**Abstract:** In this paper, a novel method for accelerating eddy currents calculation on a cell model using the finite element method (FEM) is presented. Due to the tiny thickness of cell membrane, a full-mesh cell model requires a large number of mesh elements and hence intensive computation resources and long time. In this paper, an acceleration method is proposed to reduce the number of mesh elements and therefore reduce the computing time. It is based on the principle of replacing the thin cell membrane with an equivalent thicker structure. The method can reduce the number of mesh elements to 23% and the computational time to 17%, with an error of less than 1%. The method was verified using 2D and 3D finite element methods and can potentially be extended to other thin shell structures. The simulation results were validated by measurement and analytical results.

**Keywords:** finite element method; thin shell model; β dispersion; Maxwell–Wagner effect; bio-impedance spectroscopy

#### **1. Introduction**

β-dispersion analysis has been widely investigated for medical [1,2] and industrial [3,4] applications. For example, it has been proposed for the detection of tumours [5,6] and cerebral stroke [7], and it has been applied to the quality inspection of food such as meat, fruit and vegetables [3]. The mechanism of β-dispersion was first described by Schwan [8,9]. β-dispersion takes place at radio frequency, mainly from hundreds of hertz to megahertz. This dispersion is caused by Maxwell–Wagner effect, an interfacial polarisation of cell membrane blocking the ion-flow of intra and extracellular dielectrics.

Analytical solutions are well developed to calculate the β dispersion [10–12]. However, the analytical solution is only designed to analyze the models with a regular shape (i.e., spherical cell model) [13,14]. In reality, most cell shapes are anomalous. The finite element method is a methodology that discretises an entire continuous domain into discrete domains. An irregular cell shape can be simulated in a high accuracy using the finite element method. Therefore, FEM is a feasible and effective way to simulate the dielectric dispersions for models with irregular shapes [15,16].

The difficulty of FEM is that the number of meshing elements is significantly large due to the tiny thickness of the cell membrane; it takes a long time to compute with numerous meshing elements. In this paper, a novel FEM acceleration method was proposed to simulate β dispersion. The acceleration method is to reduce the number of meshing elements and computing time by replacing the full-mesh cell membrane model with an equivalent reduced-mesh model. In this paper, the FEM

acceleration method was used to simulate the spherical and oval cell model in two-dimensional and three-dimensional versions.

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

#### *2.1. Single Spherical Cell Model*

The full-mesh cell model is the single shell spherical model introduced by Asami [15]. The spherical cell model is shown in Figure 1. The parameters of the cell model are *km* = 10−<sup>7</sup> S/m, *ka* = *kc* = 1 S/m. ε*<sup>c</sup>* = ε*<sup>a</sup>* = 80, ε*<sup>m</sup>* = 5, R = 5000 nm, *dm* = 5nm. *km* is the conductivity of the cell membrane. *ka* and *kc* are the conductivity of the extracellular and intracellular fluid, respectively. ε*<sup>c</sup>* and ε*<sup>a</sup>* are the relative permittivity of the intracellular and extracellular fluid, respectively. ε*<sup>m</sup>* is the relative permittivity of the cell membrane. R is the radius of the cell. *dm* is the cell membrane thickness [10,17,18].

**Figure 1.** Spherical cell model.

The two-dimensional simulation model is a single cell model, as shown in Figure 1. The spherical cell is put in a suspension with an applied alternating electric field. The voltage of the upper electrode was set to be 10 V, and the voltage of the lower electrode was set to be 0 V [14]. The intracellular and extracellular fluid were conductive while the cell membrane blocked current flow between intracellular and extracellular fluid at low frequency. Dielectric relaxation describes a phenomenon that the permittivity decreases when the frequency of applied electric field increases at a certain range. The relaxation occurs when the electric dipoles are not able to reverse at the same rate with the frequency of the applied field. β-dispersion was observed from 10 kHz to 1 MHz. The dispersion was described by Schwan and could be explained by Maxwell–Wagner effect [8]. The effect occurs at the boundary of layered or inhomogeneous dielectrics. Applying an electric field would build up charges at the boundaries created by different dielectrics. The double-layer dielectrics system can be represented by one capacitance in parallel with one another. Applying a static electric field would charge both capacitances and creates charges at the boundary of the dielectrics. This phenomenon is also called interfacial polarisation.

The equivalent complex relative permittivity of the cell model in Figure 1 can be calculated according to Asami [12]:

$$\varepsilon\_{\rm p}^{\*} = \varepsilon\_{\rm m}^{\*} \frac{2(1 - \mathbf{v})\varepsilon\_{\rm m}^{\*} + (1 + 2\mathbf{v})\varepsilon\_{\rm c}^{\*}}{(2 + \mathbf{v})\varepsilon\_{\rm m}^{\*} + (1 - \mathbf{v})\varepsilon\_{\rm c}^{\*}} \tag{1}$$

where <sup>v</sup> = (<sup>1</sup> <sup>−</sup> *dm* R ) n , *dm* is the cell membrane thickness, R is the radius of the cell, and n stands for the factor of dimension. When n = 2, the equation is applied to the two-dimension model. When n = 3, the equation is applied to three-dimension model. ε∗ <sup>p</sup> is the complex relative permittivity of the cell model, ε∗ <sup>m</sup> is the complex relative permittivity of the cell membrane and ε<sup>∗</sup> <sup>c</sup> is the complex relative

permittivity of cytoplasm; then the complex relative permittivity of the suspension can be calculated by Wagner's equation:

$$
\varepsilon^\* = \varepsilon\_\text{a}^\* \frac{2(1 - \text{P})\varepsilon\_\text{a}^\* + (1 + 2\text{P})\varepsilon\_\text{p}^\*}{(2 + \text{P})\varepsilon\_\text{a}^\* + (1 - \text{P})\varepsilon\_\text{p}^\*} \tag{2}
$$

where P is the volume fraction, ε∗ is the complex relative permittivity of cell suspension, ε∗ <sup>a</sup> is the complex relative permittivity of the external medium.

The conductivity of the suspension can be derived:

$$
\varepsilon^\* = \frac{\sigma^\*}{\text{j}\omega\varepsilon\_0} = \frac{\sigma\_r + \text{j}\omega\varepsilon\_r\varepsilon\_0}{\text{j}\omega\varepsilon\_0} \tag{3}
$$

where σ*<sup>r</sup>* and ε*<sup>r</sup>* are the conductivity and permittivity of the suspension, respectively. σ<sup>∗</sup> is the complex conductivity of the suspension, ε<sup>0</sup> is the permittivity of vacuum, ω = 2πf where f is the frequency of the applied field and j is the imaginary unit.

According to Schwan [9], when *dm* <sup>R</sup> 1, *km ka* = *kc*, the dielectric relaxation of the cell suspension can be represented by the Debye relaxation model:

$$
\varepsilon^\* = \varepsilon\_h + \frac{\varepsilon\_I - \varepsilon\_h}{1 + \text{j}\omega\tau} + \frac{k\_l}{\text{j}\omega\varepsilon\_0} \tag{4}
$$

where ε*<sup>h</sup>* and ε*<sup>l</sup>* are the high-frequency limit and low-frequency limit of the relative permittivity, τ is the relaxation time and *kl* is the low-frequency limit of the conductivity.

#### *2.2. D FEM Simulation*

Normally, the finite element problems could be solved by two kinds of methods which are Galerkin's method and Ritz's method [19]. In this paper, Galerkin's method is used to build up simulations.

According to the Gaussian equation,

$$
\nabla \cdot \mathbf{D} = \rho \tag{5}
$$

where <sup>→</sup> D = ε → E, → <sup>D</sup> is the electric flux density, <sup>→</sup> E is the electrical field intensity, ε = εrε0, ε<sup>r</sup> is the relative permittivity of the material, ρ is the quantity of electric charge, which is zero in the cell model. So the equation could be rewritten as:

$$\nabla \cdot (\sigma + \mathbf{j}\omega \epsilon) \stackrel{\rightarrow}{\mathbf{E}} = \nabla \cdot (\sigma + \mathbf{j}\omega \epsilon) \nabla \mathbf{U} \tag{6}$$

where σ is the conductivity of the material and ε = εrε0. ε<sup>r</sup> is the relative permittivity of the material. ε<sup>0</sup> is the permittivity of vacuum. U is the electric potential which is the objective value to be solved [20]. To obtain the equivalent conductivity of the total cell suspension, we calculate the total current flowing through the cell suspension from the top electrode to the bottom electrode [15]. According to Kirchhoff's current law, the current flow through any plane, which is in parallel with the electrodes is the same with the total current flow through the cell suspension. So the current I flow through a chosen plane in parallel with electrodes can be calculated as:

$$\mathbf{I} = \sum \mathbf{j}^{(\mathbf{i})} = \sum \sigma^{(\mathbf{i})} \overleftarrow{\mathbf{E}}^{(\mathbf{i})} \tag{7}$$

where j (i) is the current density of the element on the chosen plane, σ(i) is the conductivity of the element on the chosen plane, <sup>→</sup> E (i) is the background electrical field on the chosen plane. Then the complex conductivity of the cell suspension can be calculated as:

$$
\sigma^\* = \frac{1}{\rho^\*} = \frac{\mathcal{L}}{\mathcal{S}} \ast \frac{\mathcal{I}}{\mathcal{U}} \tag{8}
$$

where σ∗ and ρ∗ are the complex conductivity and resistivity, respectively. L is the distance between the top electrode and the bottom electrode. S is the area of the electrodes.

#### *2.3. D FEM Simulation*

An induction model was built for the three-dimension FEM simulation. As shown in Figure 2, the cylinder stands for the cell suspension, and the sphericalmodel stands for the cell. Imaginary transmitter and receiver coils are put on the top of the suspension to provide an alternating magnetic field as an exciting signal and detect induced secondary magnetic field. The single-cell model in a three-dimension situation is shown in Figure 2b. The applied field is a magnetic field excited by the transmitter in Figure 2a. The radius of the sphere cell model is R and the cell membrane thickness of the sphere cell model is *dm*. All the electric parameters of the sphere cell model are set to be the same with the parameters in Figure 1. The Mxwell–Wagner effect can be calculated by calculating the eddy current distribution induced by the applied magnetic field.

**Figure 2.** (**a**) 3D suspension model, (**b**) 3D single cell model.

According to Oszkar Biro [21], the basic formulas for the three-dimension edge elements simulation of eddy current problems, Galerkin's equation is shown as following:

$$\begin{cases} \int\_{\Omega\_{\mathrm{c}}} \nabla \times \mathbf{N}\_{\mathrm{i}} \cdot \mathbf{v} \nabla \times \mathbf{A}^{(\mathrm{n})} \mathrm{d}\Omega + \int\_{\Omega\_{\mathrm{c}}} \mathbf{j} \omega \, \sigma \mathbf{N}\_{\mathrm{i}} \cdot \mathbf{A}^{(\mathrm{n})} \mathrm{d}\Omega + \int\_{\Omega\_{\mathrm{c}}} \sigma \mathbf{N}\_{\mathrm{i}} \cdot \nabla \mathbf{V}^{(\mathrm{n})} \mathrm{d}\Omega\\ = \int\_{\Omega\_{\mathrm{c}}} \nabla \times \mathbf{N}\_{\mathrm{i}} \cdot \mathbf{v}\_{0} \nabla \times \mathbf{A}\_{\mathrm{s}} \mathrm{d}\Omega \quad \mathrm{i} = 1,2,...,6 \end{cases} \tag{9}$$

$$\underset{\Omega\_{\mathbb{C}}}{\int} \mathbf{j} \boldsymbol{\omega} \boldsymbol{\sigma} \nabla \mathcal{L}\_{\mathbb{i}} \mathcal{A}^{(\mathbb{n})} \mathrm{d}\Omega + \int\_{\Omega\_{\mathbb{c}}} \boldsymbol{\sigma} \nabla \mathcal{L}\_{\mathbb{i}} \cdot \nabla \mathcal{V}^{(\mathbb{n})} \mathrm{d}\Omega = \boldsymbol{0} \quad \text{i} = 1, 2, \dots, 4 \tag{10}$$

where Li is the elemental interpolation of i th node corresponding to the nth element of it; Ni stands for the vector interpolation of i th edge relating to the nth edge element of it; A(n) stands for the excited edge vector potential of the nth element; As represents the original edge vector potential of the nth element; V(n) represents the electrical potential excited by the nth element on the receiving coil; v stands for the reluctivity of the model; v0 is the reluctivity of the air; σ is the conductivity of the model [22–24].

The eddy current density is:

$$\mathbf{J}^{(i)} = \ ^{\mathbf{a}}\mathbf{E}^{(i)}\tag{11}$$

where E(i) denotes the vector sum of the electrical field on all the edges of each tetrahedral element, σ(i) is the complex conductivity parameter of the of each tetrahedral element.

The model is shown in Figure 2. Assuming there is an equivalent model with uniform dielectric and the background suspension has exactly the same shape with the simulation model in Figure 2. Since the normal component of the electrical field relative to each cross-section surface of suspension (shown in Figure 2) is identical, then

$$\sum \overrightarrow{\mathbf{n}} \cdot \mathbf{E}^{(i)} = \sum \overrightarrow{\mathbf{n}} \cdot \overrightarrow{\mathbf{E}}^{(i)} \tag{12}$$

where, <sup>E</sup>(i) <sup>s</sup> denotes the background electrical field of the of each tetrahedral element; <sup>→</sup> n is the normal unit vector relative to the surface of target;

Since the equivalent model is uniform, the electrical background field <sup>E</sup>(i) <sup>s</sup> is always vertical to every cylindrical cross-section surface of the suspension model shown in Figure 2. Then the equivalent complex conductivity of the original suspension (as shown in Figure 2) can be derived from (8) and (9) that

$$\sigma\_{\sf s} = \frac{\sum \sf I\_{\sf s}^{(i)}}{\sum \overrightarrow{\sf n} \cdot \sf I\_{\sf G^{(i)}}} \tag{13}$$

σ(i) is the complex conductivity of each tetrahedral element, J (i) <sup>s</sup> is the eddy current density through each cylindrical cross-section of equivalent suspension model, J (i) is the eddy current density through each cylindrical cross-section of the original suspension model with arranged cells and <sup>→</sup> n is the current flow direction.

#### *2.4. Acceleration Method*

The complicated calculation progress and the long computing time are caused by a large number of meshing elements around the thin cell membrane. The number of meshing elements could be reduced by enlarging the cell membrane thickness. However, changing the cell membrane thickness affects the conductivity spectrum of the cell suspensions. The aim of the proposed acceleration method is to keep the accuracy of the simulation result while enlarging the cell membrane thickness (reducing the number of elements). This acceleration method replaces the full-mesh model with an equivalent thicker cell membrane (reduce-mesh model) with an equivalent complex conductivity, as shown in Figure 3. The M1 in Figure 3 is the full-mesh model cell membrane with normal thickness l1. The M2 in Figure 3 is the enlarged part of the reduce-mesh model cell membrane with thickness l2, R is the radius of the cell model. The electrical parameter of the cell membrane M2 is exactly the same with the intracellular fluid. The cell membrane M1 and M2 combine as one thicker equivalent cell membrane (reduce-mesh model cell membrane), and the thickness of the reduce-mesh model cell membrane is l1 + l2. The number of meshing elements is reduced due to a thicker equivalent cell membrane. In this case, the factor v = <sup>1</sup> <sup>−</sup> *dm* R <sup>n</sup> in Equation (1) can be rewritten as <sup>v</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> (l1+l2) R n to describe the Maxwell–Wagner effect in the equivalent model. According to the Maxwell–Wagner Equation (2), as the parameter P and ε∗ <sup>a</sup> of the equivalent model are fixed, the dielectric behaviour of the suspension model ε∗ only depends on the dielectric property of the cell model ε∗ p. According to Equation (1), ε∗ <sup>p</sup> is determined by ε<sup>∗</sup> <sup>m</sup> and the factor v. In order to retain the dielectric behaviour of the suspension model, the equivalent dielectric parameters of the equivalent cell membrane ε∗ <sup>m</sup> can be calculated by Equations (14)–(18). With the equivalent complex conductivity, the behaviour of Maxwell–Wagner effect, which leads to β-dispersion remains the same.

**Figure 3.** Equivalent cell membrane model.

The parameters of the cell membrane M1 and M2 in Figure 3 are:

$$
\sigma\_1^\* = \sigma\_1 + \text{j}\omega \varepsilon\_1 \varepsilon\_0 \tag{14}
$$

$$
\sigma\_2^\* = \sigma\_2 + \text{j}\omega \varepsilon\_2 \varepsilon\_0 \tag{15}
$$

where σ∗ <sup>1</sup> and σ<sup>∗</sup> <sup>2</sup> are the complex conductivity of M1 cell membrane and M2 cell membrane in Figure 3, respectively. σ<sup>1</sup> and σ<sup>2</sup> are the conductivity of cell membrane M1 and cell membrane M2 in Figure 3, respectively. ε<sup>1</sup> and ε<sup>2</sup> are the relative permittivity of the cell membrane M1 and cell membrane M2, respectively. Considering the two cell membranes M1 and M2 are connected as two electrolytes in series. The total impedance Z is:

$$\mathbf{Z} = \mathbf{Z}\_1 + \mathbf{Z}\_2 = \frac{\mathbf{l}\_1}{\sigma\_1^\ast \mathbf{S}\_1} + \frac{\mathbf{l}\_2}{\sigma\_2^\ast \mathbf{S}\_2} \tag{16}$$

where Z1 and Z2 are the impedance of the cell membrane M1 and M2, respectively. σ<sup>∗</sup> <sup>1</sup> and σ<sup>∗</sup> <sup>2</sup> are the complex conductivity of cell membrane M1 and M2 respectively. l1 and l2 are the thickness of the cell membrane M1 and the cell membrane M2 in Figure 3, respectively. S1 and S2 are the area of cell membrane M1 and the cell membrane M2, respectively. As l1 and l2 are negligible comparing with radius R, the area S1 = S2.

Then the equivalent complex conductivity can be derived as:

$$Z = \frac{\mathbf{l\_1}}{\sigma\_1^\* \mathbf{S\_1}} + \frac{\mathbf{l\_2}}{\sigma\_2^\* \mathbf{S\_2}} = \frac{\mathbf{l\_1} \sigma\_2^\* + \mathbf{l\_2} \sigma\_1^\*}{\sigma\_1^\* \sigma\_2^\*} \frac{\mathbf{1}}{\mathbf{S\_1}} = \frac{(\mathbf{l\_1} + \mathbf{l\_2})}{\sigma\_\mathbf{e}^\*} \frac{\mathbf{1}}{\mathbf{S\_1}} \tag{17}$$

$$
\sigma\_{\mathfrak{e}}^{\*} = \frac{(\mathfrak{l}\_1 + \mathfrak{l}\_2)\sigma\_1^\* \sigma\_2^\*}{\mathfrak{l}\_1 \sigma\_2^\* + \mathfrak{l}\_2 \sigma\_1^\*} \tag{18}
$$

σ∗ <sup>e</sup> is the complex conductivity of the equivalent cell membrane (reduce-mesh model cell membrane).

The replacement of the thin cell membrane with an equivalent thicker structure appears to have reflected an equivalent Maxwell–Wagner effect from the computation results. Theoretically, less contrasting EM properties between the boundary of two materials cause less Maxwell–Wagner effect (i.e., less charge building-up). At the same time, thicker materials between the two charged surfaces would mean the equivalent capacitance is smaller. These two effects would mean the electric field E = Q/C would be kept to a constant; this is possibly the reason why we can use the replacement to carry out the calculation of the dispersion caused by the MW effect. Equations (17) and (18) considered the complex conductivity, and therefore, the physics would be masked to a certain degree, but the principle can be explained as above.

#### *2.5. Meshing Details*

The suspension model shown in Figure 1 was discretised into three regions, cell membrane and intra and extracellular fluid, for meshing. The commercial software Comsol 5.0 was used for meshing. The convergent criterion was accelerated BICGS with an optimal precondition. The initial precondition was based on optimised initial guess: the final solution from the previous frequency [25–30]. As shown in Table 1, the minimum and maximum element sizes are set to be the same to maintain the meshing accuracy. The maximum element growth rate controls the maximum changing rate of the dimension of the element among the adjacent subdomains. The narrow factor controls the mesh density on the narrow region. Increasing the cell membrane thickness reduces the narrow regions. Therefore, the meshing density is larger in the cell model with smaller cell membrane thickness. Table 2 shows the number of meshing element in each specific region, including total suspension, cell membrane, intracellular fluid and extracellular fluid. In addition, a specific region surrounding the cell membrane was created to obtain the number of meshing elements around the cell membrane. It is obvious from Table 2 that the most meshing elements concentrate at the region around the cell membrane. The narrow region is reduced when the cell membrane thickness is enlarged, thus, the number of meshing elements is decreased when increasing cell membrane thickness.



**Table 2.** The number of meshing element in different regions.


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

#### *3.1. D Single Spherical Cell Model*

The full-mesh model cell membrane thickness is 5 nm. There are two simulations with equivalent thicker cell membrane (reduce-mesh model) thickness of 10 nm and 20 nm, respectively. The ratio between the full-mesh model cell membrane (cell membrane M\_1 in Figure 3) thickness and the enlarged part of reduce-model cell membrane (cell membrane M\_2 in Figure 3) thickness is 1:1 and 1:3, respectively. The simulation results are the progress of β-dispersion which reflects dispersions on relative permittivity and conductivity, as shown in Figures 4 and 5. The main characteristics of β-dispersion are the dispersion frequency range and the magnitude of the relative permittivity and conductivity [27,30,31]. As shown in Figure 4, the dispersion frequency is ranging from 100 kHz to 10 MHz. The magnitude of relative permittivity at low frequency and high frequency are 960 and 80, respectively. The magnitude of conductivity at low frequency and high frequencies were 0.925 S/m and 1.04 S/m, respectively. This work was aimed at accelerating the computing progress with high accuracy, 15 sample points are selected in Figure 4 to describe the main feature of an entire β -dispersion and verify the accuracy of the acceleration method. The results of acceleration models are validated by the full-mesh model to verify the accuracy of the acceleration method. The errors between the acceleration model and full-mesh model at each point in Figures 4 and 5 have been calculated. The calculated error and computing time for acceleration model and full-mesh model are shown in Table 3.

**Figure 4.** Permittivity of 2D single spherical cell model.

**Figure 5.** Of 2D single spherical cell model.


1:1 35,272 0.07% 28 1:3 17,342 0.2% 13

**Table 3.** Results of the 2D spherical cell model.

As shown Table 3, the number of elements of the full-mesh model can be reduced to 24% when the equivalent cell membrane (reduce-mesh model) is four times as the thickness of the full-mesh model cell membrane (M1 cell membrane to M2 cell membrane thickness ratio is 1:3). The computing time was reduced from 73 min to 13 min, with only a 0.2% error on the simulation result. The computing time of the reduce-mesh model is 18% of the full-mesh model. The spectroscopy of conductivity and relative permittivity is shown in Figures 4 and 5, the results of the reduce-model and full-mesh model show the same magnitude over the same frequency range with an error less than 0.2%. The β dispersion of reduce-model and full-mesh model both starts at 100 kHz and both ends at 10 MHz. This shows that the reduce-mesh model can significantly accelerate the computing progress with only a tiny error on the result in the 2D-FEM spherical model simulation.

#### *3.2. D Cell Deformation Model*

This simulation is to verify that the acceleration method not only works on spherical cell model but also works on the deformation model (oval model).

The deformation model is an oval model in order to simulation the cell deformation [28]. The electrical parameters of the oval model are the same as that of the custom spherical model. The model is shown in Figure 6, where a and b are the length of the semi-major and semi-minor axis, respectively. The length of the semi-major axis is set to be a = 12 μm and the length of the semi-minor axis is set to be b = 2 μm, *ka* and *kc* are the conductivity of the extracellular and intracellular fluid, respectively, ε*<sup>c</sup>* and ε*<sup>a</sup>* are the relative permittivity of the intracellular and extracellular fluid respectively, ε*<sup>m</sup>* is the relative permittivity of cell membrane, *dm* is the cell membrane thickness. The full-mesh model cell membrane thickness is still 5 nm and the electrical properties are calculated using equation 18. The simulation result of β dispersion is shown in Figure 7, Figure 8 and the parameters of the computing time and error are shown in Table 4.

**Figure 6.** Deformation model.

**Figure 7.** The relative permittivity of the 2D single deformation cell model.

**Figure 8.** The conductivity of 2D single deformation cell model.


**Table 4.** Results of 2D Deformation Cell Model.

The number of elements of full-mesh model can be reduced to 24.1% as shown in Table 4 and the computing time is reduced from 2h 24mins to 25mins. The simulation error between reduce-mesh model and full-mesh model is only 0.11%. As shown in Figures 7 and 8, the results of the full-mesh model and reduce-mesh model exhibits β dispersion with tiny error in the same frequency range. The results show that the acceleration method is feasible to simulate not only the spherical model but also deformation model (oval model). The computing time is significantly reduced with acceptable error.

#### *3.3. D Spherical Cell Model*

The simulation model used for 3D FEM simulation is shown in Figure 2. 3-D single sphere cell model is used. The radius of the cell model is set to be 5 μm, and the cell membrane thickness is 5 nm. The top view of eddy current distributions at lower frequency 1 kHz and higher frequency 10 MHz are shown in Figure 9. The single-cell model blocks eddy currents at a lower frequency and becomes conductive at a higher frequency which meets the expectations of maxwell–wagner effect and other published works [9,11,15,29].

**Figure 9.** (**a**) Eddy current flow at low frequency 1 kHz. (**b**) Eddy current flow at high frequency 10 MHz. (**c**) Eddy current density at low frequency 1 kHz. (**d**) Eddy current density at high frequency 10 MHz.

The full-mesh model cell membrane thickness is 5 nm. There are two reduce-models with equivalent thicker cell membrane (reduce-mesh model) thickness of 10 nm and 20 nm, respectively. The ratio between the full-mesh model cell membrane thickness and the reduce-mesh cell membrane thickness is 1:1 and 1:3, respectively. The relative permittivity and conductivity result are shown

in Figures 10 and 11. β dispersion can be observed from Figures 10 and 11. The β dispersion of the full-mesh model and reduce-mesh model has an error of 2% on the magnitude. The frequency range of the β dispersion of full-mesh model and reduce-mesh model are the same, ranging from 100 kHz to 10 MHz. The results, including the number of elements, the computing time and error of full-mesh model and reduce-mesh model, are shown in Table 5. The number of elements of the full-mesh model can be reduced to 28%, as shown in Table 5, and the computing is reduced from 4 h 37 min to 71 min. The simulation error between reduce-mesh model and full-mesh model is 2%. The error is larger than the two-dimensional simulation result, but it is acceptable for β-dispersion simulations. The results show that the acceleration method is feasible to simulate a three-dimensional spherical model. The computing time is significantly reduced with acceptable error.

**Figure 10.** The relative permittivity of the 3D single spherical cell model.

**Figure 11.** The conductivity of 3D single spherical cell model.


#### *3.4. Experimental Validation by Contact-Electrode Measurement*

The four-electrode measurement system is used to measure the impedance spectroscopy of biological samples. The electrode is contacted directly to the samples to measure the impedance over the samples using an impedance analyser. The measurement system work in the following way: the electrodes generate an input signal which applies an alternating electric field and current in the samples. By measuring the current and voltage across the samples, the conductivity and permittivity of the samples can be calculated. The impedance analyser used in this work is Solatron 1260, which is efficient and accurate over the frequency ranging from hundreds of hertz to 32 MHz. The effective area of the measurement electrodes is fixed to <sup>S</sup> = 2 cm <sup>∗</sup> 3 cm = 6 cm2, and the distance between the

electric field measuring electrodes is fixed to L = 6.6 cm. The conductivity can then be calculated by σ = 1/ρ = L/RS. R is the real part of the measured impedance, ρ is the resistivity. The measurement sample is fresh potato (obtained from a local market). This measurement is to validate the FEM simulation result. All FEM simulation parameters are set to be the realistic value obtained from the measurement.

In Figure 12, the FEM and measurement results show similar β -dispersion with the same dispersion frequency range. There is some error on the magnitude of the conductivity between the measurement result and the simulation result. However, the error is acceptable for the simulation approach of the measurement result. The conductivity of measurement and simulation result is low at a lower frequency, and the dispersion starts at around 50 kHz. The dispersion ends at around 2 MHz and the conductivity of measurement and simulation result is increased to 0.15 S/m. The spectroscopy curve is flat over the higher frequency range. The measurement result validated the full-mesh and reduce-mesh FEM simulation result and agreed with other published works [9,13,15,21,25,32].

**Figure 12.** Impedance spectroscopy of FEM and measurement result.

#### *3.5. Verification by Analytical Result*

The analytical result can be calculated according to Maxwell–Wagner Equations (1) and (2). The analytical result is compared with the 3D finite element simulation result. So the factor v in equation (1) is defined as v = <sup>1</sup> <sup>−</sup> *dm* R 3 . According to Schwan [9], the impedance spectroscopy of the cell model can be described by Debye relaxation based on Equation (4). In Figures 13 and 14, both analytical and FEM results show the same magnitude and frequency range of beta dispersion with an error less than 0.1%. The volume fraction of the cell is set to be the same as P = 3.5%. The magnitude and frequency range of β-dispersion on the full-mesh model and reduce-mesh model agree with the analytical result. The result validates that the three-dimension FEM simulation on cell suspension is accurate and the improved acceleration method is also validated. The proposed acceleration method can be applied to further FEM simulations on irregular shape cell models and thin shell models.

**Figure 13.** Calculated conductivity of the analytical solution and FEM simulation.

**Figure 14.** Calculated relative permittivity of the analytical solution and FEM simulation.

#### **4. Conclusions**

This paper proposed a method to accelerate FEM calculation with bio-cell models. The idea is to replace the thin cell membrane (full-mesh model) with an equivalent thicker cell membrane (reduce-mesh model). Then the number of meshing elements of the full-mesh model is reduced, and thus the computing time is reduced.

According to the simulation results, the reduce-mesh model can be used in both 2D and 3D FEM simulations. The amount of computing time is significantly reduced with an error no more than 0.5% in 2D simulation. The error on 3D simulation result is no more than 2%. All simulation results are validated by measurements and analytical results.

The proposed acceleration method is validated to be fast and accurate on cell models and the acceleration method has potential for all other thin shell FEM models.

**Author Contributions:** J.T. and W.Y. designed, planned and analysed the major acceleration method and prepared and reviewed the manuscript. M.L. helped in building the 3D FEM simulation server. Y.X. helped in building the experimental set-up for validation and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

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

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Development of an Optical Method for the Evaluation of Whole Blood Coagulation**

**Marinos Louka and Efstathios Kaliviotis \***

Department of Mechanical Engineering and Materials Science and Engineering, Cyprus University of Technology, 3036 Limassol, Cyprus; ma.louka@edu.cut.ac.cy

**\*** Correspondence: e.kaliviotis@cut.ac.cy

**Abstract:** Blood coagulation is a defense mechanism, which is activated in case of blood loss, due to vessel damage, or other injury. Pathological cases arise from malfunctions of the blood coagulation mechanism, and rapid growth of clots results in partially or even fully blocked blood vessel. The aim of this work is to characterize blood coagulation, by analyzing the time-dependent structural properties of whole blood, using an inexpensive design and robust processing approaches. The methods used in this work include brightfield microscopy and image processing techniques, applied on finger-prick blood samples. The blood samples were produced and directly utilized in custommade glass microchannels. Color images were captured via a microscopy-camera setup for a period of 35 min, utilizing three different magnifications. Statistical information was extracted directly from the color components and the binary conversions of the images. The main advantage in the current work lies on a Boolean classification approach utilized on the binary data, which enabled to identify the interchange between specific structural elements of blood, namely the red blood cells, the plasma and the clotted regions, as a result of the clotting process. Coagulation indices produced included a bulk coagulation index, a plasma-reduction based index and a clot formation index. The results produced with the inexpensive design and the low computational complexity in the current approach, show good agreement with the literature, and a great potential for a robust characterization of blood coagulation.

**Keywords:** blood coagulation; image sensing; image classification

#### **1. Introduction**

Blood flow allows the distribution of oxygen and nutrients to the organs and tissues. The importance of blood circulation is obvious in cases where health related problems affect people's daily life. For instance, the conditions of hemophilia and hyper-coagulation have a significant impact on patients mortality rate. Hemophilia can cause major internal and external uncontrollable bleeding from minor injuries, due to the inability of proper clot formation. On the other hand, hyper-coagulation is associated with abnormal and uncontrollable formation of blood clots in the absence of injuries, or the inability of fibrinolysis to control clot growth and clot dissolution. Both conditions cause pathological problems and patients in many cases need daily monitoring in order to regulate their medication [1].

Hemostasis is a physiological process in which blood loss due to vessel injury is minimized by clotting, so that blood flow is maintained normal in the vessels. The hemostasis process involves the mechanisms of clot formation at the site of the injury, the prevention of uncontrollable clot growth, and the lysis (dissolution) of the clot. This process is characterized by a series of reactions which involve platelets, red blood cells (RBCs), coagulation factors and anticoagulants in a three-stage procedure. In the first stage of the process, named the primary hemostasis, an initial soluble white clot is created, consisting mainly of activated platelets. Activated platelets, at the site of the endothelial injury, bind with Von Willebrand factor through the GPIb-receptor [2]. Von Willebrand is a glycoprotein

**Citation:** Louka, M.; Kaliviotis, E. Development of an Optical Method for the Evaluation of Whole Blood Coagulation. *Biosensors* **2021**, *11*, 113. https://doi.org/10.3390/bios11040113

Received: 2 March 2021 Accepted: 6 April 2021 Published: 9 April 2021

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

**Copyright:** © 2021 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/).

produced in the endothelial cells and it is critical to the initial clotting as it helps the aggregation of the platelets at the site of the injury. In the next stage, the secondary hemostasis (also known as blood coagulation), the inactive coagulation factors, that exist in the plasma under physiological conditions, are activated. In a complex process the coagulation factors promote the transformation/conversion of the soluble white clot and fibrinogen into an insoluble, stable clot and fibrin mesh, that plugs the damaged vessel [3,4].

Blood coagulation is described by the coagulation Cascade (or Waterfall model), which includes the intrinsic and extrinsic pathways, both leading to the common pathway. The final stage of hemostasis regards the deactivation of the coagulation factors and the dissolution of the clot to prevent oversize and unnecessary clotting at the site of the injury, which can affect the blood flow. This stage is known as fibrinolysis and consist of natural inhibitors that terminate the growth of clots and break them down. All three of the aforementioned stages require proteins, coagulation factors and natural anticoagulants to initiate and complete hemostasis [5]. Imbalance between coagulation and anticoagulation proteins as well as deficiencies in any of them, can cause bleeding disorders like hemophilia A, hemophilia B and Von Willenbrand disease, or arterial/venous thrombosis [6–8].

When blood is removed from the body and placed in a glass test tube, red blood cells (RBCs) are able to aggregate due to low shear stress, blood clots fairly quickly, and the blood coagulation process can be assessed immediately. Red blood cell aggregation (RBCA) causes the formation of stacks of RBCs in coin-like structures (termed rouleaux), the formation of larger aggregates from rouleaux combination, and even a 3-dimensional network of aggregates [9]. The result of the contraction of RBCs to aggregates is the appearance of larger plasma regions in the blood sample. Calcium ions are required for the clotting process, whereas acidic substances like EDTA (ethylenediaminetetraacetic acid) or citrate can prevent the formation of clots by binding calcium. EDTA and citrate are used widely in the Hematology laboratory as anticoagulants. Clotting can be initiated in vitro at a later time by adding back an excess of calcium ions [10].

A variety of laboratory equipment and point-of-care (POC) devices have been developed through the years for blood coagulation assessment, exploiting a range of blood properties [11]. In those tests, platelet count, functionality, adhesion and aggregation are of primary importance. Platelet testing provides valuable information concerning patients with bleeding disorders; it allows to monitor patient's medication treatment, assess in-surgery hemostasis and guide transfusion medicine with manual or automated instruments [12]. The cost of hand-held POC devices and their disposable strips is not negligible, whereas time and complexity might be an issue for the techniques (such as thromboelastography) that exploit the mechanical properties of blood coagulation.

The extrinsic and intrinsic pathways can be used for in vitro laboratory test screening, extracting parameters like the prothrombin time (PT) and activated partial thromboplastin time (aPTT), respectively [13]. PT and aPPT indicate the time that it takes for clots to form in plasma when calcium and tissue thromboplastin are present, and their magnitudes depend on device specific parameters. The extrinsic pathway test, aPTT, is sensitive to factors, I, II, V, VII, X, whereas the intrinsic pathway test, PT, is sensitive to factors XII, XI, IX, VIII, and further, sensitive to deficiencies of prekallikrein and high molecular weight kininogen (HMWK). These indices can be derived with a variety of instruments, with each of the instrument exploiting different blood coagulation properties [14]. Also, an important widely used index for blood coagulation is the International Normalize Ratio (*INR*) which is derived from the PT test and it is used for effective monitoring patients under oral anticoagulation therapy [15]. *INR* is calculated by the formula:

$$INR = \left(\frac{PT\_{patient}}{MNP T\_{patient}}\right)^{ISI} \tag{1}$$

where *MNPTpatient* is the mean normal *PT* of a control group and *ISI* is the International Sensitivity Index, which is used to adjust sensitivity difference of a variety of reagents used by instruments.

Electrochemical techniques for analysing coagulation are widely adopted in Hematology laboratories and POC devices. These techniques are based on resistance/capacitance changes of whole blood as clot related macromolecules react to the added agents in the sample. During the coagulation process the dominant impedance is moved from the plasma impedance to the aggregated red blood cells, and finally to the fibrin network impedance [16]. As blood coagulation escalates in time, the total impedance increases and can be accurately measured and provide the *INR* index, using the formula in Equation (1). This method leads to smaller and patient friendly POC devices as it is based on integrated electronics platforms and disposable test strips, with examples like CoaguChek(*r*) XS (Roche Diagnostics, Indianapolis, IN, USA), microINR (iLine(*r*) Microsystems S.L., San Sebastian, Spain) and i-STAT (Abbott Point of Care Inc., Abbott Park, IL, USA). However as mentioned, the cost of such devices and their disposable strips is not negligible [17,18].

Another approach to assess further details of the blood coagulation process is the Thromboelastography method (*TEG*(*r*), Haemonetics, Braintree, MA, USA), which examines the mechanical characteristics of blood in terms of viscoelasticity and kinetics as blood clots form. TEG was first introduced by Hertert in 1948 [19], as a global essay that characterizes in real time the creation of clots in whole blood using a container for the sample and a constant forced rotating rod immersed into the sample. The rod is connected to a transducer, which measures the torque generated by the developed clot. Rotational thromboelastometry (*ROTEM*(*r*) TEM International, Munich, Germany) is considered the evolution of TEG, following the same fundamental principles. Both exploit the viscoelasticity on whole blood as coagulation and lysis occurs with the use of reagents, producing essentially similar indices based on characteristics of the generated curves. The basic parameters in TEG and ROTEM need 30 to 60 min and are calculated in terms of the time-evolved curve amplitude (in mm) which indicates the clot firmness. Coagulation indices are derived from linear combination of a variety of the parameters like the maximum reached amplitude, kinetic time, clotting initiation time and the angle between the initiation of coagulation tangent to the developed curve [20]. It is also suggested that Thromboelastography can potentially help to guide patients needs for transfusions in surgery and trauma induced coagulopathies, through specialized algorithms embedded in the devices [21,22].

Coagulation in blood causes intense structural changes that can be detected by image analysis techniques. Machine vision approaches offer a range of methods and techniques for image processing [23], and are used for the development of Lap-on-Chip devices [24]. These techniques can be used for real-time monitoring and characterization, in biomedical research and clinical diagnostics [24]. Methods such as image segmentation allow the identification of structures with specific characteristics [25,26]. To accommodate the blood sample to be processed, simple microfluidic configurations are often used, offering several advantages [27–30]. Other techniques, such piezocoagulography, have been also utilized for blood coagulation, showing also the possibility of studying the initial stage of fibrin formation in the clotting process [31].

In recent works presented in the literature, the optical properties of coagulating blood are exploited, by measuring the light propagation and/or scattering through a sample (please see [11] for a review on the subject). The produced time-varying curves provide valuable information about the blood coagulation process. This technique, termed photometry, requires the acquisition of sequential light patterns, usually images, utilising a light source and a camera. Photometric results can be achieved with monochromatic or polychromatic light sources, and experimental setups covering the infrared to the xray spectrum, utilizing appropriate cameras. It has been shown that light sources with a wavelength around 690 nm (near-infrared window) provide very good results, due to the minimal light absorption by hemoglobin and water [32]. During coagulation, the light patterns acquired by the camera change dynamically as clots develop in the sample.

Furthermore, light transmission through the sample is blocked by the newly formed clots and the averaged detected light is decrease due to the interaction between the platelets, the RBCs and the procoagulant plasma proteins [33,34]. The aforementioned phenomena can be characterized, after digitally process the captured images. Several studies in the literature utilize laser speckle analysis between consecutive images, and/or between a reference image (usually the initial one), in combination with cross- or auto-correlation techniques, which can be computationally demanding. The results are presented as time-varying curves with the potential of dividing them into different coagulation phases [35–37]. Indices extracted from those curves can describe blood coagulation and typically are validated by comparison with commercially available coagulometers [38,39]. Techniques, such us the aforementioned, utilize color spectra to assess blood coagulation and extract useful indices, by using a specific light wavelength. The produced information can be visualized by using pseudo-colored maps and as a result, they may lack a direct representation of the whole blood clotting process. Moreover, such techniques require specialized components and set-ups, such as laser diodes, optical components (splitters, expanders, polarizers, etc.).

In this work, an optical assessment technique for blood coagulation is presented, using whole blood from a healthy volunteer, with conventional bright-field microscopy. The tests were conducted using capillary whole blood obtained by finger prick without added reagents. Image processing algorithms have been developed to analyze a series of images acquired from a camera/microscope setup, in gray, color and binary levels, using three magnification lenses and custom-made microchannels. The produced curves display three distinct phases, the RBCA, the clot development and the completion of whole blood coagulation. A classification algorithm separates the images, by detecting the areas in which clots are developed. Results extracted from the four color channels (red, green, blue and gray scale) show good sensitivity and conclude that this method can be used to assess whole blood coagulation with low computational demands and inexpensive equipment. The current work is set in order to demonstrate the validity and the potential of the proposed approach.

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

#### *2.1. Sample Preparation*

The blood donation protocol for this work was approved by the Cyprus National Bioethics Committee (permission ref.: EEBK/EP/2019/19). A drop of blood was collected through a finger prick from a healthy volunteer (following 10 h of fasting) after the application of intermittent pressure to a finger in order to increase the blood flow. Then the area was cleaned with alcohol and a sterilized lancet was used to puncture the side of the fingertip. The first drop of blood was wiped with a sterilize cotton ball and the second drop of blood was placed at the entrance of a custom-made microchannel. The blood drop was then drawn by capillary forces in the observation area. A total of three different experiments were made and presented in this work, at three different microscope objectives. Each set of tests was completed in approximately two hours.

#### *2.2. Experimental Setup*

The microchannels were composed by two microscope glass slides (75 mm × 25 mm × 1 mm, for length, width and height, respectively) and double-faced adhesive tapes (Tesa(*r*) ) with dimensions 12.7 mm width and approximately 100 μm thickness. The tape was adhered in pairs and at a small angle between them, on one slide (Figure 1), creating an open space. The second glass slide was placed above the first one and adhered on the tapes, creating the microchannel and effectively the testing site. Two testing sites were created on the same glass slide. The glass slides overlapped for approximately 5 mm length, with entrance and exit widths of approximately 5 and 3 mm respectively, resulting in a small microchannel volume of approximately 2 mm3. The angle of the tapes and the relatively small distance from the entrance to the exit resulted in fast filling, and fast deceleration and immobilization of the blood in the channel. After the placement of the sample in the

channel, sunflower oil was applied at the entrance of the microchannel to isolate the sample from air, and restrict undesirable effects, such as air bubble formation and blood movement. The viewing window was at the center of the geometry at distances of approximately 2 mm from the entrance and exit. This region of interest in the center of the geometry, appeared to have uniform characteristics, not influenced by the boundaries of the channel.

**Figure 1.** (**a**) The design of the microchannels. (**b**) Schematic of the experimental setup, which is based on forward propagation of light into to the sample, through the objective and into a CMOS camera for the acquisition of data. (**c**) Representative images showing the initial, an intermediate and a final state of coagulated blood at the 10×, 20× and 50× magnification objectives. Scale bars equal to 200, 100 and 40 μm for the 10× 20× and 50× objectives, respectively.

For video/image acquisition a JVC TK-C1380 camera was used, integrated on an Olympus BX51 microscope (Figure 1). The microscope light illumination (100 W longlife halogen bulb) was controlled by an intensity knob which was calibrated and scaled utilising an adhered protractor around the knob. Tests were conducted to validate the linearity between the detected light intensity and the protractor's degrees, with a deviation of around ±5 degrees. So, the protractor was utilized to set the light intensity appropriately, with a fine approximation in degrees. Before sampling, the camera focus was optimized manually by adjusting the microscope stage, and by observing the camera image to avoid opaqueness and loss of detail.

The camera was connected to a computer in which the video was saved. Three experiments were conducted with different magnifications lenses of 10×, 20× and 50× of numerical apertures of 0.25, 0.40 and 0.50, respectively, using blood sample from the same volunteer. Each acquired video has a 35 min duration using the free software VirtualDub at 5 frames per second (fps), enough to detect the blood coagulation process. Representative images from the three magnifications and at different times are shown in Figure 1. The resolution of the camera/microscopy setup was 1.0 μm/pixel, 0.55 μm /pixel and 0.22 μm/pixel ±10% for the 10×, 20× and 50× respectively, the observation image area was 576 × 768 pixels, and the processed image area were 500 × 500 pixels.

#### *2.3. Whole Blood Sample Coagulation Evaluation*

Blood coagulation was evaluated by digitally processing the acquired videos utilizing the developed algorithms (The MathWorks Inc., Matlab-Natick, MA, USA). The video acquisition rate was initially set to 5 fps, however various reduced sampling rates were used in the processing stage in order to minimize storage needs and optimize the performance in terms of overall processing time. Apart from various sampling rates, different image sizes were used for performance optimization in the processing stage, while keeping the main body of the algorithm unchanged. After the evaluation of the produced time-depended indices, the minimum frame rate, preserving the details of the process, was found to be at 4 s per frame (0.25 fps), with total number of images per experiment to be 525 colored images.

The processing stage can be divided into three main sections, with each one calculating various values and parameters of interest, from the time-varying image series. First, each image is separated into their red, green, blue and gray-scale (RGBGs) color channels. For each of the four color channels the mean intensity (MI) value is calculated for each image and for all sequential images from the various tests. The MI corresponds to the average light travel through the sample at each specific moment. MI curves for all color channels were smoothed in a moving average manner in order to extract more meaningful information at the initial parts of the process (up to approximately 10 s), where the data appear noisier. The data were further normalized with the maximum values of the time series (and denoted as *MI* ∗), so that their time response of the various cases can be directly comparable.

In the second algorithmic section, the images were converted into binary data (i.e., black and white pixels), using a unique threshold for each color channel, and for each experiment. In the initial stages of the process, when no clots are present, the binary images have pixels with values one (white pixel) and zero (black pixel), representing the plasma and RBCs areas, respectively. This is due to the fact that the transmission of light through the plasma is higher than the transmission of light through the RBCs and clots (see for example the sample image at *t*0, and at 50× magnification in Figure 1). So, as RBCs block more light than the plasma, their gray scale values is lower to that of the plasma pixels. In the transformation of a gray scale image into binary the threshold value is the most important parameter. The threshold is a gray scale value that separates all the pixels into two categories, those above and those below the threshold value. The RGBG thresholds were based on the image characteristics in the peak value of their *MI* ∗ curve. That threshold value was selected for two reasons. First, at that point in time the phenomenon of RBCA is mostly complete, as it takes around three-four minutes to fully develop, and blood coagulation is the main function that affects the structure of the blood sample. Until this point an increase in the averaged light sensing was observed. Thus, this time point can be referred to as the clotting start (CS) time after which blood coagulation can be studied. Secondly, by deriving the binarization threshold value from the CS correspond image (resulting in a high threshold value), ensures the efficient detection of the growing clots in the subsequent analysis, as the distribution of pixel intensity will lean towards lower values due to blood coagulation. As the new structures, the clots, are forming they eventually assigned zero pixel values, due to their light-blocking (absorbing) nature. Thus, the area occupied by clots increase in size, by the action of coagulation factors that blood plasma contains, and the corresponding plasma area is decreased. As a result, the percentage of the white pixels correspond to the plasma area is decreased in time. Therefore, the black pixels with zero values after processing, represent RBCs (and their aggregates) and the newly formed clots. A significant advantage of the present algorithm is that RBCs and newly formed clot areas can be distinguished and properly analyzed.

In the third section of the algorithm a Boolean algebra calculation loop was developed to extract indices from the binary image series with reference to the image at the CS time, in a pixel by pixel manner, similarly to methods used in studying elastic optical networks [40]. Karnaugh mapping (KM) was used to simplify the expressions using logical operators. As mentioned above, the clot detection was classified as the change of pixels with values one to zero (i.e., from plasma pixel to clots pixel). As each binary image only consists of two

values, ones and zeros, there are four possible combinations when comparing two images, each corresponding to a specific condition/stage regarding the blood coagulation process. The first combination is the zero to zero pixel value that regards the detection of RBCs or RBCA region that remains unchanged. The second regards changes from a zero value to one, denoting a new plasma area (NPA) which can appear due to the RBC aggregation (RBCA), in which RBCs attract each other leaving an empty from cells space (plasma). This procedure however is considered complete after the initial 4–5 min as mentioned above. In a third combination, a change from one to zero corresponds to the clot formation (CF) in areas where previously were occupied by plasma. In the fourth combination, values compared from one to one represent the plasma area (PA) still unoccupied by blood clots. As in the above section of the algorithm, the change of a white pixel to a black pixel is to be expected and indicates the activity of the coagulation process in the sample. For the four pixel comparisons mentioned above, the image corresponding to CS time is used as the reference image. Using the results of the pixel comparison the algorithm can separate areas as RBCs, the newly formed clots (CF), the initial plasma area (PA) at CS, and the new plasma areas (NPA).

Table 1 shows the truth table of the aforementioned combinations, and gives information about the four outputs, which can be conceived as four new images. Each new image has binary pixel values according to the output.

**Table 1.** The four combinations of the pixel-wise comparison between the image I at reference time τ (denoted as *I*(*τ*)*i*,*<sup>j</sup>* , and the subsequent images *I*(*τ* + *t*)*i*,*<sup>j</sup>* .


*I*(*τ*)*i*,*<sup>j</sup>* denotes the value of the pixel in the *i*th row, *j*th column of the reference image and *I*(*τ* + *t*)*i*,*<sup>j</sup>* the corresponding pixel of the sequential image at a time *τ* + *t*, where τ is the CS time. The outputs RBCs, NPA, CF and PA are new binary data (images) at a given time, after CS, and they have the same size as the compared images, representing specific structural conditions in blood. The information from the Table 1 and the utilization of KM are used to derive indices describing, and quantifying, each of the aforementioned conditions in every image. The indices are developed according to Equations (2)–(5), with N noting the number of pixels of the image (total number of columns *Ni* and total number of rows *Nj* are equal to N, *I* denoting the complement pixel value of I at that position (logical NOT operation), and with the symbol ˆ denoting the logical AND operation. For each combination and for all images, the total number of ones is divided by the total number of pixels in the image, resulting in the expression of the indices as fractional quantities (denoted with the \* superscript).

$$RBCs^\* = \frac{1}{N} \sum\_{i=1}^{N} \sum\_{j=1}^{N} \overline{I}(\tau)\_{i,j} \overline{I}(\tau + t)\_{i,j} \tag{2}$$

$$NPA\ ^\ast = \frac{1}{N} \sum\_{i=1}^{N} \sum\_{j=1}^{N} \overline{I}(\tau)\_{i,\hat{j}} \hat{I}(\tau + t)\_{i,j} \tag{3}$$

$$CF^\* = \frac{1}{N} \sum\_{i=1}^{N} \sum\_{j=1}^{N} I(\pi)\_{i,j} \tilde{I}(\pi + t)\_{i,j} \tag{4}$$

$$PA^\* = \frac{1}{N} \sum\_{i=1}^{N} \sum\_{j=1}^{N} I(\tau)\_{i,j} \hat{I}(\tau + t)\_{i,j} \tag{5}$$

The addition of Equations (3) and (5) represents the total apparent blood plasma at any given time, after the CS time. Moreover, one can see that the pairs Equations (2) and (5), and Equations (3) and (4) are complementary due to their definition from Table 1. Therefore, it is of interest to further analyse the relation between *RBCs* ∗ and *CF* ∗. Figure 2a shows sections from a reference and a compared image, and the four output matrices. Figure 2b illustrates the behavior of *RBCs* ∗ and *CF* ∗. As mentioned earlier, the RBCA process was mostly complete at CS, as normally reaches a steady state after around 120 s. However, this time period could be extended due to various reasons. For example, small cell movements inside the microchannel could contribute to a prolonged RBCA time. After CS, a small activity of RBCA was observed for a small period of time, and that resulted in a small increase in *NPA* ∗ and a small decrease in *CF* ∗. The apparent increase of *RBCs* ∗ after it has reached a minimum (which denotes the completion of RBCA) seems artificial, since it is caused by the re-appearance of a zero value pixel area that should correspond to *CF* ∗, at the same location where previous *NPA* ∗ was detected. In terms of pixel values, a zero in the reference image had initially compared with one (corresponding pixel in the subsequent image) but at a later time it was compared with a zero again. However, this condition most probably corresponds to the detection of clots in the region of NPA and not to a new RBCs. In a correction process, the increase that occur from the re-appearance of *RBCs* ∗ is subtracted, and the subsequent increase is added to *CF* ∗ as Figure 2b shows.

In order to distinguish the most efficient image resolution and color channel for the image processing, a sensitivity index (*SI*) was developed (Equation (6)) to compare the resulted curves. The *SI* was applied only to the plasma reduction (*PR*) and the CF curves, as the *MI* ∗ curves are normalized using their own limits. *PR* is defined as the ratio of the total plasma area to the total image area (Equation (7)). SI takes in consideration the upper and lower limits (maximum and lower PR values) and calculates the percentage of change in that time interval. Evidently, the higher the value of SI the more susceptible is a color channel to change. SI follows the formula:

$$SI = \frac{|\min(V\_{cl}) - \max(V\_{cl})|}{\max(V\_{cl})} \times 100\tag{6}$$

*Vch* is the value of a color channel at a specific objective magnification, and at a specific point of interest (maximum or the minimum value).

Figure 3 shows a flowchart representing an overview of the algorithm. Initially the input data, an AVI file, is preprocessed in three steps: adjustment of the frame rate (reduction needed fps), define and crop the region of interest, and separate the resulted colored images into their red, green, blue and gray components. Then, the first statistical analysis is performed to produce relevant indices and other parameters (i.e., threshold values for binarization, and normalization). At the second part (second row in the flowchart) the previously processed images are converted into binary, and various statistics are calculated. Finally, in the last part of the algorithm (last row of the flowchart), Equations (2)–(5) are implemented to derive the four new regions in each binary channel and to calculate the relevant statistics and indices.

**Figure 2.** (**a**) Representative binary comparison between parts of a reference and a subsequent image, and the resulted output matrices. (**b**) The behavior of *RBCs* ∗ and *CF* ∗ after clotting start (CS) and the correction process.

**Figure 3.** Flowchart of the proposed algorithms showing the three stages.

Figure 4 below shows the result of the processing algorithm for the 20× objective and for the gray level images, as the representative of the cases. As it can be seen the clot formation is detected as negligible in the initial time (*t* = CS), following a substantial increase at the half time of the coagulation process (*t* = CT), and reaching a maximum at the end of the time (*t* = 35 min). The area occupied by the RBCs is decreasing in time as the clotted regions extent on the sample. As noted earlier the plasma area (not shown) is complementary to the RBC area at the early stages, as coagulation progresses however, the plasma and RBC areas are occupied by clots, resulting in the CF images.

**Figure 4.** Representative images and the resulted red blood cell (RBC) and clot areas (CF) in white (W) pixels, on which the relevant indices are defined.

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

Three approaches have been utilized in this work for the whole blood coagulation assessment. In the first approach the mean intensity (MI) for each color channel was calculated for the three different magnification objectives. The images were separated into their four color channels in order to detect differences between them, regarding blood coagulation. The mean image intensity was examined in order to assess the overall effect of light blocking and the characteristics of whole blood coagulation time. The MI index for all color channels in Figure 5 has been normalized with its maximum value (and denoted as *MI* ∗) in order to compare the behaviour of the different color channels. In this figure the time response of the RGBGs channels at three different magnifications is illustrated for the whole experimental period. In all cases, there is an increase in the average detected light until approximately the 10th minute, most probably due to the phenomenon of RBCA, in combination with the initial stages of fibrin formation [31]. After RBCA reaches a maximum and any cell movement has ceased, blood coagulation is the dominant process affecting blood microstructure. Indeed, after approximately 10 min from the peak of *MI* ∗, there is a significant decrease in the MI curves as light is blocked due to the structural changes in the sample, caused by the formed clots.

*MI* ∗ curves appear in the same form for the specific magnifications, with some differences between the color channels. Differences are observed for the initial parts of the curves, however it should be noted that for each magnification a different test was performed. For example, in Figure 5a a spike like behavior around the second minute is observed, most likely due to a movement of blood components. In Figure 5c for the 50× objective, an obvious difference exists at the beginning of the process between the red and green color channels. In this magnification, small displacement of blood in the microchannel may have a greater effect in the results. After the maximum of the curves (which define the CS time), all curves have a similar behavior, with small differences between the experiments (i.e., the three different magnifications) and the color channels. The CS times observed in the present work for the non-reacted normal blood sample are in the vicinity ~10 min, agreeing with times observed in other studies. For reacted blood samples the onset of the coagulation process is detected at much smaller times depending on the reaction agent [16,35,36,38].

**Figure 5.** Normalized mean intensity *MI* <sup>∗</sup> of the RGBGs responses for the three different magnifications (10× in panel (**a**), 20× in panel (**b**) and 50× in panel (**c**)). The maximum value of MI is used for the normalization.

In the second approach the images were converted to binary, as described earlier. The binarization of the data resulted in the separation of the image into regions of plasma, RBCs and clots, and subsequently to the definition of the *RBCs* ∗, *NPA* ∗, *CF* ∗ and the *PA* ∗ indices. These indices were calculated for all consecutive frames from the CS time onward. As blood coagulates, and clots increase in number and size, the total number of pixels corresponding to plasma decrease could be assessed using the *PR* ratio:

$$PR = \frac{\text{Total Plasma Area}}{\text{Total number of pixels}} \tag{7}$$

*PR* is presented in Figure 6 for the 10×, 20× and 50× objectives, showing similar behavior between channels. This is expected, according to the *MI* ∗ curves in Figure 5a,b, which have a similar behavior after CS. Figure 6c shows a greater reduction in the initial value of *PR*, probably due to a lower threshold value in the processing stage caused by the higher objective magnification.

**Figure 6.** Plasma reduction (*PR*) as the ratio between the RBCs/clots and the plasma, of sequential images in time after CS. Panels (**a**–**c**) show the results from the 10×, 20× and 50× magnifications respectively.

Information on clot formation, provided by the *CF* ∗ index, is illustrated in Figure 7. Figure 7a,b correspond to 10× and 20× objectives respectively and show that *CF* <sup>∗</sup> (clot area in the image) has similar time response and maximum values at the steady state. Differences are apparent regarding the maximum values between each color channel at steady state. In the 50× objective case (Figure 7c) the shape of the curves for the color channels remained similar to the other two objectives, however the maximum *CF* ∗ values detected are higher. All *CF* ∗ curves have an exponential characteristic, and show a rapid growth of the forming clots in the early minutes after CS.

One of the objectives of this study was to identify the most efficient approach in terms of blood coagulation assessment using brightfield microscopy. Therefore, the results for *MI* ∗, *PR*, and *CF* ∗ can be further analysed for sensitivity utilising the sensitivity index SI (Equation (6)). Tables 2–4 present the CS time, a characteristic coagulation time (CT), and SI for the three objectives used and for all color channels. CT is defined as the "half-life" time of the coagulation process after CS, i.e., the time required to reach the 50% of the *CF* ∗

growth, or the 50% in *MI* ∗ and *PR*, with maximum/minimum values defined at the 95% of the developed curves. The mean value and the standard deviation (SD) from the RGBGs channels are included for each index in the tables.

**Figure 7.** *CF* <sup>∗</sup> expressed as percentage of the image area that is occupied by clots for the three magnifications (10× (**a**), 20× (**b**) and 50× (**c**)).

**Table 2.** Clot start time (CS), characteristic coagulation time (CT) and sensitivity index (SI) for mean intensity (*MI* ∗), plasma reduction (*PR*) and clot formation (*CF* ∗) in the 10 times objective case.


**Table 3.** CS, CT and SI for mean intensity *MI\**, plasma reduction *PR* and clot formation *CF* ∗ in the 20.


**Table 4.** CS, CT and SI for mean intensity *MI\**, plasma reduction *PR* and clot formation *CF* ∗ in the 50.


Inspecting Tables 2–4 is apparent that the Mean CS times in the *MI* ∗ analysis show good agreement, with values of 9:18 ± 0:11, 10:54 ± 0:10 and 10:46 ± 0:14 min and seconds, for the 10×, 20× and 50× objectives, respectively. Similarly, for all objectives the CT of the

*MI* ∗ index shows similar Mean and SD values. Small differences in the Mean CS and CT values between the objectives could be attributed to the manual handling of the sample placement, and the sampling initiation time of the experiments. SD of the metrics CS, CT of *MI* ∗, and CT of *PR* are very close to each other for all the color channels in the individual experiments.

Despite the similarities in the general behavior of the indices, minor discrepancies are also apparent. For example, for the 10× and 20× objectives the Mean and SD of the SI for the *PR* index are −38.44 ± 9.54% and 37.93 ± 9.87%, respectively. For the 50× objective, however, this metric is equal to −92.85 ± 1.96%, a much higher mean value with smaller SD than the other objectives. This might be due to the higher magnification in the 50× lens, where there is a smaller sampling area, and therefore small cell movements, due to RBCA and/or the formation of new clots, have a greater effect in allowing/blocking the transmitted light. In the *CF* <sup>∗</sup> analysis, the difference in the SI values between the 10× and 20× with the 50× objective is smaller than previously observed, with Mean and SD values of 43.34 ± 4.22%, 39.64 ± 3.20% and 51.07 ± 6.44% for the 10×, 20× and 50× objectives, respectively.

The results of the present study are comparable with those from recent and earlier works exploiting the optical properties of blood coagulation. These works include the laser speckle contrast imaging (LSCI) methods [41], which are also non-contact techniques. LSCI methods quantify light scattering and propagation, using wavelengths in the infrared/near infrared light spectrum [36–38]. These studies have established an exponential behavior of the assessed coagulation indices, with characteristics depending on reagent and sample used. This exponential behavior is observed in the present study for the non-activated coagulation of blood in Figure 7a–c. Moreover, the present results agree both qualitatively and quantitatively with data obtained using TEG instruments, in which CS times of approximately 12 min and CT of approximately 4 min are observed [42].

A popular method of choice, in several of the aforementioned and other studies in the literature, for the extraction of plasma, or whole blood coagulation parameters, is the calculation of the (cross- or/and auto-) correlation coefficients in consecutive images. Correlation methods can quantify the similarities between two signals (images), however with a computational cost related to the number of arithmetic operations, usually scaling with the square of the sample number (image size and number of images). The total number of operations in the binary comparison used in the present study, scales directly with the sample number. The binary method utilized in this work results into specific regions of interest, depending on the desired complexity (four regions in the present case). The separation of the four regions requires a small amount of processing power due to low algorithmic complexity, as shown in the Equations (2)–(5). Also, the visible light source offers true color images that can be used to further analyze whole blood coagulation process using for example feature extraction algorithms.

#### **4. Limitations and Sources of Uncertainty**

The blood sample used in the study was provided by the same donor, therefore negligible differences in the sample are expected between the various tests. The repeatability of the technique is still to be addressed in detail, however, an indication for the robustness of the technique is given by the similar behavior of the indices, which have been obtained from different tests in each magnification. As mentioned earlier, some uncertainty is expected in the CS and CT times, due to sample handling and manipulation as the sample placement was performed manually. Environmental factors (e.g., room temperature and illumination, etc.) are not expected to have a significant effect, as they were controlled and consistent throughout the experiments. The in-house construction of the microchannels may have resulted in small geometrical differences, however no significant influence is expected in that aspect. Accurate focusing of the microscope on the sample may be an issue for the high magnification cases (50× objective), however by inspection no significant issues were identified.

In terms of image processing, the opaqueness of the sample and the probable higher concentration RBC stacks from the RBCA process, can decrease the overall light absorption by the image sensor. In the absence of the standard test, these results can be compared with data from literature [34,38]. Nevertheless, the efficiency of the image processing method adopted in this work is reflected in the analysis of the results and more specifically in the SI data.

#### **5. Conclusions**

In the present study, indices for the assessment of coagulation in a drop of blood were developed, and experiments were conducted in custom made microchannels, in order to demonstrate the validity and the potential of the approach. The different magnification lenses used showed major similarities and minor differences, regarding the behavior of the various indices. The coagulation process in whole blood was found to have an exponential behavior in agreement with other studies in the literature. The developed methodology presented in this work is very promising for application in small scale diagnostic devices for the point-of-care. Specifically, the current setup (brightfield illumination/camera imaging/low magnification microscopy/simple channel) is ideal for the construction of inexpensive POC devices. Further, the direct statistical analysis of the images, in combination with the Boolean approach can result in reduced computational needs, therefore allowing for onboard processing in a POC device.

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

**Funding:** The Project EXCELLENCE/0918/0215 is co-financed by the European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation.

**Acknowledgments:** This study was supported from the Cyprus Research and Innovation Foundation grant with ref.: EXCELLENCE/0918/0215. The help of Dimitris Pasias and Andreas Passos is also gratefully acknowledged.

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

#### **References**


## *Article* **Simultaneous Imaging of Bio- and Non-Conductive Targets by Combining Frequency and Time Difference Imaging Methods in Electrical Impedance Tomography**

**Xue Bai 1, Dun Liu 2,\*, Jinzhao Wei 1, Xu Bai 1, Shijie Sun 1,2 and Wenbin Tian <sup>3</sup>**


**Abstract:** As a promising medical imaging modality, electrical impedance tomography (EIT) can image the electrical properties within a region of interest using electrical measurements applied at electrodes on the region boundary. This paper proposes to combine frequency and time difference imaging methods in EIT to simultaneously image bio- and non-conductive targets, where the image fusion is accomplished by applying a wavelet-based technique. To enable image fusion, both time and frequency difference imaging methods are investigated regarding the reconstruction of bio- or non-conductive inclusions in the target region at varied excitation frequencies, indicating that none of those two methods can tackle with the scenarios where both bio- and non-conductive inclusions exist. This dilemma can be resolved by fusing the time difference (td) and appropriate frequency difference (fd) EIT images since they are complementary to each other. Through simulation and in vitro experiment, it is demonstrated that the proposed fusion method can reasonably reconstruct both the bio- and non-conductive inclusions within the lung models established to simulate the ventilation process, which is expected to be beneficial for the diagnosis of lung-tissue related diseases by EIT.

**Keywords:** electrical impedance tomography; frequency difference; time difference; lung imaging

#### **1. Introduction**

Electrical impedance tomography (EIT) is a non-invasive imaging technique in which the absolute value or the change of conductivity distributions are reconstructed [1–6]. The underlying rotational is based on generating a current density distribution inside the medium and measuring the resultant electromagnetic fields by means of sensors positioned circumferentially on the boundary [2]. Current excitation and voltage acquisition are repeated continuously and rapidly until all the independent electrode combinations are deployed. Attributed to its notable merits of high temporal resolution, non-radiation, low cost, and non-intrusive measurement, EIT has also gained a great deal of interest in monitoring rapid changes in multiphase flow and medical imaging [2–6].

Ever since medical interest in EIT increased in 1980s, EIT has been developed to be a promising medical imaging technique to attain useful diagnostic information where there are obvious conductive changes in the biological tissues and organs with changes in physiological conditions, such as breathing, blood flow, digestive, canceration, and nervous activity [7,8]. Up to now, the medical applications of EIT research have mainly focused on thoracic imaging. Among the organs contained in the thorax, large volumes of air are Amoved into and out of the lungs, leading to obvious conductivity changes during ventilation. Therefore, it is agreed that monitoring of breathing by imaging the lung is the

**Citation:** Bai, X.; Liu, D.; Wei, J.; Bai, X.; Sun, S.; Tian, W. Simultaneous Imaging of Bio- and Non-Conductive Targets by Combining Frequency and Time Difference Imaging Methods in Electrical Impedance Tomography. *Biosensors* **2021**, *11*, 176. https:// doi.org/10.3390/bios11060176

Received: 20 April 2021 Accepted: 28 May 2021 Published: 31 May 2021

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

**Copyright:** © 2021 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/).

most important physiological application of this technique. Monitoring ventilation by EIT has already been applied for managing mechanically ventilated patients as it is can guide the set-up of mechanical ventilator pressure, volume, and respiratory rate [9].

Most clinical and physiological research into the lung with EIT is conducted using the difference imaging method to obtain air distribution images of ventilation and perfusion in the chest. Static imaging, i.e., reconstructing the absolute conductivity distribution of interest from a single set of measurement, has already been abandoned because of the poor image quality. Time difference EIT (td-EIT) was the consensus on the performance figures of merits for EIT image reconstruction after it was proposed by Barber and Brown [10]. Meanwhile, it was tested successfully for lung imaging as it can ameliorate modeling errors contained in the measurement data [11]. By taking the merits of frequency dependence of biological tissue, frequency difference EIT (fd-EIT) aroused significant interest by injecting currents with two or more different frequencies [11–18]. Compared with the td-EIT imaging method, fd-EIT does not need the reference measurement from the past, so it tackles the issue in terms of the unavailable voltage data with respect to the reference measurement in a real clinic environment [11–13]. Subsequently, developing a multi-frequency EIT (mf-EIT) system is an important transition for the research of EIT. Remarkable work includes the Sheffield system [14], the UCL systems [15], the KHU systems [16,17], and the Edinburgh system [3].

Many lung pathologies, such as atelectasis (lung collapse), pneumothorax, pleural effusion, and pulmonary edema, explicitly indicate an abnormity of air or liquid content in the lung. For example, the "pleural effusion" model in EIT contains a high conductivity region due to the accumulation of liquid in the pleural cavity [9,18]. Generally, diagnosis of several lung pathologies, such as avoiding acute respiratory distress syndrome (ARDS) with assessment of extravascular lung water, which is also known as quantify pulmonary edema and which can be considered as a biomaterial, may benefit from the use of the fd-EIT imaging method [19]. However, fd-EIT imaging faces a main challenge. Biomaterials such as lung tissues are only shown in the reconstructed images as their conductivity varies with frequency of the injected current. Thus, for lung imaging, the reconstructed images of fd-EIT are not well characterized in terms of attaining explicit interpretation of the pulmonary monitoring and air content. For the td-image reconstruction, it can be expected that the biological materials would induce much smaller conductivity variation than that of the non-conductive materials, and thus the biomaterial distribution would be hardly seen in the reconstructed images due to the existence of the large amount of air. Another limitation is that most EIT signals from ventilation may overwhelm those from perfusion. Therefore, EIT on lung monitoring by using only fd-EIT or td-EIT cannot stratify certain clinical requirements. In medical imaging applications, the images of different modalities have already been fused to obtain reliable and accurate diagnoses [20]. Inspired by the image fusion concept, it is meaningful to investigate the possibility of improving the image quality by fusing the reconstructed images of the td-EIT and fd-EIT since it is convenient to acquire both the td-EIT and fd-EIT data simultaneously.

In this paper, to simultaneously image bio- and non-conductive targets by EIT, it is proposed to fuse the reconstructed images by fd- and td-EIT to reveal both the bio- and nonconductive materials distributions, where the wavelet-based fusion method is employed. At the beginning, the td- and fd-EIT imaging methods to image bio- or non-conductive material distributions are investigated. This provides an in-depth look at the limitations of td- and fd-EIT imaging. The wavelet transform was performed to decompose and fuse the reconstructed td- and fd-EIT images, with the performances of typical wavelet basis functions compared to attain optimized image quality. Since the inherent properties of the lungs during breathing are well matched to the application of the proposed method, the lung ventilation imaging was studied. For the lung imaging, a simulation study was carried out by considering a homogenous phantom of the human thorax with lungs simulated as real-shaped abnormal inclusions containing air. The proposed method was evaluated by imaging the simulated ventilation process of the lung. In the end, experiments were

conducted to validate the proposed method by adopting an in vitro lung model made with a white gourd.

This paper is structured as follows: In Section 2, the image reconstruction in EIT is briefly reviewed, and the wavelet-based fusion method for td- and fd-EIT images is demonstrated. In Sections 3 and 4, both simulation and experimental results are presented to validate the proposed method with simulated or in vitro lung models. Finally, conclusions are drawn in Section 5.

#### **2. Methodology**

#### *2.1. Principles of EIT*

In this investigation the boundary was *∂*Ω. Generally, for computational needs, it is assumed that Ω represents the body under the finite element model (FEM), Ω is discretized into small pixels or voxels and assigned a piecewise constant conductivity distribution, i.e., *σ* = [*σ*1,..., *σ*N] <sup>T</sup> <sup>∈</sup> <sup>R</sup><sup>N</sup> is assumed within <sup>Ω</sup>. The electric potential inside the sensing area induced by the injected AC current can be written as [5,6]:

$$\nabla \cdot (\sigma(\mathbf{x}, \mathbf{y}) \nabla u(\mathbf{x}, \mathbf{y})) = \mathbf{0}, \ (\mathbf{x}, \mathbf{y}) \in \Omega \tag{1}$$

where *u* is the electric potential inside Ω. EIT is a nonlinear problem as the excitation current is a function of unknown conductivity distribution [21,22]. Therefore, based on the nonlinear model, the relation between the conductivity inside Ω and the boundary voltages at the boundary, *∂*Ω, can be expressed as:

$$V = \mathcal{U}(\sigma) + \varepsilon \tag{2}$$

where *V* is the measured voltage, *U* is the forward operator, and *e* is the additive noise and measurement error. Regarding the difference EIT imaging, Equation (2) can be linearized given a small conductivity perturbation Δσ, which can be expressed as (3):

$$
\Delta V \approx f \Delta \sigma \tag{3}
$$

where *J* is the Jacobian matrix, and calculated by solving the EIT forward problem, i.e., Equation (1). Difference EIT calculates a vector of conductivity changes between the measured conductivity distribution and the reference conductivity distribution. One commonly employed method is the single frequency td-EIT imaging. More specifically, the AC current with selected frequency is injected into the sensing area. Hereafter, a measurement at time *t*<sup>0</sup> is acquired as the reference, and then the voltage change at time *t*<sup>1</sup> against the reference measurement is employed for image reconstruction, which is expressed as:

$$
\Delta V = V\_{t\_1} - V\_{t\_0} \tag{4}
$$

In most cases, a set of reference measurements is acquired when the EIT sensor is filled with the background medium. To improve its precision, the measured data is typically averaged over many data frames.

Recently, incorporating the advances that have been made on the hardware system, a more attractive method in biomedical applications is the frequency-difference (fd-) EIT imaging. For the fd-EIT imaging, the current source is required to output a pair of complementary currents with two different frequency components, i.e., *f*<sup>1</sup> and *f*2, which are injected to the sensor simultaneously or separately. At each frequency component, the boundary voltages are measured. The difference between the two sets of measured voltages is adopted for image reconstruction, which is calculated by

$$
\Delta V = V\_{f\_2} - V\_{f\_1} \tag{5}
$$

After attaining the voltage difference, to solve the inverse or image reconstruction problem in EIT, a variety of algorithms have been proposed in the literature [23]. In general, using the regularization technique is the predominant way. The conductivity distribution can be obtained by solving the optimization problem:

$$\arg\min \left\{ \frac{1}{2} \parallel J\Delta\sigma - \Delta V \parallel\_2^2 + \lambda R(\Delta\sigma) \right\} \tag{6}$$

where *R* is the regularization function incorporating a priori knowledge and *λ* is the regularization parameter. In this case, for lung EIT, traditional and proprietary algorithms have been an obstacle to interpretation of EIT images as the reconstructed images cannot reveal the character of the lung [2]. GREIT (Graz consensus reconstruction algorithm for EIT) is a type of regularized imaging approach for 2D linear EIT reconstruction of the lung, which is based on the experience with a large number of linear regularized EIT reconstructions. One of the advantages of GREIT is that it conforms to the performance requirements for lung application, while other reconstruction algorithms do not. One can refer to Adler [2] for the detailed description of GREIT.

#### *2.2. Wavelet-Based Image Fusion*

Image fusion has three levels: pixel level, feature level, and decision level [20]. Pixellevel fusion transcends other image fusion levels as it retains more detailed image information and can achieve higher accuracy. Among the pixel-level fusion methods, the wavelet-based method is recommended because of its non-redundancy and directionality [24–27]. Figure 1 shows the image fusion process based on the wavelet transform.

**Figure 1.** Image fusion process based on wavelet transform.

The fusion process is performed as follows: Firstly, the source images were decomposed by applying the wavelet transform with an appropriate wavelet basis function. In this procedure, reconstructed images by fd- and td-EIT can be assumed as a continuous change of two-dimensional signal *<sup>f</sup>*(*x*1, *<sup>x</sup>*2) ∈ *<sup>L</sup>*<sup>2</sup> *R*2 , *x*<sup>1</sup> and *x*<sup>2</sup> are its abscissa and ordinate, respectively. The wavelet transform of *f*(*x*1, *x*2) is expressed as:

$$WT\_f(a, b\_1, b\_2) = \left< f(\mathbf{x}\_1, \mathbf{x}\_2), \psi\_{a, b\_1, b\_2}(\mathbf{x}\_1, \mathbf{x}\_2) \right> \tag{7}$$

$$
\psi\_{a,b\_1,b\_2}(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{1}{a} \psi(\frac{\mathbf{x}\_1 - b\_1}{a}, \frac{\mathbf{x}\_2 - b\_2}{a}) \tag{8}
$$

where *ψ*(*x*1, *x*2) is a two-dimensional wavelet basis function, *a* is the scaling factor, and *b* is the displacement factor. By applying the basis function, the low-frequency and highfrequency components of the reconstructed image can be obtained, respectively. The lowfrequency components of the image give approximate features, while the high-frequency components provide detailed features [27]. Secondly, appropriate fusion rules are employed to fuse the high frequency and low frequency components. The choices of fusion rules are critical in this procedure. There are many fusion rules, such as larger absolute value of coefficients, weighted average method, and local variance criterion [23]. According to previous studies, it is suggested that the fusion rules should adapt to the features of the

source images. In this study, the fused image needs to maintain the position, size, edge, and shape of the target inclusions in the reconstructed images by fd- and td-EIT. In view of these requirements, the high frequency components were fused by taking the maximum absolute value in the source images. By combining the pixels correlation and regional variance in the image area, a hybrid fusion method was adopted for the fusion of low frequency components, as introduced by Wang et al. [27].

In the hybrid fusion rule, assuming that *L*(*X*) represents the low frequency coefficient matrix of the image *X* , let *G*(*X*, *p*) denote the regional variance of the region Q centred at p(m, n) in the low-frequency coefficient matrix of image *X*, which can be expressed as [27]:

$$G(X, p) = \sum\_{q \subseteq Q} w(q) |L(X, p) - \overline{u}(X, p)|^2 \tag{9}$$

where *L*(*X*, *p*) represents the value of the element at (m, n) in the low-frequency component coefficient matrix, and *u*(*X*, *p*) represents the average value of the elements inside Q. *w*(*q*) is the weight coefficient, which would increase if it approached p(m, n). *M*2(*p*) is the matching degree of the regional variance regarding the low-frequency coefficient matrix of the source images A and B at p(m, n), which is defined as:

$$M\_2(P) = \frac{2\sum\_{q \subset Q} w(q) |L(A, p) - \overline{u}(A, p)| |L(B, p) - \overline{u}(B, p)|}{G(A, p) + G(B, p)}\tag{10}$$

It is noted that the value of *M*2(*p*) varies from 0 to 1. The smaller the *M*2(*p*), the lower the correlation between the low frequency coefficient matrices. Let *t*2(0.5 < *t*<sup>2</sup> < 1) be the threshold of the matching degree. If *M*2(*p*) < *t*2, the option fusion strategy is adopted:

$$L(F, p) = \left\{ \begin{array}{ll} L(A, p), & G(A, p) \ge G(B, p) \\ L(B, p), & G(A, p) < G(B, p) \end{array} \right\} \tag{11}$$

Otherwise, the average fusion strategy is used:

$$L(F, p) = \begin{cases} \begin{array}{cc} \mathcal{W}\_{\max} L(A, p) + \mathcal{W}\_{\min} L(B, p), & G(A, p) \ge G(B, p) \\ \mathcal{W}\_{\min} L(A, p) + \mathcal{W}\_{\max} L(B, p), & G(A, p) < G(B, p) \end{array} \end{cases} \tag{12}$$

$$\mathcal{W}\_{\rm min} = 0.5 - 0.5 \left( \frac{1 - M\_2(p)}{1 - t\_2} \right) \tag{13}$$

$$\mathcal{W}\_{\text{max}} = 1 - \mathcal{W}\_{\text{min}} \tag{14}$$

#### *2.3. Image Quality Evaluation*

Quantitative measures, e.g., relative coverage ratio (*RCR*), relative image error (*RIE*), and structural similarity (*SSIM*) index are defined to evaluate the quality of reconstructed images, which are the differences and correlations between the true image and the reconstructed image [28]:

$$RCR = \frac{CR}{CR\_{Tunc}}\tag{15}$$

$$RIE = \frac{\parallel \hat{\mathbf{g}} - \mathbf{g} \parallel}{\parallel \mathbf{g} \parallel} \tag{16}$$

$$SSIM = \frac{2\mu\_x \mu\_y \cdot 2\sigma\_x \sigma\_y}{\left(\mu\_x^2 + \mu\_y^2\right)\left(\mu\_x^2 + \mu\_y^2\right)}\tag{17}$$

where *CR* denotes the coverage ratio defined as the ratio of the size of the inclusions to the total size of the sensing area. Correspondingly, *CRTrue* is *CR* of the true target distribution. *g* and *g*ˆ are the true and reconstructed images, respectively, *μ<sup>x</sup>* and *μ<sup>y</sup>* are the mean grey values in the true and reconstructed images, respectively, with *σ<sup>x</sup>* and *σ<sup>y</sup>* being the corresponding standard deviations, respectively. Quantitative evaluation of the reconstructed images of biological objects under different frequencies is challenging because the conductivity changes with frequencies and the ground truth is unknown, especially for the fd-EIT imaging where different images can be obtained with different frequency intervals. Therefore, the popular criteria, such as image errors and correlation coefficients, cannot be employed in this case.

#### **3. Characterization of td- and fd-EIT Imaging by Experiment**

A series of static experiments were conducted to characterize the td- and fd-EIT imaging and validate the feasibility of the proposed method.

#### *3.1. Experimental Setup*

An experimental system with a cylindrical tank sensor was established to characterize the fd- and td-EIT imaging, as shown in Figure 2. The diameter of the tank was 38 cm with four layers of EIT electrodes mounted axially. In each layer, 32 metallic electrodes (stainless steel screws, radius = 5 mm) were uniformly mounted on the inner surface of the tank. Sixteen of them (choosing one in the other two electrodes) in the second layer were employed for EIT imaging. A National Instruments (NI) PXI-based EIT system developed by Beihang University was employed in the experiment. The system mainly comprises three modules: the AC current source, the data collection system, and the PC with system control functions for image reconstruction. In the experiment, the cylindrical tank was filled with tap water as the background medium with a conductivity of 0.034 S/m. The pairwise injected current had a peak-to-peak magnitude of around 1 mA, while the frequency was set to 1 KHz, 5 KHz, 10 kHz, 20 kHz, and 50 KHz. With the adjacent excitation strategy, the excitation current was injected into a pair of adjacent electrodes. By multiplexing, the potential differences between each possible pair of adjacent electrodes were measured by the data collection system. The measurements were then sent to the PC for image reconstruction via a USB interface. Finally, an image of the target inclusions was reconstructed in MATLAB using the data received. For a single measurement, 20 frames of data were averaged to mitigate the measurement noise.

**Figure 2.** Experimental set-up.

Td- and fd-EIT were used to image three different phantoms, which were set up by inserting a conductive inclusion (a peeled radish in phantom 1), a non-conductive inclusion (a nylon rod in phantom 2), and both of them (in phantom 3) at a certain location in the tank. The diameters of the peeled radish and nylon rod were 60 mm and 55 mm, respectively. Their length was sufficiently long to extend beyond the water surface so that all tests were conducted with target inclusions homogeneous in the vertical direction. Figure 3a–f show the phantom pictures and corresponding distributions of the inclusions, respectively.

**Figure 3.** Experimental setup of three different phantoms. (**a**) Phantom 1, peeled radish; (**b**) distribution of phantom 1; (**c**) phantom 2, nylon rod; (**d**) distribution of phantom 2; (**e**) phantom 3, peeled radish and nylon rod; and (**f**) distribution of phantom 3.

According to previous work [29,30], electrical impedance spectroscopy (EIS) can reveal the frequency response of electrical impedance of bio-tissues, which can be employed to demonstrate the frequency-dependencies of the conductivities of the peeled radish and nylon rod. Figure 4 shows the EIS measurements of the peeled radish and nylon rod by an impedance analyzer (MFLI Lock-in Amplifier-500 kHz/5 MHz, Zurich Instruments, Zurich, Switzerland) over the selected frequency range, i.e., 0–100 kHz. The results of EIS indicated that the conductivity of the radish increased linearly with the increase in frequency, and the linear fitting formula was:

**Figure 4.** Electrical impedance spectroscopy of peeled radish and nylon.

Its RMSE (root mean squared error) was 0.001834, and the value of R-square (coefficient of determination) was 0.9877.

In particular, it is noted that the conductivity of the radish at 10 kHz was about 0.03 S/m, which is nearly the same as the conductivity of the background medium. With respect to the nylon rod, the EIS data illustrate that its conductivity stayed almost unchanged with frequency. This is consistent with the previous analysis and verifies the feasibility of the use of the radish and nylon rod for the dedicated study on td- and fd-EIT imaging.

#### *3.2. Experimental Results and Analysis*

By applying the GREIT algorithms for image reconstruction, the td- and fd-EIT imaging results of the three phantoms were attained, and the size of the images w 64 × 64 pixels. In fd-EIT, 1 kHz was selected as the reference frequency. The voltage differences between the measured voltages at all others frequencies and at the reference frequency were procured for image reconstruction. In td-EIT, the measured voltage when the EIT sensor was filled with the background medium was taken as the reference, i.e., the measurement at *t*1, and the voltage differences were obtained when perturbations were present in the

background medium under the tested frequencies. The td- and fd-EIT imaging results are shown in Figures 5–7.

**Figure 5.** Td- and fd-EIT imaging results of phantom 1. (**a**–**e**) are td-EIT images of phantom 1 at 1, 5, 10, 20 and 50 kHz, respectively; (**f**–**i**) are fd-EIT images of phantom 1 at 5, 10, 20 and 50 kHz, respectively.

**Figure 6.** Td- and fd-EIT imaging results of phantom 2. (**a**–**e**) are td-EIT images of phantom 2 at 1, 5, 10, 20 and 50 kHz, respectively; (**f**–**i**) are fd-EIT images of phantom 2 at 5, 10, 20 and 50 kHz, respectively.

**Figure 7.** Td- and fd-EIT imaging results of phantom 3. (**a**–**e**) are td-EIT images of phantom 3 at 1, 5, 10, 20 and 50 kHz, respectively; (**f**–**i**) are fd-EIT images of phantom 3 at 5, 10, 20 and 50 kHz, respectively.

Figures 5 and 6 show the td- and fd- imaging results of phantoms 1 and 2 at 1 kHz, 5 kHz, 10 kHz, 20 kHz, and 50 kHz. Figures 5–7 have the same color bars, and the gray value range was 0–255, corresponding to the conductivity value of 0–0.08 S/m, respectively. The same applied to the relationship between the gray value and the conductivity in Figure 8, despite of the use of different colors for presentation. As indicated, td-EIT can reconstruct

the position and shape of the single radish or nylon rob clearly (except for at 10 kHz where the conductivity of radish was close to the background medium as indicated by td-EIT imaging results). Fd-EIT could image the radish in phantom 1, but it failed to indicate the accurate position of the nylon rod in the reconstructed image since the conductivity of non-conductive nylon rod was not frequency-dependent. Regarding the phantoms 1 and 2, td-EIT provided better imaging results of the inclusions, as indicated in Tables 1 and 2 for quantitative evaluation, where td-EIT imaging results have lower RIE, as well as RCR and SSIM closer to 1. In Figure 5, td-EIT imaging results indicate that the color scale of the reconstructed inclusion changed with the increase in excitation frequency, which was due to the frequency-dependency of the conductivity of peeled radish. More specifically, the conductivity of peeled radish was lower than the tap water when the frequency was lower than 10 kHz, and larger than the tap water at frequencies >10 kHz. On the other hand, since the conductivity of the radish was related to frequency, the image quality was sensitive to the frequency intervals. Besides, size overestimation and image distortion are more significant with increase in the frequency intervals, which is attributed to the decrease of SNR with the increase in frequency. The image distortion was the most serious at 10 kHz due to the conductivity of peeled radish being very close to the tap water and the measured voltages had a much smaller signal-to-noise ratio (SNR), which would lead to wrong imaging result since the EIT reconstruction was ill-conditioned.

**Figure 8.** Fused images by applying the wavelet-based fusion method: (**a**) Ture image, (**b**) TD imaging, (**c**) FD imaging, (**d**) fused image by coif3, (**e**) fused image by fk14, (**f**) fused image by dmey, (**g**) fused image by bior3.7, (**h**) fused image by sym6, and (**i**) fused image by db8.


**Table 1.** RCR, RIE, and SSIM regarding td- and fd-EIT imaging results of phantom 1.

**Table 2.** RCR, RIE, and SSIM regarding td- and fd-EIT imaging results of phantom 2.


For further demonstration, imaging results of phantom 3 by td- and fd-EIT are shown Figure 7. In Figure 7, the position and shape of the nylon rod can be well reconstructed by td-EIT, while due to the existence of the non-conductive nylon rod, the radish can hardly be seen as it causes much smaller responses at the sensor boundary. However, regarding fd-EIT imaging, only the peeled radish is shown in the reconstructed images. The reason is that the conductivity of the radish varies with excitation frequency while that of the nylon rod does not. These findings are verified by the evaluation metrics in Table 3.


**Table 3.** RCR, RIE, and SSIM regarding td- and fd-EIT imaging results of phantom 3.

#### *3.3. Image Fusion with Respect to td- and fd-EIT Imaging*

To evaluate the performance of the above-mentioned fusion strategy, the reconstructed images of phantom 3 by td- and fd-EIT at 5 kHz were taken as the source images. Meanwhile, to investigate the impact of wavelet basis function, various basis functions available in MATLAB were employed in the wavelet transform, including "coif3", "fk14", "bior3.7", "dmey", "sym6", and "db8".

As shown in Figure 8, the source images by td- and fd-EIT cannot either demonstrate the conductive inclusion (i.e., the radish) or the non-conductive inclusion (i.e., the nylon rod). By applying the wavelet-based image fusion strategy, the fused image can successfully recover both the conductive and non-conductive inclusions, which resolves the issues of tdand fd-EIT imaging. It is noted that different basis functions have significant impact on the quality of the fused images, especially for the radish. By employing the "coif3" or "dmey" basis function in the wavelet transform, the shapes of the nylon rob in the fused images are obviously distorted besides the shapes of the radish. Overall, the "fk14", "bior3.7", "sym6", and "db8" basis functions can provide fused images of reasonable quality, which are much better than those reconstructed by either td-EIT or fd-EIT.

#### **4. Investigation for Improving Lung Ventilation Imaging with Thorax Model**

To investigate the performance of the proposed method in medical EIT imaging of lung ventilation, simulations and in vitro experiments were conducted. A real-shaped thorax model, as shown in Figure 9, was built for FEM simulation. The thorax domain and the sub-domains of the lungs were created according to the build-in function of the thorax model in EIDORS [31]. It is well known that the whole respiratory cycle contains two stages: inspiration and expiration. In the simulation, the lung was simulated as varying states by setting different air contents in the given 2-D lung model. The model is assumed ideal as it does not consider the modelling errors regarding the non-homogenous background, the lung shape movement and the change of the chest due to the breathing in the real clinical situations. The following two test cases were designed:

**Figure 9.** 2-D view of thorax model built in simulation. (**a**) Thorax model with lung, (**b**) inspiration model, and (**c**) expiration model.

> Case 1: Lung imaging in the inspiration (Figure 9b) stage, during which the breathing would increase the air content, and air would almost fully fill the lung in the end of inspiration. For this reason, it is suggested to set up a sufficiently large amount of air content inside the lung model to simulate the inspiration.

> Case 2: Lung imaging in the expiration (Figure 9c) stage, during which the air content decreases with the breathing and becomes negligible in the end of expiration. In a similar fashion, a small amount of air content was set up inside the lung model to simulate the expiration.

> In the simulation, lung conductivity was considered as a function of the excitation frequency, according to literature [32,33]. The maximum and minimum conductivity values were set to 0.079 S/m and 0.11 S/m, respectively. The conductivity of background tissues was set to 0.2 S/m, as suggested by Liu et al. [6]. Sixteen electrodes were attached on the boundary of the thorax domain. The same current injection and measurement strategy employed in the previous experimental studies was also adopted in the simulation. The simulation was conducted with 2-D FEM approximation in COMSOL. To produce a reference for fd- and td-EIT imaging, the reference frequency was chosen to be 1 kHz, and the reference distribution for td-EIT imaging was assumed to be the thorax model only containing the background medium and lung tissues, as shown in Figure 9a.

> Simulation results are shown in Figures 10 and 11. In the reconstructed images, whose size is 84 × 128 pixels, the outermost dark solid line and the innermost solid dark line denote the true lung shape and air region, respectively. In the simulated cases 1 and 2, the fd- and td-EIT imaging can successfully reconstruct the lung-shaped inclusion or the air inclusion, respectively, but not both of them. In addition, the air inclusion is overestimated by td-EIT, while the lung-shaped inclusion is underestimated by fd-EIT. In fact, the end of expiration would reduce the content of air in the lung, and hence result in a higher RCR (the air domain is significantly overestimated) in the case 2, as shown in Tables 4 and 5. Due to the limitation of the reconstruction algorithm, the sharp edge and accurate shape of the inclusions were not well preserved in the reconstructed images. Note that the fd-EIT images at 1 kHz cannot be obtained since it is the chosen reference frequency.

**Figure 10.** Td- and fd-EIT imaging results in case 1 with simulated lung. (**a**–**e**) are td-EIT images in case 1 at 1, 5, 10, 20 and 50 kHz, respectively; (**f**–**i**) are fd-EIT images in case 1 at 5, 10, 20 and 50 kHz, respectively.

**Figure 11.** Td- and fd-EIT imaging results in case 2 with simulated lung. (**a**–**e**) are td-EIT images in case 2 at 1, 5, 10, 20 and 50 kHz, respectively; (**f**–**i**) are fd-EIT images in case 2 at 5, 10, 20 and 50 kHz, respectively.




The reconstructed images by td- and fd-EIT at 5 kHz and 50 kHz were used to demonstrate the improvement of the lung imaging by the proposed fusion method. As shown in Figures 12 and 13, the image fusion for lung imaging was fairly successful. Compared with the imaging results by either td-EIT or fd-EIT, the fused images cannot only provide the rational changes of air content, but also the position and shape of the lung-shaped inclusion. It is not trivial in practice to obtain the change of air content and lung shape information simultaneously since the latter can serve as a reference and there may, in addition to the air current, be tissue changes related to certain diseases such as pulmonary edema. Meanwhile, the size and shape of the air inclusion or lung-shaped inclusion in the fused image agree well with the actual ones. It is noted that the fused images at 50 kHz indicate a much larger image contrast in terms of the lung shape. This is in accordance with the finding that the conductivity variation of the lung tissues between the

measurement frequency (5 kHz and 50 kHz) and the reference frequency (1 kHz) increases with frequency.

**Figure 12.** Fused images by applying wavelet-based fusion method in case 1. (**a**–**f**) are fused images using different wavelet bases in case 1 at 5 kHz, respectively; (**g**–**l**) are fused images using different wavelet bases in case 1 at 50 kHz, respectively.

**Figure 13.** Fused images by applying wavelet-based fusion method in case 2. (**a**–**f**) are fused images using different wavelet bases in case 2 at 5 kHz, respectively; (**g**–**l**) are fused images using different wavelet bases in case 2 at 50 kHz, respectively.

To further validate the feasibility and effectiveness of the proposed fusion method, static in vitro experiments were conducted. These experiment was conducted in the same tank as in the previous experiments. The tank was also filled with tap water as the background medium, and two lung-shaped inclusions were made of white gourd. Two pairs of air-containing inclusions were inserted into the two lung-shaped inclusions to simulate the lung-filling conditions in the expiration and inspiration stages, respectively, as in the two simulation cases. Figures 14 and 15(a1,b1) show the phantom pictures. To quantify the air content inside the lung-shaped inclusion, the volume ratio was measured. The volume ratio was 3/20 and 6/20 for the expiration and inspiration setup, respectively.

Based on the visual inspection of Figures 14 and 15, the fused images provide a reasonable estimation of the shape and location of the lung-shaped inclusion as well as the air inclusion inside it at different frequencies. The td-EIT performs well in reconstructing the air inclusion, while the fd-EIT is better in imaging the lung-shaped inclusion. It is noted that the contrast variation caused by the change of air content is clearly shown. Therefore, it can be concluded that the lung imaging based on the proposed fusion method is feasible and effective.

**Figure 14.** Reconstructed and fused images of lung-shaped white gourd containing a varied content of air at 10 kHz. (**a**) 3/20 volume ratio: (**a1**) is the actual experimental image; (**a2**,**a3**) are fd-EIT and td-EIT images, respectively; (**a4**–**a9**) are fused images using different wavelet bases, respectively. (**b**) 6/20 volume ratio: (**b1**) is the actual experimental image; (**b2**,**b3**) are fd-EIT and td-EIT images, respectively; (**b4**–**b9**) are fused images using different wavelet bases, respectively.

**Figure 15.** Reconstructed and fused images of lung-shaped white gourd containing a varied content of air at 20 kHz. (**a**) 3/20 volume ratio: (**a1**) is the actual experimental image; (**a2**,**a3**) are fd-EIT and td-EIT images, respectively; (**a4**–**a9**) are fused images using different wavelet bases, respectively. (**b**) 6/20 volume ratio: (**b1**) is the actual experimental image; (**b2**,**b3**) are fd-EIT and td-EIT images, respectively; (**b4**–**b9**) are fused images using different wavelet bases, respectively.

#### **5. Conclusions**

This paper characterizes and compares fd- and td-EIT for the imaging of target distributions composed of bio-material inclusions and non-conductive inclusions, and the lung ventilation imaging is a potentially applicable scenario. Through an experimental study

with a peeled radish or a nylon rod as targets, it was found that td-EIT reconstructions can reconstruct the bio-material inclusion (the radish) and the non-conductive inclusion (the nylon rod) if either the former or the latter exists in the sensing domain. Fd-EIT could only image the bio-materials inclusion (the radish) as the conductivity of non-conductive materials is not frequency-dependent, but has significant distortions with an appropriate frequency interval. However, both the td-EIT and fd-EIT failed to simultaneously reconstruct the radish and nylon rod if both of them exist in the sensing domain. This is because the response to the non-conductive rod (i.e., a much larger conductivity contrast to the background medium) would overwhelm that to the radish for td-EIT imaging, which would occur when imaging certain lung diseases such as pulmonary edema. To tackle with this dilemma, it is proposed to fuse the td- and fd-EIT imaging results using a dedicated wavelet-based fusion strategy, which can integrate the merits of td- and fd-EIT. Based on the simulation and experiment with lung-shaped inclusions or in vitro lung models, it was found that the proposed method has much better reconstruction accuracy than either td- or fd-EIT, which can simultaneously reconstruct both the lung-shaped inclusions and air inclusions in the fused image. Meanwhile, the proposed method can reduce the image distortions by td- or fd-EIT, and provide the relative shapes, sizes, and positions of those inclusions. In summary, the proposed method can improve the performance of EIT in lung ventilation imaging, so that the diagnosis functions of EIT can be expanded to deal with more lung-tissue-related diseases. Although such a fd- and td-EIT system certainly requires the hardware to accommodate to both the EIT modes, which can be readily achieved in practical applications, this paper aims to determine how much improvement can be expected by applying such a fusion method and understanding and identifying the compromises involved in fd- and td-EIT imaging. Future work will focus on optimizing the image reconstruction and fusion methods to further enhance the image quality. On other hand, the proposed method should be tested in more realistic cases of lung imaging.

**Author Contributions:** Methodology, formal analysis, software, and validation, X.B. (Xue Bai); methodology, original draft preparation, review and editing, D.L. and W.T.; methodology and formal analysis, J.W. and X.B. (Xu Bai); funding acquisition, review and editing, S.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Beijing Advanced Innovation Center for Biomedical Engineering.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


*Review*

## **EMG-Centered Multisensory Based Technologies for Pattern Recognition in Rehabilitation: State of the Art and Challenges**

### **Chaoming Fang 1,**†**, Bowei He 2,**†**, Yixuan Wang 1, Jin Cao <sup>3</sup> and Shuo Gao 1,4,\***


Received: 27 June 2020; Accepted: 22 July 2020; Published: 26 July 2020

**Abstract:** In the field of rehabilitation, the electromyography (EMG) signal plays an important role in interpreting patients' intentions and physical conditions. Nevertheless, utilizing merely the EMG signal suffers from difficulty in recognizing slight body movements, and the detection accuracy is strongly influenced by environmental factors. To address the above issues, multisensory integration-based EMG pattern recognition (PR) techniques have been developed in recent years, and fruitful results have been demonstrated in diverse rehabilitation scenarios, such as achieving high locomotion detection and prosthesis control accuracy. Owing to the importance and rapid development of the EMG centered multisensory fusion technologies in rehabilitation, this paper reviews both theories and applications in this emerging field. The principle of EMG signal generation and the current pattern recognition process are explained in detail, including signal preprocessing, feature extraction, classification algorithms, etc. Mechanisms of collaborations between two important multisensory fusion strategies (kinetic and kinematics) and EMG information are thoroughly explained; corresponding applications are studied, and the pros and cons are discussed. Finally, the main challenges in EMG centered multisensory pattern recognition are discussed, and a future research direction of this area is prospected.

**Keywords:** multisensory; electromyography; pattern recognition; rehabilitation

#### **1. Introduction**

Monitoring and analyzing the patient's physiological information are of significance in the process of physical rehabilitation in order to evaluate the rehabilitation effect [1,2] and control auxiliary devices during the physical rehabilitation process. Conventionally, the physiological information is divided into physical and psychological information, e.g., muscle force information and the intention of the patient. To detect these two categories of information, various types of sensors, including electromechanical sensors (such as accelerometers [3,4], gyroscopes [5,6] and force sensors [7,8]), and biosensors (such as electromyography (EMG) [9–11] magnetoencephalography (MEG) and electroencephalogram (EEG)) have been utilized. Electromechanical sensors are capable of detecting the physical information effectively. For example, in [12], flex sensors and force sensitive sensors (FSRs) are utilized to detect the bending angle of the finger and the grasp force, respectively. The sensing information is then

used for the control of robotic fingers in performing ten different tasks. In [13], the integrated use of accelerometer sensors, pulse sensors, and ambient sensors is reported for recognizing sedentary behavior with an accuracy of 95%. Nevertheless, biosensors are not only proven to reflect physical information, but convert the bioelectric signal, which is broadly treated as a direct link to human psychological information, into interpretable voltage amplitudes, and therefore have unique advantages in psychological detection over the electromechanical sensors [11,14]. Note that in the rehabilitation area, human psychological information normally refers to human intention information, which can be used to examine the synchronization level between human movements and human thinking. Hence, in this article, the detection of human intention, instead of human psychology, is focused.

Among all biosensor-captured information, EEG,MEG, and EMG are the three most relevant signals to human intention [15–17]. Among them, the EEG signal is of weak robustness due to the shortage of noninvasive electrodes in collecting a surface EEG signal, failing to provide high signal quality. MEG techniques, e.g., MRI, can offer accurate and rich information, but the time-consuming issue and huge machine volume block its use in rehabilitation. In contrast, EMG techniques enjoy relatively higher signal-noise ratio (SNR) and robustness than EEG means (especially during movements), and process information much faster than MEG techniques. Therefore, EMG is a more preferred choice for intention detection in the rehabilitation field.

Two widely used techniques to interpret the EMG signal are signal intensity registration and pattern recognition [18,19]. The former approach is based on the correlation between EMG signal intensity and quantified muscle force level [20,21]. However, the EMG amplitude of a muscle is usually determined by many uncontrollable factors, e.g., EMG crosstalk [22] and variation in muscle force [23], which will not only generate misinterpretations of the muscle activities, but also result in potential safety risks when such information is directly utilized in the exoskeleton for the motion functions impaired patients [18]. Alternatively, the latter approach offers higher accuracy in decoding EMG signals owing to the higher level of information, e.g., a set of motions or the movement intention can be extracted via raw data to provide deep insights in examining user's body condition. EMG PR was not popularized for its high computational cost compared to its counterpart, until being promoted by the rapid development of electronics and information technology in recent years. Now EMG PR has been broadly applied to disease diagnosis [24], intelligent prosthesis control [19], and gains continuing attentions from relevant researchers.

Generally, the recognition scenarios can be divided into upper limb PR and lower limb PR. In early years, research on EMG PR mainly focused on the recognition of the upper limb movements with great differences, e.g., largely changed arm movements [25]. This is because the EMG patterns of such movements are significantly different and therefore is easier to be recognized. Nevertheless, fine upper movements and lower limb movements are also important in many neuromuscular disease analysis, such as Parkinson's disease and myasthenia. Therefore, it is also important to create a reliable relationship between EMG signal and patients' intentions and physical conditions for these two kinds of movements.

However, using the EMG signal to interpret fine upper limb and lower limb motions is challenging. In terms of the former, small movements in the upper limbs like the finger or wrist only result in slight differences in EMG signals, giving rise to difficulties in effectively distinguishing the EMG signal of one pattern from others. For example, a hand movement recognition method using single-channel sEMG is presented in [26]. This work reached an accuracy of 86.7% in classifying nine finger movements, and the recognition accuracy and the pattern types are all lower than the multisensory approach [27]. As to the latter, interactions between lower limbs and unpredictable environments during movements can generate unexpected EMG patterns. For example, the muscle activity of lower limbs is often related to the flatness of terrains [28], or the different locomotion modes such as obstacle crossing [29] and stairs ascending/descending [30]. Therefore, utilizing only the EMG signal cannot satisfy the demand for reaching robust recognition performances at these complicated scenarios without learning environmental knowledge. A possible solution is to concurrently acquire the external information of

the environment as well as the EMG signals, so that information from different dimensions can be integrated for further analysis and complementing each other.

In light of this, EMG-centered multisensory-based technologies emerge in recent years, and increasing demonstrations and attempts in three main rehabilitation scenarios have been reported globally. First, multisensory fusion EMG PR plays an important role in accurate and early diagnosis for neuromuscular diseases so that the patients could start the rehabilitation process before the diseases further develop. For example, in [31], the acceleration information of the muscles and forearm of the biceps brachii of the Parkinson's patient are collected simultaneously, 12 kinds of signal features were extracted for cluster analysis, and successfully achieved a diagnostic accuracy rate of more than 90%. Second, multi-sensor fusion can evaluate the rehabilitation process of patients and the effectiveness of the treatment in more dimensions. For example, in [32], the fusion analysis of the plantar force information of the Parkinson's patient and the myoelectric information of the anterior tibial muscle. The effect of a plantar sensory stimulation therapy on improving gait and motor output in Parkinson's patients were evaluated. Third, multisensory EMG PR is used as the control strategy for assistive devices such as intelligent prosthesis and exoskeleton [33–36], these auxiliary devices can effectively help patients in physical rehabilitation, and to a certain extent reduce the burden of the physical therapy. In addition, multi-sensor fusion also allows these devices to obtain higher security and a more comfortable user experience. For example, in [37], an exoskeleton hand is designed to help paralyzed people in rehabilitation training. It uses myoelectricity as an input to represent the patient's intention and uses angle information and gripping force information as feedback to provide reference and correction for EMG signals in real-time, ensuring that the rehabilitation action can be completed accurately. The above three application scenarios are conceptually depicted in Figure 1.

**Figure 1.** The conceptual depiction of the application scenarios of the electromyography (EMG) pattern recognition technique.

To provide a timely and systematic overview of the EMG based pattern recognition techniques using multisensory fusion strategies in the context of rehabilitation, this review article is composed. The paper is structured as follows: we first introduce the basic physiology knowledge relevant to EMG-based PR in Section 2. Then, in Section 3, a detailed explanation of the processing procedure and comparison of the state-of-the-art EMG PR is presented. Afterward, Section 4 explains the collaborations between EMG and multi-sensory information and presents diverse multisensory strategies successfully applied in different rehabilitation scenarios. Finally, the challenges faced by EMG PR and the prospect in resolving these challenges are summarized in Section 5.

#### **2. Physiology Background**

For reader convenience in understanding the working principle for EMG pattern recognition, some relevant physiological backgrounds, including the mechanism of EMG generation and the target patterns, are essential prior knowledge. This chapter will explain these physiological backgrounds, which helps readers to understand the characteristics of EMG signals and the design of recognition systems for the target patterns.

#### *2.1. EMG Signal Overview*

EMG is a type of technique used to evaluate and record a series of electrical signals that emanate from body muscles [38]. The principle of the generation of the EMG signal can be described as follows. Motor nerve cells produce electrical pulses under the control of the central nervous system in the cerebral cortex. These neural signals are transmitted to muscle fibers through axons and cause pulse sequences, which activate them to contract and produce muscle tension. Meanwhile, a current is generated in the human body that brings about transmembrane potential [39]. The transmembrane potential is the difference between the internal and external potentials of muscle cell membranes. When muscle cells are in a quiet state, the cell membrane potential is polarized. The potential difference between the inside and outside of the cell membrane under the polarization state is called resetting potential. Depolarization occurs when the cell is excited, and this trend will spread to around [40]. The corresponding action potential is defined as an electromyography signal. The above procedure is demonstrated in Figure 2.

**Figure 2.** The principle of generation of the EMG signal. (**a**) the structure of the neuro-muscular system. (**b**) the schematic of the EMG signal transduction in the nerve and muscle system.

#### *2.2. Human Movement Patterns*

As the EMG signal is a bioelectrical signal controlled by the mutual effect of the receptors and nerve system of human beings, the pattern of the EMG signal relies on both the user's subjective intention and the interactive environment conditions while performing specific actions. In this section, we will introduce some frequently studied upper and lower limb patterns in relevant studies.

The upper limbs of the human body include a forearm, elbow, rear arm, wrist, and hand. The muscle mass of the upper limb is generally small and slender, and the movement range of the upper limb is not extensive. Compared with the lower limbs, the muscle strength of the upper limbs is generally weaker because the body does not need to be supported.

The main functions of the upper limbs are grasping, stretching, and expressing information with hand gestures. Therefore, the patterns in the upper limb can be divided into two categories: limb position and hand gestures. The classification of Limb positions has been elaborated comprehensively in [25], including P1: arm hanging at the side, elbow bent at 90 degrees; P2: straight arm reaching up (45 degrees

from vertical); P3: straight arm hanging at side; P4: straight arm reaching forward. Hand gestures account for the dexterous movements of the hand, including wrist and finger movements. Wrist gestures are the hand movements that rotate the whole hand around the wrist joint, but finger gestures only involve the movements of the fingers. Reference [27] proposed a challenging set of 15 gestures in which 5 (thumb, index, middle, ring, pinky) are finger gestures, 6 are wrist gestures. The above patterns are demonstrated in Figure 3.

**Figure 3.** Two types of upper limb patterns. (**a**) the limb position, adopted from [25]. (**b**) hand gestures, adopted from [27].

The lower limbs of the human body include the hip, thigh, calf, foot, hip, knee, ankle, and other joint parts. The lower limb has the characteristics of larger muscle tissue and a larger joint movement angle [41]. Its function is mainly to support the human body and change the position of the human body. Depending on these physiological functions of the lower limbs, the activities most closely related to the lower limbs are often activities with the core of walking behavior, so the activity recognition based on the lower limbs can be roughly summarized in two aspects, namely, the locomotion mode and the gait phase.

The first one is locomotion mode, which aims at the process that the brain actively coordinates the limbs to make different actions when the human body is facing different environments in the real world. Specifically, up and downslope, walking on flat ground, up and downstairs are the most widely studied in the relevant research [35,42,43]. This is not only because these five sports modes are the most common and widespread scenes in life, but also because these modes have the process of overcoming the influence of gravity and changing the center of gravity of the body, so auxiliary equipment is needed to adjust the lower limbs [44,45]. In addition to these five kinds of sports modes, some work has also been done to study such sports modes as crossing obstacles [29,46], turning [29], standing, and sitting transition [47]. These locomotion modes are a more in-depth interpretation than simple actions, such as leg swinging, touchdown, etc., which is conducive to the auxiliary equipment to have a better understanding of the user's environment and generate adaptive control based on this.

Gait phase is another widely concerned lower limb mode, which is a scientific decomposition and utilization of human periodic walking movement. There are different ways to define the gait phases according to different usage scenarios. The simplest one is to divide the gait into two phases: namely, the stance phase and swing phase [3]. Generally, the heel strike (HC) and toe-off (TO) are used as the starting point of the stance phase and the swing phase, respectively [48]. However, this division is not accurate enough and is not applicable in scenes that require high continuity. In [42], a definition of the gait phase consisting of four phases is proposed and used for real-time recognition. In the most precise definitions, the gait cycle can be divided into eight phases [49].

The function of the gait phase is to provide auxiliary equipment to refine the human walking process and provide different auxiliary control in different gait segments [42,50]. The characteristics of the gait phase, such as continuity, duration, and so on, to describe the whole walking cycle of the human body, is a standard method in gait analysis and has been proved to be related to the diagnosis of some diseases [51]. Figure 4 shows the five most commonly studied locomotion modes, as well as a detailed definition of the gait phases.

**Figure 4.** The two commonly studied lower-limb patterns, i.e., locomotion modes and gait phase. Five locomotion modes, along with the definition of the eight-phases gait cycle, are presented here.

#### **3. EMG Pattern Recognition Pipeline**

The differences in human body motion patterns can be reflected in EMG signals. However, the raw EMG data is usually very noisy and cannot be directly classified. Therefore, a series of processing procedures is necessary to distinguish these differences accurately. Generally, the PR procedure can be summarized in the following four steps (also depicted in Figure 5):


There is no universal standard for these processing procedures yet, owing to the variability for EMG PR in different rehabilitation scenarios. Therefore, in this section, we will explain the EMG PR protocols and approaches adopted by relevant studies in detail.

**Figure 5.** The pipeline of performing EMG pattern recognition.

#### *3.1. Data Acquisition*

Acquiring a high-quality EMG signal is strongly desired to ensure a satisfying recognition accuracy level, so properly designing the front-end EMG acquisition approaches is very important. In this section, we will compare the methods in front-end EMG acquisition from two aspects: the design of the EMG sensing system and the selection of muscle measurement points.

#### 3.1.1. EMG Sensing System

The sensing system for collecting EMG signals is mainly composed of electrodes, amplifiers, microprocessors, and transmission devices. The electric signal generated by the muscle is picked up by the electrodes, amplified by the amplification circuit, and then transmitted to the host computer via the transmission device. Generally, the connection between electrodes and processor is wired while the transmission is typically wireless.

In the collection of EMG, the design of the sensor module is typically highly similar, except for the electrode. Common EMG electrodes include two main kinds of pole placement methods: monopolar electrode and bipolar electrode. Both placing methods measure the potential with reference to the electrode placed at the locations without EMG response (e.g., ankle or knee), while the Monopoles method directly measures the potential difference and bipolar method apply the differentiate amplification method [52]. Bipolar electrodes have a higher usage frequency because the common mode noise can be suppressed in real application but lack of setup flexibility when compared to monopoles. Besides, the EMG signal quality is also influenced by the distance between each electrode pole and their diameters and widths when bipolar electrodes are applied [38,52].

In addition to electrode placement, the implantation of the electrode is also reviewed here. Intramuscular EMG and surface EMG are two mainstream Electromyography signal acquisition schemes. Usually, intramuscular EMG signals are recorded using percutaneous fine wire electrodes or others made of similar materials. These electrodes need to be inserted into several muscle tissues using hypodermic needles [53]. Insertion locations are identified by palpation [54] and verified by electrical stimulation and EMG channel activity during corresponding test contractions. Different from intramuscular EMG signal collection, surface EMG (sEMG) electric potentials are acquired with electrodes placed on the skin just above the target muscle [52]. Electrodes occurring in previous experiments are usually made of silver. Typically, there are two ways for operators to contact the electrode with the skin. One is by using a silver-chloride gel to achieve wet contact while the other applies dry electrodes with microneedles to record EMG signals [43].

According to the literature research, surface EMG hosts a higher usage frequency compared with intramuscular EMG [52,55–57]. The main reasons for its wide application are regarded to be non-invasive and convenient. Modern surface EMG electrodes are stuck to the surface of muscle and avoid the danger of causing potential muscle injuries to tissue, which is considered to be the critical concern in intramuscular EMG. Besides, during the use of sEMG, operators can freely select measurement points on the surface of muscle tissue and reduce the time cost to do the preparation work like disinfecting the probe points with alcohol. Nevertheless, some bold attempts in the application of intramuscular EMG have also been provided recently [53,54]. Intramuscular EMG provides another signal source to detect the human body behavior and handles some difficulties related to sEMG-based control, such as the unstable contact of the electrode with the skin. Some other additional benefits of intramuscular EMG are given in [53], for example, the ability to record deep muscle signals with little EMG crosstalk. Although the experimental results in [53,54,58] all prove that intramuscular EMG could not bring with the reduction in classification error from sEMG for single classifiers while the significant decrease of classification error for parallel classifiers can be achieved. According to the practical efforts in [53,54], the parallel configuration for simultaneous control becomes more promising via the use of intramuscular EMG, which has long been an unsolved problem in sEMG studies.

#### 3.1.2. Muscle Site Selection

In the thesis of movement recognition, reasonable muscle group selection for EMG signal obtainment is of considerable significance for robust recognition. Various muscle groups in both the lower and upper limbs are functioning and positioned differently, and Figure 6 demonstrates the anatomical structure of the muscles.

Real muscle selection needs to consider the previous experience and specific demand in the application scenario. For example, as for prosthetic control, which is mainly designed for the amputees, experimenters commonly choose the proximal muscle groups like thigh muscle group [29,35,42] and the upper arms [59] because the distal ones are considered to be amputated. While, as for the studies oriented for non-amputees subjects, the distal muscle groups such as the shank muscle tissues and forearm and wrist muscles are often selected as targeted muscles due to their lower inter-subject variability than the proximal [43,60,61]. Gao Shuo [43] chose a pair of distal antagonistic muscles, tibialis anterior (TA) and soleus (SL) to conduct the terrain identification experiment and got satisfying results. Besides, muscle characteristic is also an important aspect to be taken into account. For example, the proximal hip muscle groups (AM, GM, and PRF) are taken into consideration in [9] due to their sensitivity to different walking speeds. It has been proved in [61] that rectus femoris (RF) and Soleus (SL) are responsible for pulling the body forward at the stance phase of the gait and therefore studies related to the recognition of uphill/upstairs behaviors demanding to overcome the gravity are more likely to choose such muscles for EMG detection [62]. The selected criterion and the specific muscle groups selected by relevant studies have been summarized in Table 1.

**Figure 6.** Longitudinal representation of the muscle groups of (**a**) the lower limb (shank and thigh). (**b**) the upper limb (elbow and hand).

#### *3.2. Signal Preprocessing*

The original EMG signals obtained through the front-end acquisition equipment are generally noisy, and pattern recognition requires the use of signals of limited length for further processing. Therefore, this section will introduce filtering and windowing preprocessing methods to reduce the noise in raw data and segmenting the EMG signals.

Real EMG signals are usually mixed with various kinds of noise due to environmental disturbance and other uncontrollable factors in the experiments. Reasonable preprocessing is necessary for extracting useful features in further analysis. Filtering, normalization, and windowing are the three most common EMG data preprocessing methods. The filtering technique derives from classical signal process schemes of signal preprocessing in the frequency domain. Most of the energy of the EEG signal are mainly within the frequency range 0–500 HZ [65]. In [66,67], it is stated that high-frequency (500–1000 HZ) EMG signal is likely to be interfered with by the aliasing. As for the low-frequency region (1–10 HZ), it is mainly regarded as a kind of noise caused by the cable movement and the interface of the measurement electrode and the skin. Actually, the most commonly used EMG signal filter is essentially a finite impulse response bandpass filter with a cutoff point of 10 HZ (low) and 500 HZ (high) [68–70].



Muscles: TA = Tibialis Anterior; SL = Soleus; SAR = Sartorius; RF = rectus femoris; VL = vastus lateralis, VM = vastus medialis; GRA = gracilis; BFL = biceps femoris long head; SEM = semitendinosus; BFS = biceps femoris short head; ADM = adductor magnus.

Due to the difficulty when dealing with the long-term sensor data sequence, windowing operation is necessary to slice it into short clips for the real-time prediction task. Generally, the properties of the signal decide the window length when analyzing the EMG data. According to [71], EMG signal can be considered a wide-sense Gaussian random process. Therefore, proper window length is of great significance for information extraction. The too-long window will lead to clips with an unbearable variance, while too shot sequence may not be able to contain enough useful information for classification [72]. Huang, He et al. [29] state that the analysis window length should not exceed 200 ms, which is an ideal upper bound to control the signal variation. Research in [8,72] suggests that windows of 150–250 ms for EMG are the optimal tradeoff between the classification accuracy and the delay, while 100–250 ms windows can be used for mechanical sensors. Fast time response is also desired for continuous classification tasks in addition to high classification accuracy [29]. In [42,72], a method called overlapping analysis windows was applied to accelerate the decision update. The key aspect of this scheme lies in the careful window increment setting. The smaller the window increment

step, the faster the classification result can be achieved. Furthermore, overlapping can help obtain better classification performance when processing the transition phase, such as from the sitting phase to the standing phase [42].

#### *3.3. Feature Extraction*

After the signal was preprocessed, the noise was attenuated and the window segments were obtained. However, using the entire segment as the input data is not an ideal choice due to the high computation cost and poor correlation between the input data and the target patterns, which would negatively affect the recognition performance. Therefore, it is necessary to extract some features from each window sequence. In this section, the two most commonly applied feature extraction approaches, including time-domain features, frequency domain features, and autoregressive features, are reviewed and compared.

The time-domain features of the EMG signal are based on the characteristic index of the statistical method, which regards the EMG signal as a function of time. This is an intuitive interpretation of the EMG signal. Its calculation is directly related to the magnitude of window segments, which occupies small calculation resources. In [71,73,74], Hudgins et al. proposed several time-domain features that are widely used in the following researches, including mean absolute value (MAV), mean absolute value slope, slope sign changes (SSC), waveform lengths (WL) and zero crossings (ZC). In [75], time-domain features like root mean square (RMS), waveform length (WL), number of zero-crossings(ZC), variance (VAR) and maximum(MAX) value are combined to represent the feature space.

The frequency-domain features of EMG signal are mainly to transform the time-domain signal of EMG signal into a frequency-domain signal through Fourier transform, and then analyze the spectrum characteristics or power spectrum characteristics of the signal. The advantage of this method is that it overcomes the characteristics of a time-domain signal, which is greatly affected by noise, is relatively weak, and it is easier to extract stable characteristic indexes by analyzing EMG signals in frequency-domain. Common frequency domain features are comprehensively investigated in [76]. Lots of experimental results demonstrate that mean frequency (MNF), median frequency (MDF), mean peak frequency (PKF), mean power (MNP), frequency ratio (FR), power spectrum ratio (PSR), and variance of central frequency (VCF) are not good features for EMG signal-based locomotion mode classification. Joint time-frequency domain features have been proved to be capable of effectively representing transient EMG patterns resulting from dynamic contractions [77]. The experiments comparing the time domain methods and time-frequency methods have also been conducted, and results show that the wavelet packet transform (WPT) is the best method to increase the EMG information density.

Apart from the conventional time and frequency domain feature extraction approaches, there are other ways of extracting effective features. For example, autoregressive (AR) features are usually taken into consideration in lower limb locomotion mode analysis. Huang. He et al., applied autoregression features (three-order autoregression coefficients) and time-domain features (MAV, ZC, and WL) for locomotion identification for their computation efficiency and fast time-response property.

In recent years, there are also handcrafted features that do not rely on statistical techniques. Instead, they use the prior knowledge of biomechanics to select the specific response point at the experimental curve. For example, in [43], twenty-one biomechanical features such as the peak EMG interval between two selected muscles and the duration of the EMG activation time are carefully selected. Such an extraction technique obtained an accuracy of 96.8% and is proven to outperform the traditional feature extraction technique. Although this kind of method seems lacking the general capability dealing with different scenarios, the classification performance is satisfying for the only specific task and shows less generality than the conventional features. There are also researches dedicating to improving the efficiency of the conventional feature extraction methods by designing new features. For example, in a most recent research [78], three time-domain features, including ASS,

MSR, and ASM, are specially proposed to improve the performance of EMG-PR based strategy in arm movement classification task. Experimental results showed that the new-designed features could achieve the accuracy 92.00% ± 3.11%, which is 6.49% higher than that of commonly used time-domain features [72].

#### *3.4. Dimensionality Reduction*

When the extracted feature set occupies a high-dimensional space, it may lead to the failure to correctly classify the patterns because the discriminative features are hidden in the complex feature space. Efficient data dimensionality reduction technique not only reduces the computation cost but also proves to be powerful in increasing the inter-class distance, enhancing recognition accuracy. Two commonly used dimensionality methods, i.e., feature projection and feature selection, are being discussed in this section.

#### 3.4.1. Feature Projection

The feature projection method means projecting the original high dimension feature space into a lower-dimensional feature space. The most important characteristic of such a method is that it uses all the original features to get the reduced feature space, which will not bring any information loss. Many algorithms are proposed to accomplish this task, including approaches like principal components analysis (PCA) [77,79], Non-negative matrix factorization (NMF) [80], averaging [75], independent components analysis [81], nonlinear projection [82] and. We will mainly introduce the PCA and NMF approaches in this section.

PCA is regarded as one of the most commonly used linear dimensionality reduction techniques. It essentially takes the direction with the largest variance as the main feature and off-correlates the data in each orthogonal direction, which makes them irrelevant in different orthogonal directions. Therefore, PCA also has some limitations. For example, it can remove linear correlation very well, but there is no way for higher-order correlation. For data with high-order correlation, Kernel PCA can be considered, and the nonlinear correlation can be achieved through the Kernel function turning to linear correlation. The research in [83] illustrates the importance of having relevant embedded muscle activity features in a low-dimensional space. It is also demonstrated in this paper that the PCA technique can efficiently capture features from EMG signals in an unsupervised manner. Due to the parameterless characteristic, PCA is convenient for universal implementation but having trouble for personalized optimization itself. Besides, some other methods about the strategic combination of features have also been leveraged in recent papers, and they are improving towards the direction with higher speed and accuracy. For example, FastICA is proposed in [84] to overcome the difficulty of previously used ICA for high-density EMG signal decomposition.

Non-negative Matrix Factorization (NMF) is another dimensionality reduction technique that has also been widely studied in recent years [80,85]. D.D.Lee et al. point out in [85] that NMF reduces data dimensionality by approximately factorizing the source data into two matrices. The factorization aims to extract the latent structure from the source data to a lower-dimensional space by minimizing the predefined cost function whose value is inversely proportional to the quality of the approximation. Common cost functions include conventional least square error and Kullback-Leibler divergence. As for its application in EMG signal recognition, G.R.Naik and H.T.Nguyen first studied the identification of EMG finger movement with a nonnegative matrix factorization technique [80]. Recently, non-negative matrix factorization is also leveraged to select the subject-specific signal channel for improving the lower-limb motor imagery recognition accuracy. Apart from this, NMF is also found to be put into practice in the determination of muscle synergies while utilizing EMG signals [86].

Actually, lots of other data dimensionality reduction methods have been investigated in recent years due to the quick progress in data science, like Uncorrelated Linear Discriminant Analysis (ULDA) [87], the locality preserving projections (LPP), neighborhood preserving embedding (NPE), discriminant analysis (OFNDA) and so on. These methods provide different dimensions to

increase the distinguishable information density, which will help improve recognition accuracy and processing speed.

#### 3.4.2. Feature Selection

Feature selection means reasonably selecting a subset of features to form the new feature space, while the rest of the original features would be abandoned [88,89]. It has been proved by numerous empirical results that feature reduction can reduce data dimensionality and the computational complexity meanwhile.

A key problem in the feature selection method is that what factors should be taken into consideration to define the best EMG feature space, Zardoshti-Kermani, Mahyar, et al. [90] came up with some standards to evaluate the EMG feature space including maximum class separability, robustness, and computational complexity.

Diverse approaches have been proposed to help select more representing features from the dataset containing too much useless data. Some work derives from combining the specific physical or statistical properties. In the task like limb posture classification, deterministic approaches have been proposed to select the optimal feature-channel pairs [91]. In this work, a distance-based feature selection to determine a separability index is utilized. Besides, to measure the amount of mutual information between features and classes, a correlation-based feature selection method was applied, and both schemes were proved effective to boost the posture classification accuracy.

Evaluating relevant features is another novel approach to make feature selection proposed in recent years [92–97]. One newest work formulates the feature selection problem as the optimization problem and proposes a personal best guide binary particle swarm optimization (PBPSO) to solve the feature selection, which works by evaluating the most informative features from the original feature set [93]. In another recent paper [97], a feature selection algorithm called ReliefF is given from the heuristic angle. This algorithm first selects the random instance of one of the database classes and then searches for the k nearest neighbor instances.

#### *3.5. Classification Algorithms*

Once the features with the reduced dimensionality are determined, classification algorithms should be deployed to distinguish the different categories of the extracted feature vectors. Numerous algorithms have been proposed, and in this section, we will review some of the most frequently utilized algorithms in EMG PR, along with the comparison between them.

Traditional classification methods include linear discriminant analysis (LDA) [98], support vector machine (SVM) [99,100], k nearest neighbor (KNN) [101,102], Bayesian analysis, fuzzy logic (FL) [103] and hidden Markov models(HMM) [104] while some modern methods have also been given like artificial neural network (ANN) [105] and convolutional neural network (CNN) along with the recent progress in deep learning research. All of the above methods are effective in classifying the extracted feature space while they each have specific characteristics due to their intrinsic difference at the algorithm level. Therefore, we will briefly explain the basic principle of each type of algorithm along with the comparison of their classification performance at different rehabilitation scenarios.

LDA is one of the most simple but effective methods which has attracted great attention in recent years. Linear discriminant analysis is utilized in for lower limb prosthesis movement mode recognition and achieves competitive accuracy of 97.45% even compared with neural network classifiers while deploying a PCA reduced feature set. In another recent work [105], LDA in a One-Vs-One topology was used for non-weight bearing lower-limb movement recognition immediately after training. Nevertheless, a fatal disadvantage for the LDA model is that it is not capable of handling the linear inseparable problem even if given perfect data. At this point, ANN enjoys the ability to describe nonlinear class boundaries among different categories. MLP and Cascade are the two main structures of artificial neural networks. Identification and classification of different gait phases according to the collected EMG data is a classical linear inseparable problem. In new research focusing

on this problem [104], three different MLP is deployed for the classification of two main gait phases while directly using the raw EMG data without feature extraction. The experimental comparison demonstrates that the performance of three MLP models achieves an average higher accuracy over all the popular models. On the other hand, in scenarios that some dimensionality reduction techniques like PCA have linearized the discrimination task, LDA can behave even better than the ANN methods. Besides, due to the well-defined classifier inner architecture, LDA possesses better performance stability. Nevertheless, for a given problem, it may be not easy to determine the optimal size and structure of an ANN. At last, in spite of the advantage of extracting latent features from unprocessed EMG data, the long training period is another thing that we need to take into consideration when intending to apply the ANN.

The SVM is estimated as one of the most popular approaches utilized in the movement mode classification. It is an effective data classification approach that projects the low-dimension data to the high-dimension feature space via kernel function. SVM works by finding a hyperplane to distinguish different categories in high-dimension space through the training process. It is believed that proper kernel function can help reduce the indivisible linear data to linear separable set in high-dimension space. From this, it is obvious that the most crucial problem when building the SVM model is to determine the optimal kernel function and its parameter values. Pires, Ricardo, et al. [99] uses the SVM algorithm to classify lower-extremity EMG signals during running shod/unshod with different foot strike patterns. It is also mentioned that kernel function selection, and parameter setting should depend on the specific task scenario. In [100], the SVM algorithm is utilized to develop an automatic classification system for lower limb hemiparetic patients. The empirical results demonstrate that SVM has the highest accuracy (95.2%) than KNN (89.2%) and ANN (92.3%). Particle swarm optimization (PSO) has long been a classical heuristic optimization algorithm. According to our literature review, the related works hybridizing the PSO and SVM to detect movement patterns are growing continuously in recent years. Zheng, Jiajia, et al. [106] employ PSO-optimized SVM to classify four gait phases, including initial contact, mid stance, terminal stance, and swing phase. The experimental results show that PSO-SVM exhibits the distinctive advantages on gait phase classification and improves the classification accuracy up to 32.9%–42.8% compared with the classifier based on vanilla SVM. Recently, some work compares the performance of common classifiers in different scenarios to provide references in practical application. The comparative analysis among NLR, MLP and SVM shows that for either classification performance and for the number of classification parameters, SVM attains the highest values followed by MLP, and then by NLR [107]. It should also be noted that SVM can obtain the highest classification performance while utilizing the lowest sampling rate.

Fuzzy Logic (FL) is another technique used in the classification of EMG signals, which achieves definite conclusions just from imprecise data input in a simple manner. It has been studied that FL exhibits the advantage of control techniques in biosignal processing [103]. FL can extract unrepeatable EMG features and mimic user's intent to make a decision. On the other hand, FL requires more system memory and processing time because of the use of fixed geometric-shaped membership functions in fuzzy logic limits system knowledge more in the rule base than in the membership function base [103].

The application of the hidden Markov model (HMM) is described in [108], where it is used for recognition of gait mode based on electromyographic signals. A modified Baum-Welch was used to estimate the parameter of HMM, and the Viterbi algorithm achieved the recognition of gait mode by finding the best HMM and state to assign corresponding phases to the given segments.

KNN classifier is another kind of common biosignal classifier based on traditional machine learning techniques. The authors in [101] use the K-nearest neighbor method to classify the EMG signals recorded from lower limb muscles during standing in individuals with complete spinal cord injury implanted with spinal cord epidural stimulation. Actually, most KNN-related researches in this area are correlated with other classification algorithms to find the most suitable classifiers for some certain problem. In [102], weighted KNN and SVM are both utilized as the classifier to distinguish the sEMG signals for controlling a prosthetic foot while the comparison of the accuracy of two classifies shows that weighted KNN obtains higher efficiency.

With the rapid development of deep learning technology in recent years, increasing attention has been paid to the deep feature classifier. For example, a backpropagation neural network was used in [109] to map the optimal surface EMG features to the FE angle angles. The experimental results show that the features extracted from multichannel surface EMG signals using deep belief network method proposed in this paper outperform principal components analysis (PCA), and the RMS error between the estimated joint angles and calculated ones during human walking is reduced by about 50%.

#### **4. Multisensory Fusion**

Through a series of processing, the change of movement patterns creates a stronger correlation with the EMG signal. However, due to the interactions between the lower limb movement and the external environment, as well as the low resolution of the upper limb hand fine movement in the EMG signal, the sensing strategy relying only on the EMG signal has great limitations in these two types of rehabilitation scenes. To address these issues, multisensory fusion techniques are proposed and proved to be capable of boosting pattern recognition accuracy by collecting more efficient data from different dimensions.

Generally, the EMG-centered fusion strategy can be divided into two types, i.e., fusion with kinematic sensors and fusion with kinetic sensors. The former strategy is mainly to supplement the three-dimensional motion information like human body acceleration, angular velocity, angle, etc. While the latter one is mainly to supplement the force information during the movement process. In this section, we will explain the detailed procedures of these two fusion strategies and demonstrate some successful applications using these techniques in different rehabilitation scenarios. Figure 7 demonstrates a block diagram of the above EMG-centered fusion strategy.

**Figure 7.** EMG-centered multisensory strategies and sensing blocks.

#### *4.1. Fusion with Kinematic Sensors*

Kinematic sensors are widely used in obtaining the locomotion properties of targets. In recent years, kinematic sensors have also attracted enough attention in human behavioral study combined with biosignal sensors (e.g., EMG). Kinematic information of human lower extremity describes the movements of the joints and limbs. Generally, the kinematic sensors can be divided into two main categories, inertial sensors and motion capture systems. Inertial sensors, including accelerometer, gyroscope, and magnetometer, are often used to collect kinematic data like joint angles, joint angular velocity, body orientation, and limb acceleration, etc. Motion capture systems became popular only in recent years, their biggest advantage being that they can collect locomotion information in a contactless manner.

Kinematic sensors' advantage of intuitively reflecting the 3D motion information of the limb helps compensate for the defect of EMG sensors, which only describes the whole physical process from the human nerve activity angle. Combing bioelectrical signal and motion information like acceleration has also been investigated to improve locomotion classification accuracy. Researchers in [110] especially investigate the improvement in classification error while utilizing EMG and kinematic data together rather than pure kinematic data. The empirical results demonstrate that the higher recognition rates can be obtained with combined data compared with vanilla kinematic information from the hip joint of non-dominant/impaired limb and an accelerometer.

The most successful use of the kinematics-EMG fusion strategy lies in the control of powered lower-limb prostheses. Such devices help people get rid of the disability and restore normal walking ability. Nevertheless, the adaptability of the prosthesis to other parts of the human body is still a challenging subject in which the kernel issue can be regarded as inferring human locomotion intent. Due to the restriction of the traditional recognition approach only relying on kinematic data, some recent works are making the steps to fuse EMG and mechanical information to improve the performance. The locomotion recognition system used for lower limb prosthesis control in [35] fuse the multichannel EMG signals and accelerometer measurements in feature level. The physical test shows that the empirical accuracy is over 95%, which is significantly better than traditional schemes, only applying surface Electromyography for identifying locomotion modes [29,35]. In [111], A combination of kinematic and electromyographic (EMG) signals recorded from a person's proximal humerus was used to evaluate a novel transhumeral prosthesis controller. Most especially, they trained a time-delayed artificial neural network to predict elbow flexion/extension and forearm pronation/supination from six proximal EMG signals, and humeral angular velocity and linear acceleration. Young, A. J. [15] studied the contribution of EMG data in combination with a diverse array of mechanical sensors to locomotion mode intent recognition in transfemoral amputees using powered prostheses. And this can indeed significantly reduce intent recognition errors both for transitions between locomotion modes and steady-state locomotion.

In addition, EMG centered kinematic fusion technique is also seen in applications of the hand gesture recognition. This is because the dexterous movement of the finger movements varies more significantly in the position information than the EMG information. Therefore, it is an ideal choice to combine them to obtain intention information and position information at the same time. For example, inertial measurement unit(IMU) and myoelectric units are utilized to implement the hand gesture-based control of an omnidirectional wheelchair in [112]. The system component, sensor placement, and the characteristic collected signals are demonstrated in Figure 8. In this system, the classification, which involves recognizing the activity pattern based on the periodic shape of triaxial wrist tilt angle and EMG-RMS from the two selected muscles, helps the accuracy improve to 94%.

**Figure 8.** The fusion strategy of combining kinematics data (IMU) and EMG data, adopted from [112].

Utilizing vision technique is another effective and direct way to detect the kinematic information of human. Recently, with the vigorous development of computer vision, some researchers are trying combining the vision and EMG information for limb motion classification task. In a newly published work [113], researches implemented a framework that allows the integration of multi-sensors, EMG and visual information, to perform sensor fusion and to improve the accuracy of hand gesture recognition tasks. For embedded applications, even-based cameras were utilized to run on the limited computational resources of mobile phones. The online results of hand gesture recognition using fusion approach reached 85%, which is 13% and 11% higher than utilizing EMG and vision individually.

#### *4.2. Fusion with Kinetic Sensors*

Kinetic sensors are a kind of instrument used to measure forces and moments that are directly connected with the movement of body segments like lower limbs. Various types of force transducers, including piezoelectric [114,115], strain gauged [116,117], and capacitive transducers [118–120], have been widely used in the design of force-sensitive sensors. Joint use of kinetic sensors and electromyography sensors is believed to be able to extract human muscle reflex related information, which will help the recognition of specific body behaviors. Generally, two types of kinetic information, interaction force, and ground reaction force can both be integrated with EMG information. The following part will introduce the fusion examples with interaction force and ground reaction force, respectively.

Interaction force sensors are typically embedded in a fixed structure and then used to detect the extension or flexion of limb, which are important for interpreting the intention of the human movement. The utilization of interaction force and EMG signals are also seen in the application of the upper limb hand gesture recognition in recent years. The research in [27] showed that combining EMG and pressure data sensed only at the wrist could support the accurate classification of hand gestures. Especially, the EMG is suited to sensing finger movements, the pressure is suited to sensing wrist and forearm rotations, and their combination is significantly more accurate for a range of gestures than either technique alone. Figure 9 demonstrates the implementation of this multisensory approach in hand gesture recognition.

**Figure 9.** EMG-FSR fusion strategy and hardware diagram used to classify different hand gestures, adopted from [27].

According to our literature review, ground reaction force (GRF) is also a commonly used kinetic information, especially in lower limb rehabilitation scenarios. A typical combination pipeline focus on lower limb prosthesis is shown in Figure 10. Such integration strategy is mainly applied in the lower limb exoskeleton or prosthesis because the lower limb is mostly responsible for interaction with the ground, and many locomotion modes are closely related to the terrain changes, which is directly reflected in the GRF information. Liu, Ming. et al. [35] applied the EMG measurement results and mechanical ground reaction forces/moments recorded using a six-degree-of-freedom load cell to build an adaptive classification strategy, which further enhanced the reliability of neutrally-controlled prosthetic

legs. The researchers in [60] presented a smart sensing system utilizing flexible electromyography sensors and ground reaction force sensors for locomotion mode recognition. Here, EMG and GRF information were collected from ten healthy subjects in five common locomotion modes in daily life. Rehabilitation robots is an attractive research zone which mainly helps the patients like stroke regain the locomotion ability. In [121], ground reaction forces and lower extremity electromyography are both utilized to compare inclined treadmill walking and turning conditions with its emulation on a developed balance assessment robot in order to investigate the feasibility of the existing approaches in rehabilitation.

**Figure 10.** EMG-GRF fusion strategy and the fusion pipeline to recognize different locomotion modes adapted from [27].

#### *4.3. Fusion with Both Kinematics and Kinetic Sensors*

In addition to the information fusion work introduced above, we also notice some recent work combing EMG with both kinematic and kinetic data. With the improvement of processor performance and the development of machine learning algorithms to handle high dimensional data, this triple-type sensor information fusion, which used to be regarded as a difficult approach, has been implemented and outperforms existing approaches in some motion classification tasks.

Rehabilitation exoskeleton is an important application of kinematic-kinetic fusion strategy. Because this kind of equipment needs various dimensions of information to maintain the highest stability and safety, it can achieve the goal of using EMG as the subjective control input and kinematic and dynamic data as the feedback. The work in [122] used human-robot interaction force, ground reaction force, and EMG sensors to distinguish the walking environment and gait period to provide crucial information for exoskeleton control. The results obtained using an individual mechanical sensor together with sEMG showed improvement compared to the case of only using force sensors, reaching the highest accuracy of 97.8% at gait phase recognition scenarios. The sensing blocks, as well as the location of the sensors, have been demonstrated in Figure 11.

**Figure 11.** Multisensory fusion strategy using EMG, kinematic sensors, and kinetic sensors together for the application of rehabilitation exoskeleton, adopted from [122].

Gait pattern analysis has also been widely used to assist in the diagnosis of muscle diseases. By multi-sensor data fusion, gait analysis gives objective and quantitative information about the gait pattern and the deviation due to the muscular situation of these patients, which are of great significance to the clinical research. Cui, Chengkun, et al. [121] conducted the simultaneous recognition and assessment of post-stroke hemiparetic gait by fusing kinematic, kinetic, and electrophysiological data. In [123], researchers investigated the gait pattern of 10 patients with myotonic dystrophy (Steinert disease) compared to 20 healthy controls through manual muscle test and gait analysis in terms of kinematic, kinetic and EMG data.

Table 2 summarizes the fusion strategy, targeted recognition classes, extracted features, proposed classifiers, and the final accuracy results of the multisensory application in different scenarios of rehabilitation.


**Table 2.** Summary of the pattern recognition techniques applied in relevant studies.


**Table 2.** *Cont.*

Classes: LW = level-ground walking; RA = ramp ascent; RD = ramp descent; SA = stairs ascent; SD = stairs descent; SO = stepping over an obstacle; IT = ipsilateral turning; CT = contralateral turning; SS = standing still; LP = loading response; MST = midstance; TST = terminal stance; PS = pre-swing; IS = initial swing; MS = mid-swing; TS = terminal swing. Feature: MAV = mean absolute value; ZC = zero crossing; SSC = slope sign changes; WL = waveform length; SL = signal length.

#### **5. Challenges and Future Development**

Currently, fruitful results have been achieved by EMG-centered multisensory pattern recognition methods in the rehabilitation field. However, there are still limitations in terms of accuracy, robustness, data volume, and flexibility when current approaches are applied to practical applications. In this section, we will explain the main issues and potential solutions and research directions of multisensory EMG PR.

#### *5.1. Low Data Quality*

Low data quality is a pressing problem in EMG detection, especially under a multisensory fusion scenario. Because the data quality of each dimension can significantly affect the final recognition accuracy. Since bioelectrical signals are generally more susceptible to noise interference, in EMG centered multisensory fusion, the data quality of EMG signals is the most important influential factor.

In Section 3, we explain the preprocessing of the EMG signal, whose purpose is mainly to reduce the noise in the raw EMG signal. Although by the preprocessing technique, the deterministic noise, e.g., charger noise, can be filtered and stochastic noise such as white noise can be compressed, interference comes from the source of the bioelectric signal that is challenging to be removed by data processing technique. First, the EMG crosstalk effect is an inevitable interference in EMG detection [22]. It means that the amplitude of the EMG signal is not only contributed by the targeted muscle but also influenced by the adjacent and surrounding muscles, hindering precise detection of the muscle activity. Second, electrode shift also results in severe EMG signal misregistration, since the traditional electrochemical gel electrode suffers difficulty in firmly attaching to the skin; a shift of merely 1 cm would reduce the classification accuracy by 5% to 20% [125]. Third, even when the patient performs the same locomotion under the same environment, the applied muscle force may slightly change due to the impossibility in perfectly controlling the body, resulting in the variations of the EMG signals of the selected muscles, and reducing the data quality [19].

Advanced sensing and detecting means could be possible solutions to the problem of low-quality data. With the integration and miniaturization of sensory circuits, high-density EMG (HD-EMG) has become a promising sensing approach. This type of sensing technique increases the coverage and density of the electrodes. Hundreds of electrodes are integrated with a small size overlying a restricted area of the skin [126,127]. The most significant benefit of HD-EMG is that the possibility of subtracting information at the motor unit level greatly improved compared with conventional EMG, minimizing the effect of EMG crosstalk at this precision level. The multiple electrodes are arranged in the form of a two-dimensional array, for example, in the shape of 8 × 24 [127] or 8 × 16 [126], which will produce an "EMG image". The EMG image presents a comprehensive view of EMG signals both in the time domain and spatial domain, providing far richer information than the conventional EMG electrode placement method of implementing one or two electrodes at a single muscle. For example, the HD-EMG arrays attached to the elbow area were used for the precise hand gesture recognition, reaching the accuracy as high as 99.0% in classifying twenty-seven gestures [126]. A typical procedure of utilizing the "EMG image" generated by the HD-EMG to conduct pattern recognition is demonstrated in Figure 12. Even though abundant information can be extracted with the help of HD-EMG, it also increased the burden of the acquisition and processing hardware. Because the sampling rate of EMG signals requires at least 1000 Hz at each channel, an 128 channel HD-EMG system will produce over 128,000 observations in one second, which would make the data scale enormous and lead to challenges in real-time processing.

**Figure 12.** The procedure of applying High-Density EMG (HD-EMG) in discriminating hand gestures, adopted from [126].

#### *5.2. Inadequate and Undisclosed Data*

High-quality data acquisition is the premise of high-precision multisensory pattern recognition. Meanwhile, for multi-sensor fusion, due to the large variability in the patient's muscle and physical condition and the difficulty in obtaining high-quality multisensory data, adequate and open-source data is of equal importance. This section will explain the necessity of establishing a big and open sensor data library in detail, and review some currently developing databases.

Currently, most of the current EMG researches were recruiting a relatively small volume of the volunteers to collect EMG signals, such as 5–20 healthy individuals or no more than five patients. The relatively small size is owing to the limited access to experimental resources due to the insufficient funding budget. Therefore, recruiting tens of participants is a balanced choice for verifying the proposed algorithms for most research groups. Nevertheless, there are a few concerns about this choice.

Firstly, the variation of the EMG signal is usually significant in different individuals. A strong correlation between EMG signals and the age, height, weight, exercise habits, and health status of individuals are found in relevant studies [128]. The sample size of tens of experimenters is often difficult to cover all these differences between various populations. Besides, researchers usually adopt their experimental protocols and different sensing equipment, and it is not easy to investigate the generality of the method proposed by one work across these influential factors.

Furthermore, the requirement of sensing equipment for robust EMG acquisition is quite high, and the collection system is typically expensive. It may potentially hinder some labs or researchers who do not have the experimental conditions from investigating EMG signals on the algorithm level.

An important approach for solving the above problem is to promote the concept of "big data-based EMG," which means that the verification method of EMG pattern recognition should shift from the verification in small sample experimental data to the performance of the proposed algorithm based on big data in different populations. Well-established datasets not only allows researchers to examine the quality and correctness of their acquired data more conveniently but testing the generality of their proposed methods in different datasets. In addition, larger datasets enable the application of some advanced algorithms such as deep learning in digging the underlying information in EMG signals [129].

In the past, the experiment data of EMG were mostly stored within the laboratories of individual researchers rather than publishing online, which made it hard to gather multiple datasets. In recent years, with the occurrence of "data paper," researchers are encouraged to publish their experimental data as an individual publication [130]. Such scientific publications promote the transparency of the raw data, and thus more EMG datasets are being available. In addition, there are also studies dedicated to establishing a standard dataset, trying to provide a standard experimental procedure and cover as many different muscles, movements, and population as possible. A few influential datasets such as Surface Electromyography for the Non-Invasive Assessment of Muscles (SENIAM), The Ninapro project [131] have emerged in recent years. However, there is still a long way to go in reaching the consensus of a widely accepted protocol and establishing a thorough dataset.

#### *5.3. Discrete Interpretation of Continuous Movements*

The development of the data quality and quantity offers the possibility to develop more effective processing techniques, especially in improving the simultaneity and continuity of the EMG-PR-based control strategy, which are two key drawbacks at the state of the art technique. In this section, we will explain how simultaneity and continuity problem is affecting the EMG PR performance, and provide possible solutions and research directions in the light of the multisensory technique.

As mentioned earlier, the EMG PR technique is widely used in the real-time Human-Machine Interface system (HMI). However, the broad use of clinical viable EMG based robots are still not reported yet. The potential problem is that a complete pipeline, including sensing, signal transmission, and processing, has to be done in order to recognize the current pattern. There is a time delay in this procedure, and it is acceptable in some scenarios, such as sign language recognition, when only static patterns like hand gestures are to be recognized. However, the patterns relevant to lower limbs are mostly dynamic and require real-time and instant recognition, and therefore the time delay may greatly influence the application of conventional PR technique in lower limb HMI devices.

There are several attempts to resolve the problem of discrete interpretation. One of them is to design the transition logic between two consecutive patterns skillfully, and some pioneers studies are trying to design the transition strategy. In [42], overlapped windows are utilized to overcome the time delay problem, and the comparison between the overlapped windows and non-overlapping windows are depicted in Figure 13. The increment of two consecutive and the overlapped window is only 12 ms, generating approximate continuous decision flow. The majority vote method is used to prevent misjudging problem. In addition, because the EMG signal is generated before the actual motion, there are also studies using the currently collected signals to predict the patterns in the next moment so that the robotic actuators could react in advance to adapt to the changing patterns [132,133]. However, the prediction method is only applicable in the periodical movement like gait, and still have difficulties in predicting more complicated motions. So far, there is still no consensus on the best transition method among EMG patterns, and the robust transition is a research interest in the EMG PR field.

**Figure 13.** A conceptual depiction of: (**a**) the overlapping windowing method. (**b**) the non-overlapping windowing method.

With the development of multisensory fusion and the increasingly thorough study of human body motion, the method of combining the biomechanical model with pattern recognition is also attracting increasing attention. Some open-source software, like OpenSim, is developed to provide present biomechanical models. With the help of the multisensory fusion, both the EMG signals and other mechanical signals collected from the patient's body can be used to calibrate and specialize the raw biomechanical models. After that, with proper mapping strategies, these models can be driven by EMG signals and, therefore, intuitively simulate the movement of human extremities, providing crucial information for continuous control of HMI systems. For example, a model-based control strategy is proposed in [134], using the EMG signal as the input for the carefully tuned and calibrated hand-finger biomechanic model, and managed to continuously control a prosthesis hand to finish the task of gripping in the response time of 16.2 ms per loop. The control schemes and multisensory strategy is demonstrated in Figure 14. Nevertheless, the model-based method requires tedious procedures to calibrate and tune the model, and it is proved difficult to separate the EMG signals into distinct functional movements for the subject with small arms. At present, it is still in the laboratory research stage.

#### *5.4. Future Analysis*

In the above three sections, we explain the challenges of multisensory EMG PR and the emerging technologies potentially capable to solve these challenges. We believe that in the foreseeable future, multisensory EMG PR has high possibility to develop in these three directions for the following reasons. Firstly, the application of the novel EMG sensing technique, e.g., High-Density EMG, will greatly improve the EMG signal's quality. Besides, the spatial resolution of the EMG signals will also be enhanced from the level of distinguishing different muscle tissues to the level of distinguishing single cells, making it possible for EMG signals to apply in detection of very fine movements. Secondly, more open source databases with high data quality, complete muscle measurement points and rich application scenarios will be built, which will further reduce the obstacles of self-designed experiments and attract more researchers to join the field of EMG signal analysis. Thirdly, with the reduction of EMG research barriers, the research of EMG analysis algorithm will continue to increase, and more advanced

algorithms to solve the traditional PR problem, such as the overlapping window and biomechanical model mentioned in this paper, will be proposed and verified.

**Figure 14.** The schematic of utilizing biomechanical based and EMG centered multisensory technique in controlling a robotic hand, adopted from [134].

#### **6. Conclusions**

Employing electromyography (EMG) centered multisensory integration pattern recognition (PR) techniques demonstrates strong potential in precisely interpreting patients' intentions and locomotion modes, which are of significance for analyzing the progress of rehabilitation status. However, commercialized products widely applied in hospitals or home-based internet of health things (IoHT) have not been reported yet. This article aims to reveal the reason behind this phenomenon by thoroughly reviewing the state-of-the-art research outcomes and providing a comprehensive analysis of limitations. The content provided here not only allows readers to have a full picture view of this area but also shows potential directions in further enhancing the development of applying EMG centered multisensory fusion PR techniques in rehabilitation associated applications.

**Author Contributions:** Conceptualization, C.F., and S.G.; methodology, C.F. and B.H.; validation, C.F., B.H., and S.G.; investigation, B.H.; resources, S.G., J.C.; data curation, B.H.; writing—original draft preparation, B.H.; writing—review and editing, C.F.; visualization, Y.W.; supervision, S.G.; project administration, C.F.; funding acquisition, S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Natural Science Foundation under Grant 61803017 and Grant 61827802, and in part by the Beihang University under Grant KG12090401 and Grant No. ZG216S19C8.

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

#### **References**


© 20201 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 (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Biosensors* Editorial Office E-mail: biosensors@mdpi.com www.mdpi.com/journal/biosensors

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com