Exploring heterogeneous cell population dynamics in different microenvironments by novel analytical strategy based on images | npj Systems Biology and Applications

npj Systems Biology and Applications volume 10, Article number: 129 (2024) Cite this article

Metrics details

Understanding the dynamic states and transitions of heterogeneous cell populations is crucial for addressing fundamental biological questions. High-content imaging provides rich datasets, but it remains increasingly difficult to integrate and annotate high-dimensional and time-resolved datasets to profile heterogeneous cell population dynamics in different microenvironments. Using hepatic stellate cells (HSCs) LX-2 as model, we proposed a novel analytical strategy for image-based integration and annotation to profile dynamics of heterogeneous cell populations in 2D/3D microenvironments. High-dimensional features were extracted from extensive image datasets, and cellular states were identified based on feature profiles. Time-series clustering revealed distinct temporal patterns of cell shape and actin cytoskeleton reorganization. We found LX-2 showed more complex membrane dynamics and contractile systems with an M-shaped actin compactness trend in 3D culture, while they displayed rapid spreading in early 2D culture. This image-based integration and annotation strategy enhances our understanding of HSCs heterogeneity and dynamics in complex extracellular microenvironments.

Cellular heterogeneity is a fundamental characteristic of biological systems, reflecting the diverse functional states and behavior of cells even within a seemingly homogeneous population1. Understanding these diverse states and their dynamic transitions is crucial for assessing population phenotype and addressing fundamental questions in biology, such as how cells respond to microenvironmental changes. Despite its importance, the temporal dynamics of heterogeneous cell populations in response to different microenvironments remain poorly characterized.

Current advancements in single-cell technologies and high-content imaging have enabled the acquisition of large-scale datasets, offering insights into cellular heterogeneity2,3. However, most conventional approaches, such as single-cell sequencing, primarily capture molecular snapshots of cellular states at specific time points and are unable to track the spatial-temporal distribution of molecules within the cell and the dynamic aspects of cell behavior2. These methods often disrupt cell integrity and lack time-resolved longitudinal studies on the same cells. While high-content imaging offers a non-invasive alternative tool to generate extensive imaging data with rich information4,5, integrating and interpreting these high-dimensional and time-resolved data to fully characterize the dynamic states of heterogeneous cell populations remains challenging6,7. High-content imaging provides detailed information including cell morphology, intracellular structures, and molecular distributions, representing the highly interactive molecular networks at the genomic, transcriptomics, and proteomic levels8,9,10, which is essential for properly governing cellular behavior and ultimately determining the cellular state responses to microenvironmental cues or other stimuli11,12. However, it remains difficult to integrate this complex information to interpret the dynamic phenotypes such as membrane dynamics and cell contractility of heterogeneous cell populations in different microenvironments. Machine learning techniques, both supervised and unsupervised, have been applied to process high-content imaging data for phenotypic screening13,14. While supervised machine learning offered substantial improvements in predicting cellular responses to various treatments, it required extensive labeled training datasets and might struggle with rare cell states not represented in the training set. Unsupervised machine learning has been applied to group subpopulations by protein location15 or drug effect16. Additionally, most of these image-based profile studies used static immunofluorescence images as input4, which overlook the temporal dynamics and the full spectrum of cellular behavior. Although existing studies have developed live-cell imaging systems to analyze biological processes, including the dynamic processes of epithelial-to-mesenchymal transition17, cell responses to drugs or ligands18,19, and cell trajectories in distinct cell-cycle phases20,21, capturing the dynamic states of heterogeneous cell populations in response to 3D collagen gel microenvironments remains underexplored6.

To address these questions, we proposed a novel analytical strategy to integrate and annotate these high-dimensional and time-resolved datasets to profile the dynamics of heterogeneous cell populations in different microenvironments. First, we acquired extensive image datasets using live-cell high-content imaging combined with F-actin fluorescent labeling. We then extracted high-dimensional features from these images, including cell shape, F-actin texture, and movement features. Next, we applied principal component analysis (PCA) and unsupervised clustering to systematically identify distinct cellular states based on their feature profiles and characterize their temporal dynamics. In addition, we identified distinct temporal patterns in cellular shape and texture features using unsupervised time-series clustering to provide insights into the dynamic processes underlying cell morphology changes and actin cytoskeleton reorganization.

