First, well-conditioned HeLa cells were spread in six-well plates and cultured in adherence at an inoculum density of 1 × 10⁶ cells/mL per well for 24 h. Then, cells were stimulated with Diniconazole (DIN) at concentrations of 0, 2.5, 5, 10, and 20 μmol/L for 24 h. The culture environment was maintained at 37 °C and 5% CO₂. After stimulation, bright-field imaging was performed using a 10× objective of an optical microscope (OLYMPUS CKX53), and conditions such as exposure time and light intensity were kept consistent to ensure uniform imaging conditions for all sample images. One hundred images were acquired for each concentration, and the resolution of the images was 3072 × 2048 pixels, totaling 500 samples. For HUVEC, and HaCaT cells, stimulation was performed using DIN at concentrations of 0, 2.5, 5, 10, and 20 μmol/L, and a total of 150 samples were collected. For A549 cells, stimulation was performed using DIN at concentrations of 0, 12.5, 25, and 50 μmol/L, and a total of 120 samples were collected. For PC12 cells, stimulation was performed using DIN at concentrations of 0, 5, 10, 20 and 40 μmol/L, and a total of 150 samples were collected. It should be noted that this study adopted a phased image acquisition strategy. Initially, HeLa cells were used as the model system for method development. To ensure statistical robustness and provide a solid data foundation for feature analysis, 100 images were collected at each concentration. After successful method validation, to assess its cross-cell line generalization ability and feasibility in scenarios with limited resources, 30 images were collected per concentration in the subsequent experiments with A549, HUVEC, PC12, and HaCaT cells.

In order to further analyze the grayscale characteristics of the cellular regions, we performed batch automatic segmentation of the cellular images by calling Cellpose’s API and exported the corresponding masks (Stringer et al., 2021). Cell segmentation was performed using the Cellpose pre-trained model (cyto3). The parameters were set as follows: automatic diameter estimation, grayscale image channels, a flow threshold of 0.4, and a minimum cell size of 15 pixels. Tiled processing and network averaging were enabled to balance the processing of large images with segmentation accuracy. All segmentation tasks were accelerated using GPU processing. To visually demonstrate the image segmentation results, we overlaid the mask generated by Cellpose onto the original image during the image processing stage in Figure 1. This clearly shows the extracted cell regions (as illustrated by the grayscale example in the image processing module), thereby providing an intuitive basis for quality verification in subsequent grayscale feature extraction. Then, we batch-calculated the grayscale values of cell regions and their corresponding pixel numbers based on mask, and adjusted the data with normalization based on the average occupancy of all samples to ensure comparability among different images. Ultimately, we constructed five datasets containing the cell grayscale distribution and DIN concentration for subsequent analysis and modeling.

We performed a visual analysis of the samples at each concentration, using the sample stimulated with 10 μmol/L DIN in HeLa cells as an example. Figure 2 shows the grayscale distribution of 100 representative samples, where a small number of outliers were observed. To address this, we applied the interquartile range (IQR) method to remove outliers (Swami et al., 2023). When the scaling factor k = 3, it effectively eliminated apparent outlier curves while retaining most of the samples, thereby enhancing data reliability and preventing outliers from negatively affecting the model. Observing the gray scale distribution curve after removing the outliers in Figure 3, we found that the data still has a certain degree of dispersion. Therefore, we standardized the data to make the data obey the standard normal distribution, which in turn improves the training effect of the model (Muhammad and Faraj, 2014).

SR is a feature selection method based on statistical significance (P-value). It iteratively filters out the most relevant features to the target variable, thus reducing redundant features and improving the efficiency and predictive power of the model (Zhou and Jiang, 2016).

The SR-LR model is a combination of SR and LR, which is a simple and interpretable method of linear modeling through screened features (James et al., 2023). The SR-RR model introduces L2 regularization on this basis, i.e., it is a combination of SR and RR. RR reduces multicollinearity through L2 regularization, limits model weights, and reduces the risk of overfitting, thus improving the stability and predictive ability of the model (Wang X. et al., 2021). The advantage of RR is that it can better handle the situation of correlation between features and control the model complexity while maintaining the linear modeling framework.

