Unsupervised hyperspectral data mining and bioimaging by information entropy and self-modeling curve resolution

Unsupervised estimation of the dimensionality of hyperspectral microspectroscopy datasets containing pure and mixed spectral features, and extraction of their representative endmember spectra, remains a challenge in biochemical data mining. We report a new versatile algorithm building on semi-nonnegativity constrained self-modeling curve resolution and information entropy, to estimate the quantity of separable biochemical species from hyperspectral microspectroscopy, and extraction of their representative spectra. The algorithm is benchmarked with established methods from satellite remote sensing, spectral unmixing, and clustering. To demonstrate the widespread applicability of the developed algorithm, we collected hyperspectral datasets using spontaneous Raman, Coherent Anti-stokes Raman Scattering and Fourier Transform IR, of seven reference compounds, an oil-in-water emulsion, and tissue-engineered extracellular matrices on poly-L-lactic acid and porcine jejunum-derived small intestine submucosa scaffolds seeded with bovine chondrocytes. We show the potential of the developed algorithm by consolidating hyperspectral molecular information with sample microstructure, pertinent to fields ranging from gastrophysics to regenerative medicine.

Herein we developed an unsupervised algorithm for (i) determining the number distinct and spectrally separable biochemistries contributing to the observed spectral variation in hyperspectral datasets, which we have termed Chemical Dimensionality (CD-DE), and (ii) extraction of endmember spectra according to the estimated dimensionality, also available as a stand-alone extraction algorithm (CD-EEA). We demonstrate widespread applicability of our method, and notable performance as a biochemical data mining and image reconstruction technique, by (i) benchmarking against pseudo-randomly generated synthetic datasets with known dimensionality, (ii) comparison of estimated Chemical Dimensionality (k CD ) to established Virtual Dimensionality estimation methods developed for remote sensing applications, and (iii) benchmarking performance influenced by varying signal-to-noise ratios and sample size. We further demonstrate versatility by (i) imaging samples ranging from liquid immiscibility in oil-in-water emulsions to extracellular matrix regeneration on biomaterial scaffolds for tissue engineering in vitro and (ii) using a plurality of imaging techniques including Coherent Anti-stokes Raman Scattering (CARS) and Fourier Transform Infrared microspectroscopy (FTIR). Last, we compare the hyperspectral endmembers and their associated image reconstructions extracted by CD-EEA against K-means clusters and endmembers extracted by N-FINDR 19 . Given that input to the Chemical Dimensionality framework can range from hyperspectral datacubes to sample feature vectors, we see the algorithm as a promising candidate for any type of multidimensional systems such as chromatography or gene expression data.

Description of the Chemical Dimensionality Algorithm Chemical Dimensionality algorithm
The algorithm is divided into two parts, a model selection algorithm that determines the Chemical Dimensionality (CD) of the input dataset, and an associated endmember extraction algorithm which uses the CD framework to extract prospect endmembers using a least square-error/maximum entropy approach. Here we first treat the model selection algorithm that determines the chemical dimensionality of a hyperspectral dataset. A flowchart is shown in Supplementary Figure S1. Figure 1: (a) Left: common method of data mining in hyperspectral microspectroscopy datasets. Approximate number of endmembers is determined, typically using thresholding methods. Endmember extraction is often repeated in an exploratory fashion, making it prone to over-interpretation of returned endmembers or overlooking distinct biochemistries. Right: Algorithmic overview of the Chemical Dimensionality framework. Detailed algorithm flowchart is shown in Supplementary Figures S1 and S2, along with a full mathematical description in the methods section. Estimation of dimensionality and extraction of endmembers is carried out using information entropy and sum of squared residuals. Returned endmembers can be used to generate image reconstructions for the separated biochemistries. (b) Schematic of the developed dimensionality estimation algorithm. A hyperspectral microspectroscopy image is subset to a candidate endmember matrix by repeated hyperdimensional simplex maximization. A number of blind signal separation models are built assuming an increasing number of source signals, with subsequent estimation of the Chemical Dimensionality based on information entropy and sum of square residuals of model reconstructions.
The CD algorithm considers both datamatrices from spectroscopy experiments, and datacubes from microspectroscopy imaging experiments. Standard datamatrices has dimensions (n n , n ν ), where subscript n denotes samples rather than pixels with an explicit spatial connotation, and ν denotes an axis of wavenumbers. As the algorithm works on datamatrices, no further preprocessing is needed prior initialization, and the dataset is directly assigned as an input datamatrix Z. Hyperspectral datacubes (HSDC) has dimensions (n x , n y , n ν ), for which each pixel (x, y) in the scene contains the associated spectrum p(x, y) of length n ν . We use the term unfolding, confer earlier studies 17 , to describe the conversion from 3D hyperspectral datacubes to a datamatrix ready for computation, by reshaping the hyperspectral datacubes such that the new input datamatrix Z has dimensions (n x n y , n ν ). For the purpose of image reconstruction, it is important to note the original location of each pixel from the hyperspectral datacube in the new datamatrix Z. In the following, we consider the input datamatrix Z (n×p) , from an unfolded HSDC if necessary, containing n samples/pixels each with p wavenumbers/channels. For optimum performance the input dataset Z (n×p) should 1) contain at least one instance of a spectrally "pure" endmember, a prerequisite of most major endmember extraction algorithms 31,19,32 , 2) contain two or more endmembers, and 3) not contain endmembers that are linearly dependent, i.e. no endmembers should be expressible as a linear combination of any other set of endmembers.

