A method for screening mutation sites of enzyme enantioselectivity based on conformation flux

Abstract

The application provides a method for screening mutation sites of enzyme enantioselectivity based on conformation flux. The method quantifies enzyme molecular channel recognition, substrate inlet sampling, active pocket conformation pre-organization, near attack conformation recognition, non-productive pocket stagnation analysis, ligand escape path statistics and gate switching detection, and outputs enantioselectivity mutation site priority and recommended mutation direction.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and relates to molecular dynamics simulation and computer-aided protein design, and more specifically to an enzyme enantioselective mutation site screening method based on conformational flux. Background Technology

[0002] Enzyme-catalyzed asymmetric resolution and asymmetric synthesis are important technical routes for the preparation of chiral compounds. Compared with traditional chemical catalysis, enzyme catalysis usually has advantages such as mild reaction conditions, high regioselectivity and stereoselectivity, and environmental friendliness. Therefore, it is widely used in the preparation of pharmaceutical intermediates, fine chemicals, fragrances, pesticide intermediates, and functional material monomers.

[0003] In enantioselective enzyme catalysis, the same enzyme often needs to select between a pair of enantiomers with highly similar structures. Since the two enantiomers have the same molecular composition and only opposite spatial configurations, the enzyme's selective recognition depends not only on the static binding capacity of local residues in the active site, but also on multiple factors such as substrate entry pathway, channel bottleneck, side cavity shunting, conformational pre-organization, near-attack conformational formation, and gating dynamic switching.

[0004] For many esterases and lipases, their active sites are not fully exposed on the protein surface, but rather located inside the protein or in a semi-buried pocket. After the substrate enters the enzyme surface from the solution phase, it must pass through the surface capture region, the channel entrance, the middle section of the channel, and the vestibular region of the catalytic cavity before it can finally form a near-attack conformation that meets the catalytic requirements. Simultaneously, the substrate may also be captured by the non-productive side cavity, forming a stagnant state far from the catalytic residues, thus preventing effective catalysis. Therefore, the essence of enzyme enantioselectivity is not simply a difference in binding energy, but a difference in conformational flux distribution between productive and non-productive pathways.

[0005] Existing rational design methods for enzymes mainly include static structure-based active pocket analysis, molecular docking, molecular dynamics simulation, channel geometry analysis, energy calculation, and machine learning prediction. While each method has its advantages, they still have significant limitations when used for screening enantioselective mutation sites.

[0006] First, traditional molecular docking is mainly used to obtain the possible binding conformations and scoring results of substrates in the active pocket. However, it is usually based on static structure and cannot reflect the process of substrates entering the protein surface from the solution phase, entering the active site through the channel, and forming or losing productive conformations during dynamic gating. For a pair of enantiomers with highly similar conformations, docking energy or static binding posture alone is often insufficient to accurately explain the enantioselectivity differences in experiments.

[0007] Second, while conventional molecular dynamics simulations can describe the dynamic behavior of proteins and substrates, substrate entry into channels, release from side cavities, gating switching, or escape from active cavities are typically rare events. In conventional molecular dynamics simulations at the hundreds of nanosecond scale, partial conformational changes and localized motion patterns may be observed, but it is difficult to stably and sufficiently replicate the complete substrate entry or product release pathway. Relying solely on microsecond-scale or even longer-scale molecular dynamics simulations results in high computational costs, making them unsuitable for scenarios involving multiple mutants, multiple substrates, or high-throughput screening.

[0008] Third, channel analysis tools can identify potential substrate entry and exit channels in enzyme structures and output channel length, bottleneck radius, centerline, curvature, and channel clustering results. These methods can answer the question of which channels might exist in an enzyme structure, but they cannot directly answer questions such as which channel is productive for a particular enantiomer, which channel causes non-productive shunting, or which residues determine conformational preferences and selective flipping. In other words, the geometric presence of a channel does not equate to its functional productivity.

[0009] Fourth, methods such as stochastic accelerated molecular dynamics and stochastic external force escape simulation can accelerate ligand escape from the buried pocket and be used to infer possible exit paths or ligand residence times. However, existing methods usually focus on the direction and time of ligand exit from the pocket, and do not integrate entry sampling, exit statistics, near-attack conformation, unproductive side cavities, and gating switching into a unified enantioselectivity prediction framework, nor can they directly output the priority of candidate residues that can be used for experimental mutation construction.

[0010] Fifth, machine learning methods can be used to predict enzyme activity, mutation effects, or kinetic parameters, but they rely on a large amount of high-quality training data. Publicly available data are often insufficient for specific enzymes, specific chiral substrates, and specific channel gating mechanisms; at the same time, model predictions usually lack clear structural and mechanistic explanations, making it difficult to explain why a mutation at a certain site increases R conformation selectivity or causes selectivity reversal.

[0011] Therefore, existing technologies lack a unified quantification of substrate entry, channel migration, active pocket pre-organization, near-attack conformation formation, non-productive side cavity capture, exit path and gating dynamic switching, and a technical solution to transform these dynamic behaviors into candidate mutation site ranking results. Summary of the Invention

[0012] This invention aims to address the following problems in existing technologies: the inability of static molecular docking and pocket analysis to describe the dynamic migration, diversion, and residence of enantiomeric substrates in enzyme channel networks; the limitation that channel geometry analysis can only identify candidate channels and cannot distinguish between productive and non-productive channels; the difficulty of conventional molecular dynamics simulations in fully capturing rare events such as substrate entry, side cavity retention, gating switching, and exit escape within a limited time; the difficulty of directly converting stochastic accelerated molecular dynamics results into recommendations for enantioselective mutation sites; and the problems of large mutation libraries, high experimental screening costs, and insufficient explanation of candidate site mechanisms in enantioselective mutation screening. To solve these technical problems, this invention provides a method for screening enzyme enantioselective mutation sites based on conformational flux. This method unifies and quantifies enzyme molecular channel identification, substrate entry sampling, active pocket conformation pre-organization, near-attack conformation identification, non-productive pocket stagnation analysis, ligand escape path statistics, and gating switching detection, and outputs the priority of enantioselective mutation sites and recommended mutation directions accordingly.

[0013] The technical solution adopted in this invention is: a method for screening enzyme enantioselective mutation sites based on conformational flux, comprising the following steps:

[0014] Obtain the three-dimensional structure of the target enzyme, the enantiomeric substrate to be compared, and the reference residues of the active site of the target enzyme;

[0015] Based on the three-dimensional structure of the target enzyme, candidate channels leading from the active site to the protein surface are identified, and channel data for each candidate channel is generated.

[0016] Based on channel data, ingress sampling simulation, binding state repeated sampling simulation, and piecewise random external force escape simulation are performed for each enantiomer to obtain ingress sampling simulation data, binding state repeated sampling simulation data, and piecewise random external force escape simulation data, respectively. State discretization is then performed to construct a conformational state diagram and calculate the conformational flux. The conformational state diagram includes at least one or more of the following: bulk state, surface capture state, channel ingress state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and exit state. The conformational flux includes at least one or more of the following: ingress flux, near-attack conformational flux, non-productive residence penalty term, exit flux, and gating score.

[0017] Based on the difference in conformational flux between enantiomers, an enantioselectivity score and a site priority score are calculated, and at least one of the following is output: candidate mutation site, recommended mutation direction, and enantioselectivity prediction.

[0018] This invention proposes the concept of conformational flux. Conformational flux refers to the frequency, probability, and directional combination of substrate transfer between different functional states, such as the enzyme molecule surface, channels, active cavity, and lateral pockets. Unlike traditional observations based on static binding energy, single channel radius, or single trajectory, conformational flux emphasizes the dynamic flow of substrate along the path of "bulk solvent—surface capture—channel inlet—active pocket pre-organization—near-attack conformation or non-productive pocket—outlet." For a pair of enantiomers R and S, even if both can enter the same enzyme pocket, their conformational flux distributions may be completely different. For example, the R-configuration substrate may have a higher main channel inlet flux and a higher near-attack conformation flux, while the S-configuration substrate may have a higher non-productive lateral cavity retention flux. In this case, the enzyme exhibits R-selectivity not because the S-configuration substrate cannot enter the protein at all, but because the S-configuration substrate, after entering, is diverted to a non-productive state and cannot effectively form a catalytic conformation. Therefore, this invention no longer solely determines whether the substrate can bind or whether a channel exists, but rather determines whether the substrate reaches the catalytic state along the productive conformation flux, and further identifies which residues control this flux allocation. This is the core innovation of this invention, distinguishing it from traditional molecular docking, single-channel analysis, and simple RAMD escape statistics.

