Based on [12], we supported the notion that tinnitus heterogeneity influences the observed variability in treatment response after an analysis of collected data of 5017 tinnitus bearers where participants reported which treatments they tried, the duration and the outcome of the given treatment, alongside with the demographic and tinnitus characteristics. Sound therapy can effectively suppress tinnitus, at least in some patients [40], but there is still a lack of research on the efficacy of sound therapy. It is necessary to analyze the characteristics of individual tinnitus patients and to unify the assessment criteria of tinnitus [24]. In Tables 6 and 7, *p*-values above 0.05 and DID results suggest all the adult patients had a positive effect after the ADT-based treatment, whilst the elderly patient, under the same experimental conditions, had a negative effect. Furthermore, the subject who faced a significant improvement having the highest DID estimator and a similar statistical

distribution to the control groups before and after the sound-based treatment is the one with the lowest tinnitus intensity registered, alongside with low hearing loss in both ears.

Regarding treatment duration, it should be interpreted with caution, as it is wellknown that certain treatments require some time for adaptation, whereas other treatments require longer periods to be effective [12]. There is still uncertainty about the duration of treatment that may be required to achieve an improvement [25]. During this study, ADT-based treatments lasted 8 weeks. However, they were not applied for all patients even though 2 months is the minimum necessary time that has been empirically reported to find changes [12].

Tinnitus impairment can be quantified by various validated questionnaires such as THI. However, a recent analysis revealed a high variability in the outcome instruments used in clinical trials, indicating the need to standardize outcome measurement [9]. Furthermore, the outcome measures carried out through the THI in [12] were retrospective and subjective, which could have biased the results. This is why questionnaires are considered a subjective metric. According to [25], a further limitation of the current tools for assessing tinnitus impact is the reliability and repeatability of such measures: self-report measures of tinnitus have an associated risk of variability, as they supply a momentary snapshot, whereas the experience of tinnitus changes with time and context. Based on the previous evidence, it was proposed a first quantitative approach to objectively measure and evaluate the effects of ADT using ERD/ERS techniques along with the extraction of deep spectrum features. Significant relationship responses between the 'tinnitus group' after the sound-based treatment versus the 'control group' (Tables 5 and 6), positive DID estimators (Table 7), and close distance measures (Table 8) indicate the existence of neural modifications, which could explain why this treatment is so effective in some scenarios. Results from this research might help point out ADT as a potential solution for certain patients, but it is not a viable treatment for many others.

According to [24], patients with more severe initial tinnitus respond better to sound therapy; however, in the current study, the opposite results were observed. In Tables 6 and 7, *p*-values above 0.05 and positive DID estimators suggest that the subject who faced a better performance is the one with the lowest tinnitus intensity registered, alongside with low hearing loss in both ears. The elderly patient who did not benefit from acoustic therapy was due to the time he had suffered from tinnitus: around 30 years.

Our study comes with some inherent limitations. First, although we started analyzing 11 tinnitus patients, this number was reduced to 5 tinnitus subjects due to one of the following reasons: the rest did not show auditory material recognition responses in the initial monitoring session before receiving the ADT-based treatment or during the preprocessing stage, and the channel Fz was eliminated due to the transient or large amplitude artifacts. The final sample was insufficient, so it might not be representative of all patients with tinnitus. Second, the improvement trend is inevitable; however, it would be interesting to carry out a deep spectrum features analysis by theta, alpha, and beta bands to know exactly which cognitive demands are increasing or decreasing in terms of semantic, working memory, and attentional processes in each tinnitus subject compared with control subjects.

#### **5. Conclusions**

In conclusion, a new methodology based on ERD/ERS analysis and deep spectrum features extraction was successfully implemented to measure changes in attentional processes in tinnitus patients treated with ADT. Based on the previous implementation, our results pointed out that tinnitus attention was significantly reduced after the ninth week of an ADT-based treatment in adult patients. Furthermore, the therapy reported significant improvements in the patients with the lowest intensity recorded of tinnitus, alongside with low hearing loss in both ears. It is worth mentioning that this acoustic therapy is based on redirecting the attention that the patient has his tinnitus, this attention is focused on the deviant pulse of the oddball paradigm that is different from the frequency of the tinnitus. After eight weeks of treatment, the patient reports a reduction in the perception, but beyond

the reduction in the level of tinnitus perception, there is a reduction in the attention level, which results in the improvement of the patient.

Future work will entail measuring the EEG signals over the whole frontal lobe (Fp1, Fp2, F7, F3, Fz, F4, and F8). Furthermore, different neural network architectures could be applied to ensure the increase of the accuracy percentage in the classification stage to make the deep feature extraction stage more reliable.

**Author Contributions:** Conceptualization, I.G.R.-L., L.M.A.-V. and S.T.-R.; methodology, I.G.R.-L., L.M.A.-V., I.R.-G., S.T.-R. and R.A.S.-R.; validation, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G., R.A.S.-R. and D.I.I.-Z.; formal analysis, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G., R.A.S.-R. and D.I.I.-Z.; investigation, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G. and R.A.S.-R.; data curation, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G. and R.A.S.-R.; writing—original draft preparation, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G. and R.A.S.-R.; writing—review and editing, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G., R.A.S.-R. and D.I.I.-Z.; visualization, I.G.R.-L., L.M.A.-V., S.T.-R., I.R.-G., R.A.S.-R. and D.I.I.-Z.; supervision, I.G.R.-L., L.M.A.-V., S.T.- R., I.R.-G., R.A.S.-R. and D.I.I.-Z. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** The database used was created by following a protocol formerly approved by the Ethics Committee of the National School of Medicine of the Tecnólogico de Monterrey, described, published, and registered under the trial number: ISRCTN14553550.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The database used in this study is available at https://data.mendeley. com/datasets/kj443jc4yc/1 (accessed on 20 November 2021).

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

#### **Appendix A**

**Table A1.** Rejected channels, the percentage of bad data periods with transient or large-amplitude artifacts, and the independent components distinguished as non-brain sources.

**Figure A1.** Training and validation curves to display the classification performance of the MobileNet-V2 model used in the current research study.

#### **References**

## *Article* **Differences in Physiological Signals Due to Age and Exercise Habits of Subjects during Cycling Exercise**

**Szu-Yu Lin 1,†, Chi-Wen Jao 1,2,† , Po-Shan Wang 1,3 , Michelle Liou 4 , Jun-Liang Wu 5 , Hsiao Chun 1 , Ching-Ting Tseng 1 and Yu-Te Wu 1,6, \***

**Abstract:** Numerous studies indicated the physical benefits of regular exercise, but the neurophysiological mechanisms of regular exercise in elders were less investigated. We aimed to compare changes in brain activity during exercise in elderly people and in young adults with and without regular exercise habits. A total of 36 healthy young adults (M/F:18/18) and 35 healthy elderly adults (M/F:20/15) participated in this study. According to exercise habits, each age group were classified into regular and occasional exerciser groups. ECG, EEG, and EMG signals were recorded using V-AMP with a 1-kHz sampling rate. The participants were instructed to perform three 5-min bicycle rides with different exercise loads. The EEG spectral power of elders who exercised regularly revealed the strongest positive correlation with their exercise intensity by using Pearson correlation analysis. The results demonstrate that exercise-induced significant cortical activation in the elderly participants who exercised regularly, and most of the *p*-values are less than 0.001. No significant correlation was observed between spectral power and exercise intensity in the elders who exercised occasionally. The young participants who exercised regularly had greater cardiac and neurobiological efficiency. Our results may provide a new exercise therapy reference for adult groups with different exercise habits, especially for the elders.

 Introduction**

Regular physical exercise is associated with health benefits and is a crucial element of preventive strategies for promoting health. During exercise, moving the body requires a substantial degree of brain activity, necessitating the activation of numerous neurons to generate, receive, and interpret repeated, rapid-fire messages from the nervous system [1]. However, the neurophysiological mechanisms underlying the effects of exercise are poorly understood and require further investigation. Cycling is a common exercise, and daily cycling can enable a large proportion of the population to meet their recommended physical activity levels [2]. Several studies have reported that cycle ergometers are suitable for measuring physiological signals emitted during exercise. Studies on cycling exercise have reported that such exercise can induce specific changes in cortical activity. These changes are measured through various methods, including electroencephalography (EEG), the aim of which is to study the modulation of brain activity associated with performing cycling tasks [3–7].

**Citation:** Lin, S.-Y.; Jao, C.-W.; Wang, P.-S.; Liou, M.; Wu, J.-L.; Chun, H.; Tseng, C.-T.; Wu, Y.-T. Differences in Physiological Signals Due to Age and Exercise Habits of Subjects during Cycling Exercise. *Sensors* **2021**, *21*, 7220. https://doi.org/10.3390/ s21217220

Academic Editors: Vahid Abolghasemi, Hossein Anisi and Saideh Ferdowsi

Received: 30 September 2021 Accepted: 27 October 2021 Published: 29 October 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/).

Hottenrott et al. reported that cortical brain activation could be measured during cycling exercise; they suggested that higher cortical brain activation is necessary to increase muscle strength at higher cadences [4]. Enders et al. recently revealed that EEG power increased significantly in the frontal cortex and parietal cortex as fatigue accumulated throughout high-intensity cycling exercise activities. Notably, they observed a broadband increase in EEG power, in contrast to other studies that investigated various exercise conditions and observed changes that were limited to the alpha and beta bands [5]. Brummer et al. localized the exercise-induced changes in brain cortical activity by using the active-EEG/low-resolution electromagnetic tomography analysis and demonstrated that motor cortex activity increased with additional exercise intensity on a cycle ergometer [6]. Although Brummer et al. used different methodologies from other, earlier research, all of the aforementioned studies have focused on the effects of exercise intensity on cortical activity in young people or athletes [7]. Few studies have examined the activity of the cerebral cortex during exercise in other segments of the population, especially in older adults. Moreover, few studies have investigated the neurobiological differences between regular and occasional exercisers during physical exercise.

Accordingly, the aim of the present study was to investigate the changes in brain activity during exercise in elderly people and young adults. Previous studies have proposed the use of heart rate as a measure of exercise intensity [8]. They have described a positive linear correlation between increasing exercise intensity and changes in heart rate. However, because of age-related factors, the heart rate should not be directly used as an index for measuring exercise intensity. Santos reported that the aging process significantly alters the mean heart rate, which decreases with advancing age [9]. Therefore, the mean heart rates of young adults and elderly people at rest differ. In the present study, we used the average maximum heart rate ratio (AMHRR) [10,11], which can reduce the effect of age on the resting heart rate and maximum heart rate in response to exercise, and hypothesized that the AMHRR would facilitate the comparison of EEG and electromyography (EMG) readings between elderly people and young adults at the same exercise intensity.

We measured cardiac, cerebral, and muscular activity levels in elderly people and young adults in response to cycling exercise and investigated the differences between physiological signals obtained from four study groups: regularly exercising elderly people, occasionally exercising elderly people, regularly exercising young people, and occasionally exercising young people. In general, under a constant cycling period and intensity, regularly exercising young adults could achieve higher exercise efficiency with lower brain activation compared with the other participants. However, occasionally exercising young adults and elderly people may need to recruit more muscle units and increase the activation of the motor cortex during cycling compared with regularly exercising young adults. We hypothesized that physiological signal patterns would be similar between the occasionally exercising young adults and elderly people. We also anticipated that as age increases, the significant differences of physiological signals between occasional and regular exercisers may be more obvious in elderly adults than in young adults.

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

#### *2.1. Participants and Data Acquisition*

This study included 36 healthy young adults (18 men and 18 women aged 22.39 ± 3.56 years) and 35 elderly people (20 men and 15 women aged 64.65 ± 2.21 years) as participants. The elderly participants as well as the young participants were subdivided into 2 groups according to the time spent on exercise per week; specifically, participants who exercised for a total time of more than 3 h every week were considered as regularly exercising individuals, and those exercised for a total time of less than 3 h every week were regarded as occasionally exercising individuals [12]. All participants provided informed consent after receiving a detailed explanation of the purpose and potential benefits, and risks involved in the study. This study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of National Yang Ming Chiao

Tung University (YM106115E-1, 7 March 2019). Moreover, all participants were confirmed by physicians that their body mass index (BMI) was not overweight and without any lower limb or pelvic injuries, and had no brain-related diseases such as stroke, epilepsy, neurodegenerative diseases, orthopedic, or cardiovascular diseases. ECG, EMG, and EEG signals were recorded using V-AMP (Brain Products GmbH, Munich, Germany) with a 1 kHz sampling rate. The EEG channels included 10 wired wet electrodes, namely F3, F4, Fz, C3, C4, Cz, P3, P4, Pz, and A1, and were used according to the international 10/20 system (Figure 1a) [13]. The ground electrode was positioned at FPz. The EEG impedance level was maintained at <20 kΩ during the recording. The A1 channel was used as the reference for all electrodes. Ω