In this study, we chose human hepatic stellate cells (HSCs) LX-2 as a model to track the cellular status change in response to different microenvironments. HSCs play a central role in liver fibrosis22, where they undergo phenotypic change in response to sustained liver injury23,24,25, making them a relevant and informative model for studying cellular dynamics in varying conditions. The activation of HSCs, which is characterized by morphological alterations, occurs within the context of a complex microenvironment rich in ECM proteins, particularly collagen, which provides chemical and physical signaling cues essential for cellular function26,27. Due to the essential role of HSCs activity in liver fibrosis pathogenesis, the dynamic behavior of HSCs including contracting and membrane protruding and the potential transition of activated HSCs phenotypes in response to microenvironment cues is important, however, remained incompletely understood. Here, we cultured HSCs LX-2 cells either on 2D substrates or within 3D matrices with various collagen densities. In 2D cell cultures, the LX-2 cells were grown on the flat glass surface, whereas in 3D cell culture, they were embedded in 3D type I collagen gels to highly simulate the fibrotic liver microenvironment, providing a platform for dissecting the effects of collagen gels with various densities on HSCs behavior28,29.

This novel high-dimensional and time-resolved image-based integration and annotation strategy enhances our comprehensive understanding of cellular heterogeneity and dynamic behavior of HSCs, particularly in how their morphology and actin cytoskeleton organization adapt to varying conditions. Finally, by combining advanced imaging and analytical techniques, we pave the way for future studies aimed at deciphering the complex phenotypes and behavior of cell populations, ultimately contributing to the development of more comprehensive research strategies and a better understanding of biological systems.

Our study aimed to elucidate the distinct cellular phenotypes and their change patterns in response to different microenvironments. We developed an image-based heterogeneous cell population dynamic profile analytical platform to capture the potential heterogeneity of phenotypes and dynamic transitions between different phenotypic states at different time points (Fig. 1). LX-2 cell line was used as a model, which was known to exhibit potential phenotypic plastic states in response to microenvironments cues30,31. LX-2 cells were cultured either on 2D substrates or within 3D collagen matrices with varying collagen concentrations (ranging from 1 to 9 mg/ml) for 48 h. We acquire extensive image datasets using live-cell high-content imaging combined with live-cell F-actin fluorescent labeling. High-dimensional features including shape, texture, and movement features of each cell were extracted from the large dataset of time-scale images.

Firstly, extensive image datasets were acquired using live-cell high-content imaging combined with F-actin fluorescent labeling, from which high-dimensional features were extracted. Next, distinct cellular states were identified by PCA and unsupervised clustering based on their feature profiles. In addition, distinct temporal patterns in cellular shape and texture features were identified using unsupervised time-series clustering to provide insights into the dynamic processes underlying LX-2 cell morphology changes and actin cytoskeleton reorganization. Elements of this figure was created with Microsoft Powerpoint using figures and plots generated in this manuscript.

Firstly, the single-cell datasets were aggregated as the mean value for each feature of LX-2 on each culture condition at each time point. Distinct dynamic patterns of typical features including cell roundness, compactness, cortical profile, and speed of cells cultured on 2D substrates or within 3D collagen gel (3 mg/ml) were shown (Fig. 2a–d). We found that cell roundness in 2D culture initially decreased within the first 24 hours before increasing, while in 3D culture, cell roundness showed a relatively stable trend and maintained smaller compared to 2D culture (Fig. 2a). Cell compactness in 2D culture exhibited a gradual increase over 48 hours, whereas cells in 3D culture displayed larger compactness compared to 2D culture, with a trend of initial increase followed by a slow decrease (Fig. 2b). Regarding cortical profile, cells in 2D culture initially showed an increase within the first few hours followed by a decrease, with fluctuations observed after 24 hours, whereas in 3D culture, cells exhibited smaller cortical profiles compared to 2D culture, with less variation (Fig. 2c). The speed of movement showed a similar change pattern as cortical profile in 2D and 3D cultured cells (Fig. 2d).

