The work presented in this paper shows that ML models combined with the feature importance method SHAP are able to detect proteins that are relevant for DED and MGD. Because one aim of the study was to explore proteomic changes in subjects with increasing severity of disease, the models were trained on patients with MGD only. Among the top 15 ranked features in each of the four ML models, there were several proteins that are interesting with respect to DED and MGD. These proteins are PRP4, S100 calgranulin A8, lysozyme C, prostaglandine reductase 1, psoriasin, serotranferrin, ADH7, the polymeric immunoglobulin receptor, and the immunoglobulins IGKV2-24, IGLV8-61, IGLV6-57, IGA2 and IGG1. Their relationships to DED and MGD are discussed below.

Jackson et al.⁸ performed a literature review about proteins in DED and identified some proteins that are promising biomarkers for DED. Three of the proteins detected with the method proposed in the current study were found to be potential biomarkers that are upregulated or downregulated in DED according to the review: PRP4, S100 calgranulin A8 and lysozyme C. These three proteins are marked in Fig. 2a–d. Compared to healthy controls, the expression of PRP4 and S100 calgranulin A8 were significantly different in both aqueous-deficient and combined aqueous-deficient and evaporative DED, while lysozyme C was only significantly differentially expressed in combined aqueous-deficient and evaporative DED⁶.

Further on, prostaglandin reductase 1 is involved in the breakdown of prostaglandins⁴⁴. Prostaglandins trigger inflammatory responses on the ocular surface and were considered as a contributor to the pathogenesis of DED²⁷. Prostaglandin reductase 1 might therefore be a key protein with respect to DED and MGD.

Psoriasin, being associated with the immune system and inflammation, has been found to be present at high levels in MGs from individuals without any diseases affecting or involving the lacrimal glands³⁴. Moreover, the gene for psoriasin was upregulated in diseased MGs³⁵ as well as in saliva from patients with primary Sjögren syndrome (pSS)³⁶. It should, however, be noted that all included subjects in the latter study were females. Taken together, psoriasin could be of value for the grading of MGD and serve as a potential biomarker for DED.

Serotransferrin binds and transports iron and exhibits antimicrobial effects by keeping iron unavailable to pathogens⁴⁵. MGD is associated with microbial infections. Serotransferrin was increased in aqueous deficient and combined aqueous deficient and evaporative DED⁶. Moreover, levels of serotransferrin increased with increasing age⁴⁶. Since the risk of experiencing DED also increases with increasing age, this indicates a potential relationship between serotransferrin and DED. Also considering its antimicrobial effects and that the protein levels were altered in several types of DED, serotransferrin might play a role in the pathogenesis of of MGD.

ADH7 is part of a class IV alcohol dehydrogenase primarily involved in the oxidation of retinol to retinaldehyde and possibly retinoic acid synthesis^47,48. Retinoids are associated with the proliferation, differentiation, keratinization and apoptosis of corneal epithelial cells and deficiency of vitamin A can cause both DED and keratopathy³². A metabolite of vitamin A, 13-cis retinoic acid, is known to cause MGD and has been shown to alter gene expression related to inflammation and differentiation as well as inducing apoptosis of MG epithelial cells in vitro³³.

The polymeric immunoglobulin receptor was significantly downregulated in patients with aqueous-deficient DED compared to individuals without DED²³. The protein was also decreased in tears from rabbits with Sjögren syndrome-associated dry eye²⁴. Moreover, the protein was suppressed in reflex tears, that are produced in response to irritant stimulations²⁵. Huang et al. found, on the contrary, that the polymeric immunoglobulin receptor was upregulated in tears from patients with DED²⁶. Despite contradicting findings, changes in the protein expression have been detected in different types of DED, suggesting that it could play a role in the pathogenesis of the disease. Further on, IGA2 is important for the adaptive immune response and binds to the polymeric immunoglobuline receptor²⁹. Amorim et al.²⁸ found the levels of IGKV2-24 to be significantly higher in tears from patients with proliferative diabetic retinopathy compared to non-diabetic controls. In this study, \(74 \%\) of the proliferative diabetic retinopathy patients exhibited Schirmer test results \(< 10\) mm/5 min, and tear film breakup-times were significantly reduced compared to the controls. In our dataset, FBUT was significantly lower for patients with MGD level 4 compared to MGD levels 2 and 3. Consequently, the importance of IGKV2-24 in our models could also be affected by differences in FBUT between the patient groups. Still, considering their roles in the immune response and the fact that MGD facilitates microbial infections in the eye⁴⁹, IGKV2-24’s, IGA2’s and the polymeric immunoglobulin receptor’s roles in MGD could be investigated further.

Both IGLV8-61 and IGLV6-57 are immunoglobulins involved in the immune response. In a study investigating tear proteomics following laser-assisted in-situ keratomileusis (LASIK) and small incision lenticule extraction (SMILE) surgery, IGLV8-61 and IGLV6-57 were both upregulated one week after LASIK surgery³⁰. The patients had a mean tear film break-up time of 6.1 and Schirmer test of approximately 9 mm/5 min, which are clinical signs of DED. Even though the observed alterations in IGLV8-61 and IGLV6-57 probably arose from the surgery, they could also be connected to dry eye-related pathology.