**Figure 1.** Locations of electrodes for EMG, ECG, and EEG. (**a**) location of electrodes for EEG (**b**) location of electrodes for ECG and EMG.

Electrocardiography (ECG) signals were recorded using 2 bipolar lead electrodes. The lead 1 (negative) electrode was situated below the right clavicle, on the mid-clavicular line within the rib cage frame; the lead 2 (positive) electrode was placed on the lower left abdomen, also within the rib cage frame. The surface EMG (sEMG) electrodes were placed on the quadriceps muscle (Figure 1b).

#### *2.2. Experimental Protocol*

We conducted an experiment to record EEG, ECG, and EMG signals while the participants performed the cycling exercise. Each of the participants sat on an electronically braked cycle ergometer in the upright position, with electrodes attached to their body. The study involved a pretest session and an experimental session. The pretest session involved 10 40 s stages of increasing workload with 20 s of rest between stages. For every participant, the workload ranged from 1 to 10. After the pretest session, the participants took a 5 min rest. The root mean square (RMS) amplitudes of EMG signals recorded for each stage were calculated, and the maximum RMS amplitude was considered the subject-specific maximum workload. For a participant, a workload corresponding to 40% of the maximum RMS amplitude was defined as the suitable workload for this participant. For safety, we assigned lighter exercise loads to the elderly participants to avoid injury or muscle damage due to over-load, considering the effects of declining physiological function with aging. Hence, the pretest session was considered to be excessively strenuous for the elderly participants, their suitable workload was set to 3.

In the experimental session, the participants were asked to ride the bicycle in 3 5-min exercise stages, resting for 30 s between stages. These 3 stages corresponded to relatively light, suitable, and relatively heavy workloads. EEG, ECG, and EMG signals were recorded simultaneously while the participants performed the exercise. Signals were also recorded for 5 min before the exercise (pre-exercise period) and for another 5 min after the exercise (post-exercise period). In this study, we required subjects to minimize their head and upper body movement as much as possible during the experiment. The participants were also asked not to move during the resting period. Figure 2 illustrates the overall experimental protocol.

**Figure 2.** Overall experimental protocol.

#### *2.3. Data Analysis*

2.3.1. ECG Analysis

ECG signals were detrended to remove low-frequency shifts, and the peak-to-peak R waves were identified to calculate RR intervals. The RR interval is the time elapsed between 2 successive R waves of the QRS signal on the ECG. We further used the AMHRR to monitor the status of each participant during the experiment [14]. The AMHRR can be defined as follows:

$$\text{AMHRR} = \frac{\text{averaged heat rate in each stage} - \text{RHR}}{\text{predicted maximal heat rate} (220 - \text{age} - \text{RHR})} \times 100\% \tag{1}$$

where RHR is the average heart rate during rest [10,11].

#### 2.3.2. EEG Analysis

For each participant, EEG signals recorded during the 5 min rest and during the exercise sessions were subjected to band-pass filtering between 1 and 45 Hz. Although participants were advised not to blink their eyes, clench their teeth, tense their muscles, or move their heads, these activities occasionally occurred and introduced artifacts into the EEG data. All signals with these artifacts were discarded during offline data processing. We further applied a moving average to the remaining signals for artifact suppression. Subsequently, each signal was divided into non-overlapping 1 min segments and then subjected to the wavelet transform [15].

The wavelet transform is based on small wavelets with a limited duration. The wavelet transform of a continuous-time signal *x*(*t*) can be defined as follows:

$$\text{WT}(a,b) = \int\_{-\infty}^{\infty} x(t)\psi\_{(a,b)}(t)dt\tag{2}$$

where

$$
\psi\_{(a,b)}(t) = \frac{1}{\sqrt{|a|}} \psi\left(\frac{t-b}{a}\right) \tag{3}
$$

ඥ|| is called the mother wavelet. The notations *a* and *b* denote the dimensionless frequency scale variable and time-like translation variable, respectively. The Wavelet transform enables the achievement of excellent localization both in the time domain through translations of the mother wavelet and in the scale (frequency) domain through dilations.

In this study, we used the Morlet wavelet [15] to transform each 1 min non-overlapping segment of an EEG signal (Figure 3b) in the 9 channels into temporal-spectral maps (Figure 3c). Each of these maps had 60,000 samples on the horizontal axis and 7 passbands namely 1–4 (delta), 4–8 (theta), 8–10 (low alpha), 10–12 (high alpha), 13–21 (low beta), 21–30 (high beta), and 31–45 Hz (gamma) Hz—on the vertical axis (Figure 3d). The spectral power

levels in each frequency band were averaged to obtain a frequency-averaged temporal power curve, which was again averaged across time to derive a frequency-time-averaged value. Thus, the average power per minute per frequency band was calculated. Each exercise stage was 5 min. Thus, the average power was calculated for 3 different workloads (Figure 3e). Subsequently, to normalize the average power for each exercise stage, this power was divided by the power at rest before exercise, thus yielding the normalized power (Figure 3f).

**Figure 3.** EEG signal analysis procedure. (**a**) Five-minute EEG signals during exercise. (**b**) Fiveminute signals divided into 1-min segments. (**c**) Temporal–spectral maps after the application of the Morlet wavelet transform on 1-min segments. (**d**) Temporal–spectral map divided into seven bands: delta, theta, low-alpha, high-alpha, low-beta, high-beta, and gamma bands. (**e**) Average power of each frequency band in each exercise stage. (**f**) Normalized average power of each frequency band in each exercise stage.

#### 2.3.3. EMG Analysis

The EMG signals were detrended to remove low-frequency shifts caused by the position fluctuations produced during the cycling exercise. The EMG signals were then subjected to band-stop filtering between 55 and 65 Hz for the removal of noise effects. After preprocessing, the signals were divided into 5 s segments (5000 sample points). RMS is usually used to predict muscle activity. Generally, a higher RMS value means higher muscle activity. RMS can be derived as follows:

$$\text{RMS} = \sqrt{\frac{1}{N} \sum\_{n=1}^{N} x\_n^2} \tag{4}$$

where *x* 2 *n* represents the EMG signal and *N* represents the length of the signal.

#### 2.3.4. Statistical Analysis

 ଶ

Pearson correlation analysis was used to evaluate the linear relationships between normalized power of EEG and the AMHRR or RMS of EMG. The AMHRR was considered an indicator of heart load for the various exercise stages. Thus, we could observe EEG and EMG changes with different exercise loads. In addition, paired-sample *t*-tests were used to examine for significant within-group changes before and after exercise (stage 3 and rest 2) to determine the post-exercise recovery status. In this study, MATLAB R2013b software (Mathworks, Natick, MA, USA) was applied for data analysis. Figure 4 illustrates a summary of the analysis procedures of EEG, ECG and EMG used in this study.

RMS = ඩ

1 ଶ

ே

ୀଵ

**Figure 4.** EEG, ECG, and EMG procedures in this study.

#### **3. Results**

#### *3.1. Changes in Heart Rate and AMHRR with Exercise Stages*

The mean heart rate and the AMHRR of the young and elderly participants during the different exercise stages are presented in Tables 1 and 2. Figure 5 illustrates the ECG analysis results for mean heart rate and AMHRR.The results revealed that in all groups, the heart rate and the AMHRR increased gradually with each exercise stage. The heart rates of the elderly participants were lower than those of the young participants. However, the AMHRR values of the elderly participants were not significantly different from those of the young participants, indicating that the cardiac load conditions of both the young and elderly participants were similar. The AMHRR was derived by normalizing the heart rate and excluding the effects of basal heart rate and age. Therefore, the AMHRR was suitable for observing the physiological state of the heart. We used Pearson correlation coefficient analysis to estimate the association between normalized EEG power and AMHRR per minute.

**Table 1.** ANOVA results for heart rate in young and elderly participants during different exercise stages. \* Indicates *p*-value < 0.05.

**Table 1.** *Cont.*

**Table 2.** ANOVA results for AMHRR in young and elderly participants during different exercise stages.

#### *3.2. Changes in EEG during Exercise in Young Participants with and without Exercise Habits*

We used Pearson correlation analysis to estimate the correlation between normalized EEG power and the AMHRR. Table 3 presents a summary of the regression coefficients of normalized EEG power and the AMHRR for all frequency bands. According to this table, a moderately strong correlation was observed, with the normalized coefficient ranging from 0.4 to 0.6. The results demonstrated that changes in EEG at the most frequency bands at C3, C4, and Cz were significantly and positively correlated with the AMHRR in both the young and elderly participants (*p* < 0.001). Moreover, the effect of exercise on EEG was mainly observed in the alpha band.

**Table 3.** Regression coefficients for the correlation between normalized EEG power and AMHRR during exercise. Boldface values represent moderate positive correlation between normalized EEG power and the AMHRR (correlation coefficient > 0.4).

**Figure 5.** ECG analysis results for mean heart rate and AMHRR.

According to Table 3, we could also observe the effect of exercise habits on normalized EEG power in the young participants. The correlation between normalized EEG power at C3 and the AMHRR was higher in young participants who exercised regularly, and the correlation between normalized EEG power at C4 and the AMHRR was higher in young participants who exercised occasionally.

#### *3.3. Changes in EEG during Exercise in Elderly Participants*

As presented in Table 3, the regression coefficients revealed a moderate or high correlation between normalized EEG power and the AMHRR in the elderly participants who exercised regularly. However, the correlation observed for the elderly participants who exercised occasionally was low and nonsignificant. Combining the results for elderly participants and young participants revealed that maintaining adequate exercise habits was more imperative for older adults than for younger adults. As illustrated in Figure 6, the elderly participants who exercised regularly demonstrated consistent EEG power changes. As the AMHRR increased, the normalized EEG power also increased. By contrast, no clear trend was observed for the elderly participants who exercised occasionally. The changes in EEG power were more dispersed. These results indicated that adequate exercise habits may lead to more stable brain wave changes in elderly people during exercise.

**Figure 6. Figure 6.** Scatter plots of correlation between normalized EEG power and AMHRR in low-alpha band (C3).

#### *3.4. Paired t-Test Results Observed during and after Exercise*

Figure 7 displays the normalized power values and statistical analysis results observed at C3 (low beta) at stage 3 and during post-exercise rest. Accordingly, the normalized power value during post-exercise rest would decrease to 1 if the power value during the pre- and post-exercise rest periods were identical. According to the plots in Figure 7, we observed the recovery speed of EEG power after exercise. The results revealed a significant difference in the change in normalized power between the exercise stage and post-exercise rest state in the young participants, regardless of their exercise habits. By contrast, in the elderly participants, the difference in the change in normalized power between stage 3 and post-exercise rest states was nonsignificant. Table 4 presents a summary of the results of the paired *t*-test for normalized EEG power in stage 3 and in the post-exercise rest state. In particular, the difference between the young and elderly participants was clearly observed in the beta band. The young participants recovered faster after exercise; therefore, a significant difference in the change in normalized power was observed. By contrast, the elderly participants recovered more slowly after exercise; hence, the difference in the change in normalized power was nonsignificant.

**Figure 7.** Histogram of EEG normalized power in low-beta (C3) band during exercise stage 3 and post-exercise rest.

Overall, the alpha and beta bands could reflect changes in brain wave power during and after exercise. The alpha band can be used to observe changes in brain wave power during exercise, and the beta band can be used to observe recovery in rest states after exercise. However, further understanding of the effect of age and exercise habits on EEG changes is warranted.

#### *3.5. Relationship between EMG RMS and AMHRR for the Four Test Groups*

Figure 8 displays the results of the linear regression on the differences in EMG RMS between stages 2 and 1 (i.e., ∆EMG RMS) and the difference in AMHRR between stages 2 and 1 (i.e., ∆AMHRR). Because of the extensive individual differences in EMG RMS values and the location of the EMG bipolar electrodes, we normalized the EMG RMS values; that is, we divided the RMS values for stages 2 and 3 by those for stage 1. This can be used to observe the increase in ∆EMG RMS with exercise load. The results revealed a more significant trend of increasing ∆EMG RMS with ∆AMHRR in the young participants than in the elderly participants. Additionally, the regression coefficients for the young participants who exercised occasionally were higher than those for the young participants who exercised regularly. However, for the elderly participants, a low correlation was observed between the ∆EMG RMS values and ∆AMHRR, regardless of their exercise habits, and their ∆EMG RMS values were more clustered. This low correlation may be because in this study, the elderly participants were assigned a fixed cycling load that was lower than those assigned to the young participants.

**Figure 8.** Linear regression between △△EMG RMS and △△AMHRR.

#### **4. Discussion**

