1.

Introduction

Skin cancer is the most common type of cancer in the United States, with melanoma being the fifth most prevalent among men and women.¹ The 5-year survival rate of patients with early-stage skin melanoma is $\sim 94\%$; however, $\sim 13\%$ of skin melanoma patients are diagnosed with lesions already at intermediate or advance stages,¹ which are associated with 5-year survival rates of $\sim 61\%$ and $\sim 27\%$, respectively. The most common diagnosis strategy for skin cancer is clinical evaluation of suspicious lesions followed by biopsy for histopathological evaluation to confirm diagnosis and tissue staging. One major drawback of this practice is the inability to clinically distinguish between similar lesions; in particular, melanoma is often mistaken for other benign pigmented lesions such as seborrheic keratosis (pSK). In addition, it is known that the accuracy of melanoma diagnosis with unaided eye is only about 60%.² Therefore, clinical tools that could provide objective, in situ, and accurate noninvasive discrimination between malignant and benign skin lesions during clinical examination could significantly improve early detection of skin cancer, reduce the risk of adverse events, and lead to improved cost-conscious patient care.

One of the most common tools used by physicians to diagnose skin cancer lesions is the dermoscope^2–4 which helps the unaided eye by magnifying the features on the skin. This allows doctors to examine the morphological features of concerning lesions at a significantly more detailed level. Although dermoscopy is known to improve the diagnostic sensitivity of skin lesions by $\sim 10\%$ to 30%, its performance largely depends on both the level of experience of the dermatologist and the type of lesions.² The highly subjective nature and poor reproducibility of this method have led to the emergence of several proposed computer-aided diagnostic (CAD) systems.^3,5–7

CAD systems are becoming largely popular in both diagnosis and prognosis of various diseases as they allow automated and noninvasive analysis of the tissue conditions. Table 1 summarizes some of the published works that reports the diagnosis and classification of pigmented skin lesions. Most of the works used dermoscopic images that were either collected by the authors or from publicly available datasets (e.g., ISIC archive, ISBI, Atlas, HAM10000, or PH2). Celebi and Zornberg⁸ explored the clinically significant colors in dermoscopic images using K-means clustering and employed symbolic regression to classify the lesions. Ramlakhan and Shang⁹ designed a melanoma recognition system using smart phone photographs that are classified using k-nearest neighbor (kNN) algorithm. Satheesha et al.¹⁰ examined computerized three-dimensional (3D) dermoscopy features of skin cancer lesions to develop multiclass classifiers using Adaboost, bag of features (BoF), and support vector machine (SVM) techniques. Khristoforova et al.¹¹ used logistic regression to classify benign and malignant skin lesions using spectral features from Raman and autofluorescence spectroscopy measurements.

Table 1

Summary of previously reported works on pigmented skin lesion classification.

Classification of dermoscopy images of benign and malignant skin lesions using different deep learning approaches has also been reported. Harangi¹³ used an ensemble of different convolutional neural network (CNN) classifiers, while Romero Lopez et al.¹⁴ used transfer learning with pretrained VGGNet CNN architecture. Majtner et al.¹⁵ combined CNN with SVM classifier using handcrafted RSurf features and local binary patterns to classify melanomas from other benign skin lesions. Lee et al.¹⁶ developed the WonDerm pipeline that segments the skin cancer dermoscopic images using neural network architectures and classifies it using an ensemble approach. Amin et al.¹⁷ extracted features using pretrained AlexNet and VGG16 deep learning architectures, performed feature selection using principal component analysis, and applied traditional machine learning models including SVM, kNN, and discriminant analysis. Jojoa Acosta et al.¹⁸ utilized transfer learning with ResNet-152 architecture to classify benign and malignant skin lesions using dermoscopic images.

It has been widely established that autofluorescence responses of intrinsic fluorophores vary significantly between normal and neoplastic tissues.^19–23 Neoplastic progressions in the epithelial tissue are associated with morphological, biochemical, and functional alterations which can cause changes in the autofluorescence responses from the tissue.^22–24 The skin has several intrinsic fluorophores, such as nicotinamide adenine dinucleotide (NADH), flavin adenine dinucleotide (FAD), collagen, elastin, keratin, melanin, and porphyrins.^19,25–27 The levels of two metabolic cofactors and endogenous fluorophores in the epidermis, the reduced-form NADH and FAD, can change as skin cancer develops.^19,24 The optical redox ratio, typically defined as the ratio of fluorescence intensity of NADH to FAD, is sensitive to changes in the cellular metabolic rate. Increased cellular metabolic activity, a hallmark of neoplastic cell transformation, is usually attributed to a decrease in the optical redox ratio. In addition, the fluorescence lifetimes of these metabolic cofactors are sensitive to protein binding, thus to cellular metabolic pathways involving NADH and FAD. As a result, carcinogenesis process has been shown to cause changes in both NADH and FAD fluorescence lifetimes. Finally, cancer development also leads to extracellular matrix remodeling occurring within the dermis, which together with concurring epidermis thickening, result in a decrease in connective tissue autofluorescence that can be measured. Therefore, interrogation of NADH, FAD, and collagen autofluorescence could provide optical biomarkers of skin epithelial cancer.

Preferential excitation of these endogenous fluorophores in the tissue by multiple excitation sources could shed light on the biochemical changes in the target lesion area.²⁸ The broad emission spectral bandwidth of the fluorescence intensity signal has an intrinsic disadvantage in that it is difficult to differentiate between the intensities of overlapping emissions from multiple fluorophores. Time-resolved technique such as multispectral autofluorescence lifetime imaging (maFLIM) overcomes this challenge by quantifying the fluorescence lifetime in addition to the emission spectrum. Alex et al.²⁹ demonstrated fluorescence lifetime-based imaging of minipig skin and human skin to specifically target the endogenous fluorophores: keratin, NADH, melanin, elastin, and collagen under 725-nm multiphoton excitation. The capability of such optical biopsy techniques to serve as promising tools for dermatological research to facilitate preclinical and clinical translation is also highlighted. Huck et al.³⁰ demonstrated the effectiveness of the combined modality, multiphoton-based intravital tomography and fluorescence lifetime imaging to monitor the progression of inflammatory skin diseases. The two modalities studied the biochemical changes induced by the redistribution of mitochondria at different stages of inflammatory skin conditions.