Simplex maximization
While hyperspectral image maps from spontaneous Raman and FTIR microspectroscopy are typically relatively small, in the order of hundreds to thousands of pixels in the image, other techniques, such as Coherent Anti-stokes Raman Scattering frequently yield larger images, in the order of hundreds of thousands pixels in the image. We therefore start by subsetting the input datamatrix Z (n×p) to a reduced datamatrix V (m×p) containing m candidate endmembers where m n, subject to m k where k is the expected/ground true number of endmembers in the dataset. The N-FINDR 19 algorithm operates on the principle that k potential endmembers maximizes the volume of a simplex with vertices given from a dimensionality reduced hyperspectral dataset. However, given that the initial simplex is formed from a random set of samples/in the data, there is no guarantee that the algorithm consistently returns the same set of endmembers 33 . Furthermore, even if the dataset is known to contain k endmembers, it may require extracting more than k endmembers to return ground true endmembers. Nonetheless, returned endmembers from a number of repeated simplex maximizations can be considered candidate endmembers, and therefore used to generate the reduced candidate endmember matrix V (m×p) . In principle most endmember extraction techniques can be used to generate the endmember candidate matrix V, or if the dataset is small, the input datamatrix Z (n×p) can be used directly as a candidate endmember matrix. However, here we provide a full implementation, and adapt the simplex maximization technique described by Winter 19 .
Let i be the number of endmembers extracted using simplex maximization. The dimensionality of Z (n×p) is reduced to i˘1 components using Singular Value Decomposition. Using the dimensionality reduced matrix, the endmember matrix, E i is initialized with i random selection of pixels/samples.
Each pixel in the dataset is then substituted in the place of each prospect endmember e i , and retained if it further maximizes the simplex volume given by: This process is repeated for i = 2, 3, . . . , g, where the resulting endmembers in E i found from maximization of V olume(E i ) is continuously placed in the endmember candidate matrix V (m×p) . Owing to later comparison between the endmember candidates, it is required that g − 1 > k. We found that repeating the simplex maximization up to extraction of g = 20 endmembers generally is suitable for microspectroscopy, yielding a candidate endmember . For large input data matrices Z (n×p) , generation of the candidate endmember matrix by repeated simplex maximization can advantageously be implemented using parallel computing.

Self-modelling Curve Resolution using semi-nonnegativity constrained matrix factorization
Having repeated simplex maximization to generate the candidate endmember matrix V (m×p) , the k ground true endmembers are to be found among the m candidates in the matrix. However, with no preconception as to the magnitude of k and the corresponding endmember spectra, we can treat it as a blind signal separation problem. A canonical assumption is that endmembers E = [ e 1 e 2 . . . e k ], are considered vectors in a real vector space and describe mixture voxels, pixels or samples through convex mixing. This implies that, for each voxel, pixel or sample with spectrum p (1×p) in the unfolded data input matrix: Where x can be thought of as a vector of mixing weights for the convex mixing: While the sum-to-one constraint on linear mixing models is reasonable, it is an unnecessary constraint when determining the Chemical Dimensionality. Furthermore, not all microspectroscopy techniques lend themselves equally well to quantitative, idealized mixing models, for instance due to an inherent non-resonant background such as in conventional Coherent Anti-stokes Raman Scattering (CARS) imaging. If the sum-to-one constraint is dropped, estimating p i (1×p) becomes a conical mixing problem, retaining the non-negativity constraint on x i ≥ 0. The mixing weights in x can then be obtained by solving a non-negative constrained minimization problem in a least-squares sense: Where · 2 denotes the 2 norm. To determine the Chemical Dimensionality we let E be of size (p × u), i.e. u feature vectors with data collected at p wavenumbers/channels. As the signature of the u feature vectors is unknown, determining E becomes a blind signal separation problem. Here, blind signal separation was carried out using nonnegative matrix factorization (NMF). At its core, non-negative matrix factorization seeks to express a non-negative matrix X ∈ R p×n in terms of W ∈ R p×k and H ∈ R k×n , such that: X ≈ WH, s.t. W, H ≥ 0 (6) In Eq. 6, X can be thought of as the input spectral data matrix, W a set of u feature vectors i.e. "quasi-endmembers", and H containing the weights of each quasi-endmember. Depending on the type of spectroscopy used for data acquisition, candidate endmembers may contain slightly negative arbitrary intensities, for instance owing to prior preprocessing with baseline correction algorithms, challenging the non-negativity constraint on X. We therefore adopted the more lenient semi-NMF method by Ding et al. 34 , in which W is unconstrained and H is strictly non-negative. The semi-NMF optimization problem is thus described by: Where · F denotes the Frobenius norm. In this Chemical Dimensionality algorithm we used the semi-NMF method as implemented by Li and Ngom 35 , using the fast combinatorial non-negative least squares algorithm developed by van Benthem and Keenan 36 . From Eq. 7 we generate g models consisting of u feature vectors where u = 1, 2, . . . , g, as for the repeated simplex maximization, and X is the candidate endmember matrix V ( g n=2 n) × p , such that:

Error reduction
To select the final model of W ∈ R p×u and H ∈ R u×n , where u becomes the Chemical Dimensionality, we adapted the Shannon Entropy 20 , a concept from information theory, to calculate a total entropy of the reconstruction residuals. We also calculated an error reduction term for each of the u blind signal separation models, and together, the two metrics were used to determine the Chemical Dimensionality. Firstly, for both metrics, let R u contain the reconstruction residuals from the semi-nonnegative matrix factorization with u feature vectors: We defined the error reduction for model u as the ratio between the reduction in normalized sum of square error for models u − 1 to u, and u to u + 1 ∀u = 2, 3 , . . . , g − 1. Alternatively, the error reduction in this sense can be thought of as measure of "sudden" improvements in how well a semi-NMF model with u feature vectors describes the candidate endmember matrix V (m × p) . Denoting the matrix elements of R u as R u,mp , the sum of square error for the u th semi-NMF model, s u , is given in Eq. 10, which for interpretability is normalized, such that ε u ∈ [0, 1]. Noting that s u ∈ R + ∪ {0} ∀u ∈ 1, 2, . . . , g should be monotonically decreasing as u goes from 1 to g: u = 1 max 1≤i≤g s i s u ∀u ∈ {1, 2, . . . , g} As u varies from 1 to g in steps of 1, we calculate the error reduction, ρ ∈ R + ∪ {0}, for model u by:

Entropy of reconstruction
Subsequently, we calculated the information theoretic metric based on the Shannon entropy 20 , which hereinafter will be referred to as the reconstruction entropy. The reconstruction entropy is given by: Inspired by prior use of Shannon entropy 37 , we calculated b i such that b i ∀i sums to unity. Let r c,u = [r c,u,1 r c,u,2 · · · r c,u,p ] T , where c ∈ {1, 2, . . . , m}, be a residual vector of R u with dimensions (p × 1) across h∈{2, 3, ..., p} When the elements of b( r c,u , v) i are numerically equivalent, H is maximized. Correspondingly, the reconstruction entropy will tend to increase as the residuals tend to identically distributed noise, and decrease if the residuals carry spectral information. For a given semi-NMF model with u feature vectors, we calculate a total entropy score S u ∈ R + by summing the entropy from Eq. 13 for all c ∈ 1, 2, . . . , m residual vectors in R u , which is then repeated for ∀u ∈ {1, 2, . . . , g}.

Estimation of the Chemical Dimensionality
The reconstruction entropy is expected to increase, as the reconstruction residuals tend to identically distributed random noise. Consequently, maximizing the reconstruction entropy may not (necessarily) estimate a meaningful dimensionality. Instead, the chemical dimensionality, k CD , is estimated by searching for an increase in reconstruction entropy at or in the vicinity of the u th semi-NMF model, providing the maximum error reduction.
If model z was also associated with an increase in reconstruction entropy, S z − S (z−1) > 0 then k CD = z. Otherwise, we seek the nearest model z ± i, i ∈ {1, 2, 3, ...}, subject to (z ± i) ∈ N | 2 ≤ z ± i ≤ g − 1, associated with an increase in reconstruction entropy is the nearest model associated with an increase in reconstruction entropy, then k CD = z + i, and conversely In case both S (z−i) − S (z−i−1) > 0 and S (z+i) − S (z+i−1) > 0 for a given i, then we seek the model that resulted in the greatest increase in reconstruction entropy, i.e.

Chemical Dimensionality endmember extraction
After estimation of the Chemical Dimensionality, k CD , the algorithm proceeds to extract k CD samples/pixels from the endmember candidate matrix identified as actual endmembers of the dataset, according to the 2 norm of the reconstruction residuals and the reconstruction entropy. First, we define the Chemical Dimensionality endmember matrix E CD , containing k CD endmembers, such that: (17) It should be noted that, as a consequence of the repeated simplex maximization in Eq. 3, the endmember candidate matrix may contain multiple occurrences of the same endmember candidate. For the initialization in Eq. 17, E CD is consequently populated with the k CD samples/pixels in the endmember candidate matrix V m × p with the lowest 2 norm, subject to e 1 = e 2 = . . . = e kCD . This condition must also be checked when the algorithm proceeds to iteratively swap endmembers in E CD with new prospect endmembers.
Having populated the initial endmember matrix E CD , the non-negativity constrained least squares minimization problem in Eq. 5 is solved using E = E T CD and p i as the i th candidate endmember in V (m × p) , ∀i ∈ 1, 2, . . . , m.
Having solved for x, the performance of endmembers initially present in E CD in terms of the 2 norm is then given by: To calculate the reconstruction entropy of the endmembers, let: be a residual vector from the minimization problem in Eq. 18, where c ∈ {1, 2, . . . , m} and r c has dimensions (p × 1) across the associated wavenumbers v = [ν 1 ν 2 · · · ν p ] T , then b( r c , v) i is calculated using Eq. 14, and the reconstruction entropy becomes: The algorithm loops through each of the m candidate endmembers in V, which is then substituted in place of each e j , j ∈ 1, 2, . . . , k CD . For each swap, Eq. 18-21 is calculated. If swapping one or more endmembers in E CD with a given candidate endmember from V results in a reconstruction 2 norm lower than P 2 , then the candidate endmember is swapped in place of the endmember in E CD for which the reconstruction entropy, P S , was maximized compared to the other swaps, regardless of whether the swap also resulted in the minimum reconstruction 2 norm or not. Subsequently, P 2 is updated with the reconstruction 2 norm of the new endmember matrix E CD . The procedure is continued until looping through the candidate endmember matrix V no longer results in any swaps. A flowchart of the endmember extraction algorithm can be found in Supplementary Figure S2.

Generated datasets
To benchmark the Chemical Dimensionality algorithm we generated a series of datasets with 1) a known number of endmembers ranging from 2 to 10 endmembers, each with 2) random convex mixtures of the 9 endmember sets, to a total of 300, 1000 and 5000 samples, and 3) random convex mixtures of endmember sets to a total of 5000 samples with signal-to-noise ratios (SNR) ranging from 1000:1 to 10:1. New endmembers were generated for each endmember sets k = 2, 3, . . . , 10. For an endmember set with k endmembers, each of the k pseudo-spectra were generated by mixing a random number f = 1, 2, . . . , 6, of Gaussian functions each with random height (a), width (s) and peak center on the spectral axis (m) (Eq. 22). To test the algorithm specifically on its performance on datasets of various sample sizes and signal-to-noise ratios, all endmembers were vector normalized to unit length prior mixing, to avoid confounding with scale. Although arbitrary, we let the x-axis range from ν 1 = 900 to ν 1001 = 1900, i.e. consisting of 1001 channels/wavenumbers. Where: Once the 9 endmember sets were prepared, the endmembers of each set were mixed such that k j=1 w j = 1, where weights w j were pseudo-randomly generated to mix a total of 300 − k, 1000 − k and 5000 − k spectra. The k original Table 1: Overview of synthesized datasets with varying endmember constituents, sample mixtures and signal-to-noise ratios.