IGG1, also involved in the adaptive immune response, could play a role in the development of dry eyes. Mackie et al. reported increased tear levels of immunoglobulin gamma in patients with DED compared to healthy controls³¹. A study comparing allergic conjunctivitis in mice with and without DED observed significantly increased levels of IGG1 regardless of DED, suggesting that DED did not affect the protein expression⁵⁰. However, since allergic conjunctivitis activate the adaptive immune response, this might mask possible alterations during DED. A reduction of IGG1 in Descemet’s membrane, which is a part of the cornea, was observed for patients with Fuchs endothelial corneal dystrophy⁵¹. Taken together, alterations of IGG1 expression are observed for several diseases on the ocular surface, including DED.

When weighting the SHAP importance values by the model uncertainty, two additional proteins that are potentially important for MGD were identified. These proteins were glutathione peroxidase 1 and dynactin subunit 2. Moreover, psoriasin, IGG1 and ADH7 were ranked among the top 15 features for the multiclass model, which was not the case without uncertainty weighting. The relevance of glutathione peroxidase 1 and dynactin subunit 2 with respect to MGD and DED is discussed below.

Studies have shown that oxidative stress plays a role in the mechanism of DED⁵². Glutathione peroxidase 1 controls the levels of reactive oxygen species (ROS) and is important for the normal function of the MGs³⁷. Hyperosmolarity, which is often observed with DED and pSS are both associated with reduced levels of this protein^38,39. Consequently, glutathione peroxidase 1 stands out as a potentially promising biomarker for MGD.

Dynactin subunit 2 is a part of the dynactin complex, which binds to and activates dynein. Further on, this facilitates exocytosis from acinar cells, which reside in the lacrimal glands and secrete proteins onto the ocular surface⁴⁰. Studies of mouse models for pSS indicate that the protein secretion from acinar cells are altered independently of inflammatory responses^53,54. Future research should therefore look into whether human lacrimal gland secretion is affected by changes in levels of the dynactin subunit 2 and whether such changes are associated with MGD and DED.

Regarding the significant proteins detected by the PEAKS X Pro software, thymidine phosphorylase is potentially relevant for DED and MGD. Thymidine phosphorylase is associated with endothelial cell survival and function⁴¹. Loss-of-function has been reported to give dry eyes in a case report, although as part of a complex systemic clinical picture⁴². Further research is required to investigate the potential role of thymidine phosphorylase in MGD.

Remarkably, only one of the significant proteins detected by PEAKS Pro X, mannose-1-phosphate guanyltransferase alpha, was also found by the feature importance approach. Consequently, it seems like these two approaches complement each other, and applying both of them can result in a higher number of identified potential biomarkers for diseases.

When using the feature importance approaches, most of the proteins identified as promising biomarkers for MGD are involved in the immune response and inflammation. Other proteins, such as PRP4, ADH7 and dynactin subunit 2, contribute to the normal function of the eye. These findings confirm that MGD is a complex condition associated with disturbance of processes in the healthy eye and activation of several inflammatory and immunologic pathways. There is most likely not one single biomarker, but rather a panel of biomarkers, that characterize MGD.

Even though several of the important features identified in this study represented proteins that were previously found to be related to DED, the majority of the features did not. Highly ranked proteins such as mannose-1-phosphate guanyltransferase alpha, 60S ribosomal proteins L6 and L13 and corticosteroid-binding globulin have not been extensively studied with respect to DED and MGD. However, mannose-1-phosphate guanyltransferase alpha is associated with disease in general, and the other proteins are involved in the regulation of inflammation. Investigating the relevance of these and other highly ranked proteins with respect to DED and MGD might give rise to novel medical discoveries.

Some of the detected proteins that are potential biomarkers for MGD, for example psoriasin and the polymeric immunoglobulin receptor, are altered in pSS. It is not unreasonable that some proteins can play key roles in both MGD and pSS, especially because pSS has dry eyes as a primary symptom and has been associated with MGD^55,56. Still, the most upregulated tear protein in pSS, neutrophil gelatinase-associated lipocalin⁵⁷, was not among the detected proteins in the current study. One possible explanation is that patients with MGD make up a more diverse group, also including patients without pSS. Even though MGD and pSS can occur simultaneously, these are two separate diagnoses, and it is expected that there are differences in how they affect the ocular surface and tear composition.

The feature importance plots differed between the MGD models. This is not surprising, as we expect different protein levels to be of higher or lower importance according to the severity of the MGD. Because MGD is associated with inflammation, one would for example expect proteins related to inflammation to be ranked differently for the different levels of MGD. Still, many of the same proteins were highly ranked for all models. This means that there were also similarities between the different levels of MGD.

