Participants

94 participants were recruited to take part in this study, following the ethical approval from the Institutional Ethics Committee (IEC) of the Indian Institute of Technology Madras, Chennai, India (IITM—IEC Protocol No. IEC/2018–03/MP/01). The participants were students of the Indian Institute of Technology Madras. Their age ranged from 21 to 27 years. Data were collected only once (one set of 10 breath samples) per participant. Volunteers with epileptic disorder were excluded from participation. The experimental data collection was carried out between 8th and 17th January, 2019. All volunteers who participated in this study had given their written informed consent. The recorded time series data were analyzed anonymously.

Data collection and analysis

A schematic of the experimental setup is shown in Fig 2A. It consists of a mouthpiece assembled into an aluminium circular cross-section channel which housed the hot-wire probe aligned to its axis to measure the streamwise component of the turbulent exhaled flow velocity. The human subjects were allowed to exhale through their mouths into the experimental measurement setup. The nose was clipped during data recording to ensure that all the exhaled air passes through the oral cavity before entering the experimental setup. Each human subject was provided with a new disposable plastic mouth-piece to wrap their mouth around, through which the subjects exhaled. The obstruction of the tongue to the flow was avoided by placing the mouth-piece above the tongue. Data were obtained in each exhalation trial lasting about 1.5 seconds, with 10 trials recorded per subject. Each time series was recorded by sampling the voltage response at 10kHz. This effectively gives us 15000 data points in a time series, the relevance of which would be discussed in the following sections. The time series signal from a typical exhalation trial is shown in Fig 2B. In our study, we investigated the multifractal properties of the time series, since interestingly, human exhaled breath has been found to display multifractality, based on our analysis which is discussed in detail in Part 1 of S1 Text. This was performed using the well-known technique called multifractal detrended fluctuation analysis (MFDFA) developed by [30].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Experimental setup and recorded time series.

(A) Depiction of the experimental setup for data collection. It consists of a disposable mouth-piece, a mouth-piece mount housing a hot wire anemometer and a data acquisition system. (B) A typical human exhalation velocity signal measured using a standard hot wire anemometer. The time signals were sampled at 10kHz for 1.5 seconds.

https://doi.org/10.1371/journal.pone.0301971.g002

Given a set of time series signals from a library of users, our algorithm comprises of segmentation, normalization, feature extraction and subdivision of feature set into training and testing sets. The training dataset became part of the enrolled database, whereas the testing dataset was used for testing the performance of the authentication algorithms. The enrollment and algorithm testing depends on the type of algorithm being used. More details of user authentication systems are discussed in section titled User confirmation algorithms (page 8).

Time series segmentation, normalization and selection

Segmentation of time series is a standard practice in many data analysis techniques to obtain dividing points on a signal with or without stationarity. In machine learning problems with limited availability of time series samples, segmentation is of vital importance. By performing an efficient segmentation on the basis of certain statistical measures, we can obtain sufficient number of samples to train and test machine learning models. Fig 2B is a plot showing the instantaneous voltage response from the hot wire probe for 1.5 seconds. It was obtained by sampling at a frequency of 10kHz, giving us a sufficiently resolved long series to perform segmentation without losing any significant information on the flow physics.

We will now discuss the segmentation process. Each time signal was divided into 19 overlapping segments using a window size equal to 1/10th the length of the signal and a sliding width of half the segment size. Remember the machine learning models may tend to overfit the training data when there are large number of overlapping segments. The purpose of using overlapping windows was to capture the end effects of the time series segments during feature extraction. So, the chosen segment width and sliding width are justified as each part of the time signal appears only in two segments. This effectively gives 1500 data points to each segment making it sufficiently long to reliably extract features using tools discussed in this manuscript. As a result, a maximum of 190 representative time blocks become available for the analysis for each user. Each of the time signals were normalized before feature extraction, making the time series comparable across realizations. This would also make all signals independent of the sensor being used for the measurement, since these features only rely on the temporal correlation structure in the series and not on the raw data values. This approach can be termed as being sensor-agnostic. Regardless of whether the time series signal is measured using a hot-wire/film probe, or a laser-based technique, the performance of the algorithm will not be affected, as long as there are sufficient data points to properly capture the temporal structure in the flow. We then build an algorithm which works with these features which are invariant to the absolute value of the time series. z-score normalization which is popularly known as standardization was used to normalize the time series. To perform z-score normalization, the mean of the entire time series is subtracted from each data point in the time series. Then, the resulting values are divided by the standard deviation of the time series. This scales the data so that it has a mean of zero and a standard deviation of one. The resulting normalized time series will have values that represent the number of standard deviations away from the mean. The z-score normalization has the form shown in Eq 1. (1) where z(i) is the normalized time series, x(i) is the original time series of length N, (μ_t) is the mean of the time series, and (σ) is the standard deviation of the time series. The time series becomes unitless after normalization.