Samples in mixture dataset
Signal-to-noise ratio endmembers were then appended to the data matrix, and the entire data matrix was shuffled. The final locations of the k original endmembers in the data matrix were therefore unknown.
Prior running the Chemical Dimensionality algorithm on the synthesized datasets, white Gaussian noise was added to each sample in each mixture dataset, to avoid strict linear dependency. For studying the effect of sample size on the CD algorithm, all mixture datasets for each endmember set, was added noise to an SNR of 1000:1, while the mixture datasets containing 5000 samples were also added white Gaussian noise to an SNR of 100:1, 50:1, 20:1 and 10:1, to test the performance of the algorithm with increasing noise levels ( Table 1).
In this work we calculated signal/spectral strength as the root-mean-square of the signal. Consider a pixel/sample p i with the associated spectrum Signal strength E RMS is then given by: The signal strength E RMS was calculated for each of the samples in each of the mixture datasets. The weakest signal was used to determine the amount of white Gaussian noise to be added to obtain the desired signal-to-noise ratio. The same amount of noise was added to all remaining samples in the dataset, such that the noise level in each sample at least matches the stated SNR.

Benchmarking Chemical Dimensionality to established Virtual Dimensionality methods
The estimated Chemical Dimensionality, k CD , was compared to the dimensionality estimated by five different methods from analytical chemistry, information theory, and Neyman-Pearson detection theory, all of which has previously been applied to hyperspectral imagery from satellite remote sensing 38,39 in estimation of the Virtual Dimensionality.
Benchmarking was based on well-defined synthesized datasets and was carried out for cases where 1) p < n, p ≈ n, and p > n, where p is the number of channels, and n is the number of pixels/samples in the dataset, and 2) signal-tonoise ratios as follows; 1000:1, 100:1, 50:1, 20:1 and 10:1.
On the basis of Neyman-Pearson detection theory, the dimensionality of the generated datasets was estimated by the Harsanyi-Farrand-Chang method(Chein-I Chang, 2013) (HFC) using a false-alarm probability of 10 −5 , and from information theory the dimensionality was estimated by the Akaike Criterion (AIC) and the Minimum Description Length (MDL), as detailed in prior works 22,39 . These methods ideally require that n p, and especially for the AIC and MDL criteria must noise be independent and identically distributed 39 . The Chemical Dimensionality was also compared to the dimensionality predicted by Malinowski's Factor Indicator Function (FIF), which was developed and specifically tested on data from spectroscopy, magnetic resonance and chromatography 23 . Last, we implemented an explained variance thresholding method which frequently appears in literature. Here, the dimensionality was estimated by finding the number of principle components for which the explained variance would describe at least 95 % and 99 % of the total variability, here denoted λ 95 and λ 99 respectively, by Principal Component Analysis (PCA) using Singular Value Decomposition (SVD). As these methods cannot extract the estimated number of endmembers from a dataset, the Chemical Dimensionality algorithm was benchmarked separately on its ability to extract endmembers from a hyperspectral dataset.

Benchmarking Chemical Dimensionality endmember extraction
Endmembers extracted by the Chemical Dimensionality algorithm were benchmarked to those extracted by N-FINDR 19 , a simplex maximization-based algorithm that has previously been extensively deployed for unmixing of hyperspectral data from various microspectroscopy techniques 15,14,16,17 . Extraction was conducted for both the synthesized datasets as well as experimental datasets, and the number of endmembers extracted was dictated by the estimated Chemical Dimensionality, k CD . For the synthesized datasets and the spontaneous Raman spectroscopy dataset of seven compound standards, the endmembers extracted by N-FINDR and the CD-EEA were compared to the ground true, known, endmembers. For the synthesized dataset, performance was evaluated with 2, 5 and 10 known endmembers mixed to 5000 samples with SNR 1000:1. The experimental hyperspectral imaging datasets enabled image reconstruction by solving the non-negativity constrained least squares (NNLS) problem in Eq. 5. Here, image reconstruction was carried out by NNLS of endmembers extracted by the N-FINDR and CD-EEA again using the estimated k CD , and these were further compared with their respective cluster images obtained from K-means clustering using the K-Means++ algorithm 40 .

Computational setup
The algorithm was implemented in a x64 version of MATLAB R2018b Update 3 (MathWorks Inc., Natick, MA, USA) on a desktop computer running Windows 10. The system consisted of 10 cores Intel® Core™ i9-9820X CPU @ 3.30GHz and 128 GB of RAM. The algorithm was implemented exploiting the benefits of parallel computing where appropriate. We ran the algorithm using j max = 5 for g = 2, 3, . . . , 20 simplex maximization repeats (see Supplementary Figure S1), yielding a candidate endmember matrix of 209 candidate endmembers. All 10 cores were deployed for parallel computing.

