Task-Specialized Protein Language Models Decode the Sequence Grammar of Post-Translational Modification Sites

Post-translational modifications (PTMs) regulate protein signaling, localization, degradation, and cellular decision-making, yet the sequence determinants that distin- guish modified from chemically eligible but unmodified residues remain difficult to de- code at proteome scale. Here, we examine whether adapting a general protein language model to PTM-site prediction can reveal the biochemical logic underlying residue-level modification. We fine-tune ESM2, a protein language model trained on tens of mil- lions of evolutionarily diverse protein sequences, for phosphorylation, acetylation, and ubiquitination-site prediction. To address the pronounced class imbalance inherent in proteome-wide PTM annotation, we combine parameter-efficient fine-tuning with focal-loss training. The resulting task-specialized models show that PTM recognition depends on model capacity, annotation depth, and modification chemistry: phospho- rylation benefits from larger models, whereas acetylation and ubiquitination peak at intermediate scale. More importantly, the fine-tuned phosphorylation model exposes three layers of biological organization: it recovers canonical kinase-recognition motifs without kinase-label supervision, resolves pathway-level functional relationships among 1 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint proteins from sequence-derived embeddings, and preserves evolutionary signatures of homologous phosphorylation sites across 200 eukaryotic species. These results estab- lish task-specialized protein language models as interpretable instruments for probing PTM-site biochemistry, kinase specificity, functional organization, and evolutionary conservation.

Post-translational modifications (PTMs) are pivotal biochemical processes that regulate pro- tein function and enable dynamic cellular responses to diverse physiological and environmen- tal signals. By covalently attaching chemical groups such as phosphate, acetyl, or ubiquitin to specific amino acid residues, PTMs expand the structural and functional diversity of pro- teins far beyond the genomic blueprint. With over 400 distinct PTM types identified,1 these modifications transform an estimated 20,000 protein-coding genes into a proteome exceeding a million unique protein species.2 PTMs orchestrate critical biological processes, including signal transduction, enzymatic regulation, protein localization, and degradation, thereby serving as molecular switches that fine-tune cellular behavior.3–5 Their dysregulation is im- plicated in numerous diseases, from cancer and neurodegenerative disorders to metabolic syndromes, underscoring their significance in health and disease.6,7 Accurate prediction of PTM sites is essential for decoding protein function, understanding disease mechanisms, and developing precision therapeutics. Yet this prediction problem is biologically subtle: most serine, threonine, tyrosine, and lysine residues are chemically eligible for modification, but only a small subset are selected in vivo by the combined influence of local sequence context, enzyme specificity, protein function, and evolutionary constraint. The importance of PTMs lies in their ability to introduce functional versatility without altering the genetic code. Phosphorylation regulates protein activity in pathways controlling cell division and apoptosis;8 acetylation modulates chromatin structure and gene expres- sion;9 and ubiquitination tags proteins for degradation, maintaining cellular homeostasis.10 2 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint The combinatorial nature of PTMs, where multiple modifications coexist on a single protein, further amplifies their regulatory complexity, resembling a molecular language in which se- quence context, structural environment, and modification state jointly determine regulatory meaning.11 This analogy is useful because PTM-site selection is also context dependent: the same residue can be modified or ignored depending on its surrounding sequence, structural environment, and regulatory role.12 Protein language models therefore offer a route to ask whether the biochemical grammar of PTM recognition can be inferred directly from protein sequence. The foundation of computational PTM prediction rests upon curated databases that aggregate experimentally verified modification sites from high-throughput mass spectrome- try and literature mining. General-purpose repositories such as PhosphoSitePlus,13 Phos- pho.ELM,14 PHOSIDA,15 and the Human Protein Reference Database16 provide curated PTM annotations across model organisms, while Swiss-Prot/UniProtKB17 and RESID18 serve as authoritative cross-PTM knowledgebases. Modification-specific resources include UbiProt19 and UbiNet 2.020 for ubiquitination, while integrated platforms such as SysPTM 2.0,21 iPTMnet,22 and the Eukaryotic Phosphorylation Site Database (EPSD)23—which alone houses over 1.6 million phosphorylation sites across 68 species—consolidate multi- PTM data with pathway and ontology annotations. Structure-resolved information is avail- able through Phospho3D24 and PTM-SD,25 and the dbPTM resource,26 recently updated to dbPTM 2025,27 remains among the most comprehensive integrated platforms. Collectively, these databases form the training and validation backbone for downstream computational prediction tools. Computational prediction of PTM sites has progressed from motif- and profile-based scoring to classical machine-learning classifiers and, more recently, deep neural networks that learn modification-specific sequence features directly from data. Early methods used kinase motifs, position-specific scoring, support vector machines, random forests, and physicochem- ical encodings to capture the local sequence environments associated with phosphorylation, 3 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint acetylation, ubiquitination, and related PTMs.28–36 Deep learning subsequently enabled con- volutional, recurrent, and transformer-based models to learn longer-range sequence depen- dencies and context-dependent PTM signatures, improving residue-level prediction across multiple modification classes.37–42 Protein language models (PLMs) have further transformed the landscape by providing contextual residue embeddings learned from millions of natural protein sequences. These embeddings encode information about evolutionary conservation, structural constraints, do- main organization, and functional context without requiring explicit supervision for each property.43–45 Recent PTM predictors increasingly exploit this representational capacity by combining PLM embeddings with task-specific classifier heads, multi-task learning, prompt tuning, structural features, evolutionary profiles, motif-aware modules, or unified multi- label formulations.46–52 These approaches have substantially improved PTM coverage and proteome-scale applicability. Despite this rapid progress, the current generation of PTM predictors leaves several ques- tions open. First, the move toward broad-spectrum, multi-modal architectures has come at the cost of a clear understanding of how the underlying protein language model itself contributes to performance: when ESM2 is held fixed and combined with structural, evo- lutionary, or motif-aware modules, it becomes difficult to attribute gains to the language model versus to the additional modules. Second, although adapting large models for spe- cific tasks while training only a small fraction of parameters is now well-established for natural-language models, this kind of efficient adaptation has not been systematically char- acterized for PTM prediction across model sizes—existing efforts largely use a single model size.47 Third, proteome-wide PTM prediction is intrinsically imbalanced, with phosphoryla- tion accounting for the bulk of curated annotations,53 and even within a single PTM type, the proportion of modified residues is typically below 10% of all candidate sites—a severe imbalance that standard training procedures handle poorly.54 Finally, most predictors are evaluated narrowly on classification metrics, with little attention to whether the learned rep- 4 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint resentations recover known biological patterns such as kinase consensus motifs, functional pathway organization, or evolutionary conservation55—a question that becomes especially important when the goal is not just prediction but mechanistic understanding. Acetylation Phosphorylation Ubiquitination ESM2-8M to 3B parameter models Finetuned Acetylated model Finetuned Phosphorylated model Finetuned Ubiquitinated model Dataset Base Model LoRA Finetuning for predicting PTM sites Biological Significance Kinase Motifs GO : Embedding derived BP Evolutionary Conservation Figure 1: A schematic overview illustrating the development of post-translational modification (PTM)-specific predictive models. ESM2 models, ranging from 8M to 3B parameters, are LoRA- fine-tuned on datasets corresponding to three major PTMs—acetylation, phosphorylation, and ubiquitination. The resulting specialized models are used for accurate prediction of PTM sites. Downstream analyses interrogate the learned representations for three biological readouts: kinase- recognition motifs, pathway-level functional organization, and evolutionary conservation across eu- karyotes. Motivated by these gaps, we ask whether adapting a general protein language model to PTM-site prediction can convert it into an interpretable biological instrument for probing residue-level modification logic. We achieve this specialization through parameter-efficient fine-tuning of ESM2 family of protein language models,44,56 ranging from 8 million to 3 billion parameters and adapting each variant to PTM site classification using Low-Rank 5 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint Adaptation (LoRA).57 LoRA is an efficient adaptation strategy that updates only a small fraction of model parameters and thereby makes adapting billion-parameter models practical on a single GPU. The severe imbalance between modified and unmodified residues is handled by a focal-loss training objective that gives the model stronger feedback on the rare positive cases. Within this controlled setting, we make three contributions. First, we characterize how predictive performance scales across phosphorylation, acetylation, and ubiquitination, finding that the optimal model size depends on the PTM: phosphorylation continues to improve through the 3B-parameter model, while acetylation and ubiquitination peak at 650M parameters and decline at the largest scale—a finding directly relevant to anyone choosing a model size for new PTM tasks. Second, we show that the internal patterns of attention learned by the adapted phosphorylation model recover canonical kinase consensus motifs, with different parts of the model specializing in different kinase families—providing direct mechanistic interpretability that is hard to achieve when the language model is held fixed. Third, we validate the adapted representations along two independent biological axes: similarity-based protein networks recover Gene Ontology functional organization that complements(ratherthanreproduces)theexperimentallycuratedSTRINGinteractome, and similarity at homologous phosphorylation sites correlates strongly with phylogenetic distance across 200 eukaryotic species (Figure 1). Together, these analyses establish that an adapted protein language model is valuable not only as a predictor but as an instrument whose internal representations can be read out as testable biological hypotheses—a complementary axis of value to the coverage-oriented predictors that currently define the field. 6 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint

and Discussion Task Specialization of Protein Language Models for Residue-Level PTM Recognition Each PTM targets a characteristic set of amino acid residues. Phosphorylation targets serine (S), threonine (T), and tyrosine (Y), wherein kinases catalyze the addition of a phosphate group to regulate protein activity and signaling. Acetylation predominantly modifies lysine (K), most prominently in histone proteins, where acetyltransferases add acetyl groups to alter chromatin accessibility and gene expression. Ubiquitination also occurs on lysine, attaching ubiquitin to mark proteins for degradation or to alter their cellular localization. Given a protein sequence, our task is therefore to ask, for every S, T, Y, or K residue, whether it is modified or not. The biological challenge is that only a small fraction of these candidate residues are actually modified under physiological conditions; the model must learn to distinguish true modification sites from the much larger pool of chemically eligible but unmodified residues based on the local sequence context surrounding each site. We approach this problem using ESM2, a protein language model trained on tens of millions of evolutionarily diverse protein sequences. During its pre-training, ESM2 learns to predict amino acids that have been hidden from each sequence—a procedure that forces the model to internalize the recurring patterns by which residues co-occur in nature. As a result, themodel’sinternalrepresentationofeachresiduecapturesinformationaboutitsbiochemical environment, including features related to secondary structure, residue co-evolution, and functional domains, without the model having been explicitly taught any of these concepts. The ESM2 family spans models from 8 million to 3 billion parameters (Figure 2 (a)), with larger variants showing increasingly fine-grained awareness of protein structure—a property that is useful for tasks such as PTM prediction, where the local biochemical environment is often the decisive cue. A pre-trained ESM2 model is, however, a general-purpose protein representation; it cap- 7 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint ESM2- 35M ESM2- 150M ESM2- 650M ESM2- 3B ESM2- 8M Transformer Models(a) (b) EPSD 2.0/dbPTM PTM Database Extract Primary Sequence Filter canonical sequences Phosphorylation Acetylation/Ubiquitination Sequence : Labels : 10 0 0 0 0 0 1 1 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 P... ...N Y G H L T D Q L R G D S S V R D P W T A Q P T Q G Y L T S A P R G P G S L G C P L PTM Residue NON-PTM ResidueA B FFNN with Focal Loss for Binary Classification Base Model LoRA ESM2 (8M-3B) MADELASYHKLHGVYPT Primary Sequence Last Hidden layer Sequence : Labels : 10 0 0 0 0 0 1 1 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 -100 M... ...G K P L H K N Q S R G P K K C Q H A A K D Q R K Q M K S K T N H R A P MK L D C T L Figure 2: Pipeline for PTM site prediction. (a) The ESM2 family of protein language models, ranging from 8M to 3B parameters, serves as the starting point. (b) PTM datasets (EPSD 2.0 for phosphorylation; dbPTM for acetylation and ubiquitination) are processed by extracting primary sequences, filtering for canonical entries, and assigning per-residue labels (modified vs. unmodified) for each candidate residue. Each sequence is passed through the ESM2 model, whose internal representation of each residue is then used as input to a small feed-forward classifier that predicts whether the residue is modified. The ESM2 backbone is fine-tuned using LoRA, which adds a small number of trainable parameters to the model’s attention layers without altering the original weights. A focal-loss training objective addresses the pronounced class imbalance arising from the rarity of true modified residues among chemically eligible candidate sites. tures broad sequence biology but does not by default emphasize the specific features—kinase- binding motifs, local electrostatics, sorting signals—that distinguish a true modification site. 8 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint Adapting the model to PTM prediction therefore requires further training on labeled ex- amples. Doing so by retraining the entire 3-billion-parameter network is computationally prohibitive, while training only a small classifier on frozen ESM2 features can fail to capture the deeper sequence patterns specific to PTMs. We use a middle path: Low-Rank Adapta- tion (LoRA), which inserts a small set of additional trainable parameters into the model’s attention machinery while leaving the original weights untouched. This adjusts how the model attends to neighboring residues without disturbing its general protein knowledge, and reduces the training cost enough that even ESM2-3B can be adapted on a single GPU. A small classifier on top of the adapted model then produces a per-residue prediction (Figure 2 (b)). We therefore treat fine-tuning as a process of task-specialization: the model retains its broad protein-sequence knowledge but is redirected toward the residue-level biophysical fea- tures that distinguish modified from unmodified candidate sites. This proteome-wide class imbalance,where true modified residues constitute only a small minority of chemically eligi- ble candidate sites is handled with focal loss, which gives greater weight to rare and difficult positive examples during training. Architectural and hyperparameter details are given in Methods. Model Capacity, Annotation Depth, and PTM Chemistry Jointly Shape Prediction Performance The fine-tuned ESM2 models reveal a striking pattern: the benefit of fine-tuning depends strongly on PTM chemistry and data depth. In particular, model size interacts differently with each PTM, and increased capacity does not uniformly translate into improved perfor- mance. We evaluate predictions using four standard metrics, each with a clear biological interpretation. Recall measures what fraction of true modification sites the model success- fully identifies—the higher the Recall, the fewer biologically real sites are missed. Precision measures what fraction of the model’s predictions are correct, capturing the rate of false positives that would mislead downstream experiments. F1-score is their harmonic mean. 9 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint The Matthews Correlation Coefficient (MCC) is a single, balanced summary statistic that performs well when modified residues are heavily outnumbered by unmodified ones, as is the case here. We report each metric averaged across both modified and unmodified residues to ensure that performance on the rare modified class is not masked by the much larger unmod- ified class. As a baseline, we compare against PTMGPT2,46 the most directly comparable transformer-based PTM predictor. We first evaluated the performance of pre-trained baseline models on the phosphory- lation site prediction task using the ESM2 models and compared them with PTMGPT2 model. As shown in Figure S1, all baseline ESM2 models (labelled asB-ESM2 in figure S1 ) demonstrated poor predictive performance across all metrics. The ESM2 models achieved F1-scores ranging from 0.39 to 0.49, with AUROC values between 0.44 and 0.54, barely exceeding random chance. Precision and recall metrics were similarly underwhelming, with precision hovering around 0.48-0.61 and recall between 0.44-0.55. Most concerning were the MCC values, which ranged from -0.064 to 0.047, indicating minimal to no correlation between predictions and true phosphorylation sites. Notably, scaling to larger model sizes (650M and 3B parameters) did not yield substantial improvements, with the 650M model actually showing the worst performance (MCC: -0.064, F1: 0.44, AUROC: 0.44). In contrast, PTMGPT2 demonstrated markedly superior performance with an F1-score of 0.59, AUROC of 0.61, precision of 0.61, and MCC of 0.18, substantially outperforming all baseline ESM2 variants. However, even this best-performing baseline model exhibited limited discriminative power for phosphorylation site identification, with an MCC of only 0.182 suggesting weak predictive capability. These results indicate that general-purpose protein language models, despite being pre-trained on large-scale protein sequence data, lack the specific representa- tions necessary for accurate phosphorylation site prediction. The poor zero-shot performance across all baseline models, particularly the near-random AUROC scores for ESM2 variants motivated our decision to fine-tune these models. Fine-tuning the ESM2 variants on the phosphorylation dataset produced substantial 10 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint improvements across all evaluation metrics and model sizes. Performance climbs from fine- tuned ESM2-8M through ESM2-150M (referred asFT-ESM2), with progressive gains in Re- call and MCC indicating better sensitivity to true sites, and ESM2-3B delivers the strongest