The reliability of the estimated feature importance values will be affected if the ML models are not performing well. In this work, the balanced accuracy for the multiclass ML model was 72%, while the balanced accuracies for the three binary ML models were 85% or higher. The high model performances indicate that the predicted important features should be robust. At the moment, there exists no method that takes the model uncertainty into account when determining the feature importance. In the current work, the estimated importance values were weighted based on the predicted probabilities as a mean of including uncertainty into the feature ranking. When the ground truth annotations are available, the importance values for correct predictions with high certainty from the model can be assigned more weights than incorrect and/or uncertain predictions. By following this, a change in the feature ranking was observed, and new potentially MGD-related proteins were identified. However, when the ground truth is not available, the correct model predictions cannot be identified. In this case, only the original importance values for the most certain model predictions were included. Using this technique, no new proteins were present among the 15 most important features. This is probably due to the low number of uncertain predictions. The results indicate that the weighting is most effective when the ground truth is known.

The dataset applied in this study has several strengths. First, it includes a high number of patients. Moreover, the comprehensive proteomic analyses provide detailed information about the tear composition for the patients and were combined with several clinical parameters including level of MGD. Finally, because all patients were examined at the same eye clinic, variations between the performance of the various procedures are expected to be minimal. Taken together, these strengths are adding reliability to the reported results.

This study has some potential limitations. First, the ML models were not evaluated using an external test set. As a result, the models’ abilities to generalize to new patients are not known. The main reason for not dividing the dataset into training and test sets was because the aim was to describe patterns in the data rather than developing ML models for predictive purposes. This approach is motivated by a method called ‘microscope artificial intelligence (AI),’ where models are developed and explained in order to gain a deeper understanding of the data used to train the models^58,59. Similar to applying a microscope for studying our surroundings in more detail, explaining the ML model can enable us to view the data from a different angle, potentially leading to new discoveries. According to microscope AI, the goal is to extract knowledge out of the data rather than creating models for automatic decision-making⁵⁸. Still, the ML models from the current work should not be used for diagnostic purposes since they have not been externally evaluated.

Another potential limitation is the chosen tear collection method. The tears in the present study were collected using Schirmer strips, adhering to previous protocols presented by our group^57,60. However, the optimal methodology for collecting tear samples remains an area of debate. Collection methods such as microcapillary tubes with or without saline flush, Schirmer strips as well as the intrinsic individual variability may contribute to both composition and concentration of collected proteins. Some studies found similar results regarding protein expression following collection with microcapillary tubes and Schirmer strips^61,62, while another study indicates that the collection method might impact the proteins detected in the sample⁶³. In the present project, unanesthetized Schirmer tests were conducted, which can stimulate reflex tearing. Indeed, at least 15 proteins have been noted to differ between reflex and basal tears²⁵, where basal tears might be regarded as more relevant when studying MGD. However, it can be argued that the previous studies mentioned above indicate that the applied collection method using Schirmer strips ensures a representative collection of biological material as well as a low degree of contamination. Moreover, although an anesthetized Schirmer test might be more representative concerning basal tears, the topical anesthesia might serve as a contaminant and alter tear composition.

DED is broadly divided into aqueous deficient DED and evaporative DED⁶⁴. The latter is the more common of the two and MGD is the most common cause of evaporative DED. These subdivisions of DED are not mutually exclusive. Instead, they often overlap along a spectrum referred to as mixed DED, which is commonly seen in the clinic. In the MGD dataset used in the current work, some of the patients with MGD also exhibit signs of aqueous deficient DED. This might be a potential limitation because some proteins have been noted to be altered as a result of lacrimal tear production⁶⁵, which further on can give rise to differences in protein measurements between individuals with and without aqueous deficient DED. On the other hand, the fraction of patients with Schirmer test values \(< 5\) mm was relatively low and stable for each of the MGD levels 2 to 4, ranging between 16 and \(22\%\). Consequently, if the protein measurements were affected due to reduced tear production, it would most likely have similar effects across all these MGD levels. The inclusion of patients with mixed DED probably had little effect on the reported results. Still, it would be interesting to conduct studies specifically targeting tear proteomic patterns in patients with mixed aqueous deficient DED and evaporative DED. While the dataset in the present work only includes 46 out of 233 patients with mixed aqueous deficient DED and MGD, data collection from a larger cohort is necessary to get reliable results for this subgroup of DED patients. Future work should look into this topic.

This study explored which proteins are regarded as important for ML models predicting different levels of MGD. Several of the detected features represented proteins that from earlier research are known to be altered in DED. What might seem unusual with the presented work is that a control group without MGD was not included for training the ML models. However, the aim was not to distinguish healthy eyes from eyes diagnosed with MGD, but rather to explore the expression of proteins in tears for different levels of MGD. Individuals without MGD were not regarded as relevant to include in the training dataset because they might obscure the proteomic differences within the MGD patient group. However, future work should look into how the current results compare to a control group without MGD. Investigating the protein rankings in healthy controls can strengthen the findings in the present study.

In conclusion, this study successfully combined ML and proteomics to explore relationships between MGD severity and protein expression in tears. Examination of which proteins the ML models regarded as important for predicting levels of MGD showed that several of the highest ranked proteins are known to be upregulated or downregulated in DED. Other proteins that might be relevant in MGD were also detected. Future work should explore these proteins further with the aim to increase the understanding of MGD and develop improved treatment options. Moreover, the proposed method for detecting potentially relevant proteins could be applied to other types of medical conditions and diseases with the aim to discover new medical knowledge.