[0019] Preferably, the channel data includes at least centerline data, bottleneck radius data, and ROI residue set.

[0020] Preferably, the inlet sampling simulation data includes at least the surface capture pit, the main inlet channel, and the intermediate metastable region.

[0021] Preferably, the binding state repeated sampling simulation data includes at least the trajectory data of the substrate in the active pocket pre-organization, near-attack conformation state, and non-productive side cavity state.

[0022] Preferably, the segmented random external force escape simulation data includes at least the escape path of the substrate or product from the active cavity to the protein surface, the escape time, and the probability of use of the exit channel.

[0023] Preferably, the inlet sampling simulation is a ligand replication flooding molecular dynamics simulation, which includes: arranging multiple substrate copies in a predetermined spatial shell outside the target enzyme solvent-exposed surface, and statistically analyzing events of substrate capture by the surface, entry into the channel ROI, or entry into the active pocket during molecular dynamics. In some embodiments, the number of substrate copies is 10 to 100, preferably 20 to 60; the closest distance between the initial substrate placement region and the protein surface is 0.8 to 3.0 nm, preferably 1.2 to 2.0 nm; and the minimum distance between any two initial heavy atoms of the substrate is 0.6 to 1.2 nm, preferably 0.8 nm.

[0024] Preferably, the segmented random external force escape simulation uses a distance-contact number dual criterion as the escape criterion. The distance-contact number dual criterion simultaneously satisfies the following conditions: the robust statistical value of the distance between the substrate or product centroid and the active site reference centroid in the tail window is greater than a preset distance threshold, and the robust statistical value of the number of heavy atom contacts between the substrate or product and the active site reference centroid in the tail window is less than a preset contact number threshold. The robust statistical value in the tail window is selected from the median, truncated mean, or estimated value under median absolute deviation constraints, preferably the median of 3 to 11 sampling points in the tail. In some embodiments, the segmented random external force escape simulation includes: for each repeated trajectory, applying a random directional external force to the ligand centroid in segments of 50 to 500 ps; updating the direction after each segment and using the last frame of the previous segment as the starting point of the next segment, until the exit criterion is met or the maximum number of segments is reached; the number of repetitions is 10 to 200, preferably 20 to 50.

[0025] Preferably, the criteria for the near-attack conformation state include at least two or more of the following four: a distance threshold between the catalytic nucleophile and the substrate reaction center atom; a nucleophilic angle of attack threshold; a hydrogen bond threshold between the substrate carbonyl oxygen and the oxygen anion hole-donating residue; and a spatial orientation threshold of the substrate relative to the catalytic triplet. In some embodiments, the distance threshold is preferably 2.5 to 3.5 Å, the angle of attack threshold is preferably 90 to 120 degrees, and the hydrogen bond donor-acceptor distance threshold is preferably no greater than 3.5 Å.

[0026] Preferably, the criteria for the non-productive pocket state include: the ligand center is located within a predefined lateral cavity ROI, the continuous residence time reaches a threshold, the number of contacts with the shell residues of the lateral cavity ROI is satisfied, and at least one of the near-attack conformation criteria is not satisfied.

[0027] Preferably, the gating score is determined by combining at least two types of information: dynamic cross-correlation changes between the catalytic residue region and the channel ROI region; time-series changes in the bottleneck radius of the candidate channel; state transitions of the side chain conformation angle, main chain dihedral angle, or side chain conformation cluster of the gating residue; and changes in the occupancy rate of pocket water molecules or the accessibility of local solvents.

[0028] Preferably, the enantioselective scoring is a weighted combination of at least the following or a monotone equivalent transformation thereof: inlet flux ratio, productive flux ratio, nonproductive penalty ratio, gating score difference, and main exit channel usage probability ratio.

[0029] Preferably, the site priority score is a weighted combination of at least the following indicators or their monotone equivalent transformations: the frequency with which the site participates in channel bottlenecks; the contact asymmetry between the site and different enantiomers in NAC or NPV states; the difference in the associated motion of the site with catalytic residues or gated residues; the sensitivity of the ESS after applying local geometric perturbation or in situ virtual mutation to the site; and the structural role factor of the site, which indicates that the site belongs to a channel ROI, a lateral cavity ROI, near an oxygen anion hole, a gate frame site, and / or a gated site.

[0030] Preferably, the method includes the following steps:

[0031] S1. Obtain the three-dimensional structure of the target enzyme, one or more pairs of enantiomers of the substrate to be compared, and the reference residue set of the target enzyme active site; perform hydrogen supplementation, protonation state determination, missing atom completion, small molecule parameterization, solvation, and ion balance processing on the target enzyme structure;

[0032] S2. Based on the three-dimensional structure of the target enzyme, candidate channels leading from the active site to the protein surface are identified, and channel data for each candidate channel is generated; the channel data includes the centerline, bottleneck radius, channel length, channel curvature, exit region, and channel ROI residue set;

[0033] S3. Based on the channel data, perform an inlet sampling simulation for each enantiomer to obtain inlet sampling simulation data; the inlet sampling simulation data includes at least the surface trapping pit, the main inlet channel, and the intermediate metastable region.

[0034] S4. Based on channel data, perform binding state resampling simulation for each enantiomer to obtain binding state resampling simulation data; the binding state resampling simulation data includes at least the trajectory data of the substrate in the active pocket pre-organization, near-attack conformation state, and non-productive side cavity state; the trajectory data includes at least the occupancy probability, transfer events, residence time, and key interactions;

[0035] S5. Based on the channel data, perform a segmented random external force escape simulation for each enantiomer to obtain segmented random external force escape simulation data; the segmented random external force escape simulation data includes at least the escape path of the substrate or product from the active cavity to the protein surface, the escape time, and the probability of use of the exit channel;

[0036] S6. Discretize the obtained inlet sampling simulation data, combined state repeated sampling simulation data, and segmented random external force escape simulation data into a state discretization and construct a conformational state diagram; the conformational state diagram includes at least one or more of the following: bulk state, surface capture state, channel inlet state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and outlet state.

[0037] S7. Calculate the conformational flux of each enantiomer based on the conformational state diagram; the conformational flux includes at least one or more of the following: inlet flux, near-attack conformational flux, non-productive residence penalty, outlet flux, and gating score;

[0038] S8. Calculate enantioselectivity scores and site priority scores based on the differences in conformational flux between enantiomers;

[0039] S9. Based on site priority scoring, output at least one of the following: candidate mutation site, recommended mutation direction, enantioselectivity prediction, site priority, and site functional type.

[0040] This invention also provides an enzyme enantioselective mutation site screening system based on conformational flux, comprising:

[0041] The structure input module is used to obtain the three-dimensional structure of the target enzyme, the enantiomer of the substrate to be compared, and the reference residues of the active site of the target enzyme.

[0042] The channel recognition module, based on the three-dimensional structure of the target enzyme, identifies candidate channels leading from the active site to the protein surface and generates channel data for each candidate channel.

[0043] The ingress sampling module is used to perform ingress sampling simulation for each enantiomer based on channel data to obtain ingress sampling simulation data;

[0044] The associative state resampling module is used to perform associative state resampling simulation for each enantiomer based on channel data, and obtain associative state resampling simulation data;

[0045] The segmented random external force escape module is used to perform segmented random external force escape simulation for each mapper based on channel data, and obtain segmented random external force escape simulation data.

[0046] The state diagram construction module is used to discretize the obtained inlet sampling simulation data, combined state repeated sampling simulation data, and segmented random external force escape simulation data and construct a conformational state diagram; the conformational state diagram includes at least one or more of the following: bulk phase state, surface capture state, channel inlet state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and outlet state.

[0047] The conformation flux calculation module is used to calculate the conformation flux of each enantiomer based on the conformation state diagram; the conformation flux includes at least one or more of the following: inlet flux, near-attack conformation flux, non-productive residence penalty, outlet flux, and gating score;