overall (MCC 0.35, F1 0.66, highest AUROC; Figure 3 (a)). This is biologically rea- sonable: phosphorylation is by far the best-annotated PTM, with millions of validated sites across hundreds of kinase families, so larger models have enough labeled examples to learn the diverse sequence contexts that different kinases recognize. Acetylation tells a different story. Performance peaks at intermediate scale rather than continuing to improve. The smallest model behaves conservatively, predicting few sites but predicting them confidently (high Precision, lower Recall). Through ESM2-35M and ESM2-150M, F1 improves modestly while Recall remains stable; the best balance between sensitivity and specificity emerges at ESM2-650M (MCC 0.26). At ESM2-3B, all metrics drop including AUROC, suggesting that the largest model has begun to overfit to the com- paratively limited acetylation training set (Figure 3 (b)). Ubiquitination follows a similar trajectory but with a sharper drop at the largest scale. Performance improves through ESM2-150M, peaks at ESM2-650M (MCC 0.34, F1 0.66), and then degrades substantially at ESM2-3B (Figure 3 (c)). As with acetylation, this likely reflects the smaller and biologically narrower set of validated ubiquitination sites available for training—ubiquitination is harder to capture by mass spectrometry than phosphorylation, leading to less complete annotation databases. 11 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (a) (b) (c) Phosphorylation Acetylation Ubiquitination Figure 3: Comparison of the classification metrics for different fine-tuned ESM2 models (labelled as FT-ESM2 )for (a) Phosphorylation, (b) Acetylation and (c) Ubiquitination. Also shown the comparison with baseline model PTMGPT2 model.46 The PTMGPT2 baseline consistently underperforms across all three PTM types when 12 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint compared with fine-tuned models, with notably lower Recall indicating limited sensitivity to true modification sites and lower MCC reflecting weaker overall classification robustness compared to the LoRA-fine-tuned ESM2 models. More broadly, these results show that the relationship between model scale and PTM pre- diction performance is task-dependent. Phosphorylation, with the largest and most diverse training set, benefits monotonically from increased capacity. Acetylation and ubiquitina- tion, in contrast, peak at ESM2-650M and degrade at 3B—a pattern most consistent with overfitting in the larger model when labeled data are limited. Thus, model scaling in PTM predictionisnotgovernedbycapacityalone, butbytheinterplaybetweenmodificationchem- istry, annotation depth, and the diversity of biochemical recognition contexts represented in the training data. 13 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint Fine-tuned Attention Heads Recover Kinase-Specific Phosphoryla- tion Motifs A P M G Y EpSR T M S R L D NpSK P V G T A P HpSM R Input Embedding Positional Encoding Multi-head Attention 0.96 0.19 0.65 0.92 0.31 0.840.57 0.36 0.63 0.25 0.85 0.05 0.01 0.21 0.48 0.07 0.88 0.54 0.49 0.91 0.62 0.06 0.25 0.370.59 0.76 0.81 0.89 0.23 0.17 0.32 0.39 0.93 0.80 0.53 0.30 0.39 0.26 0.43 0.10 0.19 0.44 0.73 0.41 0.61 0.62 0.79 0.99 0.260.96 0.19 0.65 0.92 0.31 0.840.57 0.36 0.63 0.25 0.85 0.05 0.01 0.21 0.48 0.07 0.88 0.54 0.49 0.91 0.62 0.06 0.25 0.370.59 0.76 0.81 0.89 0.23 0.17 0.32 0.39 0.93 0.80 0.53 0.30 0.39 0.26 0.43 0.10 0.19 0.44 0.73 0.41 0.61 0.62 0.79 0.99 0.260.96 0.19 0.65 0.92 0.31 0.840.57 0.36 0.63 0.25 0.85 0.05 0.01 0.21 0.48 0.07 0.88 0.54 0.49 0.91 0.62 0.06 0.25 0.370.59 0.76 0.81 0.89 0.23 0.17 0.32 0.39 0.93 0.80 0.53 0.30 0.39 0.26 0.43 0.10 0.19 0.44 0.73 0.41 0.61 0.62 0.79 0.99 0.260.96 0.19 0.65 0.92 0.31 0.840.57 0.36 0.63 0.25 0.85 0.05 0.01 0.21 0.48 0.07 0.88 0.54 0.49 0.91 0.62 0.06 0.25 0.370.59 0.76 0.81 0.89 0.23 0.17 0.32 0.39 0.93 0.80 0.53 0.30 0.39 0.26 0.43 0.10 0.19 0.44 0.73 0.41 0.61 0.62 0.79 0.99 0.260.96 0.19 0.65 0.92 0.31 0.840.57 0.36 0.63 0.25 0.85 0.05 0.01 0.210.48 0.07 0.88 0.54 0.49 0.91 0.62 0.06 0.25 0.370.59 0.76 0.81 0.89 0.23 0.17 0.32 0.39 0.93 0.80 0.53 0.30 0.39 0.26 0.43 0.10 0.19 0.44 0.73 0.41 0.61 0.62 0.79 0.99 0.26 -10 -6 -2 +2 +6 +100 -10 -6 -2 +2 +6 +100 Probability to Information Bits = Log2 Fine tuned ESM2 Model Last layer attention heads of the finetuned ESM2-3B model Phosphorylation by MAPK Motif : ********P*SP********* Phosphorylation by CDK Motif : ********P*SP*K******* Motif : ********P*SP********* Phosphorylation by GSK Phosphorylation by PKA Motif : *******RR*SL********* Phosphorylation by PKC Motif : *******RR*S*K******** Phosphorylation by PKG Motif : *******RK*SLA******** (a) (b) -10 -1 0-5 +1 +10+5 Figure 4: Kinase consensus motifs recovered from internal attention patterns of the LoRA-fine- tuned ESM2 model. (a) Schematic of the analysis: protein sequences containing phosphorylated residues are passed through the fine-tuned model, and attention values from each phosphorylation site to its surrounding residues (−10 to +10 positions) are extracted from the model’s final layer (40 attention heads). These values are aggregated across all annotated sites in the dataset to produce position-specific amino acid preferences, which are then converted into information-content sequence logos. (b) Sequence logos derived from these patterns reveal the canonical kinase consensus motifs (e.g. MAPK, CDK, GSK, PKA, PKC, PKG) that the model has implicitly learned to associate with phosphorylation, despite never having been told which kinase modifies which site. 14 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint Having shown that fine-tuning improves phosphorylation prediction, we next ask whether the adapted model has learned the biochemical rules of kinase recognition rather than dataset-specific correlations. A predictor can achieve good metrics for the wrong reasons— memorizing dataset artifacts, exploiting compositional biases, or over-relying on sequence neighborhoods that happen to correlate with annotation. If the fine-tuned, task-specialized phosphorylation model has internalized biologically meaningful determinants of phosphory- lation, its attention patterns should recover kinase-recognition motifs without kinase-label supervision. To distinguish biochemical learning from such shortcuts, we interrogated the model‚Äôs attention patterns around annotated phosphorylation sites. Phosphorylation is carried out by kinases, and each kinase family recognizes a charac- teristic short sequence motif around the phosphorylated residue—for example, the proline- directed sites preferred by MAPKs and CDKs, or the basic residues upstream of PKA and PKC sites. These motifs were never given to the model during training; the labels indicate only whether a site is phosphorylated. If the adapted model has internalized real biology, the way it weights the residues surrounding a phosphorylation site—through the model’s internal “attention” mechanism, which encodes how much each residue’s representation depends on each other residue in the sequence—should reproduce these consensus motifs without any supervision over kinase identity. We therefore aggregated these attention values from the final layer of the adapted ESM2-3B model across all annotated phosphorylation sites, within a±10-residue window, and converted the resulting positional preferences into information- content sequence logos. The model has 40 such attention patterns running in parallel (the “heads”), and each head can be inspected separately to see what it has learned to focus on (Figure 4 (a); details in Methods). Acrossthemajorkinasefamiliesknowntoleavedistinctsequencesignatures, theresulting logos cleanly recover the canonical motifs (Figure 4 (b)). The CMGC family—which includes CDKs,54 MAPKs,54 and GSKs58—is well known for proline-directed phosphorylation, and the model recovers this through a single dedicated attention head that captures the canonical 15 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint P*SP* motif (proline at−2 and +1) for both MAPK and GSK substrates, while a separate head captures the CDK-specific extension P*SP*K. The AGC family—PKA,59 PKG,59 and PKC60—is characterized by basic residues upstream of the phosphorylation site, and each member is captured by a distinct head with a refined version of this pattern (RR*SL for PKA, RK*SLA for PKG, RR*S*K for PKC). The CAMK family61 similarly relies on an upstream basic residue, recovered as the R**S motif for both CAMKL and CAMK2. Beyond these family-level patterns, individual kinases also receive dedicated representa- tions: DMPK,62 MAPKAPK,63 and AKT each show characteristic basic motifs, and casein kinase 2 (CK2) recovers its acidic SD*E pattern (Figure S2). The one notable exception is casein kinase 1 (CK1), for which the model did not recover a clear consensus, consistent with CK1‚Äôs known reliance on prior phosphorylation at ‚àí3 rather than a fixed primary- sequence motif. Two points are worth highlighting. The model recovers known kinase motifs without ever being told which kinase modifies which site, and different parts of the model (different heads) specialize in different kinase signatures, suggesting that the model learns a divided-up representation of phosphorylation determinants rather than a single universal motif. To- gether, these observations show that fine-tuning exposes kinase-specific biochemical speci- ficity within the model: distinct attention heads recover distinct phosphorylation-recognition motifs without kinase-label supervision. Although attention weights do not provide a com- plete mechanistic explanation of model decisions, their recovery of known kinase consensus motifs provides a biologically grounded readout of the sequence features emphasized by the fine-tuned model. Task-Specialized Embeddings Resolve Pathway-Level Functional Or- ganization The previous section showed that the adapted model captures kinase-specific motifs at the residue level. We next asked whether task-specialization also reorganizes the protein-level 16 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint embedding space in a biologically meaningful way. In particular, we therefore asked whether the same fine-tuned representation also encodes higher-order functional organization: specif- ically, whether proteins annotated to related biological processes occupy proximal regions of the model-derived embedding space. To test this, we represented each protein as a single embedding vector and constructed a graph in which edges connect pairs of proteins with sufficiently similar embeddings. If the fine-tuned representations encode biologically meaningful functional information, pro- teins participating in the same biological process should be preferentially connected within this graph. This analysis requires both functional annotations and a reference network for comparison. For functional labels, we used the Gene Ontology (GO) Biological Process namespace,64 a hierarchical vocabulary in which proteins are annotated according to the biological processes in which they participate. As a reference network, we used the STRING protein–protein interaction database,64 restricted to high-confidence interactions, which cap- tures physical interactions as well as curated functional associations. Our comparison rests on three metrics from the framework of Chagoyen and Pazos.64 For a given GO term (e.g. “cell cycle”), we extract the subgraph induced by proteins annotated to that term and compute: (i) the density (D), defined as the fraction of possible edges among those proteins that are present; (ii) theclustering coefficient(CC), which quantifies the extent to which neighbors of a protein are themselves connected, thereby measuring lo- cal triangle-richness; and (iii) thesegregation ratio(Sr), which measures whether a protein’s connections are preferentially concentrated within its annotated GO term relative to chance expectation. Thus, density and clustering coefficient quantify local cohesiveness within a functional category, whereas segregation ratio quantifies the specificity with which that cat- egory is separated from the rest of the network. To build the embedding-derived network, we represented each protein by averaging the residue-level representations from the adapted ESM2 3B model into a single vector per protein and connected pairs of proteins whose vectors exceeded a calibrated cosine- 17 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint similarity threshold. The threshold was chosen so that the embedding-derived network and the STRING network contained comparable numbers of edges over the same set of proteins. This calibration controls for trivial differences in global connectivity and allows the com- parison to focus on how each network organizes proteins with respect to biological function (Figure 5 (a); details in Methods). 18 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint go-basic.obo goa-human.gaf Human Protein Annotations Gene Ontology (GO) GO filtering and expansion Filter : Experimental evidence, exclude IEA and NOT Expand : BP terms, is_a & part of relationships Leaf definition : ≥ 3 proteins Extract Unique Proteins EPSD 2.0 PTM Database STRING (confidence score ≥ 700)

Graph Nodes : Proteins Edges : High confidence interactions Finetuned ESM2-3B Model Protein Embedding Generation Residue level embeddings Average and L2 Norm Embedding vector ∑ Pairwise cosine similarity comparison Threshold calibration Sample ≤3000 shared proteins → compute STRING subgraph density → select cosine similarity quantile → constrain to [0.5, 0.99] Embedding derived network (using proteins shared with STRING) Embedding derived network

