Method for predicting protein oxidation risk based on multiscale spatial topology and biochemical interaction

Abstract

The application discloses a protein oxidation risk prediction method based on multi-scale space topology and biochemical interaction, and relates to the cross technical field of food science and computational biology. Through molecular dynamics simulation, the original physical characteristics of the solvent accessible surface area, the center of mass spatial distance and the root mean square fluctuation of the protein full residue under different macromolecular crowded environments are obtained. Secondly, a sliding window is constructed with each residue as the center, the local environment smooth descriptor and the environment relative mutation descriptor are calculated, the spatial local topology descriptor is obtained, the macroscopic crowded concentration is multiplied by the selected spatial local topology descriptor, and the cross-dimensional biochemical interaction descriptor is constructed to simulate the modulation effect of the environment on the microstructure. The above characteristics are integrated to form a multi-scale feature set, a machine learning model is trained by using the known oxidation state site label, and SMOTE sample balance and regularization pruning are combined. The multi-scale features of the to-be-predicted site are input into the model, and the oxidation risk probability and 3D space mapping atlas are output.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of food science and computational biology, specifically to a method for predicting protein oxidation risk based on multi-scale spatial topology and biochemical interactions. Background Technology

[0002] Protein oxidation is one of the main causes of protein quality deterioration during food processing and storage. In food systems, the electrophilic attack of reactive oxygen species (ROS) on protein side chains leads to covalent modification of amino acid residues, causing protein cross-linking, aggregation, decreased solubility, and loss of functional properties, severely affecting the nutritional value and sensory quality of food. Therefore, accurately predicting protein oxidation risk sites is of great significance for optimizing processing technology, extending shelf life, and ensuring food safety.

[0003] Currently, the detection of protein oxidation sites mainly relies on experimental methods, including mass spectrometry, immunoblotting, and fluorescent labeling. However, these methods suffer from high costs, long processing times, and low throughput, and can only perform post-oxidation detection, failing to provide early warning. The limitations of these methods are particularly pronounced for complex systems commonly found in food processing, such as the impact of macromolecular crowding in high-concentration mixtures on protein oxidation.

[0004] In recent years, machine learning-based computational prediction methods have gradually become a new approach for predicting protein oxidation sites. Existing computational methods mainly rely on the single-point physicochemical properties of amino acid residues (such as hydrophobicity and isoelectric point) as input features, and classify them into oxidation / non-oxidation categories using a classifier. However, these methods have the following technical drawbacks: First, single-point physicochemical properties only reflect the static chemical properties of the residues themselves and cannot capture the synergistic regulatory effect of the local microenvironment on oxidation risk; second, they do not consider the nonlinear modulation of micro-oxidation reactions by macro-environmental factors (such as macromolecular crowding effects), resulting in insufficient generalization ability of the model under different environmental conditions; third, under extreme sample imbalance conditions (such as a ratio of oxidized to non-oxidized sites as high as 1:46), traditional classifiers are prone to majority class bias, resulting in serious underreporting of rare oxidation targets.

[0005] Therefore, there is an urgent need for a protein oxidation risk prediction method that can integrate multi-scale spatial topological information and macro-micro biochemical interaction effects to address the shortcomings of existing technologies in feature expression and environmental adaptability. Summary of the Invention

[0006] The purpose of this invention is to overcome the technical problems existing in the prior art, and to propose a protein oxidation risk prediction method based on multi-scale spatial topology and biochemical interactions, comprising the following steps: (1) By molecular dynamics simulation, the molecular dynamics trajectory of the target protein under different macromolecular crowding environment parameters is obtained, and the original physical characteristics of each amino acid residue in the whole sequence are extracted from it. The original physical characteristics include solvent-accessible surface area, centroid spatial distance and root mean square fluctuation. (2) Based on the original physical characteristics, for each target amino acid residue in the sequence, a sliding window is constructed with it as the center, and the local environment smoothing descriptor and the environment relative mutation descriptor within the window are calculated to obtain the spatial local topological descriptor; (3) Multiply the environmental parameter representing the macro-crowding concentration with the environmental relative mutation descriptor and the local environmental smoothing descriptor in the spatial local topology descriptor respectively to construct a cross-dimensional biochemical interaction descriptor; (4) Integrate the original physical features, spatial local topological descriptors and cross-dimensional biochemical interaction descriptors to form a multi-scale feature set, and use protein site data with known oxidation states as labels to train the machine learning model. (5) Input the multi-scale feature set of the protein site to be predicted into the trained machine learning model and output the oxidation risk prediction result of the site.

[0007] Further, the local environment smoothing descriptor in step (2) is obtained by calculating the average value of the original physical characteristics of the effective residues within the calculation window; the environment relative mutation descriptor is the original physical characteristic value of the target residue minus the local environment smoothing descriptor.

[0008] Further, the cross-dimensional biochemical interaction descriptor mentioned in step (3) includes: a crowding-forced perturbation descriptor, which is obtained by multiplying macroscopic environmental parameters with a relative abrupt change environmental descriptor of root mean square fluctuations; and an environment-mediated exposure descriptor, which is obtained by multiplying macroscopic environmental parameters with a local environment smoothing descriptor of solvent-accessible surface area.

[0009] Furthermore, the machine learning model described in step (4) is the XGBoost model, and the training set is balanced with positive and negative samples using the SMOTE algorithm or SMOTETomek composite sampling algorithm before training. Under the hierarchical cross-validation framework, the intelligent hybrid pipeline ensures that the sampling only applies to the training set.

[0010] Furthermore, the multi-scale feature set described in step (4) is filtered by feature importance, and the core features retained include: polarity, local environment smoothing descriptor of solvent-accessible surface area, local environment smoothing descriptor of centroid spatial distance, environmental relative abrupt change descriptor of root mean square fluctuation, hydrophobicity, isoelectric point, original features of root mean square fluctuation, and original features of solvent-accessible surface area.

[0011] Furthermore, the oxidation risk prediction results output in step (5) include: the oxidation probability value of each site, the whole sequence oxidation risk heatmap under different crowding concentrations, and the temperature factor distribution map that maps the oxidation probability to the three-dimensional structure of the protein.

[0012] Furthermore, after training in step (4), the SHAP interpretation framework is used to perform global feature importance analysis on the trained machine learning model, and output the contribution ranking of each feature to oxidation risk and the direction of directional regulation.

[0013] Furthermore, the protein site tags with known oxidation states mentioned in step (4) are obtained by liquid chromatography-tandem mass spectrometry. Positive samples are amino acid residues that have been confirmed to have undergone oxidation modification, and negative samples are residues that have not undergone oxidation modification.

[0014] Furthermore, in step (2), when calculating the sliding window descriptor, dynamic truncation adjustment is adopted for the N-terminal and C-terminal residues of the protein sequence, and the calculation is based on the actual number of effective residues existing in the window.

[0015] Furthermore, the training process in step (4) adopts the repeated hierarchical five-fold cross-validation protocol, which is repeated 5 times, for a total of 25 independent training and validation cycles.