Several animal and human tissue studies have been published on the autofluorescence properties of skin cancer lesions.^31–35 Pastore et al.³⁶ conducted experiments with mouse models to study the autofluorescence response from melanoma skin lesions using multiphoton excitation at 740 and 900 nm and emission spectral bands at 447 and 540 nm. A significant difference in the bound and free NADH ratio between cancerous and noncancerous sites was observed, while the fluorescence decay obtained from targeting FAD remained almost the same between the two regions. It was also mentioned that the presence of melanin in the deeper layers of the skin tissue could interfere with the overall fluorescence response from the lesions. Miller et al.³² studied the autofluorescence emission properties between squamous cell carcinoma (SCC) bearing and normal mice skin under 480-nm excitation, and a decrease in the short lifetime component for SCC in comparison to normal skin was observed for 535-nm emission band. Drakaki et al.³⁷ studied the autofluorescence responses from mouse, chicken, and pig skins under ultraviolet (UV) excitation, and the structural differences and the variations in tissue constituents were investigated between the different animal species for an emission spectral band between 340 and 950 nm. De Beule et al.³⁴ investigated the autofluorescence response from ex vivo biopsy skin lesions under 355- and 440-nm excitations, and the average fluorescence lifetime was found to be useful in discriminating basal cell carcinoma (BCC) from normal skin tissues at the emission band between 390 and 600 nm. Galletly et al.³¹ imaged unstained human biopsy samples using maFLIM under a 355-nm pulsed laser excitation, and significant differences in the mean fluorescence lifetimes for the emission wavelengths 375 and 455 nm were observed between the autofluorescence responses from BCC skin lesions and healthy skin. Lohmann and Bodeker³⁸ analyzed the fluorescence intensities at the emission wavelength 470 nm, from human skin with melanoma, nevi, and dysplastic nevi lesions under 365-nm excitation, and a significant difference in fluorescence intensities was observed for melanoma and nevi lesions, while melanoma and dysplastic nevi lesions did not show much difference. Fast et al.³⁹ investigated the autofluorescence response from human skin at 780-nm frequency doubled excitation and two emission channels at 535 and 720 nm corresponding to red and green channels. Red channel collects fluorescence emission from melanin, while the green channel collects emission from keratin, NAD(P)H, FAD, and elastin.

In this work, we developed and validated a label-free and objective image-guided strategy for the clinical evaluation of suspicious pigmented skin lesions based on maFLIM dermoscopy. In addition, a computationally efficient frequency-domain deconvolution of maFLIM data is explored, and three different pools of global image-level maFLIM features were evaluated for machine-learning (ML)-based objective discrimination between malignant and benign pigmented skin lesions.

2.

Methods

A summary of the complete methodology performed in this study is shown in Fig. 1.

Fig. 1

Summary of methodology showing maFLIM image acquisition, preprocessing, feature extraction, and classification. maFLIM, multispectral autofluorescence lifetime imaging.

2.2.

maFLIM Data Preprocessing

Pixel-level preprocessing: The maFLIM data measured at each image pixel $(p,q)$ are composed of three fluorescence intensity temporal decay signals $y_{m,\lambda }(p,q,t)$ measured at the three targeted emission spectral bands ($\lambda$). The preprocessing steps applied to each pixel maFLIM temporal signal is shown in Fig. 2(a). First, offset subtraction was applied to the raw maFLIM signal, $y_{m,\lambda }(p,q,t)$, followed by spatial averaging (order $5\times 5$) to increase the signal-to-noise ratio (SNR) of the time-dependent signal. The offset was estimated by fitting a straight line on the first and last five time points in each channel. The baseline was then subtracted from the entire time vector to obtain the corrected signal. Since the background fluorescence was significantly lower than the sample fluorescence, the background correction of the signals was not performed. Second, the duration of the temporal decay signals for all emission bands was adjusted to the length of the longest signal among the three emission channels, which is 149 temporal samples (59.6 ns) by applying zero padding to the short signals. Finally, the signals from the three emission channels, $y_{\lambda }(p,q,t)$, are concatenated to form $y(p,q,t)$ as shown in Eq. (1). Signal concatenation is essential for cluster analysis in image level preprocessing, explained later in this section, as well as for frequency-domain deconvolution explained in Sec. 2.3.2. The concatenated signal at each pixel location can be represented as

Eq. (1)

$$y(p,q,t)=\sum _{n=0}^2{y}_{{\lambda }_{n+1}}(p,q,(t-M.n)),$$

where $y(p,q,t)$ is the preprocessed concatenated maFLIM decay signal; $y_{{\lambda }_1}(p,q,t)$, $y_{{\lambda }_2}(p,q,t)$, and $y_{{\lambda }_3}(p,q,t)$ are the preprocessed maFLIM decay signals from each of the three spectral channels, ${\lambda }_1=390\mathrm{nm}$, ${\lambda }_2=452\mathrm{nm}$, and ${\lambda }_3=500\mathrm{nm}$; $M$ is the temporal spacing between the signals from the three channels; $(p,q)$ indicates the pixel locations. The value of $M$ is equal to 149, which is the length of the fluorescence emission decays in each channel.

Fig. 2

(a) Transformations in a single pixel multispectral maFLIM data during pixel-level preprocessing. (b) Example maFLIM image with $K$-means cluster mask and the two separated regions. The images map the total integrated intensity of the maFLIM signals at each pixel location. maFLIM, multispectral autofluorescence lifetime imaging.

Image-level preprocessing: Pixels presenting either signal saturation or low SNR ($<15\mathrm{dB}$) were detected and masked. The majority of the acquired maFLIM images contain pixels from within and outside the skin lesion region; thus, cluster analysis was performed to group pixels based on their region of origin. The concatenated signal $y(p,q,t)$ is used as the feature vector for cluster analysis at each pixel location to simultaneously include the information from the three emission channels. The steps for cluster analysis are as follows: First, an unsupervised K-means clustering algorithm was applied to generate two cluster masks. Then, each cluster mask is applied to the maFLIM image to define two regions within the FOV of the maFLIM image. It should be noted that since the K-means clustering algorithm involves random initialization of cluster centroids, it is difficult to identify which cluster mask belongs to within or outside the skin lesion region; thus, the identified regions were taken as two arbitrary regions: region-1 and region-2. Figure 2(b) shows an example of the cluster mask and the two separated regions generated from a representative maFLIM image.