In the SR-MLP model, we employ the MLP for nonlinear modeling. The MLP is a feedforward neural network consisting of multiple hidden layers, which captures complex nonlinear relationships in data through weighted summation and activation functions (Liu et al., 2022). The structure of the MLP—such as the number of hidden layers, the number of neurons, and the activation functions—can be flexibly adjusted according to the specific task requirements. In this study, we constructed an MLP network with three hidden layers, containing 72, 48, and 24 neurons, respectively. The hidden layers use the Rectified Linear Unit (ReLU) activation function to enhance the modeling capability for complex nonlinear relationships (Alkhouly et al., 2021; Fu et al., 2020), while the output layer employs a linear transformation (Linear) to directly output continuous predictions for the regression task (Kılıçarslan et al., 2021).

To prevent overfitting and enhance the model’s generalization capability, we employed the following regularization strategies during training: First, L2 regularization was introduced by adding a weight penalty term to the loss function, thereby constraining the magnitude of weights and suppressing overfitting (Xie et al., 2022). Second, the Dropout method was applied to randomly drop a portion of hidden layer neurons during training, reducing model complexity and further improving its generalization ability (Salehin and Kang, 2023). We monitored the trends of training loss and validation loss to guide hyperparameter optimization. Model convergence was considered satisfactory when both losses decreased smoothly, remained closely aligned, and eventually stabilized. A scenario where training loss continued to decline while validation loss plateaued or increased prematurely indicated potential overfitting, whereas persistently high losses with slow reduction suggested underfitting (Wang et al., 2022). Accordingly, we systematically adjusted hyperparameters including the number of training epochs, batch size, regularization strength, learning rate, and Dropout rate to achieve optimal convergence performance.

To fully evaluate the performance of the hybrid model, we use the following commonly used metrics:

Coefficient of determination (R²): a measure of the goodness of fit of the model, the closer the value is to 1 the better the model fits the data and the better the prediction (Chicco et al., 2021).

We evaluated the performance of the three hybrid models (SR-LR, SR-RR, and SR-MLP) and the baseline models (LR, RR, and MLP) using R², RMSE, and MAE metrics, as summarized in Table 1. The results indicate that all hybrid models incorporating SR outperformed their non-SR counterparts across most datasets, underscoring the efficacy of SR in feature selection for enhancing model performance. And among these three hybrid models, the R² values of SR-MLP on all cellular datasets were significantly better than those of SR-LR and SR-RR and higher than 0.9, indicating a better fitting ability. In the case of the HeLa cell dataset, for example, the R² of SR-MLP was as high as 0.9919, while SR-LR and SR-RR were 0.9707 and 0.9708, respectively. Moreover, the RMSE and MAE of SR-MLP are lower than those of SR-LR and SR-RR on all datasets, which indicates that its prediction error is smaller and its stability is higher. In contrast, the R² of SR-LR and SR-RR was higher than 0.9 on the HeLa, PC12 and HaCaT cell datasets, but lower than 0.9 on the A549 and HUVEC cell datasets, which indicated a poorer fit and limited generalization ability.

To further assess the predictive power of SR-MLP, we compared the predicted and true values of SR-RR (e.g., Table 2) with those of SR-MLP (e.g., Table 3) models. Taking the HeLa cell dataset as an example, the predicted values of SR-MLP and SR-RR were 19.44 μmol/L and 19.07 μmol/L, respectively, at a high concentration of 20 μmol/L, whereas the predicted values of SR-MLP and SR-RR were 2.65 μmol/L and 3.54 μmol/L, respectively, at a low concentration of 2.5 μmol/L. It can be seen that the predicted values of SR-MLP are closer to the true values, and the errors are smaller, indicating that the model has the ability to predict the true values at different concentrations. MLP’s prediction value is closer to the real value with smaller error, indicating that the model has strong prediction accuracy under different concentration conditions.

In order to reduce data dimensionality and enhance model stability, we employed the SR method to screen gray-scale features based on statistical significance. During the screening process, differentiated p-value inclusion thresholds (ranging from <0.001 to <0.05) were set according to the data characteristics of different cell datasets and the need to control model complexity. Specifically, for the HeLa dataset, which exhibited clearer feature-concentration relationships and higher data quality, a stricter threshold (p < 0.001) was applied to construct a highly concise initial feature set. For the other datasets, relatively lenient thresholds (p < 0.05 or p < 0.01) were adopted to retain more potentially relevant feature information for the subsequent MLP to perform nonlinear learning. The setting of all thresholds followed the statistical principle of controlling Type I errors while also considering the connection with the subsequent nonlinear modeling stage (Aguinis et al., 2021).