[0016] Compared with the prior art, the present invention has the following beneficial technical effects: This invention addresses the technical bottlenecks of high cost and long cycle time in protein oxidation detection during food processing and storage by proposing a multi-scale spatial topology and biochemical interaction-based method for predicting protein oxidation risk. Validated in the oat globulin (12SOG) system, under extreme imbalance conditions with a positive-to-negative sample ratio as high as 1:46, this invention, through XGBoost combined with a multi-scale feature model, achieved excellent performance with an average AUC of 0.9112 ± 0.0281 and an accuracy of 0.9185, significantly outperforming traditional machine learning models.

[0017] By constructing a sliding window spatial local topological descriptor, the synergistic regulatory effect of 25 residues before and after the protein oxidation site was successfully captured, which increased the model AUC from 0.67 to 0.83. After further introducing a cross-dimensional biochemical interaction descriptor, the AUC exceeded 0.91 and the standard deviation was extremely small, proving that the model has extremely high generalization ability in complex and crowded environments.

[0018] By using the SHAP framework to decode the model decision logic, we quantitatively reveal that hydrophobicity and local average solvent-accessible surface area are the two core factors driving oxidation. We also discovered for the first time that crowded-forced perturbation features can accurately identify structural tension regions, providing clear targets for the design of anti-oxidation strategies.

[0019] The model predictions revealed that the oxidation risk changes non-linearly with increasing crowding concentration. Moderate crowding (D40) creates a shielding effect, reducing the number of high-risk sites to 97; while extreme crowding leads to the collapse of the defense, causing the number of high-risk sites to surge to 130. This finding provides a quantitative basis for optimizing food processing techniques.

[0020] By mapping the predicted probability to a 3D temperature factor, an oxidation topological fingerprint spectrum was constructed, which intuitively presents the spatial evolution trajectory of high-risk regions and can directly guide targeted modification and formulation design. Attached Figure Description

[0021] Figure 1 The configurations of the 12SOG molecular dynamics simulation system under different macromolecular crowding environment parameters are shown, where AF corresponds to D0, D20, D40, D60, D80 and D100 gradient crowding environments, respectively. Figure 2 The flowchart shows the sample imbalance handling and optimization test based on the SMOTE algorithm. Figure 3 A flowchart of microenvironment perception feature engineering based on sequence-scale spatial exhaustive search; Figure 4 The trend of the impact of different sliding window radii (w) on the average AUC of each machine learning model; Figure 5 This is a radar chart showing the performance of the upgraded XGBoost model based on SMOTETomek sampling and repeated validation. Figure 6 Ranking of contributions to 51-dimensional original topological features and selection of the Top 8 golden core features; Figure 7 Radar chart for multidimensional performance evaluation of the XGBoost model after feature engineering and Occam's razor pruning; Figure 8 A schematic diagram illustrating the construction principle and macro-micro collaborative response logic of the cross-dimensional biochemical interaction descriptor (CBI); Figure 9 To establish a feature importance ranking for the final champion model after incorporating biochemical interaction terms and performing feature pruning; Figure 10 Radar chart for multidimensional performance evaluation of the XGBoost model after feature engineering 2 and Occam's razor pruning; Figure 11 Radar chart for multidimensional performance evaluation of the XGBoost model in Feature Engineering 3; Figure 12 Radar chart for multidimensional performance evaluation of the XGBoost model in feature engineering 4; Figure 13 Different feature engineering approaches to improve XGBoost prediction performance; Figure 14 The ROC curve is used to comprehensively verify the predictive performance of the feature engineering model 2. Figure 15 Decoding the global feature importance of SHAP; Figure 16 The core driving feature attribution beehive diagram; Figure 17 This is a macroscopic regulation of the protein oxidation probability by the concentration of macromolecular crowding. Figure 18 Statistics on the number of high-risk oxidation sites under different crowding concentrations; Figure 19 Heatmap of the oxidation risk evolution of the full sequence of oat globulins under different crowding environments; Figure 20 3D structure mapping of the oxidation risk of oat globulin under macromolecular crowding. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0023] I. Experimental Design 1. Original feature set The chemical properties of amino acid residues themselves, including isoelectric point, polarity, side chain molecular weight and hydrophobicity index, are shown in Table 1, which contains the core physicochemical properties of 20 amino acids.

[0024] Table 1 alanine Ala Nonpolar (0) 89.09 6.02 1.8 Arginine Arg The polarity carries a positive charge (2) 174.20 10.76 -4.0 Asparagine Asn Polarity is uncharged (1) 132.12 5.41 -3.5 Aspartic acid Asp The polarity carries a negative charge (2) 133.10 2.97 -3.5 Cysteine Cys Polarity is uncharged (1) 121.15 5.02 2.5 glutamic acid Glu The polarity carries a negative charge (2) 147.13 3.22 -3.5 glutamine Gln Polarity is uncharged (1) 146.15 5.65 -3.5 glycine Gly Polarity is uncharged (1) 75.07 5.97 -0.4 Histidine His The polarity carries a positive charge (2) 155.15 7.59 -3.2 Isoleucine Ile Nonpolar (0) 131.17 6.02 4.5 Leucine Leu Nonpolar (0) 131.17 5.98 3.8 Lysine Lys The polarity carries a positive charge (2) 146.19 9.74 -3.9 Methionine Met Nonpolar (0) 149.21 5.75 1.9 Phenylalanine Phe Nonpolar (0) 165.19 5.48 2.8 proline Pro Nonpolar (0) 115.13 6.30 -1.6 Serine Ser Polarity is uncharged (1) 105.09 5.68 -0.8 threonine Thr Polarity is uncharged (1) 119.12 6.53 -0.7 Tryptophan Trp Nonpolar (0) 204.23 5.89 -0.9 Tyrosine Tyr Polarity is uncharged (1) 181.19 5.66 -1.3 Valine Val Nonpolar (0) 117.15 5.97 4.2 Physical characteristics: Molecular dynamics simulations were used to obtain trajectory files and other data of 12SOG under different crowding conditions. Physical characteristics of neighboring residues were extracted, including solvent accessible surface area (SASA), centroid spatial distance (Distance), and root mean square fluctuation (RMSF).

[0025] 2. Spatial Topological Descriptors The identification of protein oxidation sites depends not only on the physicochemical properties of the residues themselves, but also on the local spatial microenvironment in which the site is located. To quantify the synergistic regulatory effect of neighboring residues on oxidation risk, this invention introduces a spatial topological descriptor based on a sliding window w. This project aims to explore the spatial correlation range of oxidation susceptibility at the protein sequence scale. Let the index of the target residue in the protein sequence be i, and a sliding window of length 2w+1 is established. For the dynamic physical feature f∈{SASA,Distance,RMSF}, the following two types of derived descriptors are constructed: (1) Local Mean Descriptor ) This descriptor represents the area surrounding the target residue. The average physical background within the range. It can smooth out random noise at a single point and characterize the overall topological features of the region. Its calculation formula is: ; in, The generated local environment smoothing descriptor represents the mean local feature value of the i-th residue at a radius of w; i is the index position of the target amino acid residue in the protein sequence (Site). ID f represents the selected original physical feature type, including SASA, Distance, and RMSF; w is the sliding window radius parameter (WindowRadius), which determines the radius range of spatial perception; j is the residue index falling within the sliding window region [iw, i+w]. This represents the original observation of the j-th residue on physical feature f; The total number of effective residues actually involved in the calculation within the sliding window is automatically truncated and adjusted at the sequence edges.