[0048] The site priority scoring module is used to calculate enantioselectivity scores and site priority scores based on the differences in conformational flux between enantiomers.

[0049] The mutation site output module is used to output at least one of the following based on site priority scoring: candidate mutation sites, recommended mutation directions, and enantioselectivity predictions.

[0050] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of any one of the methods described herein.

[0051] In studies on the enantiomeric resolution of methyl 3-cyclohexene-1-carboxylate (CHCM) catalyzed by the esterase Est13, insufficient screening methods for mutation sites were also identified. Model and simulation results showed that Est13 preferred to catalyze R-configuration substrates; R-CHCM could reside within the Ser201 catalytic cavity and form a productive conformation favorable for nucleophilic attack; S-CHCM, on the other hand, tended to enter the other cavity and form a non-productive stagnation. Further analysis revealed that R-configuration substrates exhibited a gating switching phenomenon regulated by both channels and pockets during catalysis, while S-configuration substrates easily formed a "double-anchoring" effect in the Leu227-related side cavity. These phenomena indicate that the enantioselectivity of Est13 cannot be explained solely by static docking conformations or the distance to a single active site, but requires a systematic analysis from the perspective of substrate conformation flux allocation. Therefore, the enantioselective mutation sites of Est13 esterase were screened using the method described above in this invention. Specifically, Est13 esterase was used as the target enzyme, and R / S-CHCM was used as the substrate. The enantioselective output was used to guide the mutation design of Leu227, Ala230, Phe233, Ile254, Val257, Thr258, Thr137, and His139 sites. The results showed that Leu227 was identified as a shunt valve site, and Val257, Ala230, Phe233, Ile254, and Thr258 were identified as gate or bottleneck regulation sites. The mutation direction was recommended according to the different contributions of the site in productive or non-productive branches, namely, disrupting non-productive clamping, adjusting bottleneck volume, changing local polarity, or changing gating conformation preference.

[0052] The present invention also provides the application of shunt valve sites and / or channel gate sites in the construction of Est13 esterase mutants, wherein the shunt valve sites include at least Leu227; and the channel gate sites include at least Ala230, Phe233, Ile254, Val257, and Thr258.

[0053] This invention also provides an Est13 esterase mutant, obtained by any one of the following mutations into the amino acid sequence shown in SEQ ID NO.1:

[0054] (1) Y228C / D253A;

[0055] (2) Y228C / D253A / A230Y;

[0056] (3) Y228C / D253A / I254V;

[0057] (4) Y228C / D253A / F233K;

[0058] (5) Y228C / D253A / V257A;

[0059] (6) Y228C / D253A / V257M;

[0060] (7) Y228C / D253A / V257S;

[0061] (8) Y228C / D253A / T258E;

[0062] (9) Y228C / D253A / L227M;

[0063] (10) Y228C / D253A / L227F;

[0064] (11) Y228C / D253A / L227W;

[0065] (12) Y228C / D253A / L227R;

[0066] (13) Y228C / D253A / L227N;

[0067] (14) Y228C / D253A / L227H;

[0068] (15) Y228C / D253A / L227K;

[0069] (16) Y228C / D253A / L227S.

[0070] The beneficial effects of this invention are:

[0071] (1) This invention integrates channel geometry identification, inlet sampling, bound state molecular dynamics, escape path statistics, near-attack conformation identification, non-productive side cavity analysis and gating switching detection into the conformational state diagram, forming a complete enzyme enantioselectivity dynamic analysis process.

[0072] (2) This invention quantitatively distinguishes between productive and non-productive pathways by using conformation flux index, which can identify the key problem that substrates can enter enzymes but cannot form catalytic conformations, and avoids the traditional docking method from misjudging binding as catalytic.

[0073] (3) This invention improves the reliability of ligand escape path statistics by using a dual-criteria escape judgment and channel classification algorithm. The dual criteria simultaneously consider the centroid distance between the ligand and the pocket reference group and the number of heavy atom contacts, and use robust statistical values ​​of the trajectory end window to reduce misjudgments caused by instantaneous conformational fluctuations.

[0074] (4) This invention incorporates dynamic cross-correlation, channel bottleneck radius, dihedral state of gated residues and changes in local solvent accessibility into the analysis by gating score, transforming the originally qualitative gating switching phenomenon into a calculable, comparable and applicable quantitative indicator for site ranking.

[0075] (5) The present invention can output the functional type of candidate mutation sites, including main channel gate site, channel bottleneck site, side cavity anchoring site, diversion valve site and remote gate coupling site, and give recommended mutation direction, thereby significantly reducing the size of experimental mutation library.

[0076] (6) This invention has been validated in the asymmetric resolution system of CHCM by Est13 esterase. Experimental results show that mutations at channel gate sites A230, F233, I254, V257, and T258 can improve R conformation selectivity to varying degrees; the hydrophobic large-volume mutation of Leu227 maintains R selectivity and increases the eep value, while the polar or charged mutation of Leu227 can lead to a reversal of conformation selectivity towards the S conformation. These results are consistent with the mechanism for regulating the productive flux of the main channel and the non-productive flux of the lateral lumen proposed in this invention. Attached Figure Description

[0077] Figure 1 The overall flowchart of the method of the present invention includes structural input, channel identification, entry sampling, binding state repeated sampling, segmented random external force escape, state diagram construction, conformation flux calculation, site priority scoring and mutation site output.

[0078] Figure 2 This is a schematic diagram of the conformational state diagram. The nodes include the B bulk state, C surface capture state, T channel inlet state, P active pocket pre-organized state, N near-attack conformational state, V non-productive pocket state, and X outlet state. The edges represent the transition events between states and the corresponding conformational flux.

[0079] Figure 3 This is a schematic diagram of the R-CHCM productive pathway and S-CHCM non-productive pathway in the Est13 esterase examples, showing the T3 main channel, Pocket-R, Pocket-S, Ser201 catalytic chamber, Leu227 shunt valve site, and V257 gate sites. Detailed Implementation

[0080] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0081] refer to Figure 1 This invention provides a method for screening enzyme enantioselective mutation sites based on conformational flux, mainly including structural input, channel identification, entry sampling, repeated sampling of binding states, segmented random external force escape, state diagram construction, conformational flux calculation, site priority scoring, and mutation site output. The specific steps are detailed below.

[0082] S1. Data Input and Structure Preparation

[0083] The method described in this invention first obtains the three-dimensional structure of the target enzyme. The target enzyme structure can be derived from crystal structures, cryo-electron microscopy structures, NMR structures, homology modeling structures, AlphaFold predicted structures, or other structure prediction models.

[0084] For the target enzyme structure, the following pretreatment is preferred: removal of irrelevant ligands, water molecules, buffer salt ions, or crystallization additives; if the water molecules participate in catalysis or stabilize oxygen anion pores, key structural water can be retained; missing residues, missing side chains, and hydrogen atoms are completed; and the protonation state of amino acids is determined according to experimental pH or target reaction conditions.

[0085] Three-dimensional structures of the substrate in both the R and S configurations are established and parameterized using appropriate small molecule force fields. Preferably, the substrate parameters can be generated using GAFF, GAFF2, CGenFF, OpenFF, or other small molecule force fields suitable for molecular dynamics simulations; the charge can be determined using RESP, AM1-BCC, or equivalent methods.

[0086] Establish the initial conformation of the enzyme-substrate complex. The initial conformation can be derived from molecular docking, manual placement, existing molecular dynamics equilibrium conformation, or experimental structure. For esterase systems, it is preferable to orient the substrate carbonyl carbon toward the catalytic serine Oγ atom and the carbonyl oxygen toward the oxygen anion pore region.

[0087] The reference set of active sites for a target enzyme may include catalytic residues, oxyanion pore residues, active cavity shell residues, or residues near the initial substrate binding site. For example, for Est13 esterase, Ser201, Asp200, His325, Gly129, Gly130, and adjacent active cavity residues can be defined as the reference set of active sites.

[0088] S2. Candidate Channel Identification and ROI Definition