Experimental procedures Spontaneous Raman spectroscopy of compound standards
Seven compound standards were measured using spontaneous Raman spectroscopy. The standards were from a Bruker Standard Reference kit, and included; Polystyrene, Tylenol, Silicon, Calcium carbonate and Naphthalene. Stearic acid was obtained from Sigma Aldrich and β-Tricalcium phosphate was a kind donation from Particle3D. Raman spectra were collected using a Fiber Optic Probe-based setup as previously described 41 . Briefly, spontaneous Raman spectra were collected using a near-infrared laser (B&W TEK Inc., Newark, DE, USA) with excitation at 785 nm, and a power output at 50 mW. The Filtered Fiber Optic Probe (EmVision LLC., Loxahatchee, FL, USA) contains seven low OH collection fibers with NA = 0.22 and 300 µm core, arranged in a donut-shape around the 300 µm core (NA = 0.22) excitation fiber. Filtered Raman signals in the collection fibers were led to an Ibsen Eagle Raman-S OEM spectrometer (Ibsen Photonics A/S, Farum, Denmark). The Ibsen Eagle Raman-S spectrometer is based on a transmission grating design with light projected onto an Andor iVac 316 (Andor, Belfast, Northern Ireland) NIRoptimized back-illuminated Low Dark Current Deep-depletion CCD (2000 × 256 pixels) cooled to −70°C using the in-built thermoelectric Peltier cooler. Spectra were collected in full-vertical binning mode with 1s integration and 5 co-additions. The collected Raman spectra for each of the seven standards was subsequently baseline-corrected using a combined asymmetric least squares Whittaker smoothing method 42,43 , with input parameters 10 −6 ≤ p ≤ 10 −3 and 10 2 ≤ λ ≤ 10 4 . Here explicitly p is an asymmetry parameter, and λ is a smoothing parameter 43 .
The compound standards were used to test the performance of the algorithm with experimentally obtained spectra. Therefore, contrary to the synthesized endmember sets, the spontaneous Raman standards were not vector normalized prior mixing. The seven standards were mixed as-is, such that such that k j=1 w j = 1 to a set of 5000 samples, and added noise to an SNR of 1000:1 as previously described to avoid strict linear dependency.
FTIR Imaging of Polylactic acid (PLLA) and Small Intestine Submucosa (SIS) tissue scaffolds PLLA scaffold was prepared by dissolving PLLA in 1,1,1,3,3,3-hexafluoroisopropanol (HFIP) at 8 % w/v overnight. Dissolved PLLA was subsequently electrospun into a randomly oriented scaffold consisting of 1 µm PLLA fibers, to a mean scaffold thickness of 120 µm. PLLA and HFIP was purchased from Sigma Aldrich. Small Intestine Submucosa (SIS) biomaterial scaffolds were derived from porcine jejunum by mechanical delamination and subsequently decellularized for two hours in 0.1 % (vol/vol) peracetic acid, 4 % (vol/vol) ethanol, and 96 % (vol/vol) DI water. Prior final preparation of the multi-layered scaffolds, electrospun PLLA was sterilized in 70 % ethanol and extensively rinsed in PBS. SIS Scaffolds were freeze-dried and sterilized with ethylene oxide, and rehydrated in sterile PBS prior layering. Five layers of each scaffold were used for layering, and these were each seeded with isolated bovine chondrocytes with 200K cells per layer. After seeding and cell adhesion, the media was replaced by a chondrogenic differentiatiation medium containing DMEM (4.5 g l −1 l-glucose), 50 mg ml −1 L-Proline, 50 mg ml −1 ascorbic acid, 0.1 mM sodium pyruvate, 0.1 % v/v ITS premix and 10 ng ml −1 TGF-β3. All materials were procured from Sigma Aldrich, except ITS premix and TGF-β3 which were procured from BD Biosciences (HQ, NJ, USA) and Lonza (HQ, Basel, Switzerland) respectively. The multi-layered scaffolds were cultured for four weeks prior sectioning and image acquisition. Detailed descriptions on scaffold preparation and isolation of bovine chondrocytes can be found in earlier works [44][45][46][47][48] .
Hyperspectral maps of the incubated PLLA and SIS scaffolds were collected using a Perkin Elmer Spectrum 100 FTI-IR connected to a Perkin Elmer AutoImage microscope. Sections of the scaffolds were cut in a cryostat microtome by embedding in O.C.T., and subsequently placed on MgF 2 slides. Hyperspectral data was collected with 10 integrations across the spectral region from 1000 cm −1 to 4000 cm −1 , at a spatial resolution of 20 µm. Spectral data in the region from 1000 cm −1 to 1860 cm −1 was used for further analysis.

Coherent Anti-stokes Raman Scattering (CARS) microspectroscopy of immiscible lipid-water phases in mayonnaise
Immiscibility study of the lipid and aqueous phases in mayonnaise using the Chemical Dimensionality and endmember extraction algorithm was conducted using Coherent Anti-stokes Raman Microspectroscopy. The imaging setup included a Leica TCS SP8 inverted microscope (Leica Microsystems GmbH, Wetzlar, Germany) equipped with a pi-coEmerald pico-second pulsed laser (APE GmbH, Berlin, Germany) delivering a 300 mW Stokes beam at 1064 nm, and a tunable optical parametric oscillator (OPO) delivering a pump/probe beam in the range 780 nm to 940 nm with 300 mW. The Stokes and pump/probe lasers were focused onto the sample using an HC PL IRAPO 40x/1.1 objective, and a spectral stack was collected using the non-descanned detectors in forward mode by tuning the pump/probe laser in the range from 794.6 nm to 826.6 nm in steps of 1 nm (∼3200 cm −1 to 2713 cm −1 ). Handmade mayonnaise was prepared at the SDU FoodLab, and kindly donated by Mathias Porsmose Clausen.

Overview of the CD framework
A detailed flowsheet of the Chemical Dimensionality framework, CD-DE and CD-EEA, is given in Supplementary  Figures S1 and S2. Briefly, hyperspectral datasets are subset to candidate endmember matrices by repeated simplex maximization in an increasing number of dimensions. Consequently, the same pixel/sample may be considered an endmember candidate in several such hyperdimensional spaces. The candidate endmember matrices are used as input to multiple semi-nonnegativity constrained matrix factorization models, a type of blind signal separation technique, where each model assumes a successively increasing number of source signals (Fig. 1b). Dimensionality is estimated from these blind signal separation models based on the information entropy (Shannon entropy 20 ) of the reconstruction residuals and the reduction in sum of squared residuals for each model (error reduction). Individually these methods may be used with thresholding to estimate the dimensionality. However, the selection of threshold is often subjective. For both CD-DE and CD-EEA, information entropy and error reduction provides complimentary information on model performance based on metrics describing residual information in the reconstruction residuals (Shannon entropy) and a measure of model fit to the data (reduction in sum of squared residuals).