2.3.

Feature Extraction

2.3.1.

Features based on time-domain deconvolution parameter estimation

In the context of time-domain maFLIM data analysis, the fluorescence decay $y_{\lambda }(p,q,t)$ measured at each emission spectral band ($\lambda$) and spatial location $(p,q)$ can be modeled⁴² as the convolution of the fluorescence impulse response (FIR) $h_{\lambda }(p,q,t)$ of the sample and the measured IRF $u_{\lambda }(t)$:

Eq. (2)

$$y_{\lambda }(p,q,t)=u_{\lambda }(t)*h_{\lambda }(p,q,t).$$

The standard method for time-domain maFLIM data analysis proceeds by first deconvolving the IRF of each spectral band ($u_{\lambda }(t)$) from the corresponding measured time-resolved fluorescence signal $y_{\lambda }(p,q,t)$ to estimate the sample FIR for each image pixel, $h_{\lambda }(p,q,t$),³⁹ which is usually modeled as a multiexponential decay. The model order (number of exponential components) can be selected by analyzing the model-fitting mean squares error (MSE) as a function of the model order. For the maFLIM data of this study, a model order of two was selected, since the addition of a third component did not reduce the MSE. The variations in error for one, two, and three exponential components during fitting is shown in Fig. S1 in the Supplemental Material. The FIR was modeled as

Eq. (3)

$$h_{\lambda }(p,q,t)={\alpha }_{\text{fast},\lambda }{e}^{\frac{-t}{{\tau }_{\text{fast},\lambda }(p,q)}}+{\alpha }_{\text{slow},\lambda }{e}^{\frac{-t}{{\tau }_{\text{slow},\lambda }(p,q)}},$$

where ${\tau }_{\text{fast},\lambda }$ and ${\tau }_{\text{slow},\lambda }$ represent the time-constant (lifetime) of the fast and slow decay components, respectively; while ${\alpha }_{\text{fast},\lambda }$ and ${\alpha }_{\text{slow},\lambda }$ represent the contribution of the fast and slow decay components, respectively. The average fluorescence lifetime for each spectral band at each pixel location is computed as

Eq. (4)

$${\tau }_{\mathrm{avg},\lambda }(p,q)=\frac{\int t{h}_{\lambda }(p,q,t)\mathrm{d}t}{\int h_{\lambda }(p,q,t)\mathrm{d}t}.$$

The parameters of the biexponential decay model are estimated for each pixel by nonlinear least squares iterative reconvolution.⁴² After deconvolution, the biexponential parameters estimated at each pixel can be used as features representing the temporal dynamics of the fluorescence decays at each emission spectral band: ${\alpha }_{\text{fast},\lambda }(p,q)$, ${\alpha }_{\text{slow},\lambda }(p,q)$, ${\tau }_{\text{fast},\lambda }(p,q)$, ${\tau }_{\text{slow},\lambda }(p,q)$, and ${\tau }_{\mathrm{avg},\lambda }(p,q)$. Since the sum of ${\alpha }_{\text{fast},\lambda }(p,q)$ and ${\alpha }_{\text{slow},\lambda }(p,q)$ is equal to one, only one of them is kept as a feature.

In addition, the following spectral intensity features can be also estimated from the deconvolved FIR, $h_{\lambda }(p,q,t)$. Absolute fluorescence intensities $I_{\lambda }(p,q)$ for each emission spectral bands are simply computed by time integrating the FIR $h_{\lambda }(p,q,t)$:

Eq. (5)

$$I_{\lambda }(p,q)=\int h_{\lambda }(p,q,t)\mathrm{d}t.$$

The normalized fluorescence intensities $I_{\lambda ,n}(p,q)$ can be also computed from the multispectral absolute fluorescence intensities $I_{\lambda }(p,q)$ as follows:

Eq. (6)

$$I_{\lambda ,n}(p,q)=\frac{I_{\lambda }(p,q)}{{\sum }_{\lambda }{I}_{\lambda }(p,q)}.$$

Lastly, the ratio of absolute intensities from the three spectral channels is computed at each pixel location resulting in three additional features:

Eq. (7)

$$I_{390,n}/I_{452,n}(p,q)=\frac{I_{390,n}(p,q)}{I_{452,n}(p,q)},$$

Eq. (8)

$$I_{452,n}/I_{500,n}(p,q)=\frac{I_{452,n}(p,q)}{I_{500,n}(p,q)},$$

Eq. (9)

$$I_{390,n}/I_{500,n}(p,q)=\frac{I_{390,n}(p,q)}{I_{500,n}(p,q)}.$$

In general, the features can be computed at the pixel or the image level. In this study, image-level global features were explored, whereby one set of features, a single feature vector, is estimated to represent the whole image. Each global feature was computed from the corresponding pixel-level maFLIM feature map as follows. Based on the two regions identified per image using the cluster analysis described in Sec. 2.2, the feature median value for each region was computed, and the absolute value of their difference was taken as the global feature:

Eq. (10)

$${\text{Feature}}_{\text{Global}}=|\text{median}({\text{Feature}}_{\text{pixel}-\text{level},\text{Region}-1})-\text{median}({\text{Feature}}_{\text{pixel}-\text{level},\text{Region}-2})|.$$

Since ${\mathrm{Feature}}_{\text{Global}}$ is the absolute difference between the feature median values from the two clustered regions, it is independent of labeling the regions as either lesion or healthy. This is particularly beneficial as it is difficult to label the clustered regions due to the unavailability of the pixel-level ground truth labels as well as the randomness in the K-means algorithm. These defined global features have the advantage of reducing patient-to-patient variability in the extracted features, which is particularly important as the color and texture of skin vary considerably with ethnicity and age. This feature extraction approach based on time-domain deconvolution of the maFLIM data generates a total of six intensity and 12 biexponential global maFLIM features, as summarized in Table 2.

Table 2