[0089] After structural preprocessing, candidate channels connecting the active cavity and the protein surface are identified starting from the active site reference point. Candidate channel identification can employ CAVER, MOLE, HOLLOW, HOLE, or other channel identification algorithms. Preferably, channel search is performed on both static structures and molecular dynamics trajectory frame structures to obtain a dynamic channel set.

[0090] For each candidate channel, at least the following information should be recorded: channel ID, coordinates of the channel centerline node, local radius of each centerline node, channel length, channel bottleneck radius and bottleneck location, channel curvature, channel exit region, and set of channel ROI residues.

[0091] A channel ROI set refers to a set of residues located at the channel inlet, bottleneck, middle, or outlet, or in frequent contact with the substrate migration path. ROI residues can be defined by factors such as being less than a preset threshold distance from the channel centerline, constituting a channel bottleneck, having frequent contact with the substrate, being located at a pocket shunt location, or being at the inlet of a side cavity.

[0092] For Est13 esterase, candidate channels such as T1, T2, T3, and T4 can be preferentially recognized. Among them, the T3 channel can be defined as the main channel. T3-related ROI residues may include, but are not limited to, Leu227, Ala230, Phe233, Ile254, Val257, Thr258, etc.; non-productive lateral ROI residues may include, but are not limited to, Phe125, Val142, Leu146, Leu225, Leu227, Ser345, etc.

[0093] S3. Inlet Sampling Simulation

[0094] To identify the ability of substrates to enter the enzyme surface capture region and channel entrance from the solution phase, this invention employs entrance sampling simulation. The entrance sampling simulation preferably uses ligand-copy flooding molecular dynamics simulation.

[0095] In this step, multiple substrate copies are placed in a solvent environment surrounding the protein. The number of substrate copies is preferably 10 to 100, more preferably 20 to 60. The initial substrate placement region is preferably a spatial shell 0.8 to 3.0 nm from the protein surface, more preferably 1.2 to 2.0 nm. The initial minimum distance between any two substrate copies of heavy atoms is preferably not less than 0.6 to 1.2 nm, more preferably not less than 0.8 nm, to avoid initial overlap and unreasonable interactions.

[0096] Establish separate systems for each enantiomer. For R-configuration and S-configuration substrates, use the same protein structure, solvation method, ion concentration, and simulation parameters to ensure fair comparison.

[0097] In the inlet sampling simulation, the time of the first sustained contact between the substrate and the protein surface, the time of the first entry into the surface capture pit, the time of the first contact with the ROI residue of the channel inlet, the time of the first entry into the region near the center line of the candidate channel, the time of the first entry into the pre-organized region of the active pocket, and the residence time of the substrate in different inlet regions and channels were recorded.

[0098] Surface capture state C can be defined as: continuous contact between the substrate and the ROI region on the protein surface for more than [time period missing]. And the number of contact residues is greater than .in, Preferably, it is 50 to 500 ps, ​​more preferably 100 ps; Preferably 2 to 8, more preferably 3.

[0099] Channel entrance status T j This can be defined as: the shortest distance from the substrate center to the centerline of candidate channel j is less than d. T Furthermore, the substrate comes into contact with the ROI residues of channel j. T Preferably, it is 3.0 to 8.0 Å, more preferably 5.0 Å.

[0100] S4. Repetitive molecular dynamics sampling of bound states

[0101] To obtain conformational pre-organization, productive near-attack conformation, and non-productive side-cavity stagnation behavior of the substrate in the active pocket, this invention performs binding-state repetitive conventional molecular dynamics simulations for each enantiomer.

[0102] Each enantiomer preferably executes at least three independent trajectories, more preferably five to ten independent trajectories. The length of each trajectory can be 100 to 500 ns, preferably 300 ns. For systems with slow substrate diffusion and side cavity migration events, the length can be extended to the microsecond scale or supplemented by enhanced sampling methods.

[0103] Each trajectory should undergo energy minimization, gradual heating, and equilibration before formal production simulation. The preferred process includes: energy minimization to eliminate initial unreasonable contact; gradual heating of the NVT ensemble to raise the system from a low temperature to a target temperature, which can be 298–310 K; NPT ensemble density equilibration to stabilize the system pressure at 1 atm or 1 bar; gradual release of protein backbone and substrate constraints; and unconstrained or weakly constrained production simulation.

[0104] All trajectories should preserve substrate coordinates, coordinates of key protein residues, conformation of channel bottleneck residues, dihedral angles of gating residues, and necessary energy terms for subsequent state labeling and conformational flux calculation.

[0105] S4.1. Near-Attack Conformation NAC Recognition

[0106] Near-attack conformation (NAC) refers to a substrate conformation that is geometrically close to a state in which a catalytic reaction can occur. In esterase or lipase systems, NAC should typically meet the geometric requirements for catalyzing the nucleophilic attack of serine on the carbonyl carbon of the ester, and the carbonyl oxygen should be stabilized by an oxygen anion pore.

[0107] For serine hydrolase systems, NAC is preferably defined by at least two of the following conditions: the distance between the catalytic SerOγ atom and the carbonyl carbon atom of the substrate ester. Not greater than 3.5 Å, more preferably not greater than 3.2 Å; nucleophilic attack angle Located in the range of 90° to 120°, more preferably in the range of 95° to 115°; at least one hydrogen bond is formed between the substrate carbonyl oxygen and the oxygen anion pore hydrogen-donating atom, more preferably two hydrogen bonds are formed; the hydrogen bond donor-acceptor distance is not greater than 3.5 Å; the orientation of the substrate leaving group or hydrophobic framework is compatible with subsequent catalytic steps.

[0108] NAC occupancy rate can be defined as the proportion of trajectory frames that satisfy the NAC criterion to the total number of trajectory frames. NAC formation frequency can be defined as the number of times a pre-organized state P transitions to NAC state N per unit time.

[0109] Productive near-attack conformational flux It can be defined as:

[0110]

[0111] Where e represents the enumeration type, e∈{R,S}; This represents the number of times enantiomer e transitions from the active pocket pre-organized state P to the near-attack conformation state N; This represents the total observation time of enantiomer e.

[0112] S4.2 Non-Productive Pocket NPV Recognition

[0113] Nonproductive pocket state (NPV) refers to a state in which the substrate enters the enzyme interior or side cavity, forms a stable contact with the protein, but cannot meet the conformational requirements for catalytic proximity attack.

[0114] NPV status can be defined by the following conditions: the substrate core is located within a predefined nonproductive cavity ROI; and the number of substrate contacts with cavity shell residues is greater than a preset threshold. The substrate remains continuously in this cavity for a period of time greater than [a certain duration]. The substrate does not satisfy at least two of the NAC criteria.

[0115] Preferably 2–50 ns, more preferably 5–20 ns; The preferred values ​​are 3 to 8.

[0116] Non-productive residency penalties It can be defined as:

[0117]

[0118] in, This represents the rate at which enantiomer e transitions from the pre-organized state P to the unproductive pocket state V; This represents the probability of the enantiomer e occupying the unproductive state V; The weight is non-negative, preferably 0 to 1, and more preferably 0.5.

[0119] For the Est13 esterase system, S-CHCM can form a non-productive luminal retention near residues such as Leu227 and Ser345. This is characterized by the substrate being "double-anchored" by the luminal residues, preventing it from forming a stable nucleophilic attack geometry with Ser201. This state can be defined as the NPV state of the S-configuration substrate.

[0120] S5. Piecewise random external force escape simulation

[0121] To address the challenge of conventional molecular dynamics in repeatedly observing ligand escape paths within a finite timeframe, this invention employs a piecewise stochastic external force escape simulation.

[0122] The segmented random external force escape simulation includes the following steps: starting from the equilibrium conformation of the enzyme-substrate complex; applying an external force in a random direction to the substrate or product core in each simulation segment; regenerating a random unit vector as the direction of the external force for the next segment after each simulation segment ends; continuing the simulation with the last frame of the previous segment as the starting point of the next segment; terminating the repeated trajectory when the escape criterion is met or the maximum number of segments is reached; and statistically analyzing the escape time, escape direction, and escape path for multiple repeated trajectories.

[0123] The duration of each simulation segment is preferably 50 to 500 ps, ​​more preferably 200 ps; the time step is preferably 1 to 4 fs, more preferably 2 fs; the maximum number of segments is preferably 20 to 200, more preferably 60; and the number of repetitions is preferably 10 to 200, more preferably 20 to 50.