a–d The line represents mean values of typical features (roundness, compactness, cortical profile, and speed) extracted from cell populations on 2D surfaces or within 3D collagen gels (3 mg/ml) changes over time. Roundness was used to measure how closely the cell shape matches a perfect circle and was calculated by the formula Roundness = sqrt(4*pi*Area) / Perimeter. Compactness was used to describe how compact the actin filaments inside the cell. Cortical profile was used to characterize the location of the actin intensity in cortical regions. Speed was measured by the equation Speed = Step size / Time intervals. The shaded regions above and below each line show ± SEM. e PCA visualization of LX-2 cells culture on 2D surfaces and within 3D collagen gels (n = 180,000). Clusters were colored according to the three clusters identified by k-means clustering. f The proportion of distinct clusters on 2D surfaces vs. in 3D collagen gels (3 mg/ml). g The proportion of cell states changes over time on 2D surfaces. h The proportion of cell states changes over time within 3D collagen gel (3 mg/ml).

After calculating the mean values for these features over time, we observed diverse trends in their temporal changes as mentioned above. To further comprehensively characterize different cellular states, we employed the k-means unsupervised clustering on the high-dimensional feature data (n = 180,000). The features selected for clustering analysis are listed in Supplementary Table 2. We categorized cell objects into distinct groups based on the similarity between the single-cell high-dimensional feature profiles. We identified three clusters (Cluster 0, Cluster 1, and Cluster 2), representing the three cell states mentioned below.

Comparing the relative proportions of each cluster within the LX-2 cell population under different conditions, we found that Cluster 2 accounted for the highest proportion of the cell population in both 2D and 3D cultures, represented the most common cell state, while the other two clusters accounted for a relatively smaller proportion (Fig. 2e–f). Additionally, we found that the relative proportion of cell states changed over time, even under the same conditions (Fig. 2g–h). This indicated that cell states exhibit distinct levels of plasticity in both 2D and 3D cultures. Cells might transition between different states at varying rates and times, suggesting temporal dynamics and complex cellular behavior influenced by microenvironments cues. Therefore, to understand how cellular populations under different conditions exhibit various cell states and undergo cell state transition, we need to systematically and comprehensively characterize each cell state in detail.

Distinct cell states identified based on the shape, texture, and movement features of LX-2 represent the global manifestation of the structural organization of the cytoskeleton, reflecting the specific behavior of the cell in response to the microenvironment cues. To further elucidate the different phenotypic representations of cellular populations under 2D and 3D conditions, we characterize each cluster in detail based on the typical features selected above.

Specifically, we observed that cells within Cluster 0 exhibited higher roundness, while cells in Clusters 1 and Cluster 2 displayed relatively lower roundness, which meant irregular cell shapes (Fig. 3a). Regarding the compactness of intracellular actin filaments, cells in Cluster 1 displayed higher intracellular actin density, while cells in Clusters 0 and Cluster 2 displayed relatively lower actin density (Fig. 3a). The actin filaments located at the cortical region within Cluster 2 were higher than the other two clusters. The instantaneous speeds of LX-2 within Cluster 2 were a little bit faster than those in other clusters. In brief, Cluster 0 was characterized by a rounded shape and low actin density, Cluster 1 exhibited irregular shapes with high actin density, while Cluster 2 displayed irregular shapes with more actin located at the cortical region. By elucidating these distinct cellular features, our study aimed to shed light on the subtle phenotypic diversity and functional implications under different cellular states, thereby contributing to a more comprehensive understanding of HSCs dynamics.

a Violin plots of roundness, compactness, cortical profile, and speed across 3 cell clusters. Mann–Whitney-Wilcoxon test two-sided with Bonferroni correction (*: 0.01 < p < = 0.05, **: 0.001 < p < = 0.01, ***: 0.0001 < p < = 0.001, ****: p < = 0.0001). In all violin plots, central white dots represent the median values, black boxes delineate the interquartile range, and thin black lines indicate a maximum of 1.5 times the interquartile range beyond the box. Each side of the black lines consists of a kernel density estimation curve to show the distribution of the feature. b Proportion of distinct cell states in 3D matrices with varying collagen concentrations (ranging from 1 to 9 mg/ml).

