MPEPE, a predictive approach to improve protein expression in E. coli based on deep learning

Graphical abstract

a b s t r a c t
The expression of proteins in Escherichia coli is often essential for their characterization, modification, and subsequent application. Gene sequence is the major factor contributing expression. In this study, we used the expression data from 6438 heterologous proteins under the same expression condition in E. coli to construct a deep learning classifier for screening high-and low-expression proteins. In conjunction with conserved residue analysis to minimize functional disruption, a mutation predictor for enhanced protein expression (MPEPE) was proposed to identify mutations conducive to protein expression. MPEPE identified mutation sites in laccase 13B22 and the glucose dehydrogenase FAD-AtGDH, that significantly increased both expression levels and activity of these proteins. Additionally, a significant correlation of 0.46 between the predicted high level expression propensity with the constructed models and the protein abundance of endogenous genes in E. coli was also been detected. Therefore, the study provides foundational insights into the relationship between specific amino acid usage, codon usage, and protein expression, and is essential for research and industrial applications.

Introduction
High-level production of soluble recombinant proteins at low cost is crucial for their application in many fields. However, there are barriers to heterologous expression in many hosts, including in the unicellular bacterial Escherichia coli expression system. In order to minimize the effect of heterologous recombinant proteins on cell growth, a smart expression host could degrade the heterologous proteins by the proteolysis system or form an inclusion body. It is estimated that <50% of bacterial and <15% of nonbacterial proteins could be expressed in soluble form in E. coli [1]. Studies of expression system modification have examined, for example, the optimization of expression conditions [2][3][4][5], expression of co-soluble tags fused to target proteins [6], as well as co-expression of the molecular chaperone [7]. While these methods work for a small number of proteins, the level of improvement for most proteins is still limited and some strategies to increase expression reduce the catalytic activity.
Additional studies show that there are many different strategies to increase the soluble and functional expression of foreign proteins. According to Deng et al., alanine-or leucine-scanning mutagenesis increased soluble expression, and leucine could increase the protein helices and their stability against degradation by proteinase K [8]. In recent years, the well-established directed evolution approach has been employed to optimize the coding sequence based on soluble expression phage-assisted continuous evolution (SE-PACE). This method has been used to evolve some antibody fragments and maltose-binding protein (MBP) to increase their expression [9]. Michal Jamroz et al. conducted a promising rational design study in which they proposed the AGGRESCAN 3D structural aggregation predictor. This method was used to modify the green fluorescent protein and the human single-domain VH antibody, effectively reducing the aggregation propensity of the protein and increasing its expression [10]. In summary, those results indicated that the sequence is also a vital factor to regulate its expression.
One promising approach for uncovering the key factors affecting high-level and soluble expression is deep learning, a subfield of machine learning that uses multi-layered Deep Neural Networks (DNNs) to extract novel features from input data [11]. The first predictive tool for protein solubility based on amino acid sequence was proposed by Wilkinson and Harrison in 1991 [12]. Based on 81 soluble or insoluble protein sequences, the authors found that the charge average and turn-forming residue fraction correlated with inclusion body formation. In recent years, deep learning algorithms for enhancing the expression of exogenous proteins in E. coli have been developed. Approaches based on logistic regression [13] and deep learning architecture (EPSOL) [14] were proposed to learn the features related to soluble expression from the expressed sequence fragments and construct the classifier to distinguish the soluble expressed proteins, but they lacked the experimental validation. A model based on the bidirectional long-short-term memory conditional random field was constructed to optimize the codon and enhance protein expression [15]. A feed-forward artificial neural networks were also constructed based on the ribosome profiling [16]. All of the results indicated that the gene sequence could affect the protein expression. A complicated model should be trained based on the real gene sequence and design the gene to improve its expression.
In this study, we propose a predictive model, mutation predictor for enhanced protein expression (MPEPE). The 6438 proteins that were experimentally validated the expression yields in E. coli [17] were selected to train and validate the prediction model based on multi-layered deep neural networks (DNNs). The evolutionary method was incorporated into the model to virtually screen the mutant sites that might make positive contributions to protein expression but not disrupt its function. When the strategy was applied on two enzyme proteins, the laccase 13B22 and the flavin adenine dinucleotide-dependent glucose dehydrogenase (FAD-AtGDH), the expression and activity of these two enzymes in E. coli were significantly increased. This study will help researchers to understand the relationship between amino acids and soluble protein expression and is important for the industrial application of enzymatic proteins.