[0124] The preferred escape criterion is a dual criterion of "distance-number of contacts," using robust statistics from the tail window.

[0125]

[0126] in, This represents the distance between the substrate centroid and the reference centroid of the active pocket; This indicates the number of heavy atom contacts between the substrate and the reference set in the active pocket; The distance threshold is preferably 1.0–5.0 nm, more preferably 2.0 nm; The contact number threshold is preferably 0 to 10, more preferably 3; MedianTail means taking the median of 3 to 11 sampling points at the end of the trajectory, more preferably taking 3 or 5 sampling points at the end.

[0127] Using logical "AND" instead of logical "OR" as the escape criterion can avoid misjudging a conformation where the substrate is temporarily away from the pocket but still in contact with the pocket, or where the contact is temporarily reduced but has not actually left the pocket, as an escape, thereby improving the authenticity and robustness of escape event judgment.

[0128] S5.1. Escape Channel Classification Algorithm

[0129] For each successful escape trajectory, this invention uses a combined criterion of "nearest distance to the centerline + ROI crossing" to classify the channels.

[0130] For the ligand position at the end of the escape trajectory and the center line of candidate channel j The shortest distance from the ligand to the centerline of the channel is defined as:

[0131]

[0132] Meanwhile, the proportion of ligands that come into contact with the ROI residue set of channel j during the escape process is statistically analyzed. And the contact ratio between the end of the trajectory and the set of residues in the ROI. .

[0133] Channel classification score can be defined as:

[0134]

[0135] Wherein, β1, β2, and β3 are non-negative weights, preferably 0.4–0.8, 0.1–0.4, and 0.1–0.4, respectively; It is a small constant, preferably 10. -6 ~10 -4 .

[0136] Score j The channel with the highest score exceeding the threshold is defined as the channel to which the escape trajectory belongs. If all channel scores are below the threshold, the trajectory can be marked as unclassified or an unknown exit.

[0137] For each enantiomer, count the number of times n escapes through each channel.exit,j And calculate the export flux:

[0138]

[0139] S6. The obtained inlet sampling simulation data, combined state repeated sampling simulation data, and segmented random external force escape simulation data are discretized into states and a conformational state diagram is constructed; the conformational state diagram includes the bulk state, surface capture state, channel inlet state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and outlet state.

[0140] S7. Conformation flux calculation

[0141] S7.1 Gating Switching Rating

[0142] Gating switching refers to the phenomenon where conformational changes of residues around a channel or pocket alter the substrate migration path, pocket openness, or the probability of catalytic conformation formation.

[0143] This invention uses a multi-index combined approach to define the gating score G. gate The gating score may include at least two types of information: dynamic cross-correlation changes between the catalytic residue region and the channel ROI region; time-series changes in the bottleneck radius of the candidate channel; the χ angle of the side chain of the gating residue, the dihedral angle of the main chain, or the state transition of the side chain conformation cluster; and changes in local solvent accessibility or pocket water molecule occupancy.

[0144] In a preferred embodiment, the gating score can be defined as:

[0145]

[0146] in, This represents the average dynamic cross-correlation between the catalytic region and the ROI of the T3 channel; This represents the average dynamic cross-correlation between the catalytic region and the ROI of the T2 channel; and These represent the average effective radii of the bottlenecks in channels T3 and T2, respectively. This represents the change in the probability that the gated residue is in an open rotor state; and The weights are non-negative, preferably 0.1 to 2.0; Norm represents normalization.

[0147] For the Est13 system, if the correlation between the T3 channel and the Ser201 catalytic region in the R-CHCM trajectory is enhanced, and the T3 bottleneck radius increases or the gated residues are converted to open conformations, it can be considered that the R-configuration substrate induces or utilizes a gated switching that is conducive to productive pathways.

[0148] S8.1 Enantioch-selective scoring

[0149] This invention defines the substrate state transition rate as:

[0150]

[0151] in, Let e ​​be the number of times the enantiomer e transitions from state a to state b; The total observation time is for enantiomer e.

[0152] The inlet flux can be defined as:

[0153]

[0154] in, Indicates the main productive entry channel; The weight is 0 to 1, and more preferably 1.

[0155] The enantioselectivity score (ESS) can be defined as:

[0156]

[0157] Here, w1 to w5 are non-negative weights. (Default values ​​are acceptable.) Alternatively, it can be recalibrated based on known experimental data; It is a small constant, preferably 10. -6 ~10 -4 .

[0158] When ESS > 0, the prediction system is biased towards the R configuration; when ESS < 0, the prediction system is biased towards the S configuration; the larger the absolute value of ESS, the stronger the prediction enantioselectivity.

[0159] In some implementations, the ESS can be further mapped to a signed eep prediction value:

[0160]

[0161] in, and These are calibration parameters. This formula can be used for sorting and direction prediction, but is not a necessary limitation for implementing this invention.

[0162] S8.2 Mutation Site Priority Scoring

[0163] For any candidate residue i, this invention defines a site priority score M. i :

[0164]

[0165] Among them, B iLet A be the bottleneck participation frequency, representing the normalized frequency of residue i appearing in the bottleneck shell of the candidate channel; i Contact asymmetry represents the difference in contact between residue i and its R / S enantiomer in the NAC or NPV states; C i U is the gating coupling index, representing the difference in the correlation motion between residue i and catalytic or gating residues; i Local mutation sensitivity indicates the degree of change in ESS after applying local geometric perturbation, volume change, or in situ virtual mutation to residue i; R i λ1 is the structural role factor, indicating whether residue i belongs to the channel ROI, the side cavity ROI, the vicinity of the oxygen anion pore, the gate frame site, the shunt valve site, or the gated site; λ1 to λ5 are non-negative weights, which can be taken as (0.25, 0.25, 0.20, 0.20, 0.10) by default.

[0166] S9. Mutation site output

[0167] According to M i The scores can be used to classify candidate sites into main channel gate frame sites, channel bottleneck sites, side cavity anchoring sites, diversion valve sites, and remote gated coupling sites.

[0168] Based on the site type, recommended mutation directions can be further output. If the residue mainly stabilizes the non-productive side cavity stagnation of the non-dominant configuration substrate, it is recommended to weaken the error clamping by changing the polarity, charge, hydrogen bonding ability, or spatial volume. If the residue is located at the bottleneck of the productive main channel of the dominant configuration substrate, it is recommended to optimize the productive flux by adjusting the side chain volume, hydrophobicity, or flexibility. If the residue is a shunt valve site, hydrophobic large-volume mutations can be designed to enhance the original selectivity, or polar / charged mutations can be designed to induce enantioselectivity flipping.

[0169] Table 1. System Modules and Input / Output File Formats

[0170]

[0171] Table 2. Parameters and Recommended Range

[0172]

[0173] Example 1. Conformational flux analysis of asymmetric resolution of CHCM by Est13 esterase

[0174] 1. Est13 structural preparation and active site definition

[0175] This example illustrates how the method of the present invention, based on enzyme structure, substrate enantiomers, channel trajectories, and experimental validation data, completes the screening of enantioselective mutation sites using the asymmetric resolution of 3-cyclohexene-1-carboxylate methyl ester (CHCM) catalyzed by the esterase Est13, to achieve enantioselective mutation site screening. CHCM includes two configurations: R-CHCM and S-CHCM. This example uses the differences in entry points, channel migration, active pocket conformation pre-organization, near-attack conformation formation, non-productive lateral lumen retention, and escape pathways of the R / S conformation substrates within the same enzyme channel network as the objects of conformational flux calculation. The amino acid sequence of the Est13 esterase used in this example is SEQ ID NO.1. A His tag is added to the C-terminus of sequence SEQ ID NO.1.

[0176] When performing structure numbering, mutation site labeling, and catalytic triplet localization, the Est13 esterase sequence corresponding to SEQ ID NO.1 was used as the benchmark. For the C-terminal His tag region, it can be retained or truncated in structural modeling and molecular dynamics simulations depending on the model quality; if the His tag does not participate in the calculation of catalytic pockets, substrate channels, or conformational flux, it is preferable to truncate the flexible tag before structural simulation to avoid irrelevant interference with solvent regions and trajectory analysis.