graph GO term specific subgraph extraction Extract Metrics Aggregate leaf terms into BP categories • Average D, CC, Sr • Compare STRING vs embedding network (a) (b) (c) Figure 5: Methodological workflow and comparative analysis of functional coherence. (a) Pro- tein function annotations were taken from the Gene Ontology (go-basic.obo, goa_human.gaf), restricted to terms supported by experimental evidence and to proteins present in EPSD 2.0. A high-confidence STRING interaction network (score≥ 700) was constructed alongside a similarity network built from the adapted ESM2 3B representations. The similarity threshold for the second network was calibrated to match the STRING network’s edge density on the shared protein set, so that the two networks are directly comparable. Functional coherence within each Gene Ontology biological process category was quantified using three metrics: density (D), clustering coefficient (CC), and segregation ratio (Sr). (b)–(c) Differences ( ∆) between the similarity network and STRING for representative biological process categories. The trend toward∆Sr > 0 and ∆CC < 0 indicates that the similarity network is more functionally specific (proteins in the same category are more selectively connected) but locally sparser (fewer triangles among them) than STRING. 19 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint The two networks reveal distinct but complementary modes of biological organization (Figure 5 (b)–(c)). Across nearly every biological process category, the network derived from model representations shows substantially higher segregation than STRING, indicating that proteins annotated to the same biological process are more selectively connected in the fine-tuned embedding space. The largest gains are observed for transport- and localization- related categories, includingestablishment of localizationand localization, with several-fold increases. Broader categories such asimmune system processand cellular homeostasis also show large gains (Table S1). These trends suggest that the model’s sequence representations encode determinants of biological function such as signal peptides, transmembrane domains, sorting motifs, conserved family signatures, and other sequence-level features that may define functional categories without necessarily producing direct physical interactions. The increase in functional segregation is accompanied by a systematic reduction in local cohesiveness, with average decreases of roughly0.27 in clustering coefficient and0.21 in den- sityacrosscategories. Thisbehaviorisconsistentwiththedifferentinformationcontentofthe two networks. STRING integrates physical interactions, curated pathway memberships, and experimental evidence, and therefore tends to form dense, triangle-rich modules correspond- ing to protein complexes and pathway neighborhoods. In contrast, the embedding-derived network connects proteins on the basis of similarity in the fine-tuned sequence representation space. Such edges can reflect shared functional or evolutionary signatures, but they need not correspond to direct interaction partners or closed local interaction neighborhoods. The metabolic processcategory illustrates this distinction: STRING places metabolic proteins in dense pathway-like modules, whereas the embedding-derived network links them more sparsely while still preserving the overall functional grouping. The cell cycle category provides an informative boundary case, as it is the only major category in which the segregation ratio decreases in the embedding-derived network relative to STRING. This result is consistent with the biological nature of cell-cycle regulation, which depends heavily on transient assemblies, cyclin–CDK partnerships, checkpoint complexes, 20 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint and context-dependent interaction states. Such organization is more directly captured by experimental interaction networks than by primary-sequence-derived representations alone. Thus, rather than indicating a general failure of the embedding-derived network, the cell- cycle case delineates the boundary between functional similarity encoded in sequence and cellular organization that depends on dynamic interaction architecture. Taken together, the embedding-derived network captures a complementary axis of bi- ological organization—one that emphasizes shared sequence, evolutionary, and functional signatures over physical interaction architecture. Rather than reproducing curated interac- tiondatabases, itprovidesanalternativerepresentationinwhichfunctionalgroupingsemerge directly from the fine-tuned protein language model embedding space. Thus, the fine-tuned representation does not simply reproduce curated interaction databases. Instead, it reveals a sequence-encoded layer of pathway organization that is most apparent for processes shaped by shared localization signals, family signatures, and regulatory sequence features. 21 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint AdaptedRepresentationsPreserveEvolutionaryConstraintsonPhos- phorylation Sites Cicer arietinum Pisum sativum Phaseolus vulgaris Arachis hypogaea Manihot esculenta Ricinus communis Cucumis meloCucumis sativusArabidopsis thalianaArabidopsis lyrata subsp lyrataCapsella rubellaRaphanus sativusBrassica oleracea var oleraceaGossypium raimondiiSolanum tuberosumSolanum pennelliiNicotiana tabacumHelianthus annuusSetaria viridisPanicum virgatumMiscanthus sinensisZea maysOryza sativa Indica GroupTriticum aestivumMusa acuminata subsp malaccensisPhyscomitrium patensMarchantia polymorpha subsp ruderalisChlamydomonas reinhardtiiVolvox carteri f nagariensisPlasmodium yoelii yoeliiPlasmodium vivax Sal-1Plasmodium falciparum 3D7Theileria annulataBabesia microti strain RICryptosporidium hominisCryptosporidium parvum Iowa IIToxoplasma gondiiIchthyophthirius multifiliisTetrahymena thermophila SB210Paramecium tetraureliaStylonychia lemnaeThalassiosira pseudonanaTrypanosoma cruzi marinkelleiTrypanosoma brucei gambiense DAL972Leishmania donovaniLeishmania infantum Naegleria gruberi Dictyostelium discoideum Entamoeba histolytica HM-1:IMSS Acanthamoeba castellanii str Neff Trichomonas vaginalis G3 Giardia lamblia ATCC 50803Emiliania huxleyi CCMP1516 Saccharomyces cerevisiae S288C Saccharomyces kudriavzevii IFO 1802 Kluyveromyces marxianus DMKU3-1042 Kluyveromyces lactis NRRL Y-1140 [Candida] glabrata CBS 138 Candida albicans Candida tropicalis MYA-3404Debaryomyces hansenii CBS767 [Candida] auris Yarrowia lipolytica CLIB122 Neurospora crassa OR74APyricularia oryzae 70-15Fusarium graminearum PH-1 Blumeria graminis f sp tritici Aspergillus niger CBS 51388 Aspergillus fumigatus Af293Aspergillus nidulans FGSC A4Schizosaccharomyces pombe 972h- Ustilago maydis Ustilago hordei Malassezia globosa CBS 7966 Cryptococcus neoformans var neoformans JEC21 Cryptococcus gattii VGII Ram5Puccinia graminis f sp tritici Nematostella vectensis Aedes aegypti Culex quinquefasciatus Anopheles gambiaeAnopheles darlingiDrosophila melanogaster Drosophila simulans Drosophila yakuba Lucilia cuprina Ooceraea biroiHarpegnathos saltator Apis melliferaBombus terrestris Nasonia vitripennis Bombyx moriBombyx mandarina Manduca sextaHelicoverpa armigera Papilio xuthus Tribolium castaneum Acyrthosiphon pisum Rhodnius prolixus Zootermopsis nevadensisDaphnia pulex Daphnia magna Limulus polyphemus Ixodes scapularis Caenorhabditis briggsae Caenorhabditis elegans Strongylocentrotus purpuratusPoecilia reticulata Xiphophorus maculatusOryzias latipes Maylandia zebra Takifugu flavidus Takifugu rubripes Gasterosteus aculeatus aculeatus Seriola dumerili Cynoglossus semilaevis Gadus morhua Oncorhynchus mykiss Salmo salar Esox lucius Danio rerio Astyanax mexicanus Lepisosteus oculatusLatimeria chalumnae Xenopus laevis Xenopus tropicalis Corvus brachyrhynchosTaeniopygia guttata Amazona aestiva Dryobates pubescens Aptenodytes forsteri Egretta garzetta Opisthocomus hoazin Calypte anna Columba liviaGallus gallus Meleagris gallopavo Numida meleagris Callipepla squamata Anas platyrhynchos platyrhynchos Tinamus guttatus Struthio camelus australis Alligator sinensis Alligator mississippiensis Pelodiscus sinensis Chelonia mydas Anolis carolinensis Python bivittatus Ornithorhynchus anatinus Vombatus ursinusMonodelphis domesticaLoxodonta africanaTrichechus manatus latirostrisTursiops truncatusLipotes vexilliferDelphinapterus leucasPhyseter catodonCapra hircusOvis ariesBos taurusSus scrofaUrsus maritimusAiluropoda melanoleucaOdobenus rosmarus divergensLeptonychotes weddelliiCanis lupus familiarisMustela putorius furoFelis catusMyotis lucifugusPteropus alectoEquus caballusErinaceus europaeusCricetulus griseusMesocricetus auratusRattus norvegicusMus musculusFukomys damarensisHeterocephalus glaberCavia porcellusDipodomys ordiiIctidomys tridecemlineatusOryctolagus cuniculusPropithecus coquereliOtolemur garnettiiCarlito syrichtaSaimiri boliviensis boliviensisCallithrix jacchusAotus nancymaaeMacaca fascicularisMacaca mulatta Chlorocebus sabaeus Cercocebus atys Papio anubis Mandrillus leucophaeus Rhinopithecus bieti Nomascus leucogenys Pongo abelii Pan paniscus Pan troglodytes Gorilla gorilla gorilla Homo sapiens Figure 6: The circular phylogenetic tree illustrates the broad evolutionary relationships and genetic distances across a diverse range of eukaryotic species, with root at Homo Sapiens 22 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (a) (b) Serine Step 1 : Data curation and clustering Step 2 : Homology Search Step 3 : Multiple Sequence Alignment (MSA) Step 4 : Evolutionary Conservation Analysis FASTA CD-HIT Canonical Proteins Clustering EPSD2 1 2 3 FASTA Human + Homologs S S M S K S T R T T T P Y Y Y L Y Y -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- MAFFT Aligned FASTA Filter Hits Best hit per organism Max Bit score Secondary E-value BLASTP Searches Pan troglodytes FASTA FASTA FASTA Mus musculus Zea mays One FASTA from each cluster of Homo Sapiens 1 2 3 S S M S K S T R T T T P Y Y Y L Y Y -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- -------- MSA Cosine Similarity Extract embeddings from finetuned model Human Embedding Eukaryote Embeddingθ Figure 7: Overview of the evolutionary analysis pipeline for phosphorylation site conservation across eukaryotes. (a) Pipeline. Experimentally validated human phosphorylation sites from EPSD2 are filtered to canonical proteins and clustered with CD-HIT to remove redundancy. Representative sequences from each cluster are used to identify homologs across the panel of eukaryotic species via BLASTP, followed by selection of the best hit per organism and construction of multiple sequence alignments with MAFFT. The fine-tuned model is used to produce its internal representation for both the human residue and its homologous counterpart at each aligned phosphorylation site, and the similarity between the two is quantified by the cosine of the angle between the two representation vectors. (b) Average similarity at serine phosphorylation sites between human and homologs across the panel of eukaryotic species, with bars colored by phylogenetic distance from human and ordered along the eukaryotic tree. 23 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint A natural concern about any deep learning model trained on a fixed dataset is whether it has learned biology or whether it has merely memorized statistical patterns specific to the training distribution. The previous two analyses have already shown that the fine-tuned model recovers known kinase motifs and meaningful functional groupings. Here we apply a third, independent test: does the model encode features that are conserved across evolution? The logic of the test is straightforward. Phosphorylation is an ancient regulatory mech- anism, and many phosphorylation sites are functionally conserved across distantly related eukaryotes even when the surrounding sequence has diverged. If task-specialization has captured the biophysical determinants of phosphorylation rather than memorizing human- specific patterns, homologous phosphorylation sites should remain proximal in the fine-tuned embedding space in a manner that decays with phylogenetic distance. If, on the other hand, the model has merely memorized human-specific patterns, similarity should drop off sharply as soon as we leave the training distribution, with no smooth dependence on phylogeny. To test this, we assembled a panel of 200 eukaryotic species spanning 41 major taxo- nomic groups—vertebrates, invertebrates, fungi, protists, and plants—representing over a billion years of evolutionary divergence (Figure 6). For each human phosphorylation site in a curated, non-redundant set, we identified the homologous residue in each organism via BLASTP and multiple sequence alignment. We then extracted the fine-tuned model’s internal representation (a numerical vector) for both the human residue and its organism counterpart, and measured how similar the two vectors point in their high-dimensional space (the cosine of the angle between them, with1 meaning identical and0 meaning unrelated). We averaged these similarities at the organism level and asked whether the resulting score correlates with phylogenetic distance from human (Figure 7, S3, S4; details in Methods). This analysis reveals a strong, clean negative correlation between phylogenetic distance and similarity for all three phosphorylatable residues: serine (r =−0.86), threonine (r = −0.85), and tyrosine (r =−0.84). The consistency of these correlations across residue types is itself meaningful—the same evolutionary signal is recovered regardless of which amino 24 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint acid gets phosphorylated. As organisms diverge from humans, the model’s representations of corresponding phosphorylation sites become smoothly less similar, mirroring the gradual accumulation of sequence and structural changes at these regulatory positions. The fact that the trend extends from closely related primates all the way to deeply divergent plant and protist lineages—rather than collapsing to noise outside the training distribution—indicates that the model has learned features of phosphorylation that predate the human lineage by hundreds of millions of years. This provides an independent evolutionary validation that the task-specialized representation captures conserved determinants of phosphorylation-site recognition rather than merely human-specific or dataset-specific patterns.