[0026] The code uses `groupby('Condition')` to ensure that calculations are performed within the same experimental group, and uses `rolling(window=2*w+1,center=True,min)` to ensure that calculations are performed within the same experimental group. periods =1).mean() is implemented.

[0027] Boundary handling: For the N-terminal and C-terminal residues of the sequence, the algorithm sets a minimum... periods =1, ensuring that even in the event of a neighboring cell being missing, the actual number of valid residues can still be used. Calculate the accurate regional mean.

[0028] (2) Local Delta Descriptor ) This descriptor is used to quantify the deviation of a target residue from its local microenvironment, aiming to identify structural defect sites, such as sites exhibiting abnormal flexibility in globally rigid regions. The calculation formula is as follows: ; in, This is the generated environmental relative mutation descriptor, representing the degree of abnormal deviation of site i relative to the background; The original physical observation value of target site i on feature f; This is the local environment smoothing descriptor calculated at the corresponding window size w.

[0029] df in the corresponding code w [delta_col_name]=df w [col]-df w [mean_col_name]. This item can capture local conformational mutations of residues and is a key high-frequency feature for predicting oxidation risk.

[0030] 3. Cross-dimensional Biochemical Interaction Descriptors (CBI) Topological structures undergo forced reshaping and perturbation. To quantify the synergistic effect between macroscopic environmental factors and microscopic structural traits, this invention constructs a cross-dimensional biochemical interaction descriptor. Macroscopic experimental parameters (Condition, i.e., the concentration of the congestant dextran) are used as feature weights and dynamically applied to the local spatial descriptor, thereby simulating the amplification or inhibition of the oxidation risk at local exposed or flexible protein sites under different congestion environments. The following two types of descriptors are constructed: (1) Crowding-induced Perturbation Descriptor ) This descriptor quantifies the modulating effect of a crowded environment on the overall exposure level of a local area, reflecting the amplification effect of the macro-environment on the background of solvent accessibility. Its calculation formula is as follows: ; in, The generated crowding-forced perturbation descriptor represents the dynamic damage risk of the i-th residue under environmental pressure; This is a macroscopic environmental parameter (Condition), representing the number of dextran molecules in the molecular dynamics simulation system. In this embodiment, its specific value corresponds to six set gradients (D0, D20, D40, D60, D80, D100), which respectively represent 0, 20, 40, 60, 80, and 100 dextran molecules added to the simulation box. is the environmental relative mutation descriptor of the i-th residue at the window scale w; f is the corresponding physical feature, which in this embodiment mainly refers to the root mean square fluctuation (RMSF).

[0031] The implementation instruction is: in Python code, it is df['Crowding_Perturbation']=df['Condition']*df['RMSF_Delta_w25'].

[0032] (2) Environment-mediated Exposure Descriptor ) This descriptor quantifies the modulating effect of a crowded environment on the overall exposure level of a local area, reflecting the amplification effect of the macro-environment on the background of solvent accessibility. Its calculation formula is as follows: ; in, The generated environment-mediated exposure descriptor is used to assess overall oxidative susceptibility in an environmental context; For macro-environmental parameters (Condition); is the local environment smoothing descriptor for the i-th residue at window scale w; f is the corresponding physical feature, mainly referring to the solvent accessible surface area (SASA).

[0033] The implementation instruction is: In Python code, it is represented as df['Dynamic_Hydro_Exposure']=df['Hydrophobicity']*df['SASA_Mean_w25'].

[0034] The introduction of cross-dimensional interaction descriptors enables a paradigm shift in models, moving from static feature recognition to dynamic response sensing and capturing synergistic effects, rather than relying on a single... or It can often only explain a portion of the variance, while It can effectively identify deep sites that fluctuate drastically due to environmental changes, leading to changes in oxidation risk.

[0035] 4. Environmental Boundary Probing and Heterogeneity Descriptors To clarify the boundaries of Occam's Razor in predicting protein oxidation risk and prevent the model from experiencing false booms or overfitting due to disordered feature stacking, this invention designs a system of environmental boundary detection descriptors. This system serves as a negative control for the model, aiming to test the following three scientific hypotheses from mathematical and physicochemical perspectives: 1. Is oxidation driven by drastic fluctuations in local physical states (heterogeneity)? 2. Is oxidation sensitive to scale deviations (structural faults) between the macroscopic environment and microscopic local conditions? 3. Is oxidation dominated by a macroscopically continuous static chemical environment (chemical halo)? By introducing these detection indicators, the effectiveness boundaries of feature engineering can be scientifically and clearly defined. Three types of boundary detection descriptors are dynamically constructed within the scale space: (1) Microenvironment heterogeneity descriptor ) This descriptor uses the unbiased standard deviation within a sliding window to quantify the rate of drastic change (fluctuation) of the local physical state of the target residues, in order to examine whether the tension at the rigid-flexible interface of the structure will induce oxidation. Its calculation formula is as follows: ; in, The generated microenvironment heterogeneity descriptor represents the local feature standard deviation of the i-th residue at a radius of w; i is the index position of the target amino acid residue in the protein sequence (Site). ID ); f is the selected type of original physical feature ( w is the window radius parameter (WindowRadius); j is the area within the sliding window. Residue index within; This represents the original observation of the j-th residue on physical feature f; The calculated local environment smoothing descriptor (regional mean); This represents the total number of effective residues actually used in the calculation within the sliding window.

[0036] Implementation instructions: In the code, use `groupby('Condition')` to limit the experimental groups, and use `rolling(window=2*w+1,center=True,min)` to... periods =2).std().fillna(0) is used to implement this. Setting min...periods =2 because calculating the standard deviation requires at least two valid data points, and it is automatically filled with 0 when a single edge residue cannot be calculated.

[0037] (2) Multi-scale concentric ratio descriptor ) This descriptor captures relative disparities or structural defects in three-dimensional conformations (such as abnormally exposed residues deeply embedded in a hydrophobic core) by calculating the ratio of physical quantities between a minimal microscopic window and a maximal macroscopic window. The calculation formula is as follows: ; in, To generate multi-scale nested ratio descriptors, quantifying the degree of deviation between the micro-environment and the macro-background; To achieve the microscopic sliding window radius The local environment smoothing descriptor calculated below (in this embodiment, it is taken as...) ); To the macroscopic sliding window radius The local environment smoothing descriptor calculated below (in this embodiment, it is taken as...) f represents the corresponding original physical feature type; This is a small value for regularization (RegularizationTool), set in the code as follows. This is used to prevent program division by zero errors when the denominator is 0.