We observed that the proportions of relatively rare cell states with rounded shapes and low actin density (Cluster 0, 23.9%) and cells with irregular shapes and high actin density (Cluster 1, 8.4%) in 3D culture (3 mg/ml) were relatively higher than those in 2D cultures, while the proportion of the relatively common cell states with irregular shapes and high cortical actin intensity (Cluster 2, 67.7%) in 3D culture was smaller compared to those in 2D culture (Fig. 2f). This suggested that in 3D culture, the heterogeneity in cellular phenotypes was more pronounced, as observed by the increase in the proportions of cell states with low probability compared to 2D culture.

Next, we further compared the effects of different collagen concentrations (from 1 to 9 mg/ml) on cell state transitions in 3D culture. We found that the proportions of cells with irregular shapes and high actin density (Cluster 1) exhibited a unimodal trend with increasing collagen concentration, peaking at 3 mg/ml (Fig. 3b). And the proportions of cells with rounded shapes and low actin density (Cluster 0) were highest at 9 mg/ml (Fig. 3b). These findings suggest that varying collagen concentrations influence the distribution of cellular states in 3D culture, with distinct effects on different cell states.

The temporal changes in cell state proportions described above illustrated that cells might undergo transitions between different cellular states over time even under the same conditions. These cellular states are defined based on a comprehensive combination of cell shape, texture, and movement features. To further elucidate the independent contributions of these features and their temporal variation patterns to cell phenotypes, we applied k-means unsupervised time-series clustering methods to classify the dynamic patterns of individual features. Firstly, we employed two key shape features, roundness, and area, for time-series clustering analysis because cell roundness and area change patterns may indicate membrane deformation32,33 and cell behavior including cell spreading, protruding, and membrane trafficking.

For roundness, four trend clusters were identified (Roundness-trend-cluster 1 to 4). Cell roundness first decreased and then remained constant in Roundness-trend-cluster 2, while it remained stable in the early stage and then increased after about 36 hours in Roundness-trend-cluster 3. The roundness in the other two patterns showed an unsynchronized increase and then a tendency toward decrease (Fig. 4a). The proportion of cell states within each roundness trend cluster changes over time (Fig. 4b). It showed a similarity in the roundness patterns and the trend of proportion for cell state of Cluster 0 (Fig. 4b), which was identified as cell states with rounded shapes and low actin density above. We further compared the proportions of four roundness patterns of LX-2 on the 2D surface with those in 3D cultures with different collagen concentrations (ranging from 1 to 9 mg/ml). We found that the proportion of cells changed from round to irregular shape (Roundness-trend-cluster 2) in 3D culture was higher compared to that in 2D culture, with the highest proportion found at 9 mg/ml 3D collagen culture (Fig. 4c). The proportion of cell roundness remained stable in the early stage and then increased (Roundness-trend-cluster 3) in 3D culture was lower than that in 2D culture, with relatively higher proportions observed at 3 mg/ml 3D collagen culture (Fig. 4c).

a, d Identification of 4 roundness or area trend clusters, respectively. The k-means time-series clustering method is proposed for trend clustering based on 2668 time-series data. b, e The proportion of cell states at each time point for each roundness or area change pattern, respectively. The color legends (C0, C1, and C2) represented the three cell state clusters identified above in Fig. 2e. c, f The proportion of roundness or area trend clusters in 2D cultures and 3D cultures with different collagen concentrations (ranging from 1 to 9 mg/ml), respectively.

Additionally, four area trend clusters were identified as well (Area-trend-cluster 1 to 4). We found that the cell area remained stable in the initial stage and then spread in Area-trend-cluster 1 and 4, while cells showed a fast spread and then shrunk asynchronously in Area-trend-cluster 2 and 3 (Fig. 4d). The area trend was similar to the trend of proportion in cell state of Cluster 2, which referred to irregular cell morphology and high cortical actin intensity (Fig. 4e). This suggests that the main driving force during cell spreading is the polymerization of actin filaments that push the cell membrane forward. We further revealed that Area-trend-cluster 2 and 3 accounted for the main proportion of cells on the 2D surface (Fig. 4f), which suggests a fast spread mode in the initial stage of the LX-2 cultured on 2D substrates. However, LX-2 in 3D culture with different collagen concentrations showed diverse patterns of spreading and shrinking.