This study used ECG, EMG, and EEG to explore changes in physiological signals transmitted during cycling exercise in young and elderly participants with different exercise habits. We assigned lighter exercise loads to the elderly participants to avoid muscle damage or injury from overload, considering the effects of declining physiological function with aging. According to previous research, exercise intensity (workload) is reflected in the response of many physiological processes, including heart rate [16]. Therefore, we defined

the exercise load according to the AMHRR and further observed changes in EEG and EMG with gradually increasing exercise loads.

#### *4.1. Spectral Power of EEG Increases with AMHRR during Exercise*

We observed that during exercise, the normalized power of each frequency band of the EEG signal was positively and linearly correlated with the AMHRR. We also determined that an increase in normalized EEG power was consistent with an increase in AMHRR. This consistency was observed in most EEG frequency bands, including the delta, theta, low-alpha, high-alpha, low-beta, and high-beta bands. Furthermore, these phenomena were more evident in the low-alpha, high-alpha, low-beta, and high-beta frequency bands. Earlier research reported that cortical activity increased with fatigue during exercise in order to maintain a constant physical output [1]. Schillings et al. also reported that the energy loss associated with fatigue during exercise may cause increased brain activation in the motor cortex [17].

Previous studies determined that during exercise, EEG cortical activation was most affected in the alpha and beta frequency bands [18–21]. Therefore, most experiments and literature reviews on the effects of exercise on EEG cortical activity were limited to these two frequency bands. Several previous studies involving ergometer cycling revealed that incremental graded exercise tests resulted in increased alpha power in the central and parietal regions as well as increased EEG current density in the primary motor region. Bailey et al. showed an increase in alpha and beta power after sustained intensity bicycle ergometer exercise with a progressively increasing workload [3]. Lin et al. reported increased EEG power in the alpha and beta bands in the frontal and central areas during high-resistance pedaling exercise [22]. They further proposed that the fatigue situation would be accompanied by an increase in α and β power. However, increased EEG beta activity may be associated with attentional demands and higher levels of arousal. Other studies demonstrated that the effect of exercise on EEG cortical activity was not limited to the alpha and beta bands [3,5]. Our results demonstrated that the alpha band was more suitable for observing changes in brain activation during exercise. However, the beta band was more appropriate for determining the differences between brain activation observed during exercise and that observed during post-exercise rest.

#### *4.2. Young People Who Exercise Regularly Have a More Coordinated Use of Their Dominant Leg*

Our results reveal that the EEG differences between young participants who exercised regularly and those who exercised occasionally were in the activation of motor cortical areas in the left and right hemispheres (i.e., C3 and C4). A higher correlation was observed between normalized power changes at C3 and exercise load in the young participants who exercised regularly. However, the normalized power at C3 and that at C4 in the young participants who exercised occasionally were moderately correlated. The concept of limb dominance was based on the fact that the two hemispheres of the brain function differently and tend toward activities that use one limb under voluntary control [23]. Bhise et al. observed that when for an inherently manipulative task, most participants used the dominant leg [24]. Young people who exercise regularly have greater coordination in the use of the dominant leg, meaning that they require only the dominant leg to complete the exercise. However, young people who exercise occasionally must use both legs to compensate for the deficiency of the dominant leg [25–27]. The RMS of EMG signals is often used as a concise quantitative indicator of muscle activity; we found that the young participants who exercised occasionally had significantly higher EMG RMS values than did those who exercised regularly. Our results indicate that the dominant legs of young people who exercise occasionally require more force output to perform a given task. However, that the young participants who exercised occasionally had lower EEG activation in the C3 region than did those who exercised regularly.

#### *4.3. Regular Exercise in Elderly People Induces Significant Cortical Activation during Exercise*

We observed that the highest increase in EEG normalized power occurred when the participants were at their highest AMHRR (exercise workload). This phenomenon was particularly notable in the elderly participants who exercised regularly. The results reveal that the normalized EEG power increased with the AMHRR in the elderly participants who exercised regularly, with the corresponding correlation being moderate to high. The heart rate increases with the delivery of oxygenated blood around the body and into the brain. Muscles require relatively high energy during exercise. Similarly, the brain consumes glucose or other carbohydrates when the body is in motion [28]. Therefore, the brain becomes more active during exercise. This suggests that elderly people who exercise regularly require relatively high exercise performance and muscle strength during exercise, which may induce considerable activation of the cerebral cortex. However, the change in normalized EEG power with respect to exercise load was less consistent in the elderly participants who exercised occasionally. Accordingly, the results reveal no significant correlation between normalized EEG power and the AMHRR. Although this phenomenon could also be observed in the young participants, the results were less pronounced than those observed in the elderly participants. Our results show that the difference in EEG signal changes between the elderly participants who exercised occasionally and those who exercised regularly was more significant than that between the young participants who exercised occasionally and those who exercised regularly. For elderly people, regular exercise can help reduce the functional decline associated with aging.

#### *4.4. EEG Recovery after Exercise Is Slower in Elderly People*

The paired *t*-test revealed significant beta band activation in the young participants in stage 3 and during post-exercise rest. By contrast, this phenomenon was not observed in the elderly participants. These results indicate that the young participants returned to a resting state more quickly after exercise, whereas the elderly participants required a longer time to recover. Aging affects the post-exercise recovery process. Several studies have revealed a functional decrease in the replenishment of energy supply before and after exercise. Research has presented evidence of differences in acute recovery of physiological parameters after fatiguing exercise between younger and older participants. For similar exercise stimuli, elderly people require a longer recovery period when returning to baseline levels after exercise [29]. Although this study did not reveal a significant difference in exercise recovery between the elderly participants with and without exercise habits, the EEG results demonstrate that the elderly participants who exercised regularly had superior brain regulation of exercise load than did those who exercised occasionally.

However, there are still some limitations in this study. First, the muscle artifacts occurring in the head and neck musculature during cycling exercise may be recorded in EEG signals. In this study, we asked subjects to minimize their head and upper body movement as much as possible during the experiment. Unfortunately, experimental protocols are still sensitive to physiological and non-physiological artifacts, including motion artifacts that may contaminate the EEG recordings. Following the procedure of artifact suppression, we applied a simple cleaning noise method, moving average, to remove the noise caused from motion artifacts in EEG signals. Although these procedures can eliminate motion artifacts but may also decrease the sensitivity in EEG signals. Second, the strength of muscle will decrease with aging, and there exist an individual difference in this aging effect. In this study, we did not take the muscle strength decay of elderly and individual muscle ability into consideration in the experiment setup.

#### **5. Conclusions**

This study revealed the AMHRR to be a suitable indicator of exercise intensity and that the physiological indicators of ECG and EEG in elderly people are different from those in young people because of aging. We found that the EEG spectral power of elders who exercised regularly revealed the strongest positive correlation with their exercise intensity. The results demonstrate that exercise-induced significant cortical activation in the elderly participants who exercised regularly, and most of the *p*-values are less than 0.001. No significant correlation was observed between spectral power and exercise intensity in the elders who exercised occasionally. The young participants who exercised regularly had greater cardiac and neurobiological efficiency. Therefore, appropriate exercise habits may benefit brain responsiveness and improve the efficiency of cardiac and neurobiological responses to exercise. Our results may provide a new exercise therapy reference for adult groups with different exercise habits, especially for the elders.

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

**Funding:** This research was funded by National Yang Ming Chiao Tung University, grant number 110BRC-B701, Ministry of Education (MOE),Taiwan, grant number 110W219 and the APC was funded by National Yang Ming Chiao Tung University and Ministry of Education, Taiwan.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of National Yang Ming Chiao Tung University (YM106115E-1, 7 March 2019).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study, and written informed consent has been obtained from the patients to publish this paper.

**Data Availability Statement:** The data are not publicly available due to the privacy concern raised by our IRB.

**Acknowledgments:** The authors thank all the subjects participated in this study and thank Wallace academic editing company for editing this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

## *Article* **Safe** *Hb* **Concentration Measurement during Bladder Irrigation Using Artificial Intelligence**

**Gerd Reis 1,\*, Xiaoying Tan 1 , Lea Kraft 2 , Mehmet Yilmaz 2 , Dominik Stephan Schoeb 2 and Arkadiusz Miernik 2**

**Abstract:** We have developed a sensor for monitoring the hemoglobin (*Hb*) concentration in the effluent of a continuous bladder irrigation. The *Hb* concentration measurement is based on light absorption within a fixed measuring distance. The light frequency used is selected so that both arterial and venous *Hb* are equally detected. The sensor allows the measurement of the *Hb* concentration up to a maximum value of 3.2 g/dL (equivalent to ≈20% blood concentration). Since bubble formation in the outflow tract cannot be avoided with current irrigation systems, a neural network is implemented that can robustly detect air bubbles within the measurement section. The network considers both optical and temporal features and is able to effectively safeguard the measurement process. The sensor supports the use of different irrigants (salt and electrolyte-free solutions) as well as measurement through glass shielding. The sensor can be used in a non-invasive way with current irrigation systems. The sensor is positively tested in a clinical study.

**Keywords:** hemoglobin sensor; bladder irrigation monitor; absorption near infrared; artificial intelligence; bubble detection

## **1. Introduction**

In urology, continuous bladder irrigation (CBI) is an important standard care procedure [1–4] after transurethral resection of the bladder (TURB) or the prostate (TURP). The dominant goal of CBI is to prevent the formation of blood clots and consecutive bladder tamponade, a medical condition requiring an additional and foremost avoidable follow-up surgery [5]. The purpose of CBI in the given application scenario is, therefore, to keep the blood concentration in the bladder at a very low level. Technically, CBI provides a continuous dilution of the bladder content with fresh irrigation fluid (often saline) and thus, prevents clot formation. Although nearly trivial from a pure technical point of view, CBI is involved when applied in clinical practice. Improper CBI may trigger bladder spasms by irritating the bladder, cause undesired bleeding [5], lead to bladder rupture or perforation [6], and might even become life-threatening [7]. These and other possible complications result primarily from increased pressure, due to a high flow rate. Consequently, there are two optimization goals for an optimally adjusted CBI: On the one hand, the flow should be high enough to dilute the blood sufficiently, and on the other hand, the flow should be as low as possible to avoid pressure-related complications. Since the amount of bleeding after surgery cannot be controlled, the optimal flow speed changes over time. Hence, CBI demands extensive and continuous supervision and management by medical personnel, imposing a heavy burden on nurses responsible for a whole urological station [8,9].

CBI supervision is, on the one hand, comprised of rather technical aspects, such as caring for a filled fluid reservoir and an empty waste reservoir, ensuring a continuous flow of fluids into and out of the bladder. On the other hand, there are aspects with high medical

**Citation:** Reis, G.; Tan, X.; Kraft, L.; Yilmaz, M.; Schoeb, D.S.; Miernik, A. Safe *Hb* Concentration Measurement during Bladder Irrigation Using Artificial Intelligence. *Sensors* **2021**, *21*, 5723. https://doi.org/10.3390/ s21175723

Academic Editors: Vahid Abolghasemi, Hossein Anisi, Saideh Ferdowsi and Stefano Bettati

Received: 25 June 2021 Accepted: 24 August 2021 Published: 25 August 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/).

relevance, in particular, choosing the right irrigation flow. Too high a flow will keep the blood concentration low, but may cause severe negative effects to the patient. In the case of the urethral catheter becoming clogged by a clot, continued irrigation may even lead to bladder rupture [5]. Too low a flow will not suffice to dilute the blood and will result in an ineffective procedure.

Today, nurses inspect the coloring of the waste fluid and estimate the flow speed accordingly. Obviously, this approach has significant drawbacks. Estimating the blood concentration in the outflow demands a high level of experience and is observer-dependent. Furthermore, it is affected by external conditions, such as the current illumination. Additionally, nutrition can have a significant effect on urine color, such as betanins leading to beeturia [10], easily confused with hematuria. Other foods known to change urine color include blackberries and rhubarb, which turn urine pink or red, while fava beans and aloe turn it reddish brown. Medications can also affect the color of urine, such as phenazopyridine, a drug used to numb urinary symptoms, the antibiotic rifampin, and laxatives containing senna, which can turn urine reddish-orange. The anti-inflammatory drugs sulfasalazine and phenazopyridine, as well as certain chemotherapy drugs turn urine orange-red, while the antimalarials chloroquine and primaquine, the antibiotics metronidazole and nitrofurantoin, the muscle relaxant methocarbamol, and laxatives containing cascara turn urine reddish-brown. In addition to the drug itself, the food coloring contained in the coating, e.g., in the case of tablets, pills, and dragees, can also discolor the urine. The discoloration of the urine and, thus, of the waste fluid is easily mistaken for an increased amount of blood, which in turn leads to an increased irrigation flow. The color of the excretory fluid is, therefore, not a reliable indicator for the adjustment of the flow. Detection of acute bleeding might become obscured, again putting patients at risk.