Despite the differences, the screening results also indicate some similarities among key features across different cell datasets, as shown in Table 4. For instance, feature 244 was selected in the HeLa, A549, and HUVEC cell datasets, while key feature 204 for PC12 cells and key feature 248 for HaCaT cells, along with feature 244, all fell within the gray-scale range of 200–255. This suggests that these features may possess strong generalizability across different cell types and can stably reflect the effects of pollutants on cells. The gray-scale range of 200–255 typically corresponds to abnormal cells, and features within this range exhibit more significant variations. Moreover, the features screened in this range demonstrated better p-values, further validating the importance of these features in the cellular damage response. Therefore, during the feature selection process, we focused on the key features of cell damage and controlled the number of features within 20 to optimize the computational efficiency and generalization capability of the model.

To elucidate the contribution of key features in the SR-MLP model to pollutant concentration prediction, we applied the SHapley Additive exPlanations (SHAP) method. Derived from game-theoretic Shapley values, SHAP provides a consistent and interpretable quantification of feature contributions. Using SHAP’s Explainer, we computed SHAP values for key features on the test set and visualized the results with beeswarm plots. In these plots, the X-axis represents SHAP values (indicating directional contribution), the Y-axis lists features in descending order of mean absolute SHAP value (importance), and point color reflects the original feature magnitude (red for high, blue for low), thus enabling an intuitive assessment of how feature values relate to their predictive influence (Ekanayake et al., 2022).

Analysis results are highly consistent with the key features selected based on statistical significance, strongly confirming the experimental finding that “features in the high grayscale range (200–255) are the core drivers for predicting pollutant concentration.” The SHAP importance ranking reveals that across all cell types, features located within this high-brightness interval dominate (The SHAP beeswarm plot for HeLa cells is shown in Figure 6, while the corresponding results for the remaining cell types (A549, HUVEC, PC12, HaCaT) are provided in the Supporting Information as Supplementary Figures S1–S4). Specifically, in HeLa cells, the top three features by importance are 242, 245, and 244; in A549 cells, they are 244, 221, and 218; in Huvec cells, they are 244, 198, and 48; in PC12 cells, they are 204, 76, and 200; and in HaCaT cells, they are 192, 37, and 172. It is particularly noteworthy that feature 244 ranks highly in HeLa, A549, and HUVEC cells, indicating its stable predictive ability across different cell types. The results for HaCaT cells show certain particularities: the top-ranked feature in SHAP importance is 192, which, although close to, does not fully fall within the 200–255 high-brightness range. This may suggest that under DIN stimulation, HaCaT cells exhibit differences in brightness response thresholds or grayscale distribution patterns related to morphological damage compared to other cell lines—for example, a smaller increase in brightness during apoptosis or interference from optical properties of other subcellular structures. This phenomenon warrants further investigation in future studies, considering cell type specificity.

In Figure 6, key high-grayscale features (such as 242, 245, and 244) exhibit a distribution where blue points are concentrated on the left side of the X-axis, indicating that their lower feature values are associated with negative SHAP values. This suggests that the model has learned not a simple linear positive correlation but a more complex nonlinear response relationship: it is likely that only when these feature values exceed a certain threshold do they contribute positively and clearly to the prediction of high concentrations. This precise capture of nonlinear dependencies and feature interactions is precisely the intrinsic key to the SR-MLP model’s ability to achieve high-precision predictions.

Overall, although individual features exhibit complex nonlinear influences in their SHAP values, from the perspective of overall feature importance ranking and dominant patterns, these observations align highly with morphological findings: after DIN stimulation, cells contract, become rounded, and gradually detach, leading to a significant increase in brightness in the cell regions, with grayscale values concentrated in the 200–255 range (as shown in Figure 4). Therefore, SHAP analysis confirms, at the level of model decision-making mechanisms, that the model ultimately relies on and integrates the information carried by these high-grayscale features, thereby successfully mapping the biological causal relationship between “the degree of cellular damage (manifested as an increase in high-brightness pixels) and the rise in pollutant concentration.”

By monitoring the trends of training loss and validation loss, we systematically optimized the hyperparameters of the SR-MLP model. The optimized hyperparameter configurations are presented in Table 5. As shown in Figure 7, after tuning, both the training loss and validation loss exhibited a smooth downward trend and eventually stabilized, maintaining a small gap between them without evident divergence or premature convergence. These results indicate that the model did not suffer from overfitting or underfitting during training, demonstrating good convergence behavior and generalization ability (Wang et al., 2022). This validates the effectiveness of the adopted regularization strategies and hyperparameter configurations.