To further explore the relationship between roundness and area trends, we identified a Spearman correlation of -0.72 in 3D cultures and -0.66 in 2D cultures. These negative correlations indicate that as cell area increased, cells tended to become less round (i.e., more elongated or irregular in shape), and as cell area decreased, they tended to become more rounded. We performed an overlap analysis between roundness trend clusters and area trend clusters and found more overlapped cells between clusters where roundness and area exhibited opposing trends (Supplementary Table 1). This further demonstrates that most cells exhibited a negative correlation between changes in roundness and area.

Distinct temporal patterns of cellular shape features described above showed the different shape, spread and contraction patterns of cells in response to 2D and 3D microenvironments. In addition to cell shape, actin cytoskeleton dynamics play a central role in driving cellular state changes and regulating multiple cellular behavior in response to microenvironmental cues34,35. Here, we applied k-means time-series clustering analysis on the compactness and cortical profile of actin filaments, to access the dynamic pattern of the distribution and location of cytoskeletal actin filaments.

Four compactness trend clusters were identified (Compactness-trend-cluster 1 to 4). The compactness showed a bimodal trend in Compactness-trend-cluster 1, a unimodal trend in Compactness-trend-cluster 4, a slow upward trend in Compactness-trend-cluster 2, and a slow downward trend in Compactness-trend-cluster 3 (Fig. 5a). The compactness trend was consistent with the trend of proportion in cell state of Cluster 1, which referred to irregular cell shapes and high actin density (Fig. 5b). We also found that Compactness-trend-cluster 2 accounted for the highest proportion of cells on 2D culture (Fig. 5c), indicating that most 2D-cultured cells showed a slow increase in actin density. The proportion of compactness showed a bimodal trend (Compactness-trend-cluster 1) in 3D culture was higher than that in 2D culture (Fig. 5c). This suggests a higher active actin cytoskeleton of LX-2 in the 3D collagen microenvironments.

a, d Identification of 4 compactness or cortical profile trend clusters, respectively. The k-means time-series clustering method is proposed for trend clustering based on 2668 time-series data. b, e The proportion of cell states at each time point for each compactness or cortical profile change pattern, respectively. The color legends (C0, C1, and C2) represented the three cell state clusters identified above in Fig. 2e. c, f The proportion of compactness or cortical profile trend clusters in 2D cultures and 3D cultures with different collagen concentrations (ranging from 1 to 9 mg/ml), respectively.

As for cortical profile trend clustering (Cortical-profile-trend-cluster 1 to 4), the cortical actin intensity decreased slowly and then increased asynchronously in Cortical-profile-trend-cluster 1 and 2, increased slowly and then decreased in Cortical-profile-trend-cluster 3, and increased and then remained stable in Cortical-profile-trend-cluster 4. (Fig. 5d). We found that the Cortical-profile-trend-cluster 4 was consistent with the trend of proportion in cell state of Cluster 2, which referred to irregular cell morphology and high cortical actin intensity (Fig. 5e). The proportion of cortical profile showed an increase trend (trend cluster 4) in 2D culture and 3D culture with low density of collagen was higher than that in 3D culture with high density of collagen (Fig. 5f), which indicated an active actin cortex of the LX-2 cultured on the 2D substrates and in the low-density 3D collagen gel.

We developed an image-based heterogeneous cell population dynamic profile integration and annotation strategy to elucidate the distinct cellular states and the dynamic phenotypic transitions of HSCs in response to different microenvironments (Fig. 6). Our study highlighted the importance of temporal dynamic analysis techniques in capturing the full spectrum of cell population phenotypes.

Through comprehensive analysis of the shape, texture, and movement features of HSCs LX-2, we identified three cell states in both 2D and 3D cultures, including round shape and low actin density, irregular shape and high actin density, and irregular shape and high cortical actin density. The identification of dynamic patterns of specific features related to cell shape and cytoskeletal organization showed that cells exhibited higher complexity in their membrane dynamics and contractile systems in 3D microenvironments compared to 2D microenvironments. Elements of this figure was created using Microsoft Powerpoint.