[0177] A three-dimensional structural model of the Est13 esterase was obtained. This structural model was derived from the prediction results of the AlphaFold3 structure prediction software. Models with high pLDDT and reasonable conformations of the catalytic and channel regions were preferred; if there were obvious missing ring regions or unreasonable side chain orientations, they should be corrected through homology modeling, energy minimization, or short-range confinement molecular dynamics.

[0178] The structural model was preprocessed, specifically including: deleting heterologous small molecules and irrelevant water molecules that do not participate in catalysis; retaining key water molecules that may participate in oxygen anion pore stability, catalytic proton transfer, or pocket structure stability; completing missing hydrogen atoms and missing side chains; determining the protonation state of residues such as His, Asp, Glu, Lys, and Arg according to the pH of the target reaction system; examining the spatial relationship between the catalytic triplet Ser201, Asp200, and His325; and examining the positions where oxygen anion pore residues such as Gly129 and Gly130 may form hydrogen bonds with the substrate carbonyl oxygen.

[0179] Ser201, Asp200, His325, Gly129, Gly130, and pocket residues within 8 Å of the Ser201 Oγ atom are defined as the active site reference residue set POCKET; residues that may participate in channel entry, channel bottleneck, side cavity shunting, and gating switching are defined as the candidate ROI residue set. This definition is used for subsequent calculations of substrate-pocket centroid distance, heavy atom contact number, NAC criterion, NPV criterion, and escape criterion.

[0180] 2. Substrate configuration construction and initial model of enzyme-substrate complex

[0181] Construct the three-dimensional structures of R-CHCM and S-CHCM separately. Initial conformations can be generated using SMILES with RDKit, OpenBabel, SchrodingerLigPrep, or equivalent software, and geometry optimization can be performed for each conformation. Small molecule force fields can be GAFF, GAFF2, CGenFF, OpenFF, or equivalent force fields; charges can be determined using AM1-BCC, RESP, or RESP2 methods. The same small molecule parameterization procedure should be used for both R and S conformations to ensure that the comparison results originate from conformational differences rather than parameter differences.

[0182] Initial complexes of Est13-R-CHCM and Est13-S-CHCM were established using molecular docking or manual placement, respectively. The docking cassette center was preferably located near the potential reaction centers of Ser201 Oγ, His325 Nε2, and the substrate carbonyl group, with the cassette covering the Ser201 catalytic cavity, the end of the T3 channel, and the Leu227 side cavity region. The docking conformation screening was not based solely on the lowest binding energy; instead, conformations that could enter the active cavity, had the carbonyl group facing the oxygen anion pore, and did not undergo significant atomic collisions were preferentially selected as the initial molecular dynamic conformations for the bound state.

[0183] For esterase systems, if the initial conformation is used for near-attack conformation analysis, the distance between Ser201 Oγ and the carbonyl carbon of the substrate ester, the angle of attack of Ser201 Oγ-C=O, and the probability of hydrogen bonding between the substrate carbonyl oxygen and the NH in the Gly129 / Gly130 main chain should be checked. If the initial conformation clearly does not meet the reasonable catalytic orientation, short-range confinement equilibration or re-docking can be performed first to avoid distortion of subsequent state labeling due to an unreasonable initial model.

[0184] 3. Simulation of solvation, heating, equilibration, and repetitive production.

[0185] Each complex is placed in an explicit water tank, with the tank boundary preferably at least 1.0 nm from the outermost atoms of the protein. Na+, Cl-, or equivalent ions are added to neutralize the system charge, and the ionic strength is set to 0.10-0.15 mol / L as needed. The protein force field can be Amber ff14SB, Amber ff19SB, CHARMM36m, OPLS-AA / M, or an equivalent force field; the water model can be TIP3P, SPC / E, or a water model matched to the selected protein force field.

[0186] Each system was processed according to standard molecular dynamics procedures: first, energy minimization was performed until the maximum force was below a preset threshold or a preset number of steps were reached; then, the system was gradually heated to the target temperature in the NVT ensemble, preferably 300 K, which could be completed within 100-500 ps; during the heating phase, positional constraints were preferably imposed on the heavy atoms of the protein backbone and the substrate to avoid sudden shifts in the initial conformation. Afterwards, density equilibration was performed in the NPT ensemble, with a pressure preferably of 1 bar or 1 atm and an equilibration time preferably of 0.5-2 ns, and the positional constraints could be reduced in stages. Once the system temperature, pressure, density, and total energy stabilized, unconstrained or weakly constrained production simulations were initiated.

[0187] To avoid the randomness of single trajectories, this embodiment performs 10 independent, repeated production simulations for both the R-CHCM and S-CHCM systems. Each production trajectory is 300 ns long, and the cumulative sampling time for each substrate configuration is 3 μs. The 10 repeated trajectories are started with different random velocity seeds; if starting from the same initial complex, velocities can be regenerated for each trajectory after equilibration. The production simulation preferably saves coordinate frames every 10-100 ps for subsequent conformational state annotation, contact analysis, DCCM calculation, channel bottleneck radius statistics, and NAC / NPV event identification.

[0188] 4. Candidate Channel Identification and T3 Master Channel Definition

[0189] Using the Ser201 Oγ atom or the centroid of the Ser201, Asp200, and His325 catalytic triplet as the starting point for channel search, CAVER, MOLE, or equivalent channel identification programs were used to search for channels in the static structure and production trajectory frame structure of Est13. During the search, the centerline node, local radius, average radius, bottleneck radius, channel length, curvature, exit position, and channel shell residues of each candidate channel were recorded.

[0190] Based on the channel centerline and shell residues, candidate channels in Est13 connecting the active cavity to the protein surface are labeled as T1, T2, T3, and T4. The channel ROI residues are defined as residues that occur frequently within 5 Å of the channel centerline or less than 4 Å from substrate heavy atoms in any trajectory frame. For the T3 channel, ROI residues include, but are not limited to, Leu227, Ala230, Phe233, Ile254, Val257, and Thr258; these residues are distributed in the later part of the T3 channel, the bottleneck region, and the active cavity entrance region, and are therefore defined as T3 channel gate frames or bottleneck candidate sites.

[0191] If a channel exists in a static structure but its bottleneck radius remains too small in most production trajectories, and neither inlet sampling nor escape sampling shows that the substrate passes through the channel, then this channel can be labeled as a geometric candidate channel rather than a functional main channel. If a channel simultaneously satisfies optimal geometric parameters, high-frequency inlet sampling, substrate entry into the active chamber along the channel in the bound state trajectory, and can serve as a main outlet or reversible channel in escape sampling, then it can be labeled as a main functional channel. In this embodiment, considering channel geometry, inlet sampling, and trajectory behavior, T3 is defined as the main productive channel for Est13-catalyzed CHCM separation.

[0192] 5. Flooding Inlet Sampling and Surface Capture Event Recognition

[0193] To identify the differences between R / S substrates entering the protein surface capture region and channel entrance from the solvent phase, entry sampling systems were established for Est13 with multiple R-CHCM copies and multiple S-CHCM copies, respectively. The preferred substrate copy number in each system was 20-60, with the substrate initially distributed within a solvent shell 1.2-2.0 nm from the protein surface, and the initial heavy atom distance between any two substrate copies not less than 0.8 nm to avoid initial overlap.

[0194] The inlet sampling system also performs energy minimization, NVT heating, NPT equilibration, and production simulation. The inlet sampling production trajectory length is preferably 100-300 ns. Each substrate configuration is established and simulated separately to avoid statistical confounding of inlet events caused by competition between R / S substrates.

[0195] In the ingress sampling trajectory, events such as the first sustained contact between the substrate and the protein surface, the first entry into the surface capture pit, the first contact with T3 ingress ROI residues, the first entry into the 5 Å proximity region of the T3 centerline, and the first entry into the pre-organized region of the Ser201 active pocket were recorded. Surface capture events were defined as continuous contact between the substrate and the same surface ROI region for more than 100 ps and with at least 3 contacting residues; T3 ingress events were defined as substrate centers less than 5 Å from the T3 centerline and having heavy atom contact with T3 ingress ROI residues.