[0038] Implementation instructions: The corresponding code calculates w separately. micro =2 and w macro =25 rolling average, perform micro mean / (macro mean The +epsilon operation is performed and dynamically injected into the feature space.

[0039] (3) Static Chemical Halo Descriptor ) To verify whether oxidation reactions depend on specific single-point biochemical properties or on the surrounding continuous regional chemical environment, this invention extends the sliding window smoothing mechanism to static chemical features, constructing a chemical halo descriptor. Its calculation formula is as follows: ; in, The generated static chemical halo descriptor represents a continuous chemical microenvironment with radius w around the i-th residue; The intrinsic property value of the j-th residue in static chemical property c ( ).

[0040] Implementation instructions: In the corresponding code, the static chemical feature columns are traversed, and df.groupby('Condition')[col].transform(lambdax:x.rolling(...).mean()) is executed to dynamically generate the result.

[0041] II. Experimental Procedure 1. Raw data collection and congestion environment model construction (1) Molecular dynamics simulation The initial structure of 12SOG was obtained from UniProt, and molecular dynamics simulations were performed using the GROMACS software package under different concentrations of dextran. Topological files for 12SOG and dextran were constructed using the CHARMM36 force field. The dynamics simulations were conducted in the TIP3P water model with an appropriate amount of Na added. + and Cl - Ions (1M NaCl). After system construction, energy minimization (rapidest descent method) and equilibrium simulations (NVT and NPT ensembles, 1 ns each) were performed sequentially to ensure the simulated system reached a stable state at 310.15 K and 1 bar. After removing positional constraints, productive molecular dynamics simulations were performed on six independent crowded systems (D0, D20, D40, D60, D80, D100) in the NPT ensemble, with each system simulated for 100 ns and a time step of 2 fs. Trajectory files and other data for 12SOG under different crowded conditions were obtained, and physical characteristics of neighboring residues were extracted, including solvent accessible surface area (SASA), centroid spatial distance (Distance), and root mean square fluctuation (RMSF).

[0042] By analyzing trajectory files (D0, D20, D40, D60, D80, D100) under different congestion environments, the simulation architecture is shown in the diagram. Figure 1 The dynamic physical feature space (X) of all residues in the whole sequence under six-dimensional environments, including SASA, Distance and RMSF, is extracted as the benchmark input for subsequent multidimensional feature engineering evolution.

[0043] (2) Identification of actual oxidation sites Sample preparation: Weigh 1.00 g of oat globulin and dissolve it in two 100 mL solutions of 6% NaCl: one containing 60 mg / mL of 200 kDa dextran (crowded group), and the other without dextran (dilution group). Add AAPH to each solution to a final concentration of 40 mmol / L and react at 37°C in the dark for 6 h. After the reaction, cool the samples in an ice bath, and then dialyze them in deionized water at 4°C for 48 h using a dialysis bag with a molecular weight cutoff of 3.5 kDa to remove salt ions and AAPH. Freeze-dry the dialyzed samples for 48 h, and use the resulting powder for subsequent analysis.

[0044] Modification site identification: The lyophilized oat globulin sample powder was dissolved in a reconstitution solution (8M urea, 100mM Tris-HCl, pH 8.5). After determining the protein concentration using the BCA method, an equal volume of protein sample was taken and incubated at 37°C for 1 hour with the addition of TCEP and CAA to complete the reduction and alkylation treatment. Subsequently, the urea concentration was diluted to below 2M using 100mM Tris-HCl solution. Trypsin was added at a ratio of 1:50 (w / w) between enzyme and protein, and the mixture was incubated overnight at 37°C with shaking for enzymatic digestion. The reaction was terminated the next day by adding TFA. After centrifugation at 12000g, the supernatant was collected and desalted using a self-made SDB desalting column. After vacuum drying, the sample was stored at -20°C for later use. The Bruker timesTOFPro liquid chromatography-mass spectrometry system was used. Mass spectrometry analysis was performed using a timsTOFPro captured ion mobility mass spectrometer (Bruker, Germany) coupled with an UltiMate 3000 RSLC nano system (Thermo Fisher Scientific, USA). Peptide samples were automatically injected and bound to a C18 trapping column (75 μm × 2 cm, 3 μm), followed by separation on a reversed-phase analytical column (75 μm × 15 cm, 1.7 μm, Ion Opticks). Mobile phase A was 0.1% formic acid aqueous solution, and mobile phase B was 0.1% formic acid in acetonitrile solution. The mass spectrometer was operated in parallel cumulative tandem mass spectrometry (PASEF) mode. The capillary voltage was set to 1500 V, and the scan ranges for both the first stage (MS1) and second stage (MS2) were 100–1700 m / z. The reciprocal of ion mobility (1 / K0) was set to 0.75–1.4 V·s / cm. 2 Ion accumulation time and ramp time were both set to 166 ms. The single acquisition cycle time was 1.88 s, including one MS1 full scan and ten PASEFMS2 scans. Single-charge precursor ions were excluded by detecting ion mobility characteristics. Fragmentation was triggered when the precursor ion signal intensity exceeded 1000, with a target intensity set at 20000. Collision energy was linearly adjusted with ion mobility: when 1 / K0 = 1.6 V·s / cm 2At 59 eV, it linearly decreases to 1 / K0 = 0.6 V·s / cm. 2 The time was 20 eV. The dynamic exclusion time was set to 0.4 min. The raw data obtained by mass spectrometry acquisition were searched and analyzed using MaxQuant (v2.2.0.0) software and the built-in Andromeda algorithm. The search database used was a custom Fasta file constructed based on the oat globulin (12SOG) sequence in UniProt and related literature reports.

[0045] After the raw mass spectrometry data were retrieved and identified, amino acid residue sites that were confirmed to have undergone oxidation modification were defined as positive samples (Label=1), and sites that had not undergone oxidation modification were defined as negative samples (Label=0).

[0046] In this invention, mass spectrometry identification was performed under two typical conditions: 0 mg / mL (dilution group) and 60 mg / mL (crowded group), providing core supervised truth labels (Ground Truth) for subsequent machine learning algorithms.

[0047] 2. Extraction of the original feature set and establishment of the benchmark dataset (1) The physical and biochemical connotations of the seven major benchmark characteristic systems To characterize the inherent properties and spatial state of each amino acid site in the entire protein sequence from multiple dimensions, this invention constructs a baseline feature space containing seven single-point original features. Based on its acquisition channels and physicochemical properties, these seven characteristics can be divided into two core components.

[0048] A. Static biochemical and physicochemical characteristics (4 items, derived from classic bioinformatics databases) Hydrophobicity: determines whether residues tend to be hidden in the hydrophobic core inside the protein or exposed to the aqueous phase, directly affecting the probability of contact between the site and oxidative free radicals.

[0049] Isoelectric point (pI): The charged state of the reactive residues under a specific pH environment, which regulates the local charge microenvironment.

[0050] Molecular weight Weight ): Characterizes the spatial volume occupancy effect of residue side chains.

[0051] Polarity: Characterizes the non-uniformity of electron distribution in the side chains of residues and is the core essence that determines the stability of chemical covalent bonds.

[0052] B. Dynamic spatial topological features (3 items, derived from molecular dynamics simulation trajectories) Solvent-accessible surface area (SASA): quantifies the actual absolute area of ​​a single residue exposed to solvent. The larger the area, the higher the probability of direct attack by oxidants in the macroscopic environment.

