Unsupervised clustering of multiparametric fluorescent images extends the spectrum of detectable cell membrane phases with sub-micrometric resolution

: Solvatochromic probes undergo an emission shift when the hydration level of the membrane environment increases and are commonly used to distinguish between solid-ordered and liquid-disordered phases in artiﬁcial membrane bilayers. This emission shift is currently limited in unraveling the broad spectrum of membrane phases of natural cell membranes and their spatial organization. Spectrally resolved ﬂuorescence lifetime imaging can provide pixel-resolved multiparametric information about the biophysical state of the membranes, like membrane hydration, microviscosity and the partition coeﬃcient of the probe. Here, we introduce a clustering based analysis that, leveraging the multiparametric content of spectrally resolved lifetime images, allows us to classify through an unsupervised learning approach multiple membrane phases with sub-micrometric resolution. This method extends the spectrum of detectable membrane phases allowing to dissect and characterize up to six diﬀerent phases, and to study real-time phase transitions in cultured cells and tissues undergoing diﬀerent treatments. We applied this method to investigate membrane remodeling induced by high glucose on PC12 neuronal cells, associated with the development of diabetic neuropathy. Due to its wide applicability, this method provides a new paradigm in the analysis of environmentally sensitive ﬂuorescent probes.


Introduction
Cell membranes are composed of many different types of lipids, including saturated and unsaturated phospholipids, cholesterol, and fatty acids. Membrane lipids can interchange among several membrane phases, which are supramolecular structures having relatively uniform chemical and physical properties. At low temperatures, artificial lipid bilayers composed by a limited number of phospholipids are laterally ordered and well-packed together in the membrane, in a solid ordered, gel-like phase (s o ). Raising temperatures over a transition value (melting temperature, Tm), characteristic of the phospholipid composition, induces the transition to a fluidlike, liquid-disordered phase (l d ), which is characterized by a larger area occupied by polar heads of phospholipids. Sufficiently high membrane sterol concentration triggers the emergence of the liquid-ordered phase (l o ), which possesses solid-like qualities similar to s o though maintaining the high rate of lateral diffusion of the l d phase. Cell membranes display more complex behavior, consisting of a mosaic of multiple phases coexisting in the lipid bilayer [1,2]. The heterogeneity of cell membranes cannot be simplified to mere differentiation of coexisting lipid domains since the number of lipid species is significantly higher with respect to artificial bilayers, and high local concentrations of membrane-associated proteins create lateral pressure that renders lipid domain separation more thermodynamically favorable, possibly leading to phase transitions or to the emergence of new phases [2]. The resulting membrane phase dynamic equilibrium solution of Laurdan (Molecular Probes, Inc., Eugene, OR, USA) is added per milliliter of medium and cells are treated with 11 mM, 50 mM and 75 mM glucose for 48 h before labeling. For organelle labeling experiments, 1 µl of 1 mM stock solution of Nile Red (Molecular Probes, Inc., Eugene, OR, USA) is added per milliliter of medium and incubated for at least 30 minutes in the dark for staining intracellular lipid droplets (LD). 40 µl of different fluorescent protein-based markers, CellLight ER-GFP BacMam 2.0, CellLight Golgi-GFP BacMam 2.0 and CellLight Mitochondria-GFP BacMam 2.0 (Molecular Probes, Inc., Eugene, OR, USA), are used for the staining of the endoplasmic reticulum (ER), Golgi apparatus and mitochondria, respectively. Incubation is performed for 16 hours at 37°C.