MFDFA was performed on all normalised time series, and it revealed that not all spectra exhibit the expected shape. The general shape of a multifractal spectrum is convex or more precisely an inverted parabola, with the peak occurring at the central moment. This convex shape signifies the presence of multifractal scaling, indicating that different parts of the time series exhibit distinct scaling behaviors. Certain time segments were observed to result in a spectrum with folds or distortions. Fig 3 shows an example of such a distortion. The multifractal spectrum for a time signal and three randomly chosen segments X, Y and Z from the same time series are displayed. Fig 3A shows the entire time signal and the chosen segments. Out of the three segments, X and Z show a typical spectral shape, whereas segment Y consists of a fold towards the left hand side of the spectrum (see Fig 3B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Multifractal spectra for different segments of a time signal.

The multifractal spectra corresponding to the entire time signal (maroon) and time segments X, Y and Z (black, bounded by gray band) in (A) are shown in (B). It is evident that few segments exhibit an inverted parabola shape and spectrum B has a distortion.

https://doi.org/10.1371/journal.pone.0301971.g003

There could be several reasons for the appearance of folds in the multifractal spectrum. (i) They could occur due to irregularities or data artifacts in the time series itself, such as noises, outliers, etc. which may arise due to inconsistent exhalation by the user during data acquisition. For example, during the period of 1.5 seconds, if the user exhales abruptly for the first 1 second of the trial, and then the breath velocity steadily decays for the remaining 0.5 seconds. The segment which falls between these two regions might contain irregularities within it. Such irregularities could introduce inconsistencies in the scaling behaviour. (ii) The spectrum may be affected by the non-stationarity of the time series, which is when the statistical properties change with time, such as due to change in breath velocity. (iii) Spectral folds might even arise due to the finite size of the time segment. Limited number of data points may not capture the scaling properties at different scales. Investigating the type of distortions or the reason behind this behaviour of the spectrum for certain time segments fell outside the scope of this work. Instead, we made use of this behaviour as an indicator to judge whether a segment is valid or not. All segments which showed non-convex singularity spectra were discarded in our analysis. Also, the segments which produce a spectral width less than 0.05 were rejected, since they exhibit a very low degree of multifractality. These two strategies effectively make MFDFA a tool for time series selection, for further feature extraction and analysis. Any time signal which contains significant number of segments with inconsistent scaling behaviour can be rejected using this tool during the data recording step itself. A numerical example discussing how a multifractal singularity spectrum can have non-convex shapes can be found in [31].

Feature extraction

Features were extracted from normalized time signals using various time series feature extraction techniques. Unlike other physiological biometric systems where image-based patterns or features are used as templates to match an individual’s identity, our input data is a time series from an individual, which requires feature extraction. Several features of the time series were studied in order to develop insights into the data. The multifractal spectral information was incorporated into our analysis by including them in the set of features. The fact that the time series contains information pertaining to the correlation structure becomes relevant to machine learning algorithms. In keeping with this principle, we extract a set of three important features from the spectrum: (i) β, the abscissa corresponding to the spectral maxima, (ii) ω, the width of the spectrum, and (iii) ϵ, the bias or asymmetry parameter of the spectrum. The parameters β, ω and ϵ are dimensionless. These features are visualised on the multifractal spectrum of an exhaled breath time signal in Fig 4. It was also noted from our analysis that the spectra showed clear differences in their temporal structure; i.e., parameters such as β, ω and ϵ were different for different time signals. Several other multifractal spectral features have also been considered in the literature [31–33]. We chose these three features for simplicity and also they encompass the most important descriptions of a multifractal spectrum. Investigating how unique these features behave is of interest to this work.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. The multifractal spectrum.

Plot of the spectrum of singularities f(α) against the singularity strength α, computed for an exhalation time series segment. The parameters β, ω and ϵ are the features that characterize a multifractal spectrum.

https://doi.org/10.1371/journal.pone.0301971.g004

In addition to the use of MFDFA as a feature extraction algorithm, we also use an automated time series feature extraction algorithm named tsfresh (Time Series FeatuRe Extraction on the basis of Scalable Hypothesis tests) developed by [34]. The tool generates over 700 time series features using 63 different time series characterization methods. The following discussion pertains to the preparation of dataset for model building, training and testing of the algorithms. A consolidated pipeline of the algorithm towards model library building including time series normalization, and selection, followed by feature extraction and reduction, is shown in Fig 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Flow chart of the algorithm.

Flow chart showing the algorithm pipeline, including time series normalization, filtering, feature extraction, feature reduction, and data splitting into training and testing. The time signal shown here is one of the segments of the original time series. Note that the representation of blue bar for training dataset and green bar for testing dataset will be consistent in further discussions in this manuscript. The training data of all users were used for building ⁿC₂ binary classifier models, which becomes the process known as enrollment.

https://doi.org/10.1371/journal.pone.0301971.g005

Features extracted by these algorithms from all available time series are concatenated and passed through a low-variance filter. This was done to remove those feature columns with a variance value below a given threshold, which in our case was 1%. The rationale behind applying this low variance filter was to eliminate features that exhibit very little variation across instances. Such low-variance features may not provide useful insights for classification tasks. Furthermore, highly correlated features were removed from the feature set. A correlation threshold of 80% was chosen for this purpose. Removing features by these techniques reduce the dimensionality, simplifies the model, and potentially improves model performance by focusing on more informative features. All features which were derived from the absolute values of the time series, such as maximum/minimum values, quantile information, etc., were disregarded. For example, inclusion of mean value of a signal will bias the algorithms and allow them to classify on the basis of the mean values itself, which was undesired. It was observed that different human subjects were able to exhale in different velocity bands depending on their lung capacity. The filtered feature matrix thus obtained is a stack of vectors from each time series sample available, and it consisted of approximately 450 time series features. This feature space is high dimensional and may contain redundant features that can be excluded. The reduced feature set will also reduce the computational complexity of the algorithms. We adopted a feature selection method using binary random forest classifiers. Binary classifiers were built on pairwise combinations of the users’ feature datasets. The importance of the features can be quantified for every random forest binary classifier by estimating how much the random forest’s performance would suffer if a given feature were to be eliminated. This impurity-based feature importance developed by [35] was used for picking the top features. The top 10 most prevalent features among users were chosen as the feature space after computing the top 10 features from each classifier. In the later sections of this manuscript, the methods used for model construction and the physical insights of these features will be described. The reduced feature matrix thus obtained contains features of all the users in the database. For each user, the dataset was split into training (60%) and test (40%) sets. It is important to note that this splitting was done after shuffling between groups of features corresponding to the 19 time blocks for each subject. We know that there were 190 time signals for each user in the database with each set of 19 signals coming from a single recorded time series (see subsection titled Time series segmentation, normalization and selection (page 5)). Shuffling without grouping would result in the same information being spread across the training and testing dataset, which was undesired. By doing this, we made sure that out of 10 exhaled breath samples, 6 became part of the training set and 4 became part of the test set. The training feature set was used to build the model library and the test feature set was used for user confirmation/identification tests.

Building of model library

We have formulated the multi-class classification problem into a series of binary classification problems. Several studies have explored the application of pairwise binary classifiers for handling multi-class problems. A description of this technique which is also known in literature as class binarisation or round robin classification can be found in [36]. [37–39] are few others who have studied class binarisation for multi-class classification. In order to perform tests with a machine learning based algorithm, it was necessary to build binary classifier models using binary combinations of the training datasets and these models were stored in a model library. Computational simulations were setup to evaluate the performance of the user confirmation and identification algorithms. Let us briefly see how the model library grows with the addition of users to the existing database of users. This is known as enrollment mode of the biometric system. Say, there are n disjointed users U₁, U₂, …U_n in the current state of the users’ database. ⁿC₂ binary classifier models can be built, which makes up the complete model library. With the addition of a user, the updated size of the users’ database becomes n + 1. Therefore, the size of the model library increases by n and becomes ⁿ⁺¹C₂. This growth can be expressed as (2)

This means that when a new user is added to the users’ database, n additional binary classifier models are to be built and stored in the model library. Expectedly, this follows a second-order power-law variation of the form y = ax^m with the multiplication factor a ≈ 0.5 and exponent m ≈ 2.

Confirmation algorithm based on hypothesis testing

The use of hypothesis testing as an instance-based binary classifier has been attempted in the literature. [41] compared the machine learning approach and the statistical testing based on p − variations; and the idea of instance-based classification by hypothesis testing was investigated by [42]. [43] provided a detailed description on how binary decision problems can be formulated as hypothesis testing and/or binary classification. In a system based on hypothesis testing, the library comprises of the training datasets of all the users. Since we are building an algorithm which is intended to work alongside a machine learning algorithm, we formulate the hypothesis test based algorithm to work on binary pairs of users. To be more precise, the library will comprise of training datasets of pairs of users. It will be referred to as user-pair data in further discussions. Fig 6 shows a flow chart of the user confirmation algorithm which is based on hypothesis testing principles. The equality-of-means test was performed between a test data and each training data in pairs present in the library to infer whether the null hypothesis is to be rejected or not, as depicted in Fig 6. Here, the null hypothesis states that the two samples come from the same distribution (H₀: μ_a = μ_b), and the alternate hypothesis states that the samples come from different distributions (H₁: μ_a ≠ μ_b). A detailed description on the test statistic and formulation of the Hotelling’s T² test can be found in the original work by [40].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. User confirmation algorithm based on hypothesis testing.

A flow chart of the user confirmation algorithm based on hypothesis testing. The user confirmation block will be made use in the user identification algorithm later in this manuscript. An example of the hypothesis test against user-pair is illustrated inside the dotted box, directed from the user confirmation block by the red asterisk. Given a user i, the user confirmation block’s output was reposed to answer the question “Are you indeed User i?” based on a threshold.

https://doi.org/10.1371/journal.pone.0301971.g006

When a test user, say ‘User i’ was to provide the input, the pairwise Hotelling’s T² tests are performed between the test user’s data and the training data of n − 1 pairs of users which include ‘User i’, where n is the number of users in the database. Let us look at one of those tests as shown inside dotted box in Fig 6. By performing a hypothesis test against a user-pair, for example, (1, 2), we get a pair of p − values, (p₁, p₂). The tests were performed with a confidence level of 99.9%, and therefore a p − value of 0.001 or less was sufficient to reject the null hypothesis. At least one of the two p − values need to be above 0.001 for the algorithm to accept the null hypothesis. The predicted user is then the user corresponding to a higher p − value. If both p − values are either equal to or below 0.001, no predictions were made. After the test, the predictions made here are reposed as an answer to the question “Is it User i? (Yes/No)”. The pipeline discussed so far becomes the ‘User Confirmation Block—HT’ for the hypothesis testing based algorithm. The output of this block is a scalar v which is equal to the count of model predictions which says ‘Yes’. Here, a threshold of 50% of the predictions was used for defining the minimum confidence of confirmation. This means that HT(i, i) accepts the null hypothesis and HT(i, j) ∀j = 1, 2, …n and i ≠ j rejects the null hypothesis in at least 50% of the cases. Then, the User i is so confirmed. Here, HT(i, j) stands for the hypothesis test between a User i and User j.

The equality-of-means test can actually be viewed from two perspectives: (a) Testing the distribution of test data against the distribution of n training data; (b) Testing the distribution of test data against the distribution of training data in pairs as discussed so far. The former strategy produces n test results and the algorithm would face one of three scenarios: (i) If only one test accepts the null hypothesis, the user identity is presumed to be of the user corresponding to that particular test; (ii) If more than one tests accept the null hypothesis, the user corresponding to the test which corresponds to highest p-value is presumed to be of the user identity (predicted user). In either case, if the predicted user matches with the test user, the user is confirmed, otherwise not; (iii) If all tests reject or no test rejects the null hypothesis, then the user is not confirmed. Although the former case (procedure (a)) is a computationally simpler formulation, the latter case (procedure (b)) becomes more relevant in our study since we are trying to build a multi-model approach for user identification. It was also noted that the latter approach gave better confirmation results (for UCA.HT) compared to the former approach.

Confirmation algorithm based on machine learning

Following the discussions from subsection titled Building of model library (page 8), generating ⁿC₂ binary classifiers is necessary to handle the multiclass problem. The choice of a classifier depends on the specific characteristics of the dataset. A detailed discussion on the model-building procedure and the choice of a binary classifier can be found in Part 2 of S1 Text. Based on this analysis, we chose random forest as the appropriate binary classifier model for the model library. For the rest of this work, we will employ random forest as our machine learning algorithm and report results from this tool for both user confirmation and user identification.

Once the model building was complete and the entire library was stored, the test user data were given as input, say ‘User i’. The algorithm selects those models from the library which were built using the same test user and makes predictions using each model as depicted in the flow chart in Fig 7. The predictions made here are answers to the question “Is it User i? (Yes/No)”. The pipeline discussed so far becomes the ‘User confirmation block—ML’. The output of this block is a scalar v which is equal to the count of model predictions which says ‘Yes’. Here, a threshold of (again) 50% of the predictions was used for defining the minimum confidence of confirmation. This means that if the algorithm confirms the user in more than half the classification trials, i.e., when v > (n/2), the user is confirmed, else not.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. User confirmation algorithm based on machine learning.

A flow chart of the user confirmation algorithm based on machine learning. The user confirmation block will be made use in the user identification algorithm later in this manuscript. Given a user i, the user confirmation block’s output was reposed to answer the question “Are you indeed User i?” based on a threshold.

https://doi.org/10.1371/journal.pone.0301971.g007

User confirmation system

Confirmation tests were performed for all users (n) available in the database. Each set of confirmation tests were repeated 66 times by shuffling training and test data split-up. The results of the algorithm from each of these trials can be interpreted as follows: number of confirmed users denoted by c, and number of unconfirmed users denoted by u. In order to quantify the performance of the algorithms, we define a metric called the true confirmation rate (TCR) which is a ratio of the confirmed users and total number of users as shown in Eq 5. (5)

The confidence of confirmation (η) for a user confirmation algorithm is the percentage prediction of the favourable user during a confirmation test. It directly quantifies how confident the algorithm is while attempting to confirm a user i. It can be defined as, (6) where, v_i is the favourable user predictions as seen in Figs 6 and 7, i.e., the total number of model predictions that matches the user that the algorithm is attempting to confirm, and n is the total number of users in the database. The value of η_i ∀i = 1, 2, …n has to pass a threshold confidence of confirmation, say η_t, for a user to be confirmed. A comparison of the histogram of η_i is shown in Fig 9 for one trial of n confirmation tests. The study revealed that the machine learning based algorithm performs better than the hypothesis testing based algorithm. This validates the ability of a random forest classifier to capture the decision boundary better, when compared to its hypothesis testing based counterpart. For the UCA.HT, the TCR was 50±9.6%, whereas, for the UCA.ML, the TCR was 97±2.5%. This implies that almost every user was able to pass the threshold of 50% in the machine learning based algorithm. This signifies that the algorithm achieves a greater level of confidence while confirming a user using UCS.ML.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Comparison of the confidence of confirmation η_i.

Histograms of confidence of confirmation η_i compared between (A) a machine learning based approach (random forest classifiers) and (B) a hypothesis testing based classification approach, for one trial of n confirmation tests. In the example shown here, the predictions from ML classifiers give a range of η_i values distributed between ≈38% to 100%, whereas the predictions from HT based classifiers produce η_i values only close to 0% and 100%.

https://doi.org/10.1371/journal.pone.0301971.g009

We shall now investigate why the machine learning based classification algorithm performs better in comparison with a hypothesis test based classification. In the case of hypothesis testing, we know that the rejection of null hypothesis is based on the confidence level chosen. The confidence level can be visualized as a demarcating hyper-surface between two n-dimensional normal distributions. For simplicity, let us have a look at the decision boundaries captured by the random forest classifier and the hypothesis test based classifier in a chosen two dimensional feature space. Fig 10 shows a visualisation in the (β, ω) plane for a randomly chosen user-pair. The blue and red markers are the training data points corresponding to two user classes, respectively. The class regions are computed using a structured synthetic dataset in the feature space.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Comparison of the decision boundaries in (β, ω) plane.

Decision boundaries captured by (A) random forest classifier and (B) hypothesis testing based classifier for a randomly chosen user-pair. The scattered points are the training data points with red and blue labels denoting their true classes respectively. The line separating the two contour regions is the decision boundary. Accuracy of each model against the test data is displayed at the top right corner of their respective plots. The RF classifier captures a complex decision boundary compared to the HT based classifier.

https://doi.org/10.1371/journal.pone.0301971.g010

For the purpose of visualisation of a hypothesis test based classifier’s decision boundary, z − tests were performed in each dimension separately, for every data point from the synthetic dataset against one of the user’s training data. The tests were performed under the null hypothesis that the data point belongs to the distribution of the training data, under a confidence level of 99.9%. The overall null hypothesis is accepted only if the null hypothesis in both the dimensions are accepted. Comparing the decision boundaries captured by a hypothesis test based algorithm and a random forest model for the same pair of users, one can observe that the random forest model has the ability to capture a more complex decision boundary between two user classes. This lets the random forest classifier to achieve a test data accuracy of 90.9%, whereas the hypothesis testing based classifier achieves only 73.9%. Now that we have established that the machine learning based algorithm is better than the hypothesis test based algorithm for user confirmation, we will now investigate how these two algorithms perform for user identification in the following section.

User identification system

The identification algorithm discussed in Fig 8 shows that we obtain a vector V of favourable user predictions. Based on the values of vector V_j with j = 1, 2, 3, …n, we can obtain the following outcomes:

• True positives (t)—Number of users who were identified correctly.
• False positives (f)—Number of users who were identified incorrectly.
• Not identified (h)—Number of users who the algorithm was unable to identify.

We shall define the following performance metrics to evaluate the user identification algorithm:

1. Precision (P) or Positive Predictive Value (PPV), which quantifies the percentage of users who were identified correctly among all the identified users. (7)
   This parameter quantifies the probability of correct predictions given a judgement (identification) by the algorithm.
2. Accuracy (E), which quantifies the percentage of users who were identified correctly among all the users n. (8)

The precision and accuracy values computed using Eqs 7 and 8 respectively, were 35±10.5% and 29±9.1% respectively, for the hypothesis test based algorithm. The results reported in this section are in the format ‘μ_p±2σ_p’ where μ_p and σ_p are mean and standard deviation of the performance metrics respectively. For the random forest based algorithm, we were able to observe precision and accuracy values of 26±7.2% and 22±6.4% respectively. These values were computed on the basis of the maximum votes received by a user among n confirmation trials, as described previously in Fig 8. When we combine the results from both the algorithms using Eq 3 with w₁ = 0.3 and w₂ = 0.7, we get precision and accuracy values of 32±8.5% and 31±8.5% respectively. Note that the values reported here are also influenced by the threshold η_t which in this case was set to 55%. The parameters w₁, w₂, and η_t can be tweaked to make the algorithm behave on both extremes—(i) to be very liberal (low precision, low accuracy); (ii) to be very conservative (high precision, low accuracy). Taking the example of a particular trial with n = 94, for a weights setting of w₁ = 0.3 and w₂ = 0.7, η_t = 50% produces the outcomes (t, f, h) = (31, 58, 5), giving a precision of 34.8% and accuracy of 33.0%. For the same weights, η_t = 96% produces the outcomes (t, f, h) = (18, 6, 70), giving a precision of 75.0% and accuracy of 19.1%. The former case allows for a lot of false positives by making judgements on most of the instances, whereas the latter case of the algorithm makes judgements stringently.

With the right set of hyperparameters (w₁, w₂, …w_r (in the general case from Eq 4) and η_t), a multi-modal approach is expected to improve the robustness of the overall algorithm. If one classifier produces incorrect predictions for certain trials, other classifiers in the ensemble can compensate for it and provide correct predictions. The contribution of each algorithm can be controlled by the weights. This robustness helps in improving the generalization of the ensemble model. The following discussion is based on results produced from this combined algorithm. We know that the highest voted user becomes the identified user from the algorithm. Based on the 66 shuffle trials, we have the following understanding of the user database. 21.3% to 42.6% of the users can be correctly identified by them being the highest voted users, 39.4% to 57.4% of the users can be correctly identified as at least the second highest voted users, and 50.0% to 66.0% of the users can be correctly identified as at least the third highest voted users. This is remarkable given that it is the first attempt in the literature to classify and uniquely identify individuals based solely on the fluid physics of the exhaled breath. We believe that this is conclusive evidence that the fluid dynamic structure of the exhaled breath contains uniquely identifiable information.

This algorithm holds tremendous potential for future use in the area of personalised medicine and also as a novel way to store biological data. This can be achieved by careful model selection and generalisation of classifier models. Advanced models such as deep neural networks can be made use to enhance the multi-model approach discussed in this manuscript.

Physical insights: Understanding the defining features

In order to make a physics-based argument for the uniqueness of human exhalation, it is important to investigate the physical significance of the most important features that result in robust classification. These would be the set of features or attributes which inherently differentiate the classes for a given training data. As we have seen, the importance of the features were quantified for every random forest binary classifier for choosing a reduced feature set in subsection titled Feature extraction (page 6). These features are to be investigated to understand their physical meaning in the context of the current problem in hand. A description of the most important classifying features (in the decreasing order of importance) are as follows.

1. The singularity strength or Hölder exponent corresponding to the maximum (β) of the multifractal spectrum of the exhaled breath time series: This is a feature extracted using the MFDFA. β explains the long range correlation present in the time series. A low value indicates that the underlying process becomes correlated and loses fine structure, becoming more regular in appearance [30]. This, in our case, would relate to the organised motion of vortical structures in the turbulent exhaled air flow. For some subjects the vorticity pattern might be more irregular than the others, which could be attributed to the extrathoracic morphology.
2. The sum over the absolute value of consecutive changes in the velocity time series: This feature represents the total magnitude of absolute differences between successive data points. In the context of our study, a higher value of this feature indicates a greater overall change in velocity between consecutive data points, i.e., the velocity changes rapidly and frequently. In contrast, a low value of this feature indicates that the velocity is smooth and consistent. It provides a quantitative measure of how much the velocity values fluctuate over successive time intervals, which in our case is 0.1 milliseconds. The detection of distinctive patterns in these fluctuations can provide insights into the presence of vortical structures in exhaled breath flow, contributing to the uniqueness of these patterns for individual subjects and enabling their classification by the algorithm.
3. Third coefficient of the autoregressive AR(r) model with order parameter r = 10: The parameter r is the maximum lag of the autoregressive process. The AR model generally predicts future behavior based on past data. The importance of the third as well as fourth (point 8) coefficients show that there is some correlation between successive values in the time series for most of the users.
4. The number of peaks in the time series with a support (s) of at least 1: A peak of support s is defined as a sub-sequence in the time series where a value occurs that is greater than its s neighbors to the left and to the right. When s is set to 1, this feature computes the number of peaks in the time series where a value is greater than its immediate neighbors. This feature can provide insights into the presence or intensity of localised fluctuations in the flow.
5. The number of different Continuous Wavelet Transform (CWT) peaks present in the signal for smoothing width of 1: This feature was extracted from the time series by applying CWT using Ricker wavelet with width, w = 1. This method simultaneously evaluates the signal in the temporal and frequency domains. In the context of our study, the identified CWT peaks represent distinctive features in the breath signature. Physically, these peaks may correspond to specific events or patterns that are characterised by rapid changes in both time and frequency domains. For instance, a CWT peak could signify the presence of a sudden, localized change in the breath velocity with a particular frequency content. The number of distinct peaks across the considered width scales provides a quantitative measure of the breath signature’s complexity. It can be utilized to compare the signals based on their peak characteristics.
6. The value of partial autocorrelation function at a lag of 3: The partial autocorrelation is a statistical measure that quantifies the linear relationship between a time series variable and its lagged values. In the context of our exhaled breath flow, the partial autocorrelation can provide insights into the temporal dependence and correlation structure of the breath velocity. This means that this feature can be useful in understanding the persistence or memory of the signal. It suggests that a strong linear relationship between the current flow state and its state 3 time steps ago have been important for the classification of human subjects. In our analysis, a ‘time step’ corresponds to the original sampling rate of 10kHz. Therefore, when we refer to a lag of 3 time steps, it signifies a duration of 0.3 milliseconds.
7. Width of the multifractal spectrum (ω) of the exhaled breath time series: ω describes the richness of the multifractality present in the time series, i.e., wider the range of singularity strength, richer the structure of the signal. The spectral width can implicitly represent the intensity or the level of turbulence present in the flow of exhaled breath. Turbulence is characterized by fluctuations in velocity at different scales. A wider range of turbulence scales is reflected by a wider spectral width, indicating a more turbulent flow. This might be attributed to factors such as extrathoracic constriction, or increased turbulence due to specific breath patterns or breath dynamics.
8. Fourth coefficient of the autoregressive AR(r) model with order parameter r = 10.
9. The number of different continuous wavelet transform (CWT) peaks present in the signal for smoothing width of 5. This feature was extracted using the same technique as discussed in point 5, but with a width of w = 5. A larger smoothing width typically leads to a broader wavelet. A wider wavelet provides a smoother analysis that might emphasize broader features and lower-frequency components in the signal. Conversely, a smaller smoothing width of w = 1 (point 5) would result in a narrow wavelet, allowing for a more detailed examination of rapid changes in the signal (sensitive to high-frequency components).
10. Kurtosis of the velocity time series calculated with the adjusted Fisher-Pearson standardized moment coefficient, g2: We know that Kurtosis is a higher-order statistical attribute of velocity signals. The heaviness of the tails of the probability density functions of normalized time series could be distinct for each user. This feature will help us in assessing the degree of deviation from the Gaussian distribution and provides evidence of skewed behaviour of the time series.

Computational complexity of the algorithm

Run-time of an algorithm is an extremely important factor for a real-time biometric system. It was generally observed that the size of the input feature set affects the amount of computational resources required to run an algorithm. It was observed that the hypothesis test based algorithm performs predictions faster than the machine learning based algorithm which is because the former is an instance-based classifier. Since the user identification algorithm depends on the number of users and in turn the number of models in the model library, the identification time per user was expected to scale up with the size of the library. The identification time was observed to show a linear relationship with the size of the library (of the form y = ax, with slope a ≈ 1) as seen in the Fig 11. The error bars show 95% confidence interval at every data point.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Dependence of user identification time on the size of model library.

Plot showing the linear relation of user identification time with the growth of model library. This is applicable to the ML based algorithms which include building of binary classifier models (also known as enrollment in the context of biometrics). The error bars show 95% confidence interval at every data point.

https://doi.org/10.1371/journal.pone.0301971.g011

One of the advantages of building an algorithm which uses ⁿC₂ binary classifiers instead of a single multi-class classifier is that it is massively parallelisable. As long as we have sufficient number of cores to run model loading and prediction, the parallelisation is possible. This significantly improves the computational time by several orders.