Through comprehensive analysis of the shape, texture, and movement features of HSCs LX-2, we identified three cell states in both 2D and 3D cultures, including round shape and low actin density, irregular shape and high actin density, and irregular shape and high cortical actin density. The study revealed that cells exhibit heterogeneous responses to different microenvironments, with varying proportions of cell states and dynamic transitions over time.

We first observed differences in dynamic patterns of various cellular features, including cell roundness, compactness, cortical profile, and speed, between 2D and 3D cultures which highlighted the influence of microenvironmental cues on HSC behavior (Fig. 2a–d). 3D culture introduced a more physiologically relevant tissue-like architecture, including spatial constraints and cell-ECM interactions. In this context, the observed changes in cell morphology may indicate adaptations to the 3D microenvironment, including cell spreading and cytoskeleton organization within the 3D collagen matrix. LX-2 in 3D culture exhibited more irregular morphology and larger compactness compared to 2D culture might reflect differences in the activation state of LX-2 cells (Fig. 2a). The altered morphology of LX-2 cells in 3D culture could be influenced by mechanical forces exerted by the surrounding ECM36. Cells might respond to these forces by adopting a more elongated or compact cytoskeleton morphology (Fig. 2b), which could have implications for their functional behavior, including ECM remodeling and fibrogenesis. Additionally, the temporal changes of cortical profile and speed (Fig. 2c–d), showed similar patterns in both 2D and 3D cultures, suggesting correlative regulatory mechanisms underlying these dynamic processes37.

While traditional statistical methods based on the average level of individual features at various time points above provided valuable information, they might tend to focus on overall trends and could overlook the potential presence of diverse cell phenotypes, especially rare ones, which might exhibit spatially and temporally heterogeneous responses to external stimuli. To address these limitations, we employed unsupervised clustering to categorize cells into distinct groups based on their feature profiles (Fig. 2e). We identified three clusters representing different cell states. The first cell state (Cluster 0) with a round shape and low actin density was often associated with a less activated state of HSCs38 reduced capacity for mechanical support or a high extracellular space constraint. It also suggested cells with more stable and less dynamic actin cytoskeleton, which might indicate cells in a resting state or one that was not actively engaged in activated processes. The second cell state (Cluster 1) with an irregular shape and high actin density was associated with a stable cytoskeleton with mechanical support and force generation capacity, indicating a typically more plastic phenotype39. The third cell state (Cluster 2), with an irregular shape and high cortical actin density, indicating a high degree of plasticity and the ability to interact dynamically with the surrounding microenvironments, counted for the highest proportion relative to the others. Actin filaments were concentrated in the cortical region, potentially corresponding to processes for force generation, including the formation of stress fibers40, lamellipodia41, or filopodia42, which were essential for cell adhesion and interaction with the extracellular matrix. The observed differences in cell state proportions between 2D and 3D cultures shed light on the complex interplay between cell morphology and microenvironmental cues. The increase in heterogeneity in cellular phenotypes in 3D culture, as shown by the higher proportion of those relatively rare cell states occurring less frequently (Fig. 2f), underscored the importance of considering the dynamic transitions in cellular states that represent the spatial organization and cell-matrix interactions in different contexts.