[0053] Centroid distance: Characterizes the distance from the geometric centroid of the residue side chain to the geometric center of the protein, used to quantitatively characterize the degree of site burial / exposure.

[0054] Root mean square fluctuations (RMSF): characterize the local flexibility and structural amplitude of a single amino acid residue during kinetic simulation. Violent fluctuations usually indicate conformational instability.

[0055] (2) Topological alignment of original data across groups Raw data collected from multiple sources (mass spectrometry, MD simulation, bioinformatics databases) is first broadly segregated by macroscopic experimental group (Condition) before being directly fed into the classifier. Then, within each group, data is strictly classified according to amino acid position index (Site). ID The residues are arranged in ascending order. This mechanism cuts off cross-contamination of residue features between different concentration groups, laying a solid matrix structure foundation for subsequent calculations.

[0056] 3. Machine Learning Model Training and Implementation Steps (1) Selection of training sample set and construction of benchmark matrix This invention first extracts seven established single-point original physicochemical features and integrates them with the macroscopic crowding environmental variable Condition to jointly construct an 8-dimensional benchmark input feature space. To enable supervised learning and training, the program retrieves mass spectrometry identification tags. Unlabeled blank sites are filtered out, and only valid samples that are confirmed to be oxidized (labeled 1) and unoxidized (labeled 0) are retained. Finally, a supervised benchmark training matrix with fixed dimensions is constructed as the initial input to the machine learning model.

[0057] (2) Model training process To address the unique nonlinear and multidimensional feature interaction characteristics of bioinformatics data, this invention establishes a parallel evaluation process that includes four mainstream classifiers: Tree ensemble model: Gradient boosting decision tree (XGBoost) and random forest (RandomForest) are selected to capture spatial features by taking advantage of their natural non-linearity.

[0058] Traditional statistical and kernel function models: Logistic Regression and Support Vector Machine (SVM) are selected. Since these two algorithms are extremely sensitive to the dimensions (absolute numerical magnitude) of features, this invention constructs a StandardScaler standardization pipeline before the algorithm to achieve unbiased de-dimensionalization of features. Simultaneously, the SVM is configured with probability=True to enable it to output continuous predicted probabilities for subsequent robustness evaluation.

[0059] (3) Stratified 5-FoldCV protocol Given the significant class imbalance in mass spectrometry experimental data, this invention employs a Stratified 5-Fold Cross-Validation strategy to prevent false booms caused by the randomness of data partitioning. This method rigorously divides the dataset into five mutually exclusive subsets and mandates that the ratio of positive to negative samples within each fold (training and testing sets) strictly conforms to the overall distribution. The training process is executed 5 times, alternating between supervised training with four folds and blind testing with one fold, thus avoiding data cross-contamination.

[0060] (4) Benchmark test results and analysis of algorithm selection Based on the original 8-dimensional baseline features, the baseline performance of each machine learning model under five-fold cross-validation is shown in Table 2.

[0061] Table 2 RandomForest 0.6760±0.0582 0.1077 0.1610 0.4757 0.8276 XGBoost 0.6675±0.0961 0.1846 0.2536 0.4893 0.8251 LogisticRegression 0.6476±0.0511 0.0000 0.0000 0.0000 0.8374 SVM 0.4944±0.1188 0.0000 0.0000 0.0000 0.8399 Based on the five-fold cross-validation baseline data in Table 2, the model has the following problems: Falling into the trap of imbalanced samples (majority class bias): In linear and kernel function models (Logistic Regression and SVM), the average accuracy reached 83.74% and 83.99%, respectively. However, the average recall and F1-score for both were zero (0.0000). This proves that under the extreme imbalance of 46:1, traditional classifiers, in order to minimize global loss, adopted a completely passive strategy—blindly predicting all test sites as the majority class (not oxidizing 0). This phenomenon of high accuracy and zero recall irrefutably proves the complete failure of traditional models on the original data.

[0062] Local resistance and limitations of tree ensemble models: Thanks to their nonlinear spatial partitioning architecture, RandomForest and XGBoost initially demonstrated some resistance to imbalance, achieving average AUCs of 0.6760 and 0.6675, respectively. However, both exhibited extremely low recalls (0.10–0.18), and XGBoost had a larger standard deviation (0.0961), indicating that the system was highly unstable, with over 80% of real oxidation targets being severely underreported.

[0063] (5) Sample imbalance handling and optimization test based on SMOTE algorithm To address the curse of majority class bias exposed in benchmark testing, this invention introduces SMOTE (Synthetic Minority Over-sampling Technique). This algorithm analyzes the neighborhood distribution of rare oxidation samples in a multidimensional topological feature space, performs nonlinear interpolation between pairwise positive samples, and synthesizes reasonable virtual oxidation samples, thereby mathematically forcibly bringing the ratio of positive to negative samples back to balance. To avoid the most fatal beginner trap in machine learning—data leakage—when deploying SMOTE, this invention constructs an ImbPipeline mechanism based on the imbalanced-learn toolchain. This invention encapsulates SMOTE within the ImbPipeline pipeline. During stratified five-fold cross-validation, it ensures that SMOTE is only activated and generates virtual samples on the training set of the current fold, while keeping the original, true imbalanced distribution of the validation set (used as a blind test control) undisturbed. The process is as follows: Figure 2 As shown.

[0064] Based on the original 8-dimensional baseline features, after sample equalization processing using the SMOTE algorithm, the baseline performance of each machine learning model under five-fold cross-validation is shown in Table 3. After introducing the SMOTE algorithm, the major models successfully broke the bias of blindly betting on the majority class, which greatly improved the recognition sensitivity (recall) of rare oxidation sites.