[0196] Ingress flux was calculated using the number of ingress events, first arrival time, and residence time distribution. If R-CHCM enters the T3 ingress more frequently and reaches the pre-organized region of the active pocket faster than S-CHCM, it indicates that the R-configuration substrate has a higher productive ingress flux. If S-CHCM can enter the T3 ingress but then deflects to the side cavity, it is further determined by NPV flux.

[0197] 6. Combining state trajectory labeling with NAC / NPV determination

[0198] The states of 10 R-CHCM binding state trajectories and 10 S-CHCM binding state trajectories were labeled frame by frame. The states included at least: B bulk state, C surface capture state, T channel inlet state, P active pocket pre-organized state, N near-attack conformation state, V non-productive pocket state, and X exit state. This embodiment focuses on analyzing the P, N, and V internal states.

[0199] The P-state is defined as the substrate located in the Ser201 catalytic cavity or its vestibular region, with the substrate centroid less than a preset threshold from the POCKET reference set centroid, and in contact with the shell residues of the active cavity, but not yet fully satisfying the NAC criterion. The N-state is defined as the substrate satisfying the near-attack conformation criterion: the distance between the Ser201 Oγ and the substrate ester carbonyl carbon is no greater than 3.5 Å, the nucleophilic attack angle is within the range of 90°–120°, and the substrate carbonyl oxygen forms at least one hydrogen bond with the hydrogen atom donated by the Gly129 / Gly130 oxygen anion pore. The robustness of these conditions can be determined by satisfying at least two of them, depending on the simulation system and substrate size.

[0200] The V state is defined as the substrate residing in the Pocket-S or other non-productive cavities, forming stable contact with cavity residues, residing continuously for at least 5 ns, and not satisfying at least two of the NAC criteria. For the Est13 system, the Pocket-S can be defined by Leu227, Ser345, and neighboring residues such as Phe125, Val142, Leu146, and Leu225. When S-CHCM forms a "double anchor" in this region, although it can be stably bound by the protein, its carbonyl orientation and Ser201 attack geometry are unfavorable, therefore it is labeled as an NPV non-productive state.

[0201] In the R-CHCM representative trajectory, a productive near-attack conformation located within the Ser201 catalytic cavity can be captured at approximately 257.834 ns. In this conformation, the substrate carbonyl group interacts stably with the oxygen anion pore, and the distance between the Ser201 Oγ and the substrate carbonyl carbon is within the nucleophilic attack allowable range. In contrast, in the S-CHCM 300 ns representative conformation, the substrate is located in a non-productive cavity on the other side, held by the Leu227 / Ser345 neighborhood, and cannot form an effective near-attack conformation. This comparison provides structural evidence for the high Φ_NAC of the R configuration and the high Ω_NPV of the S configuration in the conformational state diagram.

[0202] 7. Segmented random external force escape simulation and channel classification

[0203] To address the insufficient sampling of channel migration and ligand escape events in conventional molecular dynamics, segmented stochastic external force escape simulations are performed starting from the fully equilibrated conformation of the enzyme-substrate or enzyme-product complex. In each repetition, a constant external force in a random direction is applied to the center of mass of the substrate or product. After each segment of the simulation, a new random unit vector is generated, and the last frame of the previous segment is used as the starting point for the next segment.

[0204] In this embodiment, the segment duration is preferably 200 ps, ​​the time step is preferably 2 fs, the maximum number of segments is preferably 60, and each system preferably executes 20 independent repeating trajectories. The escape criterion uses a distance-contact number dual criterion: the median distance in the final window between the centroid of the substrate or product and the centroid of the POCKET reference group is greater than 2.0 nm, and the median number of heavy atom contacts between the substrate or product and the POCKET reference group in the final window is less than 3. The final window preferably consists of the last 3 or 5 sampling points. Successful escape is determined only when both conditions are met simultaneously.

[0205] In practice, the LIG and POCKET groups are first defined using an index file. The LIG group contains all heavy atoms of the CHCM substrate or product, while the POCKET group contains Ser201, Asp200, His325, Gly129, Gly130, and active cavity residues within 8 Å of the Ser201 Oγ atom. Each system uses a representative conformation after energy minimization, NVT heating, NPT equilibration, and a 300 ns production simulation as the starting point for RAMD, avoiding direct application of external forces to the unequilibrated conformation.

[0206] At the start of each repetition, the starting structure is copied as the 0th segment structure; before the start of each segment, a three-dimensional random unit vector is generated by random numbers, and the current segment's MDP file is automatically generated based on the MDP template. Then, grompp and mdrun are executed sequentially; if a previous segment checkpoint file exists, it is used as the input for the current segment to inherit velocity information and dynamic continuity.

[0207] After each segment, distance analysis was used to calculate the time series of distances between COM(LIG) and COM(POCKET), and minimum distance or contact analysis was used to calculate the time series of heavy atom contact numbers between LIG and POCKET. For the last three consecutive sampling points at the end of each segment, the median distance and the median contact number were calculated respectively; only when the median distance was greater than 2.0 nm and the median contact number was less than 3, the repeat was marked as a successful escape.

[0208] For repeats that successfully escape, all segmented trajectories are concatenated into ramd_full.xtc; for repeats that do not escape within 60 segments, they are recorded as not having escaped. The repeat, total_time_ps, segments, last_dist_nm, last_contacts, and exit_flag of all repeats are written to summary_ramd.csv for subsequent calculation of the proportion of escape channels, average escape time, and exit flux.

[0209] Successful escape trajectories are further categorized into channels. First, the minimum distance from the ligand centroid to the centerline of each channel (T1, T2, T3, and T4) at the end of the escape phase is calculated. Second, the contact frequency, contact duration, and crossing sequence between the ligand and the ROI residues of each channel during escape are statistically analyzed. Finally, the escape trajectory is categorized into the corresponding channel using a joint score of the closest centerline distance and the ROI crossing event. If a trajectory does not meet any channel threshold, it is marked as an unknown exit and is not used for main channel flux calculation or separate classification statistics.

[0210] 8. Conformation flux calculation and site priority output

[0211] For each enantiomer, the ingress sampling, binding state repeating trajectory, and escape trajectory are jointly mapped to the conformational state diagram. The number of state transition events such as B→T3, C→T3, P→N, P→V, and P→X are calculated and normalized to the state transition rate over the total observation time. The ingress flux is jointly determined by events entering the T3 main channel and entering the pre-organized region of the active pocket; the near-attack conformational flux is determined by the rate of transition from the P state to the N state; the non-productive residence penalty term is jointly determined by the rate of transition from the P state to the V state and the probability of V state occupancy; the exit flux is determined by the proportion of escape trajectories released via T3 or other channels.

[0212] Enantioselectivity score (ESS) was calculated based on the conformational flux difference between R / S configurations. When ESS is greater than 0, the predicted system is biased towards the R configuration; when ESS is less than 0, the predicted system is biased towards the S configuration. A site priority score (M_i) was calculated for each residue, which integrates residue bottleneck participation frequency, R / S contact asymmetry, gating coupling index, local virtual mutation sensitivity, and structural role factor. Residues with higher scores were preferentially recommended for site-directed or combinatorial mutagenesis.

[0213] In the Est13 system, A230, F233, I254, V257, and T258 were identified as gate or bottleneck sites in the T3 main channel, which mainly affect the geometric matching, spatial permeability, and conformational convergence of R-CHCM after entering the catalytic chamber along the T3 channel; L227 was identified as a diversion valve site, which mainly affects whether the substrate enters the Ser201 catalytic chamber from the end of T3 or is deflected into the non-productive side chamber of Pocket-S.

[0214] 9. Mutation Construction and Experimental Validation Procedure

[0215] Based on the site priority score Mi, Y228C / D253A was selected as the starting background. The candidate sites for Est13 were ranked, yielding the following priority sites: The first category consists of T3 main channel gate sites, including A230, F233, I254, V257, and T258. These sites primarily regulate the geometric matching, spatial permeability, and conformational convergence of R-configuration substrates entering the catalytic cavity along the T3 channel. The second category consists of split-valve sites, including L227. This site determines the split direction of the substrate between the T3 productive channel and the Pocket-S non-productive side cavity, and is a crucial site controlling the R / S selectivity direction and selectivity reversal. The third category consists of secondary gating or dynamically coupled sites, such as T137 and H139.