Collection of the protein expression dataset
We collected the protein expression dataset from a published study [17]. All of the 6438 proteins in the dataset have been classified to six classes (Class1, Class2, Class3, Class4, Class5, and Class6) based on the protein expression level under the identical conditions and from the same promoter in E. coli in the reference [17]. The protein level was scored on the integer scale from class1 (lowest) to class6 (highest) based on the visual inspection of whole cell lysates in Coomassie-blue-stained SDS-PAGE gels [17][18]. The proteins in Class1 have the lowest protein levels while the proteins in Class6 have the highest protein levels. For the aim to construct the relative balanced training dataset in terms of the size, the low expression dataset consisted of 2308 proteins in Class1, Class2, and Class3, while the 1973 proteins in Class6 comprised the high expression dataset. Therefore, the ratio of low expression proteins to high expression proteins was $1:1. The independent validation dataset was composed of the proteins in Class4 and Class5 and its size was 2067. Table 1 shows detailed information on the size of each dataset. The protein abundance data were downloaded from the paxdb database and the weighted average of WHOLE_ORGAN-ISM between the different protein abundance values was taken (http://pax-db.org/) [19].

Data processing of protein expression datasets and coding schemes
In order to use synonymous codon number, the specific amino acid, and specific nucleotide combination to construct three Deep Neural Network (DNN) models, the nucleotide (codon) sequences in the protein expression dataset were translated into amino acid sequences and codon number sequences, the detailed coding information was shown in Table S1.
The negative and positive samples are randomly partitioned into k (k = 10) equal size subsamples. Of the k subsamples, a single subsample is retained as the validation data for testing the model, and the remaining k-1 subsamples are used as training data. The 10cross-validation process is then repeated k(k = 10) times (the folds), with each of the k subsamples used exactly once as the validation data. Based on this method, all samples are used for both training and validation, and each sample is used for validation exactly once. The detailed architectures and optimized hyperparameters showed in Fig. S1 and three DNN models were trained on the TensorFlow platform [20] based on Keras 2.1.5 in the Python 2.7.15 programming environment.
Three models were successfully constructed. The architecture of each one started with an embedding layer, followed by two convolutional, maximum pooling layers, an LSTM layer, and a final prediction layer. In addition to the prediction layer, each layer was followed again by a Batch Normalization layer and a final dropout layer. In this model, the rectified linear unit (ReLU) activation function was used (except for the final prediction layer) as follows: x; else where x denotes the feature map from the convolution operation (the weighted sum of a neuron). The softmax activation function was used in the final prediction layer.
To avoid over-training and the resultant over-fitting of the prediction model, we used an early stopping technique to detect the prediction accuracy achieved the high score on the test dataset in the training process. We optimized the various hyper-parameters in the DNN architecture, including the number of layers, number of kernels, kernel size, fully connected layer size, dropout rate, learning rate, batch size, activation functions, number of nodes, and optimizers using the optimization package HYPERAS (https:// github.com/maxpumperla/hyperas). The final parameters used for the prediction models are shown in Table S2. We also verified the performance of the proposed method by a 10-fold crossvalidation method. The code for the above DNN architecture with Python v2.7.15 environments is available from GitHub (https:// github.com/BRITian/MPEPE).

Evaluation of prediction performance
A 10-fold cross-validation strategy was used to train and evaluate the performance of the prediction models. The performance of the prediction models was evaluated with metrics including accuracy, recall, precision, F1-score, and the area under the receiver operating characteristic curve (AUROC) score, which were each calculated based on 10-fold cross-validation. The area under ROC (AUROC) and the area under PRC (AUPRC) of the prediction models were greater than 0.5, which indicated that the performance of the constructed models was better than that of a random classifier. The metrics were calculated using the Keras package as follows: where TP, TN, FP, and FN indicate true positive, true negative, false positive, and false negative, respectively.

The PSAP and entropy analysis
Using BLASTP, we identified and downloaded sequences of homologous proteins of laccase 13B22 and the glucose dehydrogenase FAD-AtGDH (identity >40% and coverage >50%), from NCBI (https://www.ncbi.nlm.nih.gov/). The homologous protein sequences were aligned with their respective wild type and spliced to the same length. Then the position-specific amino-acid probability (PSAP) matrices for both enzymes were calculated to determine the conservation of amino acids at each locus. The entropy of each residue was then determined using the PSAP calculation tool [4,21,22]. In this study, we used the cutoffs of PSAP value of mutations greater than 0.1 and the entropy value of the residue greater than 0.65 to select the non-conserved residues and natural selected mutations.

Simulation of protein three-dimensional structure
The AlphaFold 2.0 tool was used to simulate the threedimensional structure of laccase 13B22 and FAD-AtGDH from its amino acid sequence with the publicly available code and default parameters [23]. All five structures were simulated and the structure with the least structural energy was selected for the following analysis.

Construction of mutants
We generated single-point mutations in the residues that differed from the sequence of the wild type in the mutants using the two-step PCR mutagenesis strategy [24]. Mutant primers were designed using Oligo software ver.7.0 (Table S3) and were then used to generate the single-point mutants. Using wild-type plasmid DNA template, we amplified the mutant site by T7-F primer and each downstream primer containing a mutation sequence or T7-R and each upstream primer containing a mutation sequence under the action of Phanta Max Super-Fidelity DNA Polymerase (Vazyme, Nanjing, China). PCR was performed as follows: 95°C for 5 min for the preheat, 95°C for 30 s, 58°C for 30 s, 72°C for 1 min, 32 cycles, and 72°C for 10 min. We used an AxyPrep DNA gel extraction kit (Axygene, California, USA) to recover the target gene fragment, and then it was used as the primer and the wildtype plasmid as the template for the second round of PCR amplification. The mutant plasmid was obtained by PCR amplification. The second round PCR was performed as follows: 95°C, 5 min for the preheat, 95°C for 30 s, 72°C for 5 min, 32 cycles, and 72°C for 10 min. The wild-type plasmid was eliminated by DpnI (NEB, Ipswich, United Kingdom) and the mutant plasmids were recovered using a purification and recovery kit (TransGen, Beijing, China) and then transformed into E. coli Top10 competent cells (TransGen, Beijing, China) using standard procedures [25]. Through bacterial liquid PCR, the single clone with the correct size of the target gene was preliminarily identified and sequenced by TSINGKE biological technology (Beijing, China), and then the correct mutant plasmid was transformed into competent E. coli BL21 (DE3) cells (TransGen, Beijing, China). Multi-site mutations were obtained by adding two mutation sites in each round according to the above method.

Detection of the expression of laccase 13B22
E. coli BL21 (DE3) cells harboring the recombinant plasmid were cultured in a 50 mL LB (Lysogeny Broth) medium supplemented with kanamycin (50 lg/mL) at 30°C with shaking at 120 rpm to an OD 600 of 0.7$0.75. Then, 0.2 mM isopropyl-b-d-thiogalactopyra noside (IPTG) and 0.4 mM CuCl 2 were added to the culture medium, and the temperature was reduced to 25°C. Incubation was continued for a further 4 h, during which microaerobic conditions were achieved by switching off the shaking function [26,27]. Cells were harvested after a further 20 h of growth by centrifugation at 8000g for 10 min [28]. The pellets were resuspended in 5 mL buffer containing 20 mM Tris-HCl (pH 8.0). The cells were disrupted by sonication on ice, and debris was removed by centrifugation at 4°C and 8000g for 30 min. After transferring the crushed supernatant to a new pre-cooled 10 mL EP tube and resuspending the crushed pellet in 5 mL of 20 mM Tris-Cl buffer (pH 8.0), the laccase activity was determined. 50 mL protein were taken concurrently to prepare SDS-PAGE samples.
SDS-PAGE was performed using standard procedures and gels were stained with Coomassie brilliant blue R250. For Western blotting analysis, SDS-PAGE gels were transferred to polyvinylidene fluoride (PVDF) membranes (Amersham, Piscataway, NJ, USA), and these were blocked with 5% non-fat milk (Applygen, Beijing, China) and incubated with mouse anti-His monoclonal antibody (TransGen, Beijing, China), followed by horseradish peroxidase HRP-conjugated goat anti-mouse IgG (TransGen, Beijing, China). Proteins were visualized using a BeyoECL Plus Chemiluminescence Detection Kit (Beyotime, Shanghai, China) and Chemiluminescence Touch Imaging System (e-BLOT, Shanghai, China).

The detection of the activity of the laccase 13B22
Laccase activity was assayed at 37°C using 2,2-azino-di-(3-et hylbenzthiazoline-sulfonate) ABTS (Sigma-Aldrich, St. Louis, USA) as substrate. An assay mixture containing 200 mL 5 mM ABTS and 750 mL 50 mM citrate/phosphate was preheated at 37°C for 2 min, and then 50 mL crushed supernatant was added to react accurately for 3 min. This mixture was then terminated in an icewater bath for 1 min. The increase in absorbance due to the oxidation of ABTS at 420 nm was measured (e420 = 36,000 M À1 cm À1 ). One unit was defined as the amount of enzyme that oxidized 1 lmol of substrate per minute [29].
2.10. Detection of the expression of glucose dehydrogenase activity of FAD-AtGDH E. coli BL21 (DE3) cells harboring the recombinant plasmid were cultured in a 50 mL LB medium supplemented with kanamycin (50 lg/mL) at 37°C with shaking at 200 rpm to an OD600 of 0.7$0.75. Then 0.3 mM isopropyl-b-d-thiogalactopyranoside was added to the culture medium, which was grown for 18 to 20 h in LB medium at 16°C and 200 rpm to obtain FAD-AtGDH as soluble and active proteins. Cells were harvested by centrifugation at 8000g for 10 min. The pellets were resuspended in 5 mL buffer containing 20 mM Tris-HCl (pH 8.0). The cells were disrupted by sonication on ice, and debris was removed by centrifugation at 4°C and 8000g for 30 min. The crushed supernatant was then transferred to a new pre-cooled 10 mL EP tube, and the crushed pellets resuspended in 5 mL of 20 mM Tris-Cl buffer (pH 8.0), and the glucose dehydrogenase activity was determined. At the same time, 50 mL protein were taken to prepare SDS-PAGE samples.
2.11. Detection of the activity of glucose dehydrogenase activity of FAD-AtGDH FAD-AtGDH activity was assayed spectrophotometrically using 2,6-dichloroindophenol (DCIP, e600 = 16.3/mM/cm) (Solarbio, Beijing, China) and phenazine methosulfate (PMS) (Solarbio, Beijing, China) as electron acceptors. The reaction was followed for per min measured and the reaction continued for 5 min at 600 nm using the SpectraMax M2 microplate reader (Molecular Devices, Silicon Valley, USA). The DCIP-based assay contained 50 mM Potassium phosphate buffer (pH 6.5), 0.06 mM DCIP, 0.6 mM PMS and 100 mM D-glucose. One unit of FAD-AtGDH activity was defined as the amount of enzyme required for the reduction of 1 lmol glucose or electron acceptor per min under the assay condition [30].

Overview of the MPEPE prediction strategy
In order to predict mutations that could enhance the heterologous expression of proteins in the soluble fraction in E. coli, we proposed the MPEPE strategy for rational optimization of gene sequence related to translation rate and protein surface charge. To this end, we constructed DNNs, the predicted accuracy of which was higher than the other machine learning methods (Fig. S2, Table S4), to select mutations with a high potential to increase protein expression levels based on published, experimentallydetermined expression data generated in E. coli under the identical expression condition and expression system including the promoter ( Fig. 1A; and more details in the next sub-section) [17]. In addition, evolutionary analysis was employed to select nonconserved residues for mutagenesis that appeared unlikely to disrupt protein function (Fig. 1B). Based on the above two design components ( Fig. 1A and B), we screened single-point mutations that might positively contribute to expression levels, but not disrupt the catalytic function of the protein of interest (Fig. 1C). We experimentally verified the contribution of each screened mutation by SDS-PAGE and enzymatic activity analysis, then generated mutants with stacked mutations to determine the combined effects of these confirmed positive mutations (Fig. 1D). Through this discovery pipeline we obtained several candidate mutations for exploration of the determining factors related to increased protein expression in the soluble fraction in E. coli (Fig. 1E) and to subsequently design a mutant with significantly enhanced expression in E. coli.

Protein expression data-based selection of charged amino acids and rapidly translated codons in E. coli
To investigate relationships between protein expression level and amino acid charge properties or synonymous codon number for specific amino acids, we identified preferentially selected amino acids/codons in 6438 published proteins that were experimentally classified into six categories based on their heterologous protein expression under identical conditions and promoter in E. coli, collected from diverse phylogenetic sources and provided broad sampling of codon and amino acid space due to variations in codon-usage frequency in the source organisms [17]. The expression indexes for the six classes ranged from 1 to 6 (from lowest to highest), and consisted of n = 1754 proteins in group 1, n = 131 in group 2, n = 423 in group 3, n = 896 in group 4, n = 1171 in group 5, and n = 1973 in group 6. Due to substantial bias in the number of proteins allocated among the six classes, groups 1, 2, and 3 were combined to represent low expression proteins, while group 6 represented highly expressed proteins. The resulting dataset was relatively balanced, with a low:high ratio of 1:0.85. Proteins in groups 4 and 5 were selected as independent categories to validate the performance of the prediction model.
An initial analysis of the data revealed significantly different usage of amino acids between the low-and high-expression proteins ( Fig. 2A). Among the eleven amino acids (E, K, D, Q, H, T, M, N, Y, V, and F) preferred in the highly expressed proteins (Fig. 2B), four were charged amino acids (E, K, D, and H) and eight amino acids (E, K, D, Q, H, N, Y, and F) were encoded by only two codons. Additionally, the amino acid E, K and D were also preferred in the endogenous proteins of E. coli (Fig. S3). However, among the other nine amino acids (R, L, A, P, G, S, C, W, and I) preferentially found in low-expression proteins (Fig. 2C), only one was a charged amino acid (R), and six of these amino acids (R, L, A, P, G, and S) were encoded by four or six codons.
In addition, five of the nine amino acids preferentially used in low-expression proteins (R, S, G, C and P) contained the codons AGG, AGT, GGT, TGT, and CCC which have been established to significantly affect in vivo translation speed in E. coli (Loss of attenuation >3) [31]. Here, the loss of attenuation of the codon was relative score that was the b-galactosidase activity of each codon construct divided by the activity of the wild-type construct that contains the CAC His codon [31]. The high loss of attenuation of the codon indicated that the codon could decrease the translation speed, and the vice versa [31]. Notably, in high-expression proteins, the codons encoded any of the preferred eleven amino acids that did not greatly affect the translation rate [31]. The average attenuation loss of preferential high-expression codons was 1.57 ± 0.53, significantly (p < 0.01) lower than the average attenuation loss (2.75 ± 2.12) of preferred codons in the low expression proteins [31] (Table S5). Therefore, the average translation rate of high expression proteins was significantly greater than that of the low expression proteins. Additionally, the correlation coefficients between the codon usage vector of the endogenous genes in the E. coli and the genes for the high or low expressed heterologous proteins (Fig. S3B and S3C) were 0.58 and 0.45, respectively. These results revealed that highly expressed proteins were more likely to use optimal synonymous codons than low-expression proteins, which could facilitate faster translation of those proteins. In light of these findings, we hypothesized that the selection of charged residues [32] and optimal synonymous (i.e., penaltyfree) codons in highly expressed proteins likely enhanced their solubility and translation rate [33,34], respectively.

Evaluation of deep learning model performance
For training the classifiers to accurately predict the effects of a given mutation on expression, we used the dataset constructed in the above section [17]. Since each codon in a given gene sequence could also be categorized according to synonymous codon number, i.e., number of codons encoding the same amino acid (1 to 6), the specific amino acid it called (1 to 20), or its specific nucleotide combination (1 to 61) (Table S1), three classifiers were constructed based on these three coding patterns, respectively. Receiver Operating Characteristics curve (ROC) and Precision-Recall curve (PRC) analyses were used to evaluate the predictive accuracy of each classifier through 10-fold cross validation. Notably, AUROC and AUPRC scores of the classifier that used specific nucleotide combinations were both markedly higher than those of the other two classification strategies, which were 0.764 and 0.751, respectively (Fig. 3B). Additionally, we also considered several alternative metrics, including the prediction, accuracy, F1 score, precision score, and recall score; all of these indexes supported classification based on the 61 combinations of nucleotides as the highest performing model for predicting protein expression levels (Fig. 3C). This result demonstrated the importance effect of codon usage on protein levels.
When the three prediction models were applied to groups 4 and 5, the results showed that group 5 proteins had a significantly greater propensity for high expression in E. coli than that of proteins in the 4th group (p < 0.001) (Fig. 3D). The ratios of predicted average highly expressed propensity of the genes between groups 5 and 4 were 1.01, 1.09, and 1.13 for the synonymous codon number-, amino acid-, and nucleotide combination-based models, respectively. Thus, these collective results all indicated that specific nucleotide combinations were most informative aspect of codon Fig. 2. Comparison of the amino acid or codon frequencies in lowly and highly expressed proteins. A. Amino acids usage differences between the lowly and highly expressed proteins. B-C. Codon frequencies in the highly and-lowly expressed proteins, and in the endogenous E. coli genome. D-E. Codon frequencies in the highly and-lowly expressed proteins, and in the exogenous E. coli genome. In addition, ''ns" denotes no significant, ''*" denotes 0.01 < p-value 0.05, ''**" denotes 0.001 < p-value 0.01, and ''***" denotes p-value 0.001. frequency for predictive modeling of soluble protein expression. We therefore selected the model trained with specific nucleotide combination frequencies to design subsequent mutations.

Correlation between the predicted high-level expression propensity and protein abundance
To evaluate the relationship between the predicted high-level expression propensity and the protein abundance, we collected the protein abundance (PA) data of seven species of unicellular organisms from Paxdb and predicted the high-level expression propensity of those endogenous genes with the three coding schemes. As shown in Table 2 and Fig. S4, most of the spearman rank correlations of the predicted high-level propensity with the coding scheme of 61 combinations of nucleotides were higher than the other two coding schemes. Additionally, among all spearman rank correlations of species, the correlation between the PA of E. coli and the predicted high-level propensity with the coding scheme of 61 combinations of nucleotides exhibited the highest score which was 0.4581. The results indicated that the constructed models based on the analysis of heterologous gene expression propensity in the E. coli expression system gave significant results when predicted endogenous gene expression in various organisms across the tree of life. This suggest that there are universal gene expression codes that our model was able to detect.

Mutation prediction to increase expression of laccase 13B22 and AtGDH
To further validate the MPEPE model, we selected two proteins, laccase 13B22 (GenBank number: MZ817083) and the glucose dehydrogenase, AtGDH (GenBank number: XM_001216916), both of which were reported to have high value in commercial applications and were not included in the training, validation and test dataset [35,36]. Laccase 13B22 can be used in oxidative bioremediation of toxic xenobiotic compounds, while glucose dehydrogenase is an oxidoreductase that catalyzes the oxidation of glucose into gluconic acid d-lactone. Although, the coding sequences of laccase 13B22 and glucose dehydrogenase FAD-AtGDH were both optimized for expression in E. coli, the two proteins still exhibited poor efficiency in their expression in E. coli, which we speculated was related to their codon composition.
In order to ensure that both enzymes remained functional after modification for high, soluble expression, we first examined their conserved residues through protein sequence alignment and residue entropy analysis. Using BLASTP, we identified and downloaded sequences of 1335 and 2413 homologous proteins of 13B22 and FAD-AtGDH, respectively. Based on those collected sequences from NCBI, the position-specific amino-acid probability (PSAP) matrices for both enzymes were calculated. The PSAP cutoff value was set to 0.05, which was the average PSAP value of each amino acid at a given target residue (Figs. S5 and S6). The entropy of each residue was then determined and the higher the entropy, the lower the conservation at that site. For our purposes, we set the demarcation point to 0.65, and higher entropy values indicated that residues did not have conserved function at that position. Finally, the laccase 13B22 and FAD-AtGDH amino acid sequences were individually used as inputs for the MPEPE model to predict all single-point mutations that could improve their soluble expression in E. coli (see more details in the Methods section). As a result, there were 21 and 30 predicted mutations that were satisfied the above three conditions, located at the non-conserved sites and the predicted Fig. 3. Evaluation predictive performance of the three constructed models. A-B. Receiver operator characteristic and precision recall curves for the three models output based on the results of 10-fold cross-validation. C. Model evaluation metrics. D. Prediction results of the three models on the independent test set class4 and class5. In addition, ''ns" denotes no significant, ''*" denotes 0.01 < p-value 0.05, and ''**" denotes p-value 0.01. high-level expression propensity by the constructed model higher than that of the wild-type (Tables S5 and S6). Notably, most of these mutations were located at the C-terminal end and exposed at the protein surface ( Fig. 4A and C). The AlphaFold 2.0 tool was used to simulate the three-dimensional structures of laccase 13B22 and FAD-AtGDH. The location of the mutation sites on the 3D structure were analyzed and their percent solvent accessibilities (PSA) were calculated using Discovery Studio (2016), respectively. The average PSA of the mutations of 13B22 and FAD-AtGDH were 43.67% and 32.12%, respectively ( Fig. 4B and D, Table S6).

Expression levels and enzyme activity validation of the predicted point mutations
To validate the effects of the predicted single point mutations on soluble protein expression in E. coli, we individually constructed the 21 and 30 single-point mutations into laccase 13B22 and FAD-AtGDH, respectively. The proteins of all of the variants was expressed in E. coli under the same induction and expression conditions. Since SDS-PAGE analysis could not clearly distinguish differences in expression between the 13B22 single point mutants and wild type (Fig. S7), we screened for laccase positive mutants using laccase activity assays with culture supernatant to characterize their soluble expression levels. By contrast, SDS-PAGE clearly indicated that FAD-AtGDH single point mutants were expressed at significantly higher levels than that of wild type (Fig. S8). The enzymatic activities of FAD-AtGDH variants and wild type in crude enzyme solution were also examined. The results showed that, among the 21 laccase 13B22 single point mutants, 16 mutants exhibited higher enzymatic activity than that of wild type. In particular, L508F-TTC showed the highest activity of 841.73 ± 33.25 U/ L, 2.854 times higher than that of wild type (Fig. 5A, Table S7). Among the 30 FAD-AtGDH single point mutants, 9 mutants had higher enzymatic activity than the wild type, and of the 15 variants that retained activity after mutagenesis. A454V-GTA had the highest activity (555.49 ± 4.51 U/L), which was 3.05 times greater than that of the wild type (137.02 ± 4.92 U/L) (Fig. 5B, Table S8).
Sequence analysis of these positive mutants revealed that glutamic acid (E) and lysine (K), encoded by two synonymous codons, appeared at higher frequency in the 13B22 and FAD-AtGDH variants, while arginine (R) or leucine (L), encoded by six synonymous codons, were less abundant in the variants than in the wild type ( Fig. 5C and D). These findings suggested that most of variants with improved soluble expression likely benefitted from the charged amino acids and acceleration of translational rate in the protein C-terminal.

Characterization of variants with multiple point mutations
In order to determine whether combining the point mutations could further increase the soluble expression levels of the two proteins, we next generated variants harboring 3, 5, or 7 point mutations (selected based on the magnitude of their effects on soluble expression), and expressed them in E. coli. Soluble expression of 13B22 in culture supernatant gradually increased with increasing numbers of mutated sites ( Fig. 6A and B). In particular, laccase mutant L7M (harboring L508F/L507M/H436E/R431D/R355K/T294 E/N196D) showed the highest expression, approximately 3.49 ± 0.90 times that of wild type (Fig. 6A). In addition, laccase 13B22-L7M activity reached 3954.03 ± 74.03 U/L, approximately 12.40 times that of the wild type (294.938 ± 7.962 U/L) (Fig. 6C). The soluble protein expression of FAD-AtGDH by E. coli BL21 (DE3) also increased significantly with increasing numbers of mutation sites. Notably, AtGDH-A7M (carrying mutations A454V/ Y457M/E493D/S495D/N490D/K373E/E469E) was apparently expressed almost entirely in the soluble fraction ( Fig. 6D and E). The enzymatic activity of the AtGDH-A7M variant was 1213.48 ± 2 07.07 U/L, or 7.86 times greater than wild type (137.02 ± 4.92 U/L) (Fig. 6F). In summary, these results demonstrated that introducing mutations predicted by the preferential usage of charged amino acids in non-conserved sites and codons with fewer synonymous codons significantly improved the expression of these two proteins in E. coli, thereby validating the predictive accuracy of MPEPE.

Discussion
In this study, we constructed DNN-based models for predicting mutations conducive to higher soluble protein expression in E. coli. In particular, we developed a predictive strategy (MPEPE) to screen for mutations that confer these properties and trained the model with gene sequence data calculated from publicly available, experimental, protein expression data. The dataset analysis showed the preferential use of four charged amino acids (E, K, D, and H) in the high-level expressed dataset. These results support previous studies that the charged amino acid can affect translation speed initiation and elongation [37][38][39].
In addition, we evaluated the high-level expression propensity of a human acetylcholinesterase (hAChE) variant bearing 51 mutations, which was approximately 2000-fold higher than that of the wild type in E. coli [40]. Based on our predictive scoring with MPEPE, the 51-point mutant of hAChE was 0.482, higher than the score of 0.463 for the wild type, which suggested an increased propensity for soluble expression. Moreover, the amino acids were mutated to the charged amino acids and the codons with high translational rate (Fig. S9).  Our results support the fact that expression levels are affected by the gene translation rate. In 2004, Wernisch and co-workers introduced the concept of tAI, which measures that adaptation of a coding region to the tRNA pool in the cell [41]; this measure was later generalized to fit various organism [42]. It was suggested in various additional studies that codons with higher adaptation to the tRNA pool (i.e. that are recognized by tRNAs with higher concentration in the cell) tend to have fast translational speed [43,44]. For example, E. coli has four copies of the tRNA gene for codon GAA, encoding glutamate (E), and six copies of the tRNA gene for codon AAA, encoding lysine (K). By contrast, the E. coli genome contained one or zero copies of the vast majority of tRNA genes matching codons for arginine (R) and leucine (L), which both have several synonymous codons, coinciding with the rule we found through development of MPEPE [45]. In addition, the mutation sites that we found could positively affect 13B22 and FAD-AtGDH expression in this study were mainly concentrated in the C-terminus of the protein. These results were consistent with the findings reported in [41] (see also [46,47]) which showed the Nterminus of proteins are translated at a slower rate than the downstream region, resulting in increased translation efficiency and ribosomal allocation.
The predictors designed here are based only on the coding sequence and where trained based on the analysis of the expression genes with the same promoter. Thus, the fact that there is a correlation of 0.46 between the predicted protein levels by our model and the protein abundance of endogenous genes in E. coli suggest that may suggest that at least 20% of the gene expression variability in E. coli can be explained by coding region features. This is a non-negligible value that should be considered when deigning heterologous genes and when studying genome evolution.
Compared to the relatively scattered data that might apply different expression vectors or different expression conditions used by Yu et al., which was only collected from the PDB database [14,48], we used high-quality data published by Letso et al. strictly derived from E. coli under the same expression vector and conditions [17], to construct the MPEPE deep learning model. Although these experimental data cannot be clustered by traditional methods (Fig. S10), our model can effectively extract the features relevant to our protein engineering research question. In addition, our model was validated through mutagenesis and activity assays in laccase 13B22 and FAD-AtGDH, suggesting that the MPEPE provides highly accurate predictions, with broad potential applicability.

Data availability
The gene sequence of laccase 13B22 was uploaded to GenBank under the accession number MZ817083. The details of PCR primers are available in the Supplementary Data.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


extfieldsofstudy:
Computer Science
Engineering

text:
Analyzing the Effects of Handling Data Imbalance on Learned Features from Medical Images by Looking Into the Models

One challenging property lurking in medical datasets is the imbalanced data distribution, where the frequency of the samples between the different classes is not balanced. Training a model on an imbalanced dataset can introduce unique challenges to the learning problem where a model is biased towards the highly frequent class. Many methods are proposed to tackle the distributional differences and the imbalanced problem. However, the impact of these approaches on the learned features is not well studied. In this paper, we look deeper into the internal units of neural networks to observe how handling data imbalance affects the learned features. We study several popular cost-sensitive approaches for handling data imbalance and analyze the feature maps of the convolutional neural networks from multiple perspectives: analyzing the alignment of salient features with pathologies and analyzing the pathology-related concepts encoded by the networks. Our study reveals differences and insights regarding the trained models that are not reflected by quantitative metrics such as AUROC and AP and show up only by looking at the models through a lens.

Introduction
Medical imaging datasets often appear in imbalanced distribution where the frequency of the samples between the different classes of the training dataset is not similar or balanced. The low amount of training samples for infrequent classes or tailed distribution makes it challenging to learn optimal classification boundaries in the representation and can lead to a biased model. Existing methods to tackle class imbalance problem either modify data distribution or learn denotes equal contribution 1 https://github.com/CAMP-eXplain-AI/imba-explain appropriate costs to re-weight class errors. At the data level, the objective is to balance the data distribution through re-sampling techniques which often are prone to over-fitting [18,3,6]. On the other hand, the cost-sensitive approaches modify the learning algorithm to alleviate the bias towards frequent classes or head of distribution [5,26]. The efficacy of these approaches is demonstrated by conventional metrics such as precision and recall and their derivatives. However, the effect of these approaches on the learned features is not well studied. The study of learned features not only improves our understanding of what happens within the models, is specifically insightful when the conventional evaluation metrics do not reflect the effect of applying such cost-sensitive approaches.
This study explores what happens to the learned features when it is trained by cost-sensitive approaches to handle the data imbalance. To understand the effect on learned features, we analyze the internal units of neural networks. Specifically, we analyze the feature maps (outputs of convolutional layers) in deep layers by class activation maps and network dissection [1]. First, we optimize models (ResNet [7] and DenseNet [9]) with Binary Cross-Entropy (BCE) [23]. Several recent and popular cost-sensitive losses such as Weighted BCE (WBCE) [25], Focal loss [21], and Class-Balanced Focal loss (CB-Focal) [5] on NIH Chest Xray [31] dataset and report classical metrics: area Under ROC curve (AUROC), average precision (AP) and predicted probabilities. Then, we visually analyze the impact on salient learned features using class activation maps [33,28] and provide quantitative evaluations to validate the visual observation. We then proceed to analyze the learned features using network dissection [1,16] which quantitatively identifies the concepts encoded (learned) by the model. Statement of Contribution By placing models trained with different learning strategies under the lens, we observe that while metrics such as AUROC and AP report equivalent performance, the models trained with cost-sensitive losses encode more pathology-related concepts. Moreover, we observe an increased alignment between salient learned features and pathology-related features.

Related Works
Handling Data Imbalance: Much of the recent works on imbalanced learning focused on alleviating this problem using novel objective function. Lin et. al. [21] introduce a Focal loss for dense object detection where class-specific weights were automatically learned. Cui et. al. [5] re-weight the loss by the inverse effective number of examples to learn balanced representations. Similarly, [18] modified the weights according the uncertainty of predictions. Others address this problem by multi-task learning [10,26] that used selective instances for training on imbalanced sets for each task. Interpreting Neural Networks: Two principal neural network interpretation approaches are feature attribution [29,32,13,30,22,17,12] (i.e. saliency methods [4]) and analyzing internal units (e.g. feature visualization [24] and dissection [2,16]). Here, we look into the models from both perspectives. We are interested in both the contribution of input features to the output (i.e., feature attribution) and concepts encoded by the network (via analyzing units).
The are many attribution methods, however the identified important features are different for different methods [20,32,15,14]. This is a caveat for researchers using the attribution toolkit. Within the various feature attribution methods, CAM [28] being a classic and intuitive method, it satisfies benchmarks in terms of faithfulness [8,32,15]. Moreover, regardless of its advantages for attribution, the method provides a summary of activation maps in a certain layer by performing a weighted sum of the attributions. We use network dissection [1,16] to identify concepts associated with individual convolutional feature maps.

Method
The question we explore in this paper is what happens to the learned features within the model when we apply cost-sensitive approaches to handle the class imbalance during training. Specifically, we study the intra-class data imbalance, i.e. the imbalance between the number of positive and negative samples for each class. We first introduce the cost-sensitive methods (see Section 3.1) that consider the data imbalance between positive and negative samples. Then, we explain our approaches for analyzing the learned features (see Section 3.2).

Handling Data Imbalance
, our goal is to learn the optimal boundary θ * obtained by empirical loss minimization L D (θ) on the training set D as: θ * = arg min θ L D (θ). The class imbalanced problem exists when the frequency of the samples among different categories are extremely mismatched. Therefore θ learned on D using conventional empirical loss can be biased towards the low frequent classes and significantly different from the ideal boundary θ * . In other words, because of the imbalanced proportion between classes, the optimal boundary is more likely to afford a higher empirical error than an alternative hypothesis based on an empirical loss. This is due to the nature of imbalanced class distribution that forces the classifier to shift θ closer to low-frequent classes because it reduces empirical error. A principal strategy to handle data imbalance is through the cost-sensitive loss functions [21,5,10,27].
Here, we study data imbalance problem in multi-label and binary classifica- m of class m for sample x (i) given by a neural network is defined with multiple outputs as cross entropy error between predicted value and true label using Binary Cross Entropy (BCE) loss. The cost-sensitive losses applied to class can be summarized into a unified framework: where w (i) − are the positive and negative class weights for sample x (i) , respectively. We omit the class index m for simplicity. The four losses considered in this work can be derived as follows: where N + , N − are the number of positive and negative samples of class m, and α, γ and β are hyper-parameters. As we are discussing intra-class imbalance, each loss is considered individually and we apply the cost-sensitive methods on each loss to strike a balance between positive and negative samples for each class.

Looking into the Models
The objective is to analyze the features learned by the model in different training scenarios. In this work, we analyze the internal units of the model, specifically the activation pattern of convolutional feature maps in deeper layers. The deep convolutional layers are chosen due to the following two reasons: 1) Final layers in principle encode higher-level concepts and we are interested to know if these concepts align with pathological features. 2) We focus on convolutional layers instead of deeper fully connected layers, since convolutional layers keep the spatial correspondence with input features, and therefore, we can leverage image annotations from experts to study if the activations align with pathology related features. We analyze the internal features maps from two different perspectives: Class Activation Maps [33,28] can be deemed as a tool that summarizes the contributions of feature maps of the final convolutional layer in one tensor (for the networks under investigation, GradCAM [28], and CAM [33] are equivalent). The method combines the feature maps using their associated weight (or gradient of output with respect to feature map in GradCAM) of their connections to the output, thus summarizing how feature maps contribute to the final output. The method, in essence, shows which areas are contributing to the output for each sample input. Thus using annotations from experts, we can check whether the salient regions are aligned with the pathologies [31,11].
Analyzing Encoded Concepts [1,16]: We use an approach inspired by network dissection [1] which is a tool for individually analyzing the convolutional feature maps of the neural network. In essence, it identifies the concepts associated with each feature map, thus allowing for understanding what concepts are encoded by the network during training. In this work, we use a conceptually similar approach. For a chosen convolutional layer, let A c denote the activation (feature) map of channel c. Similar to [1], we first apply a threshold τ c with P(A C ≥ τ c ) = q to activation values of each channel c to filter the dataset-wise significant activation values, where q is a pre-defined hyper-parameter. For each image in the annotation dataset we compute the corresponding A c and apply the threshold τ c . Subsequently, we identify the connected components in the threshold activations. If the connected component overlaps the bounding box region, we consider the connected component as a pathology related concept. Using this approach we can report two values that reflect the concepts encoded in the network quantitatively, which are explained in Sec. 4.

Experiments
Experimental Setup In this paper, we examine and perform two different datasets and learning tasks: multi-label classification and binary classification on the NIH Chest X-ray dataset [31]. For binary classification. We define samples with the "No Finding" label as "Healthy" samples, and samples with other labels as "Unhealthy" samples. Random crop, horizontal flip, and color jitter are used as the augmentations. In the end, the augmented images are resized to 224 × 224 before being fed into an image classifier.
We use two ImageNet pre-trained classifiers: ResNet50 [7] and DenseNet [9]. The last feature maps of ResNet50 have the size 7 × 7. In order to have a finegrained analysis on smaller features such as Nodule (for which 7 × 7 is a coarse resolution) we use a truncated version of DenseNet. We discard the last two dense blocks and the associated transition blocks, so that the last feature maps have size 28 × 28. We refer to this truncated variant of DenseNet as T-DenseNet. We adopt the same training configurations for all models. Specifically, we train the classifiers on 4 GPUs (DGX-A100) for 50 epochs using the Adam [19] optimizer with weight decay 10 −6 and initial learning rate 0.0004. The learning rate is decayed following the cosine-annealing policy. In addition, we set the batch size to 512. For the Focal loss, the γ and α set to 2.0 and 0.25, respectively; for the class-balanced focal loss, we set γ to 2.0 and β to 0.9999.

Feature Alignment Analysis
The objective is to have a metric for measuring the alignment between the contributing features (from CAMs) for a single prediction and the pathologies. This can be achieved by classical object detection metrics IoBB(Intersection over Bounding Box) and IoR(Intersection of detected Region). IoBB and IoR can be considered as the visual counterparts of recall and precision. Computing these metrics requires thresholding the CAMs. In order to reduce the sensitivity of the results to the chosen threshold, we implement soft IoBB and IoR computation. This is performed by normalizing the CAMs to [0, 1] and applying a summation of values within the regions under consideration (intersection, bounding box, detected region). for all images and normalize by the number of images. The provided metrics reflect how individual feature maps are aligned with pathologies.

Results and Discussion
Performance Analysis We evaluate the classification performance using AU-ROC and AP (Fig. 2). The performance differences between BCE and weighted losses (WBCE, Focal, and CB-Focal) are marginal. Therefore the effect of handling data imbalance via losses is not visible via these metrics. We also report the predicted probabilities of different losses. The predicted probabilities of the weighted losses are substantially higher than that of the BCE. This is expected as more weight is given to positive samples. Moreover, as more weight is given to the positive weights, the recall is higher. However, it is not clear if the increased recall is due to model's inclination towards identifying samples as positive or due to the model having learned and leveraging more predictive features.

Feature Alignment Analysis
In this section, we analyze the models' internal feature maps to see if the increased recall and probability values are related to predictive features relevant to the pathologies. This is realized by comparing CAMs and the salient features for individual predictions with annotations from radiologists. In Fig. 1 we observe that in models trained with cost-sensitive losses, the additionally activated regions are aligned with the pathologies (bounding boxes). This shows that the increased recall is not due to the models' propensity towards identifying samples as positive but due to having learned meaningful Task  Multi-label  Binary  Classifier T-DenseNet  ResNet50  T-DenseNet  ResNet50  Metric IoBB  IoR  IoBB  IoR  IoBB  IoR  IoBB  IoR  BCE Table 1: Feature Alignment Analysis. The average IoBB and IoR results for test sets images that have bounding box annotations. These metrics are provided to quantitatively validate the observation in Fig. 1   predictive features. Note that metrics such as AUROC and AP are roughly equivalent for all models (Fig. 2). In order to quantify the observed behavior in Fig. 1 we report IoBB and IoR results in Tab. 1. The increased IoBB for models trained with cost-sensitive losses shows that the features cover more areas of bounding boxes (pathologies). We also observe a decrease in IoR, which is also evident in Fig. 1 as falsely positive active regions have also increased.

Analysis of Concepts
This section analyzes the effect of handling data imbalance on the concepts encoded by the model. This is performed by counting the number of detected pathology-related concepts. In Tab. 2 we observe that the models trained with cost-sensitive losses consistently have more detectors. Thus it can be concluded that the models encode more pathology-related concepts when the data imbalance is considered in the loss.

Conclusion
In this work, we study the effect of handling data imbalance using cost-sensitive losses during training on the learned feature maps. The feature maps are analyzed from two perspectives: class activation maps and the concept encoded by each feature map. We observe that although classical metrics such as AUROC and AP report equivalent performance for the studied models, the learned features are different in these models. When data imbalance is handled, the models encode more pathology-related concepts. Moreover, we observe that overall increased recall and predicted probability of these models are accompanied by an increased alignment between learned features and pathology-related features.


extfieldsofstudy:
Computer Science
Medicine

text:
Distance-Based Detection of Cough, Wheeze, and Breath Sounds on Wearable Devices

Smart wearable sensors are essential for continuous health-monitoring applications and detection accuracy of symptoms and energy efficiency of processing algorithms are key challenges for such devices. While several machine-learning-based algorithms for the detection of abnormal breath sounds are reported in literature, they are either too computationally expensive to implement into a wearable device or inaccurate in multi-class detection. In this paper, a kernel-like minimum distance classifier (K-MDC) for acoustic signal processing in wearable devices was proposed. The proposed algorithm was tested with data acquired from open-source databases, participants, and hospitals. It was observed that the proposed K-MDC classifier achieves accurate detection in up to 91.23% of cases, and it reaches various detection accuracies with a fewer number of features compared with other classifiers. The proposed algorithm’s low computational complexity and classification effectiveness translate to great potential for implementation in health-monitoring wearable devices.

Introduction
Acoustic signal processing is a growing area of focus in healthcare and biomedicine. The characteristics of acoustic signals, especially that of abnormal breath sounds, provide clinicians with valuable information on respiratory diseases [1]. Traditionally, clinicians use stethoscopes as an indispensable tool for healthcare delivery for their reliability and efficiency in aiding the investigation of bodily sounds. However, the interpretation of breath sounds heard through a stethoscope is user dependent, and its single examination mode of measurement disallows continuous monitoring of daily symptoms and fluctuations [2]. With the progression of sensing technology development, wearable devices are now available as promising solutions. They facilitate automatic and continuous acoustic analysis and assist clinicians in evaluating the effectiveness of a prescribed intervention [3].
Amongst the various breath sounds, wheeze is of particular interest as its presence and duration are a significant reference for clinicians to diagnose and monitor various pulmonary pathologies, including chronic obstructive pulmonary disease (COPD), bronchiolitis, and asthma [4]. In addition to wheezing, cough is another common respiratory sound and significant symptom of interest. Its pattern and changes over time can indicate disease evolution in more than 100 different diseases, and its resolution reflects the effectiveness of therapy [5].
The biological features of cough, wheeze, and breath sounds can be characterized by both the spectral and temporal analysis [6,7]. For example, the wheezing sounds during respiration occur when the airway is obstructed or narrowed; hence, the fundamental frequency of the respiratory sound was of a higher frequency when compared with a normal breath. Furthermore, wheeze is often perceived as musical tone where the harmonics can be observed in the spectral domain. Coughs, although having similar properties as wheezes, such as having a higher fundamental frequency and harmonics, are explosive in nature, whereas wheezes are continuous and more commonly found in the expiration phase rather than the entire respiratory cycle. Given these biological characteristics of wheeze and cough sounds, various classification algorithms were developed to differentiate wheeze and cough signals using their spectral and temporal features. Applying hidden Markov models on the extracted cepstral features from a database of 2155 cough events, Matos et al. [8] achieved 82% detection accuracy in the binary classification of cough and noncough sound samples. Aydore et al. [9] investigated 246 wheeze and non-wheeze epochs using Fisher's ratio and Neyman Pearson hypothesis and achieved 93.5% sensitivity. Jain et al. [10] extracted spectral features from 19 wheeze samples and 21 other lung sounds and applied a multi-stage fixed threshold classifier for 85% overall accuracy. Lin et al. [11] developed a bilateral filter-based edge-processing technique for 90 wheeze samples and detected 87 out of them. Time-series deep learning models were also examined in previous studies. For example, Justice Amoh [12] used the convolutional neural network and recurrent neural network for cough detection among 14 subjects and achieved 82.5% and 79.9% overall accuracy, respectively. Arati Gurung [13] reviewed 12 recent studies of automatic abnormal lung sound detection, and the average accuracy was 85%. However, all the above studies focused on the detection of a single symptom; the robustness of such methods is not addressed if applied to multiple concurrent symptoms, a feature of 'real-world' clinical practice. To the best of the authors' knowledge, the only study that explored the overall detection of multiple concurrent symptoms including cough and wheeze was performed by Himanshu S. Markandeya [14]. This system employed a two-stage classification structure and achieved 90% accuracy from 74 recordings. It identified discrete wavelet transform (DWT) coefficients for symptoms and applied a specific processing method for each identified DWT coefficient to determine symptoms. This classification algorithm, however, assumes that the DWT coefficient of wheeze sounds are resolved as D 5 , but another study [15] illustrated that wheeze data exhibits useful information throughout D 6 to D 2 . Therefore, using exclusive processing methods for specific DWT coefficients would be unsuitable when the actual symptom is not from its expected category.
Health monitoring wearable devices are designed to simultaneously detect multiple respiratory sounds (and perhaps symptoms) to provide comprehensive health and medical insights to the wearer and healthcare providers. However, previous studies mainly focused on a binary classification problem in cough or wheeze detection. Extending them to multiclass classification requires different model signs and higher model complexity to attain satisfactory detection performance. Moreover, none of these existing models consider the unique patterns (for example, distribution of each class in each feature) of various respiratory sounds. Researchers have long been calling for simple yet effective methods to diagnose and monitor patients on wearable devices [16].
This study addresses the challenges of existing cough or wheeze detection studies, attempting to find an optimal acoustic signal processing method for implementation in a wearable device. The proposed method includes a kernel-like input mapping strategy by first transforming the original sound signals in time series to a higher dimensional feature space that characterizes the temporal and spectral patterns of respiratory sounds, followed by a tailored dimension reduction for the algorithm to achieve robustness and generalizability by distilling the prediction information from only a few features. Since data processing is one of the main consumers of battery power in wearable devices [17], the superior prediction performance and minimized feature extraction in the proposed method would make such devices more suitable in clinical applications. This paper is organized as follows. Section 2 describes the collected dataset, the proposed data-processing methods, and its unique features. It also includes a short discussion on its implementation on an embedded system. In Section 3, the results from the proposed data processing methods are presented and analyzed. Finally, the findings are discussed in Section 4.

Data
For the generalizability of the study, raw data were collected from a variety of sources, as illustrated in Table 1. These sound signals were recorded using different devices with a large sampling frequency and noise-signal ratio variation. Hospital patients' data were collected either in the outpatient clinic or in-patient ward of a hospital in Singapore, with informed consent obtained from all participants. Respiratory sounds were recorded over the right side of the chest wall using an electret microphone with a sampling rate of 8 kHz and a dedicated cavity for air-coupling purposes. Other data types were included in the final analysis. The open-source datasets were selected from readily available sound libraries or published works [18]. Participant subjects were volunteers who contributed breathing or coughing sounds by recording with their own devices. In this project, 604 segments of acoustic data were acquired and converted to WAV digital format. Each segment might contain one or several respiratory cycles. The total duration of the segmented data was 1935s.

Pre-Processing
Each segment was manually trimmed from its original recording to contain one complete bout of breathing, wheezing, or coughing event. Due to the variability of recording devices, the number of channels in each segment was different. These segments were standardized to mono signals by averaging all channels. Following this, a band-pass filter was implemented, as different respiratory sounds only exhibit frequency peaks in certain regions, and signals outside the range are likely noise. For wheeze, a frequency range of spectral peaks was reported as 80-1600 Hz [19]; whereas for cough sounds, the frequency components were found to occur up to 3.0 kHz, with the majority of peaks appearing at 100-900 Hz for various cough types [15]. In this study, the band-pass filter was designed to only keep 80-3000 Hz signals.

Proposed Distance Function
In the choice of machine learning algorithms, it is essential to achieve both accuracy and computational simplicity to be applicable for wearable devices. Considering a sound sample of length t, x = ( x 1 , x 2 , . . . x t ), the detection task can be generalized as: where k ∈ (1, 2, 3) are three detection classes (cough, wheeze, and breath), c k are the representative signals of each class, and l k are the respective likelihood functions. In this study, the negative distance between the sample and cluster center was used to assess Note that D can be any distance metric. However, as x is a time series with unequal length, directly measuring the distance is time consuming and computationally expensive [20,21]. Instead, a kernel-like trick was proposed to transform the segmented frames from the time series into well-established high-dimensional acoustic features, x i = x i,1 , . . . , x i,j , where i is the ith segment of x. Various temporal and spectral features were extracted as cough and wheeze sounds have significant features in both time and frequency domains [22,23]. In the spectral domain, extracted features included melfrequency cepstral coefficients (MFCC), spectral roll-off, spectral entropy, etc. Temporal features included amplitude and change to the energy level. In total, j = 72 features were extracted in the sequence, as shown in Table 2. This feature set constructed the feature mapping from the time series to the static high-dimensional feature space. Equation (2) is a simple mathematical model, with essential importance on the choice of distance metric D, as it needs to capture the distinct characteristics of three classes of acoustic signals in each feature while satisfying the restriction of low computational power available on wearable devices. To capture the unique characteristics of different signals in each feature, the distance function D should consider not only the absolute distance between the feature value x i and cluster center c k = c k,1 , . . . , c k,j , but also the feature importance, variance, and distribution (see details in Supplementary Part A). For any feature j, such a distance metric between a sample x and an output class center c k,j is defined as: where σ k,j is the standard deviation of jth feature of class k, and ρ x j , c k,j is the probability density matrix that the feature value x i,j belongs to class k evaluated using the histograms in training data. Incorporating σ k,j into the denominator avoids the distance scaling issue where different classes and features have different variances. ρ x i,j , c k,j is similar to the concept of the value difference matrix in memory-based reasoning [24], defined as the probability that a sample belongs to class k when only considering feature j. As this probability density only considers a single feature at a time, this assumption seems to be similar to a naïve Bayesian classifier. However, the introduction of feature weight w j mitigates the impact of such an assumption: when the feature is highly correlated with other features, w j decreases to minimize the correlation bias. Hence, w j can be evaluated using the game theory approach by comparing the detection accuracy before and after the addition of feature j: where p −j and p are the detection accuracies with/without feature j.

Proposed Feature Selection Algorithms
As there were 72 features extracted, using all features on wearable devices is power consuming and unnecessary. Furthermore, the classifier is more likely to overfit the training data. For robustness and generalizability, we want the classifier to be as simple as possible, using the fewest number of features (Occam's razor). There are various feature reduction methods available to evaluate the significance of the extracted features. In this study, we proposed a tailored feature reduction/selection method using the proposed distance function.
A better feature selection method should reflect the real discriminant power of features with respect to all output classes; meanwhile, it should consider the interactions between combined features. To achieve this, a recursive feature selection algorithm was proposed by integrating the proposed distance function in Equation (2). To reflect the real discriminant power of features with respect to all output classes and consider the interactions between combined features, the proposed method selects the feature subsets that minimize the distance from the training sample to the center of its class while maximizing distance from other classes. Let SS d be the collection of all feature subsets with cardinality d, containing any feature subset S d ∈ SS d , and x i be the feature vector of the ith training data segment . The proposed algorithm selects the best feature subset S d from SS d in four steps. First, it takes one feature subset S d and calculates the distances D k (x i , c k , S d ) from each feature vector x i to output classes c 1 , c 2 , and c 3 according to Equation (3). Secondly, it labels each feature vector x i with the class that has the minimum distance, denoted by argmin Thirdly, it compares the labeled classes to the ground truth for all samples and determines the accuracy of the features subset S d . Finally, it repeats the first three steps for other candidate feature subsets S d in SS d and detects the feature subset with the highest accuracy. This process can be expressed as: where c x i is the real output class of feature vector x i and 1 D( is an indicator function that returns 1 if two distances are equal and returns 0 otherwise. When selecting the best feature subset S d containing d features, the above calculation must be repeated for 72 d times to exhaust all feature subsets from SS d . This could be time-consuming. To shorten this process, a recursive algorithm was proposed to start with a small number of features and expand the subsets by gradually adding more features into feature subsets. The stepwise algorithm (Algorithm 1) is as follows: In this study, l is set to be 500. 6: Step 3: For every selected subset S m from Step 2, find the absolute complement of S m , 7: represented by S m . For each feature f i ∈ S m , create a new subset S m+1 by adding f i 8: to S m ; and keep this new subset S m+1 only if its performance is better than S m . 9: Step 4: Insert all remaining S m+1 from step 3 to new feature set SS m+1 . Increase m by 1. 10: Step 5: Repeat step 2 to 4 until m = d. With the designed distance metric and the feature selection algorithm, the kernel-like minimum distance classifier (K-MDC) can be written as: In essence, K-MDC transforms the original one-dimensional acoustic time series into a high-dimensional static feature space (like the kernel trick in support vector machines). Finally, the feature space is further reduced from high dimensional to d-dimensional to optimize the classification accuracy and robustness.

Proposed Embedded Architecture for K-MDC Implementation
After identifying a d-optimized feature space from the original one-dimensional acoustic training data samples, K-MDC can be implemented on an embedded solution using only the necessary d features extraction. Concretely, we proposed the following embedded architecture where we split the processing into two steps: (a) feature extraction, as shown in  In essence, K-MDC transforms the original one-dimensional acoustic time series into a high-dimensional static feature space (like the kernel trick in support vector machines). Finally, the feature space is further reduced from high dimensional to -dimensional to optimize the classification accuracy and robustness.

Proposed Embedded Architecture for K-MDC Implementation
After identifying a -optimized feature space from the original one-dimensional acoustic training data samples, K-MDC can be implemented on an embedded solution using only the necessary features extraction. Concretely, we proposed the following embedded architecture where we split the processing into two steps: (a) feature extraction, as shown in

Results and Discussion
To evaluate the proposed K-MDC, experiments were designed and conducted the study as follows. First, detection accuracy was assessed in various output classes by increasing the model complexity. Next, K-MDC was compared with the state-of-art algorithms in the mapped high-dimensional feature space. Last, the impact of the feature reduction algorithm was explored both quantitatively and qualitatively. In this section, the quantitative results were reported based on five-fold cross validation. Additionally, a demonstration of implementation in an embedded wearable based on the Renesas synergy internet-of-things (IoT) platform was evaluated.
In essence, K-MDC transforms the original one-dimensional acoustic time series into a high-dimensional static feature space (like the kernel trick in support vector machines). Finally, the feature space is further reduced from high dimensional to -dimensional to optimize the classification accuracy and robustness.

Proposed Embedded Architecture for K-MDC Implementation
After identifying a -optimized feature space from the original one-dimensional acoustic training data samples, K-MDC can be implemented on an embedded solution using only the necessary features extraction. Concretely, we proposed the following embedded architecture where we split the processing into two steps: (a) feature extraction, as shown in Figure 1

Results and Discussion
To evaluate the proposed K-MDC, experiments were designed and conducted the study as follows. First, detection accuracy was assessed in various output classes by increasing the model complexity. Next, K-MDC was compared with the state-of-art algorithms in the mapped high-dimensional feature space. Last, the impact of the feature reduction algorithm was explored both quantitatively and qualitatively. In this section, the quantitative results were reported based on five-fold cross validation. Additionally, a demonstration of implementation in an embedded wearable based on the Renesas synergy internet-of-things (IoT) platform was evaluated.

Results and Discussion
To evaluate the proposed K-MDC, experiments were designed and conducted the study as follows. First, detection accuracy was assessed in various output classes by increasing the model complexity. Next, K-MDC was compared with the state-of-art algorithms in the mapped high-dimensional feature space. Last, the impact of the feature reduction algorithm was explored both quantitatively and qualitatively. In this section, the quantitative results were reported based on five-fold cross validation. Additionally, a demonstration of implementation in an embedded wearable based on the Renesas synergy internet-of-things (IoT) platform was evaluated.

Performance of the K-MDC in Various Output Classes
To evaluate the performance of K-MDC in three types of acoustic signals, prediction accuracies were recorded while increasing feature subset size of S d from d = 1 to 10. Feature subsets were gradually expanded by the proposed recursive feature selection algorithm in Section 2.4. The detection accuracies for each output category as well as overall performance were plotted in Figure 3. K-MDC was able to achieve overall detection accuracy of more than 90% and it delivered the best detection performance for cough data (≥ 95% when d ≥ 4). The detection accuracies of all classes increased more drastically when d ≤ 4, and tended to taper as the number of features further increased. When d = 10, the overall detection accuracy was 91.23%.

Performance of the K-MDC in Various Output Classes
To evaluate the performance of K-MDC in three types of acoustic signals, predictio accuracies were recorded while increasing feature subset size of from = 1 to 10. Fe ture subsets were gradually expanded by the proposed recursive feature selection alg rithm in Section 2.4. The detection accuracies for each output category as well as overa performance were plotted in Figure 3. K-MDC was able to achieve overall detection accu racy of more than 90% and it delivered the best detection performance for cough data ( 95% when ≥ 4). The detection accuracies of all classes increased more drastically whe ≤ 4, and tended to taper as the number of features further increased. When = 10, th overall detection accuracy was 91.23%.

Comparison with State-Of-Art Algorithms
Next, the proposed K-MDC was compared with state-of-art machine learning alg rithms which include naïve Bayes (NB), support vector machine (SVM), and artificial neu ral network (ANN), which are widely used classifiers in acoustic signal processing. T study the effects of deep learning in this detection problem, two variants of ANN mode were implemented: a shallow neural network (SNN) with 1 hidden layer of size 10 and deep neural network (DNN) with three hidden layers of sizes (128, 64, and 32). MATLA Statistics and Machine Learning Toolbox and MATLAB Neural Network Pattern Reco nition Toolbox were used for implementing DN, SVM, and SNN, and Python Pytorch [2 was used for implementing DNN. The comparison of K-MDC with baselines was re orded at various feature dimensionalities in Figure 4. It shows that the proposed K-MD method achieved superior performance compared with the implemented machine-learn ing classifiers in most cases, implying its capability of distilling predictive informatio from limited feature dimensionality. When comparing SNN with DNN models, bot models achieved similar prediction accuracy with a larger number of input feature When given a smaller number of input features, DNN was more capable of distilling in formation and produced better prediction than SNN when ≤ 3. However, such compl cated architecture was less stable and was outperformed by SNN when 4 ≤ ≤ 9.

Comparison with State-of-Art Algorithms
Next, the proposed K-MDC was compared with state-of-art machine learning algorithms which include naïve Bayes (NB), support vector machine (SVM), and artificial neural network (ANN), which are widely used classifiers in acoustic signal processing. To study the effects of deep learning in this detection problem, two variants of ANN models were implemented: a shallow neural network (SNN) with 1 hidden layer of size 10 and a deep neural network (DNN) with three hidden layers of sizes (128, 64, and 32). MATLAB Statistics and Machine Learning Toolbox and MATLAB Neural Network Pattern Recognition Toolbox were used for implementing DN, SVM, and SNN, and Python Pytorch [25] was used for implementing DNN. The comparison of K-MDC with baselines was recorded at various feature dimensionalities in Figure 4. It shows that the proposed K-MDC method achieved superior performance compared with the implemented machine-learning classifiers in most cases, implying its capability of distilling predictive information from limited feature dimensionality. When comparing SNN with DNN models, both models achieved similar prediction accuracy with a larger number of input features. When given a smaller number of input features, DNN was more capable of distilling information and produced better prediction than SNN when d ≤ 3. However, such complicated architecture was less stable and was outperformed by SNN when 4 ≤ d ≤ 9.
As the requirement of smaller features implies a shorter computation time and power for feature extraction and classification (hence faster detection speed and working hours) on wearable devices, it is interesting to investigate the number of features required to achieve certain levels of prediction accuracy for each machine-learning classifier. Moreover, simpler models also lead to more robustness and better generalizability, according to Occam's razor. The minimum required feature dimensionality to achieve prediction accuracy of 80%, 85%, and 90% for each classifier, respectively, was measured. In addition to the abovementioned classifiers, this test was extended to K-nearest neighbors (KNN) using MATLAB Statistics and Machine Learning Toolbox. The optimal number of K in KNN was found to be K = 18, and the results of all classifiers were recorded in Table 3. As can be seen, K-MDC reached various accuracy thresholds with the smallest number of required features. This implies the best distillation and utilization of feature information in K-MDC as it takes variance, density distribution, and mutual information of each class and each feature into consideration. The sensitivity analysis of each feature in the K-MDC model was conducted in Supplementary Part B using the game theory approach.  As the requirement of smaller features implies a shorter computation time and pow for feature extraction and classification (hence faster detection speed and working hour on wearable devices, it is interesting to investigate the number of features required achieve certain levels of prediction accuracy for each machine-learning classifier. More ver, simpler models also lead to more robustness and better generalizability, according Occam's razor. The minimum required feature dimensionality to achieve prediction acc racy of 80%, 85%, and 90% for each classifier, respectively, was measured. In addition the abovementioned classifiers, this test was extended to K-nearest neighbors (KNN) u ing MATLAB Statistics and Machine Learning Toolbox. The optimal number of K in KN was found to be K = 18, and the results of all classifiers were recorded in Table 3. As ca be seen, K-MDC reached various accuracy thresholds with the smallest number of r quired features. This implies the best distillation and utilization of feature information K-MDC as it takes variance, density distribution, and mutual information of each cla and each feature into consideration. The sensitivity analysis of each feature in the K-MD model was conducted in Supplementary Part B using the game theory approach.   Table 3. Minimum required number of features to achieve specific detection accuracies. '-' indicates that the classifier could not achieve the required accuracy threshold using our dataset.

Performance of Feature Selection
Performance of the proposed feature selection algorithm was both quantitatively and qualitatively evaluated. Table 4 lists the average number of features required for each machine-learning algorithm to converge to its optimal accuracy, where the convergence was identified when the accuracy was less than 3% away from the highest accuracy. Note that the highest testing accuracy may not be attained with the highest dimensionality because the higher dimensionality leads to more model complexity and overfitting. To exemplify this, we calculated the accuracy difference between the model constructed with the proposed feature selection method and the model constructed using all features. As shown in Table 4, these classifiers converge to their optimal performance with fewer than 8 features out of 72 features with the proposed feature selection method. The best performance improvement (in the case of K-MDC) was up to 13.24% when compared with the performance using randomly selected features with the same dimensionality. However, the ANN models (both SNN and DNN) with feature selection did not perform better than without feature selection, even though feature selection reduced total computation effort. This is due to the nature of the ANN architecture with early stopping, and such a finding is in accordance with other studies using different datasets [26]. Therefore, great caution must be taken when performing dimension reduction techniques with neural networks to consider the trade-off between computation time and detection accuracy. Computation complexity analysis of the proposed model and other baselines was discussed in Supplementary Part C. Table 4. Performance difference with or without feature selection methods and required number of features for the detection accuracy to converge, where the convergence is defined as less than 3% from the optimal accuracy.

Classifiers
Performance Compared with other feature selection methods, for example, principal component analysis (PCA), the proposed method also performed better. PCA selected spectral roll-off (f 40 ) and energy variation (f 54 ) as the best two-feature combinations, contributing to a total of 99.57% of its explained variance. When using these two features for classification, it yielded an overall detection accuracy of 72.28%, which was inferior to the selected two-feature combination (f 2 and f 53 ) in the proposed algorithm with an accuracy of 85.20%.
To qualitatively visualize the superiority of the proposed method, Figure 5 illustrates the clustering of samples by recursively expanding the feature space from one dimensional to three dimensional. Training samples belonging to different classes were initially hard to distinguish when only one feature was chosen from the feature subset. In subsequent steps, the proposed algorithm selected one feature each time to optimize the objective function and the separation between different classes became increasingly distinct.
out of 72 features with the proposed feature selection method. The best performance improvement (in the case of K-MDC) was up to 13.24% when compared with the performance using randomly selected features with the same dimensionality. However, the ANN models (both SNN and DNN) with feature selection did not perform better than without feature selection, even though feature selection reduced total computation effort. This is due to the nature of the ANN architecture with early stopping, and such a finding is in accordance with other studies using different datasets [26]. Therefore, great caution must be taken when performing dimension reduction techniques with neural networks to consider the trade-off between computation time and detection accuracy. Computation complexity analysis of the proposed model and other baselines was discussed in Supplementary Part C. Compared with other feature selection methods, for example, principal component analysis (PCA), the proposed method also performed better. PCA selected spectral rolloff ( ) and energy variation ( ) as the best two-feature combinations, contributing to a total of 99.57% of its explained variance. When using these two features for classification, it yielded an overall detection accuracy of 72.28%, which was inferior to the selected twofeature combination ( and ) in the proposed algorithm with an accuracy of 85.20%. To qualitatively visualize the superiority of the proposed method, Figure 5 illustrates the clustering of samples by recursively expanding the feature space from one dimensional to three dimensional. Training samples belonging to different classes were initially hard to distinguish when only one feature was chosen from the feature subset. In subsequent steps, the proposed algorithm selected one feature each time to optimize the objective function and the separation between different classes became increasingly distinct. (a)

Implementation of the Proposed Algorithm in Embedded Systems
Finally, we implemented the K-MDC using the proposed architecture described in Section 2.5 onto an embedded system based on Koh et al. [27] running on the Renesas S5D9 120 MHz Arm ® Cortex ® -M4 CPU in an IoT environment. The classification results from K-MDC are published on an IOT cloud, and the results are extracted through a smartphone connected to the same IOT cloud. Specifically, in Figure 6, we show the results for the classification of (a) cough, (b) breath, and (c) wheeze. In all cases, the classified sounds are preloaded in the embedded system as the main purpose is to determine the feasibility of K-MDC in a wearable device. The feature space is constructed by recursively adding features into feature subsets according to Algorithm 1. Each dot represents an acoustic segment of cough (red), breath (blue), or wheeze (green) sound.

Implementation of the Proposed Algorithm in Embedded Systems
Finally, we implemented the K-MDC using the proposed architecture described in Section 2.5 onto an embedded system based on Koh et al. [27] running on the Renesas S5D9 120 MHz Arm ® Cortex ® -M4 CPU in an IoT environment. The classification results from K-MDC are published on an IOT cloud, and the results are extracted through a smartphone connected to the same IOT cloud. Specifically, in Figure 6, we show the results for the classification of (a) cough, (b) breath, and (c) wheeze. In all cases, the classified sounds are preloaded in the embedded system as the main purpose is to determine the feasibility of K-MDC in a wearable device.

Conclusions
This study proposed a kernel-like detection algorithm for classifying acoustic time series into cough, wheezing, and breathing. It has a novel distance measure that considers each feature's unique properties, including the class variance, probability distribution, and feature importance. Results showed that the proposed method achieved better detection accuracy than existing algorithms, and the proposed feature reduction method effectively reduced the feature dimensions to 4 from 72 and, at the same time, improved the overall classification by 13.24%. In wearable device applications, feature dimensionality is of concern as a smaller number of features requires shorter computation time for feature extraction and classification (faster detection speed) and hence less consumption of battery power (longer working hours). From our findings, the proposed method potentially tackles the current challenges of limited processing power by minimizing the feature extraction efforts and detection accuracy in health-monitoring wearable devices by capturing the unique distribution characteristics in different features and sound classes. This is also exemplified in this work through the implementation of the proposed algorithm in an embedded platform based on the IOT chipset from the Renesas synergy series. Furthermore, the findings in this study may be applied to other health-monitoring applications, including the quantification and severity analysis of clinical symptoms.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Figure S1: Variances of features by class. For easy comparison, the variances of breath are scaled to 1 for all features; Figure S2. The distribution patterns of features by classes. Bar shades are the relative frequencies for each class and curves are the fitted probability density distributions. Figure S3. Bimodal and multi-modal distributions of features in three classes. Bars are the relative frequencies for each class and curves are the fitted probability density distributions. Figure S4. Sensitivity analysis of features in each outcome and overall accuracy. Performance change is measured by the improvement (if any) when the current feature is included in the model. Table S1: Computation complexity of machine learning classifiers using Big O notation. is the size of training data, k is the number of selected features, and K is the number of nearest neighbors in KNN. '-' implies that the computation complexity is not applicable or cannot be directly expressed [28,29].

Conclusions
This study proposed a kernel-like detection algorithm for classifying acoustic time series into cough, wheezing, and breathing. It has a novel distance measure that considers each feature's unique properties, including the class variance, probability distribution, and feature importance. Results showed that the proposed method achieved better detection accuracy than existing algorithms, and the proposed feature reduction method effectively reduced the feature dimensions to 4 from 72 and, at the same time, improved the overall classification by 13.24%. In wearable device applications, feature dimensionality is of concern as a smaller number of features requires shorter computation time for feature extraction and classification (faster detection speed) and hence less consumption of battery power (longer working hours). From our findings, the proposed method potentially tackles the current challenges of limited processing power by minimizing the feature extraction efforts and detection accuracy in health-monitoring wearable devices by capturing the unique distribution characteristics in different features and sound classes. This is also exemplified in this work through the implementation of the proposed algorithm in an embedded platform based on the IOT chipset from the Renesas synergy series. Furthermore, the findings in this study may be applied to other health-monitoring applications, including the quantification and severity analysis of clinical symptoms.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/s22062167/s1, Figure S1: Variances of features by class. For easy comparison, the variances of breath are scaled to 1 for all features; Figure S2. The distribution patterns of features by classes. Bar shades are the relative frequencies for each class and curves are the fitted probability density distributions. Figure S3. Bimodal and multi-modal distributions of features in three classes. Bars are the relative frequencies for each class and curves are the fitted probability density distributions. Figure S4. Sensitivity analysis of features in each outcome and overall accuracy. Performance change is measured by the improvement (if any) when the current feature is included in the model. Table S1: Computation complexity of machine learning classifiers using Big O notation. n is the size of training data, k is the number of selected features, and K is the number of nearest neighbors in KNN. '-' implies that the computation complexity is not applicable or cannot be directly expressed [28,29].