Accelerated discovery of highly active enzyme nanohybrids with parallelized Bayesian optimization in hybrid space

Domain: chemistry
Published: 2026-03-07
Editorial outcome: major_revision

Abstract
Artificial intelligence (AI) has significantly advanced protein engineering, enabling rapid enzyme evolution for diverse applications. However, the fragile nature of biomacromolecules requires enzyme immobilization to preserve catalytic activity under harsh industrial conditions, which often restricts substrate diffusion and reduces enzymatic activity. This challenge demands extensive trial-and-error experiments to optimize immobilized carriers with high activity for different enzymes. Here we show a machine-learning-guided workflow along with an algorithm named parallelized hybrid-space Bayesian optimization (PHBO) to accelerate the discovery of nanocarriers for specific enzymes and reactions. Leveraging prior knowledge, machine learning and iterative feedback, within limited number of experiments, this workflow explores the reaction space of over 107 experiments and achieves activity recovery of 100%, 90%, and 79% for glucose oxidase, catalase, and Candida Antarctica lipase B, respectively. These results demonstrate that data-efficient optimization can substantially accelerate the discovery of enzyme nanohybrids with high catalytic activity across diverse enzymatic systems.

Introduction
From photosynthesis and nitrogen fixation to phosphorylation and carbohydrate metabolization, enzymes catalyze the transformation of energy and substance with strict precision and high efficiency1,2. Inspired by nature’s repertoire, scientists have employed various computational and chemical methods, including directed evolution3 and de novo design4, to engineer enzymes with enhanced functionality or activities not observed in naturally occurring enzymes. Breakthroughs in artificial intelligence (AI) and rapid accumulation of protein structure database also facilitate more sophisticated understanding of the protein sequences-functions relationship5, from which researchers have established more reliable libraries, more efficient high-throughput screening techniques, and more complicated generation mechanisms to explore a broader space for new enzymes.

When exposed to the real industrial conditions, these biocatalysts still suffer from the fragile nature of biomacromolecules that undergo the dynamic conformational changes with high frequency and delicacy6. Immobilizing enzymes through bounding to a solid support represents an effective and widely used chemical approach for enzyme engineering, which can reduce the extent of free diffusion of enzymes in solution and provide a confined environment for the three-dimensional structures. This approach can improve enzyme stability in industrial applications but often leads to a serious loss of activity due to mass transfer resistance and/or synthetic conditions7. The diversity of enzymes and catalyzed reactions requires customized carriers to ensure activity when improving the enzyme stability and reusability8. However, current immobilization approaches rely on large amounts of trial-and-error experiments based on the scientific intuition from biochemists. Therefore, it is highly desirable to develop an efficient decision-making method that can simultaneously recommend multiple synthetic parameters and accelerate the discovery of carriers which suit specific enzymes, achieving the holistic evolution of immobilized biocatalysts.

AI-driven decision-making methods have emerged as effective tools to assist experimental design, guide data collection, and interpret scientific insights in the era of intelligent synthesis9,10. As for the immobilized enzymes, due to the high value of enzyme molecules, sophistication of enzymatic assays, large amount of carrier materials, and variation of immobilizing experimental protocols, it is difficult to generate a large database of structure-activity relationship with high consistency and accuracy to train a delicate AI model. Compared with expert systems and evolution algorithms, Bayesian optimization (BO), as an efficient global optimization algorithm for expensive black-box objective functions11,12, effectively mitigates the small sample problem by significantly enhancing sampling efficiency, and thus accelerates the discovery of new materials including quantum dots13, polymers14, copolymers15, photocatalyst16, and photosensitizer17. For the intelligent synthetic condition recommendation of immobilized enzymes, two challenges remain to be addressed: hybrid variable spaces and parallelized sampling. The variety of synthetic precursors is coupled with their concentration, which generates a hybrid variable space. BO is originally designed for pure continuous variable spaces. While existing studies have extended the application of BO to pure categorical variable spaces18,19, BO for categorical-continuous hybrid variable spaces is rarely explored. Gryffin20 performs well in hybrid spaces but relies heavily on prior domain knowledge, as its effectiveness strongly depends on the quality and relevance of user-provided descriptors. If these descriptors are poorly correlated with the optimization objective, the advantages of Gryffin may diminish significantly. Conducting multiple sets of experiments simultaneously can significantly reduce the cost of a single trail of synthesis and evaluation, which requires multiple recommendations from the algorithm at the same time. The basic multiple acquisition function method for parallelized sampling tends to cause sample points to cluster, resulting in highly homogeneous evaluations and inefficient utilization of high-throughput experimental devices21.

In this study, an AI-assisted workflow has been established to accelerate the evolution of immobilized enzymes with high activity targeting specific enzymes and reactions. Metal-organic frameworks (MOFs) synthesized by the coordination of metal nodes and organic ligands22,23 were chosen as the starting point for designing the nanocarriers for enzyme immobilization. Our group24,25 and fellow colleagues26 pioneered the development of a one-pot co-precipitation method to immobilize enzymes by MOFs in aqueous solution inspired by biomineralization in nature, where enzyme molecules, metal ions and organic ligands were easily mixed in aqueous solution. This method is convenient and applicable to a wide range of enzymes, and the nanocomposites can maintain good stability and enzymatic activity. To discover the most suitable nanocarriers for enzyme encapsulation, 7 soluble zinc salts and 17 soluble organic ligands were screened to generate a large scope of reaction space, which on the other hand geometrically increased the number of reactions required to explore all the possibilities, far from what human researchers could endure (Supplementary Tables 1, 2, Supplementary Note 1). Therefore, there were four crucial issues to be solved: (i) small sample problem, (ii) categorical-continuous hybrid variable space, (iii) parallelized recommendations of experimental conditions, and (iv) specificity of different enzymes. Aiming at addressing these challenges, we integrated probabilistic reparameterization, customized kernel, and the Nearby Liar method to propose the parallelized hybrid-space Bayesian optimization (PHBO) algorithm, accelerating the discovery of high-performance immobilized enzyme nanohybrids for specific enzymes.

Results
Chemical space and workflow design
Two categorical variables (zinc salt type and ligand type) and three continuous variables (zinc salt concentration, ligand concentration, and reaction time) were chosen (Fig. 1a). 7 water-soluble zinc salts and 17 water-soluble imidazolyl or similar ligands were screened and listed in Fig. 2, whose solubilities at room temperature and pressure were tested and used as the upper limits for precursor concentrations (Supplementary Tables 1, 2). During the synthetic process, the volumes of precursor solutions were fixed, and the concentrations were varied in the limits of respective solubility. The volumes and concentrations of enzyme solutions were fixed for each type. In this context, over 107 experiments were required to build an extensive database for these precursors, which was apparently impractical due to the limitation of resources and time (Supplementary Note 1). For the proposed workflow, PHBO can suggest 8 sets of experimental conditions in parallel based on the prior knowledge database, which were then artificially implemented and fed back to the model. PHBO would recommend again in these iterations, until the requirements of researchers were reached (Fig. 1b). Activity recovery (AR) was selected as the objective to optimize the overall immobilized enzymes, which contained comprehensive information during the immobilization process, including whether the reaction can form precipitates, encapsulation efficiency (EE), relative activity (RA), and separation and reuse properties (Fig. 1c). Once ideal carriers were found, detailed enzymatic properties including encapsulation efficiency and relative activity were investigated to get better understandings of the newly discovered nanohybrids.

Parallelized hybrid-space Bayesian optimization construction
To address the challenges of modeling and optimization in hybrid categorical-continuous reaction spaces, we developed the PHBO algorithm, which integrated three key innovations: (i) a surrogate model built in a reparameterized continuous distribution space, which allowed customized GP modeling; (ii) an acquisition function defined as an expectation over categorical variable distributions, enabling efficient gradient-based optimization via Monte Carlo gradient estimation; and (iii) the Nearby Liar parallelized sampling method, which incrementally updated the surrogate model to diversify batch recommendations and reduce redundant sampling. A comparison of PHBO with representative hybrid-space or categorical-space BO algorithms was summarized in Table 1.