We have shown that parameter-efficient fine-tuning can convert a general protein language model into an interpretable instrument for PTM-site prediction and biological discovery. By fine-tuning the ESM2 family from 8M to 3B parameters with LoRA and a focal-loss objective that handles the rarity of true PTM sites, we obtain strong residue-level predictions for phosphorylation, acetylation, and ubiquitination on a single GPU. More importantly, this design provides a clean setting in which the contribution of the PLM backbone, model scale, and PTM-specific data depth can be examined directly. A central practical finding is that the optimal model size depends on the PTM. Phospho- rylation, with the largest and most diverse training set, benefits all the way to ESM2-3B, which gives the best metrics. Acetylation and ubiquitination peak at ESM2-650M and de- cline at 3B, a pattern most consistent with the larger model overfitting when labeled data are limited. The practical take-away is that backbone size should be matched to the PTM and the amount of labeled data available rather than defaulting to the largest available model. Beyond classification, three independent biological readouts show that the fine-tuned phosphorylation model learns interpretable determinants of PTM-site recognition. Inside 25 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint the model, fine-tuned attention patterns recover canonical kinase consensus motifs—proline- directed for CMGC kinases, basic upstream for AGC and CAMK families, acidic for CK2— without the model ever being told which kinase modifies which site. Networks built from similarity between fine-tuned protein representations recover Gene Ontology Biological Pro- cess organization with substantially higher functional coherence than the experimentally curated STRING interactome, while remaining sparser locally; the one major exception— cell cycle, dominated by transient assemblies that primary sequence does not encode—marks a sensible boundary of what sequence-based representations can recover. Finally, similarity between human and homologous phosphorylation sites correlates strongly and consistently with phylogenetic distance (r≈−0.85) across 200 eukaryotic species, indicating that the model has learned phosphorylation features conserved across more than a billion years of evolutionary divergence. More broadly, this study establishes fine-tuning as a route for transforming general pro- tein language models into biological instruments. Once adapted to PTM-site recognition, their embeddings and attention patterns can be interrogated to reveal kinase specificity, pathway-level functional organization, and evolutionary conservation. Task-specialization therefore provides a framework not only for residue-level PTM prediction, but for extracting mechanistic hypotheses about protein regulatory biochemistry from sequence-trained mod- els.