Benchmarking CD-DE to dimensionality estimation techniques from remote sensing
First we benchmarked the developed CD-DE to established methods in analytical chemistry and satellite remote sensing, hereinafter, collectively denoted Virtual Dimensionality methods; the Harsanyi-Farrand-Chang (HFC) method 21 based on Neyman-Pearson detection theory, the Akaike Criterion (AIC) and Minimum Description Length (MDL) 22 based on information theory, the Malinowski Factor Indicator Function (FIF) derived from factor analysis 23 , and last, dimensionality estimated from the number of principal components required to explain at least 95 % and 99 % of the total variance, denoted λ 95 and λ 99 respectively. Initially we sought to benchmark CD-DE on its ability to handle datasets in which the number of samples/pixels n is less than the number of channels/wavenumbers p, approximately equal n ≈ p, and greater than the number of channels n > p. We generated pseudo-spectral datasets with a known number of endmembers, ranging from two to ten endmembers and p = 1001 channels, vector normalized to unit length, convexly mixed in silico to n = 300, 1000 and 5000 samples and added white Gaussian noise to a signal-tonoise ratio (SNR) of 1000:1 to avoid strict linear dependencies ( Table 1 in Methods).
We found for n = 300 samples and p = 1001 channels, the developed CD-DE was the only method estimating the correct dimensionality (Supplementary Table S1). AIC, MDL and FIF were unable to accurately estimate the dimensionality for datasets where the number of samples was approximately equal or less than the number of channels (Supplementary Tables S1 and S2), while the explained variance thresholding measures, λ 95 and λ 99 generally underestimated the dimensionality for n = 300, 1000 and 5000, but was consistent in its estimations (Supplementary Tables S1, S2 and S3). In all three cases, n = 300, 1000, 5000 the developed CD-DE method correctly and consistently estimated the correct dimensionality based 30 repeats of CD-DE algorithm (Supplementary Tables S1, S2 and S3), as the only one of the tested methods (full confusion matrix in Supplementary Table S4). Based on Neyman-Pearson Detection Theory using a false-alarm probability of P F = 10 −5 , the HFC method was able to estimate the dimensionality for all three cases, but generally overestimating the dimensionality for n = 300 and 1000 (Supplementary Tables S1 and S2).
We then tested whether the accuracy of CD-DE was contingent on the size of the candidate endmember matrix. Using the mixture dataset containing n = 5000, p = 1001 at an SNR of 1000:1, we varied the number of hyperdimensional simplex maximization repeats used to populate the candidate endmember matrix. We found no effect of varying size of the candidate endmember matrix by volume maximization of hyperdimensional simplexes up to a maximum of 14, 19 and 29 dimensions (Supplementary Table S5), translating to matrices of 119, 209 and 464 endmember candidates.

Benchmarking CD-DE to Virtual Dimensionality methods with respect to signal noise
We then examined the robustness of CD-DE in response to increasing noise levels, benchmarked to existing Virtual Dimensionality methods. Ground-true endmembers remained the same as in prior analysis, however, we generated four new noised datasets, each consisting of nine constituent datasets containing between two and ten endmembers, that were mixed by convex mixing in silico to n = 5000 samples, and added white Gaussian noise to a signal-to-noise ratio of SNR = 100:1, 50:1, 20:1 and 10:1 respectively (see Methods).
CD-DE was unaffected by increasing noise levels for SNR=1000:1, 100:1, and 50:1 (Supplementary Tables S3, S6 and S7), correctly and consistently estimating the ground true number of endmembers across all 30 repeats of the algorithm confer confusion matrices in Supplementary Tables S4, S10 and S11. The same was true for AIC, MDL, FIF and HFC at all tested noise levels, however, HFC only yielded the ground true number of endmembers in four out of nine cases. λ 95 and λ 99 were consistent, but generally underestimated the number of endmembers for SNR 100:1 and 50:1.
For SNR 20:1 the CD-DE estimated the correct dimensionality for all endmember datasets except the dataset containing three endmembers, and with minute variability for the dataset containing nine endmembers, for which 29 of 30 repeats estimated the correct dimensionality (Supplementary Table S8). At SNR 20:1 the λ 99 began to significantly overestimate the dimensionality, and for SNR 10:1 both λ 95 and λ 99 greatly overestimated the number of endmembers in all but one case ( Table 2, and Supplementary Table S8). CD-DE remained generally consistent, only with variability in the dataset containing three endmembers, confer confusion matrices in Supplementary Tables S12 and S13. The estimated dimensionality was correct in 7 out of 9 cases, and never deviated more than one endmember from the ground truth ( Table 2).

Benchmarking of CD-EEA endmember extraction
We then investigated if the developed framework could function as an endmember extraction algorithm (CD-EEA). The goal of an endmember extraction algorithm (EEA) is to retrieve the "purest" feature vector/spectrum for each distinct signal source, that combined describes the remaining mixture samples in the dataset. Using the dimensionality estimated by CD-DE, k CD , CD-EEA extracts endmember spectra that 1) lowers the squared 2 norm, and 2) if possible, maximizes the entropy of the reconstruction residuals (Fig. 3a, Methods). To evaluate the performance of the CD framework as an endmember extraction method, we initially benchmarked CD-EEA against N-FINDR, based on synthetic datasets containing 2, 5 and 10 endmembers vector normalized to unit length, and convexly mixed in silico to 5000 samples with SNR = 1000:1. Out of the 5000 spectra, both CD-EEA and N-FINDR succeeded at recovering all ground true endmembers in all three cases (Supplementary Figures S3 to S5), suggesting that CD-EEA indeed can function as an endmember extraction algorithm.
To further substantiate the ability of CD-EEA as an endmember extraction technique, we collected Raman spectra of seven reference compounds; polystyrene, tylenol, silicon, calcium carbonate, naphthalene, stearic acid and β-tricalcium phosphate, varying especially by the complexity of their spectral signatures (Figure 3c). The seven measure references were mixed convexly in silico to 5000 samples and added white Gaussian noise to SNR = 1000:1 to avoid strict linear dependencies. The Chemical Dimensionality algorithm correctly estimated k CD = 7, as indicated by both the error reduction and associated increase in reconstruction entropy (Figure 3d-f). Although the CD-EEA and N-FINDR returns endmembers unordered, they both uncovered the ground true endmember spectra out of the 5000 mixture spectra (Figure 3a-b).