Feature set showing both intensity and biexponential global maFLIM features.

Intensity features		Biexponential features
$I_{390,n}$	$I_{390,n}/I_{452,n}$	${\alpha }_{\text{fast},390}$	${\alpha }_{\text{fast},452}$	${\alpha }_{\text{fast},500}$
$I_{452,n}$	$I_{452,n}/I_{500,n}$	${\tau }_{\text{slow},390}$	${\tau }_{\text{slow},452}$	${\tau }_{\text{slow},500}$
$I_{500,n}$	$I_{390,n}/I_{500,n}$	${\tau }_{\text{fast},390}$	${\tau }_{\text{fast},452}$	${\tau }_{\text{fast},500}$
—	—	${\tau }_{\mathrm{avg},390}$	${\tau }_{\mathrm{avg},452}$	${\tau }_{\mathrm{avg},500}$

2.3.2.

Phasor-based features from frequency-domain deconvolved signals

As mentioned in Sec. 2.3.1, the parameters of the biexponential decay model are estimated for each pixel by nonlinear least squares iterative reconvolution, which is computationally expensive and time consuming. This brings about the need to develop a much simpler algorithm for extracting maFLIM features with comparatively similar discriminative capability. An alternate fitting-free strategy explored in this work is inspired by Campos-Delgado et al.,⁴³ where a model-free representation of maFLIM data was developed utilizing the frequency-domain properties of the phasor representations. Here, we aim to replace the iterative time-domain deconvolution process by a simple division operation in the frequency domain.⁴⁴ The computational overload is further reduced by processing the fluorescence decays of all the three spectral channels together, unlike in the traditional method where the maFLIM signal for each spectral channel must be processed separately. Subsequently, several features can be extracted from the frequency-domain phasor representation of the maFLIM data.^45,46 The presented method proceeds in three steps: (1) performing frequency-domain deconvolution of the instrument response from the concatenated pre-processed fluorescence decays from all the three spectral channels, (2) constructing phasor plots for the maFLIM data, and (3) extracting global features from the phasor plots representing each maFLIM image. A detailed description of this method is presented as follows.

In this approach, the preprocessed and concatenated maFLIM signals at each pixel location, $y(p,q,t)$ are normalized to sum one for further processing and feature extraction. Similar to the concatenated, preprocessed signal, $y(p,q,t)$, the concatenated IRF from all the three spectral channels can be mathematically represented as

Eq. (11)

$$u(t)=\sum _{n=0}^2{u}_{{\lambda }_{n+1}}(t-M.n),$$

where $u(t)$ is the concatenated IRF; $u_{{\lambda }_1}(t)$, $u_{{\lambda }_2}(t)$, and $u_{{\lambda }_3}(t)$ are the IRF signals from the three spectral channels; $M$ is the temporal spacing between the signals from the three channels; and $(p,q)$ indicates the pixel positions.

The first step of the algorithm is to compute the Fourier transform (FT) of both the preprocessed concatenated signal and the concatenated IRF. The FT of the signals $y(p,q,t)$ and $u(t)$ can be represented as $Y(p,q,\omega )$ and $U(\omega )$, respectively, where $\omega$ is the angular frequency:

Eq. (12)

$$Y(p,q,\omega )=\mathrm{FT}\{y(p,q,t)\},$$

Eq. (13)

$$U(\omega )=\mathrm{FT}\{u(t)\}.$$

If the effective FIR from all the three fluorescence emission channels is denoted as $h(p,q,t)$, the effective fluorescence frequency response $H(p,q,\omega )$ is obtained from the FT of $h(p,q,t)$ as

Eq. (14)

$$H(p,q,\omega )=\mathrm{FT}\{h(p,q,t)\}.$$

Therefore, the convolution in Eq. (1) can be represented as a multiplication in the frequency domain as

Eq. (15)

$$Y(p,q,\omega )=H(p,q,\omega )U(\omega ).$$

To uniformly scale all the frequency components of $Y(p,q,\omega )$, the normalization with respect to the DC response $Y(p,q,0)$ can be applied as follows:

Eq. (16)

$$\frac{Y(p,q,\omega )}{Y(p,q,0)}=\frac{H(p,q,\omega )U(\omega )}{H(p,q,0)U(0)}.$$

Subsequently, the normalized fluorescence frequency response $P(p,q,\omega )$ can be estimated as

Eq. (17)

$$P(p,q,\omega )=\frac{H(p,q,\omega )}{H(p,q,0)}=\frac{Y(p,q,\omega )U(0)}{U(\omega )Y(p,q,0)}.$$

A phasor representation for the normalized frequency response $P(p,q,\omega )$ at specific values of the frequency $\omega$ can be generated by plotting the real $\mathrm{Re}[P(p,q,\omega )]$ versus the imaginary $\mathrm{Im}[P(p,q,\omega )]$ components of $P(p,q,\omega )$. Therefore, each pixel of the maFLIM image is mapped to a point in the corresponding phasor plot generated for a specific frequency. This transformation is shown in Fig. 3, where a representative maFLIM image is mapped to its corresponding phasor plot for an arbitrary frequency. Region-1 and region-2 marked on the maFLIM image represent the regions obtained after clustering. The pixels of each region are mapped into a two-dimensional (2D)-histogram distribution on the phasor plot, as shown in Fig. 3.

Fig. 3

Transition of a sample maFLIM image to the corresponding 2D histogram distribution on the phasor plot. Figure also shows the transformation of pixels from both regions 1 and 2 on the maFLIM image into corresponding points on the phasor plot computed at an arbitrary frequency component.