Dataset Curation and Handling Class Imbalance The model is trained and evaluated on curated datasets for three major post-translational modifications (PTMs): phosphorylation, acetylation, and ubiquitination. Phosphorylation data were obtained from the Eukaryotic Phosphosite Database (EPSD 2.0),65 which provides experimentally validated sites across 223 eukaryotic species. Acetylation and ubiquitination 26 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint data were retrieved from dbPTM,27 a comprehensive repository of experimentally verified PTM annotations. Phosphorylation Dataset.The initial dataset comprised 362,707 proteins, which were refined to 327,516 unique sequences after resolving duplicate groups and removing redundant entries. A quality control step further excluded 48 invalid PTM annotations that did not cor- respond to serine (S), threonine (T), or tyrosine (Y) residues. In total, 35,760,439 candidate residues were identified, including approximately 19.1 million serines, 11.2 million threonines, and 5.3 million tyrosines. Among these, 2,573,592 (7.2%) were experimentally validated phosphorylation sites, while the remaining 33,186,847 (92.8%) were non-phosphorylated, re- sulting in a negative-to-positive ratio of 12.9:1. This imbalance is particularly severe for tyrosine residues (3.6% phosphorylated), compared to serine (9.3%) and threonine (5.4%). The dataset was partitioned into training, validation, and test sets in an 80:10:10 ratio, stratified by protein sequence to prevent data leakage. The training set contains 262,012 sequences (28,584,069 residues, 2,055,505 positives, 7.19%), while the validation and test sets each contain 32,752 sequences with comparable residue distributions and positive rates (7.22%). Residue-specific proportions are preserved across splits, ensuring consistent evalu- ation. Acetylation Dataset.The acetylation dataset consists of 35,572 protein sequences with 21,269,017 total residues, including 1,480,397 lysine (K) residues. Among these, 110,386 (7.5%) are acetylated, yielding a negative-to-positive ratio of 12.4:1. The data were split into training (28,457 sequences), validation (3,557 sequences), and test (3,558 sequences) sets using the same stratified 80:10:10 scheme. The class distribution remains consistent across splits, with approximately 7–8% of lysines labeled as modified. Ubiquitination Dataset.The ubiquitination dataset includes 30,874 protein sequences comprising 17,481,793 residues and 1,084,833 lysine residues. Of these, 118,293 (10.9%) are ubiquitinated, resulting in a comparatively less severe imbalance (8.2:1). The dataset is similarly partitioned into training (24,699 sequences), validation (3,087 sequences), and test 27 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (3,088 sequences) sets, maintaining consistent class proportions across all splits. Akeychallengeacrossalldatasetsisthepronouncedclassimbalanceinherentinproteome- wide PTM prediction: true modified residues are rare relative to the much larger pool of chemically eligible but unmodified candidate sites.66 Standard loss functions such as cross- entropy tend to bias models toward the majority class, leading to poor sensitivity for rare but biologically significant modification sites. To address this, we employ a weighted focal loss function67 with focusing parameterγ = 2.0: Lfocal =−αt(1−pt)γ log(pt), (1) wherept denotes the predicted probability of the true class andαt is a class-specific weighting factor. Class weights are set to 1 and 12.9 for negative and positive samples in phosphoryla- tion, 1 and 12.4 for acetylation, and 1 and 8.2 for ubiquitination, reflecting their respective imbalance ratios. Padding tokens and masked residues are excluded from loss computation to ensure that only valid positions contribute to training. Overall, the framework combines the evolutionary and structural representations cap- tured by ESM244 with parameter-efficient fine-tuning via LoRA57 and imbalance-aware op- timization through weighted focal loss. This design enables robust residue-level prediction performance, maintaining both precision and sensitivity in large-scale, highly imbalanced PTM datasets. Model Architecture and Fine-Tuning Strategy Pre-trained ESM2 models (8M, 35M, 150M, 650M, and 3B parameters) were used as back- bone encoders for contextual residue-level embeddings and fine-tuned with LoRA. In this parameter-efficient fine-tuning scheme, the original backbone weights remained frozen while trainable low-rank updates were introduced into the transformer layers. LoRA modules were inserted into the query, key, value, and output projection layers of 28 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint each attention block.56 Each adapter consisted of a rank-r decomposition withr = 16 and scaling factorα = 32. The original backbone weights were kept frozen, while LoRA parame- ters and all bias terms were trainable. The per-residue classification head was implemented as a feed-forward neural network (FFNN), with architecture dependent on the embedding dimension of each ESM2 variant: • ESM2-8M / 35M: 256→ 128→ 32→ 16→ 8→ 4→ 2 • ESM2-150M: 512→ 256→ 128→ 32→ 16→ 8→ 4→ 2 • ESM2-650M: 1024→ 512→ 256→ 128→ 32→ 16→ 8→ 4→ 2 • ESM2-3B: 2048→ 1024→ 512→ 256→ 128→ 32→ 16→ 8→ 4→ 2 Each layer consisted of a linear transformation followed by layer normalization68 and Tanh activation. Linear weights were initialized using Xavier uniform initialization69 and biases were initialized to zero. Protein sequences were preprocessed by replacing rare amino acids (O, B, U, Z, J) with “X”. Tokenization was performed using the ESM2 tokenizer with a maximum sequence length of 1024 tokens. To ensure alignment with special tokens, residue- level labels were truncated to 1022 positions and padded with ignore indices (−100) at both termini. A custom DataCollatorForTokenClassificationESM was used to maintain con- sistency between input IDs, attention masks, and label tensors. To address class imbalance, a weighted focal loss was used withγ = 2.0. Loss computation excluded padding and masked positions (−100). Training was performed using mixed-precision (FP16) to improve com- putational efficiency70 on NVIDIA A6000 GPUs. The AdamW optimizer71 was used with a learning rate of 5× 10−5 and weight decay of 1.0. A batch size of 2 was used without gradient accumulation. Models were trained for up to 10 epochs with early stopping72 based on validation loss, using a minimum improvement threshold of 0.005 and patience of 5 eval- uation steps. Evaluation was conducted every 1250 training steps. Metrics were computed only on valid (non-padded) residue positions. Regularization was provided by the intrinsic 29 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint dropout layers in ESM2 (dropout rate = 0.1)73 and layer normalization in the classifier. No additional explicit regularization was applied. All models were implemented in PyTorch74 using the HuggingFace Transformers library.75 Custom components included the focal loss function, data collator, and early stopping mechanism. All experiments were performed with fixed random seeds for PyTorch, NumPy, and Python to ensure reproducibility. Model performance was evaluated using Precision, Recall, F1-score, and Matthews Correlation Co- efficient (MCC),76 computed from true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN): Positive class metrics : Precisionpos = TP TP +FP (2) Recallpos = TP TP +FN (3) F1pos = 2· Precisionpos·Recallpos Precisionpos +Recallpos = 2TP 2TP +FP +FN (4) Negative class metrics : Precisionneg = TN TN +FN (5) Recallneg = TN TN +FP (6) F1neg = 2· Precisionneg·Recallneg Precisionneg +Recallneg = 2TN 2TN +FN +FP (7) Matthews Correlation Coefficient (MCC): MCC = TP·TN−FP·FN√ (TP +FP )(TP +FN )(TN +FP )(TN +FN ) (8) These metrics were computed for all models to quantitatively assess classification accu- racy and enable fair comparison between the LoRA-fine-tuned ESM2 and PTMGPT2 model. 30 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint Attention-based Motif Analysis The optimal fine-tuned ESM2 model was used to extract attention weights by performing forward passes with attention outputs for each protein sequence in the dataset. Sequences were tokenized, and the full attention tensors were obtained in evaluation mode. From these, only the final transformer layer was considered, comprising 40 attention heads. To character- ize phosphorylation-associated sequence patterns, attention values were aggregated across all annotated phosphorylation sites. For each attention head independently, a positional window spanning -10 to +10 residues relative to each phosphorylation site was considered. For every position within this window, attention scores from the phosphorylation site to surrounding residues were accumulated conditioned on the amino acid identity, yielding a 21× 20 matrix representing positional amino acid preferences across all sites. The aggregated attention matrices were converted to probability distributions by row-wise normalization, where each row sum was divided into individual elements to obtain relative frequencies of amino acids at each position. Positions with zero total attention were assigned a uniform distribution (1/20 for each amino acid) as a fallback. The probability matrix was then transformed into an in- formation content matrix using logomaker’s77 probability-to-information conversion, which computes bits as log2(probability/background) where background frequencies were assumed uniform (0.05 per amino acid). The information content matrices were visualized as sequence logos78 with amino acid letters stacked proportionally to their information content at each position, using a custom color scheme to distinguish the 20 amino acids. Evolutionary analysis For evolutionary analysis, human phosphorylation data were obtained from the EPSD2 database, including protein sequences and annotated phosphorylation sites. Only canon- ical protein sequences were retained by excluding isoforms, resulting in a curated FASTA dataset of human proteins. To reduce sequence redundancy and group homologous proteins, clustering was performed using CD-HIT79 with a sequence identity threshold of 70%. This 31 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint procedure yielded 20,933 clusters comprising a total of 48,570 sequences, of which 12,451 clusters (59.5%) were singletons. The resulting file was parsed to extract cluster-level statis- tics, including clustersize, representativesequence, representativesequencelength, andmean pairwise identity of non-representative members to the representative. A diversity score for each cluster was defined as1−mean identity/100. To obtain a high-quality and diverse sub- set suitable for downstream analysis, clusters were filtered based on the following criteria: (i) cluster size≥ 3 to ensure sufficient intra-cluster representation, (ii) mean sequence identity < 90% to retain sequence diversity, and (iii) representative sequence length between 100 and 1500 amino acids to exclude very short or excessively long proteins. After filtering, 1,369 clusters were retained. The selected clusters exhibited a median size of 4 sequences (mean 6.2, maximum 410) and an average diversity of 0.162. Cluster sizes were predominantly in the range of 3-5 sequences, with a smaller number of larger clusters. For each selected cluster, the representative sequence (as defined by CD-HIT) was extracted and used for subsequent analysis, yielding a non-redundant and evolutionarily diverse set of protein sequences. The representative sequences from the 1,369 selected clusters were used as queries for large-scale homology search against a diverse set of eukaryotic proteomes. A concatenated protein database comprising proteomes from 200 eukaryotic organisms was constructed by merging individual FASTA files into a single reference file. From this concatenated dataset, organism identifiers were parsed from FASTA headers and a non-redundant list of organisms was generated. A BLAST80 protein database was created from the concatenated proteome. To enable scalable computation, the query FASTA file containing the 1,369 representative human se- quences was split into individual single-sequence FASTA files, each corresponding to one query protein. BLASTP searches were performed independently for each query sequence against the concatenated database using an E-value cutoff of1× 10−5.81 For each query, all potential hits were retrieved, and results were processed to identify the best match per or- ganism. Specifically, for every organism in the dataset, hits were filtered based on organism 32 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint annotation in the subject title field, and the optimal hit was selected based on the highest bit score, with E-value used as a secondary criterion in the case of ties. This procedure resulted in a per-query mapping of the best homologous sequence across all 200 eukaryotic organisms. To enable comparison of phosphorylation sites between human proteins and their ho- mologs across 200 eukaryotic organisms, multiple sequence alignments were performed for each set of homologous sequences using the MAFFT82 alignment program. For each repre- sentativehumanprotein, acorrespondingFASTAfilecontainingthehomologsidentifiedfrom BLASTP was used as input. This procedure resulted in a collection of multiple sequence alignments for each human protein and its homologs, enabling positional correspondence analysis of phosphorylation sites across species. To quantify evolutionary conservation of phosphorylation-related features in the learned representation space, residue-level embeddings were extracted from the fine-tuned model for serine, threonine, and tyrosine residues in both human proteins and their aligned homolo- gous sequences. For each human protein, embeddings corresponding to the aligned phos- phorylation sites were identified using the multiple sequence alignments, ensuring positional correspondence between the human sequence and each homolog. For each aligned residue pair, similarity between the human and organism embeddings was computed using cosine similarity. Given two embedding vectorsu and v, the similarity was calculated as similarity = u· v ∥u∥∥v∥, where u· v denotes the dot product and∥·∥ represents the Euclidean norm. To ensure numerical stability, cases where either vector had zero norm were assigned a similarity score of zero. For each human protein and corresponding organism, cosine similarity values were com- puted across all aligned phosphorylation sites and stored. These values were then aggregated at the organism level by pooling similarity scores across all proteins and computing the mean 33 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint similarity for each organism. This resulted in a single representative similarity score per or- ganism, reflecting the overall conservation of phosphorylation-related embedding features relative to human proteins. To relate embedding similarity with evolutionary divergence, phylogenetic distances for each organism were obtained from an external dataset. Organism-wise average similarity scores were matched with their corresponding phylogenetic distances and sorted in ascending order of distance from human. The results were visualized using a bar plot, where each bar represents an organism, the height corresponds to the average cosine similarity, and the color encodes phylogenetic distance using a continuous colormap. This analysis enabled systematic assessment of how embedding similarity varies as a function of evolutionary distance. Construction of Quantitatively Comparable Embedding and Inter- action Graphs Gene Ontology (GO) data were obtained from thego-basic.obo ontology file,83,84 and hu- man protein annotations were retrieved from thegoa_human.gaf dataset.85 Only biological process (BP) annotations supported by experimental evidence were retained by excluding entries with theIEA evidence code and those annotated with theNOT qualifier. Protein iden- tifiers were restricted to the set of unique human UniProt accessions present in the EPSD 2.0 database. For each protein, direct GO annotations were expanded by traversing the GO hierarchy using is_a and part_of relationships, thereby incorporating ancestral terms. A mapping from GO terms to associated proteins was constructed, and leaf terms were de- fined as those without annotated children within the filtered dataset. Only leaf GO terms associated with at least three proteins were retained for subsequent analysis. Protein-protein interaction data were obtained from STRING v12.0,86,87 retaining only interactions with a confidence score of at least 700. STRING identifiers were mapped to UniProt accessions using the provided alias file, and an undirected graph was constructed in which nodes represent proteins and edges correspond to high-confidence interactions. This 34 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint graph serves as a reference network reflecting experimentally supported functional relation- ships. Protein representations were derived from embeddings generated using the fine-tuned ESM2-3B language model. For each protein, residue-level embeddings were averaged to pro- duce a fixed-length vector, and all vectors were L2-normalized such that the cosine similarity between two proteins corresponds to the dot product of their embeddings. To construct an embedding-based similarity network, pairwise cosine similarities were computed and edges were defined by thresholding these similarities. The similarity threshold was calibrated to match the edge density of the STRING network as closely as possible, following principles used in prior network comparison studies.64 Specifically, a random subset of up to 3,000 pro- teins common to both datasets was sampled with a fixed random seed for reproducibility, and the threshold was chosen as the cosine similarity quantile corresponding to the edge density of the STRING subgraph on those proteins. The calibrated threshold was subsequently con- strained to the interval[0.5, 0.99] to ensure numerical validity. The final embedding-derived network was constructed using only proteins shared between the STRING and embedding datasets. To quantify functional coherence, GO term-specific subgraphs were extracted and ana- lyzed using three graph-based metrics.88 The density (D) was defined as the ratio of observed edges to the maximum possible number of edges within the induced subgraph. The cluster- ing coefficient (CC) was computed as the average local clustering coefficient across nodes within the induced subgraph of each GO term,89 considering only edges between proteins annotated to that term. Nodes with degree one within the induced subgraph were assigned a local clustering coefficient of one, as their single neighbor trivially satisfies the clustering condition, while isolated nodes were assigned a value of zero. Functional segregation was quantified using the segregation ratio (Sr), defined as Sr = ⟨ kin i /ki ⟩ n/N , (9) 35 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint wherekin i denotes the number of neighbors of nodei within the GO term subgraph,ki is its total degree in the full graph,n is the size of the GO term, andN is the total number of nodes. Proteins with zero degree in the full graph were excluded from the segregation ratio computation. Leaf GO terms were grouped into higher-level biological process categories based on the GO hierarchy, and all three metrics were averaged across leaf terms within each category. The number of leaf terms, average term size, and mean values ofSr, CC, and D were reported for both the STRING network and the embedding-derived network, enabling a direct and controlled comparison of their ability to capture biologically meaningful functional organization. Code and Data A vailability The code for PTM site prediction models is publicly available athttps://github.com/ JMLab-tifrh/PTM_prediction. The repository includes pre-trained models for acetylation, phosphorylation, and ubiquitination site prediction, along with Jupyter notebooks for run- ning predictions on new protein sequences. The datasets used for model training and evalua- tion were obtained from publicly available databases: EPSD 2.0https://epsd.biocuckoo. cn/Download.php#1 and dbPTM https://biomics.lab.nycu.edu.tw/dbPTM/download. php#dataset. Supporting Information The Supporting Information (SI) provides supplemental figures, supplemental table. Figure S1 shows the comparison of base ESM2 models with PTMGPT2 model. Figure S2 shows the Logo map for kinase specific motifs. Figure S3 and S4 shows the cosine similarity of threonine and tyrosine phosphorylation across eukaryotic species. Table S1 shows the comparison of STRING and ESM2 network derived metrics for different GO terms. 36 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint

We acknowledge support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4007. We sincerely acknowledge Tata Institute of Funda- mental Research Hyderabad, India for providing the support of computing resources. JM acknowledges the core research grant approved by the Department of Science and Technology (DST) of India (CRG/2023/001426).