Image acquisition settings
Laurdan emission intensity images are acquired with a Nikon A1-MP confocal microscope equipped with a 60× oil-immersion objective (1.4 NA), a 2-photon Ti:Sapphire laser (MAI TAI DEEP SEE EHPDS -007 (2017), Spectra Physics, Newport Beach, CA) producing <70-fs pulses at a repetition rate of 80 ± 1 MHz. A PML-SPEC 16 GaAsP (B&H, Germany) multi-wavelength detector (grating part is the No. 77414) coupled to a SPC-830 TCSPC/FLIM device (B&H, Germany) working in FIFO mode is used to collect the decay data. Calibration was performed through the registration of second harmonic generation (SHG) signals registered at several excitation wavelengths, as recommended by the manufacturer. The software used for acquisition is SPCM-64 (B&H, Germany). Laurdan fluorophore is excited at 780 nm. This value constitutes a good compromise between the maximization of absorption and the reduction of autofluorescence excitation. Signals are integrated into the time window of width 12.5 ns and into the wavelength region 400-600 nm. For image acquisition, the pixel frame size is set to 512 × 512 and the pixel dwell time is 60 µs. Imaging is performed at 37°C with 5% CO2 (cage incubator OKOLAB). The average laser power at the sample is maintained at the mW level (10 mW). For FLIM, data are acquired, without performing time integration Signals are integrated into the wavelength region 420-460 nm for the blue channel (I blue (t,x,y)) and 500-560 nm for the green channel (I green (t,x,y)). FLIM acquisitions provide for the decay curves, for each pixel.

Fluorescence intensity microscopy for the partition coefficient determination:
construction of the P-image Fluorescence intensity of a lipophilic probe in the membrane phase is proportional to the partition coefficient P, defined as the ratio of the emission intensity of the probe between the phase at pixel (x,y) and a reference phase: Log[P(x, y)] = log I x,y I ref (1) where I x,y is the emission intensity of Laurdan at the pixel(x,y) integrated on the whole emission spectrum, and I ref is the emission intensity of Laurdan in a reference phase. The reference phase is the phase characteristic of the core of LD, which is characterized by the highest uptake of the Laurdan probe (Fig. S1), as already qualitatively observed in [25]. Therefore, the log P value we defined herein is a measure of the hydrophobicity of the membrane taking the LD as a reference. The use of LD's Laurdan intensity as reference for maximum has been chosen accordingly to the high compartmentalization of the probe in these organelles, as also observed in a previous article [25]. In the Supplement 1, the detailed workflow for the determination of the reference is described. P value is dependent on the composition and the physical state of the membrane phase, and especially on the distribution, length and three-dimensional organization of the hydrophobic tails. The P value is calculated pixel-by-pixel using the equation [1] yielding the P image, which is one of the three image inputs for the unsupervised analysis. Background values are measured and subtracted for each image, and debris or other aggregates are removed to avoid biases in the analysis. Background values are removed using the software SPC-image 8.1 (B&H, Germany) through the function 'Threshold'. This function defines the minimum number of photons in the peak of a fluorescence curve. Pixel with lower photon number are not analyzed by the fitting procedure, therefore suppressing dark pixels. This value is sampled on n=100 curves on the area outside cells, ranging from 5 to 10.
2.4. Fluorescence ratio imaging microscopy for membrane hydration determination: construction of the GP image For the fluorescence ratio imaging approach, the same experimental setup and the same acquisition data are used, but signals are integrated in the time window of width 12.5 ns and into the wavelength region 420-460 nm for the blue channel and 500-560 nm for the green channel to evaluate the spectral shift of the probe. The fluorescence spectrum of the probe Laurdan, which incorporates into the lipid phase in the membrane, is correlated to its physical state. Laurdan's excited-state relaxation is highly sensitive to the presence and mobility of water molecules within the membrane bilayer, yielding information on membrane hydration by a shift in its emission spectrum depending on the surrounding lipid phase state (i.e. bluish in ordered, gel phases and greenish in disordered, liquid-crystalline phases) [26]. By using this probe, coexisting lipid domains are labeled according to their distinctive fluorescence spectra and dual-wavelength ratio measurements. As a normalized ratio of the intensity at the two emission wavelength regions, the generalized polarization (GP) provides a measure of membrane order (high membrane micropolarity, low GP; low membrane micropolarity, high GP). The GP index is calculated for each pixel using the two Laurdan intensity images I blue (x,y) and I green (x,y) [27], according to the following formula: G is a correction factor calculated by using the following Eq.
where GP theo is the GP theoretical value of a standard solution of 5µM Laurdan in DMSO, which has a known value (GP theo =0.207), where-as GP exp is the GP value of the same solution measured in our confocal microscope [28]. GP index is independent of excitation intensities, probe concentrations, and other artefacts, relying on the ratiometric properties of the probe [6,15].