The specific features including cell area, roundness, actin compactness, and the cortical actin profile are well-established indicators of membrane deformation32,33 and cell contractility43,44. The interaction of cells with the surrounding microenvironment leads to trigger membrane deformation that allows cell spreading33, as well as cytoskeletal rearrangement and cell contraction that drives cell migration43. These activities are closely related to the cell ability to navigate through complex 3D matrices45 and the broader signaling pathways that govern cellular mechanosensation and response to extracellular cues46. Then we applied k-means time-series clustering methods to classify the dynamic patterns of these specific features to further elucidate their diverse temporal patterns and their independent contributions to cellular activities. We found a higher proportion of cells cultured in 3D with a tendency to change from round to irregular shape (Fig. 4c), indicating higher membrane dynamics in 3D collagen microenvironment. The rapid spread pattern of LX-2 in the early stage of culture in the 2D microenvironments illustrated a different spreading mode compared with that in the 3D microenvironments. The negative correlation between roundness and area trends was particularly significant in 3D cultures, where cells adapted to a more mechanically complex environment. The increase in cell area accompanied by decreased roundness likely reflects the need for greater membrane flexibility and contractile properties, facilitating cellular navigation and interaction with the mechanically complex extracellular matrix. These morphological adaptations reflects cytoskeletal reorganization and mechanical responsiveness to environmental stiffness and matrix density, which are critical in HSC activation and fibrosis35,47. Temporal trends in compactness and cortical profile of actin filaments correlated with shifts in the proportions of specific cell states, suggesting a link between cytoskeletal dynamics and cellular phenotypic transitions. The higher proportion of an M-shaped compactness pattern with two distinct peaks in the 3D collagen microenvironments indicated the enhanced actin cytoskeletal reorganization and cell contractility in this context. LX-2 cells cultured on 2D substrates and in low-density 3D collagen gel exhibited a more active actin cortex, suggesting enhanced cell adhesion activity. It indicates the sensitivity of HSCs to spatial constraints and physically complex extracellular matrix47, which requires a higher degree of membrane flexibility and contractile force generation to enable spreading and migration.

Our study proposed a novel analytical strategy designed for high-dimensional, time-resolved analysis of cellular morphology in both 2D and 3D microenvironments. Using the novel analytical strategy, we revealed the heterogeneous cell states and dynamic activities that were not well-characterized before. Our study demonstrates new insights into cell behavior, particularly by identifying complex membrane dynamics and contractile systems in LX-2 cultured in 3D microenvironments. These findings suggest that 3D collagen microenvironment influenced membrane flexibility and cell contractility of HSCs, which are critical for cellular processes such as spreading and migration, both of which play key roles in HSCs activation and their fibrogenic response26,35. This can help elucidate how cells regulate their mechanical properties in response to external cues and how alterations in cell-matrix interactions affect cell membrane deformation and cell contractility. Our study focused on cell morphology and dynamic behavior, which are believed to provide a foundational perspective that bridged morphological changes with regulatory systems8,9. It has been reported that the complex interactions between cells and their microenvironment, particularly in 3D collagen matrices, influence intracellular signaling pathways that regulate membrane tension and cytoskeletal contractility46,48. For example, previous studies have shown that mechanical cues from the ECM are transduced via integrin-mediated focal adhesions and activate RhoA/ROCK signaling pathways, driving changes in actin polymerization and cell contractility49,50,51. These pathways are key to understanding how HSCs respond to ECM stiffness during fibrosis progression, providing a functional link between the morphological changes we observe and underlying regulatory networks52,53. Our findings provided a critical basis for such mechanistic studies, offering new insights and hypotheses that can be further investigated in the context of systems biology. Our study provides a valuable platform and resource that can be integrated with molecular mechanistic studies in future research to elucidate the potential system biological mechanisms involved.

In conclusion, our study provides a novel analytical strategy for image-based integration and annotation to profile the dynamics of heterogeneous HSCs cell populations in different microenvironments, highlighting their dynamic nature and sensitivity to microenvironmental cues. This high-dimensional and time-resolved image-based integration and annotation solution is able to be broadly applied to other cellular systems and conditions, uncovering additional critical insights into cellular dynamics and environmental responses. Future studies should focus on exploring the molecular mechanisms underlying these phenotypic transitions and their implications for treatment.

Human HSCs cell line LX-254 was purchased from Procell Life Science & Technology Co., Ltd. (CL-0560, China). In 2D cell cultures, LX-2 cells were grown on the flat glass surface of the CellCarrier Ultra 96-well microplate (Revvity, USA) and maintained in Dulbecco’s Modified Eagle medium (DMEM, Gibco, USA), supplemented with 2% (vol/vol) fetal bovine serum (FBS, Gibco, USA) and 1% penicillin/streptomycin (Gibco, USA).