[0216] Transform the correctly sequenced mutant plasmid into a suitable expression host for expression. If using the His-tagged Est13 sequence, protein purification can be performed using Ni-NTA affinity chromatography or equivalent methods; alternatively, activity and selectivity can be determined using crude enzyme solution, cell lysis supernatant, or purified enzyme, depending on the experimental purpose. All mutants should be measured under the same expression, induction, purification, and reaction conditions to ensure the comparability of EEP data.

[0217] Racemic CHCM was used as a substrate for enzymatic resolution. The substrate concentration, enzyme dosage, buffer pH, temperature, organic co-solvent ratio, and reaction time were kept consistent throughout the reaction system. After the reaction, the product composition was analyzed using gas chromatography, liquid chromatography, or chiral column chromatography, and the enantiomeric excess value (eep) was calculated. If enzyme activity was simultaneously measured, the substrate conversion rate per unit amount of enzyme per unit time could be recorded to analyze the trade-off between increased selectivity and catalytic efficiency.

[0218] Table 3. Est13 esterase mutants

[0219]

[0220] The above experimental verifications demonstrate that this invention does not merely output residues near the static pocket, but can identify sites of different functional types based on conformational flux allocation: A230, F233, I254, V257, and T258 belong to T3 gate sites that enhance the productive flux of the R conformation; L227 belongs to a shunt valve site that regulates R / S path shunting and selective inversion. Therefore, this invention can be used to transform complex molecular dynamics trajectories into directly experimentally verifiable mutation site screening schemes.

[0221] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope of the present invention.

Claims

1. A method for screening enzyme enantioselective mutation sites based on conformational flux, characterized in that, Includes the following steps: Obtain the three-dimensional structure of the target enzyme, the enantiomer of the substrate to be compared, and the reference residues of the active site of the target enzyme; Based on the three-dimensional structure of the target enzyme, candidate channels leading from the active site to the protein surface are identified, and channel data for each candidate channel is generated. Based on channel data, ingress sampling simulation, binding state repeated sampling simulation, and piecewise random external force escape simulation are performed for each enantiomer to obtain ingress sampling simulation data, binding state repeated sampling simulation data, and piecewise random external force escape simulation data, respectively. State discretization is then performed to construct a conformational state diagram and calculate the conformational flux. The conformational state diagram includes at least one or more of the following: bulk state, surface capture state, channel ingress state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and exit state. The conformational flux includes at least one or more of the following: ingress flux, near-attack conformational flux, non-productive residence penalty term, exit flux, and gating score. Based on the difference in conformational flux between enantiomers, an enantioselectivity score and a site priority score are calculated, and at least one of the following is output: candidate mutation site, recommended mutation direction, and enantioselectivity prediction.

2. The method as described in claim 1, characterized in that, The channel data includes at least centerline data, bottleneck radius data, and ROI residue set; and / or, The inlet sampling simulation data includes at least surface trapping pits, the main inlet channel, and the intermediate metastable region; and / or, The combined state repeated sampling simulation data includes at least the trajectory data of the substrate in the active pocket pre-organization, near-attack conformation state, and non-productive side cavity state; And / or, The segmented random external force escape simulation data includes at least the escape path of the substrate or product from the active cavity to the protein surface, the escape time, and the probability of use of the exit channel.

3. The method as described in claim 1, characterized in that, The inlet sampling simulation is a molecular dynamics simulation of ligand replication flooding. The molecular dynamics simulation of ligand replication flooding includes: arranging multiple substrate copies in a preset spatial shell outside the target enzyme solvent-exposed surface, and statistically analyzing events during molecular dynamics in which the substrate is captured by the surface, enters the channel ROI, or enters the active pocket.

4. The method as described in claim 1, characterized in that, The segmented random external force escape simulation uses the distance-contact number dual criterion as the escape criterion; the distance-contact number dual criterion simultaneously satisfies the following conditions: the robust statistical value of the distance between the centroid of the substrate or product and the centroid of the active site reference group in the tail window is greater than a preset distance threshold, and the robust statistical value of the number of heavy atom contacts between the substrate or product and the active site reference group in the tail window is less than a preset contact number threshold.

5. The method as described in claim 1, characterized in that, The criteria for the near-attack conformation state include at least two or more of the following four: the distance threshold between the catalytic nucleophile atom and the substrate reaction center atom; the nucleophilic attack angle threshold; the hydrogen bond threshold between the substrate carbonyl oxygen and the oxygen anion hole hydrogen-donating residue; and the spatial orientation threshold of the substrate relative to the catalytic triplet.

6. The method as described in claim 1, characterized in that, The gating score is determined by combining at least two types of information: dynamic cross-correlation changes between the catalytic residue region and the channel ROI region; time-series changes in the bottleneck radius of the candidate channel; state transitions of the side chain conformation angle, main chain dihedral angle, or side chain conformation cluster of the gating residue; and changes in pocket water molecule occupancy or local solvent accessibility.

7. The method as described in claim 1, characterized in that, The site priority score is a weighted combination of at least the following indicators or their monotone equivalent transformations: the frequency with which the site participates in channel bottlenecks; the contact asymmetry between the site and different enantiomers in NAC or NPV states; the difference in the associated motion of the site with catalytic residues or gated residues; the sensitivity of the ESS after applying local geometric perturbations or in situ virtual mutations to the site; and the site's structural role factor, which indicates whether the site belongs to a channel ROI, a lateral cavity ROI, near an oxygen anion hole, a gate frame site, and / or a gated site.

8. A system for screening enzyme enantioselective mutation sites based on conformational flux, characterized in that, include: The structure input module is used to obtain the three-dimensional structure of the target enzyme, the enantiomer of the substrate to be compared, and the reference residues of the active site of the target enzyme. The channel recognition module, based on the three-dimensional structure of the target enzyme, identifies candidate channels leading from the active site to the protein surface and generates channel data for each candidate channel. The ingress sampling module is used to perform ingress sampling simulation for each enantiomer based on channel data to obtain ingress sampling simulation data; The associative state resampling module is used to perform associative state resampling simulation for each enantiomer based on channel data, and obtain associative state resampling simulation data; The segmented random external force escape module is used to perform segmented random external force escape simulation for each mapper based on channel data, and obtain segmented random external force escape simulation data. The state diagram construction module is used to discretize the obtained inlet sampling simulation data, combined state repeated sampling simulation data, and segmented random external force escape simulation data and construct a conformational state diagram; the conformational state diagram includes at least one or more of the following: bulk phase state, surface capture state, channel inlet state, active pocket pre-organized state, near-attack conformational state, non-productive pocket state, and outlet state. The conformation flux calculation module is used to calculate the conformation flux of each enantiomer based on the conformation state diagram; the conformation flux includes at least one or more of the following: inlet flux, near-attack conformation flux, non-productive residence penalty, outlet flux, and gating score; The site priority scoring module is used to calculate enantioselectivity scores and site priority scores based on the differences in conformational flux between enantiomers. The mutation site output module is used to output at least one of the following based on site priority scoring: candidate mutation sites, recommended mutation directions, and enantioselectivity predictions.

9. The application of shunt valve sites and / or channel gate sites in the construction of Est13 esterase mutants, characterized in that, The diversion valve location includes at least Leu227; The channel gate location includes at least Ala230, Phe233, Ile254, Val257, and Thr258.

10. An Est13 esterase mutant, characterized in that, Obtained from the amino acid sequence shown in SEQ ID NO.1 through any of the following mutations: (1) Y228C / D253A; (2) Y228C / D253A / A230Y; (3) Y228C / D253A / I254V; (4) Y228C / D253A / F233K; (5) Y228C / D253A / V257A; (6) Y228C / D253A / V257M; (7) Y228C / D253A / V257S; (8) Y228C / D253A / T258E; (9) Y228C / D253A / L227M; (10) Y228C / D253A / L227F; (11) Y228C / D253A / L227W; (12) Y228C / D253A / L227R; (13) Y228C / D253A / L227N; (14) Y228C / D253A / L227H; (15) Y228C / D253A / L227K; (16) Y228C / D253A / L227S.

Images