Fluorescence lifetime microscopy for microviscosity determination: construction of the η image
To determine the lifetime, especially when the decay time exceeds 30% of the 12.5 ns time window, a good estimate of τ green and τ blue is retrieved by fitting the data with an incomplete-decay model which includes the fluorescence remaining from the previous laser pulses in the model function. This is the case of Laurdan lifetime as in [29]. The model is a sum of exponential terms, taking into account the fluorescence that does not fully decay within a single laser pulse period. Based on 'Repetition Time' the fluorescence leftover from all previous excitation pulses is included in the model [30]. The 'offset' is kept fixed and set to zero. With the PML-16-GaAsP detector, this situation is achieved by strictly avoiding pickup of room light. This was obtained by using obscuring panels applied to the OKOLAB cage incubator, thus realizing total dark in the acquisition chamber. The value of the background was measured for our setup and the average photon count for each pixel on a 512 × 512 image is 0.48 ± 0.8. In membranes, the solvent relaxation (decrease in GP), analyzed so far, depends on the number of the surrounding water molecules, and thus increases with membrane hydration levels [31]. This effect can be also revealed by a decrease of the fluorescence lifetime due to the enhanced emission from the relaxed state, which is characterized by a lower lifetime. Another quantity that can be measured analyzing the fluorescence decay is the rate of solvent relaxation, obtained by measuring the rate of the spectral shift [31,32]. The speed of solvent relaxation is related to the rotational mobility of the water molecules within the membrane and is referred to as membrane microviscosity [29]. When the rate of spectral relaxation is in the picosecond scale, the temporal resolution of the time-correlated single photon counting devices is not able to resolve any changes. However, in very viscous environments the time scale of relaxation can increase up to the nanosecond time scale. In this context, the change in membrane polarity due to water hydration from the effects due to viscosity can be uncovered by analyzing the decay time of the blue channel (emission: 450/50 nm) and the green channel (emission: 540/50 nm) [32]. The excited-state decay in the blue channel presents an apparent decay time that reduces the standard decay time of the fluorescence emission from the relaxed state because the depopulating process from the locally excited state to the relaxed state contributes to the overall decay. Conversely, the excited-state decay in the green channel increases due to the populating process from the locally excited state to the relaxed state. Overall, the populating and depopulating processes due to the solvent relaxation lead to a contextual decrease of the fluorescence lifetime in the blue channel, and an increase of the measured lifetime in the green channel. Thus, the ratio between the blue and green apparent lifetimes provides for a measure of the extent of the solvent relaxation. The rationale of using the ratio is to summarize a common effect due to solvent relaxation in the two channels: from previous studies [29,[31][32][33], as already explained upwards, the main effect of solvent relaxation is to increase the lifetime of the green channel and decrease concomitantly the lifetime of the blue channel.
The ratio allows to reduce the extent of lifetime variations due to fluidity (which are already accounted by GP). Amplitudes are more strongly affected by fluidity variations. We added this explanation and this reference to justify this choice. This parameter, independent from the probe concentration, uneven illumination, and other artefacts gives a further characterization of the membrane phase, thus providing information about its microviscosity. An increase in η(x, y)Z indicates a slower spectral relaxation, which is related to slower motility of water molecules due to the formation of hydrogen bonds with a membrane phase that presents a high grade of structured water, increasing the microviscosity of the phase. In several works, an increase in structured water is linked to the presence of cholesterol [14,32,33]. The η value is calculated pixel-by-pixel using the equation [4] yielding the η-image, which serves as one of the inputs for the unsupervised analysis. Background values are measured and subtracted for each image, and debris or other aggregates are removed to avoid biases in the analysis.

Confocal microscopy of intracellular organelles
For the detection of intracellular organelles, petri dishes with GFP or Nile Red labeled cells are placed on the inverted confocal microscope (Nikon A1 MP) equipped with an on stage incubator (T = 37°C, 5% CO2, OKOLAB) and a 32 channel spectral detector [21]. The 32 channel spectral images are obtained using a 60× objective (NA 1.4) under 488 nm excitation in the emission range 510-700 nm (6 nm bandwidth). Internal photon multiplier tubes collected 16 bit, 1024×1024-pixel images at 0.125 ms dwell time. To isolate intracellular organelles from the background an intensity-based segmentation is applied. For the intracellular GFP-stained organelles segmentation, the intensity image is integrated into the spectral range 510-520 nm.
For the Nile Red labeled LD, the intensity image is integrated into the spectral range 534-557 nm. Organelles segmentation is obtained through the open-source software ImageJ (NIH) by applying the Max Entropy method (built-in function 'Threshold'). Plasma membranes (PM) are isolated by applying the 'Find Edges' tool.
The quantification of the fraction of area covered by a specific organelle (central sections of the cells are sampled as reference) is evaluated as the ratio of the number of pixels covered by the organelle (N pixel ORG ) and the total number of pixels covered by the cell (N pixel CELL ) Because of the possibility of an overlap between the organelles (especially ER and Mitochondria), the values given in graph represent upper estimates.

Statistics
Statistical tests for sets of biological/biophysical data ( Fig. 2(F), 2(G), 5(B), 7) are performed by R Studio (https://www.rstudio.com/). For image quantification, a minimum of n=40 cell per sample is analyzed. Baseline characteristics among groups have been compared with ANOVA for parametric variables. Quartile-quartile plots checked data normality graphically. A comparison between groups' couples has been made with Tukey's Test.

Analysis of P maps, GP maps, and η maps of PC-12 cells
The pixel resolved parameters extracted from the hyperspectral lifetime acquisitions are visualized in the P maps, the GP maps and η maps, reflecting pixel-resolved differences in the partition coefficient of the probe within the phase, in the hydration levels and in the microviscosity of the phase, respectively. In Fig. 1(A)-1(B)-1(C), P maps, GP maps and η maps of a representative, untreatedPC-12 cell are reported. In Fig. 1(D)-1(E)-1(F) the spatial distributions of normalized intensity (contour plots generated with Image-J, NIH), GP and η values are visibly different since each quantity reflects a peculiar property of the membrane phase. P maps present peaks in correspondence of LD [25], which are resolved and display the highest P value (P=1, white regions), and an intermediate value (P∼0.5, yellow regions) in the intracellular regions different from LD. Based on its location and extent within the cell (Fig. 6), these internal regions correspond to ER membranes, mitochondrial membranes and the membrane of Golgi apparatus. Perinuclear and peri-PM regions display the lowest P values. The GP distribution ( Fig. 1(E)) is more homogeneous in LD and in the regions around them. In the other intracellular regions, it is possible to observe subsets of pixels characterized by a flat spatial P distribution, but with a GP distribution dominated by a wrinkled landscape ( Fig. 1(D) and 1(E), red arrows). The η spatial distribution ( Fig. 1(F)) provides an additional informative content: the corrugated landscape is very different from the one characteristic of P and GP distributions, highlighting the presence of several membrane regions characterized by a range of diverse spectral relaxation values.

P maps, GP maps, and η maps of PC-12 cells in glucose overload conditions
In Fig. 2(A)-2(B)-2(C) representative P Maps, GP Maps and η maps of PC-12 cells, respectively at 11 mM, 50 mM and 75 mM glucose are reported. The average variations of GP and η with increasing glucose concentration (averages of the image histograms over the entire dataset) are not significant (Fig. 2(F)-2(G)). However, from the GP histogram calculated for n=40 cells (Fig.  2(D)), a rearrangement of the membrane domains following the glucose treatment, consisting of a decrease of the high GP pixels (in the dark-green colored region around GP∼0.9, low hydration) and an increase of the low GP pixels (in the light-green colored region around GP∼0.5), can be seen. Concomitantly, a slight shift of the η distribution towards lower values can be observed ( Fig. 2(E)).

Artificial intelligence clustering of multiparametric fluorescent images
The separate analysis of the three maps, though indicating that some changes are occurring, lacks the specificity needed to fully describe the phases of intracellular membranes. To extend the spectrum of detectable membrane phases, we report in Fig. 3 the workflow for obtaining AI images. In the first step, multiparametric images are acquired with a hyperspectral imaging system. For each pixel, intensity values and lifetime values for the blue channel (400-470 nm) and green channel (490-600 nm) are obtained, for a total of 4 parameters. These parameters are further processed (Materials and Methods), adding additional information in the form of two factors: the G factor and the referenced normalization of intensity (Eq. (1) and Eq.3), transforming them in the P value, the GP value and the ratio η of the lifetimes. In the second step, images are stacked, creating a multiparametric map in which each pixel is associated with a real phase-vector ì ϕ = (P, GP, η). The parameters are then rescaled and the pixels of the entire dataset analyzed are sent to a k-means classifier. K-means clustering aims to partition the n observations into k (≤ n) sets S = S 1 , S 2 , . . . , S k to minimize the within-cluster sum of squares (WCSS) (i.e. variance). The obtained cluster center is the representative of its cluster. The number of observations n corresponds to the total number of pixels analyzed (n=512 × 512*m, where m is the number of images).
All the m images of the experiment (m∼ 15, resolution 512 × 512) are sent to the classifier. In each image from 6 to 10 cells are present. The k-means algorithm [34,35] represents each of the k clusters Cj by the mean (or weighted average) cj of its points (centroid). The sum of distances between elements of a set of points and its centroid expressed through an appropriate distance function is used as the objective function. We employed the L2 norm-based objective function, i.e. the sum of the squares of errors between the points and the corresponding centroids, which is k-means iterative optimization consists of two-step major iterations that: (1) reassign all the points to their nearest centroids, and (2) recompute centroids of newly assembled groups. Iterations continue until a stopping criterion is achieved (no reassignments with tolerance <10 −5 ). This version, known as Forgy's algorithm [35] works with any Lp norm and it does not depend on data ordering. The k-means algorithm also has certain shortcomings which are minimized or mitigated according to the following: 1) Dependence from the initial guess of centroids: To mitigate the effects of cluster initialization we applied the method suggested by Bradley and Fayyad [36]. First, we performed k-means on several small samples of data with a random initial guess. Each of these constructed systems is then used as a potential initialization for a union of all the samples. Centroids of the best system constructed this way are suggested as an intelligent initial guess to ignite the k-means algorithm on the full data.
2) Computation of local minima: No initialization guarantees a global minimum for k-means. This is a general problem in combinatorial optimization, which was tackled by allowing uphill movements and determines the global minimum of a given set of n-dimension functions with constraints using simulated annealing (SA) method [37].

3) Problem of k-selection:
The optimal number of clusters (k=8) is determined by the elbow method and validated with the Silhouette method (Fig. 4). The elbow method consists in Fig. 3. Workflow for the artificial intelligence clustering of multiparametric fluorescent images. In the first step, multiparametric images are acquired with a hyperspectral imaging system. For each pixel, the P value, the GP value, and the ratio η are measured to quantify the spectral properties of the solvatochromic probe. In the second step, images are stacked, creating a multiparametric map in which each pixel is associated with a 3D phase state-vector (P, GP, η). The dataset of phase state vectors is then sent to a k-means classifier. K-means clustering aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. Each classified pixel is then remapped to the pseudo-image, that we called AI image. Phases are then ranked and pseudo-colored and the computed distribution of the phases is reported in the phase-profile graph.
the determination of the within-cluster-sum of squared errors (WSS) for different values of k, and in the selection of the k value for which WSS becomes first starts to diminish. In the plot of WSS-versus-k, this is visible as an elbow (Fig. 4). Silhouette is a method to validate consistency within clusters of data. The silhouette value is a measure of how similar an object is to its cluster compared to other clusters. The silhouette ranges from −1 to +1, where a high value indicates that the object is well matched to its cluster and poorly matched to neighboring clusters. The silhouette plot displays a local maximum at the chosen k=6 and has an average Silhouette value of 0.88 (Fig. 4).