Based on the comparison of model performance, we conducted an in-depth analysis of the fundamental reasons behind the superior performance of the SR-MLP model. Analysis reveals that the core reason for the superior predictive performance of SR-MLP over SR-LR and SR-RR is its ability to model nonlinear relationships. The SR-LR and SR-RR model are based on linear regression, which assumes a simple linear relationship between the input features and pollutant concentrations, whereas SR-MLP is able to capture more complex patterns through a multilayered neural network and nonlinear activation functions. To validate the nonlinear relationship between the variables, we visualized the relationship between features and concentrations and used polynomial fit curves to observe data trends. Common manifestations of nonlinearity: the fitted curve is not a straight line, has inflection points, and exhibits significant curvature or fluctuation (e.g., U-shape (Yang et al., 2020), inverted U-shape (Yang et al., 2020; Northoff and Tumati, 2019), and S-shape (Xiao et al., 2009), etc.). By observing the variable fitting curves of HeLa cells as in Figure 8, we found that the fitting curves of most of the features showed obvious nonlinearities, such as Feature 115 and Feature 133 showed U-shape variations. The results indicate that SR-MLP overcomes the limitations of the SR-LR and SR-RR linearity assumption by modeling nonlinear relationships, and is able to more accurately describe the complex mapping between the gray-scale features of cell damage and pollutant concentrations, improving the fitting ability and prediction accuracy of the model.

Having established the core advantage of the SR-MLP model in handling nonlinear relationships, we further applied it to quantitatively predict environmental pollutant concentrations based on cellular image features. This approach establishes a direct, data-driven correlation between cellular morphological damage severity and pollutant concentration, offering a novel research perspective for environmental toxicity assessment. The efficacy of the SR-MLP framework in environmental concentration prediction is corroborated by related studies. For instance, Liu et al. successfully calibrated air quality monitoring data using a similar hybrid model (Liu et al., 2021). By transferring and validating this effective modeling strategy at the cellular level, this study proposes a novel, efficient, and reliable computational method for predicting pollutant concentrations based on cellular damage features. This expands the application scope of the hybrid model to microscopic biological response analysis.

To validate the applicability of the method in different cell types, we applied the SR-MLP model to four other cell types (A549, HUVEC, PC12, HaCaT) and evaluated the performance of the model on different cell datasets. As can be seen from the data in Table 4, the SR-MLP model showed relatively small errors between the predicted and true values at each concentration level in the HeLa, A549, HUVEC, and PC12 cell datasets, whereas the predictions of certain concentrations in the HaCaT cell dataset were slightly biased, e.g., the predicted value of 2.92 μmol/L at 5 μmol/L true value, but the predictions of other concentrations were still accurate. This indicates that the model still has good adaptability and generalization ability across different cellular datasets. We found that despite the differences in biological properties and data distribution among cell types, the overall trends were highly similar, and the SR-MLP model was able to accurately predict contaminant concentrations and compare the degree of cellular damage in different cell models.

To ensure that changes in image grayscale features primarily reflect the toxic effects of pollutants, this study implemented systematic experimental design and data processing procedures to control for potential confounding factors. Specific measures included: during the cell culture phase, all cell lines were cultured simultaneously in a standard incubator under constant temperature and humidity, using the same batch of culture medium and serum, while controlling passage number and inoculation density to maintain approximate consistency; during the image acquisition phase, all samples were imaged using the same optical microscope and a 10× objective, with fixed exposure time, light intensity, and white balance parameters, and each imaging session was preceded by background correction and uniform field calibration; during the image analysis phase, an automated segmentation pipeline based on a pre-trained model (Cellpose cyto3) was employed, with unified flow threshold and minimum cell size settings to eliminate morphological and grayscale biases due to segmentation variability; during data processing, rigorous outlier removal and data normalization procedures were applied.

Nevertheless, we recognize that certain potential confounding factors are difficult to fully eliminate, such as intrinsic cellular heterogeneity (e.g., cell cycle distribution, metabolic state fluctuations) (Kondo et al., 2021), spontaneous morphological changes unrelated to toxicity that may occur during apoptosis or stress, and potential edge illumination attenuation in microscopic imaging fields despite calibration. These unmeasured or incompletely controlled variables may lead the model to partially rely on non-toxicity-related imaging artifacts or biological noise. Therefore, the quantitative “image-concentration” relationship established in this study should be interpreted as a strong statistical association captured under highly controlled experimental conditions. While this association provides a robust computational foundation and scalable framework for image-based toxicity prediction, the underlying biological causal chain warrants further systematic validation in the future through approaches such as intervention experiments (e.g., co-imaging with oxidative stress probes), multimodal data integration (e.g., combining fluorescence labeling with transcriptomics), and causal inference modeling.