From the phasor representations of the maFLIM images at specific frequency components, $\omega =2\pi f$, the following features were computed as follows. First, a bivariate Gaussian function was fitted to the phasor distribution of each region (region-1 and region-2) of a given maFLIM image: $f_1=N({\mu }_1,{\mathrm{\Sigma }}_1)$, $f_2=N({\mu }_2,{\mathrm{\Sigma }}_2)$ (Fig. 4). The “distance” between the phasor distributions of the two regions was then estimated as the magnitude of the difference of the distribution means: $d=|{\mu }_1-{\mu }_2|$ [Fig. 4(a)]. The determinant of the covariance matrix $|\mathrm{\Sigma }|$ from the fitted Gaussian distribution provides a measure of the “spread” of the distribution. The difference in spread of the phasor distributions of the two regions was thus estimated as: $\mathrm{\Delta }\mathrm{\Sigma }=|{\mathrm{\Sigma }}_1-{\mathrm{\Sigma }}_2|$ [Fig. 4(b)]. The “angle” $\theta$ between major axes of the phasor distributions of the two regions was estimated as the acute angle between the eigenvectors of maximum variance of the multivariate Gaussian distributions (Fig. 4(b)). Finally, the “symmetry” of the Gaussian distribution can be quantified as the ratio of the variances along the orthogonal directions, $s=\frac{{\sigma }_p^2}{{\sigma }_q^2}$ [Fig. 4(c)]. The difference in symmetry of the phasor distributions of the two regions was thus estimated as: $\mathrm{\Delta }s=|s_1-s_2|$.

Fig. 4

(a) 2D histogram phasor distributions from the pixels corresponding to the two regions in an maFLIM image. The distance between the distributions is indicated by “$d$.” (b) Phasor distribution scatter plots with bivariate Gaussian fits on regions 1 and 2. The covariance matrices ${\mathrm{\Sigma }}_1$ and ${\mathrm{\Sigma }}_2$ give a measure of spread of the two regions and $\theta$ represents the angle between their major axes. (c) Phasor distribution scatter plot showing the variances ${\sigma }_p^2$ and ${\sigma }_q^2$ along the major axes. The ratio of the variances indicates the symmetry of the distribution.