Combined CD-DE and CD-EEA image reconstructions for imaging oil-in-water emulsions by Coherent Anti-stokes Raman Scattering
Having benchmarked CD-DE and CD-EEA for estimation of data dimensionality and endmember extraction, we sought to study more complex (bio)chemical system. We first studied liquid immiscibility in mayonnaise, a model oilin-water emulsion, using Coherent Anti-stokes Raman Scattering (CARS). Contrary to most other microspectroscopy techniques, CARS signals are not easily used quantitatively owing to an inherent non-resonant background arising from coherent addition of off-resonant vibrational contributions, induced by the pump and Stokes fields 24 , challenging both estimation of dimensionality and extraction of endmembers.
We compared endmembers extracted by CD-EEA to N-FINDR and cluster centroids from K-Means++ clustering, along with their associated image reconstructions, according the k CD estimated by CD-DE. For the CARS hyperspectral datacube of mayonnaise estimated dimensionality was k CD = 2 (Figure 4a-c), a reasonable estimate given that mayonnaise is an oil-in-water emulsion. However, AIC, MDL, λ 95 and λ 99 all estimated 26 endmembers or more (Table 3).
Remarkably, CD-EEA algorithm and N-FINDR both returned the same two endmembers given k CD = 2 (Figure 4d-e). Endmember spectra were generally characterized by a background spectrum corresponding to the water-phase, and an endmember with peaks around 2847 cm −1 and 2937 cm −1 , assigned to C-H and CH 3 stretching vibrations 25,26 , consistent with the presence of lipid. The cluster centroids from K-Means++ indicates lipid in both endmembers, which may occur if pixels containing lipid is ascribed to both clusters (Figure 3f). Nonetheless, the K-Means++ algorithm succeeds at revealing the lipid-droplet microstructure (Figure 4i), as expected from an oil-in-water emulsion and previously demonstrated using both CARS and fluorescence microscopy 27,28 . The same microstructure is observed in image reconstructions using endmembers returned by CD-EEA and N-FINDR (Figure 3g-h).

CD-DE and CD-EEA for monitoring tissue-engineered extracellular matrices with FTIR microspectroscopy
Label-free microspectroscopy plays a pivotal role in research and development of engineered tissue grown in vitro.
Monitoring cell growth, differentiation and migration in response to physical and mechanical cues served by the biomaterial scaffold is paramount for engineering personalized, biocompatible tissue replacements. Label-free microspectroscopy can offer a wealth of information with molecular specificity, however, the success hyperspectral imaging for monitoring is intimately related to the estimation of dimensionality and the ability of the endmember extraction algorithm to distinguish mixed spectral information from spectral information corresponding to an endmember.
We engineered a scaffold from electrospun fibers of Poly-L-Lactic acid (PLLA), and a scaffold from Small Intestine Submucosa (SIS) derived from porcine jejunum. Both scaffolds were incubated with isolated primary bovine chondrocytes and transforming growth factor β3 (TGF-β3). After 4 weeks, scaffolds were imaged using FTIR microspectroscopy and analyzed using CD-DE and CD-EEA.
The estimated chemical dimensionality for the PLLA and SIS scaffolds were k CD = 5 and k CD = 4 respectively ( Table 3). Endmember extraction on the SIS-scaffold FTIR dataset yielded a similar endmember for the collagenous extracellular matrix as for the PLLA scaffolds, using both N-FINDR (endmember 4; green) and CD-EEA (endmember 3; green) for endmember extraction (Figure 5d-e). Absorptions at 1652 cm −1 , 1548 cm −1 , and 1236 cm −1 were assigned to amide I, amide II and amide III respectively, likewise displaying absorptions at 1080 cm −1 and 1032 cm −1 , that were assigned to C-O and C-O-C vibrations stemming from carbohydrate species in the tissue 29 . Endmembers extracted are generally consistent between the two methods, however, CD-EEA discovers and extracts a spectrum containing spectral signatures at 1736 cm −1 and 1472 cm −1 (endmember 2; cyan; Figure 5d) in addition to the ECMrelated amide signatures. The spectrum is not returned by N-FINDR (endmember 1; cyan; Figure 5e). The absorptions were assigned to carbonyl stretching (C=O; 1736 cm −1 ) and methylene scissoring (C-H; 1472 cm −1 ), suggesting the presence of lipid 30 in the submucosal layer of the small intestine. The associated endmember reconstructions indicate little co-localization between CD-EEA endmember 2 and the remaining endmembers (Figure 5g). Thus, this endmember resembles a lipid-rich region of the tissue. From the viewpoint of simplex volume maximization, N-FINDR returns a background spectrum (endmember 1; cyan; Figure 5b), confirmed by image reconstruction with the endmember in question (Figure 5h). Inspection of the cluster centroids from clustering with K-Means++ mainly indicates the presence of the three amide bands, confer findings with CD-EEA and N-FINDR, but are also indicative of weak signatures at 1736 cm −1 in all four endmembers (Figure 5f). The associated cluster membership images do not provide detailed microstructure of the tissue scaffold, and cluster centroids generally does not convey the molecular specificity offered by CD-EEA and N-FINDR (Figure 5d-e). More elaborate pre-processing would typically be required, such as background removal and extended multiplicative scattering corrections including spectral interference subtraction, to provide enhanced chemical and structural information. To visualize the contrast in structural and chemical informa-