As mentioned above, the objective of the synthetic optimization was to maximize AR, that was, to identify the sample point with the highest AR within the reaction space using the fewest possible experiments:
(x*,z*) = argmax f(x,z) (1)
(x,z) ∈ X × Z
The BO framework utilized Gaussian processes (GP) as the surrogate model of the reaction space27,28, and PHBO followed this approach. Given the absence of an explicit ordinal relationship between the levels of categorical variables, it was necessary to customize a dedicated GP kernel for them. We first defined the kernel for categorical variables as:
kCat(z,z') = { 1, if z = z'; e^-1/l, if z ≠ z' } (2)
For continuous variables, PHBO used the commonly adopted Matern 5/2 kernel. The complete hybrid kernel was derived from the sum and product of these two kernels, weighted accordingly:
kHyb((x,z),(x',z')) = (1-λ)(kCon(x,x') + kCat(z,z')) + λ(kCon(x,x')kCat(z,z')) (3)
The kernel weights were optimized using maximum likelihood estimation, allowing data-driven adjustment of the contributions from the categorical and continuous components.

The hybrid kernel better captured the characteristics of the variables and avoided the pathological issues of the kernel for categorical variables due to the insufficient data in the early stages of iterations.

Based on the surrogate model, the BO framework iteratively constructed and optimized the acquisition function to suggest a sample point for each iteration round. Taking round t as an example:
(xt,zt) = argmax α(x,z) (4)
(x,z) ∈ X × Z
To utilize efficient gradient-based optimization algorithms in the discontinuous optimization problem of acquisition function, PHBO utilized the idea of probabilistic reparameterization. Specifically, PHBO replaced the categorical variables in the acquisition function with probability distributions (Fig. 3a, b), and accordingly replaced the acquisition function with its expectation (Fig. 3c, d):
(xt,θt) = argmax Ez[α(x,z)] (5)
(x,θ) ∈ X × Θ
where z was drawn from distribution p(Z|θ) throughout Eqs. (5)-(8).

The gradient of the acquisition function with respect to the continuous variables could be directly derived:
∂/∂x Ez[α(x,z)] = ∑ p(z|θ) ∂/∂x α(x,z) (6)
z∈Z
The gradient with respect to the categorical variables could be converted into an expectation form by variable substitution and then estimated through Monte Carlo sampling29. Specifically, we referred to the REINFORCE strategy gradient estimation method from reinforcement learning30:
∂/∂θ Ez[α(x,z)] = ∑ α(x,z) ∂/∂θ p(z|θ) ≈ 1/NMC ∑ α(x,z̃i) (I - p(z̃i|θ)) (7)
z∈Z i=1
To reduce the variance of the REINFORCE estimator in (Eq. (7)), we adopted a control variate technique where the acquisition function α was adjusted by subtracting a proportional baseline β, defined as B times the running mean of α from previous iterations. This yielded a modified gradient estimator (Eq. (8)) with significantly lower variance. In practice, we used 128 Monte Carlo samples for each REINFORCE estimation, which represented a practical trade-off between variance reduction and computational efficiency.
∂/∂θ Ez[α(x,z)] ≈ 1/NMC ∑ (α(x,z̃i) - β) (I - p(z̃i|θ)) (8)
i=1
With the gradients with respect to these two parts, PHBO employed the L-BFGS algorithm for acquisition function optimization to suggest the next sample point21.

To parallelize sampling, inspired by the Constant Liar method31, we developed the Nearby Liar method. In Liar methods, parallelized recommendations were sequentially generated by adding each newly suggested sample point to the training set with a temporary “Liar value”, that was, a fabricated evaluation used to update the surrogate model before actual evaluation. In Nearby Liar, instead of assigning a fixed value, the Liar value was set to the average objective value of the closest previously evaluated sample points in the reaction space, improving estimation accuracy and reducing redundant sampling. Specifically, PHBO selected a series of sample points that are closest to the suggested sample point from the evaluated dataset:
select {(xi,zi) | (xi,zi) ∈ D sorted by Dis((xi,zi), (xt,zt))} (9)
According to the objective function values of these points, the Liar objective function value of the suggested sample point was calculated and merged into the dataset (Fig. 3e):
Dt,1 = Dt ∪ {(xt,zt), 1/NNL ∑ y i=1} (10)
Then the surrogate model was updated with the new dataset. The suggestion was repeated until the number of sample points to evaluate reaches the capacity of experimental devices. Then, batch experiments were conducted to replace all Liar values with exact results, thus completing round t of iterations. PHBO continued iterations until the results converged, or the budget was exhausted.

The computational complexity of PHBO was primarily composed of two parts: GP model training and acquisition function optimization. The GP model training incurred a complexity of O(n3), where n was the number of observed samples. The acquisition function, estimated via Monte Carlo sampling over reparametrized distributions, had an optimization complexity of O(TMn2), where T was the number of optimization iterations, and M was the sample size. The parallelized sampling using Nearby Liar introduced negligible additional overhead. Given that experimental evaluations were typically more time-consuming than these computations, PHBO offered a favorable trade-off between decision-making efficiency and computational cost.

PHBO for the accelerated discovery of GOx nanohybrids
We first conducted ablation study of PHBO and chose glucose oxidase (GOx) as the model enzyme. One variable at a time (OVAT) by experienced chemists was employed to generate the prior knowledge database and used as a control optimization as well. By the alteration of zinc salt type, ligand type, reaction time, molar ratio of precursors, and dilution ratio of precursors, the AR was increased from 5.2% of ZIF-8 to 86.1% of Zn(2,4-dmIM)2 through 47 experiments in 6 iterations (Fig. 4a, Supplementary Figs. 1, 2). Despite the nearly 17-fold increase in AR, it was difficult for OVAT to further suggest new sample points purely based on chemist's prior experience.

LocalSearch, a BO variant that combined confidence bounds and local search32, was a commonly used method for handling hybrid variable spaces. During the optimization of the acquisition function, LocalSearch alternated between fixing the categorical and continuous parts of the operating variables, optimizing the other part using gradient-based optimization algorithm for continuous variables and local search algorithm for categorical variables, thus suggesting the final sample point. As a comparison, the only difference between LocalSearch and PHBO was limited to the optimization of the acquisition function; and all other configurations were identical, including the Nearby Liar method. The AR optimized by LocalSearch exceeded OVAT in rounds 7 and 8 and reached nearly 100%, proving the effectiveness of BO framework and Nearby Liar technique (Fig. 4b). PHBO, however, outperformed OVAT and LocalSearch throughout the first 6 rounds, and overtook LocalSearch in round 10 with an AR of around 110%, demonstrating the validity of probabilistic reparameterization (Fig. 4c).

The proportions of immobilized GOx with AR over 80% generated by OVAT, LocalSearch, and PHBO were 6.4%, 12.1%, and 18.0%, respectively (Fig. 4d). The reaction conditions screened by the different methods were further visualized by the precursor type, concentration, and reaction time (Fig. 4e-h, Supplementary Fig. 3). Limited by experimental experience, the chemical space explored by OVAT was regular and concentrated, rarely involving high precursor concentrations or long reaction time. LocalSearch and PHBO explored the entire variable space more evenly, with attention to long reaction time and abundant combinations of different precursors. PHBO especially focused on the low concentrations, which may produce highly active immobilized enzyme at a lower cost. Overall, PHBO was able to efficiently recommend a rich set of experimental conditions with high enzymatic activity in parallel.

The obtained reaction conditions with high activity were repeated and summarized in Supplementary Table 3. The combination of Zn(Ac)2 and 2-tIM in the vicinity of a ligand/zinc molar ratio of 1.0 yield a concentrated region of high AR, which were labeled as group B, and the other combinations were labeled as group A. Scanning electron microscopy (SEM) images showed diverse morphology ranging from micron-sized lumpy precipitates to tens of nanometer-sized spherical-like particles (Fig. 4i, Supplementary Fig. 4). We discovered that enzyme nanohybrids synthesized at low ligand/zinc molar ratios presented amorphous X-ray diffraction (XRD) curves (Fig. 4j). Attenuated total reflection flourier transformed infrared spectroscopy (ATR-FTIR) confirmed the immobilization of enzymes by the appearance of the amide band (1700-1500 cm-1) (Fig. 4k). The nitrogen adsorption and desorption curves showed that GOx@ZIF-8 presented a typical microporous adsorption profile with a peak at 1.43 nm and a BET surface area of 1581 m2 g-1 (Fig. 4l, m). Whereas, GOx@A1 and GOx@A6 presented weak mesoporous hysteresis loops that may be consistent with flaky granular particles. The pore size distribution indicated the presence of mesopores of GOx@A1 and GOx@A6, which can be attributed to particle stacking according to SEM images. BET surface areas of these two groups were 181.1 and 357.4 m2 g-1, respectively. The nitrogen adsorption and desorption curves and the pore size distribution of GOx@B6 indicated that there were almost no porous structures, with a BET surface area of only 50.3 m2 g-1.

For GOx, the optimal carrier with highest AR was B6, which was synthesized by Zn(Ac)2 and 2-tIM with a ligand/zinc molar ratio of 1.4. The encapsulation efficiency and relative activity of GOx@B6 achieved nearly 100% as well. We further immobilized different types of enzymes using ZIF-8 and B6, including lactate oxidase (LOx), horseradish peroxidase (HRP), catalase (CAT), Candida Antarctica lipase B (CALB), and keratinase (Ker). Supplementary Fig. 5 showed that the encapsulation efficiency and relative activity of the same carrier for different enzymes varied very much, but overall B6 outperformed ZIF-8 in five groups. B6 can achieve an encapsulation efficiency of more than 90% for five enzymes while ZIF-8 can barely immobilize three kinds of enzymes. As for relative activity, B6 outperformed ZIF-8 except for the case of CAT. This suggested that on the one hand, B6 discovered by PHBO had significant advantages of enzyme immobilization performance compared to ZIF-8. On the other hand, the most suitable carriers for different enzymes may have different compositions and structures, which required dedicated optimization with the overall enzyme nanohybrids as the target.

Model expansion for CAT and CALB nanohybrids
We then incorporated transfer learning into PHBO to utilize the transcendental relationship of synthetic conditions and properties from GOx to accelerate the discovery of other immobilized enzymes. The core idea of transfer learning was to extract knowledge from the source data and migrate it to the target domain to improve learning-based algorithms’ performance33-35. Specifically, we employed two techniques: GP prior parameter update and BO warm up. Based on the optimization process of GOx, GP prior parameter provided a reference for the mean parameter of the GP surrogate model, while BO warm up expanded the initial dataset of PHBO. The first round of reaction conditions for CAT were recommended by PHBO based on the results of GOx, where the maximum of AR reached over 90%. After 6 rounds of iterations, structures with AR of 100% were discovered (Fig. 5a). Compared to 11 rounds of iterations from 11.4% to 106.6% of GOx immobilization, the validity of transfer learning was confirmed. The visualization of experimental conditions showed that PHBO explored the region of low precursor concentrations and medium reaction times to achieve high AR (Fig. 5b, c). The preferred synthetic conditions were repeated and summarized in Supplementary Table 4, some of which generated large lumpy precipitates and others small amorphous particles of tens of nanometers (Supplementary Fig. 6). Similarly, CALB was set as the target. The first round achieved an AR over 50%, compared to 1.2% of ZIF-8 as the starting point. After 10 rounds of iterations, structures with AR more than 80% were discovered (Fig. 5d). Low-to-medium precursor concentrations and long reaction times tended to produce structures with high activity as shown in Fig. 5e, f. Supplementary Table 5 showed the results of group D (ZnBr2 and 2,4-dmIM) and group E (Zn(NO3)2 and 2,4-dmIM) with a generally broad ligand/zinc molar ratio. The morphology varied from precursor types. Group D were overall irregular particles ranging from 100 nm to micrometers with rough surfaces, while the amorphous particles in group E were much smaller, with sizes of only 20-50 nm (Supplementary Fig. 7).

Furthermore, to quantitatively assess the effect of transfer learning, we compared PHBO to a baseline version initialized without GP prior updates nor warm-up data. For both CAT and CALB, one batch (8 parallel conditions) was recommended by the control baseline and experimentally evaluated. The results presented in Supplementary Tables 6 and 7 showed that in the first iteration, the transfer-enabled PHBO achieved average AR of 71.31% for CAT and 28.49% for CALB, while the control baseline only achieved average AR of 28.74% for CAT and 4.39% for CALB (Supplementary Notes 2, 3). Such improvement validated the effectiveness of transfer learning and the benefit of incorporation prior knowledge from GOx optimization. These results suggested that with the assistance of transfer learning, PHBO was able to effectively utilize the prior information of existing enzymes in the database to accelerate the carriers design for other types of enzymes.

Response surface modeling for given precursors
The optimal carriers for GOx, CAT, and CALB discovered by the PHBO-assisted workflow were different: B6: Zn(2-tIM)2 (2-tIM : Zn(Ac)2 =1.4) for GOx, C2: Zn(HTM)2 (HTM : Zn(Ac)2 =2.0) for CAT, and E6: Zn(2,4-dmIM)2 (2,4-dmIM : Zn(NO3)2 =11.1) for CALB. Combinations near the optimal molar ratio also had relatively high AR. Therefore, we conducted the response surface modeling (RSM) using PHBO surrogate model based on the existing database to provide a thorough understanding of the selected reaction space. Fig. 6a-c and Supplementary Fig. 8 presented a reference for required precursors, and highly active enzyme nanohybrids can be formed over a range of concentration. We then selected several experimental conditions to verify the accuracy of RSM. As shown in Supplementary Table 8, the PHBO predictions agreed with the general trend of the experimental results. In this way, the RSM workflow can be utilized to directly study precursor combinations of interest and provide guidance to chemical engineers in their experiments.

We then set the optimal molar ratios of three enzymes as the starting points and expanded the exploration space to study the properties of enzyme nanohybrids (Fig. 6d-f). As shown in Supplementary Figs. 9-11, the size of the particles gradually increased when increasing ligand concentrations. When the ligand/zinc molar ratio was too low, it was usually difficult to form sufficient precipitates, resulting in low encapsulation efficiency. When the molar ratio was too high, the strong alkaline environment by the hydrolysis of ligand resulted in a more severe loss of enzyme activity, leading to low relative activity and AR. The increase of ligand concentration also decreased the formation of more particles, and CAT and CALB were unable to precipitate when the ligand/zinc molar ratio exceeded 40.

Investigation of high activity-structure relationships
We further performed the detailed structural characterizations of the three optimal enzyme nanohybrids discovered by PHBO to uncover the origin of high activity. Take CALB as the first example. As the starting point, CALB@ZIF-8 achieved an encapsulation efficiency of 100% but suffered from very low relative activity of only 7.6%. Meanwhile, the relative activity of CALB@Zn(2,4-dmIM)2 was 80.4%, resulting in a high AR of 85.1% (Fig. 7a). A scaled reaction was performed to determine the catalytic properties of both free CALB and CALB@Zn(2,4-dmIM)2, where the relative activity and activity recovery showed similar results as above (Supplementary Table 9, Supplementary Note 4).

Energy-dispersive X-ray spectroscopy (EDS) mapping confirmed the dispersed distribution of Zn and N elements from metal salts and ligands, and the distribution of S elements confirmed the immobilization of enzymes (Fig. 7b). The stretching vibration peaks of NC-N at 1420-1400 cm-1 and bending vibrational peaks at 900-650 cm-1 of secondary amine appeared in the curve of CALB@Zn(2,4-dmIM)2, which verified the presence of 2,4-dmIM in the particles. The amide band of 1700-1500 cm-1 in FTIR curves confirmed the immobilization of enzymes (Fig. 7c). XRD patterns verified the crystallinity of CALB@ZIF-8, while CALB@Zn(2,4-dmIM)2 exhibited typical amorphous curves with broad peaks at 10°-20° (Fig. 7d). CALB@ZIF-8 possessed micropores of 0.93 nm and exhibited typical microporous adsorption curves, whereas CALB@Zn(2,4-dmIM)2 had almost no porous structures and very weak adsorption capabilities (Fig. 7e). As discussed earlier in Figure 6d-f and Supplementary Figs. 9-11, by adjusting the ligand-to-metal molar ratio under the same precursor combination, we obtained enzyme nanohybrids with different particle sizes. When the ligand-to-metal molar ratio was too low (e.g., 0.1-4.0 for CALB), it was difficult to form well-defined particles. Once stable particles were obtained, AR decreased as particle size increased. Under optimal conditions, the resulting particles exhibited negligible porosity, suggesting that enzymes were loosely confined within the solid matrix rather than encapsulated within a porous framework. In this context, smaller particles offered a higher surface area-to-volume ratio, enhancing enzyme exposure to substrates and thus contributing to the observed improvement in catalytic activity.

In addition, high-resolution Zn 2p X-ray photoelectron spectroscopy (XPS) spectra revealed a shift in the Zn 2p binding energy of CALB@Zn(2,4-dmIM)2 toward lower energy compared to CALB@ZIF-8 (Figure 7f), indicating a lower partial positive charge for zinc attributed to the additional methyl group of 2,4-dmIM. This would result in the relatively loose and unstable interactions between Zn and N, providing a more flexible environment to maintain enzyme structures and activity. As shown in CD diagrams, CALB released from CALB@Zn(2,4-dmIM)2 showed the well-preserved secondary structures, while CALB released from CALB@ZIF-8 exhibited clear blue-shift of α-helix peaks (Fig. 7g). Moreover, amorphous particles usually suffered from instability, easy breakage, and low reusability. To observe the tolerance of particles to mechanical shear force, the immobilized enzyme suspension was treated at high centrifugation of 16,060× g, 3 min, then resuspended, sampled, and tested for activity. After 10 cycles, the activity of CALB@Zn(2,4-dmIM)2 remained nearly 70%, demonstrating an efficient recovery and excellent stability of the enzyme nanohybrids (Fig. 7h). GOx@Zn(2-tIM)2 and CAT@Zn(HTM)2 were analyzed as well which showed similar improvement compared to GOx@ZIF-8 and CAT@ZIF-8 (Supplementary Figs. 12, 13).

In summary, the proposed integrated workflow effectively enabled the rapid design of high-efficiency immobilization carriers for various enzymes. By loosely confining enzymes within a small solid region, the structural integrity and catalytic activity were well preserved. Leveraging transfer learning techniques, this workflow can be applied to a wide range of other enzymes for different industrial biocatalysis. Moreover, by modifying the reaction space and optimization objectives, the PHBO algorithm can facilitate the efficient design of other structures, thereby accelerating the development and iteration of functional materials.

Methods
Synthesis of hybrid biocatalysts
Hybrid biocatalysts were synthesized as below. Typically, 100 μL protein solution and 160 μL zinc salt solution were simultaneously added into 1.6 mL ligand solution, and then stirred (600 rpm) at room temperature for 5 min. The precipitate was collected by centrifugation of 16,060× g for 5 min and washed 3 times with deionized water. For the synthesis of ZIF-8, 0.31 M Zn(NO3)2 and 1.25 M 2-methylimidazole were used. Concentrations of proteins were 10 mg mL-1 for GOx, 30 mg mL-1 for CALB, and 7.5 mg mL-1 for CAT, respectively. Composites were freeze-dried and stored at 4 ℃ until further characterization.

Characterizations of hybrid biocatalysts
Composites were sputtered with Au in a vacuum for 60 s before being observed with JOEL scanning electron microscopy (SEM) 7900f. Samples were dispersed in ethanol and then dropped onto a carbon-coated copper grid. Transmission electron microscopy (TEM) images were taken at 200 kV by a JEOL 2100 plus equipped with an Oxford energy-dispersive X-ray (EDX) analyzer. Powder X-ray diffraction (XRD) patterns were recorded using a Bruker D8 Advance X-ray diffractometer with a 5° min-1 scanning speed and 5°-50° diffraction angle. Circular dichroism (CD) was performed by Chirascan from Applied Photophysics Ltd. Attenuated total reflection flourier transformed infrared spectroscopy (ATR-FTIR) was performed on Horiba Bruker Vertex 70 with ATR accessories. All characterization data were directly exported from instrument software without additional processing.

Enzymatic activity assays
GOx activity was measured using glucose as substrate and ABTS (2,2′-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) as the chromogenic agent. Typically, horseradish peroxidase (HRP) and ABTS were dissolved in water to make a mixed solution containing 1.5 mg mL-1 HRP and 1.5 mg mL-1ABTS. 30 μL of such mixed solution was added into 940 μL potassium phosphate buffer (K-PB, 10 mM, pH=7.5) containing 50 mM glucose, and 30 μL GOx solution was added at last. The absorbance at 415 nm was recorded for 15 s.

CALB activity was measured using 4-nitrophenyl butyrate (p-NPB) as the substrate. Typically, 4 μL p-NPB was dissolved in 2 mL acetone, which was then added to 40 mL 50 mM K-PB (pH=7.5). 50 μL CALB solution was added into 950 μL of such mixed solution, and the absorbance at 405 nm was recorded for 15 s.

CAT activity was measured using hydrogen peroxide as the substrate. Typically, 20 μL CAT solution was added into 980 μL potassium phosphate buffer (50 mM, pH=7.5) containing 10 mM hydrogen peroxide. The absorbance at 240 nm was recorded for 15 s.

Tyrosinase (Tyr) activity was measured using a commercial assay kit (ACMEC, AC10628). Tyr catalyzed the oxidation of L-Dopa to dopa quinone, which subsequently reacted with MBTH to produce a chromogenic product. Reactions were conducted at 37 ℃ for 20 min and the absorbance was measured at 505 nm. Enzyme activity was calculated based on the absorbance change following the manufacturer’s protocol.

β-Galactosidase (GAL) activity was measured using a commercial assay kit (Solarbio, BC2580). GAL hydrolyzed p-nitrophenyl-β-D-galactopyranoside (pNPGal, a standard galactosidic substrate) to yield p-nitrophenol. Reactions were conducted at 37 ℃ for 30 min and the absorbance was measured at 400 nm. Enzyme activity was calculated based on the absorbance change following the manufacturer’s protocol.

CA activity was measured using the Wilbur-Anderson method based on pH change. Briefly, a 20 mM TAPS buffer (solution A) and CO2-saturated water (solution B) were prepared in advance. For the blank measurement, 2.45 mL of solution A and 50 μL of DI water were added to a 10 mL tube under magnetic stirring, and the time required for the pH to drop from 8.3 to 7.0 after adding solution B was recorded (40 ± 2 s). For enzyme measurements, the same procedure was followed except that 50 μL of free or immobilized CA solution (diluted to 15 μg mL-1) replaced the water. The decrease in time relative to the blank reflects catalytic activity.

For enzymatic activity assay, each sample were tested separately for three times, and the final activity were represented as mean ± SD (n=3) unless otherwise stated.

Calculation of enzymatic properties
Initially, AR was selected as the most direct optimization target because this property contained comprehensive information of the entire immobilization process. The calculation formula was as follows:
AR = apparent activity / total activity (11)
wherein, the total activity was the specific activity of the free enzyme times the total amount of enzyme input in the immobilization reaction.

When the iteration reached a certain stage, the immobilized material with higher AR was found, so its enzymatic properties were comprehensively measured. EE was the percentage of the amount of immobilized enzyme to the total amount of enzyme input. The amount of immobilized enzyme was obtained by subtracting the amount of enzyme in the supernatant after the reaction from the total amount of enzyme input.
EE = Quantity immobilized enzymes / Quantity input enzymes (12)

RA was the proportion of residual enzyme activity in the immobilized material compared to the free enzyme at the same protein concentration. According to the encapsulation efficiency, the immobilized enzyme solution was prepared with the same concentration as the free enzyme, and the activity of both samples were tested under the same conditions. The calculation formula was as follows:
RA = apparent activity / free enzyme activity at the same concentration (13)

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript \"Accelerated discovery of highly active enzyme nanohybirds with parallelized hybrid-space Bayesian optimization\" by Jun Ge et al, is presenting an interesting combination of enzyme immobiliziation and the use of AI to efficiently sample the vast space of all possible experimental parameters to test out in order to arrive at more active catalyst. In general, the idea is excellent as there are not many studies which could address the rather \"trial-error\" area of enzyme immobilization where usually the through put of experimental set up is limited. However, the manuscript is completely lacking experimental evidence which would thoroughly prove the successful outcome of the work. I do not consider spectrophotometric assays with a surrogate/chromogenic substrate as quantitative proof for assessing activity. I would expect scale up on at least few mgs of substrate, isolating and characterizing the product. Also I did not find any experimental procedure for e.g. encapsulation efficiency, calculating activity recovery and reuse performance (how was the catalyst handled between each reaction set up). For those reasons, I can recommend the manuscript for publishing.", "reviewer_2": "The authors present a novel PHBO method to optimize enzyme immobilization carriers based on MOFs. Their approach combines a customized Gaussian-process kernel designed for mixed categorical-continuous variable spaces, probabilistic reparameterization to enable gradient-based optimization, and a batch-suggestion strategy. The authors demonstrate significant improvement in activity recovery (AR) for immobilized enzymes compared to conventional methods. Transfer learning from GOx further reduced the optimization iterations for other enzymes. I found this work exciting, particularly because the proposed PHBO method creatively integrates multiple recent advancements in BO, such as hybrid kernel design, probabilistic reparameterization, and batch sampling. I would not hesitate to offer my support for its publication in Nat Comms and would like to suggest the following additions/corrections/etc. 1. The authors may consider more clearly differentiating their method from established methods such as COMBO (Oh et al., NeurIPS 2019) and Gryffin (by Aspuru-Guzik and co-workers). Perhaps, the authors could include a comparative table summarizing algorithmic components (e.g., kernels, recommendation strategies, etc.) alongside expected computational or performance benefits. 2. I have several queries about the following essential details: 2.1 How were kernel weights (Eq. 3) determined? Were these optimized, e.g., through maximum likelihood estimation or set manually? 2.2 In the Monte Carlo sampling (Eq. 7), what was the number of samples used for the REINFORCE gradient estimator? Were variance-reduction techniques, such as baseline subtraction, employed? 2.3 Please consider discussing or briefly reporting a computational complexity analysis. 3. The authors illustrate valuable transfer-learning capabilities from GOx to CAT and CALB, but the explanation remains qualitative. Please consider providing some quantitative metrics to evaluate transfer effectiveness. 4. The manuscript seems to only associate particle size and amorphousness with high AR. I failed to catch any discussion of a direct causal relationship. If these two mentioned factors are the only/most significant ones, then wouldn’t it make sense to conduct control experiments: synthesizing enzyme nanohybrids with consistent compositions but varied particle sizes? 5. I didn’t find where \"Nearby Lair\" was explained or introduced."}, "author_response": "Responses to reviewers’ comments and changes we have made. Thank you very much for all the comments and suggestions from the reviewers, which have greatly improved our manuscript. Here we provided a point-to-point response and listed all new experiments and revisions in blue font. We believe that with all the efforts from our reviewers, the significance and novelty of our study has been robustly described and the revised manuscript would be suitable for publication. Reviewer #1 (Remarks to the Author) General Comment. The manuscript \"Accelerated discovery of highly active enzyme nanohybrids with parallelized hybrid-space Bayesian optimization\" by Jun Ge et al, is presenting an interesting combination of enzyme immobilization and the use of AI to efficiently sample the vast space of all possible experimental parameters to test out in order to arrive at more active catalyst. In general, the idea is excellent as there are not many studies which could address the rather \"trial-error\" area of enzyme immobilization where usually the through put of experimental set up is limited. However, the manuscript is completely lacking experimental evidence which would thoroughly prove the successful outcome of the work. I do not consider spectrophotometric assays with a surrogate/chromogenic substrate as quantitative proof for assessing activity. I would expect scale up on at least few mgs of substrate, isolating and characterizing the product. Also I did not find any experimental procedure for e.g. encapsulation efficiency, calculating activity recovery and reuse performance (how was the catalyst handled between each reaction set up). For those reasons, I can recommend the manuscript for publishing. Response: We sincerely thank your appreciation of the conceptual value of combining AI-based experimental design with enzyme immobilization field. Below, we would like to address your concerns point by point: 1. The method for enzyme activity assay The assays we used in our study are the standard method for testing enzyme activity, as can be seen in many previous studies (Nat. Commun. 2022, 13, 951., Nat. Catal. 2019, 2, 718., Nat. Commun. 2015, 6, 7240). During the activity measurement processes, the substrates we used are the real substrates for enzymes: ·For Glucose Oxidase (GOx), the activity is measured via the oxidation of the substrate-glucose, which can be further quantified through the chromogenic reaction with ABTS at 415 nm. ·Catalase (CAT) activity is measured by monitoring the decomposition of hydrogen peroxide. ·For CALB (lipase), the activity is assessed by tracking the release of p-nitrophenol from p-nitrophenyl esters which are also the standard substrates for lipase. These methods are both reliable and provide valuable initial insights into the catalytic activity of enzymes immobilized on various carriers. We believe these assays serve as a robust measure of enzyme performance and offer the advantage of high-throughput screening in early-stage research. 2. The catalytic reactions and product isolation/characterization Thank you for raising this point. Due to the high-throughput nature of our studies, we first focused on rapid screening of activity. Then, after finding the optimal carrier for each enzyme, we designed experiments to correlate enzyme activity with product conversion. Taking CAT as an example, define the amount of enzyme that decomposes 1 micromole of hydrogen peroxide per minute as 1 U. The measured enzyme activity data and relative enzyme activity in a scaled-up reaction system are listed in the table. Similarly, the enzyme activity data of GOx and CALB can be calculated. According to your valuable suggestions, we have added the corresponding descriptions in the revised manuscript. 3. The encapsulation efficiency, activity recovery, and reusability Thank you for your valuable suggestions. We have added the detailed calculation and test methods for the above properties in the revised Supporting Information. Your comments have helped us significantly improve the clarity and completeness of our manuscript. Revisions: Manuscript, Page 6: Once ideal carriers were found, detailed enzymatic properties including encapsulation efficiency and relative activity were investigated to get better understandings of the newly discovered nanohybrids. Manuscript, Page 18: A scaled reaction was performed to determine the catalytic properties of both free CALB and CALB@Zn(2,4-dmIM) , where the relative activity and activity recovery 2 showed similar results as above (Table S9). Manuscript, Page 19: Moreover, amorphous particles usually suffered from instability, easy breakage, and low reusability. In order to observe the tolerance of particles to mechanical shear force, the immobilized enzyme suspension was treated at high centrifugation of 13000rpm, 3 min, then resuspended, sampled, and tested for activity. Supporting Information, Page 3: Calculation of enzymatic properties At the initial stage of the work, activity recovery (AR) was selected as the most direct optimization target because this property contained comprehensive information of the entire immobilization process. The calculation formula was as follows: A apparent activity AR = ×100% A total activity wherein, the total activity was the specific activity of the free enzyme × the total amount of enzyme input in the immobilization reaction. After the iteration reached a certain stage, the immobilized material with higher AR was found, so its enzymatic properties were comprehensively measured. Encapsulation efficiency (EE) was the percentage of the amount of immobilized enzyme to the total amount of enzyme input. The amount of immobilized enzyme was obtained by subtracting the amount of enzyme in the supernatant after the reaction from the total amount of enzyme input. Q uantity immobilized enzymes EE = ×100% Q uantity input enzymes Relative activity (RA) was the proportion of residual enzyme activity in the immobilized material compared to the free enzyme at the same protein concentration. According to the encapsulation efficiency, the immobilized enzyme solution was prepared with the same concentration as the free enzyme, and the activity of both samples were tested under the same condition The calculation formula was as follows: A apparent activity RA = ×100% A free enzyme activity at same C Supporting Information, Page 10: Table S9. Comparison of prediction by PHBO and real experimental results. Enzyme activity (U) Enzyme Encapsula Relati Molar extinction concentration Activity ed in best efficiency activit coefficient (M-1 Free in the recovery carrier (%) y (%) cm-1) enzyme reaction (%) (μg/mL) 33000 (oxidized 0.118±0. 0.119±0.0 100.8 GOx 0.064 100.00 100.86 ABTS)[1] 001 02 6 5.614±0. 4.986±0.0 CAT 43.6 (H O )[2] 0.283 99.00 88.82 89.91 2 2 106 82 18000 (p- 6.150±0. 4.946±0.0 CALB 1.475 100.00 80.43 87.70 nitrophenol)[3] 130 68 Reaction: GOx: 3 mL reaction system, substrate was 47 mM glucose, containing 45 μg/mL HRP and 45 μg/mL ABTS. The final concentration of free enzyme and immobilized enzyme was 0.064 μg/mL. The degradation of glucose was determined by the change of absorbance at 415 nm within 15s of the production of oxidized ABTS, which was the activity of GOx. The amount of enzyme that decomposed 1 micromole of glucose and generated 1 micromole of oxidized ABTS per minute was defined as 1U. CAT: 3 mL reaction system, substrate was 9.8 mM hydrogen peroxide, the final concentration of free enzyme and immobilized enzyme was 0.283 μg/mL. The activity of CAT was determined by the decrease of absorbance at 240 nm (characteristic absorption peak of hydrogen peroxide) within 15s. The amount of enzyme that decomposed 1 micromole of hydrogen peroxide per minute was defined as 1U. CALB: 3 mL reaction system, substrate was 0.90 μg/mL p-NPB, the final concentration of free enzyme and immobilized enzyme was 1.475 μg/mL. The degradation of p-NPB was determined by the change of absorbance at 348 nm within 15s of the production of p-nitrophenol, which was the activity of CALB. The amount of enzyme that decomposed 1 micromole of p-NPB and generated 1 micromole of p-nitrophenol per minute was defined as 1U. [1] Cano et al. A method to measure antioxidant activity in organic media: application to lipophilic vitamins. Redox Report, 2000, 5, 365. [2] Nobel et al. The reaction of ferrous horseradish peroxidase with hydrogen peroxid. J. Biol. Chem. 1970, 245, 2409. [3] Bessey et al. Preparation and measurement of the purity of the phosphatase reagent, disodium p-nitrophenyl phosphate. J. Biol. Chem. 1952, 196, 175. Reviewer #2 (Remarks to the Author) General Comment. The authors present a novel PHBO method to optimize enzyme immobilization carriers based on MOFs. Their approach combines a customized Gaussian-process kernel designed for mixed categorical-continuous variable spaces, probabilistic reparameterization to enable gradient-based optimization, and a batch-suggestion strategy. The authors demonstrate significant improvement in activity recovery (AR) for immobilized enzymes compared to conventional methods. Transfer learning from GOx further reduced the optimization iterations for other enzymes. I found this work exciting, particularly because the proposed PHBO method creatively integrates multiple recent advancements in BO, such as hybrid kernel design, probabilistic reparameterization, and batch sampling. I would not hesitate to offer my support for its publication in Nat Comms and would like to suggest the following additions/corrections/etc. Comment 1. The authors may consider more clearly differentiating their method from established methods such as COMBO (Oh et al., NeurIPS 2019) and Gryffin (by Aspuru-Guzik and co-workers). Perhaps, the authors could include a comparative table summarizing algorithmic components (e.g., kernels, recommendation strategies, etc.) alongside expected computational or performance benefits. Response: We thank the reviewer for this insightful suggestion. To better clarify the novelty and strengths of our proposed PHBO algorithm, we have added a dedicated comparison of PHBO against existing Bayesian optimization methods designed for categorical or hybrid variable spaces, namely COMBO and Gryffin. This comparison highlights key differences in three core algorithmic components: surrogate model construction, acquisition function formulation and optimization, and parallelization. (1) Surrogate Model Construction: COMBO constructs a graph over the Cartesian product of discrete variables and defines a graph-based GP model using Laplacian eigenvectors as input features. Gryffin integrates prior expert knowledge into a weighted discrete domain using categorical priors, and then applies a one-hot encoding for GP modeling. In contrast, PHBO employs a probabilistic reparameterization strategy that maps categorical variables into a continuous distributional space. This allows the construction of a customized GP model directly in the reparameterized continuous space, reducing modeling complexity and minimizing dependence on domain-specific priors. (2) Acquisition Function and Optimization: COMBO and Gryffin both adopt standard acquisition functions (e.g., Expected Improvement, Upper Confidence Bound). COMBO performs optimization via graph traversal or enumeration, while Gryffin uses a rule-based heuristic search over the weighted categorical space. PHBO defines the acquisition function as the expectation over categorical distributions, which is estimated via Monte Carlo sampling. Its continuous and differentiable structure enables gradient-based optimization using stochastic methods. This formulation avoids the combinatorial explosion in COMBO when scaling to higher-dimensional categorical spaces and also mitigates Gryffin’s dependency on well-defined prior knowledge. (3) Parallelization: COMBO does not natively support parallelization. Gryffin achieves parallelism via diversified λ-weighted acquisition functions. PHBO leverages the Nearby Liar method, which iteratively adds pseudo-evaluated points (Liar points) with estimated objective values (Liar values) derived from nearby evaluated samples. Each recommendation of PHBO is based on an updated GP model, effectively reducing sample clustering and improving batch diversity. Revisions: Manuscript, Page 4: Gryffin20 performs well in hybrid spaces but relies heavily on prior domain knowledge, as its effectiveness strongly depends on the quality and relevance of user-provided descriptors. If these descriptors are poorly correlated with the optimization objective, the advantages of Gryffin may diminish significantly. Manuscript, Page 7: To address the challenges of modeling and optimization in hybrid categorical-continuous reaction spaces, we developed the Parallelized Hybrid-space Bayesian Optimization (PHBO) algorithm. PHBO integrates three key innovations: (i) a surrogate model built in a reparametrized continuous distribution space, which allows customized GP modeling; (ii) an acquisition function defined as an expectation over categorical variable distributions, enabling efficient gradient-based optimization via Monte Carlo gradient estimation; and (iii) the Nearby Liar parallelized sampling method, which incrementally updates the surrogate model to diversify batch recommendations and reduce redundant sampling. A comparison of PHBO with representative hybrid-space or categorical-space BO algorithms is summarized in Table 1. Table 1 Comparison of PHBO with existing BO algorithms for hybrid or categorical spaces COMBO19 Gryffin20 PHBO Expected Benefits of PHBO Surrogate Graph-based GP over One-hot encoded Customized GP on Lower complexity; model discrete product space GP with categorical reparameterized less reliance on using Laplacian priors distributional space expert priors eigenvectors Acquisition EI/UCB with graph EI/UCB with rule- UCB expectation over Scalable; enables function enumeration search based heuristics categorical variable stochastic gradient distributions optimization Parallelization Not supported Diversified λ Nearby Liar method Reduces sample heuristics clustering; efficient batch diversity Comment 2. I have several queries about the following essential details: Comment 2.1 How were kernel weights (Eq. 3) determined? Were these optimized, e.g., through maximum likelihood estimation or set manually? Response: We appreciate the reviewer’s question regarding the determination of kernel weights in (Eq. 3). In our implementation, the kernel weights were optimized through maximum likelihood estimation during the training of the GP surrogate model. This approach was adopted because it enables data-driven calibration of the relative contributions of the categorical and continuous components of the hybrid kernel. By maximizing the marginal likelihood of the observed data, this approach provides an efficient and principled way to tune the kernel weight, leading to improved optimization performance. This approach also avoids the need for manual tuning and allows the model to adapt to varying importance of variable types across different enzymes. Revisions: Manuscript, Page 8: The kernel weights were optimized using maximum likelihood estimation, allowing data-driven adjustment of the contributions from the categorical and continuous components. Comment 2.2 In the Monte Carlo sampling (Eq. 7), what was the number of samples used for the REINFORCE gradient estimator? Were variance-reduction techniques, such as baseline subtraction, employed? Response: We are grateful to reviewer’s insightful comment. As rightly pointed out, while the REINFORCE gradient estimator provided in (Eq. 7) ((Eq. 7a) in the revised manuscript) is unbiased, its variance can be substantial, especially when the underlying distribution is broad. To mitigate this issue and enhance the computational efficiency of the PHBO, we incorporated a variance-reduction strategy based on the use of control variates. Specifically, we employed the REINFORCE gradient estimator itself as a control variate by adjusting the score function with a proportional baseline derived from the historical mean of the acquisition function value, denoted as 𝛽. This results in the modified gradient estimator shown in Eq. (7b), where 𝛽 is defined as 𝐵 times the running average of 𝛼, and 𝐵 is a tunable hyperparameter controlling the extent of variance reduction. To balance accuracy and efficiency, we used 128 Monte Carlo samples per call to the REINFORCE estimator, which was found sufficient to achieve stable gradient estimation in practice. Revisions: Manuscript, Page 9: We have relabeled the original Eq. (7) as Eq. (7a), and added the following paragraph to the revised manuscript: “To reduce the variance of the REINFORCE estimator in (Eq. 7a), we adopted a control variate technique where the acquisition function 𝛼 is adjusted by subtracting a proportional baseline 𝛽, defined as 𝐵 times the running mean of 𝛼 from previous iterations. This yields a modified gradient estimator (Eq. 7b) with significantly lower variance. In practice, we used 128 Monte Carlo samples for each REINFORCE estimation, which represented a practical trade-off between variance reduction and computational efficiency.” 𝑁𝑀𝐶 𝜕 1 𝐸 [𝛼(𝑥,𝑧)] ≈ ∑(𝛼(𝑥,𝑧̃ )−𝛽)(𝐼 −𝑝(𝑧̃ |𝜃)) （7b） 𝜕𝜃 𝑧~𝑝(𝑍|𝜃) 𝑁 𝑖 𝑖 𝑀𝐶 𝑛=1 Comment 2.3 Please consider discussing or briefly reporting a computational complexity analysis. Response: We thank the reviewer for this valuable suggestion. We have now added a brief computational complexity analysis of PHBO to better contextualize its computational cost relative to its experimental efficiency. PHBO primarily involves two computational components: (1) training of the GP surrogate model and (2) optimization of the acquisition function. The complexity of each component is as follows: (1) GP model training. As described on Page 7 of the manuscript, PHBO employs GP as the surrogate model. The computational complexity for training the GP model is O(n3), where n is the number of observed sample points. This is consistent with the standard BO framework. (2) Acquisition function optimization. In PHBO, acquisition function optimization is performed in the reparameterized continuous space of probability distributions. The acquisition function is defined as an expectation over this distribution, which is estimated using Monte Carlo sampling. Let M be the number of samples per iteration, T the number of optimization iterations, the time complexity of acquisition function optimization is O(TMn2), as each Monte Carlo sample requires a GP-based prediction with complexity O(n) for mean and O(n2) for variance. To support batch recommendations, PHBO adopts the Nearby Liar strategy for parallel sampling. This method imposes minimal additional computational cost, as it leverages existing evaluations to simulate new observations and efficiently update the surrogate model. It is important to note that in typical experimental optimization scenarios, especially in laboratory-based synthesis and evaluation, the cost of running wet experiments far exceeds the computational cost of PHBO. Therefore, PHBO remains computationally efficient and well-suited for practical applications involving expensive experimental evaluations. Revisions: Manuscript, Page 10: The computational complexity of PHBO is primarily composed of two parts: GP model training and acquisition function optimization. The GP model training incurs a complexity of O(n3), where n is the number of observed samples. The acquisition function, estimated via Monte Carlo sampling over reparameterized distributions, has an optimization complexity of O(TMn2), where T is the number of optimization iterations, and M is the sample size. The parallelized sampling using Nearby Liar introduces negligible additional overhead. Given that experimental evaluations are typically more time-consuming than these computations, PHBO offers a favorable trade-off between decision-making efficiency and computational cost. Comment 3. The authors illustrate valuable transfer-learning capabilities from GOx to CAT and CALB, but the explanation remains qualitative. Please consider providing some quantitative metrics to evaluate transfer effectiveness. Response: We thank the reviewer for this constructive suggestion. We fully agree that assessing the effectiveness of transfer learning in a quantitative manner is important. However, due to the diversity of transfer learning strategies and the lack of a universally accepted metric to quantify source-target domain relatedness, there is currently no standardized evaluation protocol for transfer effectiveness in BO. To address this, we adopted a controlled comparative strategy to provide an intuitive evaluation of transfer learning performance in our system. Specifically, we compared the performance of PHBO with and without transfer knowledge for both CAT and CALB. The control setup used a “blank” PHBO initialization without GP prior parameter update nor BO warm-up, and generated one batch (8 conditions) of experimental recommendations for each enzyme. These results were then compared with the first iteration round of PHBO that included transfer knowledge from the GOx optimization. The results showed that the transfer-enabled PHBO achieved average AR of 71.31% for CAT and 28.49% for CALB, while the control baseline only achieved average AR of 28.74% for CAT and 4.39% for CALB. Therefore, the AR achieved in the first iteration with transfer was 148% and 549% higher than that of the blank control for CAT and CALB, respectively. These quantitative comparisons clearly demonstrate the positive effect of transfer learning in accelerating the identification of high-performance enzyme nanohybrids. Revisions: Manuscript, Page 15: Furthermore, to quantitatively assess the effect of transfer learning, we compared PHBO to a baseline version initialized without GP prior updates nor warm-up data. For both CAT and CALB, one batch (8 parallel conditions) was recommended by the control baseline and experimentally evaluated. The results presented in Table S6 and S7 showed that in the first iteration, the transfer-enabled PHBO achieved average AR of 51.16% for CAT and 9.77% for CALB, while the control baseline only achieved average AR of 28.74% for CAT and 4.39% for CALB. Such improvement validated the effectiveness of transfer learning and the benefit of incorporation prior knowledge from GOx optimization. Supporting Information, Page 9-10: Table S1. Comparison of the effectiveness of transfer learning between PHBO and control baseline: in the case of immobilized CAT biocatalysts. Ligand conc Molar ratio Reaction Activity Group Zinc type Zn conc (M) Ligand type (M) of ligand/Zn time (min) recovery (%) Zn(Ac)2 0.1550 2,4-dmIM 0.6250 40.3 30 85.7±6.1 Zn(NO3)2 1.0616 2,4-dmIM 0.5358 5.0 30 90.8±1.6 ZnBr2 1.0616 2,4-dmIM 0.5358 5.0 30 87.0±1.5 Zn(Ac)2 0.6134 PR 0.1658 2.7 30 79.7±1.3 PHBO Zn(Ac)2 0.2601 2-eIM 0.0100 0.4 30 54.9±2.2 Zn(Ac)2 0.6333 HTM 0.0341 0.5 30 85.8±1.2 Zn(Ac)2 0.5817 2-tIM 0.0814 1.4 30 47.9±0.3 Zn(Ac)2 0.5817 2-tIM 0.2327 4.0 30 38.7±0.8 Control ZnSO4 3.9591 2-mtIM 6.2296 15.7 50 14.7±0.5 baseline Zn(C3H5O3)2 0.0039 2-IMdn 3.4604 8872.8 105 1.9±0.5 ZnSO4 1.9194 2,4-dmIM 8.8332 46.0 55 0.0±0.0 Zn(H2PO4)2 0.2756 2-eIM 0.3678 13.3 40 54.7±2.3 Zn(Ac)2 0.0959 2-tIM 0.0674 7.0 80 0.6±0.6 Zn(Ac)2 0.6512 PR 0.5681 8.7 65 71.7±1.9 Zn(NO3)2 2.2209 2-IMdn 5.5182 24.8 40 38.1±1.3 Zn(Ac)2 0.6112 4-mIM 2.9172 47.7 120 48.3±0.6 Note: The control baseline experimental group was recommended without any prior knowledge. The conditions extensively explored the entire reaction space, with precursor concentrations reaching up to 8.8 M and a ligand-to-zinc ratio of 8872.8. The exploration efficiency was relatively low, resulting in an average AR of only 28.74%. In contrast, PHBO made recommendations based on knowledge extracted from the GOx experimental data. The conditions explored more of the low precursor concentration and low ligand-to-zinc ratio space, leading to an average AR of 71.31%. Table S2. Comparison of the effectiveness of transfer learning between PHBO and control baseline: in the case of immobilized CALB biocatalysts. Ligand conc Molar ratio Reaction Activity Group Zinc type Zn conc (M) Ligand type (M) of ligand/Zn time (min) recovery (%) Zn(NO3)2 1.0616 2,4-dmIM 0.5358 5.0 30 51.4±8.5 ZnBr2 1.0616 2,4-dmIM 0.5358 5.0 30 57.2±2.5 Zn(Ac)2 0.6134 PR 0.1658 2.7 30 32.3±1.5 Zn(Ac)2 0.5817 2-tIM 0.0814 1.4 30 17.7±0.5 PHBO Zn(Ac)2 0.5817 2-tIM 0.2327 4.0 30 31.9±1.1 Zn(Ac)2 0.6174 2-tIM 0.0752 1.2 30 13.7±0.8 Zn(Ac)2 0.5817 2-tIM 0.4654 8.0 30 8.9±0.5 Zn(Ac)2 0.2067 2,4-dmIM 0.8333 40.3 30 14.8±1.3 ZnBr2 8.1262 2-tIM 0.3613 0.4 50 18.3±1.7 ZnSO4 0.1058 2-tIM 0.0531 5.0 85 0.0±9.0 Zn(H2PO4)2 0.3748 2-mIM 1.1350 30.3 30 1.0±0.9 Control Zn(NO3)2 0.1811 PR 0.8154 45.0 100 11.2±0.7 baseline ZnBr2 9.3371 2,4-dmIM 1.6148 1.7 65 0.0±0.0 Zn(C3H5O3)2 0.0214 1-mOmIM 0.525 245.3 70 0.0±0.0 Zn(Ac)2 0.7879 TZ 4.9201 62.4 45 0.0±0.0 ZnSO4 1.0340 4-mIM 7.6677 74.2 60 4.7±0.9 Note: Similar results can be observed for CALB as well. The control baseline experimental group explored the reaction space more widely, with a ligand-to-zinc ratio ranging from 0.4 to 245.3. CALB faced difficulty in forming nanohybrids and relatively low AR. In this case, 3 out of 8 conditions recommended by control baseline failed to form participates. The average AR generated by PHBO and control baseline were 28.49% and 4.39%, respectively. Comment 4. The manuscript seems to only associate particle size and amorphousness with high AR. I failed to catch any discussion of a direct causal relationship. If these two mentioned factors are the only/most significant ones, then wouldn’t it make sense to conduct control experiments: synthesizing enzyme nanohybrids with consistent compositions but varied particle sizes? Response: We appreciate the reviewer’s valuable comments. Regarding the relationship between particle size, amorphous structure, and high activity retention (AR), we have partially addressed this issue with experimental data in the manuscript and would like to further clarify as follows: Due to the limitations of the aqueous in situ co-precipitation method, it is difficult to obtain particles with different sizes under exactly the same formulation. Therefore, we fixed the optimal precursor combination and achieved particles with varying sizes by adjusting the molar ratio of ligand to metal (Fig S9, S10, S11). Taking CAT as an example, at lower ligand-to-metal ratios, smaller and moderately aggregated particles were formed, corresponding to higher AR. As the ratio increased, particle size gradually increased while AR decreased (Fig 4d-f). This can be explained by the fact that smaller particles provide a larger surface area, enhancing enzyme exposure and substrate diffusion, thereby improving catalytic activity; whereas excessive ligand concentration raises the reaction system’s pH, which may impair enzyme activity. A similar trend was observed for GOx. As for CALB: at very low ligand-to-metal ratios (<4.0), particle morphology was irregular and immobilization efficiency was poor, yielding limited meaningful data. At moderate ratios, particle sizes were smaller, and AR reached a maximum, further supporting the above mechanistic interpretation. In summary, our experimental results demonstrate the direct and significant impact of particle size and carrier structure on enzyme activity. Corresponding revisions have been made in the main text to clarify this point. Revisions: Manuscript, Page 16: Fig. 4: Activity recovery, encapsulation efficiency, and relative activity of (d) GOx (Zn(Ac) and 2-tIM), (e) CAT (Zn(Ac) and HTM), (f) CALB (Zn(NO3) and 2,4- 2 2 2 dmIM) synthesized at different molar ratios. Manuscript, Page 18: As discussed earlier in Fig 4d-f and S9-S11, by adjusting the ligand-to-metal molar ratio under the same precursor combination, we obtained enzyme nanohybrids with different particle sizes. When the ligand-to-metal molar ratio was too low (e.g., 0.1–4.0 for CALB), it was difficult to form well-defined particles. Once stable particles were obtained, a clear trend emerged: the activity retention (AR) decreased as particle size increased. Under optimal conditions, the resulting particles exhibited negligible porosity, suggesting that enzymes were loosely confined within the solid matrix rather than encapsulated within a porous framework. In this context, smaller particles offered a higher surface area-to-volume ratio, enhancing enzyme exposure to substrates and thus contributing to the observed improvement in catalytic activity. Supporting Information: Figure S9. SEM images of GOx hybrid biocatalysts synthesized with Zn(Ac)2 and 2-tIM with different molar ratios. Figure S10. SEM images of CAT hybrid biocatalysts synthesized with Zn(Ac)2 and HTM with different molar ratios. Figure S11. SEM images of CALB hybrid biocatalysts synthesized with Zn(NO3)2 and 2,4-dmIM with different molar ratios. Comment 5. I didn’t find where \"Nearby Lair\" was explained or introduced. Response: We appreciate the reviewer for pointing out this oversight. The Nearby Liar method was briefly mentioned at the bottom of Page 8 in the original manuscript (“To parallelize sampling, inspired by the Constant Liar method…”), but we agree that its mechanism was not sufficiently described. To clarify, Nearby Liar belongs to the family of Liar methods, which are designed to enable batch acquisition in BO. These methods decompose parallelized sampling into a sequence of steps: after a sample point is suggested, it is added to the training dataset as a “Liar point” along with a “Liar value”, which is a fabricated objective value. The surrogate model is then updated based on this temporary value. This process is repeated until the desired number of sample points is collected, after which physical experiments are conducted and the Liar values are replaced with actual observations. In contrast to the Constant Liar method, Nearby Liar improves the accuracy of the temporary evaluation by assigning the average objective value of several previously evaluated sample points that are nearest to the new sample point in the chemical space. This more accurate estimation helps reduce over-exploration and leads to more efficient parallel sampling. Revisions: Manuscript, Page 9: We have modified the last paragraph in this page in the revised manuscript: “To parallelize sampling, inspired by the Constant Liar method31, we developed the Nearby Liar method. In Liar methods, parallelized recommendations are sequentially generated by adding each newly suggested sample point to the training set with a temporary “Liar value”, that is, a fabricated evaluation used to update the surrogate model before actual evaluation. In Nearby Liar, instead of assigning a fixed value, the Liar value is set to the average objective value of the closest previously evaluated sample points in the reaction space, improving estimation accuracy and reducing redundant sampling. after suggesting a sample point, PHBO does not immediately conduct the experiment but uses the Nearby Lair method to continue suggesting new sample points. Specifically, PHBO selects a series of sample points that are closest to the suggested sample point from the evaluated dataset:”"}, {"round": 2, "reviewer_comments": {"reviewer_2": "My comments have been fully addressed and I don't have further suggestions.", "reviewer_3": "The authors introduce the application of AI tools to sample a set of many experimental parameters affecting the immobilization of enzymes on their activities, aiming to discover optimized experimental conditions for biocatalytic transformations. While the concept and approach are interesting, and could advance the area of biocatalysis beyond the “trial-error” method, this reviewer, being an experimentalist, rather than an AI specialist, has serious doubts regarding the overall value of the paper. 1. The selection of enzymes adapting the method is not convincing, neither are the outcoming results. The authors selected very active highly stable enzymes (almost immobilization insensitive) as test cases. Glucose oxidase, catalase and the lipase- CALB. Standard conditions given in any catalog would yield comparable activity results without any AI-predictions. If so, what are the advantages of AI? This reviewer would like to see one or two examples of enzymes that upon immobilization under “standard” conditions lead a low activity (say 20% of the native homogeneous enzyme) and a significantly enhanced activity (60-80)% upon AI-predicted conditions. Sensitive enzymes, such as carbonic anhydrase, beta-galactosidase, formate dehydrogenase, uricase, should have been tested. 2. The enzyme assays used by the authors to probe their AI predicted results implemented non-native agents (ABTS, p-nitrophenol). Monitoring “real” biological products (CO2 for formate dehydrogenase, bicarbonate for carbonic anhydrase or galactose for beta-galactosidase would be experimental convincing examples. 3. As pointed us by the first-round reviewers, the immobilization efficacy is not only valued by enhanced activities but, also by other environmental factors, such as temporal loss of biocatalytic activity or scaleability. This reviewer would like to see an example where ‘standard” immobilization condition yield an unstable time-dependent system, whereas the AI predicted conditions lead to enhanced temporal stabilities."}, "author_response": "Responses to reviewers’ comments and changes we have made. Thank you very much for all the comments and suggestions from the editor and reviewers, which have greatly improved our manuscript. Here we provided a point-to-point response and listed all new experiments and revisions in blue font. We believe that with all the efforts from our editor and reviewers, the significance and novelty of our study has been robustly described and the revised manuscript would be suitable for publication. Reviewer #2 (Remarks to the Author) My comments have been fully addressed and I don't have further suggestions. Response: We sincerely thank you for recognizing and supporting our first round of revisions. Our manuscript has been significantly improved due to your review. Reviewer #3 (Remarks to the Author) Comment 1: The authors introduce the application of AI tools to sample a set of many experimental parameters affecting the immobilization of enzymes on their activities, aiming to discover optimized experimental conditions for biocatalytic transformations. While the concept and approach are interesting, and could advance the area of biocatalysis beyond the “trial-error” method, this reviewer, being an experimentalist, rather than an AI specialist, has serious doubts regarding the overall value of the paper. 1. The selection of enzymes adapting the method is not convincing, neither are the outcoming results. The authors selected very active highly stable enzymes (almost immobilization insensitive) as test cases. Glucose oxidase, catalase and the lipase- CALB. Standard conditions given in any catalog would yield comparable activity results without any AI-predictions. If so, what are the advantages of AI? This reviewer would like to see one or two examples of enzymes that upon immobilization under “standard” conditions lead a low activity (say 20% of the native homogeneous enzyme) and a significantly enhanced activity (60-80)% upon AI-predicted conditions. Sensitive enzymes, such as carbonic anhydrase, beta-galactosidase, formate dehydrogenase, uricase, should have been tested. Response: We thank the reviewer for this comment. Regarding the use of more sensitive enzymes, we would like to respectfully clarify that GOx, CALB, and CAT should not be regarded as “insensitive” enzymes for immobilization. Many publications have demonstrated that conventional immobilization methods and commonly used carrier materials frequently lead to low relative activities or limited encapsulation efficiencies for these enzymes. For example, in Biosensors & Bioelectronics 6 (1991) 663, GOx immobilized on a glutaraldehyde-crosslinked Pt electrode retained only approximately 24.6% of the native enzyme activity (24,036 vs. 97,822 units). In Nature Communications 13 (2022) 951, even when using ZIF-8 as the carrier, different crystallization pathways yielded markedly divergent outcomes, with GOx@BZIF-8-S retaining only 3-5% of the free-enzyme activity. In Biochemical Engineering Journal 112 (2016) 20, GOx immobilized on multi-walled carbon nanotubes exhibited a relative activity of merely ~20%. For CALB, Journal of Nanomaterials (2021) 8812240 reported that immobilization on silica nanoparticles with varying surface areas and hydrophobicity produced relative activity retention ranging widely from 12.5% to 100%. Likewise, in Langmuir 40 (2024) 16338, CAT immobilized on unoptimized plastic nanobeads retained only ~22.2% activity. Consistent with these prior observations, our own measurements using ZIF-8, one of the most widely adopted MOF carriers, resulted in low relative activities for all three enzymes (27.9%, 62.3%, and 7.6% for GOx, CAT, and CALB, respectively, data shown in Figure S5). These results showed that these three enzymes used in our work still suffered from substantial activity loss under many conventional immobilization conditions. Therefore, further exploration was required to identify the most suitable immobilization carriers for each enzyme. Another reason for selecting GOx, CALB, and CAT was that these enzymes were among the most extensively used ones in industrial biocatalysis, with well-established enzymatic assays. These features make them ideal benchmarks for evaluating whether a new AI-based optimization method is accurate, reproducible, and generally applicable. Starting from low relative activity obtained with baseline ZIF-8, our AI-guided workflow was able to identify optimized conditions that increased the relative activity to 100%, 88.8%, and 80.4%, respectively (data shown in Table S9). This significant improvement proved the effectiveness and the general scope of the method. In addition, we fully acknowledge the reviewer’s interest in more fragile enzymes. To address this point, we have conducted supplementary experiments on tyrosinase (Tyr), β-galactosidase (GAL), and carbonic anhydrase (CA). In these cases, we performed a small number of iterative rounds (three cycles) to identify the most promising carriers for each enzyme. Only the initial and final results were reported here, as the intermediate formulations are part of ongoing studies regarding the applications of these enzymes to be published separately. Even with limited optimization, we observed substantial improvements in both encapsulation efficiency and relative activity compared with the initial ZIF-8 baseline. Specifically, the relative activity of Tyr, GAL, and CA immobilized with ZIF-8 were 19.29%, 0.16%, and -0.03%, respectively, while those immobilized with their respective optimized carrier increased to 77.57%, 47.06%, and 88.82%, respectively (Figure S14). These preliminary results further support the reliability and extensibility of this AI-guided workflow across a broader range of enzymes. At the same time, more comprehensive optimization and application-oriented studies on these enzymes are currently ongoing in our laboratory and will be reported in future publications, as they extend beyond the scope of the present methodological study. Revisions: Manuscript, Page 20: Beyond these primary model systems, preliminary validation on additional enzymes including on tyrosinase (Tyr), β-galactosidase (GAL), and carbonic anhydrase (CA) further illustrated the extensibility of the workflow to broader biocatalysts without requiring major methodological adjustments (Fig S14). Supporting Information, Page 27: Figure S14. (a) Encapsulation efficiency and (b) relative activity of different enzymes immobilized by ZIF-8 and respective best carriers. Data were represented as mean ± SD (n=3). Respective best carrier: For each enzyme, PHBO was initiated from ZIF-8 and iteratively screened 8 carriers per round for three rounds. The carrier yielding the highest activity recovery was selected as the “best carrier” and its encapsulation efficiency and relative activity were subsequently measured. Even with limited optimization, substantial improvements were observed in both encapsulation efficiency and relative activity compared with the initial ZIF-8 baseline. Specifically, the relative activity of Tyr, GAL, and CA immobilized with ZIF-8 were 19.29%, 0.16%, and -0.03%, respectively, while those immobilized with their respective best carrier increased to 77.57%, 47.06%, and 88.82%, respectively. 2. The enzyme assays used by the authors to probe their AI predicted results implemented non-native agents (ABTS, p-nitrophenol). Monitoring “real” biological products (CO2 for formate dehydrogenase, bicarbonate for carbonic anhydrase or galactose for beta-galactosidase would be experimental convincing examples. Response: We thank the reviewer for this comment. The three enzymes evaluated in our original submission (GOx, CAT, and CAL) were evaluated using their native or native-like substrates. Specifically: · GOx activity was measured using glucose (its natural substrate), with ABTS serving only as a chromogenic probe for the H O generated during the native 2 2 oxidation reaction. · CAT activity was directly measured from the decomposition rate of its natural substrate, H O . 2 2 · CALB (lipase) activity was measured using p-nitrophenyl esters (p-NPB), which were widely accepted native-like substrates that mimic natural ester hydrolysis for lipase activity evaluation in enzyme immobilization studies. As we noted in our response to Reviewer #1 in the last revision, the activity assays used in this study were standard and widely adopted methods in the enzyme immobilization literature (Nat. Commun. 2022, 13, 951., Nat. Catal. 2019, 2, 718., Nat. Commun. 2015, 6, 7240). Because our work aimed to establish a new AI-assisted optimization workflow which required screening hundreds of immobilization conditions across a high-dimensional parameter space, it was essential to employ activity assays that were high-throughput, reproducible, and readily standardized. In this context, colorimetric assays provide a reliable and literature-validated proxy for preliminary activity evaluation. After identifying the optimal carrier materials for each enzyme, we further conducted more detailed analysis including morphology, relative activity, and encapsulation efficiency. Additionally, we evaluated the catalytic performance toward each enzyme’s native substrate. These results have been provided in the Supporting Information Table S9. Similarly, the other three enzymes (Tyr, GAL, and CA) were evaluated likewise based on their native or physiologically relevant substrates. Specifically: · Tyr activity was measured using its natural substrate L-Dopa, which was enzymatically oxidized to dopa-quinone and subsequently reacted with MBTH to form a chromogenic product detectable at 505 nm. · GAL activity was measured using p-nitrophenyl-β-D-galactopyranoside (pNPGal), a chromogenic analog of lactose that retained the same β-1,4-galactosidic linkage that could be hydrolyzed by β-GAL. pNPGal was a standard galactosidic substrate whose cleavage yield p-nitrophenol with a characteristic absorbance at 400 nm. · CA activity was measured using its native substrate CO by Wilbur-Anderson assay, 2 which was a classical and widely adopted method for measuring CA activity. In this case, the accelerated pH drop reflected the rapid formation of bicarbonate and the associated release of protons, thereby directly reporting catalytic turnover toward the natural substrate. Revisions: Supporting Information, Page 3: Tyrosinase (Tyr) activity was measured using a commercial assay kit (ACMEC, AC10628). Tyr catalyzed the oxidation of L-Dopa to dopa quinone, which subsequently reacted with MBTH to produce fuchsia substance with absorbance at 505 nm. Reactions were performed according to the manufacturer’s instructions. β-Galactosidase (GAL) activity was measured using a commercial assay kit (Solarbio, BC2580). GAL hydrolyzed p-nitrophenyl-β-D-galactopyranoside (pNPGal, a standard galactosidic substrate) to yield p-nitrophenol, with a characteristic absorbance at 400 nm. Reactions were performed according to the manufacturer’s instructions. CA activity was measured using the Wilbur-Anderson method based on pH change. Briefly, a 20 mM TAPS buffer (solution A) and CO -saturated water (solution B) were 2 prepared in advance. For the blank measurement, 2.45 mL of solution A and 50 μL of DI water were added to a 10 mL tube under magnetic stirring, and the time required for the pH to drop from 8.3 to 7.0 after adding solution B was recorded (40 ± 2 s). For enzyme measurements, the same procedure was followed except that 50 μL of free or immobilized CA solution (diluted to 15 μg/mL) replaced the water. The decrease in time relative to the blank reflects catalytic activity. 3. As pointed us by the first-round reviewers, the immobilization efficacy is not only valued by enhanced activities but, also by other environmental factors, such as temporal loss of biocatalytic activity or scaleability. This reviewer would like to see an example where ‘standard” immobilization condition yield an unstable time-dependent system, whereas the AI predicted conditions lead to enhanced temporal stabilities. Response: We thank the reviewer foe this comment. We acknowledge the reviewer’s concern on the temporal and environmental stability. In this study, our optimization objective was enzyme activity, and therefore the AI-derived conditions were not explicitly directed toward maximizing stability. Nevertheless, we conducted a preliminary experiment using carbonic anhydrase (CA), a representative enzyme that required tolerance to high-temperature and strongly alkaline amine solvents in industrial chemical absorption process. N-methyldiethanolamine (MDEA) was a widely used organic amine solvent in chemical absorption of CO -capture process, 2 whose adsorption mechanism could synergize well with CA, as the bicarbonate generated by the enzyme was a key intermediate in MDEA-based carbon capture (Chemical Engineering Journal 2017, 307, 776). Therefore, we evaluated CA in 30% MDEA solution at 80 °C, which was a deliberately defined condition that better reflected real industrial process. Because the ZIF-8 formulation used in our workflow did not preserve CA activity well (Figure S14), we adopted a literature-reported ZIF-L synthesis with higher baseline activity as the “standard immobilization” control (Separation and Purification Technology 2023, 315, 123683). Under the heated and alkaline condition, the standard immobilized CA rapidly lost activity, whereas the AI-identified best carrier maintained 32.3% activity after 6 h of heating. These results showed that, even though stability was not part of the optimization target, the AI-guided carrier nevertheless conferred enhanced resistance to harsh conditions. We fully agree that stability is crucial for industrial biocatalysis. Importantly, such parameters can be directly set as optimization objectives in the PHBO workflow. This provides a clear path for future work, where long-term stability or tolerance to harsh environments can be directly set as the optimization function to find the most suitable carriers to different conditions and demands. Revisions: Manuscript, Page 20: Moreover, although catalytic activity was the optimization objective in this study, the AI-identified carriers also exhibited improved stability relative to conventionally prepared materials on the case of CA (Fig S15). In the future work, different properties could be directly defined as the optimization objectives, for example, long-term stability or tolerance to harsh environments, enabling the algorithm to discover carriers tailored to specific operational conditions and application needs. Supporting Information, Page 28: Figure S15. Thermal stability of CA immobilized on different carriers in 30% MDEA solution at 80 °C. N-methyldiethanolamine (MDEA) was a widely used organic amine solvent in chemical absorption of CO -capture process. 2 Relative activity was monitored over time to compare the stability of CA immobilized by AI-identified carriers and control formulations. The ZIF-L carrier was chosen as the control[1] and synthesized following the procedure reported in the reference. Under the heated and alkaline condition, CA@ZIF-L rapidly lost activity, whereas the AI-identified best carrier maintained 32.3% activity after 6 h of heating. Data were represented as mean ± SD (n=3). [1] Peijing Shao, et al. Shape controlled ZIF-8 crystals for carbonic anhydrase immobilization to boost CO2 uptake into aqueous MDEA solution. Separation and Purification Technology, Volume 315, 2023, 123683."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors made a serious effort to respond to the reviewer’s comments. Nevertheless, the original doubts and critical comments are still in place, and have not been convincingly answered in the revised paper. The authors have to admit that all the enzymes mentioned in the paper are highly stable, highly active, and can be easily immobilized. None of the enzymes reported on by the authors demonstrate the “revolutionary” AI optimized activity that could otherwise not have been achieved without the AI platform. In fact, AI-driven achievements reported by the authors could have been easily achieved by simple chemical steps available to chemists, who are active in enzyme biotechnology. The authors should concede that it is easy to kill the activity of many of the enzymes reported in the study and rescue the killing/degradation by AI recipes (but also by simple chemical steps). Nevertheless, the concept introduced by the authors is interesting, scientifically logical and legitimate. The authors made efforts to improve the paper, and this reviewer is convinced that they believe that their attitude can serve as a powerful tool for advancing enzyme biotechnology. Accordingly, the reviewer recommends accepting the paper. The authors, however, are advised to “tone-down” the “revolutionary” aspect of the concept in the conclusion section, while explicitly stating the weaknesses of the present study: Implementation of highly active and stable enzymes, use of artificial synthetic non-native probes to transducing the catalytic process and lack of a unique example achievable by AI that cannot be achieved by conventional biotechnological methods."}, "author_response": "Response: We thank the reviewer for the appreciation of the conceptual interest and scientific soundness of the study. In the meantime, we have revised the Discussion part according to the reviewer’s suggestion to better present the scope of the work in a more balanced manner. We now explicitly state the current limitations of the study, including enzyme selection, probe design, and the absence of AI-exclusive examples. We believe that these revisions provide a more accurate and transparent representation of the present study and better align the manuscript with the reviewer’s constructive suggestions. We acknowledge that the present study focuses on a limited set of enzymes that are relatively well studied and, in some cases, commercially relevant. Compared to enzymes that are newly discovered, derived from native sources, or exhibit highly specialized catalytic functions, these enzymes are considered more robust in an application context. However, many studies have shown that the choice of immobilization materials and structures can substantially influence the catalytic performance of these enzymes. Our method accelerates the finding process of suitable immobilization materials. The strength of the workflow is its ability to efficiently navigate high-dimensional hybrid design spaces, reduce empirical trial-and-error, and identify favorable parameter regions with substantially fewer experimental iterations. In this context, the AI workflow serves as a practical tool for accelerating optimization and improving experimental efficiency. In fact, ongoing work in our group is extending this workflow to additional enzymes and more unique examples. We hope that the present study can serve as a methodological starting point and inspire broader application of this approach to a wider range of enzymes. Revisions: Manuscript, Page 20: While the present study focused on a limited set of highly active and stable enzyme systems and employs artificial probes to transduce catalytic performance, it served as an initial demonstration of how data-driven optimization can be integrated into enzyme immobilization design. The favorable design regions identified here could, in principle, also be discovered through conventional chemical or biochemical exploration, which on the other hand typically relied on extensive trial-and-error and became increasingly inefficient as the parameter space expands. Therefore, the primary role of the PHBO-guided workflow is to accelerate the navigation of complex, high-dimensional design spaces with substantially reduced experimental effort. Future studies may define different properties, such as long-term stability or tolerance to harsh environments, enabling the algorithm to discover carriers tailored to specific operational conditions and application needs. Overall, the algorithm's adaptability to diverse materials, reaction spaces, and optimization parameters suggest its utility as an efficient method for developing high-performance biocatalytic composites."}]

Decision letter