(1) Mann, M.; Jensen, O. N. Proteomic analysis of post-translational modifications.Nature biotechnology 2003, 21, 255–261. (2) Venter, J. C.; Adams, M. D.; Myers, E. W.; Li, P. W.; Mural, R. J.; Sutton, G. G.; Smith, H. O.; Yandell, M.; Evans, C. A.; Holt, R. A., et al. The sequence of the human genome. science 2001, 291, 1304–1351. (3) Walsh, C. Posttranslational modification of proteins: expanding nature’s inventory; Roberts and Company Publishers, 2006. (4) Adhikari, S.; Mondal, J. Machine learning subtle conformational change due to phos- phorylation in intrinsically disordered proteins.The Journal of Physical Chemistry B 2023, 127, 9433–9449. (5) Sarkar, D. K.; Adhikari, S.; Surolia, A. Integrating Molecular Dynamics and Deep Learning to Elucidate Conformational Plasticity Underlying the Reduced Activity of Glycocin F.The Journal of Physical Chemistry B2026, (6) Deribe, Y. L.; Pawson, T.; Dikic, I. Post-translational modifications in signal integra- tion. Nature structural & molecular biology2010, 17, 666–672. 37 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (7) Olsen, J. V.; Mann, M. Status of large-scale analysis of post-translational modifications by mass spectrometry.Molecular & cellular proteomics2013, 12, 3444–3452. (8) Blom, N.; Sicheritz-Pontén, T.; Gupta, R.; Gammeltoft, S.; Brunak, S. Prediction of post-translational glycosylation and phosphorylation of proteins from the amino acid sequence. Proteomics 2004, 4, 1633–1649. (9) Choudhary, C.; Mann, M. Decoding signalling networks by mass spectrometry-based proteomics. Nature reviews Molecular cell biology2010, 11, 427–439. (10) Hershko, A.; Ciechanover, A. The ubiquitin system. Annual review of biochemistry 1998, 67, 425–479. (11) Dill, K. A.; MacCallum, J. L. The protein-folding problem, 50 years on.science 2012, 338, 1042–1046. (12) Alipanahi, B.; Delong, A.; Weirauch, M. T.; Frey, B. J. Predicting the sequence speci- ficitiesofDNA-andRNA-bindingproteinsbydeeplearning. Nature biotechnology2015, 33, 831–838. (13) Hornbeck, P. V.; Zhang, B.; Murray, B.; Kornhauser, J. M.; Latham, V.; Skrzypek, E. PhosphoSitePlus, 2014: mutations, PTMs and recalibrations. Nucleic acids research 2015, 43, D512–D520. (14) Dinkel, H.; Chica, C.; Via, A.; Gould, C. M.; Jensen, L. J.; Gibson, T. J.; Diella, F. Phospho. ELM: a database of phosphorylation sites-update 2011.Nucleic acids research 2010, 39, D261–D267. (15) Gnad, F.; Ren, S.; Cox, J.; Olsen, J. V.; Macek, B.; Oroshi, M.; Mann, M. PHOSIDA (phosphorylation site database): management, structural and evolutionary investiga- tion, and prediction of phosphosites.Genome biology2007, 8, R250. 38 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (16) Keshava Prasad, T.; Goel, R.; Kandasamy, K.; Keerthikumar, S.; Kumar, S.; Math- ivanan, S.; Telikicherla, D.; Raju, R.; Shafreen, B.; Venugopal, A., et al. Human protein

