2873 ncomms Nature Communications Nat Commun Nature Publishing Group PMC13494029 13494029 13494029 42469205 10.1038/s41467-026-75283-3 An enzyme-specific protein language model for catalytic property prediction Wang Chong 1 2 3 4 5 Li Mengyao 1 Geng Shaolei 2 5 Li Weidong 2 5 Zhou Xuezhi 3 4 Wang Yu Guang 6 7 Yu Yi 1 3 4 ✉ Wang Tianyun 2 5 ✉ Shen Yiqing 8 ✉ 1 School of Medical Engineering, Henan Medical University, Xinxiang, China 2 Henan Key Laboratory of Neurorestoratology and Protein Modification, The First Affiliated Hospital of Henan Medical University, Xinxiang, China 3 Engineering Technology Research Center of Neurosense and Control of Henan Province, Xinxiang, China 4 Henan International Joint Laboratory of Neural Information Analysis and Drug Intelligent Design, Henan Medical University, Xinxiang, China 5 International Joint Research Laboratory for Recombinant Pharmaceutical Protein Expression System of Henan, School of Basic Medicine, Henan Medical University, Xinxiang, China 6 Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai, China 7 Shanghai AI Laboratory, Shanghai, China 8 Department of Computer Science, Johns Hopkins University, Baltimore, MD USA ✉ Corresponding author. 17 7 2026 17 8760 8760 22 8 2026 © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . Abstract Enzymes drive cellular metabolism, yet predicting catalytic properties from amino acid sequences remains challenging. Existing protein language models (PLMs) provide powerful general-purpose representations but are often inefficient for high-throughput screening and insufficiently adapted to enzyme-specific tasks. Here, we propose EnzGFM, an enzyme-specific PLM based on a Mamba-Transformer hybrid architecture with hierarchical pre-training to capture enzyme-specific patterns. Across enzyme property prediction benchmarks, EnzGFM consistently outperforms Transformer-based PLMs with 2–5-fold acceleration, achieving relative improvements of 16.67% in kinetic parameter prediction, 15.69% in enzyme-reaction mapping, 13.19% in EC number classification, and 20.04% in mutation effect assessment. Building on EnzGFM, we develop EnzGFM-Agent, an enzyme-focused agentic pipeline. Experimental validation further suggests that EnzGFM-Agent can enrich beneficial variants within small candidate pools. Together, these results demonstrate that EnzGFM captures enzyme-specific sequence-function patterns, while EnzGFM-Agent translates these predictions into experimentally actionable candidates and can help reduce wet-lab screening burden for practical enzyme engineering. Subject terms: Computational models, Machine learning Enzyme function prediction from amino acid sequences remains a central challenge in computational biology, despite recent advances in protein language models. This manuscript introduces EnzGFM, an enzyme-specific hybrid model that improves both accuracy and efficiency across multiple prediction tasks and, together with the EnzGFM-Agent pipeline, demonstrates the ability to identify experimentally validated beneficial variants while reducing screening effort. status released display-pdf yes is-olf no is-manuscript no is-preprint no is-journal-matter no is-scanned no is-retracted no Received 2025 Sep 9; Accepted 2026 Jun 22; Collection date 2026. Introduction Enzymes drive biochemical reactions in cellular metabolism and enable the production of pharmaceuticals, biofuels, and chemicals, catalyzing transformations from molecular-scale reactions to industrial manufacturing 1 , 2 . Predicting enzyme catalytic properties from amino acid sequences is central to understanding and engineering function, as the sequence-function relationship determines key characteristics such as catalytic efficiency (e.g., k cat , K m ) and substrate specificity 3 , 4 . The vast sequence space (20 n possible sequences for an n -residue enzyme) makes exhaustive experimental characterization impractical, driving the need for computational approaches 5 . For example, Davidi et al. 6 estimated enzyme functional parameters in Escherichia coli (E. coli) by integrating computationally inferred metabolic fluxes with proteomic measurements across multiple growth conditions, yet could obtain functional parameters for only 436 enzymes even in this well-characterized organism. Deep learning (DL) methods have attempted to address these scalability limitations through supervised learning on enzyme-specific tasks, achieving improved predictions for kinetic parameters and EC classifications 7 – 9 . However, these DL approaches require extensive labeled training data and show limited performance beyond well-annotated enzyme families 10 , creating a need for more generalizable methods that can capture enzyme-specific patterns from sequence. Recent advances in self-supervised protein language models (PLMs) offer a promising framework for addressing this challenge: by pre-training on large-scale unlabeled sequence corpora, PLMs learn evolutionary patterns and structural constraints that support broad transfer to diverse downstream tasks 11 , 12 . Representative models such as ESM 13 , 14 and ProtTrans 12 have demonstrated promising performance across a range of applications, including structure prediction, functional annotation, and variant effect prediction 15 – 18 . Building on this success, existing state-of-the-art methods for enzyme-focused prediction tasks mainly use PLMs to extract representations from sequence, which are then fed into task-specific architectures. For instance, Wang et al. used the pre-trained ProtT5-XL-3B to extract sequence representations, which were then combined with graph neural networks to predict enzyme kinetic parameters 19 . Similarly, Yu et al. employed sequence features extracted by ESM-1B to predict enzyme commission (EC) numbers 20 . Despite these advances, general-purpose PLMs still face several limitations when applied to enzyme-specific prediction tasks. First, existing Transformer-based PLMs show limited computational efficiency, which constrains their application in high-throughput enzyme screening 21 . The quadratic complexity of self-attention mechanisms in these PLMs reduces both training efficiency and inference speed for large-scale enzyme datasets, resulting in processing speeds over 2 × slower than those of linear-complexity models 22 , 23 . Second, general PLMs are trained on heterogeneous protein corpora without adaptation to enzyme-specific data and may fail to learn representations that capture the unique sequence patterns associated with enzymatic function 24 – 26 . Together, these limitations underscore the importance of enzyme-domain adaptation, for example, through continued pretraining on enzyme-focused corpora 27 . To address these limitations, we present EnzGFM, an enzyme-specific PLM for enzyme property prediction built upon hierarchical pre-training and an efficient long-context backbone. We further develop EnzGFM-Agent, a dedicated multi-agent system for enzymes that enables automated code generation, analysis, and decision support for enzyme mutagenesis, retrieval, and discovery through a user-friendly interactive interface. Hybrid SSM-attention architectures have shown favorable accuracy-efficiency trade-offs for long biological sequences and have been adopted in biological foundation models, including Evo and related protein sequence models 28 – 31 . Inspired by these advances, EnzGFM adopts a Mamba-Transformer 23 , 32 hybrid architecture, combining Mamba layers for efficient long-range context modeling with Transformer layers for residue-level dependency and catalytic motif learning. EnzGFM is trained through a two-stage hierarchical pre-training strategy, first on general protein sequences and then on enzyme sequences, enabling adaptation from general protein representations to enzyme-specific features. Across four enzyme-specific prediction benchmarks, EnzGFM consistently improves over representative Transformer-based PLMs while achieving 2–5-fold acceleration, with relative gains in kinetic parameter prediction, enzyme-reaction mapping, EC number classification, and mutation effect assessment. We further evaluated EnzGFM-Agent in an experimental enzyme mutagenesis workflow using secreted alkaline phosphatase (SEAP), where all generated single-point mutants showed statistically significant activity improvements over the original sequence, with the best variant reaching a 7.13-fold increase in activity. Together, these results establish EnzGFM as an efficient and accurate enzyme-specific foundation model, and EnzGFM-Agent as a practical multi-agent system for translating enzyme-specific predictions into actionable workflows for enzyme analysis and engineering. Results Evaluation of EnzGFM pre-training performance EnzGFM addresses the computational and representational limitations of existing general-purpose PLMs through a hybrid architecture that combines Mamba structured state space models (SSMs) 23 with Transformer blocks 32 (Fig. 1 a). The Mamba processes long-range dependencies with linear-time complexity O ( n ) for sequence length n , while Transformer layers capture local catalytic motifs through self-attention mechanisms. The EnzGFM implements hierarchical pre-training that progressively adapts from general protein patterns to enzyme-specific features. We first pre-train on 48.28 million protein sequences (16.67 billion tokens) from UniRef50 33 using masked language modeling, establishing foundational protein language representations. Enzyme-specific pre-training then specializes the EnzGFM on 5.09 million enzyme sequences (1.96 billion tokens) curated from CAZy 34 , BRENDA 35 , and SABIO-RK 37 . This enzyme dataset maintains natural taxonomic distributions across all seven EC classes, with sequence lengths spanning 50 to 5000 amino acids (Fig. 1 b, Supplementary Fig. 1 ). Fig. 1 EnzGFM architecture and hierarchical pre-training strategy. a Two-stage training pipeline employing hybrid Mamba-Transformer architecture. General pre-training establishes foundational protein language understanding using protein sequences, followed by enzyme-specific pre-training that adapts representations using enzyme sequences. Architecture combines Mamba structured state space models for efficient long-range dependency modeling with Transformer blocks for local catalytic motif recognition. Four downstream tasks are used to evaluate enzyme property prediction. b Sequence length distributions across both training stages. c Scaling analysis demonstrating perplexity reduction with model size, from 30M to 1.5B parameters. d t-SNE visualization of enzyme embeddings colored by EC class. Left: general protein pre-training. Right: enzyme-specific continued pre-training. Clustering quality is assessed by Silhouette coefficient, where higher is better; Davies-Bouldin index, where lower is better; and Calinski-Harabasz index, where higher is better. EnzGFM performance improved systematically with model scale. Across models ranging from 30 million to 1.5 billion parameters, training loss decreased from 2.0881 in EnzGFM-30M to 0.9230 in EnzGFM-1.5B (Supplementary Fig. 2 a, c). Perplexity followed a power-law relationship with model size, with size-scaling exponents of α size = 0.3682 for general pre-training and α size = 0.2935 after enzyme-specific pre-training (Fig. 1 c). At matched parameter counts, enzyme-specific pre-training reduced perplexity by 5–35%, consistent with improved modeling efficiency after adaptation to enzyme sequences. Similar power-law behavior was observed with respect to computational cost, with larger models requiring proportionally less additional compute to achieve equivalent perplexity reductions (Supplementary Fig. 2 c, d). Computational efficiency measurements demonstrate that EnzGFM variants achieve 2–5 × inference acceleration compared with Transformer-based PLMs (ESM2 13 , ESM-C 36 , and ProtTrans 12 ) (Supplementary Fig. 3 a). EnzGFM-650M processes sequences in 12.42 ms on average, running 2.4 × faster than ESM2-650M (29.94 ms). Analysis of learned representations through dimensionality reduction reveals improved functional organization after enzyme-specific pre-training (Fig. 1 d). Representations from general pre-training show limited separation between enzyme classes, as reflected by low Silhouette coefficient 38 , high Davies-Bouldin index 39 , and low Calinski-Harabasz index 40 . After enzyme-specific pre-training, all three metrics improve consistently, indicating enhanced functional separability in the learned representation space, although overlap at the EC class level remains expected given shared protein folds and domain architectures, as well as the broad heterogeneity within each EC class. These functionally organized representations establish the foundation for improved performance across enzyme-specific prediction tasks. Kinetic parameter prediction performance evaluation Enzyme kinetic parameters provide quantitative measures of catalytic function that inform biocatalyst design and metabolic engineering strategies 41 . We evaluated EnzGFM’s capacity to predict three core kinetic parameters, namely the turnover number ( k cat ), Michaelis constant ( K m ), and catalytic efficiency ( k cat / K m ), which characterize maximum reaction velocity, substrate binding affinity, and overall catalytic efficiency, respectively 42 . Our evaluation employed the UniKP framework 3 , which integrates protein sequence embeddings with molecular substrate representations through SMILES Transformer encoding (Fig. 2 a). Four EnzGFM variants (30M, 150M, 650M, and 1.5B parameters) were benchmarked against state-of-the-art PLMs including ESM2 13 (35M, 150M, 650M, and 3B), ESM-C 36 (300M and 600M), and ProtTrans models 12 (ProtAlbert-224M, ProtXLNet-409M, and ProtT5-XL-3B). We assessed performance on both biased datasets with random train-test splits from UniKP 3 (16,838 k cat and 11,722 K m samples) and unbiased datasets constructed through CD-HIT 43 clustering at 40% sequence similarity 10 (27,658 k cat , 42,018 K m , and 25,831 k cat / K m samples). All experiments employed tenfold cross-validation. Fig. 2 Enzyme kinetic parameter prediction performance evaluation. a UniKP framework architecture integrating enzyme sequence embeddings from PLMs with substrate molecular representations via SMILES Transformer, SM-Trans. Concatenated features undergo ExtraTree regression to predict kinetic parameters. b , e , h Performance comparison between unbiased and biased evaluations on UniKP for k cat , K m , and k cat / K m , respectively. The x-axis reports R 2 on the unbiased test set, and the y-axis reports R 2 under random train--test splits, representing the biased setting. Colors denote models, with the same model using the same color. Marker size indicates model capacity, measured by the number of parameters. c , f , i Performance on the unbiased split for k cat , K m , and k cat / K m , respectively. Boxplots show PCC distributions from tenfold cross-validation across model-size bins ( n = 10), including < 35M, 150–300M, 400–650M, and > 650M. The center line indicates the median; box boundaries indicate the interquartile range (IQR), corresponding to the 25th and 75th percentiles. Whiskers indicate the minimum and maximum non-outlier values, defined as the smallest and largest observations within 1.5 × IQR of the lower and upper quartiles, respectively. Individual points indicate outliers beyond this range. Colors denote different PLMs. d , g , j Experimental versus predicted correlations on the unbiased split for k cat , K m , and k cat / K m , respectively, shown as density distributions. Source data are provided as a Source Data file. EnzGFM demonstrated consistent improvements across all kinetic parameters on both dataset types (Supplementary Tables 1 – 3 ). For turnover number prediction on biased datasets, EnzGFM-1.5B achieved R 2 = 0.6646, matching the performance of ProtT5-XL-3B ( R 2 = 0.6644) (Fig. 2 b, Supplementary Fig. 4 a). In K m prediction, EnzGFM-650M reached R 2 = 0.5086, exceeding the leading baseline ESM2-3B ( R 2 = 0.5023) by 1.25% (Fig. 2 e, Supplementary Fig. 4 b). For catalytic efficiency prediction, EnzGFM-1.5B attained R 2 = 0.5526, surpassing ESM2-3B ( R 2 = 0.5481) (Fig. 2 h, Supplementary Fig. 4 c). The performance advantages became more pronounced on unbiased datasets, which represent more realistic evaluation conditions with evolutionarily distant test sequences. EnzGFM-150M achieved the best performance across all parameters with R 2 = 0.2073 for k cat prediction (16.67% improvement over ESM2-150M), PCC = 0.6211 for K m prediction (2.15% improvement), and PCC = 0.3779 for catalytic efficiency (6.15% improvement over ESM2-150M) (Fig. 2 c, f, i). We further assessed potential pre-training evaluation leakage by computing the maximum sequence identity between each benchmark protein and the enzyme-specific pre-training corpus using MMseqs2. We stratified benchmark sequence into identity bins and evaluated k cat / K m prediction performance within each bin on unbiased dataset (Supplementary Fig. 5 a–b). EnzGFM retained competitive performance even in the lowest identity range (0.0–0.2), and the log-error distributions did not differ significantly across identity bins (Kruskal–Wallis p = 0.4736), suggesting that the model learns transferable representations rather than relying on sequence similarity. We also examined whether the naturally imbalanced taxonomic composition of the pre-training data leads to organism-dependent degradation. Stratifying the k cat benchmark by organism groups showed comparable correlation metrics across groups, and the per-sample log-error distributions did not differ significantly (Kruskal–Wallis p = 0.0579; Supplementary Fig. 5 c, d). To quantify the contribution of the hybrid architecture, we conducted an ablation on the unbiased split (Supplementary Table 4 ), where the hybrid Mamba–Transformer outperformed Mamba-only and Transformer-only counterparts across k cat , K m , and k cat / K m , with the largest relative gains in SCC for k cat (19.13% vs. Transformer-only), R 2 for K m (11.69%), and R 2 for k cat / K m (49.06%). To evaluate the impact of substituting task-specific PLM backbones with EnzGFM, we kept the downstream architectures of TurNuP ( k cat ) 44 , UniKP ( K m ) 3 , CataPro ( k cat / K m ) 10 , and CatPred ( K i ) 41 unchanged and replaced only the sequence encoder with EnzGFM (Supplementary Fig. 6 ). In k cat prediction with TurNuP, EnzGFM improved performance, with SCC increasing from 0.6387 to 0.6569. For K m prediction, UniKP with EnzGFM achieved a higher R 2 (0.5065 vs. 0.4958) than the original backbone. In the k cat / K m task with CataPro, PCC increased from 0.5572 to 0.6394 (a 14.75% improvement). For K i prediction with CatPred, replacing the backbone with EnzGFM improved R 2 from 0.6299 to 0.6382. Overall, these comparisons indicate that EnzGFM serves as an effective enzyme-pretrained feature extractor, improving performance when used as a drop-in replacement for task-specific PLM backbones across all diverse kinetic prediction tasks. Model scaling analysis revealed a split-dependent effect of model size. Under biased splits, where test sequences were more similar to the training distribution, larger PLMs achieved the best performance for two of the three kinetic endpoints. In contrast, under low-similarity, unbiased splits, intermediate-sized EnzGFM models performed best across all kinetic parameters. A similar trend was observed for general PLMs, where ESM2-35M outperformed ESM2-3B by 2.64–11.60% on unbiased k cat prediction. This pattern is consistent with prior observations 45 that medium-sized PLMs can generalize better in low-similarity transfer settings, potentially because they are less sensitive to training-set-specific regularities under limited data. Enzyme-reaction retrieval performance evaluation Accurate prediction of enzyme-reaction correspondences drives computational enzyme engineering and metabolic pathway design 46 . Therefore, we evaluated EnzGFM’s ability to capture biocatalytic specificity using the ReactZyme framework 47 , which employs cross-attention mechanisms to align enzyme sequence embeddings with reaction prediction (Fig. 3 a). We evaluated two complementary retrieval tasks, namely enzyme-to-reaction prediction (ranking candidate reactions for a given enzyme sequence) and reaction-to-enzyme prediction (ranking candidate enzyme sequences for a given reaction). Given prior results that EnzGFM outperforms parameter-matched PLM counterparts on enzyme kinetic parameter prediction, we benchmarked four EnzGFM variants (30M, 150M, 650M, 1.5B) against three large-scale PLMs (ESMC-600M 36 , ProtT5-XL-3B 12 , ESM2-3B 13 ). The evaluation employed the ReactZyme dataset 47 containing 178,463 enzyme-reaction pairs across 178,327 enzymes and 7,726 reactions. To ensure generalization assessment, we enforced ≥60% Levenshtein sequence difference 48 between training and test enzymes using a 95%/5% data split. Fig. 3 Enzyme-reaction retrieval prediction across PLMs. a ReactZyme architecture for bidirectional enzyme--reaction prediction. Reaction SMILES and enzyme sequences are encoded into reaction and enzyme embeddings, respectively, followed by cross-attention; the confusion matrix shows predicted probabilities for enzyme--reaction pairs. b MAP scores comparing PLM performance on enzyme-to-reaction prediction tasks. Bars show mean MAP ± SD from five independent experiments ( n = 5); and overlaid white circles indicate individual experimental values. Statistical significance between the best- and second-best-performing models was assessed using a one-tailed paired t -test. The exact p value was p = 0.0004. Full statistics are in Supplementary Table 6 . c Top- K retrieval accuracy curves for enzyme-to-reaction prediction, with K = 1, 5, 10, 20, 50. Lines represent mean accuracy from five independent experiments ( n = 5), and shaded regions indicate ± SD. The color legend shown in panel ( c ) applies to both panels ( c ) and ( d ), with the same colors denoting the same PLMs. d Top- K NDCG scores for reaction retrieval across different PLMs, with K = 1, 5, 10, 20, 50. Bars represent mean NDCG scores from five independent experiments ( n = 5); error bars indicate SD; and overlaid black circles indicate individual experimental values. e Correlation analysis between model performance and physicochemical protein properties. f MAP versus training data size, from 1M to 28M samples, for EnzGFM variants. Solid lines show mean MAP ( n = 5); shaded regions indicate ± SD. g MAP versus model size, from 30M to 1.5B parameters, under different training data regimes. Dashed curves in ( f ) and ( g ) denote power-law fits; exponents and R 2 are shown in ( g ). *** denotes p < 0.001. Source data are provided as a Source Data file. EnzGFM demonstrated superior performance in enzyme-to-reaction prediction across all model sizes (Fig. 3 b and Supplementary Table 6 ). EnzGFM-1.5B achieved MAP of 0.5156, exceeding ProtT5-XL-3B (MAP of 0.4457) by 15.69% and ESM2-3B (MAP of 0.4785) by 7.75% ( p < 0.001). Top-K retrieval analysis corroborated these improvements (Fig. 3 c, d). EnzGFM-1.5B achieved Top-5 accuracy of 0.6636 (7.10% increase over ESM2-3B’s 0.6196) and NDCG of 0.5152 (9.11% increase over ESM2-3B’s 0.4722). The EnzGFM-650M variant maintained competitive retrieval performance with Top-5 accuracy of 0.6403 and NDCG of 0.5001. The reciprocal reaction-to-enzyme task yielded concordant results (Supplementary Fig. 7 ). EnzGFM-1.5B attained MAP of 0.8211, outperforming the leading baseline ESM2-3B (MAP of 0.7897) by 3.98%. The intermediate EnzGFM-650M model (MAP of 0.8159) also exceeded ProtT5-XL-3B performance, supporting that enzyme-specific pre-training improves reaction-enzyme retrieval performance. Robustness analysis across molecular and biochemical properties revealed consistent performance (Fig. 3 e). EnzGFM maintained high MAP values between 0.7873 and 0.8651 across bins of protein length, molecular weight, hydrophobicity, net charge, and aromatic residue content, with no pronounced degradation at distribution extremes. Correlation analysis showed positive associations between performance and substrate chemical similarity, with correlation coefficients 0.2806–0.5239 (Supplementary Fig. 8 ), suggesting that learned representations encode chemically meaningful enzyme-substrate relationships. Performance was consistent across all seven enzyme classes (Supplementary Fig. 9 ), including EC7, despite its low frequency in the pre-training corpus, suggesting no systematic degradation for low-frequency classes in reaction prediction. We further investigated the scaling relationship between fine-tuning dataset size and performance (Fig. 3 f, g) to reveal the distinct performance trends with model size and training data volume. Training set growth followed approximate power-law improvement with exponents α = 0.06–0.10 across model variants. At the maximum data scale (28M samples), scaling the model from 30M to 1.5B parameters increased enzyme-to-reaction MAP from 0.3510 to 0.5156 and reaction-to-enzyme MAP from 0.7001 to 0.8211. In contrast, under limited supervision (5% training data; 1.4M samples), the absolute gains were smaller: enzyme-to-reaction MAP increased from 0.0941 to 0.1351 and reaction-to-enzyme MAP from 0.3230 to 0.3712. This pattern indicates that larger PLMs require sufficient training data to realize their representational capacity. The observed scaling relationships and performance correlates provide practical guidance for deploying enzyme-specific models across diverse biochemical prediction scenarios. EC number classification performance evaluation Enzyme classification through EC number prediction enables functional annotation of uncharacterized proteins and informs biocatalyst discovery efforts 49 . Consequently, we evaluated EnzGFM’s classification performance using the CLEAN framework 20 , which distinguish enzymes with identical EC numbers using the embedding from PLM by maximizing separation between different functional classes (Fig. 4 a). We benchmarked EnzGFM variants against state-of-the-art PLMs, including ESM-1B (the original CLEAN backbone) 50 , ESMC-600M 36 , ProtT5-XL-3B 12 , and ESM2-3B 13 . The evaluation employed a curated dataset from Yu et al. 20 containing 220,000 filtered SwissProt sequences with complete EC annotations. Performance was assessed on independent test datasets NEW-392 and Price-149 51 . Fig. 4 EC number prediction. a CLEAN architecture for enzyme EC classification. EnzGFM embeddings are optimized by contrastive learning to cluster enzymes with identical EC numbers and separate distinct classes. b Precision scores across sequence identity thresholds, from 10 to 100%. c Precision and F1-score comparisons. Each dot represents one independent experimental replicate, and horizontal lines and error bars indicate mean ± SD from five independent experiments, n = 5. Statistical significance between the best- and second-best-performing models was assessed using a one-tailed paired t -test. The exact p values were p = 0.0093 for Precision and p = 0.0260 for F1-score. df = 4, 95% CI, t , and Cohen’s d are reported in Supplementary Table 7 . The color legend shown in panel ( c ) applies to panels ( b )–( d ), with the same colors denoting the same PLMs. d Hierarchical classification accuracy across EC classification levels, from EC X.--- through EC X.X.X.X. e F1-score progression across sequence identity thresholds with corresponding training dataset sizes. Performance curves illustrate model behavior under varying evolutionary distance constraints between training and test sequences. f t-SNE embedding visualization, with points colored by EC number subclass membership; gray points denote other ECs. g Inter-category classification performance matrix across EC classes, with color intensity indicating prediction accuracy between different enzyme functional categories. * denotes p < 0.05, ** denotes p < 0.01 relative to baseline models. Source data are provided as a Source Data file. EnzGFM showed improved EC classification performance across sequence identity thresholds (Fig. 4 b). At 100% sequence identity, EnzGFM-1.5B achieved precision of 0.6358, outperforming ESM2-3B (0.6159) and ProtT5-XL-3B (0.6084). The smaller EnzGFM-650M variant also exceeded all baseline models with precision of 0.6256, representing a 1.57% improvement over ESM2-3B despite having fewer parameters. Performance advantages became more pronounced at lower sequence identity thresholds, indicating superior generalization to evolutionarily distant sequences. At 50% sequence identity, EnzGFM-1.5B achieved mean precision of 0.5424 ( p < 0.001) and F1-score of 0.4534 ( p < 0.05) on both test datasets (Fig. 4 c and Supplementary Fig. 10 ). The improvement was most substantial at the stringent 10% identity threshold, where EnzGFM achieved precision of 0.3732 compared to ESM2-3B’s precision of 0.3297, representing a 13.19% increase. In terms of accuracy across EC hierarchy levels, EnzGFM-1.5B and ProtT5-XL-3B showed similar overall performance (Fig. 4 d). Together, these results show that enzyme-specific pre-training improves EC classification robustness under stringent low-identity evaluation, supporting better functional generalization beyond close homologs. We further examined model robustness under progressive sequence similarity constraints between training and test sets. Larger EnzGFM variants maintained robust performance even under substantial distributional shift at 30% identity, with EnzGFM-650M attaining F1 of 0.6390 and EnzGFM-1.5B reaching F1 of 0.6330 (Fig. 4 e). Increasing the identity threshold to 50% yielded performance gains across all model scales, with EnzGFM-1.5B achieved F1 of 0.8770, EnzGFM-650M reached F1 of 0.8822, and EnzGFM-150M attained F1 of 0.8188. At 100% identity, EnzGFM achieved F1 scores of 0.9531, 0.9416, and 0.9314, respectively, demonstrating a narrowing performance gap between model variants at higher sequence identities. Embedding analysis revealed that EnzGFM captures fine-grained functional relationships beyond sequence similarity patterns. t-SNE visualization of EnzGFM-650M embeddings showed that enzymes from specific EC subclasses formed compact and distinct clusters (Fig. 4 f). Entries from EC 3.1.2.22, EC 4.2.1.113, and EC 3.5.4.38 each exhibited high intra-cluster density with clear separation from other subclasses. This clustering pattern indicates that EnzGFM captures detailed and biologically meaningful functional relationships at the subclass level rather than relying solely on overall sequence similarity. Performance variations across EC categories reflected the inherent structural and catalytic diversity within enzyme families (Fig. 4 g). Analysis across individual EC categories further shows that EnzGFM maintains stable predictive performance even for classes with low representation in the pre-training corpus, suggesting that enzyme-class imbalance does not materially affect EC number prediction (Supplementary Table 8 ). These analyses demonstrate that EnzGFM learns hierarchical enzyme representations encoding fine-grained functional specificity, capturing both subtle distinctions between closely related enzyme subclasses and broader catalytic diversity across enzyme families. Single amino acid substitution effect prediction performance Predicting the functional consequences of single amino acid substitutions drives rational enzyme design and protein engineering applications 18 . We therefore evaluated EnzGFM’s capacity to assess mutation effects using a curated dataset of 225 enzyme mutational studies from literature and public databases, comprising 16,607 experimentally characterized single substitutions after filtering for datasets with at least 5 mutations per enzyme. Mutation effects were quantified as log conditional probability differences between mutant and wild-type sequences, where higher scores indicate greater contextual preference for the substitution. Following the ESM-1V 18 protocol, which scores substitutions by masked-token likelihood differences, we compared the MLM-pretrained models EnzGFM-650M/1.5B with ESM-1V and ESM2-3B 13 using correlation-based metrics across the complete mutation dataset. Across the mutation-effect datasets, EnzGFM-1.5B attained the highest mean correlations, with a PCC of 0.2703, SCC of 0.2695, and Kendall’s τ of 0.2044, compared with 0.2384, 0.2245, and 0.1675, respectively, for ESM-1V, and 0.2625, 0.2497, and 0.1870 for ESM2-3B (Fig. 5 a and Supplementary Table 9 ). This represents a 20.04% relative improvement in SCC compared to ESM-1V and a 7.92% improvement over ESM2-3B. Despite substantial across-dataset variability (wide confidence intervals, likely reflecting heterogeneity across mutational scans such as sequence context and assay design), paired Wilcoxon signed-rank tests supported EnzGFM-1.5B’s improvement over ESM-1V for rank-based metrics (SCC p = 0.0150; Kendall’s τ p = 0.0155), and EnzGFM-1.5B also achieved higher mean correlation metrics than ESM2-3B. Analysis across mutation count strata revealed consistent model performance improvements as dataset sizes increased from 0–50 to >500 mutations per enzyme (Fig. 5 b). Fig. 5 Single amino acid substitution effect prediction across PLMs. a Correlation coefficients, including SCC, PCC, and Kendall, between predicted and experimental mutation effects across benchmark datasets. Violin plots show data distributions across benchmark datasets, n = 225. Embedded box plots show the median, center line; interquartile range, IQR, box limits; and 1.5 × IQR, whiskers. Dots represent individual datasets. Black circles and error bars indicate mean ± 95% confidence interval. b SCC scaling with mutation dataset size. Points or lines indicate mean SCC within each mutation-count bin, and shaded regions indicate 95% CI across benchmark datasets in each bin. Functional parameter category analysis comparing prediction accuracy for ( c ) substrate affinity parameters ( n = 11), d catalytic activity parameters ( n = 87), and e protein stability parameters ( n = 127). Box plots display correlation distributions across experimental datasets within each functional category. For each box plot, the center line indicates the median; box limits indicate the IQR (25th-75th percentiles); whiskers extend to 1.5 × IQR; and individual data points are overlaid with jitter. f SCC distributions across 25 physicochemical amino acid substitution classes. g–i Case studies on three enzyme targets, 4RZH, 8CIP, and 1H0C, comprising 6, 4, and 9 variants, respectively, with duplicate positions excluded. For each target, the left panel shows the structural mapping of experimentally measured mutation effects, the right panel presents the corresponding model predictions, and the middle panel provides a magnified view of the local environment. Mutation sites are indicated as colored spheres, with color gradients representing normalized mutation effect scores ranging from − 1 to + 1. Source data are provided as a Source Data file. We examined model performance across three distinct functional parameter categories to assess specificity for different enzymatic properties. For substrate affinity parameters, EnzGFM-1.5B achieved SCC of 0.3832, surpassing ESM-1V (0.1325) and ESM2-3B (0.1341) by 189.20% and 185.75%, respectively (Fig. 5 c). In comparison, for catalytic activity parameters, the prediction performance is slightly lower than that for substrate affinity parameters. In catalytic activity prediction (including k cat , V max , and activity), EnzGFM-1.5B reached SCC of 0.2993, exceeding ESM-1V (0.2927) by 2.25% and ESM2-3B (0.2686) by 11.43% (Fig. 5 d). This discrepancy may be attributed to differences in the biochemical determinants and data characteristics of these endpoints: substrate-affinity–related measurements are typically more directly linked to recognition features, whereas catalytic activity are more sensitive to assay conditions and reporting practices, which can introduce higher annotation heterogeneity across databases. For protein stability parameters, EnzGFM attained SCC of 0.2392, showing a 28.81% improvement over ESM-1V (0.1857) but falling 3.08% below ESM2-3B (0.2468) (Fig. 5 e). This relatively lower performance on stability-related prediction is consistent with prior reports that stability effects rely more on global structural context and long-range residue couplings, where structure-aware models typically outperform sequence-only PLMs 52 . To further examine task-dependent model behavior, we analyzed residue-level attention distributions in mutation effect prediction (Supplementary Fig. 12 ). We computed the mean attention score on annotated active-site residues and on all remaining residues and report their difference (Δattention). The active-site–vs–non–active-site attention difference is smallest for protein stability (mean Δattention = 0.1122) and largest for substrate affinity (mean Δattention = 0.1817), with catalytic activity in between (mean Δattention = 0.1578), suggesting that the active-site–vs–non–active-site attention contrast differs across parameter category. We also assessed whether self-attention highlights functionally important regions by comparing, on a per-protein basis, the mean attention scores of UniProt-annotated catalytic, active-site, and binding-site residues against all remaining positions. As shown in Supplementary Fig. 13 , all three annotated residue types have significantly higher attention scores than other residues (Wilcoxon signed-rank test). This enrichment of attention at annotated functional sites suggests that self-attention captures biologically meaningful signals. Across 25 physicochemical source to target substitution categories, the model showed higher average SCC for inter category than within category substitutions (0.3398 versus 0.2657), with strongest performance for charge or structure altering transitions (Fig. 5 f). Conversely, within category substitutions showed reduced correlations, likely due to subtle and context dependent effects that challenge monotonic ranking. Representative validation on three enzyme targets demonstrated strong experimental agreement, with correlation coefficients ranging from 0.7477 to 0.9502 (Fig. 5 g–i). As an external validation on a widely used DMS benchmark, we also evaluated mutation-effect prediction on ProteinGym 53 . We selected 75 enzyme sequences from the 217 proteins in ProteinGym by mapping each ProteinGym protein to its UniProt entry and retaining those with an assigned EC number in UniProt. On this enzyme subset, EnzGFM-1.5B achieved the best overall performance among the evaluated models (SCC = 0.4630; Supplementary Fig. 11 ), improving over ESM-1V (SCC = 0.4279) and being slightly higher than ESM2-3B (SCC = 0.4577). Experimental evaluation of enzyme mutagenesis with EnzGFM-Agent To facilitate enzyme design and analysis, we developed a dedicated multi-agent system for enzymes (EnzGFM-Agent), enabling automated code generation and analysis for enzyme mutagenesis, retrieval, and discovery through a user-friendly interactive interface (Fig. 6 a, Supplementary Fig. 14 ). We selected the SEAP reporter gene to evaluate the performance of EnzGFM in single-point mutation prediction and validated the results in Chinese hamster ovary (CHO) cells. CHO cells were chosen because of their rapid growth, ease of culture, low endogenous protein secretion, and human-like post-translational modification capacity 54 , 55 . Under normal circumstances, eukaryotic cells do not naturally produce SEAP protein, and SEAP activity cannot be detected in eukaryotic cells that have not been transfected with the SEAP gene, thereby improving the accuracy and reliability of the assay. Furthermore, compared with E. coli and yeast expression systems, recombinant therapeutic proteins (RTPs) expressed in CHO cells more closely resemble their native counterparts in molecular characteristics, three-dimensional conformation, and biological activity. Thus, the CHO-SEAP platform represents an appropriate system for functional screening of enzyme mutants and for establishing activity evaluation models of RTPs. Fig. 6 EnzGFM-Agent architecture and experimental validation of SEAP single-point mutants. a Schematic of the multi-agent system. User inputs, including bioinformatic files or natural language queries, are handled by the Orchestrator, which routes tasks to specialized agents: Chat Agent for text responses, Data Agent for enzyme information extraction, Task Agent for task parsing, and downstream code execution agents, including Code, Execute, and Analyst. b–d Measured activities of SEAP single-point mutants predicted by EnzGFM and ESM2. (b) Mutants predicted by both models. c Mutants uniquely predicted by EnzGFM. d Mutants uniquely predicted by ESM2. Activities are normalized to the wild-type (WT) mean measured in the same experimental setting, which is set to 1.0 and indicated by dashed lines. Bars represent mean ± SEM from three independent biological replicates ( n = 3), with individual data points overlaid. Statistical significance compared to WT was assessed by one-way ANOVA and is indicated above each bar: * p < 0.05, ** p < 0.01, *** p < 0.001. Exact p values and associated statistics for all mutant-versus-WT comparisons are provided in the Source Data file for ( b – d ). Using EnzGFM-Agent, all possible single-point mutations of SEAP were scored, and the top 10 candidates from EnzGFM-1.5B and ESM2-3B were selected, respectively, yielding a total of 15 unique variants, including 5 shared by both models. The shared variants exhibited significantly higher activity than the wild-type (WT) sequence (ANOVA, p < 0.05) (Fig. 6 b). Among the five variants uniquely prioritized by EnzGFM, three exceeded the commonly used directed-evolution hit threshold of a 1.5-fold increase over WT, with A46I reaching 6.52-fold activity. In addition, the remaining two EnzGFM-specific variants also showed significantly higher activity than WT ( p < 0.05), indicating that all five EnzGFM specific candidates were experimentally beneficial (Fig. 6 c). In contrast, among the five variants uniquely prioritized by ESM2, only two showed significantly higher activity than WT, whereas the other three exhibited lower activity than WT (Fig. 6 d, Supplementary Fig. 15 ). Although the EnzGFM prioritized candidate set showed a higher overall success rate in identifying variants with experimentally validated activity improvement, the absolute number of high-threshold hits was not substantially larger under stricter criteria. This finding highlights the advantage of the EnzGFM-Agent pipeline in candidate prioritization, as it integrates mutation scoring, ranking, and selection to translate model predictions into experimentally actionable candidates. By enriching beneficial mutations within a small candidate pool, this workflow can reduce the burden of wet-lab screening. Discussion General-purpose PLMs face two constraints when applied to enzyme analysis: (1) quadratic computational complexity that limits large-scale screening applications, and (2) training on heterogeneous protein corpora that dilutes enzyme-specific patterns. EnzGFM addresses both limitations through architectural innovation and targeted training strategies, while EnzGFM-Agent extends these predictions toward practical enzyme engineering workflows. The hybrid Mamba-Transformer architecture of EnzGFM reduces computational complexity from O ( n 2 ) to approximately O ( n ) behavior while preserving the representational capacity needed for enzyme function prediction. Mamba’s structured SSMs 23 can process long-range sequence dependencies efficiently, while Transformer blocks 32 capture local catalytic motifs through self-attention mechanisms. This design achieves 2 − 5 × faster inference compared to Transformer-based PLM without sacrificing predictive accuracy. The hierarchical pre-training paradigm first captures general protein knowledge from large-scale sequence data and then refines these representations using curated enzyme-specific sequences, thereby improving representation quality and functional clustering. Ablation experiments further indicate that these gains are not attributable to enzyme-specific training alone: when ESM2 and ProtT5 were continued pre-trained using the same datasets and protocols, EnzGFM still achieved higher accuracy across all kinetic parameter prediction tasks, with R 2 improvements over continued-pretrained ESM2 of 11.03%, 2.77%, and 9.91% for k cat , K m , and k cat / K m , respectively (Supplementary Table 5 ). Across the evaluated downstream tasks, EnzGFM generally achieved competitive or improved performance relative to baseline models, although the strength of statistical support varied across tasks and evaluation metrics. The clearest and most consistent statistically supported gains were observed in enzyme reaction retrieval, suggesting that representations derived from EnzGFM may better encode information relevant to enzyme functional relationships and associations between sequences and reactions. However, the improvement in Recall for EC number prediction was not statistically significant. These results suggest that general-purpose protein language models can remain competitive on certain tasks and metrics. EnzGFM nevertheless provided a more favorable balance between cross-task performance and computational efficiency, supporting its potential utility for large-scale enzyme function prediction, biocatalyst screening, and enzyme reaction matching. EnzGFM achieved improvements in predicting substrate affinity mutation effects, indicating enhanced sensitivity to substitutions affecting catalytic function, whereas its more modest gains in protein stability prediction suggest that different PLMs retain task-specific strengths and limitations. For practical enzyme engineering, mutation effect predictions must be translated into experimentally testable designs by integrating multiple predictive signals and experimental constraints. EnzGFM-Agent was designed to address this gap by integrating EnzGFM-derived predictions with complementary model outputs, including ESM2, within a unified decision framework to support mutation analysis and experimental prioritization. In the SEAP mutagenesis workflow, candidates prioritized by this integrated pipeline were enriched for experimentally beneficial variants, demonstrating that EnzGFM-Agent can translate sequence-level predictions into actionable mutagenesis decisions within a small validation workload. The observed scaling relationships provide practical guidance for model deployment. In our benchmarks, under unbiased splits designed to yield low train-test sequence similarity, intermediate-sized models (150M parameters) outperformed larger variants, whereas under random splits, larger models generally performed best. This pattern is consistent with prior observations 45 that medium-sized PLMs can transfer better under realistic low-similarity settings. We hypothesize that this interaction between model size and split difficulty reflects a data-capacity trade-off: as the task becomes more out-of-distribution, larger models may require substantially more diverse training data and/or stronger regularization for their additional capacity to translate into improved generalization, rather than fitting training-set-specific regularities. Accordingly, model selection should jointly consider (i) expected train-test similarity (task difficulty), (ii) the amount and diversity of labeled training data available for fine-tuning, and (iii) the target generalization regime (in-distribution versus low-similarity transfer). Several limitations warrant consideration for future development. The curated training datasets may introduce taxonomic or functional biases, as evidenced by underrepresentation of translocases (EC 7, 0.13%). Although our downstream evaluations did not reveal clear systematic degradation in performance for low-frequency classes, such class imbalance may still limit generalization to rare enzyme families or emerging catalytic mechanisms, particularly in settings or tasks that beyond those considered here. Future research can explore expanding taxonomic and functional coverage through active learning strategies that explicitly identify underrepresented enzyme families for targeted data collection or model fine tuning. Imbalance aware pretraining schemes, such as reweighting, targeted sampling, or data augmentation for rare enzyme classes, also represent a promising direction, especially for applications focused on low abundance or novel activities. In addition, the sequence-only representation, while computationally efficient, cannot fully capture the three-dimensional structural context that governs enzyme function. Integration of structural information through multimodal architectures could address this limitation and improve prediction accuracy for structure dependent properties like protein stability. Overall, EnzGFM demonstrates that domain-specific pre-training and architectural optimization can substantially improve PLM performance for enzyme-specific applications. The consistent improvements across diverse enzyme-specific downstream tasks, combined with enhanced computational efficiency, position EnzGFM for practical deployment in biotechnology workflows including biocatalyst discovery, metabolic engineering, and therapeutic enzyme development. Moreover, EnzGFM-Agent extends these predictive capabilities into a practical multi-agent framework, enabling users to translate enzyme-specific predictions into executable analysis, candidate design, and experimental decision-making workflows. Together, EnzGFM and EnzGFM-Agent provide an integrated foundation for efficient enzyme understanding, prediction, and engineering. Methods Dataset construction and curation Our training dataset construction followed a hierarchical two-stage approach designed to adapt model representations from general protein patterns to enzyme-specific features. The first stage established foundational protein language understanding through pre-training on 48.28 million protein sequences from UniRef50 33 , totaling 16.67 billion tokens after tokenization. This broad protein corpus enabled EnzGFM to capture evolutionary patterns across diverse protein families. The second stage specialized EnzGFM for enzymatic function through curated enzyme sequences from CAZy 34 , BRENDA 35 , and SABIO-RK 37 . To reduce potential noise from hypothetical or ambiguously annotated entries, we removed records flagged by uncertainty keywords, required each sequence to have an assigned EC number, and removed duplicate sequences. The combined dataset preserved the naturally imbalanced distribution across all seven EC classes, spanning oxidoreductases (EC 1, 16.97%), transferases (EC 2, 34.14%), hydrolases (EC 3, 20.52%), lyases (EC 4, 10.20%), isomerases (EC 5, 6.96%), ligases (EC 6, 11.08%), and translocases (EC 7, 0.13%). Across both pre-training stages, we retained sequences with lengths between 50 and 5,000 amino acids to filter out fragments and unusually long entries that may reflect annotation errors. To prevent data leakage and ensure proper evaluation on evolutionarily distant sequences, we implemented sequence redundancy removal using MMseqs2 56 . The deduplication process employed the following parameters: mmseqs easy-search -db-load-mode 2 -min-seq-id 0.8 -alignment-mode 3 -max-seqs 300 -s 4 -c 0.8. This procedure yielded a final enzyme training dataset of 5.09 million sequences (1.96 billion tokens), while preserving functional diversity across enzyme families and catalytic mechanisms. Validation dataset for the enzyme-specific pre-training stage was constructed through random sampling of 5% of sequences from the curated enzyme dataset, resulting in 264,133 validation sequences. This validation set maintained the same EC class distribution as the training set and underwent identical quality control procedures. Model architecture Tokenization Protein sequences of enzyme underwent tokenization using a 29-token vocabulary adapted from the ESM 13 , encompassing the 20 canonical amino acids plus 9 specialized tokens for sequence boundaries, padding, masking operations, and non-standard residue handling. Input sequences were converted to integer arrays with positional embeddings incorporated to preserve sequential dependencies. A maximum sequence length of 1024 tokens was enforced, with longer sequences truncated and shorter sequences padded using [PAD] tokens. Hybrid Mamba-transformer architecture EnzGFM implements a hybrid architecture that combines Mamba 23 with Transformer blocks 32 to balance computational efficiency with representational capability. The Mamba component processes long-range sequence dependencies through selective SSMs with linear time complexity O ( n ), enabling efficient modeling of variable-length enzyme sequences without the quadratic scaling limitations of traditional attention mechanisms. The Transformer component captures local catalytic motifs and sequence patterns through multi-head attention with configurable attention patterns optimized for protein sequence analysis. EnzGFM architecture employs a controlled layer arrangement strategy governed by offset and period hyperparameters. The offset parameter determines the initial placement of attention layers, while the period parameter controls their frequency throughout the network depth. Layer selection at position l follows: 1 L l ( H ( l ) ) = f attn ( H ( l ) ) if ( l − offset ) mod period = 0 , f SSM ( H ( l ) ) otherwise , where H ( l ) denotes the hidden state at layer l , while f attn and f ssm represent attention-based and SSM-based transformations, respectively. This alternating arrangement facilitates information flow between the efficient sequence modeling capabilities of Mamba and the established protein pattern recognition strengths of Transformer architectures. To enable scale-matched comparisons with ESM2, four model variants with increasing parameter counts were developed, using the offset/period schedule in Eq. ( 1 ) to control the attention–SSM interleaving. EnzGFM-30M uses a denser attention schedule (offset = 2, period = 4), whereas EnzGFM-150M/650M/1.5B use a sparser schedule (offset = 4, period = 8) to preserve sufficient attention capacity at scale under the targeted parameter budgets. All variants further incorporate expert layer arrangements with a fixed schedule (offset = 1, period = 3), ensuring consistent expert integration across model scales. Training settings General protein pre-training Initial model training employed the UniRef50 33 dataset with a masked language modeling (MLM) 57 , where 15% of sequence tokens were randomly masked and reconstructed from surrounding context. The AdamW optimizer 58 was configured with a learning rate of 1 × 10 −4 , weight decay coefficient of 0.01, and a gradient clipping threshold of 1.0. Training proceeded for 4 epochs with gradient accumulation strategies to achieve effective batch sizes of 192 across distributed computing resources. The optimization objective minimized the negative log-likelihood: 2 L MLM = − ∑ i ∈ M log P ( x i ∣ X \ i ; θ ) , where M represents the set of masked token positions, X ⧹ i denotes the input sequence with position i masked, and θ corresponds to learnable model parameters of EnzGFM. Learning rate scheduling implemented linear warmup across the initial 5% of training steps, followed by cosine annealing decay to zero. Early stopping criteria were established based on validation loss plateauing across 5 consecutive evaluation intervals, with assessments conducted every 1000 training steps. Model convergence was monitored through perplexity calculations and embedding quality metrics. All training operations employed bfloat16 mixed precision arithmetic to optimize memory consumption and computational throughput. Enzyme-specific continued pre-training EnzGFM underwent domain adaptation starting from general protein pre-training checkpoints through continued training on the curated enzyme sequence dataset. The learning rate was decreased to 1 × 10 −5 to preserve previously acquired general protein knowledge while enabling enzyme-specific pattern acquisition. We trained for 5 epochs while keeping the remaining optimization settings consistent with the initial pre-training phase. Implementation details Pre-training was executed on 8 NVIDIA A800 GPUs, implementing PyTorch version 2.1.2 with distributed data parallel processing across all available devices. Downstream task evaluations and performance assessments were conducted using 8 NVIDIA RTX 4090 GPUs. Mixed precision training was deployed through automatic mixed precision (AMP) frameworks to reduce memory overhead and accelerate computational workflows across all experimental phases. Downstream task evaluation datasets Enzyme kinetic parameter dataset We constructed both biased and unbiased evaluation sets to assess model generalization across different sequence similarity distributions. For biased evaluation, we employed the DLKcat dataset from UniKP 3 , containing 16,838 samples for k cat prediction with logarithmic transformation. The Michaelis constant ( K m ) prediction dataset comprised 11,722 samples from UniKP, while catalytic efficiency ( k cat / K m ) prediction used 25,831 samples curated by Wang et al. 10 with random partitioning for training and validation phases. All biased datasets underwent tenfold cross-validation with 90%/10% for training/validation splits. For unbiased evaluation addressing evolutionary distance constraints, we adopted the dataset constructed by Wang et al. 10 , encompassing 27,658 unique k cat entries, 42,018 unique K m entries, and 25,831 samples containing both kinetic parameters. We used the training, validation, and test partitions provided with the dataset. For each kinetic parameter, sequences were first clustered using CD-HIT 43 at 40% sequence identity. Clusters, rather than individual samples, were then assigned to ten folds to ensure that no enzyme sequence in the test fold shared more than 40% identity with any sequence in the training folds. The same fold assignments were used for all compared models. Enzyme-reaction correspondence datasets Enzyme-reaction prediction evaluation employed the ReactZyme 47 dataset containing 178,463 positive enzyme-reaction pairs spanning 178,327 unique enzymes and 7726 distinct reactions. The dataset provides three partitioning strategies: (1) temporal splitting (93%/7%), (2) enzyme similarity-based division (95%/5%), and (3) reaction similarity-based separation (91%/9%). We implemented enzyme similarity-based division using Levenshtein distance 48 calculations to enforce a minimum of 60% sequence divergence between training and test sets, yielding 169,724 training pairs and 8,739 test pairs with a 95%/5% distribution. Negative sample construction using enzyme and reaction similarity ranking. Enzyme functional annotation datasets EC number prediction evaluation used the curated SwissProt dataset from Yu et al. 20 , comprising 220,000 enzyme sequences with complete four-digit EC classifications. Sequence filtering ensured EC annotation coverage across all hierarchical levels. Data partitioning employed MMSeqs2 56 clustering at multiple sequence identity thresholds (10%, 30%, 50%, 70%) to prevent data leakage and assess generalization under varying evolutionary distances. The dataset follows established 80%/20% splits with five-fold cross-validation. Model validation employed two independent test sets, namely the NEW-392 containing 392 experimentally characterized sequences and Price-149 comprising 149 validated entries 51 . Single amino acid substitution datasets Mutation effect prediction analysis compiled 225 enzyme mutation studies from previous works 59 – 61 . Quality control required minimum 5 mutations per dataset with corresponding experimental measurements to ensure statistical reliability. The compiled dataset encompasses three functional categories reflecting different aspects of enzyme behavior: (1) protein stability parameters (127 sequences), (2) substrate affinity parameters (11 sequences), and (3) catalytic activity parameters (87 sequences). Distribution analysis revealed 189 enzymes with limited mutation sets (fewer than 50 mutations), 20 enzymes with moderate coverage (50-100 mutations), and 16 enzymes with extensive mutational characterization (exceeding 100 mutations). Downstream evaluation settings Enzyme kinetic parameter prediction Kinetic parameter prediction employed the UniKP framework 3 . Enzyme sequences were embedded through pre-trained PLM (ESM 13 , 14 , ProtTrans 12 , and our EnzGFM), where individual amino acid residues were transformed into high-dimensional representations at the final hidden layer. These embeddings were then processed via mean pooling aggregation and normalization to generate representative enzyme feature vectors. Substrate molecular structures were encoded using a pre-trained SMILES Transformer model 62 , converting chemical structures to SMILES notation before processing into 1024-dimensional feature vectors. The substrate representation combined mean and maximum pooling operations from the final layer with initial outputs from both the ultimate and penultimate layers through concatenation. An extra trees ensemble model was trained on concatenated enzyme-substrate feature representations to predict three kinetic parameters, namely the turnover number k cat (s −1 ), Michaelis constant K m (M), and catalytic efficiency k cat / K m (M −1 s −1 ). For optimal pH prediction tasks, only enzyme sequence-derived features served as model input. Training implementation used scikit-learn with 10-fold cross-validation protocols. Enzyme-reaction prediction The ReactZyme framework 47 was adapted to incorporate pre-trained PLM embeddings as input features. Cross-attention mechanisms operated with 8 attention heads and 128-dimensional hidden projections to align enzyme embeddings with reaction representations. Molecular embeddings (1024 dimensions from molecular encoders) were processed through a three-layer linear network implementing progressive dimensionality reduction (1024 → 256 → 256 → 128). Enzyme sequence embeddings underwent parallel processing through analogous architecture with matching dimensionality reduction to 128 dimensions. Bidirectional cross-attention operated between enzyme and molecular representations, followed by a Transformer encoder-decoder architecture comprising 4 layers each. Final predictions were generated through a multi-layer perceptron with sequential dimensionality reduction (128 → 64 → 16 → 1). ReactZyme training employed the AdamW optimizer with learning rate 5 × 10 −4 and weight decay 1 × 10 −3 . Training schedule implemented linear warmup across the initial 5 epochs, followed by linear decay throughout the remaining training period. Gradient accumulation and clipping maintained optimization stability with maximum gradient norm of 1.0. Binary cross-entropy loss with logits served as the objective function, incorporating class weights to address imbalanced positive and negative enzyme-substrate pair distributions. Training proceeded for up to 50 epochs with early stopping patience of 5 epochs based on validation loss monitoring. Enzyme function prediction EC number prediction utilized the CLEAN framework 20 , implementing contrastive learning to distinguish enzymes with identical EC numbers from functionally distinct variants. Protein sequences were processed through pre-trained PLMs to generate sequence embeddings, which served as input to a supervised contrastive learning network. Positive pairs consisted of enzymes sharing identical EC classifications, while negative pairs were sampled from different EC classes using hierarchical sampling strategies. The training procedure optimized distance minimization between positive sequence pairs and their corresponding anchor sequences while maximizing distances between negative sequences and anchors. EC number prediction compared query sequence embeddings against EC number representations, derived through averaging embeddings of all enzymes within each classification. Pairwise Euclidean distance served as the similarity metric for functional assessment. CLEAN training employed triplet margin loss with margin parameter 1.0. Training epochs scaled according to sequence identity splits following previous work 20 : 10% identity (2000 epochs), 30% (2500), 50% (3500), 70% (5000), and 100% (7000). The Adam optimizer operated with learning rate 5 × 10 −4 and adaptive scheduling applied every 100 epochs. Hard negative mining was executed every 100 epochs using 30 nearest neighbors to update distance mappings for improved contrastive learning performance. Single-point mutation effect prediction Mutation effect prediction employed EnzGFM without additional supervision to assess protein functional impacts. It implemented a wild-type marginal scoring strategy, computing log-probability differences between mutant and wild-type amino acids at specific positions. For each protein sequence, mutations were encoded using ESM tokenization 13 , with the model generating token probabilities through masked language modeling inference. Formally, the scoring calculation followed: 3 S ( i , a w t → a m t ) = log p ( x i = a m t ∣ x \ i ) − log p ( x i = a w t ∣ x \ i ) , where S ( i , a w t → a m t ) represents the mutation score for position i , x i denotes the amino acid at position i , a w t and a m t correspond to wild-type and mutant amino acids respectively, x \ i represents the sequence context with position i masked during inference, and p ( x i = a ∣ x \ i ) indicates the conditional probability of amino acid a at position i given the masked sequence context. Evaluation metrics Model performance evaluation employed task-specific evaluation metrics tailored to each downstream task. For enzyme kinetic parameter prediction, we measured the coefficient of determination as R 2 = 1 − ∑ i = 1 n ( y i − y ^ i ) 2 ∑ i = 1 n ( y i − ȳ ) 2 , where y i represents observed values, y ^ i denotes predicted values, and ȳ corresponds to the mean of observed values. Additional metrics included root mean square error RMSE = 1 n ∑ i = 1 n ( y i − y ^ i ) 2 , Pearson correlation coefficient PCC = ∑ i = 1 n ( y ^ i − y ^ ¯ ) ( y i − ȳ ) ∑ i = 1 n ( y ^ i − y ^ ¯ ) 2 ∑ i = 1 n ( y i − ȳ ) 2 , and Spearman correlation coefficient computed as the Pearson correlation between ranked predicted and observed values. Enzyme-reaction prediction tasks utilized ranking-based metrics including mean average precision (MAP) MAP = 1 ∣ Q ∣ ∑ q = 1 ∣ Q ∣ 1 m q ∑ k = 1 N q P q ( k ) rel q ( k ) , where ∣ Q ∣ denotes the number of queries, m q is the number of relevant items for query q , N q is the number of ranked candidate items for query q , P q ( k ) denotes the precision at rank k , and rel q ( k ) indicates whether the item at rank k is relevant. Normalized Discounted Cumulative Gain at cutoff k (NDCG@ k ) was calculated as NDCG @ k = 1 ∣ Q ∣ ∑ i = 1 ∣ Q ∣ DCG i @ k IDCG i @ k , where for query i at cutoff k , DCG i @ k = ∑ j = 1 k 2 rel i , j − 1 log 2 ( j + 1 ) , with rel i , j denoting the graded relevance of the result at rank j , and IDCG i @ k denoting the maximum possible DCG at cutoff k , obtained by sorting documents for query i in descending order of relevance. Top-K accuracy for K ∈ {1, 3, 5, 10} was defined as Top − K = 1 n ∑ i = 1 n I ( y i ∈ Y ^ i , K ) , where I represents the indicator function and Y ^ i , K denotes top-K predictions for sample i . EC number prediction employed standard classification metrics, namely the precision Precision = TP TP + FP and recall Recall = TP TP + FN , where TP, FP, and FN represent true positives, false positives, and false negatives respectively. F1-score was computed as F1 = 2 × Precision × Recall Precision + Recall , accompanied by AUC analysis, which summarizes the trade-off between true positive rate and false positive rate across classification thresholds. Single-point mutation effect prediction utilized correlation-based metrics including SCC, PCC, and Kendall’s tau correlation coefficient τ = 2 n ( n − 1 ) ∑ i < j sgn ( y ^ i − y ^ j ) sgn ( y i − y j ) , where sgn denotes the sign function, y i denotes the experimental mutation effect, and y ^ i denotes the predicted mutation effect. Compared baselines Benchmarking was conducted against established Transformer-based PLMs that have demonstrated superior performance in enzymatic prediction tasks. Baseline models included ESM2 13 , ESM-C 36 , ESM-1V 18 , and ProtTrans 12 , representing current state-of-the-art approaches in protein sequence modeling with extensive validation across diverse protein function prediction applications. All models underwent evaluation under identical experimental conditions, employing uniform computational resources, standardized data preprocessing pipelines, and consistent evaluation frameworks to eliminate potential confounding factors. EnzGFM-Agent design The EnzGFM-Agent for enzyme engineering employs a hierarchical multi-agent architecture comprising an Orchestrator, perception and task layers, execution and functional prediction layers, and a result analysis layer. The Orchestrator serves as the central routing and scheduling hub. It classifies user inputs, including text queries and uploaded files, to determine the processing pathway: pure Q&A queries are handled by the Chat Agent, whereas task-oriented inputs are processed by the Data Agent and Task Agent, which generate structured task descriptions and initiate the execution pipeline. The execution layer features an iterative code subsystem composed of the Code, Execute, and Analyst Agents, supporting automated error detection and correction to accommodate diverse data formats and computational errors. The Mutation Agent predicts single-point mutation effects and integrates kinetic parameter predictions ( k cat and K m ), using models such as EnzGFM, ESM2, or ProtT5-XL, either across all sites or for specified positions. The Function Agent supports enzyme retrieval and database-driven enzyme discovery through BLASTp and UniProt APIs, providing sequence homology and structured information. The Result Analyst Agent consolidates outputs from all execution agents, performing scientific interpretation, ranking, and visualization. The system supports multi-round iterative analysis through feedback integration. User interaction is provided via the Chainlit interface, with all agents sharing a single LLM backbone (Deepseek-R1 63 ) for specialized behavior through context injection. Computationally intensive model inference is delegated to GPU servers and executed asynchronously, supporting parallel processing of multiple sequences and enabling end-to-end directed evolution workflows within a single session. Experimental validation Plasmid construction Using EnzGFM-Agent, all possible single-point mutations of secreted alkaline phosphatase (SEAP) were scored, and the top 10 candidates from each predictive model (EnzGFM-1.5B and ESM2-3B) were selected. This yielded a total of 15 unique variants, including five variants predicted by both models. The wild-type SEAP coding sequence was derived from the pSEAP2-Control plasmid sequence, GenBank accession U89938.1 . The mutation sites of the selected SEAP variants are shown in Fig. 6 , and both the mutation sites and full sequences are provided in Supplementary Data. The 6 × His tag was fused to the C-terminus of the wild-type and selected mutant SEAP sequences. These constructs were cloned into a laboratory-optimized expression vector, pCMV-SEAP, containing the BSD resistance gene and multiple regulatory elements, as described previously 64 . All mutant plasmids were generated by General Biology (Anhui) Co., Ltd. The resultant plasmids were identified by enzyme digestion and sequencing. Plasmids generated in this study are available from the corresponding author upon reasonable request. Cell culture and transfection CHO-S cells, which are derived from Chinese hamster ovary cells, were obtained from Thermo Fisher Scientific (Cat. No. A11557-01) and maintained in DMEM/F-12 medium supplemented with 10% fetal bovine serum (ProteinEasy Biological Co. Ltd., Xinxiang, China). Cells were seeded at 2 × 10 5 cells per well in 24-well plates and transfected at 70–90% confluence using Lipofectamine 6000 (C0526, Beyotime Biotechnology, Shanghai, China) according to the manufacturer’s instructions. CHO-S cells transfected with the vector lacking the SEAP gene were used as a negative control. SEAP activity assay SEAP enzymatic activity was measured using a commercial alkaline phosphatase kit (P0321M, Beyotime Biotechnology, Shanghai, China) in 96-well plates. Briefly, 48 h post-transfection, cell culture supernatants were collected and centrifuged at 100 × g for 5 min. The supernatants were transferred to fresh tubes and centrifuged at 13,800 × g for 5 min to remove cellular debris. The clarified supernatants were then diluted 1:10, and 50 μL of each sample was added to the assay wells. Samples were then mixed 1:1 with chromogenic substrate and incubated at 37 °C for 10 min. The reaction was terminated by adding 50 μL of stop solution. Absorbance at 405 nm was measured, and enzyme activity was quantified based on a standard curve. All experiments were performed in triplicate. Western blot analysis SEAP protein expression levels were assessed by Western blot. Culture supernatants were collected, mixed with 5 × loading buffer at a 4:1 ratio, and boiled for 10 min. Equal amounts of protein were separated on 10% SDS-PAGE gels and transferred to PVDF membranes. Membranes were blocked with 5% nonfat milk and incubated overnight at 4 °C with an HRP-conjugated anti-His-tag monoclonal antibody (HRP-conjugated 6 × His, His-Tag Mouse McAb, Proteintech, Cat. No. HRP-66005, clone 1B7G5, lot 21016595; 1:10,000). Chemiluminescent signals were captured using a Tanon imaging system (Shanghai, China) and quantified with ImageJ. Antibody specificity was assessed by detection of a specific band at the expected molecular weight of C-terminal 6 × His-tagged SEAP proteins, with CHO-S cells transfected with the vector lacking the SEAP gene serving as a negative control. Statistical analysis Statistical comparisons between model performances were conducted through paired hypothesis testing protocols. Data distribution normality was first assessed using the Shapiro-Wilk test to determine appropriate statistical procedures. For normally distributed performance metrics, paired t-tests were applied to evaluate differences between model variants. Non-normally distributed data underwent analysis using Wilcoxon signed-rank tests to assess statistical differences while accounting for non-parametric distributions. All computational analyses, data processing workflows, and result visualization were implemented using Python version 3.10.14. The analytical framework incorporated NumPy version 1.26.3 for numerical operations, SciPy version 1.13.0 for statistical testing procedures, scikit-learn version 1.4.2 for machine learning implementations and clustering metrics, and Matplotlib version 3.8.4 for figure generation and data visualization. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Supplementary information Supplementary Information Transparent Peer Review file Reporting Summary Description of Additional Supplementary Files Supplementary Data Source data Source Data Acknowledgements We thank all members of our research groups for their helpful discussions and technical support. Author contributions C.W. and Y.S. designed and supervised the project. C.W. designed the methodology and performed the computational experiments. C.W., M.L. and X.Z. performed the data analysis and validation experiments. T.W., S.G. and W.L. designed and performed the wet-lab experiments. Y.G.W. provided technical guidance and theoretical framework. C.W., M.L., Y.Y., T.W. and Y.S. analyzed the results. C.W., Y.Y., T.W. and Y.S. wrote the manuscript. C.W., Y.Y. and T.W. revised the manuscript. Peer review Peer review information : Nature Communications thanks Liangzhen Zheng and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available. Funding This work was supported by grants from the National Natural Science Foundation of China [U23A20270, 82572359]; Henan Provincial Key R&D Program [251111310700]; Central Plains Leading Talent Program [264200510002]; Natural Science Foundation of Henan Province [262300422160]; Joint Fund of the Henan Provincial Science and Technology Research and Development Program [252103810371]; Henan Province Science and Technology Research [252102310363]; Young Elite Scientists Sponsorship Program by HAST [2025HYTP013]. Data availability All data analyzed in this manuscript are available from the public datasets. The UniRef50 dataset used for protein language model pre-training is publicly available from the UniProt Consortium at https://www.uniprot.org/uniref/ . Curated enzyme sequence datasets were obtained from CAZy ( http://www.cazy.org/ ), BRENDA , and SABIO-RK ( https://sabiork.h-its.org/ ) databases, each under their respective access policies. The data generated and used in this study have been deposited in the Zenodo repository at 10.5281/zenodo.19412283. Source data are provided as a Source Data file. Additional supporting figures and tables are provided in the Supplementary Information. Source data are provided with this paper. Code availability The code used in this study has been deposited in the Zenodo repository under 10.5281/zenodo.19412283 and is also publicly available on GitHub at https://github.com/DeepBxM/EnzGFM . Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Yi Yu, Email: yuyi@xxmu.edu.cn. Tianyun Wang, Email: wtianyuncn@126.com. Yiqing Shen, Email: yshen92@jhu.edu. Supplementary information The online version contains supplementary material available at 10.1038/s41467-026-75283-3. References Associated Data Supplementary Materials Supplementary Information Transparent Peer Review file Reporting Summary Description of Additional Supplementary Files Supplementary Data Source Data Data Availability Statement All data analyzed in this manuscript are available from the public datasets. The UniRef50 dataset used for protein language model pre-training is publicly available from the UniProt Consortium at https://www.uniprot.org/uniref/ . Curated enzyme sequence datasets were obtained from CAZy ( http://www.cazy.org/ ), BRENDA , and SABIO-RK ( https://sabiork.h-its.org/ ) databases, each under their respective access policies. The data generated and used in this study have been deposited in the Zenodo repository at 10.5281/zenodo.19412283. Source data are provided as a Source Data file. Additional supporting figures and tables are provided in the Supplementary Information. Source data are provided with this paper. The code used in this study has been deposited in the Zenodo repository under 10.5281/zenodo.19412283 and is also publicly available on GitHub at https://github.com/DeepBxM/EnzGFM .