In order to better estimate the blood concentration, Hageman et al. [11] developed the Hemostick, a color scale to visually compare blood color. While this device significantly contributes to standardize the estimation of blood concentration between nurses, it still is affected by illumination influences. Furthermore only a limited number of discrete reference colors are available. Therefore, as can be observed in Figure 9, the distinction of blood concentration based on observation or color comparison becomes quickly unfeasible. Ding et al. [12] evaluated a CBI control system able to adjust the irrigation flow automatically based on the estimation of the blood concentration. The system features a color monitor to estimate the blood concentration. It thus resembles an automated version of the approach presented in [11]. In [13], Chan et al. presented a device to measure the blood concentration, using the light absorption principle. In order to estimate the blood concentration, they investigated light emitting diodes (LED) with different colors (i.e., red, green, blue) and finally decided on using green LED. Unfortunately, they did not specify the particular wavelengths of the used LEDs. To measure the transmitted light, a light dependent resistor (LDR) was used. Timm et al. [14] proposed a system based on light absorption for non-invasive estimation of hemoglobin (*Hb*) concentration in human tissue. They used three different, well-defined wavelengths to estimate the *Hb* concentration. The system was tested in a technical setting. Zhang et al. [15] described a system to estimate blood loss during endoscopic surgery based on *Hb* measurements. Their system is also based on absorption since they use a photoelectric sensor. Unfortunately, they did not provide any specifics on the sensor.

The new *Hb* concentration sensor proposed in this article also exploits the light absorption principle and thus, is closely related to the work in [13–15]. The advantage of using the absorption principle compared to color monitoring is the improved accuracy of the estimates. Independence from dietary coloring effects can be achieved through the selection of appropriate light frequencies. In contrast to [13], our approach is far more rigorous. The light absorption properties of blood and other coloring components were taken into account, and an optimal light frequency was selected. In contrast to [14], our system is dedicated to monitor CBI. The measurements take the light absorption of the irrigation tubing into account, as well as influences of the irrigation fluid. In contrast to [15], our system measures the *Hb* concentration in a defined tube in contrast to being integrated in a collecting bucket. In addition, practical aspects, such as gross mismeasurements due to bubble formation in the outflow tubes and measurement through thin glass panes, were taken into account. The system was evaluated in a clinical study.

The proposed sensor is part of a comprehensive mobile CBI monitoring system. Although the overall system is outside the scope of this paper, the main reasons for its developments will be briefly presented. As mentioned earlier, continuous monitoring of CBI is important. In this regard, a single nurse can take care of about two to four patients, provided that they are in the same room and the nurse has no other duties. In reality, however, there are easily twenty concurrent CBI on a urology ward, spread across multiple rooms and cared for by a single nurse. Therefore, continuous monitoring of all patients is not feasible. In addition, the nurse must maintain the CBIs by emptying or changing the waste bags and replacing empty saline bags in a timely manner. Finally, the nurse also has to manage the daily routine of attending to the well-being of patients, dispensing medications, conferring with physicians, and much more. Regardless of what a nurse is doing, he/she always has the pressure in the background that he/she should be monitoring CBIs. This has several negative consequences. Nurses are under constant stress, even when monitoring a particular CBI because other CBIs cannot be monitored at the same time. Nevertheless, acute bleeding in one or more patients can occur at any time. On the other hand, maintenance of the CBIs also places a burden on the nursing staff since overflowing waste bags as well as empty reservoir bags should be avoided as much as possible. The monitoring system relieves the nursing staff of this pressure since the flow rate and *Hb*/blood concentration are constantly monitored by a technical system. The flow and concentration measurements can be used to effectively calculate whether irrigation is being performed optimally. As soon as one of the parameters leaves its predefined limits, alarms can be transmitted to a mobile device and/or to the ward room. In addition, the current status of all CBIs can be visualized simultaneously, giving the nurse an optimal overview of all patients. The use of the monitoring system allows nurses to fully concentrate on their respective tasks, while giving them the security of knowing that they will be informed in time if intervention is required. Overall, the system thus provides the opportunity to improve CBI monitoring and thus patient care, while at the same time reducing the workload of nurses, thus indirectly further improving patient care. The proposed sensor is an extremely important and integral component of such a monitoring system. Additional fields of application are conceivable, but will not be considered further at this point.

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

#### *2.1. Selecting the Measurement Method*

A requirement for the sensor development was the seamless integration with existing CBI systems. A refractometer-based approach or the use of other sensors that require direct contact with the medium to measure *Hb* concentration is, therefore, inappropriate. The *Hb* measurement in this work is performed based on light absorption. A known amount of light is exerted by a LED, passed through a tube filled with CBI waste fluid and captured on the opposite side. According to Beer's Law, the concentration is then related to the intensities and measurement section via the following:

$$I\_t = I\_0 \cdot 10^{-\varepsilon \Delta z c} \text{ }$$

where *It* and *I*0 denote the transmitted and incident light intensities, respectively, *ε* and *c* denote the molar absorptivity and concentration of the light absorptive compound, and ∆*z* is the optical path length through the drainage fluid. By selecting a proper power supply and wavelength for the LED, we obtain constant *I*0 and *ε*. When the outflow tube is full of drainage fluid, ∆*z* is also constant and determined by the diameter of the tube. Therefore, we can infer the *Hb* concentration *c* solely based on *It* .

#### *2.2. Selecting the Light Source*

In this simple form, Beer's Law is only valid for a single light absorbing material. The waste fluid is, however, composed of various materials, i.e., irrigation fluid and potentially additional medications, as well as urine and blood. The latter two materials are compounds on their own. We, therefore, investigated how strong the influence of the additional ingredients would be. In a work from Pegau et al. [16], we found that the influence of salt on the absorption property of water is very small in the near infrared frequency band and that salt might even lower absorptivity. From Palmer et al. [17], we learned that the absorption of water compared to *Hb* is negligible (approx. ratio 1:40,000). Additionally, urine coloring molecules, such as betanins, have their main absorption in the range of 400–600 nm and absorb only very little light in the near infrared as Gonçalves et al. [18] reported. In summary, we can state that *Hb* absorption is absolutely dominant, particularly in the near infrared but also in the visible spectral band. Hence, we apply Beer's Law in its simple form without danger of significantly degrading the measurement accuracy.

The absorption curve of *Hb* and *HbO*2 is shown in Figure 1. For the formation of blood clots, it is basically irrelevant whether *Hb* is oxygenated or not. In order to account for both *Hb* versions at the same time, isobestic points, i.e., points at which two chemical species have the same molar absorptivity, should be used for measurement. For *Hb* and *HbO*2 the main isobestic points are at frequencies of 420, 545, 570, and 800 nm. In tests with blood samples and the targeted tubing, it was found that the absorption for the lower frequency points is way too high for a CBI monitoring system. The transmitted amount of light quickly drops below values that can be reliably distinguished from background noise. Additionally, total absorption is reached very quickly, so increasing the light intensity would not help. Combined with knowledge of the absorption properties of other relevant materials, these findings led to the determination of the central measurement frequency at 800 nm.

**Figure 1.** Molar extinction coefficient e in [ L mol·cm ] of *Hb* and *HbO*2 for wavelengths between 250 and 1000 nm. Major isobestic points are located at 420, 545, 570 and 800 nm.

For an LED, the emitted light intensity as well as the light frequency depend on the applied current. Especially when using batteries in a mobile setting, care must be taken to provide a constant current supply. To this end, we implemented a small electronic component (see Figure 2), that guarantees a constant current over a wide voltage range. Implementation was done on a stripboard, using regular size electronic components.

**Figure 2.** LED needs a constant current to emit a well-defined light. (**Left**): Circuit diagram of the used electronic. (**Right**): Test results with and without load. Constant current can be guaranteed for a wide range of input voltages.

#### *2.3. Selecting the Light Sensor*

In order to measure the transmitted light, a photo resistor (LDR), or a photo diode (PD) would theoretically be sufficient, as shown in [13]. Taking practical considerations into account, i.e., the fact that the CBI tubes do have reinforcing ribs, a single measurement is not sufficient for a reliable monitoring system. Even more problematic is the fact that bubbles are traveling unpredictably in the outflow tube. Gain ribs appear darker in the detector image and suggest increased *Hb* values, while the presence of bubbles dramatically decreases light absorption and suggests lower *Hb* values. While slightly increased *Hb* readings might be tolerable, significantly lower readings would put patients at risk. As a consequence, measuring the *Hb* concentration with a single point measurement (LDR/PD) is very dangerous in practice. We propose to use a camera chip instead. At only moderately increased cost, a wealth of measurements can be performed, enabling the application of advanced image processing techniques.

For system demonstration, we chose a monochrome camera without an IR filter, which provides us with 10 bit linear intensities. While the concrete camera model is irrelevant, the three mentioned properties of the camera are important. A monochrome camera lacks a Bayer filter; hence, all pixels generate readings of the same quality. The lack of IR filtering allows for reliable measurements in the near infrared frequency band. The 10 bit color depth increases the dynamic range and provides an increased measurement range.

#### *2.4. Extending the Measurement Interval*

Considering the exponential relationship between the intensity and *Hb* concentration, the dynamic range of even a 10 bit sensor is rather limited. Starting from an exposure *ELow* that just saturates the sensor for pure irrigation fluid, the maximum measurable *Hb* concentration is less than 0.6 g/dL (equivalent to ≈4% blood). Increasing the exposure time by a factor of 8 (*EHigh* ) allows measurements of more than 3.2 mg/dL *Hb* (≈20% blood), but leads to overexposure for the *Hb* concentration below 0.1 g/dL. We, therefore, introduced a dual-exposure setup, which is controlled, using a hysteresis threshold *TLow*,*High* . We chose *TLow* (i.e., where we switch from *EHigh* to *ELow*) at 0.2 g/dL *Hb* and *THigh* at 0.6 g/dL *Hb*. In terms of sensor readings, i.e., intensities, the *TLow* is at 800 and the *THigh* is at 80. The overall measurement range of the proposed sensor, therefore, covers the interval [0–3.23] g/dL hemoglobin concentration or [0–21.5]% blood concentration.

#### *2.5. Processing of Sensor Images*

We apply three different image processing methods, namely the detection of shadow artifacts caused by the gain ribs, detection of bubbles, and automatic exposure adjustment.

Figure 3 gives an impression of the effect created by the gain ribs. The left subfigure shows a small, opened segment of the tube placed on a white paper with parallel pencil lines. Significant distortion of the lines can be observed near the ribs, which act like planoconvex lenses. The subfigure on the right shows the effect of the ribs on the detector image. The horizontal position of the artifacts depends on how the tube was inserted into the sensor. Vertical line artifacts with mainly decreased intensity can be observed. In some cases, even a slight intensity increase can occur, depending on the positioning of a rib within the optical pass. Although the positioning of the tube cannot be controlled in practice, the artifacts can be reliably detected. Since the tube does not change its position during the entire CBI process, it is practical to automatically select an artifact-free measurement ROI when setting up the system. This can be realized by inspecting the intensity profile of a horizontal line and comparing it to a known intensity profile. An exemplary profile is depicted in Figure 4. The blue regions are those that are not directly illuminated, while the yellow regions are transitional zones that are partially lit. Both regions are a constant size and determined during the sensor assembly since they solely depend on the actual camera and LED placement. The union of the red and green regions is the potential measurement area. In the example, we detected two artifact regions, marked in red, and an artifact-free region, where we can safely place the measurement ROI.

**Figure 3.** (**Left**): A segment of the waste fluid tube. Reinforcement ribs are running along the tube. The ribs significantly change the optical pass of light traveling through the tube. (**Right**): Acquired sensor image with marked artifacts caused by reinforcements. The horizontal location of the shadows depends on the tube placement and cannot be controlled. We note that the right figure was intensity adjusted for better depiction, which amplified intensity variations in the vertical direction as well as the noise level.

**Figure 4.** Exemplary intensity profile segmented into regions: (blue) background region, (yellow) transitional region, (red) artifact region, (green) measurement region.

As already mentioned earlier, bubbles pose a significant problem to the measurement process. While the gain ribs create static artifacts that only slightly and systematically alter the measurements, bubbles are dynamic by nature and severely impact the measurement. Figure 5 demonstrates the effect of a bubble in the sensor image. Please note that for all image sequences, the same specimen was used, i.e., a small tube section filled with diluted blood and sealed with a small fraction of air, similar to those depicted in Figure 9. Striking is the difference in the image sequences depicted in Figure 5. The only difference between the sequences is the location with respect to the camera and LED, where the bubble travels along the tube. Already, these examples demonstrate that a bubble detection based on classic approaches is not feasible. The situation becomes even more challenging when several small bubbles form a foam-like cluster. Furthermore, we found that also