In 3D cell culture, LX-2 cells were embedded in gels of controlled type I collagen density, as described previously55. Briefly, cells suspended in the culture medium were mixed with an appropriate volume of rat-tail collagen I (354249, Corning, USA) in a 1:4 (vol/vol) ratio to obtain the desired collagen concentration (1 to 9 mg/ml). A calculated amount of 1 M NaOH was quickly added to mix the final solution well to achieve a pH of ∼7. The mixed suspension was added to the microplate mentioned above and immediately placed in an incubator maintained at 37 °C for 40 minutes to allow polymerization. Fresh culture medium as mentioned above was added to maintain cell growth. Cells in 2D and 3D cultures were all maintained in a humidified incubator at 37 °C and 5% (vol/vol) carbon dioxide.

CellMask Deep Red Actin Stain (Invitrogen, USA) was used to label the F-actin as per the manufacturer’s instructions. Cells were incubated in culture medium with 0.01% CellMask Deep Red Actin Stain and Hoechst 33342 (Thermo Fisher Scientific, USA) for at least 30 minutes at 37 °C and 5% CO2.

Live-cell high-content imaging was performed using Opera Phenix High-Content Screening System (Revvity) at 5% CO2 and 37 °C temperature. Twenty-five fields of view were imaged for each well, with 20 z-stacks per field at 5 μm intervals using two channels including Hoechst (excitation, 405 nm; emission, 435–480 nm) and CellMask Deep Red (excitation, 640 nm; emission, 650–760 nm). Imaging was performed every 90 min for 48 hours with a 20× water immersion lens in confocal mode to track cellular phenotypic dynamics using 2×2 pixel binning. A total of 33000 images were captured for each well.

An analysis sequence was designed in Harmony Software (Revvity, v5.1). Firstly, flatfield correction was performed on all images and maximum intensity projection was made from each plane of images at different focus depths. Next, images were segmented to identify the cell nucleus and cytoplasm. Cells with low intensity in any of the channels or located at the edge of the image were excluded. Multiple phenotypic parameters were then quantitatively profiled for up to five regions, including cell, nucleus, cytoplasm, membrane and ring region. These parameters included classical features such as intensity or morphology and advanced descriptors such as STAR morphology and SER texture, which provided statistically robust measurement of F-actin distribution patterns within cells. The object tracking process was then performed on the time series measurement. The time-dependent properties such as current speed or object displacement were then calculated. Overall, approximately 500 cells were analyzed per group at a single time point, and a total of 1327 features (shape, texture, movement, etc.) were measured for individual cells at each time point.

The high-dimensional and time-resolved dataset was analyzed and plotted in Python. Firstly, normalization was performed by median absolute deviation (MAD) normalization. A feature selection procedure was then performed to eliminate features with low variance, missing values, and outliers56, and 103 features remained for further analysis (Supplementary Table 2). In detail, we applied a variance thresholding approach to remove features with low variance, which could introduce noise or provide little information for downstream clustering analyses. Features were excluded if they had a frequency threshold below 0.05 (the ratio of the second most common value to the most common value was less than 0.05) or a uniqueness threshold below 0.01 (the number of distinct values was less than 1% of the total number of samples). PCA was performed to extract components accounting for at least 95% variability of the dataset. Distinct clusters were then identified by applying the K-means clustering method (k = 3, n = 180,000). The optimal number of clusters k was determined by the maximum of the average silhouette. As for single-cell time-series profiling, cells that were not tracked successfully for every time point over 48 hours were excluded and 2668 single-cell time-series profiles remained for subsequent PCA. Distinct trend-clusters of specific cell shape and texture features were obtained by applying K-means time-series clustering (k = 4) to obtain dynamic patterns.

Experiments were performed at least twice. Statistical analysis was performed with Python packages, including SciPy57 and Scikit-learn58. The details of statistical tests can be found in the legend of each figure. We represent significance in the following way with asterisks: * (0.01 < P ≤ 0.05), ** (0.001 < P ≤ 0.01), *** (0.0001 < P ≤ 0.001), and **** (P ≤ 0.0001).

The datasets used in this study are available in Figshare, https://figshare.com/s/8ae0cfb6cbfbc306cce1.

The underlying code for this study is available in Figshare, https://figshare.com/s/8ae0cfb6cbfbc306cce1.

Altschuler, S. J. & Wu, L. F. Cellular heterogeneity: do differences make a difference? Cell 141, 559–563 (2010).