The fluorescence frequency response $H(p,q,\omega )$ is bandlimited with a bandwidth of $\sim 60\mathrm{MHz}$. To cover the bandwidth of $H(p,q,\omega )$, only the first nine frequency components of $H(p,q,\omega )$ were selected, corresponding to the frequencies 5.6, 11.2, 16.8, 22.4, 28, 33.6, 39.2, 44.8, and 50.4 MHz. These frequency values are the first nine harmonics of the signal Fourier spectrum, which was calculated at a frequency resolution of 5.6 MHz (sampling frequency/#$\text{samples}=2.5\mathrm{GHz}/(3\times 149)$). For each of these frequency components, the four phasor features were computed, resulting in a total of 36 phasor features.

2.4.

Feature Selection

For maFLIM dermoscopy-based automated classification of benign versus malignant skin lesions, a simple quadratic discriminant analysis (QDA) model was explored with the global features described in Sec. 2.3. Three different pools of global features were considered: (1) intensity ($n_{\text{features}}=6$), (2) biexponential ($n_{\text{features}}=12$), and (3) phasor ($n_{\text{features}}=36$) features, where $n_{\text{features}}$ is the number of features in each feature pool. In addition to these individual feature pools, different combinations of the feature pools were also considered: phasor ∪ biexponential ($n_{\text{features}}=48$), phasor ∪ intensity ($n_{\text{features}}=42$), intensity ∪ biexponential ($n_{\text{features}}=18$), and phasor ∪ biexponential ∪ intensity ($n_{\text{features}}=54$). Feature selection using sequential forward search (SFS)⁴⁷ was performed on each feature pool independently using a leave-one-patient-out cross-validation (LOPO-CV) strategy, whereby the data of one patient are left out at each fold [Fig. 5(a)]. This assures that the left-out patient data are not used for feature selection and model training. Unlike exhaustive search where every possible combination of features is examined, SFS is computationally simpler and provides an efficient strategy to investigate the importance of the available features. The steps involved in the SFS feature selection process are shown in Fig. 5(b). The feature subset yielding the highest area under the curve (AUC) value in receiver operator characteristics (ROC) analysis was selected ($f_{\mathrm{SFS}}$) at each iteration of the LOPO-CV. Subsequently, the QDA classifier is retrained using $f_{\mathrm{SFS}}$ at each LOPO-CV iteration and then tested using the left-out patient data. A threshold of 0.5 on the generated posterior probability classifies the left-out patient data into either benign or malignant. Upon accomplishing all folds, each patient data had been used as test, and a confusion matrix was constructed to estimate the classification model performance. The number of features selected during SFS ($n_{\mathrm{SFS}}$) is kept constant in each fold of the LOPO-CV. Therefore, to determine the number of features that can produce the best classifier, the experiments are repeated for different values of $n_{\mathrm{SFS}}$, varying from 1 to 7 (1 to 6 for intensity feature pool). The maximum number of features in the feature set is chosen based on Hua et al.,⁴⁸ where the optimal number of features for feature sets with some degree of correlation is $\sqrt{n_{\mathrm{data}}}$, where $n_{\mathrm{data}}$ is the number of data points. In this case, the maximum number of features is chosen as the closest and lower integer value to $\sqrt{n_{\text{lesions}}}=\sqrt{60}=7.75\sim 7$. The value of $n_{\mathrm{SFS}}$ that produced the best $F$-score computed from the final confusion matrix was then chosen as the number of features in the final feature set ($n_{\text{selected}}$). This is because $F$-score gives a combined estimate of both the sensitivity and specificity. However, in cases when two values of $n_{\mathrm{SFS}}$ produce the same $F$-score, a higher sensitivity is given more preference. This is because it is critical to ensure that the malignant lesions are correctly classified to provide adequate and timely treatment to the patients.

Fig. 5

Flow diagram showing (a) feature selection process using LOPO-CV along with SFS algorithm, and (b) detailed steps involved in the SFS algorithm. The number of features selected, $n_{\mathrm{SFS}}$, is varied from 1 to 7 for all feature pools and 1 to 6 for intensity feature pool. SFS, sequential forward search; LOPO-CV, leave-one-patient-out cross-validation; AUC,– area under the curve.

Since the LOPO-CV iterates $N_{\text{patient}}$ times, where $N_{\text{patient}}$ is the number of patient data, the features $f_{\mathrm{SFS}}$ selected in each iteration for a particular value of $n_{\mathrm{SFS}}$ depend on the ($N_{\text{patient}}-1$) patient data that are not left out by that iteration. Thus, there can be some variation in the features that are picked out in each iteration of the LOPO-CV. Therefore, the selection frequency is noted, which is defined as the number of times each feature becomes part of the $f_{\mathrm{SFS}}$ during all iterations of the LOPO-CV. This allowed to identify the most frequent (thus most relevant) features ($f_{\text{selected}}$) from each feature pool. It is to be noted that the number of features in each $f_{\text{selected}}$ is denoted as $n_{\text{selected}}$.

2.5.

Classification of Skin Cancer Lesions Using Selected Features

QDA classifiers trained on the best features, $f_{\text{selected}}$, from the three feature pools as explained in Sec. 2.4 were also combined in an ensemble fashion as shown in Fig. 6(a). Separate classifiers are trained on “$k$” feature pools with an LOPO-CV loop. If all the feature pools are used, $k=3$, otherwise, $k=2$. The left-out patient data from the LOPO-CV are tested on each of the individual classifiers, generating a set of posterior probabilities, $P_{\text{pool}\text{-}1},P_{\text{pool}\text{-}2},\dots ,P_{\text{pool}\text{-}k}$, corresponding to classifiers trained on each feature pool. Here, Pool-k is either phasor, intensity, or biexponential feature pool. Subsequently, a weighted average of the posterior probabilities is computed as

Eq. (18)

$$P=w_1{P}_{\text{pool}-1}+w_2{P}_{\text{pool}-2}\dots +w_k{P}_{\text{pool}-k},$$

Eq. (19)

$$w_1+w_2+\dots +w_k=1,$$

where $w_1,w_2,\dots ,w_k$ are the weights on the posterior probabilities generated from each feature pool, while the sum of weights equals one. A threshold of 0.5 on the weighted average of the posterior probabilities assigns a label for the left-out patient data. Figure 6(b) shows the process of weight optimization for the ensemble classifiers. Weight optimization is performed within another LOPO-CV loop. The classifiers are trained on the data from $(n_{\text{patients}}-2)$ patients and the left-out patient data are predicted to generate corresponding posterior probabilities. These posterior probabilities are combined using a weighted average. The weights are varied from 0 to 1 in steps of 0.1. The corresponding sensitivities and specificities are obtained using a threshold of 0.5 on the weighted-average probability. An ROC curve is constructed on these weights, and the weight closest to the ideal point ([0,1]) is selected.

Fig. 6

(a) Schematic of classification of skin lesions. The posterior probabilities from the individual classifiers are combined in an ensemble fashion. (b) Weight optimization for the ensemble classifier. The optimum weight is selected from the ROC curve. LOPO-CV, leave-one-patient-out cross-validation; QDA, quadratic discriminant analysis; ROC, receiver operator characteristics.

3.

Results

3.1.

maFLIM Dermoscopy Clinical Imaging of Skin Lesions

The distribution of patients ($n_{\text{patients}}=30$) and lesions ($n_{\text{lesions}}=60$) imaged in this study showing benign and malignant conditions is provided in Table 3. Benign lesions included solar lentigo and pSK, while malignant lesions included pigmented superficial BCC, pigmented nodular BCC, and melanoma.

Table 3

Distribution of imaged benign and malignant lesions.

	Type	No. patients	No. lesions
Benign	Solar lentigo	2	10
	Pigmented seborrheic keratosis	15	31
Malignant	Pigmented superficial BCC	2	6
	Pigmented nodular BCC	5	5
	Melanoma	6	8

Figure 7(a) shows a handheld maFLIM dermoscope imaging the forearm of a patient. The clinical photograph of a sample melanoma skin lesion is shown in Fig. 7(b), and its corresponding maFLIM feature maps are shown in Fig. 7(c). The scales of the feature maps across the three spectral wavelengths are kept the same for comparison purposes. Most of the feature maps (including ${\alpha }_{\text{fast},390}$, ${\alpha }_{\text{fast},452}$, ${\tau }_{\text{fast},452}$, ${\tau }_{\text{fast},500}$, ${\tau }_{\text{slow},390}$, ${\tau }_{\text{slow},452}$, ${\tau }_{\text{slow},500}$, ${\tau }_{\mathrm{avg},390}$, ${\tau }_{\mathrm{avg},452}$ ${\tau }_{\mathrm{avg},500}$, $I_{n,452}$, $I_{n,500}$, $I_{452,n}/I_{500,n}$ and $I_{390,n}/I_{500,n}$) clearly show two distinguishable regions: the center and the surrounding regions. The cluster masks for this sample image are also shown in the last row of Fig. 7(c). The maps showing the two clustered regions: region-1 and region-2 are plotted using the normalized integrated intensities ($I_{\text{Integrated}}$) from the three spectral emission channels.

Fig. 7

(a) Handheld maFLIM dermoscope imaging the forearm of a patient. (b) Clinical photograph of a melanoma lesion. (c) Time-domain maFLIM feature maps of a melanoma lesion. The columns show the feature maps corresponding to the three emission channels. First row shows the weight of the fast decay. Second row shows the fast lifetime maps, while the third row shows the slow lifetime maps. Average lifetime maps are shown in the fourth row. Fifth row shows the integrated intensity maps of each spectral emission channel, and the ratio of the intensities are shown in the sixth row. The last row shows the cluster mask generated for the lesion and the integrated intensities from all the channels for the clustered regions 1 and 2. The horizontal strip in the images is due to the presence of hair on the skin during imaging.

The relative contributions of fast lifetime (${\alpha }_{\text{fast}}$) for the spectral channels 390 and 452 nm exhibit a higher value in the central region compared with the surrounding regions. Fast lifetime (${\tau }_{\text{fast}}$) values for the spectral channels 452 and 500 nm are lower in the central region than the surrounding. Slow lifetime (${\tau }_{\text{slow}}$) values for all the three spectral channels are higher in the center than the surrounding regions. The average lifetimes (${\tau }_{\mathrm{avg}}$) of the pixels from the central region are higher than the surrounding parts. The relative intensity ($I_n$) values for spectral channel 452 nm is lower in the center compared to the surrounding, while that for the spectral channel 500 nm is higher in the center compared with the surrounding. The ratios of the intensities, $I_{452,n}/I_{500,n}$ and $I_{390,n}/I_{500,n}$, are lower in the central region compared with the surrounding regions.

3.4.

Classification of Skin Cancer Lesions Using Selected Features

The methodology of classifying skin lesions using the selected features is explained in Sec. 2.5. As shown in Fig. 6, “$k$” QDA classifiers trained separately on the $f_{\text{selected}}$ features from “$k$” feature pools are combined in an ensemble fashion. Four different combinations of feature pools are used for constructing ensemble classifiers, as shown in Table 6. A weighted average of the posterior probabilities generated from the “$k$” QDA classifiers is calculated to predict the label of the left-out patient data. Since the weights are optimized within an LOPO-CV loop, there can be some variation in the weights selected during the $n_{\text{patients}}$ iterations.

Table 6

Performance metrics of ensemble classifiers trained with multiple combinations of feature pools.

Feature sets for ensemble	Accuracy (%)	Sn (%)	Sp (%)	F-score (%)	Confusion matrices
					True	Predicted
						Benign	Malignant
${(f_{\text{selected}})}_{\text{biexponential}}+{(f_{\text{selected}})}_{\text{intensity}}$	71.67	84.21	65.85	65.31	Benign	27	14
					Malignant	3	16
${(f_{\text{selected}})}_{\text{phasor}}+{(f_{\text{selected}})}_{\text{biexponential}}$	88.33	84.21	90.24	82.05	Benign	37	4
					Malignant	3	16
${(f_{\text{selected}})}_{\text{phasor}}+{(f_{\text{selected}})}_{\text{intensity}}$	86.67	78.95	90.24	78.95	Benign	37	4
					Malignant	4	15
${(f_{\text{selected}})}_{\text{phasor}}+{(f_{\text{selected}})}_{\text{biexponential}}+{(f_{\text{selected}})}_{\text{intensity}}$	88.33	84.21	90.24	82.05	Benign	37	4
					Malignant	3	16

Figure 8 shows the histograms of the weights obtained during all the iterations for the ensemble classifiers using the four feature pool combinations. Since the sum of weights of “$k$” feature pools is one, it is sufficient to show the weights of ($k-1$) feature pools. It can be seen from Table 6 that the ensemble combination of QDA classifiers trained on phasor and biexponential features as well as the ensemble combination of all the three feature pools, produced the best performance with 88.33% accuracy, 84.21% sensitivity, 90.24% specificity, and $F$-score of 82.05%. The next best performance is obtained by the ensemble combination of phasor and intensity feature pools with an accuracy of 86.67%, sensitivity of 78.95, specificity of 90.4%, and an $F$-score of 78.95%. Ensemble combination of intensity and biexponential features produced an accuracy of 71.67%, sensitivity of 84.21%, specificity of 65.85%, and $F$-score of 65.31%. Table 6 also shows the confusion matrices for all the classifiers. While analyzing the weights on the feature pools during the ensemble combination, it can be seen in Figs. 8(a) and 8(b) that intensity features have a higher weightage when combined with phasor features. When phasor and biexponential features are combined, both the feature pools have similar weightage as shown in Fig. 8(c). Similarly, when biexponential and intensity feature pools are combined, both the pools have comparable weightage. When all the feature pools are combined, it can be seen from Fig. 8(d) that the phasor and biexponential features pools have similar weight distribution in the range [0.1, 0.5]. The weights on the intensity feature pool are widely dispersed in the range [0, 0.8].

Fig. 8

Histogram of weights on one of the feature pools, when combined in an ensemble fashion for (a) phasor-intensity, (b) biexponential-intensity, (c) phasor-biexponential, and (d) phasor-biexponential-intensity feature pools.

4.

Discussion

In this study, clinical widefield autofluorescence imaging of benign and malignant pigmented skin lesions was successfully performed in 30 patients using a recently developed maFLIM dermoscope.⁴⁰ The resulting maFLIM images from 60 pigmented lesions enabled exploring steady-state (intensity) and time-resolved (biexponential, phasor) autofluorescence global features. Results based on rigorous cross-validation methods demonstrate that simple ML classification models (QDA) based on selected time-resolved autofluorescence global features have the potential to provide discrimination of malignant from benign pigmented skin lesions.

To the best of our knowledge, only one published work has reported the use of machine learning models based on autofluorescence lifetime imaging features for the classification of pigmented skin lesions.¹² In that study, however, only skin melanoma lesions were imaged, and the classification task was restricted to discriminate early-stage from advanced-stage skin melanoma. In contrast, a more comprehensive set of pigmented skin lesions were imaged in this work (two benign and three malignant lesion categories). Moreover, the classification task focused on discriminating malignant from benign pigmented lesions, which might be clinically more relevant for early detection of skin cancer.

In multidimensional imaging data, such as in maFLIM, image features can be extracted at the pixel or the image level. We have recently explored pixel-level maFLIM features for the classification of oral dysplasia and early-stage cancer.⁴⁹ Pixel-level features, however, require the labeling of each pixel which is generally impractical. In this work, the maFLIM data were labeled at the lesion level based on the histopathology diagnosis obtained from the lesion biopsy samples; therefore, an image-level global feature extraction strategy was preferred. As shown in Fig. 7, two regions were frequently observed in the maFLIM images, corresponding to pixels either within or outside the lesion extension. In an attempt to reduce interpatient variability, the explored relative features were defined in terms of difference in autofluorescence properties between the two regions identified in each lesion maFLIM image. This strategy of using global (image-level) and relative features can find applications in many other classification tasks based on optical imaging data.

The performance of ML classification models needs to be carefully estimated when trained on limited sample size. To minimize overfitting and avoid overoptimistic performance estimations, a rigorous strategy was adopted for feature selection, model training, and performance estimation. First, the maximum number of features allowed (seven) was limited based on the sample size.⁴⁸ Second, a simple nonlinear classification model (QDA) was adopted. Third, cross-validation was applied at the patient-level (LOPO-CV) to ensure that data from the same patient is not used for both training and validation. Fourth, feature selection was performed together with model training to make sure that the validation data are not used during neither feature selection nor model training. It should be noted that at each fold of the cross-validation strategy, a different classification model (with different selected features and model parameters) is applied to the validation set. Thus, although a single optimal model is not necessarily defined, this approach still enables identifying relevant features and providing unbiased classification performance estimation.

Classification performance was dependent of the feature pool used in the model. The most frequently selected intensity feature was $I_{452,n}$ which is associated to NADH fluorescence contribution. Although the classification models using intensity features showed good sensitivity ($\sim 84\%$), their specificity was poor ($\sim 31\%$). This means that steady state intensity features are not sufficient to minimize false positives while discriminating benign and malignant lesions. Classification models using biexponential features showed similar sensitivity ($\sim 84\%$) than those with intensity features, but significantly higher specificity ($\sim 70\%$). These results indicate that time-resolved properties of pigmented skin lesion autofluorescence could represent biomarkers of skin cancer. While examining the selected features from the biexponential feature pool, it can be seen from Table 5 that the most frequently selected features were associated to NADH (${\tau }_{\text{fast},452}$, ${\alpha }_{\text{fast},452}$) and collagen (${\tau }_{\text{slow},390}$, ${\alpha }_{\text{fast},390}$) fluorescence temporal dynamics. The fast lifetime component of NADH (${\tau }_{\text{fast},452}$, ${\alpha }_{\text{fast},452}$) is associated with the free state of the molecule,¹⁹ suggesting that metabolic pathway changes induced by malignant transformations alter the microenvironment of NADH molecules. Unlike NADH or FAD, the collagen is not involved in cellular respiration and is not part of metabolic pathways. Changes in collagen autofluorescence response can occur in benign and malignant lesions as they mostly correspond to skin thickening, extracellular matrix remodeling, and texture changes.²³ Therefore, the selected features corresponding to collagen autofluorescence (${\tau }_{\text{slow},390}$, ${\alpha }_{\text{fast},390}$) can provide high sensitivity but may not contribute to high specificity leading to several false positives as can be seen in Table 4.

The phasor representation of fluorescence lifetime imaging data analysis is a noniterative, model-free, fast approach to visualize the lifetime components (and their distributions) of the fluorescence emission of a sample.^50–53 In this work, a modified version of this method was applied to the clinical multispectral maFLIM data, and a set of global image features were extracted from the corresponding phasor representation of the two regions present in each lesion maFLIM image (Fig. 4). The processing times taken for computing phasor and biexponential features from one sample lesion using a computer with an i7 Intel core processor and 48 GB RAM were found to be 10.5 and 251.8 s, respectively. Therefore, phasor-based features are simpler and $\sim 25$ times faster than traditional time resolved biexponential features. Unlike biexponential maFLIM features, the phasor features explored cannot be directly interpreted in terms of the skin autofluorescent constituents. Nevertheless, the classification models using phasor features showed superior specificity ($\sim 80\%$) than those using biexponential features, suggesting that these two different pools of time-resolved maFLIM features might be complementary. Phasor-based features are computed from the concatenated signal and contain information from all three emission channels. It is challenging to develop solid reasoning for the high specificity obtained from these features; however, since every phasor-based feature contains information from all three channels, each channel’s contribution is not quantifiable.

Feature selection starting with combinations of feature pools were also explored, although the resulting models showed lower classification performance overall (Table 4). On the other hand, ensemble classifiers based on models using the most frequently selected features of each pool outperformed any other models (Table 5). In particular, the ensemble classifier combining the models based on the most frequently selected biexponential and phasor features resulted in the best performance overall ($\sim 84\%$ sensitivity and $\sim 90\%$ specificity). Moreover, the optimum weights identified for these ensemble models indicate that both the phasor and biexponential features contribute equally to the weighted probability [Fig. 8(c)]. These results further indicate that these two maFLIM feature pools might be complementary, as biexponential features seem to contribute to higher sensitivity, while phasor features to higher specificity.

Some of the studies summarized in Table 1 produced superior classification performances compared with that reported in this work. However, those studies employed complex machine learning and deep learning algorithms, such as Adaboost,¹⁰ BoF,¹⁰ and CNN’s for classification,¹⁸ and feature extraction.¹⁷ This work uses a simple QDA machine learning model to analyze the performance with different feature pools. It is plausible that more complex machine learning classification models, such as neural networks, might provide superior performance. Given the limited data available for training and validation, however, a simple QDA model was selected, as it still enables nonlinear decision boundaries while reducing the chances of overfitting. In addition, by choosing a simple classification method, it is possible to analyze other aspects of classification model, such as the feature pools, number of features, and ensemble combinations of the feature pools.

Recent advances in artificial intelligence (AI) are allowing the development of CAD systems for discriminating benign from malignant skin lesions based on digital dermoscopy data. Although these CAD systems have not been translated yet to the clinic, preliminary validation studies demonstrate their potential to discriminate typical benign from malignant lesions. Dermoscopy data, however, provide information limited to the visual appearance of the lesions; thus, AI-assisted digital dermoscopy is less suitable for the discrimination of visually similar benign and malignant skin lesions. maFLIM dermoscopy data, on the other hand, can capture autofluorescence-based biochemical and metabolic biomarkers of skin malignant transformation, which is independent of lesion visual appearance. This property is especially important in distinguishing visually similar benign and malignant lesions, thus minimizing the number of false positives and, in turn, reducing the number of unnecessary biopsies. In particular, the lesion types included in this work—pSK and melanoma, are visually similar benign and malignant lesions, which can be easily misdiagnosed during dermoscopic evaluation. Therefore, maFLIM dermoscopy has the potential to complement digital dermoscopy by providing superior performance in the discrimination of visually similar benign and malignant skin lesions.

5.

Conclusion

The results of this study demonstrate the capabilities of maFLIM dermoscopy to clinically image a plurality of autofluorescence biomarkers of malignant skin pigmented lesions. Moreover, some of these autofluorescence biomarkers were identified as promising for malignant lesion identification, particularly those quantifying the time-resolved fluorescence characteristics of skin lesions. In addition, these relevant autofluorescence biomarkers were successfully used as features in ML models trained to discriminate malignant from benign pigmented skin lesions with promising accuracy ($\sim 84\%$ sensitivity and $\sim 90\%$ specificity). Further developments in maFLIM instrumentation and image analysis methods could result in clinical tools for noninvasive, label-free, accurate, and objective in situ detection of malignant from benign skin lesions, with the potential to impact the clinical management of skin cancer patients.