pure image processing of a single image is not sufficient to reliably detect bubbles: in the case that a bubble of a reasonable size remains static within the light pass, the images can be indistinguishable from images without bubbles at a lower *Hb* concentration. It is, therefore, imperative to use a system that can analyze image sequences. To this end, we trained a convolutional recurrent neural network comprised of convolution layers to extract optical features followed by long short-term memory (LSTM) layers to account for the temporal aspect. The network was trained on video sequences acquired with the sensor hardware. To this end, test specimen with various blood concentrations and bubble sizes were prepared (see Section 2.8). The network was trained fully supervised. The ground truth was generated by manually marking the images in the sequences where a bubble just entered the view or almost completely left the view. The remaining images were classified automatically. Figure 6 depicts the network architecture. For details on the network and how it was trained, we refer the interested reader to [19].

**Figure 5.** Examples of the appearance of bubbles in the sensor image. (**Upper row**): bubble passing the sensor close to the LED. (**Middle row**): bubble passing the sensor close to the camera. (**Lower row**): bubble passing near the center of the tube. For all images, the same specimen was used, so *Hb* concentration as well as bubble size are the same.

**Figure 6.** Network to detect bubbles in the optical path. The convolutional layers extract image features, such as bright and dark edges, bright spots, etc., while the LSTM layers account for the temporal components. The network is able to detect bubbles traversing in both directions and even accounts for those that remain static in the optical path for an extended time.

#### *2.6. Sensor Signal Processing*

Sensor signal processing, i.e., bubble detection, *Hb* value generation, and communication, was implemented on a Raspberry Pi 4B. The camera was connected via the USB 3 interface, which also provided power to the LED. Furthermore, the Raspberry Pi was used to communicate the data to the main system and to control additional sensors that are out of the scope of this paper. The operating system used was a completely stripped-down Ubuntu Linux with a custom compiled kernel.

The bubble detection network was run on the ARM-CPU exploiting the ONNX Runtime [20]. The camera provides images at a frame rate of 25 fps. All images are analyzed by the network. Sequences of 8 images are combined to guarantee stable outputs, even in the case that the network fails to detect a bubble or falsely detects a bubble in an image. The output frequency of the sensor was set to 10 measurements per second. Data were provided