database-2009 update.Nucleic acids research2009, 37, D767–D772. (17) UniProt Consortium, T. UniProt: the universal protein knowledgebase.Nucleic acids research2018, 46, 2699–2699. (18) Garavelli, J. S. The RESID Database of Protein Modifications as a resource and anno- tation tool.Proteomics 2004, 4, 1527–1533. (19) Chernorudskiy, A. L.; Garcia, A.; Eremin, E. V.; Shorina, A. S.; Kondratieva, E. V.; Gainullin, M. R. UbiProt: a database of ubiquitylated proteins.BMC bioinformatics 2007, 8, 126. (20) Li, Z.; Chen, S.; Jhong, J.-H.; Pang, Y.; Huang, K.-Y.; Li, S.; Lee, T.-Y. UbiNet 2.0: a verified, classified, annotated and updated database of E3 ubiquitin ligase–substrate interactions. Database 2021, 2021, baab010. (21) Li, J.; Jia, J.; Li, H.; Yu, J.; Sun, H.; He, Y.; Lv, D.; Yang, X.; Glocker, M. O.; Ma, L., et al. SysPTM 2.0: an updated systematic resource for post-translational modification. Database 2014, 2014, bau025. (22) Huang, H.; Arighi, C. N.; Ross, K. E.; Ren, J.; Li, G.; Chen, S.-C.; Wang, Q.; Cow- art, J.; Vijay-Shanker, K.; Wu, C. H. iPTMnet: an integrated resource for protein post-translational modification network discovery. Nucleic acids research 2018, 46, D542–D550. (23) Lin, S.; Wang, C.; Zhou, J.; Shi, Y.; Ruan, C.; Tu, Y.; Yao, L.; Peng, D.; Xue, Y. EPSD: a well-annotated data resource of protein phosphorylation sites in eukaryotes. Briefings in Bioinformatics2021, 22, 298–307. 39 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (24) Zanzoni, A.; Ausiello, G.; Via, A.; Gherardini, P. F.; Helmer-Citterich, M. Phospho3D: a database of three-dimensional structures of protein phosphorylation sites.Nucleic acids research2007, 35, D229–D231. (25) Craveur, P.; Rebehmed, J.; de Brevern, A. G. PTM-SD: a database of structurally resolved and annotated posttranslational modifications in proteins. Database 2014, 2014, bau041. (26) Huang, K.-Y.; Su, M.-G.; Kao, H.-J.; Hsieh, Y.-C.; Jhong, J.-H.; Cheng, K.-H.; Huang, H.-D.; Lee, T.-Y. dbPTM 2016: 10-year anniversary of a resource for post- translational modification of proteins.Nucleic acids research2016, 44, D435–D446. (27) Chung, C.-R.; Tang, Y.; Chiu, Y.-P.; Li, S.; Hsieh, W.-K.; Yao, L.; Chiang, Y.-C.; Pang, Y.; Chen, G.-T.; Chou, K.-C., et al. dbPTM 2025 update: comprehensive inte- gration of PTMs and proteomic data for advanced insights into cancer research.Nucleic acids research2025, 53, D377–D386. (28) Blom, N.; Gammeltoft, S.; Brunak, S. Sequence and structure-based prediction of eu- karyotic protein phosphorylation sites.Journal of molecular biology1999, 294, 1351– 1362. (29) Obenauer, J. C.; Cantley, L. C.; Yaffe, M. B. Scansite 2.0: Proteome-wide prediction of cell signaling interactions using short sequence motifs.Nucleic acids research2003, 31, 3635–3641. (30) Gao, J.; Thelen, J. J.; Dunker, A. K.; Xu, D. Musite, a tool for global prediction of general and kinase-specific phosphorylation sites.Molecular & Cellular Proteomics 2010, 9, 2586–2600. (31) Xue, Y.; Ren, J.; Gao, X.; Jin, C.; Wen, L.; Yao, X. GPS 2.0, a tool to predict kinase- specific phosphorylation sites in hierarchy.Molecular & cellular proteomics2008, 7, 1598–1608. 40 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (32) Khan, Y. D.; Rasool, N.; Hussain, W.; Khan, S. A.; Chou, K.-C. iPhosY-PseAAC: Iden- tify phosphotyrosine sites by incorporating sequence statistical moments into PseAAC. Molecular biology reports2018, 45, 2501–2509. (33) Yang, H.; Wang, M.; Liu, X.; Zhao, X.-M.; Li, A. PhosIDN: an integrated deep neural network for improving protein phosphorylation site prediction by combining sequence and protein–protein interaction information.Bioinformatics 2021, 37, 4668–4676. (34) Ismail, H. D.; Jones, A.; Kim, J. H.; Newman, R. H.; Kc, D. B. RF-Phos: A novel general phosphorylation site prediction tool based on random Forest.BioMed research international 2016, 2016, 3281590. (35) Deng, W.; Wang, C.; Zhang, Y.; Xu, Y.; Zhang, S.; Liu, Z.; Xue, Y. GPS-PAIL: prediction of lysine acetyltransferase-specific modification sites from protein sequences. Scientific reports2016, 6, 39787. (36) Chen, X.; Qiu, J.-D.; Shi, S.-P.; Suo, S.-B.; Huang, S.-Y.; Liang, R.-P. Incorporating key position and amino acid residue features to identify general and species-specific Ubiquitin conjugation sites.Bioinformatics 2013, 29, 1614–1622. (37) LeCun, Y.; Bengio, Y.; Hinton, G. Deep learning.nature 2015, 521, 436–444. (38) Wang, D.; Zeng, S.; Xu, C.; Qiu, W.; Liang, Y.; Joshi, T.; Xu, D. MusiteDeep: a deep-learning framework for general and kinase-specific phosphorylation site prediction. Bioinformatics 2017, 33, 3909–3916. (39) Wang, H.; Wang, Z.; Li, Z.; Lee, T.-Y. Incorporating deep learning with word embed- ding to identify plant ubiquitylation sites.Frontiers in Cell and Developmental Biology 2020, 8, 572195. (40) Liu, Y.; Liu, Y.; Wang, G.-A.; Cheng, Y.; Bi, S.; Zhu, X. BERT-Kgly: A bidirectional 41 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint encoder representations from transformers (BERT)-based model for predicting lysine glycation site for homo sapiens.Frontiers in Bioinformatics2022, 2, 834153. (41) Soylu, N. N.; Sefer, E. DeepPTM: protein post-translational modification prediction from protein sequences by combining deep protein language model with vision trans- formers. Current Bioinformatics2024, 19, 810–824. (42) Wen, B.; Wang, C.; Li, K.; Han, P.; Holt, M. V.; Savage, S. R.; Lei, J. T.; Dou, Y.; Shi, Z.; Li, Y., et al. DeepMVP: deep learning models trained on high-quality data accurately predict PTM sites and variant-induced alterations.Nature Methods2025, 22, 1857–1867. (43) Rives, A.; Meier, J.; Sercu, T.; Goyal, S.; Lin, Z.; Liu, J.; Guo, D.; Ott, M.; Zit- nick, C. L.; Ma, J., et al. Biological structure and function emerge from scaling unsu- pervised learning to 250 million protein sequences.Proceedings of the National Academy of Sciences2021, 118, e2016239118. (44) Lin, Z.; Akin, H.; Rao, R.; Hie, B.; Zhu, Z.; Lu, W.; Smetanin, N.; Verkuil, R.; Ka- beli, O.; Shmueli, Y., et al. Evolutionary-scale prediction of atomic-level protein struc- ture with a language model.Science 2023, 379, 1123–1130. (45) Jumper, J.; Evans, R.; Pritzel, A.; Green, T.; Figurnov, M.; Ronneberger, O.; Tun- yasuvunakool, K.; Bates, R.; Žídek, A.; Potapenko, A., et al. Highly accurate protein structure prediction with AlphaFold.nature 2021, 596, 583–589. (46) Shrestha, P.; Kandel, J.; Tayara, H.; Chong, K. T. Post-translational modification prediction via prompt-based fine-tuning of a GPT-2 model.Nature Communications 2024, 15, 6699. (47) Han, Y.; He, F.; Shao, Q.; Wang, D.; Xu, D. MTPrompt-PTM: A Multi-Task Method for Post-Translational Modification Prediction Using Prompt Tuning on a Structure- Aware Protein Language Model.Biomolecules 2025, 15, 843. 42 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (48) Tariq, A.; Wong, D.; Glab, A.; Iturraran, E.; Kim, I.; Burley, S.; Slochower, D. AstraPTM: Context-Aware PTM Prediction Model for Large-Scale Proteins.bioRxiv 2025, (49) Tariq, A.; Wong, D., et al. AstraPTM2: A Context-Aware Transformer for Broad- Spectrum PTM Prediction.bioRxiv 2025, (50) Liu, L.; Liu, Y.; Wong, K.-C. UniPTM: Multiple PTM site prediction on full-length protein sequence.bioRxiv 2024, (51) Lin, Y., et al. UniPTMs: The First Unified Multi-type PTM Site Prediction Model via Master-Slave Architecture-Based Multi-Stage Fusion Strategy and Hierarchical Con- trastive Loss.arXiv preprint arXiv:2506.054432025, (52) Peng, Z.; Schussheim, B.; Chatterjee, P. PTM-Mamba: a PTM-aware protein language model with bidirectional gated Mamba blocks.Nature Methods2025, (53) Smith, L. M.; Kelleher, N. L. Proteoform: a single term describing protein complexity. Nature methods2013, 10, 186–187. (54) Schwartz,D.; Gygi,S.P.Aniterativestatisticalapproachtotheidentificationofprotein phosphorylation motifs from large-scale data sets.Nature biotechnology2005, 23, 1391– 1398. (55) Radivojac, P.; Vacic, V.; Haynes, C.; Cocklin, R. R.; Mohan, A.; Heyen, J. W.; Goebl, M. G.; Iakoucheva, L. M. Identification, analysis, and prediction of protein ubiquitination sites.Proteins: Structure, Function, and Bioinformatics2010, 78, 365– 380. (56) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; Polosukhin, I. Attention is all you need.Advances in neural information processing systems 2017, 30. 43 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (57) Hu, E. J.; Shen, Y.; Wallis, P.; Allen-Zhu, Z.; Li, Y.; Wang, S.; Wang, L.; Chen, W., et al. Lora: Low-rank adaptation of large language models.ICLR 2022, 1, 3. (58) Ryu, G.-M.; Song, P.; Kim, K.-W.; Oh, K.-S.; Park, K.-J.; Kim, J. H. Genome-wide analysis to predict protein sequence variations that change phosphorylation sites or their corresponding kinases.Nucleic acids research2009, 37, 1297–1307. (59) Tegge, W.; Frank, R.; Hofmann, F.; Dostmann, W. R. Determination of cyclic nucleotide-dependent protein kinase substrate specificity by the use of peptide libraries on cellulose paper.Biochemistry 1995, 34, 10569–10577. (60) Kreegipuu, A.; Blom, N.; Brunak, S. PhosphoBase, a database of phosphorylation sites: release 2.0.Nucleic acids research1999, 27, 237–239. (61) Viengkhou, B.; White, M. Y.; Cordwell, S. J.; Campbell, I. L.; Hofer, M. J. A novel phosphoproteomic landscape evoked in response to type I interferon in the brain and in glial cells.Journal of Neuroinflammation2021, 18, 237. (62) Wansink, D. G.; van Herpen, R. E.; Coerwinkel-Driessen, M. M.; Groenen, P. J.; Hem- mings, B. A.; Wieringa, B. Alternative splicing controls myotonic dystrophy protein kinase structure, enzymatic activity, and subcellular localization.Molecular and cellu- lar biology2003, 23, 5489–5501. (63) Manke, I. A.; Nguyen, A.; Lim, D.; Stewart, M. Q.; Elia, A. E.; Yaffe, M. B. MAPKAP kinase-2 is a cell cycle checkpoint kinase that regulates the G2/M transition and S phase progression in response to UV irradiation.Molecular cell2005, 17, 37–48. (64) Chagoyen, M.; Pazos, F. Quantifying the biological significance of gene ontology biolog- ical processes-implications for the analysis of systems-wide data.Bioinformatics 2010, 26, 378–384. 44 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (65) Chen, M.; Gou, Y.; Lei, M.; Xiao, L.; Zhao, M.; Huang, X.; Liu, D.; Feng, Z.; Peng, D.; Xue, Y. EPSD 2.0: An Updated Database of Protein Phosphorylation Sites Across Eukaryotic Species.Genomics, Proteomics & Bioinformatics2025, 23, qzaf057. (66) Chawla, N. V.; Bowyer, K. W.; Hall, L. O.; Kegelmeyer, W. P. SMOTE: synthetic minority over-sampling technique.Journal of artificial intelligence research2002, 16, 321–357. (67) Lin, T.-Y.; Goyal, P.; Girshick, R.; He, K.; Dollár, P. Focal loss for dense object detec- tion. Proceedings of the IEEE international conference on computer vision. 2017; pp 2980–2988. (68) Ba, J. L.; Kiros, J. R.; Hinton, G. E. Layer normalization. arXiv preprint arXiv:1607.06450 2016, (69) Glorot, X.; Bengio, Y. Understanding the difficulty of training deep feedforward neural networks. Proceedings of the thirteenth international conference on artificial intelligence and statistics. 2010; pp 249–256. (70) Micikevicius, P.; Narang, S.; Alben, J.; Diamos, G.; Elsen, E.; Garcia, D.; Gins- burg, B.; Houston, M.; Kuchaiev, O.; Venkatesh, G., et al. Mixed precision training. arXiv preprint arXiv:1710.037402017, (71) Loshchilov, I.; Hutter, F. Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 2017, (72) Montavon, G.; Orr, G.; Müller, K.-R.Neural networks: tricks of the trade; springer, 2012; Vol. 7700. (73) Srivastava, N.; Hinton, G.; Krizhevsky, A.; Sutskever, I.; Salakhutdinov, R. Dropout: a simple way to prevent neural networks from overfitting.The journal of machine learning research2014, 15, 1929–1958. 45 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (74) Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L., et al. Pytorch: An imperative style, high- performance deep learning library.Advances in neural information processing systems 2019, 32. (75) Wolf, T.; Debut, L.; Sanh, V.; Chaumond, J.; Delangue, C.; Moi, A.; Cistac, P.; Rault, T.; Louf, R.; Funtowicz, M., et al. Transformers: State-of-the-art natural lan- guage processing. Proceedings of the 2020 conference on empirical methods in natural language processing: system demonstrations. 2020; pp 38–45. (76) Chicco, D.; Jurman, G. The Matthews correlation coefficient (MCC) should replace the ROC AUC as the standard metric for assessing binary classification.BioData mining 2023, 16, 4. (77) Tareen, A.; Kinney, J. B. Logomaker: beautiful sequence logos in Python.Bioinfor- matics 2020, 36, 2272–2274. (78) Schneider, T. D.; Stephens, R. M. Sequence logos: a new way to display consensus sequences. Nucleic acids research1990, 18, 6097–6100. (79) Li, W.; Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences.Bioinformatics 2006, 22, 1658–1659. (80) Camacho, C.; Coulouris, G.; Avagyan, V.; Ma, N.; Papadopoulos, J.; Bealer, K.; Mad- den, T. L. BLAST+: architecture and applications. BMC bioinformatics 2009, 10, 421. (81) Altschul, S. F.; Madden, T. L.; Schäffer, A. A.; Zhang, J.; Zhang, Z.; Miller, W.; Lipman, D. J. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs.Nucleic acids research1997, 25, 3389–3402. 46 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint (82) Katoh, K.; Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability.Molecular biology and evolution2013, 30, 772–780. (83) Ashburner, M.; Ball, C. A.; Blake, J. A.; Botstein, D.; Butler, H.; Cherry, J. M.; Davis, A. P.; Dolinski, K.; Dwight, S. S.; Eppig, J. T., et al. Gene ontology: tool for the unification of biology.Nature genetics2000, 25, 25–29. (84) The Gene Ontology resource: enriching a GOld mine.Nucleic acids research2021, 49, D325–D334. (85) Huntley, R. P.; Sawford, T.; Mutowo-Meullenet, P.; Shypitsyna, A.; Bonilla, C.; Mar- tin, M. J.; O’Donovan, C. The GOA database: gene ontology annotation updates for 2015. Nucleic acids research2015, 43, D1057–D1063. (86) Snel, B.; Lehmann, G.; Bork, P.; Huynen, M. A. STRING: a web-server to retrieve and display the repeatedly occurring neighbourhood of a gene.Nucleic acids research2000, 28, 3442–3444. (87) Szklarczyk, D.; Kirsch, R.; Koutrouli, M.; Nastou, K.; Mehryary, F.; Hachilif, R.; Gable, A. L.; Fang, T.; Doncheva, N. T.; Pyysalo, S., et al. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest.Nucleic acids research2023, 51, D638–D646. (88) Newman, M.Networks; Oxford university press, 2018. (89) Watts, D. J.; Strogatz, S. H. Collective dynamics of ‘small-world’ networks.nature 1998, 393, 440–442. 47 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted May 12, 2026. ; https://doi.org/10.64898/2026.05.08.723918doi: bioRxiv preprint