[0065] Table 3 RandomForest 0.6688±0.0703 0.3385 0.3235 0.7782 0.6688±0.0703 XGBoost 0.6684±0.0387 0.3077 0.3074 0.7905 0.6684±0.0387 LogisticRegression 0.6436±0.0400 0.6769 0.3677 0.6232 0.6436±0.0400 SVM 0.6203±0.0681 0.5077 0.2961 0.6207 0.6203±0.0681 4. Feature Engineering (1) Feature Engineering 1: Sequence Spatial Microenvironment Perception and Scale Exhaustive Optimization Topological space sliding window radius ( Seeking the best To systematically explore the patterns of protein local linearity and spatial topological microenvironment influencing oxidative susceptibility, this invention designs and implements a scale-space exhaustive search scheme. The covalent modification of amino acid side chains depends not only on the static physicochemical properties of the site itself, but also on the synergistic influence of surrounding local conformational oscillations, spatial exposure, and macromolecular crowding pressure. Therefore, in the first stage of feature engineering, this invention performs a sliding window transformation on three dynamic spatial descriptors: SASA, Distance, and RMSF, adjusting the window radius parameter... The value range is set to 0 to 35 for model training. Since the total window length expression is... This search path essentially covers one amino acid ( ) to 71 amino acids ( A total of 36 different microenvironmental physical spans were examined to assess the influence of the spatial conformation of residues preceding and following the oxidative target. In the underlying program execution for feature matrix construction, the algorithm introduced a conditional isolation mechanism based on experimental group constraints. Physical isolation walls were established by grouping calculations according to concentration conditions, effectively preventing cross-group contamination of data by cutting off cross-interference of neighboring residue features under different macroscopic crowding pressures. Simultaneously, to address the matrix element null value calculation crashes that are easily caused by the loss of neighboring residues when the sliding window sweeps across the two ends (N-terminus and C-terminus) of the protein sequence, the program configured dynamic truncation adjustment (min) in the rolling operator. periods =1), ensuring that the algorithm can perform unbiased topological calculations based on the actual number of valid residues when the sequence edges are constrained, thus guaranteeing the alignment of the feature matrix. See the detailed process below. Figure 3 .

[0066] The iterative trend of the average AUC of various machine learning models in scale-space exhaustive search is as follows: Figure 4 As shown in the chart. Observing the chart, it can be seen that under the baseline state without introducing spatial microenvironment characteristics ( The predictive performance of all classifiers was poor (average AUC between 0.62 and 0.67), strongly demonstrating that isolated single-point physicochemical properties cannot provide sufficient spatial mechanism discrimination information for the model. With the sliding window radius... Starting from 1 and gradually increasing, the performance curves of XGBoost, RandomForest, and SVM all show a highly significant upward trend. Among them, XGBoost exhibits the optimal upper limit for feature capture, and... It precisely reached the peak of its global performance. The value reached 0.8289. This indicates that the window value centered on the target site and spanning 25 residues before and after the linear sequence (a total of 51 amino acids) is the optimal physical boundary for macromolecular crowding-induced protein oxidation sensing in the system of this invention.

[0067] XGboost model optimization Locking the optimal window radius Subsequently, this invention, using the best-performing XGBoost architecture as its core, implemented two rounds of targeted deep optimization experiments, aiming to transform the performance model into a highly robust and interpretable model with extremely high generalization ability. Addressing the performance fluctuations caused by the extreme scarcity (approximately 46:1) of positive protein oxidation samples, the research first adopted a dual robustness strategy in the first stage. This involved upgrading the validation foundation to five repeated Stratified 5-Fold Validation (SMOTE), executing a total of 25 independent training and blind testing cycles, aiming to completely eliminate statistical bias caused by single random sampling from the algorithm's underlying layers. At the sample governance level, the program introduces the SMOTE composite sampling algorithm through the ImbPipeline intelligent hybrid pipeline. Its core logic lies in using SMOTE to perform nonlinear interpolation in the minority class feature space to expand oxidation samples, and simultaneously using the TomekLinks operator to identify and remove overlapping noise sample pairs on the decision boundary, thereby improving AUC (e.g., ...). Figure 5 (As shown).

[0068] Subsequently, in the second phase, following Occam's Razor principle, the study performed deep weight analysis and simplification on the 51-dimensional original topological features. To suppress the model runaway phenomenon induced by high-dimensional features under small sample sizes, strict regularization constraints were imposed on the XGBoost classifier: setting max_depth=4 to force simplification of decision branches, using colsample_bytree=0.8 to achieve feature tolerance and decoupling, and simultaneously injecting L1 (reg_alpha=1) and L2 (reg_lambda=5) penalty terms to smooth the weights of leaf nodes and strictly prevent overfitting. Ultimately, the feature set was streamlined and locked into a Top 8 core descriptor system consisting of Polarity, SASA_Mean_w25 (regional average exposure), Distance_Mean_w25 (spatial depth), RMSF_Delta_w25 (regional flexible mutation), Hydrophobicity, pI (isoelectric point), RMSF, and SASA (see [link to relevant documentation]). Figure 6 ).

[0069] (2) Feature Engineering II: Cross-Dimensional Biochemical Interaction Descriptor (CBI) Protein oxidation is essentially an electrophilic attack of reactive oxygen species (ROS) on specific chemical groups. This process depends not only on the physicochemical properties of the residues themselves, but also on the nonlinear modulation of spatial compression and solvent repulsion effects caused by macromolecular congestion conditions. To capture this response mechanism between macroscopic environmental factors and microscopic structural traits, a set of high-order interaction operators based on product and ratio mappings was designed to simulate the amplification or inhibition of oxidation risk at locally exposed or flexible sites under different environmental pressures.

[0070] In the specific feature construction, the study first generated a dynamic hydrophobic exposure index (Dynamic_Hydro_Exposure) using the command df['Hydrophobicity']*df['SASA_Mean_w25']. Its significance lies in coupling static hydrophobicity with dynamic regional exposure, aiming to identify core engine sites that possess both strong hydrophobic repulsion and are in a highly exposed state. Subsequently, to establish a logical bridge between the macroscopic environment and microscopic conformational fluctuations, the program introduced a crowding-forced perturbation descriptor (Crowding_Perturbation), defined as follows: This feature, by weighting the congestant concentration as a local flexible mutation term, can accurately characterize the structural breakpoints that still exhibit anomalous conformational fluctuations under high-concentration congestion environments, i.e., the so-called structural tension concentration regions. Furthermore, regarding the high charge selectivity of oxidative attack, the study utilizes... A local electromagnetic ratio (pI) microenvironment descriptor was constructed to finely modulate the top-ranked feature polarity through the isoelectric point (pI), thereby outlining the integrated charge and hydrophilic network around the site that attracts reactive oxygen species.

[0071] This feature engineering approach helped the model complete the paradigm shift from spatial perspective to reaction dynamics awareness, significantly improving the AUC performance of XGBoost to 0.8927. However, the injection of high-dimensional interaction terms also introduced some information redundancy and noise interference. Therefore, the study again followed Occam's razor to implement deep feature pruning. By stripping away the mediocre original SASA, pI, and some low-gain interaction variants, the system finally extracted the ultimate Top 8 eigenvalues ​​consisting of SASA_Mean_w25, Hydrohobicity, Crowding_Perturbation, Polarity, RMSF_Delta_w25, Local_Electromagnetic, Distance_Mean_w25, and RMSF (see [link to relevant documentation]). Figure 9This refinement process, by eliminating the risk of overfitting, enabled the model to achieve a breakthrough in the final validation, with an AUC of 0.9112 ± 0.0281. Its standard deviation reached a minimum, and the accuracy remained stable at 0.9185. This signifies that the synergistic prediction system of biochemical interaction and topological microenvironment has achieved extremely high predictive ability under the extreme imbalance of 46:1 data (see...). Figure 10 ).

[0072] (3) Feature Engineering 3: Ablation Experiment (without local environmental window) ) This feature engineering design is for windows without local environments. Rigorous ablation experiments (Window-free). At the code implementation level, the feature matrix was forcibly compressed into a zero-window version model consisting only of single-point biochemical attributes, forcibly stripping features with global scale backgrounds such as Distance_Mean_w25 and SASA_Mean_w25. This ablation version model was also deployed within a robust validation protocol based on repeated cross-validation for baseline performance stress testing.

[0073] The experimental results strongly demonstrate the necessity of spatial window background information for predicting crowding-induced oxidation of macromolecules. For example... Figure 11 As shown, when dealing with 3108 extreme imbalanced samples (oxidized:unoxidized = 65:3043), the ablation model maintained a high global Accurate (0.92) and a good Recall. This demonstrates that even in the absence of scale context, single-point attributes and CBI descriptors still possess a certain discriminative signal. However, the AUC drops to [missing value]. .

[0074] Furthermore, the Top9 eigenvalue ranking of ablation experimental models decoded the reasons for model failure from a mechanistic perspective. When lacking... In the context of window-scale conditions, the model is forced to rely excessively on static, non-smooth physicochemical parameters again. Although the Self_Crowding_Perturbation descriptor, which is a result of the coupling between the macroscopic crowded environment and single-point flexible nonlinearity, still holds the top importance (0.2387), confirming the robustness of the cross-dimensional interaction mechanism, the importance of the underlying primitive descriptors such as Hydrophobicity (0.1363), SASA (single-point exposure, 0.1151), and Polarity (0.1147) has significantly rebounded (see ranking list). This feature confirms the scientific hypothesis: once the optimal sliding window radius (w) descriptor proposed in this invention is removed, the model will lose its ability to sense the local background and be forced to revert to the fitting trap for single-point residue characteristics, thus leading to the failure to accurately distinguish between oxidized and non-oxidized signals in the context of complex macromolecular crowding dynamics. This ablation experiment strongly supports that window optimization (w=25) is not only effective in terms of data but also scientifically reasonable in terms of physical meaning.

[0075] (4) Feature Engineering 4: Microenvironment Heterogeneity / Volatility Index The core of this stage is to move away from dependence on local average states and instead explore how drastic fluctuations in physical properties induce chemical bond breaking. The principle lies in identifying the mechanical tension generated at the rigid-flexible interface in the protein sequence and the surface roughness at folding cracks to pinpoint the breeding ground for oxidation events. In implementation, this invention adheres to the principle of isolated testing, retaining only the basic eigenvalues, conditions, and newly generated volatility indices for supervised training.

[0076] Experimental results ( Figure 12 The results show that the independent feature set achieved an AUC of 0.8524 ± 0.0350 in the five-fold cross-validation. Although this is significantly better than the original feature baseline, it failed to break through the AUC value of 3.4.2. This negative validation result strongly defines the technical boundary and proves that although local physical ruggedness has an advantage, it is not the absolute main factor driving oxidation in crowded environments.

[0077] III. Experimental Results 1. Model performance evolution and algorithm superiority verification To verify the effectiveness of the constructed multi-scale spatial topological feature system, this invention records the driving trajectories of different feature engineering stages on the prediction performance of the XGBoost model. Figure 13In the control group containing only traditional physicochemical descriptors, the model's average AUC remained at a low level (~0.67), indicating that single sequence information cannot effectively capture complex oxidation patterns. When spatially local topological descriptors were introduced (feature engineering 1), the model performance improved dramatically; and when cross-dimensional biochemical interaction descriptors were further injected (feature engineering 2), the model performance reached global optimum, with an average AUC as high as 0.9112. This strongly demonstrates that protein oxidation is not an isolated amino acid mutation, but is determined by the topological tension of its surrounding microenvironment and the biochemical network.

[0078] After establishing the optimal feature space, this invention conducts a rigorous comprehensive performance evaluation of the champion model. Figure 14 The ROC curve generated through multiple stratified cross-validations is extremely smooth and steep, and its very narrow error band (light red shaded area) indicates that the model exhibits exceptional robustness and generalization ability when dealing with extreme positive-negative sample imbalances of up to 46:1. This model successfully overcomes the limitations of traditional structure prediction and possesses the potential for application in real-world experimental scenarios.

[0079] 2. Core driver feature decoding based on SHAP To break the black-box nature of machine learning and explore the underlying driving mechanism of protein oxidation, this invention introduces the game theory-based SHAP (Shapley Additive ex Planations) explanatory framework. Figure 15 The global feature importance decoding results show that Hydrodrophobicity (native hydrophobicity, 1.1398) and SASA_Mean_w25 (average solvent-accessible surface area of ​​the microenvironment, 1.0340) are absolutely dominant. This result is highly consistent with classical physicochemical theory: oxidative attacks tend to occur first in microenvironments with high exposure and rich in hydrophobic residues.

[0080] Figure 16 The core driving feature attribution beehive plot further reveals the targeted regulatory mechanism of each feature on the oxidation probability. In the beehive plot, the red scatter points (high feature values) are highly clustered in the SHAP>0 region of Hydrodrophobicity and SASA_Mean_w25 (right side of the coordinate axis), intuitively confirming that the increase in hydrophobicity and spatial exposure directly promotes the oxidation risk. It is worth noting that Crowding_Perturbation, as the third key feature, also significantly drives the increase in oxidation probability in its high value range. This profoundly indicates that in the crowded environment of macromolecules, the abnormal fluctuations and tension release of local structures induced by spatial compression are the core causes of breaking down the protein's antioxidant defense and generating oxidation hotspots.

[0081] 3. Macroscopic regulation of oxidation probability by macromolecular crowding concentration To explore the macroscopic regulatory mechanisms of the crowding effect, this invention utilizes a global prediction model to output the evolutionary distribution of protein oxidation probabilities under different macromolecular crowding environment parameters (0, 20, 40, 60, 80, 100). Figure 17 The morphological evolution and average trend line (red line) of the violin plot reveal a remarkably significant nonlinear response mechanism: when the crowding concentration is at a low to medium level (D20-D40), the overall oxidation risk decreases instead of increasing, the bottom of the violin plot widens, and the average probability drops to its lowest point. This suggests that a moderately crowded environment may encourage proteins to form a denser, more stable self-protective conformation, masking potential sensitive sites. However, when the concentration exceeds the critical value and reaches extreme crowding (D80-D100), the top region of the violin plot expands sharply, indicating a large emergence of extremely high-risk amino acid groups, and the structural defenses suffer a global collapse.

[0082] Figure 18 The absolute number of sites confirms this macroscopic pattern. Under the safe haven state of D40, the number of high-risk oxidation sites (Predicted=1) drops to a minimum of 97; however, under the extreme stress environment of D100, the number of high-risk sites surges to 130. This order-of-magnitude jump quantitatively characterizes the destructive power of an extremely crowded microenvironment on the three-dimensional topological network of proteins.

[0083] 4. Sequence heatmap and 3D topological fingerprint mapping To achieve high-precision spatiotemporal localization of oxidation risk, this invention constructs a joint mapping system of full-sequence 2D heatmap and 3D topological fingerprint. Figure 19 The full-sequence oxidation risk evolution heatmap clearly shows that crowding-induced oxidation does not occur uniformly along the peptide chain, but rather exhibits strong region specificity. As the crowding concentration gradient increases (vertically downwards), specific sequence regions that were previously within the safe threshold in dilute solutions (such as continuous blue bands) undergo a dramatic red shift after crossing a specific concentration threshold. These vertically penetrating red high-risk bands constitute the specific oxidation-sensitive fingerprint of oat globulin under crowding conditions.

[0084] After converting the above multidimensional prediction probabilities into a three-dimensional temperature factor (B-factor), Figure 20This digital twin visualization vividly presents the 3D structural evolution trajectory of oxidation risk. In uncrowded (D0) and moderately crowded (D20-D40) environments, the red regions representing high risk are only scattered sporadically in the peripheral loop regions. However, with increasing environmental pressure (D80-D100), these high-risk red hotspots rapidly spread and aggregate over large areas in specific hydrophobic regions and core structural domains. This digital twin visualization perfectly validates the mechanism inferences regarding crowding-induced perturbations from previous SHAP analyses, reproducing from a microscopic perspective of atomic spatial topology how extreme microenvironments strip away the natural folds of proteins, thereby exposing vulnerable oxidation targets.

[0085] This invention successfully constructs a digital twin model for predicting protein oxidation risk based on multi-scale spatial topological features and biochemical interaction descriptors. Research shows that relying solely on traditional physicochemical properties cannot accurately define oxidation hotspots; however, by deeply mining the microenvironmental properties of amino acid residues (such as SASA_Mean_w25 and Crowding_Perturbation) through feature engineering, and combining the XGBoost algorithm with the SMOTE sample balancing strategy, the model can accurately capture oxidation targets in an extremely imbalanced sample space with extremely high precision (average AUC = 0.9112). This model completely overcomes the lag of traditional experimental methods, providing a highly robust and generalizable computational paradigm for predicting the oxidation of food proteins.

[0086] Based on the SHAP interpretation framework and global prediction model, this invention, for the first time, quantitatively reveals the nonlinear dynamic response mechanism of oat globulin in different macromolecular crowding microenvironments at the microscale. The study found that oxidation does not increase linearly with increasing crowding concentration: in moderately crowded environments (D40-D60), the protein forms a dense shielding effect through structural remodeling, significantly reducing the overall oxidation probability; however, when the crowding pressure exceeds the critical threshold (D80-D100), extreme forced perturbation leads to an abnormal release of local topological tension, causing extreme exposure of the hydrophobic regions at the interface, thereby completely breaking down the antioxidant defenses and resulting in an exponential surge in the number of high-risk oxidation sites.

[0087] This invention innovatively maps a one-dimensional prediction probability matrix to a three-dimensional temperature factor (B-factor), constructing a visualized 3D oxidation topological fingerprint spectrum. The sequence heatmap and the 3D fingerprint spectrum perfectly corroborate the SHAP mechanism analysis, intuitively presenting the spatial evolution trajectory of high-risk oxidation regions under extreme crowding and pressure. This visualization tool not only imbues the extremely complex machine learning output with clear physicochemical intuition, but also provides new ideas for subsequent targeted design of antioxidant strategies and precise modification of vulnerable sites.

[0088] In the era of big health and precision nutrition, the dramatic changes in the microenvironment during food processing and storage pose significant challenges to protein quality. This invention breaks away from the empirical framework of traditional food chemistry, which relies on macroscopic characterization and post-hoc testing. It innovatively integrates molecular dynamics simulation, multidimensional spatial topological feature engineering, and interpretable artificial intelligence. This not only reconstructs our understanding of the oxidative evolution of proteins under macromolecular crowding at the underlying mechanism level, but also provides highly forward-looking theoretical guidance and digital twin underlying technology support for future food processing systems (such as the formulation design of high-concentration mixed systems and the intelligent control of storage conditions).

[0089] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0090] The above embodiments are merely illustrative examples and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for predicting protein oxidation risk based on multi-scale spatial topology and biochemical interactions, characterized in that, Includes the following steps: (1) By molecular dynamics simulation, the molecular dynamics trajectory of the target protein under different macromolecular crowding environment parameters is obtained, and the original physical characteristics of each amino acid residue in the whole sequence are extracted from it. The original physical characteristics include solvent-accessible surface area, centroid spatial distance and root mean square fluctuation. (2) Based on the original physical characteristics, for each target amino acid residue in the sequence, a sliding window is constructed with it as the center, and the local environment smoothing descriptor and the environment relative mutation descriptor within the window are calculated to obtain the spatial local topology descriptor; the local environment smoothing descriptor is obtained by calculating the average value of the original physical characteristics of the effective residues within the window; the environment relative mutation descriptor is the original physical characteristic value of the target residue minus the local environment smoothing descriptor. (3) Multiply the environmental parameters representing macroscopic crowding concentration with the environmental relative mutation descriptor and the local environmental smoothing descriptor in the spatial local topology descriptor respectively to construct a cross-dimensional biochemical interaction descriptor; the cross-dimensional biochemical interaction descriptor includes: a crowding-forced disturbance descriptor, which is obtained by multiplying the macroscopic environmental parameters with the environmental relative mutation descriptor of root mean square fluctuation; and an environment-mediated exposure descriptor, which is obtained by multiplying the macroscopic environmental parameters with the local environmental smoothing descriptor of solvent-accessible surface area; (4) Integrate the original physical features, spatial local topological descriptors and cross-dimensional biochemical interaction descriptors to form a multi-scale feature set, and use protein site data with known oxidation states as labels to train the machine learning model. The protein site tags with known oxidation states were obtained by liquid chromatography-tandem mass spectrometry. Positive samples were amino acid residues that had been confirmed to have undergone oxidation modification, and negative samples were residues that had not undergone oxidation modification. The machine learning model is the XGBoost model, and the training set is balanced with positive and negative samples using the SMOTE algorithm or SMOTETomek composite sampling algorithm before training. Under the hierarchical cross-validation framework, the intelligent hybrid pipeline ensures that the sampling only applies to the training set. The multi-scale feature set was filtered by feature importance, and the core features retained include: local environmental smoothing descriptors of polarity and solvent-accessible surface area, local environmental smoothing descriptors of centroid spatial distance, environmental relative abrupt change descriptors of root mean square fluctuations, hydrophobicity, isoelectric point, original features of root mean square fluctuations, and original features of solvent-accessible surface area. After training, the SHAP interpretation framework is used to perform global feature importance analysis on the trained machine learning model, and output the ranking of the contribution of each feature to oxidation risk and the direction of directional regulation. (5) Input the multi-scale feature set of the protein site to be predicted into the trained machine learning model and output the oxidation risk prediction result of the site. The output oxidation risk prediction result includes: the oxidation probability value of each site, the full sequence oxidation risk heatmap under different crowding concentrations, and the temperature factor distribution map that maps the oxidation probability to the three-dimensional structure of the protein.

2. The protein oxidation risk prediction method based on multi-scale spatial topology and biochemical interactions according to claim 1, characterized in that, In step (2), when calculating the sliding window descriptor, dynamic truncation adjustment is used for the N-terminal and C-terminal residues of the protein sequence, and the calculation is based on the actual number of effective residues in the window.

3. The protein oxidation risk prediction method based on multi-scale spatial topology and biochemical interactions according to claim 1, characterized in that, The training process in step (4) adopts the repeated stratified five-fold cross-validation protocol, which is repeated 5 times, for a total of 25 independent training and validation cycles.

Images