to the user side (C# DLL) by means of a double buffer. Overall, power consumption for the running system is below 700 mA.

To further support future integration of the sensor, we quantized the network to fixed point representation (Figure 7), using a quantization aware training (QAT) method. The resulting network was implemented on FPGA (Xilinx xc7z020-clg400-1). As shown in Table 1, the standard version of the network still leaves sufficient resources for additional hardware, such as a MIPI camera interface as well as the calibration tables to convert sensor readings into *Hb* values. Thus, the whole signal processing can be performed on a single SoC of a size less than 20 × 20 × 2 mm. This not only reduces the overall sensor size, but also significantly lowers the overall power consumption and, in turn, enables mobile CBI monitoring.

**Figure 7.** Schematic representation of the sensor. On the left side, an LED exerts light that is transmitted through the tube filled with irrigation waste fluid. The transmitted light is captured by a camera. The bright squares next to the tube represent thin glass plates that can be inserted into the optical pass to seal the LED and the camera and to enable effective cleaning and disinfection in clinical application.

**Table 1.** Listing of hardware resources used for bubble detection and MIPI interface on the FPGA. Two different implementations of the bubble detection network are given: one that complies with the current specifications regarding frame rates as well as a maximum performance version. Please note that the standard version leaves sufficient space to realize a MIPI interface to directly connect the camera.

#### *2.7. Sensor Housing*

A sensor used in clinical practice needs to be cleaned and disinfected. Neither the camera lens nor the LED are well suited for this process. Besides potential deterioration effects, camera lenses feature sharp edges where bacteria and other germs might survive. It is therefore mandatory to seal the optical components. A particularly suited material for the case at hand appears to be glass. We conducted a series of tests using glass cover slips for microscopic samples with a thickness of 0.15 mm ± 0.02 mm. Assuming circular openings for the camera lens ( ≈ 1.0 cm) and the LED light channel ( ≈ 0.5 cm), the cover slips appear to be sufficiently robust.

Since the positions of the glass plates, camera, and LED are rigid and only light frequencies in a very small band are used, we confirmed in tests that the influence of glass is a constant offset to the measurement. This was expected by altering Beer's Law in the following way:

$$I\_t - \Delta I\_t = (I\_0 - \Delta I\_0) \cdot 10^{-\varepsilon \Delta z c} \text{ \AA}$$

where ∆*It* and ∆*I*0 denote the attenuation of the exerted and captured intensity by the two glass slips.

The housing itself was 3D printed from black acrylonitrile butadiene styrene (ABS). Tests with white ABS showed significant light penetration in the sensor, which greatly affected the measurements. Although it would be possible to shield the measurement path by other means, using black material seemed to be the most direct approach.

In contrast to [13], we decided for a rigid housing instead of a clamp. While a clamp has the advantage that it can be most easily attached to a variety of tubes, it is much harder to know the true length of the measurement distance. Furthermore, the clamp pressure deforms the tube and thus introduces uncertainty. Lastly, it is much harder to avoid ambient illumination to impact on the measurement. An image of the sensor is given in Figure 8. Three versions of the sensor were assembled and provided for integration.

**Figure 8.** Image of the assembled sensor as used in the clinical validation study. Communication and power supply were realized using a USB 3 interface. A short tube segment was inserted to showcase sensor attachment to the CBI system.

#### *2.8. Preparation of Test Specimen*

Figure 9 depicts a set of test specimen used for sensor development. A complete set always covered a concentration range between 0% and 22% in 2% increments, i.e., 12 samples. Each specimen was manufactured in two different versions: one with saline irrigation fluid and one with an electrolyte-free irrigation fluid. The latter is mainly used during surgery, but opened bags are used up afterwards.

**Figure 9.** Test specimen samples. Tube segments were filled with diluted blood and sealed with hot glue. In order to test the effect of bubbles, a small fraction of air was left inside. This also prevented direct contact between the hot glue and test liquid. Samples were produced for blood concentration between 0% and 22% in increments of 2%. *Hb* concentration was measured using a hemoglobin analyzer.

To produce the test specimens, the drain tube of an irrigation system was cut into pieces approximately 6 cm long. One end of the tube segments was sealed with hot glue. Percentage indicators were written on each tube, using waterproof ink. This was done in a preparatory phase to avoid any delay when working with the blood.

The blood used for development was fresh porcine blood obtained immediately after slaughter. A total of 6 mL diluted Heparin (250 IU/mL) was added to 200 mL blood to prevent clotting (coagulation), and the mixture was carefully stirred to avoid sedimentation/separation. Since the proposed sensor actually measures the *Hb* concentration, the ground truth *Hb* value was measured using a hemoglobin analyzer (Measuring range: 0–25.5 g/dL; Imprecision (within run): CV < 1%; Calibrated against HiCN reference method.). Then, the blood was diluted with the respective irrigation fluid. The prepared liquid was poured into the tube segment, leaving approximately 0.5–1.0 cm of air, and the tube was sealed with hot glue. Leaving a small air bubble in the tube served two purposes: first, the presence of an air bubble allowed the effect of air bubbles to be studied accurately during scanning. Second, the air prevented possible adverse effects from the hot glue coming in direct contact with the liquid.

Samples prepared in this way provided valid measurements for one working day. Repeated measurements the following day showed significant deterioration, even when the samples were stored in the refrigerator.

In order to assess the effect of urine, additional test specimens were prepared. Human urine, collected immediately after waking up, was used for this purpose. With this urine, we were able to ensure a maximal content of salts and other metabolites. Tests on pure urine samples and samples composed of blood, urine and irrigation fluid revealed no significant effect on the measurements.

#### **3. Results**

In a planned self-experiment, the sensor was tested using blood from three human donors. Using human blood donations did not require ethics approval, according to the local ethics council. The full laboratory setup of the CBI monitoring system is shown in Figure 10 (left). The liquid flowed from the upper bag via a Foley catheter through a 3D printed bladder model into the urine bag. A pump controlled the rate of inflow. The *Hb* sensor was mounted to the outflow tube to detect the *Hb* level for data collection. To obtain five different blood concentrations, we prepared five bags filled with 500 mL saline solution and varied amounts of blood. As a reference measurement, the *Hb* concentration of each bag was measured via blood gas analysis (BGA). The average difference of the sensordetected *Hb* level from the BGA-detected *Hb* level was 0.29 g/dL. The *Hb* concentration in human blood is 15 g/dL on average, so the measurement error with regard to the full measurement interval was estimated as ~10%. We note that the full measurement range of the sensor was tested in the laboratory trials. Since measurement accuracy drops significantly at higher concentrations, an increased mean deviation from the BGA was to be expected. Table 2 gives a rough estimate of the measurement accuracy for some blood concentrations. The highest blood concentration that can be quantified by the sensor is roughly 20%, i.e., 3 g/dL *Hb* concentration.

Additionally, we employed random *Hb* concentrations to model the rise and fall in *Hb* level to simulate a real-life complication scenario, presented as a graph of the *Hb* level over time as depicted in Figure 10 (right). To imitate acute bleeding, blood was injected into the 3D-printed bladder model. The blood supply was replaced with a saline solution to resemble stopped bleeding. In the experiments, the sensor time constants for the CBI system, with the sensor placed near the waste bag, were 33 s for rise and 119 s for fall. Accordingly, the sensor's response times were 2.7 min for rise and 9.9 min for fall. However, we would like to point out that these values strongly depend on temperature, irrigation flow, emissivity of the material and size, i.e., length and diameter of the respective irrigation tube, including the catheter. Furthermore, bubble detection was tested. The sensor reliably recognized random air bubbles or artificially generated air bubbles in the outflow tube.

**Figure 10.** (**Left**): full system setup for clinical laboratory trials using blood from human donors: (1) blood donor, (2) display, (3) 3D bladder model, (4) *Hb* sensor, (5) urine bag, and (6) pump. The sensor was used in an extended setup, which is not in the scope of this paper. (**Right**): two exemplary measurement curves for increasing and decreasing *Hb* concentrations. Sensor time constants were estimated as 33 s for rise and 119 s for fall.

**Table 2.** Measurement errors with respect to the measured blood concentration. With increasing blood concentration the measurement accuracy drops due to the non-linearity of Beer's Law. We note that these values were acquired during development and were not part of the clinical study.

After receiving a positive ethical vote, a clinical trial involving 20 patients was conducted to test the sensor and other components of the CBI monitoring system. All patients were male (mean age: 73.1 years) and underwent TURB. CBI was administered as part of the standard treatment procedure. Further patient information was not collected to keep the data protection footprint low. For each patient, *Hb* levels were monitored for at least three hours, using the CBI monitoring system. As a reference measurement, we regularly checked the *Hb* level in the urine bag via BGA. We would like to point out that, unlike the *Hb* sensor, our BGA device did not allow to measure the *Hb* level continuously. The *Hb* level detected by BGA is the average of several minutes. Figure 11 depicts an example of the sensor-detected *Hb* level's development illustrated as a red line. As gray data points, ten BGA-detected *Hb* levels are presented. The median *Hb* values are at a comparable level. The mean deviation of the sensor-detected *Hb* level from that of the BGA-detected *Hb* level across all patients was −0.003 g/dL. Between minute 120 and 140, we noted a rise in *Hb* levels. The reason was an incorrectly configured inflow speed of saline solution. As the inflow speed was accelerated, the *Hb* level decreased to normal. Thus the sensor reliably measured *Hb* levels within a clinically acceptable deviation and with a high responsiveness across all experiments. In the patient study, the functionality of the *Hb* sensor's air bubble detection could also be positively confirmed.

**Figure 11.** Measurement curve for a real patient. The blood gas analysis agrees well with the sensor readings. At approximately minute 125, the irrigation was inadvertently left off after the urine bag was emptied, resulting in a real-life example of the importance of the continuous monitoring of CBI.

#### **4. Discussion**

As shown in the results, we were able to design a sensor to measure the *Hb* concentration in the effluent of a CBI with decent accuracy (±0.003 g/dL). While accuracy is highest for very small concentrations, it drops for higher concentrations, due to the exponential relation between concentration and absorption. However, for practical application this does not have negative implications. Already at 10% blood concentration, medical personnel should consider intervention. In case of a low irrigation flow, the flow should be increased to prevent cloth formation. In case of an already high irrigation flow, acute bleeding should be assumed. To increase measurement accuracy, we extended the measurement interval, using multiple exposure times. Further extension using even higher exposure times does not improve the situation, since total absorption occurs—given the fixed tube diameter.

We decided on a rigid sensor housing to achieve the most reliable measurements for system demonstration. This, of course, comes at the drawback that only a particular tube diameter is supported by the sensor. However, insets can be easily designed to account for various tube diameters. Vastly larger tube diameters decrease the measurement interval significantly since total absorption occurs earlier, while smaller tube diameters increase the measurement interval. In any case, an individual re-calibration of the sensor will be necessary.

The effect of salt on the measurement cannot be completely ignored. As can be found in the work of Pegau et al. [16], salt lowers the absorption of near infrared light by water and thus, can lead to an underestimation of the *Hb* concentration. As a consequence, it is definitely necessary to calibrate the sensor for the respective irrigation fluid. However, the amount of salt introduced by urine is negligible and does not have a significant effect on the measurement accuracy.

As stated in Section 2.8, we prepared the blood samples for the development from porcine blood. Although it is very similar to human blood, the calibration tables do not perfectly fit the practical application. A re-calibration with human blood was not conducted prior to the study since the developers had no access to human blood samples. A re-calibration at the hospital site was not possible since the procedure is technically involved and developers could not travel due to COVID-19 restrictions. Still, the achieved results during the clinical study were more than satisfying.

For system demonstration, we used a cheap camera without housing and other accessories. This camera was not designed for medical applications but rather for rapid prototyping. As such, it is relatively large (40 × 40 × 38 mm) and the main reason for the current overall sensor size of 44 × 44 × 85 mm. Additionally, the constant current supply electronic is relatively large (40 × 23 × 12 mm). Using a miniaturized camera and an integrated electronic, the overall size of the sensor can easily be reduced to less than 15 × 15 × 30 mm in total. Realizing a SoC version, the sensor including evaluation electronic could be realized in a package not larger than 25 × 25 × 30 mm.

Additionally, an approach similar to that of Hageman et al. [11] can be extended into a non-invasive measurement solution, such as the one presented by Ding et al. [12]. However, such a solution is inferior to the absorption-based approach for several reasons. Firstly, in the absorption approach, the emitter and sensor can be placed in line with the

test sample. This allows for a technically less-involved setup, because all geometries can be assumed to be planar and of a constant size, without significant sacrifice of accuracy. Any additional material, such as the tube or glass slides (c.f. Sections 2.5 and 2.7) only play a minor role when considering Beer's Law.

In contrast, a reflectivity-based system must consider a technically involved angled arrangement because the light emitter and sensor cannot be in the same position. More importantly, the tube has a significant, non-trivial impact on the reflectance characteristics. In particular, the coloring of the tube affects the perceived color, regardless of whether a color scale is used by personnel or a computerized evaluation is performed. Lastly, reflection spectra appear to be less specific than absorption spectra.

A completely different method for estimating blood concentration in the effluent could be developed based on the photoacoustic effect. Briefly, when a medium is irradiated with a sequence of light flashes, periodic heating (and cooling) occurs. The resulting alternation of volume expansion and contraction constitutes a source of sound [21]. Using laser sources with specific wavelengths and, for example, piezoelectric acoustic receivers, it is possible to measure *Hb* concentration with high precision and specificity. In [22] a system for in vitro measurement of *Hb* and *HbO*2 is described. Integrating the laser emitter into a catheter, a similar system could even allow concentration measurements directly in the bladder, which would significantly reduce the response time of the sensor (cf. Section 3). On the other hand, a rather expensive optical fiber would have to be integrated into the catheter, which is a disposable product. In [23], an interesting work is presented that aims at miniaturizing a photoacoustic sensing system. However, we assume that current solutions are still relatively large, complex and expensive and often require much more energy, compared to the sensor presented here.

Lastly, the clinical study only evaluated sensor functionality, i.e., the quality of the *Hb* measurement as well as bubble detection during clinical application. Further processing of sensor readings to, for example, generate alarms or treatment guidance, was neither performed nor part of the study. Sensor readings were acquired blindly and evaluated retrospectively.

#### **5. Conclusions**

Although further investigations are required as part of an approval study in accordance with the German as well as the international Medical Devices Act, the basic practical suitability of the sensor was demonstrated. Its non-invasive applicability allows a simple extension of the existing CBI systems and thus, a gapless monitoring of the irrigation procedure. Measuring the *Hb* concentration using the absorption principle (Beer's Law) provides reliable data for the application at hand. In combination with a (mobile) communication system, monitoring data can be immediately transferred to nurse guard rooms or even mobile devices. This significantly reduces the burden on nursing staff to care for an entire ward, while still ensuring effective CBI for the benefit of patients.

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

**Funding:** This research was partially funded by the German Ministry of Research and Education grant number BMBF VISIMON 16SV7861K.

**Institutional Review Board Statement:** Positive ethics vote was achieved from Ethics Committee of Albert-Ludwigs-Universität Freiburg, Freiburg im Breisgau, Germany, 20-1177 MPG §23b, 29 October 2020. The study was registered according to article 35 of the Declaration of Helsinki 2013 at Deutsches Register Klinischer Studien (DRKS) under ID DRKS00023647, 27 January 2021.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors want to thank the VISIMON project partner DITABIS AG, Pforzheim, Germany for the integration of the sensor into a system fit for practical application in a patient study.

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

#### **References**

## *Communication* **A Cell's Viscoelasticity Measurement Method Based on the Spheroidization Process of Non-Spherical Shaped Cell**

**Yaowei Liu 1 , Yujie Zhang 1 , Maosheng Cui 2 , Xiangfei Zhao 1 , Mingzhu Sun 1 and Xin Zhao 1, \***

**Abstract:** The mechanical properties of biological cells, especially the elastic modulus and viscosity of cells, have been identified to reflect cell viability and cell states. The existing measuring techniques need additional equipment or operation condition. This paper presents a cell's viscoelasticity measurement method based on the spheroidization process of non-spherical shaped cell. The viscoelasticity of porcine fetal fibroblast was measured. Firstly, we introduced the process of recording the spheroidization process of porcine fetal fibroblast. Secondly, we built the viscoelastic model for simulating a cell's spheroidization process. Then, we simulated the spheroidization process of porcine fetal fibroblast and got the simulated spheroidization process. By identifying the parameters in the viscoelastic model, we got the elasticity (500 Pa) and viscosity (10 Pa·s) of porcine fetal fibroblast. The results showed that the magnitude of the elasticity and viscosity were in agreement with those measured by traditional method. To verify the accuracy of the proposed method, we imitated the spheroidization process with silicone oil, a kind of viscous and uniform liquid with determined viscosity. We did the silicone oil's spheroidization experiment and simulated this process. The simulation results also fitted the experimental results well.

**Keywords:** robotic cell manipulation; mechanical properties; elasticity measurement; viscosity measurement; cell mechanics

## **1. Introduction**

The mechanical properties of biological cells, especially the elastic modulus and viscosity of cells, can provide an important basis for the evaluation of cell viability and cell states and the judgment of biological activity [1–3], and is crucial for the understanding of cell structure and physiological function [4–6]. To measure cell viscoelasticity, scientists have developed methods such as atomic force microscopy (AFM) [7–12], magnetic tweezers technique [13,14], optical tweezers technique [15,16], microfluidic technique [17–19], and micropipette aspiration (MA) technique [20–24]. These methods are suitable for different situations. AFM technique detects the viscoelasticity of the cell by moving the cantilever probe in vertical direction and monitoring its bending displacement. Magnetic tweezers technique and optical tweezers technique apply a certain force to the magnetic beads or silicon beads adhered to cells through magnetic field or light field to deform cells and obtain the viscoelasticity of cells. Microfluidic technique obtains the viscoelasticity of cells by detecting the deformation of cells under different microchannels and different shear forces. MA technique obtains the viscoelasticity of cells by measuring the length of the cells aspirated into the micropipette under different pressures.

Among these techniques, MA technique has become widely used due to the reasons of no need to purchase or prepare additional equipment, lower measurement cost, and easier integration into existing commercial micro-operation systems [25]. However, the

**Citation:** Liu, Y.; Zhang, Y.; Cui, M.; Zhao, X.; Sun, M.; Zhao, X. A Cell's Viscoelasticity Measurement Method Based on the Spheroidization Process of Non-Spherical Shaped Cell. *Sensors* **2021**, *21*, 5561. https://doi.org/ 10.3390/s21165561

Academic Editors: Vahid Abolghasemi, Hossein Anisi and Saideh Ferdowsi

Received: 26 July 2021 Accepted: 16 August 2021 Published: 18 August 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/).

micropipette aspiration method has high requirements for the seal between the cell and the micropipette in the measurement process [26]. Slightly improper sealing will result in ineffective MA operations, which will have a great impact on the measurement results. Meanwhile, the measuring of results is highly dependent on the accuracy of the force sensor. The viscoelasticity differences of the same cells measured by the same research groups using the micropipette aspiration method will also be very large. For example, the elasticity of human chondrocytes measured by Jones et al. was 0.65 ± 0.63 kPa [27], wherein the standard deviation was as large as the measured value. As the shape of the cell might be non-spherical, it will be more difficult to seal the cell and micropipette. In order to eliminate the influence of sealing on the measurement results, it is necessary to design a cell viscoelastic measurement method based on the micropipette aspiration platform and with low requirements for sealing.

In this paper, we proposed a cell's viscoelasticity measurement method based on the spheroidization process of non-spherical shaped cell. The spheroidization process means the process of some deformable non-spherical objects turning into spherical shapes due to surface tension. We firstly introduced the method of recording the spheroidization process of porcine fetal fibroblast and recorded the fetal fibroblast's spheroidization process. Secondly, we built the viscoelastic model for simulating non-spherical shaped cell's spheroidization process based on the fact that the capsule-like porcine fetal fibroblast will finally become spherical (Figure 1). Then, we simulated the spheroidization process of porcine fetal fibroblast and got the simulated spheroidization process. By changing the parameters in the simulations, we got the elasticity and viscosity that best fitted the experiments. The magnitude of the elasticity and viscosity of fetal fibroblast was in agreement with those measured in other literatures. To verify the accuracy of this method, we imitated the spheroidization process with silicone oil, a kind of viscous and uniform liquid with determined viscosity. We did the silicone oil's spheroidization experiment and simulated this process. The simulation results fitted the experimental results well.

**Figure 1.** Spheroidization process of the capsule-like porcine fetal fibroblast. (**a**–**d**) From capsule-like cell to spherical-like cell.

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

#### *2.1. System Setup*

The spheroidization experiment of porcine fetal fibroblast was performed on the self-developed NK-MR601 micro-operation system [28–30] (Figure 2). The system consists a microscope (CK-40, Olympus, Tokyo, Japan); a CCD camera (W-V-460, Panasonic, μ

Osaka, Japan, frame rate: 20 frame/s); a motorized X-Y stage (travel range: 100 mm, repeatability: ±1 µm/s, maximum speed: 2 mm/s); two X-Y-Z manipulators (travel range: 50 mm, repeatability: ±1 µm/s, maximum speed: 1 mm/s); a self-developed micro-injector, providing negative pressure to aspirate the fetal fibroblast and positive pressure to eject the fetal fibroblast; a self-developed motion control box, controlling the micro-platform, micro-manipulators, and micro-injector through the host computer. μ μ

μ

**Figure 2.** NK-MR601 micro-operation system.

The silicone oil spheroidization experiment was performed on NK-MR601 with the CCD replaced by a highspeed camera (C110, Miro, Wayne, NJ, USA, frame rate: 1000 frame/s).

μ The micropipettes used in the spheroidization experiments of porcine fetal fibroblast and silicone oil were made from borosilicate glass tubes with an outer diameter of 1 mm and an inner diameter of 0.8 mm. The micropipette used in the fetal fibroblast spheroidization experiments were pulled by the puller (MODEL P-97, Sutter Instrument, Novato, CA, USA), and fractured by the microforge (MF-900, NARISHIGE, Tokyo, Japan) with an inner diameter of 10 µm. The micropipette used in the silicone oil spheroidization experiments was pulled and fractured by hand by Yaowei Liu, with an outer diameter of 200 µm.

#### *2.2. Preparation and Spheroidization of Porcine Fetal Fibroblast*

We obtained the porcine fetal fibroblast from a sow at day 35 of pregnancy. After removal of head, internal organs and limbs, the remaining parts were cut into pieces at approximately 1 mm3 . We smeared the pieces evenly in a 35 mm dish and cultured in Dulbecco's modified Eagle's medium (DMEM), containing 15% fetal calf serum (FCS), 0.1 mM non-essential amino acids (NEAA), 6 µL/mL Gentamycin and 0.05 mM L-glutamine. Cells were cultured in a 37 ◦C humidified incubator containing 5% CO2. Cells were trypsinized and cryo-preserved for use when cells grown to ~90% confluence.

The spheroidization experiments of porcine fetal fibroblast were carried out in Medium 199 (Sigma). Figure 3 shows the typical images of the porcine fetal fibroblast spheroidization process:

and measured with 2 frames per minute. The initial ratio was determined by the inner diameter of micropipette and the cell volume in the experiment. The method of detecting the size of capsule-like fetal fibroblast is described in Appendix A.

**Figure 3.** Typical images in the spheroidization process of porcine fetal fibroblast. (**a**–**d**) From porcine fetal fibroblast just coming out of the micropipette to porcine fetal fibroblast becoming a sphere.

#### *2.3. Spheroidization of Silicone Oil*

Figure 4 shows the method of recording the silicone oil spheroidization process:

**Figure 4.** (**a**) Drop M199 into a petri dish and overlay M199 with silicone oil; (**b**) aspirate silicone oil into the micropipette; (**c**) move the micropipette tip into M199 and eject silicone oil.

#### *2.4. Viscoelastic Model*

We use a viscoelastic model to study the spheroidization process. The cell is modeled as homogeneous viscoelastic liquid, which is surrounded by infinitesimal thin cortical layer. We use the Jeffrey's viscoelastic fluid model (Equation (2)) because it is independent of the frame of reference and the motion as a whole in space [31]. Besides, it has only 2 additional parameters, while being able to imitate the viscoelastic behavior. More complex models (e.g., heterogeneous liquid) are hard to modify the parameters to obtain reliable results. In the simulation, the cortical layer is realized by surface tension. We made the following assumptions: ,

$$\overset{\nabla}{T} = \frac{\partial T}{\partial T} + (\boldsymbol{\sigma} \cdot \boldsymbol{\nabla})\boldsymbol{T} - \boldsymbol{\nabla}\boldsymbol{\sigma}\cdot\boldsymbol{T} - \boldsymbol{T}\cdot\left(\boldsymbol{\nabla}\boldsymbol{v}\right)^{\mathrm{T}} \tag{1}$$

( ) ()

$$
\lambda \stackrel{\text{v}}{T} + T = 2\eta E \tag{2}
$$

T

$$
\rho \frac{D\mathbf{v}}{Dt} = \nabla \cdot (-p\mathbf{I} + \mathbf{K} + \mathbf{T})\tag{3}
$$

$$\mathbf{K} = 2\mu \mathbf{E} \tag{4}$$

where ∇ *T* denotes the upper convection derivative [32] of *T* defined by Equation (1). *T* is the viscoelastic stress tensor that changes with time according to Equation (2). *v* is the velocity field. *λ* is the characteristic time. *η* is the viscosity in the viscoelastic term. *E* is the strain-rate tensor. *ρ* is the density of porcine fetal fibroblast and is assumed to be constant in the following simulation. *p* is the pressure. *Dv*/*Dt* is the material derivative of *v*. **I** is the unit tensor. *K* is the shear stress tensor which can be obtained from Equation (4). *µ* is another viscosity in the pure viscous term. Equation (3) is the Navier-Stokes equation. *λ η ρ μ*

An analogy of the model in one dimension can be illustrated as Figure 5. The total stress tensor (right terms in the bracket of Equation (3)) is composed of hydrostatic pressure −*p***I**, viscous stress *K* and viscoelastic stress *T*. In this figure, *E* is the stiffness coefficient of spring, *η* and *µ* are the viscosities of two dashpots, and *λ* = *η*/*E*. The bottom line represents a Maxwell model, for which the relationship of strain rate *e* and stress *σ* is *λ* . *σ* + *σ* = *ηe* in 1D case. By replacing the time derivative with upper convected derivative and extend the equation to 3D tensor form, we get Equation (2). The usage of upper convected derivative for continuum materials was argued by Oldroyd in [31]. − ܧ/ߟ = ߣ *μ η σ* ߪߣሶ + ݁ߟ = ߪ

**Figure 5.** Illustration of viscoelastic model.

#### *2.5. Simulation of the Spheroidization Process*

We used Ansys Student Fluent software to simulate the spheroidization processes. Two-dimensional axisymmetric was adopted for efficiency. μ μ

For the simulation of fibroblast, we set the capsule-like porcine fetal fibroblast as a cylinder in the middle and two hemispheres at both ends, as shown in Figure 6 (only 1/4 part was used by applying the axisymmetric and symmetric condition). In this paper, the length was set as 20.8 µm, the width was set as 8.6 µm. The computation domain was a rectangle of 16 × 10 µm2 which was divided into 0.1 × 0.1 µm2 structural quadrilateral grids. The pressure variance in the scale of several micrometers is of the order of 0.01 Pa, which is far less than the barometric pressure. Besides, it is balanced by the gravity, so we neglected both pressure variance and the gravity. We used volume of fluid (VOF) model to introduce the surface tension. Laminar flow was adopted because of low Reynolds number. We set the four boundaries as axisymmetric, symmetric and pressure outlet, respectively. The densities were set as 1080 kg/m3 [25] for fibroblast and 998.2 kg/m3 (the density of water at 20 ◦C) for surrounding liquid. The surface tension coefficient *T* was set as 10 µN/m [33]. The viscosity *η* and elasticity *E* (*E* = *η*/*λ*) were introduced with user-defined scalars (UDS, see Appendix B). Ansys Student Fluent software solves the momentum Equation (3) without viscoelastic stress term *T* by default. We used user-defined scalars (UDS) to insert *T* into the equation (see Appendix B for more details). To study the influence of viscosity and elasticity, we firstly set *λ* = 1 s, and changed viscosity *η* as 10 Pa·s, 20 Pa·s, 50 Pa·s, 100 Pa·s, 200 Pa·s, and 500 Pa·s. Secondly, we set *η* = 500 Pa·s, and changed *λ* as 500 s, 100 s, 20 s, 1 s, 0 s, 1 s and 0.02 s [33]. The timestep was 3 s. Based on the results obtained when the parameters selected in a wide range, we made more compact selections and compared the results with experimental data. The one that fitted best was viewed as measurement result. μ μ ° μ *η η λ λ η η λ*

**Figure 6.** Capsule-like porcine fetal fibroblast.

μ μ *μ* For the simulation of silicone oil, the initial shape was set as axisymmetric while the contour being obtained by image processing procedure (see Appendix A). The computation domain was a rectangle of 500 × 200 µm2 which was divided into 2 × 2 µm2 structural quadrilateral grids. As is shown in Figure 7, one boundary is symmetric axis and others are pressure outlet. Using the contour obtained by image processing (Appendix A), a user-defined function sets the corresponding region as silicone oil (secondary phase), while the remainder as culture medium (primary phase). The viscosity *µ* was inserted by setting the material property of silicone oil in the software. The viscoelastic stress term *T* was removed because it was considered as pure viscous liquid. Volume of fluid model and laminar flow was adopted. Then we run the simulation with 1 ms timestep. Please see the simulation procedure details in Appendix B.

**Figure 7.** Geometric shape and boundary settings.

#### **3. Results**

*3.1. Spheroidization Result of Porcine Fetal Fibroblast and Its Simulation*

The porcine fetal fibroblast was used in the experiments.

The typical images in the spheroidization process of porcine fetal fibroblast have been shown in Figures 1 and 2 (Video S1). The length and the width changing process was shown in Figure 8. The whole spheroidization process took 15 min to reach a 90% width–length ratio.

**Figure 8.** Variation of porcine fetal fibroblast length and width with time.

*λ λ η η* We got the simulated width–length ratio at the condition that *λ* = 1 and changed viscosity *η*, as shown in Figure 9. Video S2 shows the simulated spheroidization process of porcine fetal fibroblast. The results show that a larger *η* will prevent the porcine fetal fibroblast from turning into a sphere, and the spheroidization time becomes longer.

*η*

*η*

*λ η* **Figure 9.** The simulated spheroidization process when *λ* = 1 and viscosity *η* changed.

*η λ λ* We got the simulated width–length ratio at the condition that *η* = 500 Pa·s and changed *λ*, as shown in Figure 10. The experimental results showed that the spheroidization process was more intense in the initial stage, but because the elasticity was smaller when *λ* was larger, the small elasticity will bring a lag effect in the later stage of spheroidization, which would make the later stage of spheroidization slow down significantly. *η λ λ*

*η λ* **Figure 10.** The simulated spheroidization process when *η* = 500 Pa·s and *λ* changed.

*η λ η λ η* Finally, by changing the values of *η* and *λ* in the simulation experiment, different curves of the spheroidization process of the simulated porcine fetal fibroblast were obtained. By comparing with the curves of the spheroidization process obtained in the real experiment, the elasticity and viscosity could be obtained. Figure 11 shows the length and width variation of porcine fetal fibroblast with time in the experiment and simulation. The viscosity *η* obtained in this experiment is 10 Pa·s and the elasticity *E* is 500 Pa. The magnitude of the results was in agreement with the measured results in [29].

*η λ*

*η*

**Figure 11.** Variation of porcine fetal fibroblast length and width with time in the experiment and simulation.

#### *3.2. Spheroidization Result of Silicone Oil and Its Simulation*

Figure 12 shows the typical images in the spheroidization process of silicone oil. Video S3 shows this process of slowing down 100 times. Video S4 shows the simulated spheroidization process of silicone oil.

**Figure 12.** Typical images in the spheroidization process of silicone oil. (**a**–**d**) From silicone oil just coming out of the micropipette to silicone oil becoming a sphere.

We defined the time span from the release to 95% width–length ratio as spheroidization time, which is denoted as *ts* . The *ts* was 80 ms in the experiment.

The viscosity and density of the silicone oil (Sigma-Aldrich) at 25 ◦C was 9.71 Pa·s and 0.971 g/mL. To simulate the spheroidization process of silicone oil, we also need to know the surface tension coefficient with culture medium of silicone oil. We measured the surface tension coefficient between silicone oil and culture medium by Du Noüy ring method [26], the coefficient was 0.024 N/m. The detail is shown in Appendix C. The method of detecting the contour of the silicone oil is shown in Appendix A.

Figure 13 shows the simulation results of silicone oil with different viscosities and surface tension coefficient. The R-square values of the fitted curves are 0.99 and 0.98. The results revealed that spheroidization time increases linearly as viscosity and the reciprocal of surface tension coefficient increases. The results showed that the spheroidization time *ts* was 100 ms in the simulation, which was similar to the real experiment (80 ms).

**Figure 13.** Simulation results: (**a**) spheroidization time increases linearly as viscosity increases. (**b**) Spheroidization time increases linearly as the surface tension coefficient increases.

#### **4. Discussion**

We should know that the results measured in this paper were based on the bulk measurements, by which the cells were assumed as isotropic, homogeneous. However, in reality cells are very heterogeneous and contain organelles. We also need to measure the local force and dissipative gradients, as well as map them across the cell surface [34–36]. Considering the measuring efficiency, only two parameters are necessary to describe the cellular mechanics, so the bulk measurement is more appropriate.

We used the cells in the suspension state instead of adherent state in this paper. Because our method needs to aspirate to the whole cell into the micropipette, and it is difficult to aspirate the adherent cells into the micropipette because of the adhesion. Since the whole suspension cell was sucked into the micropipette, the cell spheroidization process was only related to the shape of the micropipette. The seal between the cell and the micropipette will not affect the spheroidization recording results, which can avoid the influence of seal in the micropipette aspiration method.

We performed the cell experiments three times and the simulations 12 times per cell. As the cells were collected from one batch, the experimental curves were very similar. The parameters in the simulation were set not very accurate (just integers), so the results of these three cells were the same. The measurement results may be significantly different among different cell types or different cell batches. Our future work will focus on measuring the viscoelasticity of different cell types and improving the simulation accuracy by adjusting the parameters more accurately.

We could see that the simulation results shown in Figure 10 do not overlap with the experiment results shown in Figure 7 exactly. We supposed that there were three reasons:

(1) The whole spheroidization process would take a tremendously long time in a in vitro environment, which would influence the viscoelasticity of the cell a lot. We only recorded the spheroidization process when the cells reached a 90% width–length ratio, which took about 15 min. Meanwhile, the simulation process recorded the whole spheroidization process. So, the experiment results and the simulation results could not overlap exactly.

(2) As mentioned above, the parameters in the simulation were not set very accurately (just integers), so the simulation results could not exactly fit the experiment results.

(3) The initial velocity was set to zero. This may have caused an initial acceleration stage from zero, while the non-zero initial shrink in velocity was found from the experiment. The problem was handled by running several steps in advance.

We know that if the viscous term of a viscoelastic body is increased, then it takes longer to get back to its original shape. Our results showed that a larger *η* will prevent the porcine fetal fibroblast from turning into a sphere, which could further verify the validity of our simulation results.

#### **5. Conclusions**

This paper presents a cell's viscoelasticity measurement method based on the spheroidization process of a non-spherical shaped cell. We firstly introduced the process of recording the spheroidization process of porcine fetal fibroblast. We secondly built the viscoelastic model for simulating a cell's spheroidization process. We simulated the spheroidization process of porcine fetal fibroblast and got the simulated spheroidization process. Then we got the elasticity (500 Pa) and viscosity (10 Pa·s) of porcine fetal fibroblast by identifying the parameters in the viscoelastic model. The results showed that the magnitude of the elasticity and viscosity were in agreement with those measured by a traditional method. To verify the accuracy of the proposed method, we imitated the spheroidization process with silicone oil, a kind of viscous and uniform liquid with determined viscosity. We did the silicone oil's spheroidization experiment and simulated this process. The simulation results also fitted the experimental results well.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/s21165561/s1, Video S1: The spheroidization process of porcine fetal fibroblast (0 to 15 s: slowing down 5 times; 16 to 23 s speeding up 100 times). Video S2: The simulated spheroidization process of porcine fetal fibroblast. Video S3: The spheroidization process of silicone oil. Video S4: The simulated spheroidization process of silicone oil.

**Author Contributions:** X.Z. (Xin Zhao) and Y.L. conceived the idea for the study. X.Z. (Xin Zhao), Y.L., Y.Z. and M.S. designed the experiments. Y.L., Y.Z., M.C. and X.Z. (Xiangfei Zhao) performed the experiments. X.Z. (Xin Zhao) supervised the project. Y.L. and Y.Z. wrote the manuscript. Y.L., Y.Z. and M.C. contribute equally to this paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was jointly supported by the National Key R&D Program of China (2018YFB1304905, 2020YFB1313101), the National Natural Science Foundation of China (U1813210, 62003174, 62003173), and the China Postdoctoral Science Foundation (2020M680865).

**Institutional Review Board Statement:** All the procedures were approved by the Animal Care and Use Committee of Tianjin Animal Science and Veterinary Research Institute, and were performed in accordance with the NIH Guide for the Care and Use of Laboratory Animals (No. 8023, revised in 1996).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

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

#### **Appendix A. Size Measurements of Porcine Fetal Fibroblast and Silicone Oil**

In the initialization step of the simulation, we chose the image where the porcine fetal fibroblast and silicone oil were just out of the micropipette, then extracted its contour. The flow chart of image processing is shown in Figure A1. After the processes, we got a list of contour points where x axis was the axisymmetric axis. As the contour was symmetric, we used the upper half of it to conduct an axisymmetric simulation, which saved time and computational resources, and was more convenient to show the results. For the consecutive images, length and width were measured automatically for silicone oil and manually for fibroblast because of the low contrast ratio.

**Figure A1.** Image processing: (**a**) ROI of grayscale image; (**b**) binarized image; (**c**) axisymmetric contour; (**d**) upper side of contour; (**e**) flow chart of image processing.

#### **Appendix B. Simulation Procedure**

We used Ansys Student Fluent to simulate the spheroidization process. The Fluent software provides the volume of fluid (VOF) module for multiphase simulation such that we can simulate the surface tension between cells and culture medium. Because the velocities in such processes are small (of the order of 1 nm/s), the Reynolds number is small, so the laminar model was adopted. For the cell simulation, we used user-defined scalars (UDS) to introduce the viscoelastic stress term. Details of the simulation are described below.

μ μ **Geometry:** rectangle region of 16 × 10 µm2 for fibroblast and 500 × 200 µm2 for silicone oil;

μ μ **Mesh:** structural quadrilateral grids of 0.1 × 0.1 µm2 for fibroblast and 2 × 2 µm2 for silicone oil;

#### **Boundary conditions:**

left and top—pressure outlet with 0 Pascal gauge pressure;

right—symmetric for fibroblast and pressure outlet for silicone oil;

down—axisymmetric;

#### **Fluid filed settings:**

General: 2D axisymmetric;

Models panel: multiphase—VOF (Phase Interactions—Surface Tension), Viscous (Laminar); **Materials**: The density and viscosity are set in this stage. In this paper, the values are: 1080 kg/m3 of density and 500 Pa s of viscosity for fibroblast, 9.71 kg/m3 of density and 9.71 Pa s of viscosity for silicone oil, 998.2 kg/m3 (the density of water at 20 ◦C) of density and 1.003 × 10−3 Pa s for surrounding liquid.

#### **Initialization and UDS equations:**

− Ansys Student Fluent provides user-defined functions (UDF) for customizing fluid simulation. In this paper, a C code file utilizing predefined macros sets the initial fibroblast and silicone oil area, as well as introduces equations about viscoelastic stress *T*. To add *T* into the momentum equation, four scalars were defined and added to the momentum equation by DEFINE\_SOURCE macro. The usage of UDF can be found in official ANSYS Fluent UDF Manual. The initial velocities were set as zero.

For an arbitrary user-defined scalar *φ*, Fluent solves the Equation (A1). ( )

$$\frac{\partial \rho \phi}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i} (\rho u\_i \phi - \Gamma\_k \frac{\partial \phi}{\partial \mathbf{x}\_i}) = \mathbf{S} \tag{A1}$$

where *ρ* is the density, *u* is the velocity, Γ is the diffusion coefficient, *S* is the source term. The terms in the equation represent unsteady term, convective flux, diffusion and source from left to right. Expanding Equation (2) and using Einstein summation convention, we get Equation (A2):

$$
\lambda \frac{\partial \pi\_{lj}}{\partial t} + \lambda u\_k \frac{\partial \pi\_{lj}}{\partial \mathbf{x}\_k} - \lambda \pi\_{ik} \frac{\partial u\_j}{\partial \mathbf{x}\_k} - \lambda \pi\_{kj} \frac{\partial u\_i}{\partial \mathbf{x}\_k} + \pi\_{ij} = \mu (\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}) \tag{A2}
$$

The first term was added with DEFINE\_UDS\_UNSTEADY macro. The second term was added with DEFINE\_UDS\_FLUX macro. Other terms were added with DE-FINE\_SOURCE macro. It should be noticed that for 2D axisymmetric setup in this paper, the gradients along normal direction of the geometry plane were all zero. Therefore, only four scalars *τ*11, *τ*12, *τ*21, *τ*22 were defined. Besides, a factor 2π should be multiplied for volume and area because of axisymmetric condition according to the manual. 11 12 21 22 ,,, π

**Solver**: 1 ms of timestep for silicone oil and 3 s for fibroblast. For each time step, iterate at most 100 steps with 0.001 as convergence absolute criteria.

#### **Notes:**

1. The quadrilateral grids are suggested for the meshing as it is more stable than triangular grids in the simulation. We suppose that it was because triangular grids cause larger curvature, which makes surface tension change dramatically within a local area.

2. The mesh size should not be too small, not only for the computation efficiency, but also for the convergence. When reducing the mesh size down to a certain scale, the simulation gets hard to converge. It is also supposed to be the result of large local surface tension.

#### **Appendix C. Surface Tension Coefficient Measurement of Silicone Oil**

As is illustrated in Figure A2, when the interface between silicone oil and culture medium was flat and the platinum ring was on the interface, we have

$$F\_0 + B = G \tag{A3}$$

where *F*0 is the external force in this case, *B* is the buoyancy, *G* is the total weight of platinum ring and its frame.

**Figure A2.** Du Noüy ring method to measure surface tension coefficient.

When we pull the platinum ring up, it stuck the interface, causing a circle of interface protruding (Figure A2), and an extreme thin film formed near the ring. Inner and outer interface both contributed to the surface tension, i.e.,:

$$F\_{\text{total\\_tension}} = 2\pi d\gamma\tag{A4}$$

where *γ* was the surface tension coefficient. The pull is a quasi-static process, we can get

$$F\_m + B = G + F\_{total\\_tenion} \tag{A5}$$

where *Fm* was the maximum external force in the process. The change of *B* was negligible compared with change of external force, so we used the same notation *B*. Combining (A3), (A4) and (A5), we get:

$$\gamma = \frac{F\_m - F\_0}{2\pi d} = \frac{\Delta F}{2\pi d} \tag{A6}$$

In this experiment, the diameter *d* was 14 mm, and we got that ∆*F* was 2.11 mN in average, so *γ* was 0.024 N/m.

#### **References**

## *Article* **DMAS Beamforming with Complementary Subset Transmit for Ultrasound Coherence-Based Power Doppler Detection in Multi-Angle Plane-Wave Imaging**

**Che-Chou Shen \* and Yen-Chen Chu**

Department of Electrical Engineering, National Taiwan University of Science and Technology, Taipei 106335, Taiwan<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>Tel.: +886-2-27301229; Fax: +886-2-27376699

**Abstract:** Conventional ultrasonic coherent plane-wave (PW) compounding corresponds to Delayand-Sum (DAS) beamforming of low-resolution images from distinct PW transmit angles. Nonetheless, the trade-off between the level of clutter artifacts and the number of PW transmit angle may compromise the image quality in ultrafast acquisition. Delay-Multiply-and-Sum (DMAS) beamforming in the dimension of PW transmit angle is capable of suppressing clutter interference and is readily compatible with the conventional method. In DMAS, a tunable *p* value is used to modulate the signal coherence estimated from the low-resolution images to produce the final high-resolution output and does not require huge memory allocation to record all the received channel data in multi-angle PW imaging. In this study, DMAS beamforming is used to construct a novel coherence-based power Doppler detection together with the complementary subset transmit (CST) technique to further reduce the noise level. For *p* = 2.0 as an example, simulation results indicate that the DMAS beamforming alone can improve the Doppler SNR by 8.2 dB compared to DAS counterpart. Another 6-dB increase in Doppler SNR can be further obtained when the CST technique is combined with DMAS beamforming with sufficient ensemble averaging. The CST technique can also be performed with DAS beamforming, though the improvement in Doppler SNR and CNR is relatively minor. Experimental results also agree with the simulations. Nonetheless, since the DMAS beamforming involves multiplicative operation, clutter filtering in the ensemble direction has to be performed on the low-resolution images before DMAS to remove the stationary tissue without coupling from the flow signal.

**Keywords:** delay-and-sum (DAS); delay-multiply-and-sum (DMAS); signal coherence; power doppler detection; plane-wave (PW) imaging; complementary subset transmit (CST); coherent plane-wave compounding (CPWC)

#### **1. Background**

Delay-and-Sum (DAS) beamforming is routinely adopted to produce image output in medical ultrasound imaging by compensating the time delay of the received echo according to the geometric path of propagation before coherent summation [1]. However, it suffers from intrinsic limitations such as insufficient image resolution and noticeable off-axis clutter. For plane-wave (PW) imaging which depends on the unfocused transmit wave to illuminate a wide field-of-view [2], these limitations are notably evident. In singleangle PW imaging, the received backscattered echoes from one PW transmit event are processed using DAS beamforming in the direction of receiving channel to generate the corresponding low-resolution image at frame rate on the order of kHz. Therefore, PW imaging is also referred as ultrafast imaging. The image quality of PW imaging can be improved by coherent plane wave compounding (CPWC) in which synthetic transmit focusing is achieved from multi-angle PW transmit [3,4]. Specifically, low-resolution images are firstly acquired from several PW transmit angles and then coherently combined

**Citation:** Shen, C.-C.; Chu, Y.-C. DMAS Beamforming with Complementary Subset Transmit for Ultrasound Coherence-Based Power Doppler Detection in Multi-Angle Plane-Wave Imaging. *Sensors* **2021**, *21*, 4856. https://doi.org/10.3390/ s21144856

Academic Editors: Vahid Abolghasemi, Hossein Anisi, Saideh Ferdowsi and Manuel Graña

Received: 2 June 2021 Accepted: 13 July 2021 Published: 16 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/).

to achieve the final high-resolution CPWC image. Considering an imaging depth of 75 mm, the pulse-repetition-interval (PRI) for pulse-echo imaging will be 100 µs with the sound velocity of 1.5 mm/µs. This corresponds to a frame rate of 10 kHz for low-resolution PW imaging. Assuming that every 10 low-resolution images is coherently combined to improve the image quality, the resultant high-resolution CPWC images will be produced at a frame rate of 1 kHz. With this high frame rate, CPWC imaging has been utilized to detect the motion of imaged objects in transient elastography [4,5] and Doppler flow imaging [6–13]. However, it should be noted that the image quality in CPWC imaging relies on the number of low-resolution images involved in the compounding, and thus an inevitable trade-off between the image quality and the frame rate exists. In other words, with just a few PW transmit angles, the image quality in CPWC imaging is generally unsatisfactory due to the presence of high clutter artifacts.