4) Sensitivity to outliers:
The parameters are rescaled using as an offset the mean value of the distribution, and as scaling factor αs, where s is the standard deviation of the distribution and α a tunable parameter. If the distributions are almost Gaussians, α >3 ensures that more than 99% of the values are considered. This was verified by visual inspection and qq-plots. When there are significant deviations from Gaussians the α factor is adjusted to include at least 99% of values by direct calculation.
Once the optimization is performed, each pixel is classified in 8 classes, with value going from 1 to 8, and then remapped to the pseudo-image, that we called AI image. Phases are ranked according to the || ì ϕ|| value. The distribution of phases is calculated by reporting the fraction of pixels belonging to each class k.
We selected P, GP and eta as input parameters because they provide a referenced and thus comparable measure of probe concentration (P) and phase shift (GP) among different samples, from experiment to experiment, and from lab to lab. Using these data, instead of the raw data,  1)). Eta is intrinsically independent from these detrimental fluctuations, and it was introduced to summarize a common effect due to solvent relaxation in the two channels: in several studies [29,[31][32][33], it is clear that the effect of solvent relaxation is to increase the lifetime of the green channel and decrease concomitantly the lifetime of the blue channel. This ratio allows to reduce the extent of lifetime variations due to fluidity (which are already accounted by GP). The usage of processed data allows also employing in k-means clustering three parameters instead of four, thus reducing the number of parameters employed in the model (if the number of dimensions (n) increases to four, with k=8 cluster the number of free parameters (equal to n×k) will increase from 24 to 32).
Keeping all the settings fixed, the usage of raw data shows a convergence of the algorithms nearly on the same clusters, but the 4p clustering is less optimal than the 3p clustering (Silhouette value=0.81). Therefore, we conclude that our approach based on processed quantities has to be preferred, allowing comparison on the same clusters of results coming from different experimental settings, different labs, different days, long acquisition times. Moreover, the maximum Silhouette value denotes an overall reduction of noise: while noise can be added in retrieving indirect quantities, the increase of specificity of these parameters and the removal of noise related to intensity fluctuations and other spurious fluctuations can lead to optimization of data clustering.

AI images of PC-12 cells undergoing glucose overload and characterization of the phase profiles
In Fig. 5(A), AI images, obtained as described in the previous section, are reported. The retrieved phases are labeled by a number going from 1 to 8. In Fig. 5(B), the average fraction of the phase is reported. Each class is characterized by the mean values of (P, GP, η) calculated on the respective clusters. These values are shown in Fig. 5(C), ranked in four classes (white = fourth quartile-lowest values, light blue = third quartile, blue = second quartile, dark blue = first quartilehighest values). While clusters 1 and 2 correspond to background values (BK 1 -BK 2 ), class 3 (P, GP, η) =(Q3-Q3-Q3) represents a fluid phase with low viscosity and low permeability. This phase is localized on the cell PM, the fraction occupied is ∼5% and it is organized in macro-domains. Since fluidity and low viscosity are characteristic of the liquid-disordered phase [2,3], we label this phase as l PM d . Class 4 (Q3-Q4-Q1) represents a hyper-fluid, hyper-viscous phase with low permeability. Also, this phase is mainly localized on the cell PM, organized in macro-domains, with an occupied fraction ∼5%. Since hyper-viscosity is characteristic of the liquid-ordered phase [29,32], we define phase 4 as l PM o phase. The sum of the area fractions occupied by the two phases is coherent with the area fraction occupied by PM ( §3.5 and Fig. 6-7). Class 5 (Q2-Q3-Q2) represents a fluid and viscous phase, with intermediate permeability, localized in the inner part of the cell. This phase is organized as a region formed by a limited number of connected blocks with a networked structure and represents the larger fraction of the cell (Fig.  4(B)). In §3.5 and Fig. 6-7 a comparison with the spatial distribution of the major intracellular organelles of PC-12 cells labeled with ER-GFP, Golgi-GFP and Mitochondria-GFP constructs and Nile Red is reported. The spatial organization and the area fraction covered by this phase indicates that it is characteristic of some portions of ER, the Golgi and the mitochondrial network. The phase has the same fluidity of the liquid-disordered phase [2,3] but has a higher viscosity. Thus it is defined as l NET d . Class 6 (Q2-Q2-Q1) represents probe-permeable, rigid, hyper-viscous domains. These membranes constitute a network of spread structures, which are less connected than class 5. Class 6 may constitute another phase of the above-mentioned networked structure, which includes vesicular membranes like endosomes and synaptic vesicles (Fig. 6). Apart from the different morphology which comes out from k-means analysis, further investigation is required to understand exactly which organelles take part to class 5 or class 6. This is not a straightforward since our technique requires the exploitation of a wide portion of the visible spectrum making it difficult to separate the fluorescence contributes. The hyper-viscosity of this phase is characteristic of the liquid-ordered phase, thus we define this phase as l NET o . Differently from l PM o , this phase is less hydrated, thus indicating a higher lipid packing. Class 7 represents hyper-permeable, hyper-rigid, viscous domains (Q1-Q1-Q2). This class is contiguous to class 5 and class 8 (LD). Therefore, this region can be localized in the ER, in correspondence of the LD formation sites. These quantities are characteristic of the gel-like phase; thus, we define this phase as s ER o . Though this phase was recognized by the k-means clustering as an independent phase, at the moment is it not possible to confirm the direct colocalization of these sites, for the same reason as above. Finally, class 8 (Q1-Q1-Q3) represents hyper permeable, hyper-rigid, low viscosity domains. This class is localized in LD (Fig. 6), as shown in Fig. S1 (Supplement 1). Since these quantities are characteristic of the gel-like phase [2,3], we call this phase s LD o . Glucose treatment causes an increase in the extent of phase 5 (80%) at the expenses of phases 6 and 7 (Fig. 5(A)

Spatial distribution and organelle area fraction of PC-12 cells
In Fig. 6, AI images for control PC-12 cells and a comparison with the spatial distribution of the major intracellular organelles are reported. The organelles are labeled with organelle-GFP constructs as indicated in Materials and Methods. NGF-treated PC-12 cells present a polarized, elongated shape along a major axis.
l PM d and l PM o are localized on PM, covering ∼ (11 ± 3)% of the total cell area (Fig. 5(B)). The area of the PM is about (15 ± 3)% consistent (within 2σ) with the value retrieved for the phase area fraction (Fig. 7).
The phases l NET d and l NET o and s ER o cover ∼ (76 ± 12)% of the area fraction ( Fig. 5(B)). This value is consistent with the area fraction occupied by Golgi, ER and Mitochondria and other organelles (73 ± 10)% (Fig. 7). The spatial distribution of the phases is constituted by less connected and more connected, scattered domains, as one of the organelles. s LD o phase presents the same spatial distribution of LD labeled with Nile Red. This spatial distribution is characterized by a limited number of perinuclear connected regions consisting of LD clusters which are localized along the major axis of the cell. The average organelle area fraction of LD (12 ± 3)% (Fig.  7) is consistent (within 2σ) with the value retrieved for the phase area fraction (7 ± 3)% ( Fig.  5(B)). This value slightly deviates systemically because s LD o labels only the core of LD, while the external layer consists mainly of phase s ER o .

Discussion
In this work, we introduced an artificial intelligence clustering tool that leverages the multiparametric content of spectrally resolved lifetime images, resulting in an extension of the spectrum of detectable membrane phases. This method allows to reveal the phase organization inside PC-12 cell, dissecting and characterizing up to six different phases, with pixel resolution, and to study in real-time phase separations induced by different treatments. We found that PM of PC-12 cell is divided into two principal phases: a liquid-disordered phase l PM d (class 3) and a hyper-fluid, hyper-viscous, liquid-disordered phase l PM o (class 4). The developed method enabled the imaging of liquid-ordered domains of l PM o , which are characterized by a higher cholesterol concentration [32,33]. The liquid-ordered phase induced by the presence of cholesterol in the membranes has intermediate physical properties between s o and l d . While it increases the fluidity in high packed membranes, as the ones composed by saturated fatty acid tails, cholesterol brings order in highly fluid domains [38]. The imaging of liquid-ordered domains, which are characterized by a higher cholesterol content, though possible with a visual inspection of the traditional analysis and comparison of the several parametric images (Fig. 2(A)), cannot distinguish with enough specificity these two different phases in which a deep restructuration of lipid domains induced by cholesterol is present. We point out that the phase scenario of a PM is more complex and a more specific clustering, especially with the aid of Total Internal Reflection Fluorescence microscopy (TIRF) or Supercritical Angle Fluorescence (SAF) techniques, can in principle be achieved with the presented analysis method. The inner part of the cell is divided into four different phases: the first is a liquid-disordered phase l NET d (class 5) which is characteristic of intracellular connected regions with a networked structure, with the same fluidity of the liquid-disordered phase [2,3], but with a higher viscosity, which can be functional to the three-dimensional organization of these structures [38]. The spatial organization and the area fraction covered by this phase indicates that it is peculiar to some portions of ER, Golgi and the mitochondrial network membranes. This finding is supported by recent research, showing that ER and mitochondria cannot be considered as individual organelles in the cell, as they are structurally and functionally linked through contact points defined as mitochondria-associated membranes (MAMs) [39][40][41][42][43]. The first described function of ER-mitochondria contact sites is the exchange of phospholipids between organelles [44], thus facilitating their phase homogenization. This interconnected family of organelles can undergo a transition towards a liquid-ordered phase l NET o (class 6). Differently from l PM o , this phase is less hydrated thus indicating a higher packing. Its hyper-viscosity indicates the presence of cholesterol or a peculiar membrane protein organization, which can eventually favor the formation of more scattered vesicular membranes, like endosomes and synaptic-like vesicles. These isolated structures can account for the observed fragmented distribution of this phase. Then, the unsupervised learning of hyperspectral information detected a phase at ER sites surrounding LD (phase s LD o , class 8). This phase may be characteristic of LD formation sites (phase s ER o , class 7), that are specialized ER sites in which LD formation occurs, whose presence is supported by evidence in both yeast and mammalian cells [45]. Proteins involved in triacylglycerol synthesis are enriched at these sites within the ER, suggesting that neutral lipids could be synthesized at these discrete zones and that these platforms can help the clustering of specific domains allowing the LD budding-off. The sites retrieved by this method reveal that their phase is highly non-polar (phase s ER o ), thus enriched in triglycerides [20]. The distribution and the fraction of LD formation sites can, therefore, be revealed, providing insights in alterations of the activation of lipid storage and lipolysis pathways like the one we found, following glucose overload, on PC-12 cells, a well-established model to monitor the effects of diabetic neuropathy. Diabetic neuropathy is the most common microvascular complication of diabetes [46], and hyperglycemia has a relevant role in its development [47], inducing oxidative stress, apoptosis, and dysfunctions in neurons [48], glycation [49], the lack of neurotrophins and proinflammatory processes [50,51]. We observed that glucose treatment triggers a main transition l NET o → l NET d and a minor transition s ER o → l NET d Z. The main transition indicates a decrease of the less connected, vesicular phase inside cells. Since these vesicles originate mostly from the ER, this transition could be an early indication of the reduction of vesicle biogenesis rate. The alteration of these membrane phases leads to a visible rearrangement of the inner organelles, both in shape and dimensions, which can play an important role in affecting calcium homeostasis and signaling and in altering the pattern of neurotransmitter release [52,53], by impairing vesicle fusion and impedance of the membrane. This transition can also be linked to the high glucose level-induced disruption of MAM integrity and function through the activation of the pentose phosphate-PP2A pathway [54]. Indeed, MAM plays a crucial role in organelles biogenesis, being the site of autophagosome formation [55]. The second transition reveals a reduction of the LD formation sites, indicating the trigger of a metabolic switch turning off LD formation, though the number of LD stays nearly constant. The biological interpretation of these transitions, though needing further morphological and functional investigation, underlines the role of membrane phase remodeling in organelles biogenesis, and cannot be provided by contextual examination of the single parameter images (Fig. 2). A limitation of this technique is the requirement of a complex experimental setup (time resolved spectral detectors, two -photon laser). Moreover, because this technique requires the exploitation of a wide portion of the visible spectrum, it is difficult to perform co-labeling experiments and separate fluorescence contributes. A setup exploiting detection of far-red/infrared probes could in principle overcome this limit. Also spatial resolution can constitute a limit, especially when dealing with phase boundaries. With respect to TIRF and SAF microscopy for studying membrane dynamics [56], this technique lacks the selectivity for studying plasma membrane dynamics and phase transition. An advantage, however, is that it provides the possibility to analyze and dissect phase transitions for inner membranes.
The artificial intelligence clustering of multiparametric fluorescent images, by extending the spectrum of detectable cell membrane phases with the sub-micrometric resolution, opens the possibility to deeply investigate organelles biogenesis, a central focus of cell biology, and enables contextual insights on membrane remodeling and lipid metabolism. The real-time visualization of the biophysical alterations, by fully exploiting the information content of the photons emitted by solvatochromic probes, provides for a deeper understanding of the mechanisms linking the supramolecular emergent properties to cell functions.