Discussion
We have developed an entropy-guided dimensionality estimation (CD-DE) and endmember extraction algorithm (CD-EEA) for (bio)chemical hyperspectral imagery. Our findings suggest a combination of simplex volume maximization, semi-nonnegativity constrained matrix factorization and information entropy, can both estimate the dimensionality of hyperspectral datasets across various types of microspectroscopy techniques and serve as basis for endmember extraction. To our knowledge, this is the first time that dimensionality estimation and endmember extraction is achieved using blind signal separation and Shannon entropy, across quantitative spectral information in spontaneous Raman and FTIR, as well as in CARS containing an inherent non-resonant background.
Estimated dimensionalities using the CD-DE were compared to established methods in satellite remote sensing, including the HFC, AIC, MDL, Malinowski FIF and canonical explained variance thresholding measures from analytical chemistry. CD-EEA endmember extraction for CARS, FTIR, spontaneous Raman and synthetic datasets was benchmarked to endmembers extracted by N-FINDR, and to cluster centroids from K-Means++ clustering of the data. In comparison to prior art 19 , we relaxed the mixing requirements from convex-to conical mixing, and exploited relaxation of the non-negativity constraints on the self-modelling curve resolution models to account for weak negativity in intensities, for instance occurring as a part of the preprocessing of the hyperspectral datasets with baseline removal algorithms. We especially noted that although all Virtual Dimensionality methods were able to estimate dimensionalities of all experimental datasets, the estimates were for certain experimental datasets unreasonably high, predicting the presence of more than 100 endmembers in a single hyperspectral dataset (Table 3). Consequently, we argue that CD-DE is less dependent on i) the type of imaging/spectrometry and ii) the dimensions of the unfolded hyperspectral datamatrix, i.e. the number of samples and channels available in the datacube, than other Virtual Dimensionality methods.
The generic applicability offered by entropy-guided estimation of dimensionality and endmember extraction, as demonstrated by the plurality of samples and hyperspectral imaging methods used for demonstrating algorithm, serves as illustration of how CD-DE and CD-EEA can prove to be a powerful for interrogating chemical and biochemical systems, from gastrophysics to tissue engineering, especially by consolidating hyperspectral molecular information with sample microstructure.

Supplementary tables
The CD-DE algorithm was compared performance-wise to established Virtual Dimensionality methods, under scenarios where 1) the sample size of the input datamatrix varied from 300-5000 samples (SNR = 1000:1), and 2) signal-tonoise ratios from 1000:1 to 10:1 for n = 5000 samples. For the CD-DE algorithm alone, was the effect of varying the simplex maximization repeats (g) also tested, i.e. the number of samples carried over from the input datamatrix to the reduced candidate endmember matrix V. All benchmarking was carried out on the synthesized datasets containing a known number of endmembers, from two to ten endmembers, sufficient for most uses of microspectroscopy.

FTIR microspectroscopy of PLLA scaffold
Using the PLLA and SIS scaffolds incubated with isolated primary bovine chondrocytes and TGF-β3, we estimated the Chemical Dimensionality using the CD-algorithm and observed the endmembers returned from CD-EEA and N-FINDR, as well as the cluster centroids from K-Means++. The estimated chemical dimensionality for the PLLA and SIS scaffolds were k CD = 5 and k CD = 4 respectively.
While the maximum error reduction coincided with an increase in entropy for the SIS scaffold (Figure 5b-c in the paper), this was not the case for the PLLA scaffold. Here, the semi-nonnegativity constrained matrix factorization models containing between two and four feature vectors were associated with consecutively decreasing reconstruction entropies, despite the maximum error reduction was observed for the model containing two feature vectors.
An increase in entropy was first observed for k CD = 5, that noteworthily, also resulted in a peak, albeit smaller compared to the maximum, in the error reduction ( Figure S6 in Supplementary Information). This raised the question whether there were grounds for estimating a higher dimensionality based on the reconstruction entropy. Indeed, if two endmembers were extracted to represent a single endmember species de facto, one would expect close colocalization between the two endmembers. However, as observed in Figure S7d, the brightest pixels of endmember 5 (cyan), are spatially exclusive to the pixels associating with endmember 4 (magenta). Inspecting the corresponding endmember spectra reveals spectral differences in the low-wavenumber region from 1100 cm −1 to 1000 cm −1 ( Figure S7a) across endmembers 2 (red), 4 (magenta) and 5 (cyan), assigned to C-O stretching vibration in PLLA 1 and glycosaminoglycans 2 . Nonetheless, the remaining spectral features corresponding to PLLA (1756 cm −1 , C=O; 1452 cm −1 , 1384 cm −1 , 1364 cm −1 , C-H) and the extracellular matrix (ECM; 1652 cm −1 , amide I; 1548 cm −1 , amide II; 1236 cm −1 , amide III) are present in each of the three endmembers ( Figure S7a). Taken together, these observations suggest the algorithm is sensitive toward narrow-band spectral dissimilarity, even if features in remaining parts of the spectrum correlate to other endmembers. Figure S1: Algorithm flowchart for estimating the Chemical Dimensionality. Off-page connector labelled 'A' directs to the endmember extraction sub-algorithm.         [a] Based on 30 repeats of the Chemical Dimensionality (CD) algorithm. Full confusion matrix in Table S4.

Supplementary tables
[b] False alarm probability of 10 −5 . [a] Based on 30 repeats of the Chemical Dimensionality algorithm. Confusion matrix in Table  S4.
[b] False alarm probability of 10 −5 .     [a] Based on 30 repeats of the Chemical Dimensionality algorithm. Confusion matrix in Table  S10.