Traditional Chinese medicine action mechanism analysis method based on network pharmacology

Abstract

The invention relates to the technical field of traditional Chinese medicine action mechanism analysis methods based on network pharmacology, discloses a traditional Chinese medicine action mechanism analysis method based on network pharmacology, and aims to solve the technical problems of low action mechanism analysis precision and difficulty in prediction and verification in the prior art. The method comprises the following steps: screening traditional Chinese medicine bioactive components; screening disease related targets in combination with disease tissue specificity; constructing a weighted bipartite graph recognition synergistic effect module; dynamically simulating path disturbance; and the core action mechanism is verified in combination with experimental omics data. By adopting the technical scheme, the precision and reliability of traditional Chinese medicine action mechanism analysis can be remarkably improved, the synergistic effect and the time sequence effect are comprehensively revealed, and a scientific basis is provided for modern research of traditional Chinese medicine.

Description

Technical Field

[0001] This invention belongs to the field of biomedical technology, and more specifically, relates to a method for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology. Background Technology

[0002] The field of biomedicine integrates multiple disciplines such as modern life sciences, pharmacology, bioinformatics, and computer science. The core goal of this field is to reveal the molecular basis of life activities, elucidate the mechanisms of disease occurrence and development, and develop new diagnostic and treatment methods accordingly. As an important part of this field, traditional Chinese medicine provides a valuable resource for new drug development with its long history and rich clinical practice. With the development of genomics and proteomics technologies, systematically analyzing the mechanism of action of traditional Chinese medicine at the molecular level has become a research hotspot in this field, promoting the transformation of traditional medicine towards precision and modernization.

[0003] Among them, network pharmacology analysis refers to a research strategy based on systems biology theory, which systematically predicts and elucidates the mechanism of drug action by constructing a multi-dimensional interaction network of drugs, targets, and diseases. This method aims to solve the complexity of the synergistic effects of multiple components, targets, and pathways of traditional Chinese medicine that traditional pharmacology research cannot handle. Its basic process includes screening active ingredients of traditional Chinese medicine, predicting potential targets, and constructing and analyzing the "component-target-disease" interaction network, thereby revealing the molecular mechanism of drug intervention in biological networks at the overall level. This method provides a new perspective and technical means for understanding the overall regulatory effects of traditional Chinese medicine.

[0004] Current technologies for network pharmacology analysis largely rely on manual retrieval and integration of data from multiple public databases. This process is cumbersome, has low data standardization, and is prone to introducing biases. Network construction commonly suffers from node redundancy and information overload; while the generated networks are large, crucial information is buried, making it difficult to effectively identify core functional modules. Analytical methods often focus on static calculations of network topology parameters, lacking in-depth analysis of dynamic network responses and pathway interactions, leading to inaccurate identification of major drug pathways and key targets. Furthermore, analytical results often consist mainly of listing massive numbers of targets and pathways, lacking effective extraction and visualization of core mechanisms of action. This is particularly pronounced in the analysis of complex systems like traditional Chinese medicine compound formulas, directly impacting the accuracy of the results and the depth of understanding of the core drug mechanisms. Summary of the Invention

[0005] The purpose of this invention is to provide a method for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology, aiming to solve the technical problems of low accuracy in mechanism analysis and poor consistency between prediction results and biological verification caused by the proliferation of targets, lack of biological spatiotemporal specificity, and unclear synergistic mechanism in the prior art.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a method for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology, comprising the following steps: S1: Receive traditional Chinese medicine compound data, which includes the names of all medicinal materials constituting the compound; parse the traditional Chinese medicine compound data, query a pre-set traditional Chinese medicine component database, obtain all compound components associated with each medicinal material, and extract the molecular structure identifier of each compound to generate an initial compound set; based on the initial compound set, call a pre-trained graph convolutional neural network pharmacokinetic property prediction model for each compound to calculate the oral bioavailability value and drug-likeness value of each compound, and perform dual filtering based on preset oral bioavailability thresholds and drug-likeness thresholds to screen out compounds whose oral bioavailability and drug-likeness are both higher than the corresponding thresholds, forming a bioactive component profile.

[0007] S2: Receive the disease name associated with the indication of the traditional Chinese medicine compound; based on the disease name, query a preset human tissue gene expression map database, identify the lesion tissue type associated with the disease, and extract all genes whose gene expression level is higher than the preset expression level percentile in the lesion tissue to form a tissue-specific high-expression gene set; simultaneously, query a target-compound interaction database to obtain all potential targets corresponding to each compound in the bioactive ingredient spectrum to form a compound potential total target set; perform a set intersection operation on the tissue-specific high-expression gene set and the compound potential total target set to generate a highly correlated disease-drug target set that simultaneously meets the dual conditions of tissue-specific high expression and compound ingredient targeting.

[0008] S3: Based on the bioactive component spectrum and the disease-drug target set, a weighted bipartite graph network model is constructed, where the node set on one side of the bipartite graph consists of each compound in the bioactive component spectrum, and the node set on the other side consists of each target in the disease-drug target set. If a compound interacts with a target, a connection edge is established between the two nodes, and a corresponding weight value is assigned to the connection edge based on the binding strength score recorded in the target-compound interaction database. For the constructed weighted bipartite graph network model, a weighted random block model algorithm is applied to perform community detection analysis to identify tightly connected subnetwork structures composed of multiple compounds acting on a functionally related target cluster, where the total internal connection weight is much higher than the random expectation, and each such subnetwork structure is determined as a synergistic module.

[0009] S4: For each identified synergistic module, map its target set to a pre-set biological pathway database to determine the core signaling pathway mainly regulated by the module; extract the topological structure of the core signaling pathway and construct a pathway dynamic system model described by a set of ordinary differential equations, where the state variables (such as concentration or activity) of each pathway component change over time by its upstream activation signal, inhibition signal, and its own degradation rate; input the effects (activation or inhibition) of the compounds in the synergistic module on the target as initial perturbation conditions into the set of ordinary differential equations model, and solve numerically using the fourth-order Runge-Kutta method to simulate the dynamic changes of each key downstream molecule in the core signaling pathway within a preset time series under drug perturbation, generating a quantified predictive pathway perturbation map containing data on the changes in concentration or activity of each molecule over time.

[0010] S5: Receive experimental omics data obtained after intervention of the traditional Chinese medicine compound on a biological model (cell or animal model) related to the disease. The experimental omics data is a list of differentially expressed genes or proteins, including the fold change in expression and statistical significance value of each gene or protein. Sort the experimental omics data list according to the logarithm of the fold change in expression. Take the set of perturbed downstream molecules in the predictive pathway perturbation map generated in step S4 as the gene set to be verified. Apply a pre-sorted gene set enrichment analysis algorithm to calculate the enrichment score of the gene set to be verified in the sorted experimental omics data list, and generate a normalized enrichment score and false discovery rate. Calculate a comprehensive verification score based on the normalized enrichment score and false discovery rate. When the comprehensive verification score is higher than a preset compliance threshold, confirm that the synergistic module and its regulated core signaling pathway corresponding to the score are the core mechanism of action of the traditional Chinese medicine compound verified by experimental data, and output an analysis report containing detailed information on the core mechanism of action.

[0011] As a further embodiment of the present invention, the process of generating the bioactive component spectrum in step S1 is specifically as follows: S111: Receive traditional Chinese medicine compound data composed of medicinal material names, access a locally deployed structured query language database containing hierarchical relationships between traditional Chinese medicine prescriptions, medicinal materials, and chemical components, perform a query operation, extract the names, CAS numbers, and two-dimensional molecular structure SMILES strings of all chemical components corresponding to each medicinal material, merge and remove duplicates to form an initial compound set; S112: For each compound in the initial compound set, input its SMILES string into a molecular graph conversion module to generate the graph representation data structure of the compound; convert the graph... The data structure is used as input to a pre-trained graph convolutional neural network model, which is trained on a large-scale ADMET dataset of known drugs. Its output layer consists of two independent neurons that output oral bioavailability prediction probability values ​​and drug-likeness prediction probability values ​​between 0 and 1, respectively. S113: Set the oral bioavailability threshold to 0.3 and the drug-likeness threshold to 0.18. Iterate through the prediction probability values ​​of all compounds, select compounds with oral bioavailability prediction probability values ​​greater than 0.3 and drug-likeness prediction probability values ​​greater than 0.18, and store their related SMILES strings, CAS numbers, and medicinal material information into a new data table to form the bioactive component profile.

[0012] As a further embodiment of the present invention, the process of generating the disease-drug target set in step S2 is as follows: S211: Receive the disease name represented by the International Classification of Diseases (ICD) code and access a locally mirrored Human Protein Atlas PostgreSQL database; based on the tissue type associated with the disease name, perform a database query to select genes whose transcript expression per million base pairs (TPM) value is among the top 25% of all genes in that tissue, and extract their corresponding UniProt protein IDs to form a tissue-specific high-expression gene set; S212: Access a local STITCH database containing compound-protein interaction information, use the CAS numbers of all compounds in the bioactive component spectrum as an index, batch query and extract the UniProt IDs of all protein targets with interaction records, and remove duplicates to form a compound potential total target set; S213: Compare the UniProt ID list in the tissue-specific high-expression gene set with the UniProt ID list in the compound potential total target set using an efficient algorithm for finding the intersection of hash sets, and output the UniProt IDs common to both sets. IDs constitute the disease-drug target set.

[0013] As a further embodiment of the present invention, the identification process of the synergistic module in step S3 is specifically as follows: S311: Create a graph data structure whose node attribute table includes node ID, node type (compound or target), and node name, and whose edge attribute table includes source node ID, target node ID, and weight value; fill the node attribute table with compounds in the bioactive component spectrum as one type of node and targets in the disease-drug target set as another type of node; S312: Traverse all known interaction pairs between compounds and targets, create connection edges between corresponding nodes, and extract their values ​​from the STITCH database. The combination score (range 0-1000) is used as the weight value of the edge and filled into the edge attribute table to complete the construction of the weighted bipartite graph network; S313: The weighted random block model algorithm is applied to iteratively partition the graph network. The objective function of the algorithm is to maximize the posterior probability of the graph. It is optimized by the Markov chain Monte Carlo method until the community partitioning results converge; The algorithm finally outputs several communities, each of which is a sub-network structure containing several compound nodes and target nodes. The internal weighted connection density between these nodes is significantly higher than the connection density between communities. Each such community is defined as a cooperative module.

[0014] As a further embodiment of the present invention, the process of generating the predictive path perturbation map in step S4 is specifically as follows: S411: Extract all target points within a synergistic module, batch query the API of the KEGG database to obtain the IDs and names of all signal paths in which these target points participate; perform statistical sorting based on the frequency of the returned path IDs, and select the path with the highest frequency as the core signal path of the module; S412: Download the KGML (KEGGMarkup) of the core signal path. The language file is parsed to obtain all molecular nodes in the pathway and their interrelationships (such as activation, inhibition, phosphorylation), and a directed graph is constructed to represent the pathway topology. This topology is transformed into a system of ordinary differential equations. For example, for the reaction A+B→C, the equation is d[C] / dt=k_f*[A]*[B]-k_r*[C], where k is the reaction rate constant. S413: The effect of the compounds in the synergistic module on the target site (activation or inhibition) is set as the initial condition of the equation system. Specifically, the initial concentration of the state variable or the reaction rate constant corresponding to the target site is quantitatively increased or decreased. The odeint function in the Python SciPy library deployed on the computing server is used. This function implements the fourth-order Runge-Kutta numerical integration algorithm to solve the entire equation system in the virtual time range from 0 to t_end, with a step size of dt, to obtain an N×M dimensional matrix of the change of the concentration of all molecules in the pathway over time, where N is the number of molecules and M is the number of time points. This matrix is ​​the predictive pathway perturbation map.

[0015] As a further embodiment of the present invention, the verification and confirmation process of the core mechanism of action in step S5 is specifically as follows: S511: Receive a text file containing three columns (gene / protein ID, log2FoldChange, p-value) as experimental omics data; firstly, filter out statistically insignificant data based on p-value < 0.05, and then sort the remaining genes / proteins in descending order according to log2FoldChange values ​​to generate a pre-sorted list of experimental genes; S512: Extract a set of downstream perturbed molecules corresponding to a predictive pathway perturbation map, and calculate the perturbation effect of each molecule in the model... The concentration change rate between the simulated end point and the initial point is obtained by screening molecules whose absolute change rate is greater than a preset perturbation threshold; S513: Call a locally implemented GSEA program, input the pre-sorted experimental gene list and the downstream perturbed molecule set; the program calculates an enrichment score (ES), by traversing from the top to the bottom of the list, adding an increment related to the expression change when encountering a gene in the set, and deducting a fixed decrement when encountering a gene outside the set; S514: By performing 1000 random permutations on the gene list and repeatedly calculating ES, a null distribution of ES is constructed, thereby calculating the standardized enrichment score (NES) and false discovery rate (FDR q-value); calculate the comprehensive validation score V=NES*(-10*log10(FDR)); set the compliance threshold to 2.0, when V>2.0, mark the corresponding synergistic module and its regulatory pathway as "validated", and integrate its topology map, pathway perturbation simulation curve, GSEA enrichment map and all related scores into a PDF document as the final analysis report output.

[0016] This invention also provides a system for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology, used to perform the above method, the system comprising: The bioactive ingredient screening module receives data from traditional Chinese medicine compound prescriptions, connects to a database of traditional Chinese medicine ingredients and a pre-trained graph convolutional neural network ADMET prediction model, and calculates and filters all compounds in the compound prescriptions by setting dual thresholds for oral-oral bioavailability and drug-likeness, generating a data spectrum containing only potential bioactive ingredients. The tissue-specific target enrichment module receives disease names, connects to the human tissue gene expression map database and the target-compound interaction database, and generates a set of drugs and targets that are highly correlated with the disease phenotype by querying highly expressed genes in diseased tissues and performing intersection operations with the potential target set of compound components. The synergistic effect module discovery module is used to construct a weighted compound-target bipartite graph network based on the bioactive component spectrum and drug-target set, and run a weighted random block model algorithm to identify the statistically significant tight-connected sub-network structure formed by multiple components synergistically acting on multiple related targets, i.e., the synergistic effect module. The pathway dynamic perturbation simulation module is used to map the identified synergistic modules to a biological pathway database, construct a dynamic model of the ordinary differential equation of the core regulatory pathway, and use the targeting effect of the compounds in the module as the initial perturbation. The module simulates the dynamic response process of the pathway under drug intervention through numerical solution, and generates a quantitative predictive pathway perturbation map. The experimental evidence-based verification module receives externally input experimental omics data, applies a pre-sorted gene set enrichment analysis algorithm, statistically compares the predicted pathway perturbation map with the actual biological experimental results, calculates a comprehensive verification score, and confirms or refutes the predicted synergistic mechanism based on the score, and finally outputs a verified core mechanism report.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention significantly reduces the chemical space of initial analysis by introducing pharmacokinetic property prediction as a pre-filter, eliminating a large number of ineffective components that cannot reach effective concentrations in vivo, thus improving the biological relevance of the analysis from the source. Furthermore, this invention uses disease tissue specificity as a key constraint, focusing target screening on genes highly expressed at lesion sites, avoiding ineffective searches on systemically abundant targets, resulting in stronger disease correlation of the target set and significantly improving the accuracy of mechanism prediction. Crucially, this invention abandons the traditional analysis mode of simple superposition of single targets or pathways, innovatively employing a weighted bipartite graph community detection algorithm to actively identify "synergistic modules" of multiple compounds and multiple targets. This method can reveal the combinatorial and systemic action mode of traditional Chinese medicine compound as a whole from complex interaction networks, more closely reflecting the actual effects of traditional Chinese medicine. In addition, this invention constructs a dynamic model of ordinary differential equations to simulate pathway perturbations, elevating static network topology analysis to the prediction of dynamic system behavior, revealing the temporal effects and downstream influences of drug action. Ultimately, this invention establishes a closed loop of mutual verification between computational prediction and experimental data. By using gene set enrichment analysis, the computational results are quantitatively compared and scored with real omics data, ensuring that the final output mechanism of action is not only theoretically reasonable but also supported at the biological experimental level. This provides a highly reliable and credible scientific basis for the modernization research of traditional Chinese medicine and the discovery of new drugs. Attached Figure Description

[0018] Figure 1This is a schematic diagram of the overall process of a method for analyzing the mechanism of action of traditional Chinese medicine provided in an embodiment of the present invention; Figure 2 A functional module block diagram of a traditional Chinese medicine mechanism analysis system provided in an embodiment of the present invention; Figure 3 This is a flowchart illustrating the discovery process of the synergistic effect module in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating the principle of the path dynamic disturbance simulation process in an embodiment of the present invention; Figure 5 This is a schematic diagram illustrating an application scenario combining computational prediction and experimental data verification, provided by an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In the description of this invention, it should be understood that the terms "upper," "lower," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, in the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0020] Please see Figure 1 The present invention provides a method for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology, the overall flowchart of which is shown in the figure below. Figure 1 As shown. The method includes the following steps: Step S1: Receive and parse the data of traditional Chinese medicine compound prescriptions, and screen the compounds in the compound prescriptions by calling the pre-trained pharmacokinetic property prediction model to generate a bioactive component profile.

[0021] Step S2: Receive the disease name, combine human tissue gene expression profiles and compound target data, perform tissue-specific target enrichment, and generate a drug-target set highly related to the disease.

[0022] Step S3: Based on the bioactive component profile and disease-drug target set, a weighted bipartite graph network is constructed, and a community detection algorithm is applied to identify synergistic modules. A flowchart can be found below. Figure 3 .

[0023] Step S4: For the identified synergistic modules, construct a dynamic system model of their core regulatory pathways, simulate the dynamic response of the pathways under drug perturbation, and generate a predictive pathway perturbation map. A schematic diagram of this principle can be found in [reference needed]. Figure 4 .

[0024] Step S5 involves receiving external experimental omics data, statistically validating the predicted results through pre-sorted gene set enrichment analysis, confirming the core mechanism of action, and outputting an analysis report. An application scenario diagram can be found in the provided text. Figure 5 .

[0025] In a specific embodiment, the specific execution process of steps S1 and S2 is described in detail.

[0026] Step S1, the generation of the bioactive component profile, can be broken down into the following sub-steps: Step S101: Receive traditional Chinese medicine (TCM) compound data, query a pre-set TCM component database, and obtain an initial set of compounds. Specifically, the system receives TCM compound data defined in a structured format (e.g., JSON or XML). In one instance, the input data is the "Siwu Tang" compound for "liver stagnation and blood deficiency" syndrome, with the data format {"formula_name":"Si-Wu-Tang","herbs":["Angelica sinensis","Radix Paeoniae Alba","RhizomaChuanxiong","Radix Rehmanniae Preparata"]}. After receiving this data, the system accesses a locally deployed Structured Query Language (SQL) database containing a hierarchical relationship between TCM formulas, medicinal materials, and chemical components. The database contains at least four core tables: the `formulas` table (stores prescription information), the `herbs` table (stores medicinal material information), and the `compounds` table (stores chemical component information, including CAS numbers and SMILES molecular structure strings), as well as the `formula_herb_map` and `herb_compound_map` mapping tables that link them together. The system executes a cascading query operation: first, it finds the corresponding `formula_id` in the `formulas` table based on "Si-Wu-Tang"; then, it finds all associated `herb_id`s in the `formula_herb_map` table; and finally, it extracts records of all chemical components corresponding to each `herb_id` from the `herb_compound_map` table, including the compound name, CAS registration number, and SMILES strings representing the two-dimensional molecular structure.

[0027] For example, for the medicinal material "Angelica sinensis" (Chinese angelica), the compounds that the query may return include: Ferulic acid (CAS: 1135-24-6, SMILES: COC1=C(O)C=C( / C=C / C(=O)O)C=C1), Z-Ligustilide (CAS: 4431-01-0, SMILES: C1C=C(CCC)C2=C(C1)C(=O)OC2), etc. For the medicinal material "Radix Paeoniae Alba" (white peony root), Paeoniflorin (CAS: 23180-57-6, SMILES: O=C(OC[C@H]1OC@@HC@HC@H[C@@H]1O)C4=CC=CC=C4) may be returned. The records of all the compounds queried from all four medicinal materials are merged, and duplicate entries are removed by performing uniqueness checks on the CAS numbers or SMILES strings, finally generating an initial compound set containing hundreds or even thousands of compounds. The data structure of this set is a list, where each element is a dictionary object containing {CompoundName, CAS_Number, SMILES_String, Source_Herb}.

[0028] In step S102, for each compound in the initial compound set, a pre-trained graph convolutional neural network (GCN) model is called to calculate the predicted values of its oral bioavailability (OB) and drug-likeness (DL). Specifically, the system traverses the initial compound set. For each compound, first, a molecular graph conversion module processes its SMILES string. This module parses the SMILES string into a graph representation data structure. In this structure, each atom is represented as a node, and the chemical bond is represented as an edge connecting the nodes. The node features are encoded as a high-dimensional vector, which specifically includes: the type of the atom (such as C, N, O, S, one-hot encoded), the atomic number, the node degree (the number of connected bonds), the formal charge, the hybridization type (sp, sp2, sp3, one-hot encoded), and whether it is in an aromatic ring (boolean value). The features of the edge are also encoded, including the type of the bond (single bond, double bond, triple bond, aromatic bond) and whether it is in a ring.

[0029] Subsequently, this diagram illustrates how the data structure is fed into a pre-trained graph convolutional neural network model. In one specific embodiment, this model is based on a Graph Isomorphism Network (GIN) architecture, containing five GIN convolutional layers, each followed by a batch normalization layer and a ReLU activation function. After stacking the convolutional layers, a global average pooling layer aggregates the node features into a representation vector of the entire molecular graph. This vector then passes through two independent fully connected layers, each followed by a sigmoid activation function, outputting a probability value between 0 and 1, namely the oral bioavailability (OB) prediction and the drug-likeness (DL) prediction. It should be noted that the pre-training process of this model was performed on a public dataset (e.g., ADMETLab2.0) containing over 200,000 molecules with known ADMET (absorption, distribution, metabolism, excretion, toxicity) properties. Training employed the Adam optimizer with a learning rate of 0.0001 and a binary cross-entropy loss function.

[0030] As a simulation example, for ferulic acid in the initial compound set, its SMILES string is converted into a graph structure and input into the model. The model output might be {"OB_score":0.453,"DL_score":0.218}. For another compound with a larger molecular weight and more complex structure, such as a certain natural saponin, the model output might be {"OB_score":0.076,"DL_score":0.152}.

[0031] In step S103, the predicted values ​​of all compounds are double-filtered according to preset thresholds to select qualified compounds and form the final bioactive component profile. Specifically, the system sets the oral bioavailability threshold to 0.3 and the drug-likeness threshold to 0.18. These two thresholds are based on statistical analysis of approved oral drugs in the DrugBank database, which shows that approximately 85% of successful oral drugs meet these two criteria. This setting aims to maximize the retention of potential active components while filtering out a large number of ineffective components with poor pharmacokinetic properties. The system iterates through all compounds and their OB and DL predicted values ​​calculated in step S102. A compound is determined to have potential oral activity if and only if its OB_score > 0.3 and DL_score > 0.18.

[0032] Continuing with the previous simulation example, the predicted value for ferulic acid is {"OB_score":0.453,"DL_score":0.218}. Since 0.453 > 0.3 and 0.218 > 0.18, this compound passes the screening. The predicted value for the natural saponin is {"OB_score":0.076,"DL_score":0.152}. Since 0.076 < 0.3, this compound is filtered out. All the compounds that pass the screening, along with their related information (CAS number, SMILES string, medicinal material, predicted OB and DL values), are stored in a new data table or data frame. This data structure is the bioactive component profile. This profile represents the chemical basis of the compound that is most likely to be absorbed into the bloodstream and exert biological effects after oral administration.

[0033] As an alternative, the active ingredient screening process in step S1 can employ a chemical rule-based filtering method. Specifically, instead of GCN model prediction, the system can calculate the physicochemical properties of each compound in the initial compound set and screen according to the extended Lipinski's Rule of Five. This process includes calculating the molecular weight (MW), logarithm of the lipid-water partition coefficient (LogP), number of hydrogen bond donors (HBD), and number of hydrogen bond acceptors (HBA) from the SMILES string of each compound using a cheminformatics software library (such as RDKit). The screening rules are set as follows: MW ≤ 500 Da, LogP ≤ 5, HBD ≤ 5, HBA ≤ 10. Only compounds that simultaneously meet all four conditions are retained. This alternative has lower computational costs and clearer rules, but it may not be able to identify some atypical structural molecules that do not conform to the classical rules but still have good oral activity. In contrast, the GCN model used in the embodiments of this invention can learn more complex nonlinear relationships from molecular graph structures, and its prediction accuracy and generalization ability are generally superior.

[0034] Step S2, the generation of the disease-drug target set, can be broken down into the following sub-steps: Step S201: Receive the disease name, query the human tissue gene expression map database, and construct a tissue-specific high-expression gene set. Specifically, the system receives a disease name associated with the indication of the traditional Chinese medicine compound. This name is preferably an International Classification of Diseases (ICD-10) code or a Medical Subject Headings (MeSH) ID to ensure uniqueness. In one example, assuming the indication of the compound is hepatocellular carcinoma, the received disease identifier is ICD-10:C22.0. Based on this identifier, the system first maps it to the associated lesion tissue type, namely "liver".

[0035] The system then accesses a locally mirrored Human Protein Atlas (HPA) PostgreSQL database. This database stores gene expression data in different human tissues based on RNA sequencing. The system executes an SQL query to filter for genes highly expressed in liver tissue. The query logic is: `SELECT gene_uniprot_id, tpm_value FROM rna_tissue_expression WHERE tissue_name='liver'`. The query returns a list containing all genes expressed in liver tissue and their corresponding transcript expression levels per million base pairs (TPM). The system sorts this list in descending order by `tpm_value`.

[0036] Next, a screening process is performed based on a preset percentile threshold for expression levels. In this embodiment, this threshold is set to the top 25%. This parameter is derived from the analysis of GTEx (Genotype-Tissue Expression) project data, which shows that genes defining tissue-specific functions are typically enriched in the top quartiles of the expression profile. Therefore, the system selects genes with TPM values ​​ranking in the top 25% of all liver-expressed genes and extracts their corresponding UniProt protein IDs. These IDs constitute a set, namely the set of tissue-specific highly expressed genes. For example, if a total of 19,670 expressed genes were detected in liver tissue, the 4,917 genes with the highest expression levels are selected. This set may include genes such as EGFR (P00533) and TP53 (P04637) that are actively expressed in the liver and are associated with the occurrence and development of liver cancer.

[0037] Step S202: Based on the bioactive component profile generated in step S1, the system queries the target-compound interaction database to obtain the total potential target set for the compound. Specifically, the system extracts a list of all compound CAS numbers from the bioactive component profile. Using this list, a batch query is initiated to a locally deployed copy of the STITCH (Search Tool for Interactions of Chemicals) database, which contains large-scale compound-protein interaction information. The purpose of the query is to identify all protein targets that have known interactions with these active components. To ensure the reliability of the interactions, a filter condition is added to the query, for example, requiring that the combined_score (a binding strength score that integrates experimental evidence, text mining, database co-references, etc., ranging from 0-1000) recorded in the STITCH database is greater than a preset confidence threshold, such as 400.

[0038] As a simulation example, for the active ingredient ferulic acid (CAS: 1135-24-6), querying the STITCH database might return protein targets with a combined_score > 400, including PTGS2, SOD1, and CAT. For another active ingredient, quercetin, possible targets include TP53, MAPK1, and AKT1. The system will aggregate the UniProt IDs of all protein targets retrieved from all bioactive ingredients and perform deduplication to form a single, unique set of IDs. This set represents the total potential target set of the compound. This set may contain thousands of protein targets, representing the entire target space that all active ingredients in the compound could theoretically act upon.

[0039] Step S203 involves performing a set intersection operation between the tissue-specific highly expressed gene set and the set of potential total targets for the compound drug, generating a highly correlated disease-drug target set. Specifically, to perform this operation efficiently, the system employs a hash set-based intersection algorithm. First, the tissue-specific highly expressed gene set (a list of UniProt IDs) generated in step S201 and an empty hash set A are used as input. The list is traversed, and each UniProt ID is added to the hash set A. The average time complexity of this operation is O(n), where n is the number of highly expressed genes. Next, the set of potential total targets for the compound drug (another list of UniProt IDs) generated in step S202 and an empty result list B are used as input. This list is traversed, and for each UniProt ID, it is checked whether it exists in the hash set A. Since the average time complexity of the hash set lookup operation is O(1), the average total time complexity of this step is O(m), where m is the number of potential total targets. If it exists, the ID is added to the result list B. Ultimately, the contents of list B represent the intersection of the two sets.

[0040] This calculation selects targets that simultaneously meet two conditions: first, they are potential targets of the orally administered active ingredients in the compound; second, they are significantly highly expressed in the specific tissues where the disease occurs. The intersection of these two conditions greatly narrows the target range, eliminating a large number of targets irrelevant to the disease tissue background, thus generating a smaller set of disease-drug targets that are biologically more closely related to the disease-drug mechanism of action. For example, if EGFR (P00533) is both in the set of genes highly expressed in the liver and in the list of targets for quercetin, then it will be included in the final set of disease-drug targets.

[0041] It should be noted that this step includes an exception handling mechanism. If the intersection operation results in an empty set, it indicates that no target satisfying both conditions could be found under the current parameter settings. The system will trigger a rollback process. First, the definition of tissue-specific high expression is automatically relaxed, for example, the TPM expression percentile threshold is lowered from the top 25% to the top 40%, and then steps S201 and S203 are re-executed. If the intersection is still empty, the confidence threshold for the STITCH database query is further lowered, for example, from 400 to 250, and steps S202 and S203 are re-executed. If a non-empty target set cannot be obtained after multiple adjustments, the system will generate an alarm message, prompting the user that the association between the traditional Chinese medicine compound and the direct target of the specified disease under the current analysis framework is weak, and suggesting checking the input data or adjusting the analysis strategy.

[0042] As an alternative, the disease-related target acquisition process in step S2 can be performed without relying on tissue expression profiling, instead directly querying disease-gene association databases. Specifically, after receiving the disease name, the system queries databases such as DisGeNET or OMIM that specifically record the relationship between genes and human diseases. Using the disease name or its ID, a list of known genes genetically or functionally associated with the disease is directly obtained. This list constitutes a disease-associated gene set. Subsequently, this disease-associated gene set is intersected with the set of potential total targets obtained in step S202. This alternative directly utilizes validated disease-gene knowledge to locate targets related to the genetic basis of the disease. In contrast, the tissue expression profiling-based method used in this embodiment focuses on identifying targets that are functionally active under specific pathophysiological states (i.e., in diseased tissues), better reflecting the tissue-specific phenotype of the disease, and potentially revealing some equally important functional targets not included in classic disease gene lists.

[0043] Please see Figure 3 Step S3, identification of synergistic modules, aims to discover combination patterns of multiple compounds synergistically acting on multiple functionally related targets from the target list generated in the previous step. The initial part of this step specifically involves: Step S301: Construct a weighted bipartite graph network model. Specifically, the system instantiates a graph data structure object that supports the storage of node and edge attributes. In a specific embodiment, this object can be created using Python's NetworkX library. The graph construction process is as follows: First, nodes are created. The system iterates through the bioactive component spectrum generated in step S1, adding each compound as a node to the graph. The node attributes include at least: node_id (unique identifier, which can use a CAS number), node_type (set to 'compound'), and node_label (the generic name of the compound). Next, the system iterates through the disease-drug target set generated in step S2, adding each target (represented by a UniProt ID) as another type of node to the graph. The attributes of this type of node include: node_id (UniProt ID), node_type (set to 'target'), and node_label (the gene symbol corresponding to the target). At this point, the node set of the graph is completed, consisting of two disjoint subsets of nodes.

[0044] The next step is to add connecting edges and assign weights to the graph data structure. Specifically, the system iterates through all possible combinations of compounds in the bioactive component spectrum from step S1 and all targets in the disease-drug target set from step S2. For each "compound-target" pair, the system queries the local STITCH database again to determine if there is a recorded interaction between them. If an interaction exists and its combined_score is greater than a preset confidence threshold (e.g., 400), a connecting edge is established between the corresponding compound node and target node in the graph. Subsequently, the numerical combined_score extracted from the STITCH database is used as the weight value of this connecting edge and stored in the edge's attributes. For example, if the combined_score of the compound quercetin and the target TP53 is 950, an edge is created between the node representing quercetin and the node representing TP53, and its weight attribute is set to 950. This process is repeated until all compound-target interactions that meet the conditions are represented as weighted edges in the graph. At this point, a complete weighted bipartite graph network representing the intensity of the interaction between the active ingredient and disease-related targets has been constructed.

[0045] Step S302 involves applying the Weighted Stochastic Block Model (wSBM) algorithm to perform community detection analysis on the constructed weighted bipartite graph network. Specifically, the objective function of this algorithm is to maximize the posterior probability of the graph under a given community partitioning, i.e., to find a node grouping method such that the distribution pattern of the sum of connection weights between nodes within a group and the sum of connection weights between nodes between groups best matches the generation model of wSBM. In a specific embodiment, the system uses an optimizer based on the Markov Chain Monte Carlo (MCMC) method to solve this problem. The optimization process iteratively and randomly moves nodes from one community to another, and decides whether to accept the move based on the change in the posterior probability of the graph after the move. Specifically, a Gibbs sampling strategy is used. In each iteration, for each node, the conditional probability of it being assigned to each possible community is calculated, and a new community affiliation is randomly selected for it based on these probabilities.

[0046] The algorithm's convergence is determined by monitoring changes in the overall partition entropy or modularity score of the graph. The algorithm terminates when, after multiple iterations (e.g., 1000), the community partitioning structure no longer changes significantly, or when the likelihood function describing the entire graph reaches a stable plateau. The algorithm's output is a mapping table that assigns all nodes in the graph (including compound and target nodes) to different communities.

[0047] Step S303 involves parsing and filtering the communities output by the algorithm, identifying communities that meet specific criteria as synergistic modules. Specifically, within the wSBM framework, a community's weighted density of internal connections (i.e., the sum of all edge weights within the community) is significantly higher than the connection density based on the model's random expectation. The system iterates through each community (i.e., sub-network structure) output by the algorithm. For each community, the system performs structural filtering. The filtering criteria include: first, the community must contain at least two types of nodes, i.e., it must contain both compound nodes and target nodes; second, the community's size must reach a preset minimum size, for example, containing at least 2 compound nodes and 3 target nodes. This criterion aims to filter out substructures that are too simple (e.g., a single compound acting on a single target) or too small to be statistically significant. Each community that passes the above filtering is formally defined as a synergistic module. Each synergistic module represents a potential action unit in which multiple bioactive components act together on a set of functionally related disease targets.

[0048] Please see Figure 4 Step S4, the generation of predictive pathway perturbation maps, aims to transform static synergistic modules into dynamic biological process simulations. This step specifically involves: Step S401: For each identified synergistic module, the set of targets within it is mapped to a biological pathway database to determine the core signaling pathways primarily regulated by that module. Specifically, the system extracts a list of UniProt IDs for all target nodes within a synergistic module. Subsequently, using this list, a pre-configured, localized Kyoto Encyclopedia of Genes and Genomes (KEGG) database is accessed programmatically. The system performs a batch query to retrieve the IDs and names of all KEGG signaling pathways involved by these target proteins. The query returns a list of pathways, which may contain duplicate pathway IDs because multiple targets may participate in the same pathway.

[0049] The system performs frequency statistics and sorts the returned pathway ID list, selecting the pathway with the highest frequency as the core signaling pathway of the module. In one example, assuming a synergistic module contains 5 targets, the pathway list returned after querying KEGG is: {hsa04151, hsa04010, hsa05200, hsa04151, hsa04066, hsa04151}. Statistical analysis shows that hsa04151 (PI3K-Akt signaling pathway) appears 3 times, having the highest frequency. Therefore, the system determines that the PI3K-Akt signaling pathway is the core signaling pathway regulated by this synergistic module. As an enhancement strategy, a hypergeometric distribution test can be used for pathway enrichment analysis, calculating the significance (p-value) of each pathway in the module's target set, and selecting the pathway with the lowest p-value (most significant) as the core signaling pathway. This method is more rigorous than simple frequency statistics.

[0050] Step S402 involves obtaining the topological structure of the core signaling pathway and transforming it into a dynamic system model described by a system of ordinary differential equations (ODEs). Specifically, the system downloads the KEGG Markup Language (KGML) file for the pathway based on the core signaling pathway ID (e.g., hsa04151) determined in step S401. KGML is an XML format file that details all molecular nodes (genes, proteins, compounds) in the pathway and their interrelationships (e.g., activation, inhibition, phosphorylation, dephosphorylation, binding, etc.). The system deploys a KGML parser to read the file, extract all node information and edge (i.e., interaction) information, thereby constructing a directed graph representing the pathway topology in memory.

[0051] Next, the system automatically transforms this directed graph topology into a system of ordinary differential equations. Each molecular component in the pathway (e.g., the non-phosphorylated and phosphorylated forms of a protein) is treated as a state variable, and its concentration change over time is described by a series of ordinary differential equations. The equations are constructed following the law of mass action or more complex enzyme kinetic models (such as the Michaelis-Menten equation). For example: For a simple activation reaction A→B (A activates B), the equation describing the change in concentration of the active form B_active is d[B_active] / dt=k_act*[A_active]*[B_inactive]-k_deact*[B_active], where k_act is the activation rate constant and k_deact is the deactivation rate constant.

[0052] For an inhibitory reaction C-|B (C inhibits B), an inhibition term will be added to the equation, for example, d[B_active] / dt=...-k_inhibit*[C_active]*[B_active].

[0053] The phosphorylation reaction Kinase + Substrate → Kinase + Substrate-P can be described by the Michaelis-Menten equation: d[Substrate-P] / dt = Vmax*[Substrate] / (Km+[Substrate])-k_dephospho*[Substrate-P].

[0054] The kinetic parameters required in the equation system (such as the reaction rate constant k, Michaelis constant Km, maximum reaction rate Vmax, degradation rate, etc.) are preferentially obtained from publicly available kinetic parameter databases (such as SABIO-RK or BRENDA). For parameters missing from the database, the system uses reasonable default values ​​or estimates based on literature, allowing users to make subsequent adjustments. Ultimately, a system containing dozens or even hundreds of coupled ordinary differential equations is formed, which fully describes the baseline dynamic behavior of the core signaling pathway in the absence of external perturbations.

[0055] Step S403 uses the targeting effect of compounds in the synergistic module as the initial perturbation condition, and simulates and generates a predictive pathway perturbation map by numerically solving a system of ordinary differential equations. Specifically, the system analyzes the interaction type between compounds in the synergistic module and their targets. This type information can be obtained from the STITCH database or other pharmacology databases and is usually labeled as "activation" or "inhibition / binding". The system converts this qualitative effect into a quantitative perturbation of ODE model parameters or state variables.

[0056] In a specific embodiment, if a compound activates target A, then at the start of the simulation (t=0), the initial concentration [A_active] representing the active form of target A in the model is increased by a preset factor (e.g., increased by 50%), or the rate constant k_act of its catalyzed downstream reaction is increased by a fixed proportion. Conversely, if the compound inhibits, the corresponding parameter is decreased. If multiple compounds in the module act on multiple different targets in the pathway, all corresponding parameters or initial conditions are perturbed simultaneously to simulate the synergistic effect of multiple components.

[0057] After setting the initial perturbation conditions, the system invokes a high-performance numerical integrator deployed on a computing server. In one embodiment, the `solve_ivp` function from the Python scientific computing library SciPy is used, which implements several high-precision adaptive step-size solution algorithms, including the fourth-order Runge-Kutta method (RK45). The system sets a virtual simulation time range, for example, from 0 to 48 hours, and specifies an output time sequence (e.g., outputting results every 0.5 hours). The solver performs numerical integration on the entire system of ordinary differential equations, calculating the dynamic changes in the concentrations of all molecular components in the pathway within the preset time sequence.

[0058] The solution yields an N×M dimensional matrix, where N represents the number of molecules in the pathway and M represents the number of time points. Each element C(i,j) in the matrix represents the predicted concentration or activity level of molecule i at time point j. This matrix is ​​a quantified, time-series-informed predictive pathway perturbation map. This map visually demonstrates how the core signaling pathway is perturbed under the synergistic effect of the components of the traditional Chinese medicine compound, and how this perturbation propagates downstream along the pathway, ultimately affecting the behavior of key effector molecules.

[0059] Please see Figure 5 Step S5, experimental evidence-based verification and confirmation of the core mechanism of action, aims to evaluate and confirm the reliability of the computational predictions using real biological experimental data. This step specifically involves: Step S501: Receive and process externally input experimental omics data to generate a pre-sorted gene list. Specifically, the system receives an experimental omics data file uploaded by the user, obtained after intervention with the traditional Chinese medicine compound on a disease-related biological model (e.g., the HepG2 liver cancer cell line treated with "Siwu Tang," or a corresponding animal model). This file is typically in text format (e.g., CSV or TSV) and contains a list of differentially expressed genes or proteins. Each line of the file represents a gene or protein and contains at least three columns of key information: gene / protein ID (e.g., gene symbol or UniProt ID), the logarithm of the expression change (log2FoldChange), and the p-value or false discovery rate (FDR / q-value) indicating the statistical significance of the expression change.

[0060] The system first preprocesses the data list. The first step is filtering based on statistical significance, removing gene / protein entries with p-values ​​or FDRs greater than a preset threshold (e.g., 0.05), as the expression changes of these entries may be caused by random fluctuations. The second step involves strictly sorting the filtered gene / protein list in descending order according to their log2 FoldChange values. This results in the most significantly upregulated genes at the top of the list and the most significantly downregulated genes at the bottom. This filtered and sorted list constitutes the pre-sorted experimental gene list, forming the basis for subsequent enrichment analysis.

[0061] Step S502: Based on the predictive pathway perturbation map generated in step S4, extract the set of downstream perturbed molecules to be validated. Specifically, the system analyzes the predictive pathway perturbation map (i.e., an N×M dimensional concentration matrix). To determine which molecules undergo significant changes under drug perturbation, the system calculates the concentration change rate of each molecule in the pathway from the simulation end point (t_end) to the simulation start point (t=0), using the formula (Concentration_at_t_end - Concentration_at_t_0) / Concentration_at_t_0. The system filters out molecules whose absolute concentration change rate is greater than a preset perturbation threshold (e.g., 0.5, indicating a concentration change exceeding 50%). These filtered molecules represent the effector molecules that exhibit the most significant downstream response after being perturbed by the drug synergistic module in the computational model. The gene IDs of these molecules are extracted to form a set, i.e., the set of downstream perturbed molecules to be validated. This set is the core "fingerprint" of the computational prediction.

[0062] Step S503: Apply the Pre-sorted Gene Set Enrichment Analysis (GSEA) algorithm to calculate the enrichment status of the gene set to be validated within the sorted experimental data list. Specifically, the system calls a locally implemented GSEA program module. This module receives two inputs: the pre-sorted experimental gene list generated in step S501, and the downstream perturbed molecule set to be validated generated in step S502. The core logic of the algorithm is to calculate an enrichment score (ES). The calculation process is as follows: The program starts from the top of the pre-sorted experimental gene list and traverses downwards one by one. When a gene is encountered that belongs to the molecule set to be validated, an accumulated running score is increased; when a gene is encountered that does not belong to the set, the score is decreased. The increment step is related to the gene's sorting position in the experimental list (or its log2FoldChange value), while the decrement step is a fixed penalty value. During the entire traversal, the maximum peak (or minimum trough) of the running score is the original enrichment score ES of the gene set. A positive ES value indicates that the gene set is enriched at the top of the sorting list (i.e., it tends to be upregulated in experiments), while a negative ES value indicates that it is enriched at the bottom (i.e., it tends to be downregulated).

[0063] Step S504: The enrichment score is standardized using a permutation test, and a comprehensive validation score is calculated. Specifically, to assess the statistical significance of the observed ES, the system performs a permutation test procedure. This procedure simulates the ES distribution under a null hypothesis (i.e., no association between the gene set and the phenotype) by performing multiple (e.g., 1000) random permutations on the gene tags in the pre-sorted experimental gene list. For each permutation, the ES is recalculated, thus constructing a null distribution of the ES. The standardized enrichment score (NES) is calculated by comparing the original, unpermuted ES with this null distribution. The NES is obtained by dividing the original ES by the mean of the null distribution, taking into account differences in gene set size. Simultaneously, based on the null distribution, the false discovery rate (FDR q-value) can be calculated, representing the probability that the observed current or stronger enrichment signal is a random event.

[0064] The system then calculates a comprehensive validation score (V) based on NES and FDR. In a specific embodiment, this score is calculated as: V = NES * (-10 * log10(FDR)). This formula combines the direction and intensity of enrichment (represented by NES) with the statistical confidence level (represented by FDR). The smaller the FDR, the larger the value of the -10 * log10(FDR) term, and the stronger the weighting effect on NES.

[0065] Step S505: Based on the comprehensive verification score, the predicted core mechanism of action is confirmed, and a final analysis report is output. Specifically, the system sets a compliance threshold, such as 2.0. When the comprehensive verification score V calculated for a certain synergistic module is greater than this threshold, the system determines that the mechanism of action constituted by the module and its regulated core signaling pathway has been supported and verified by real experimental data. This synergistic module is marked as a "verified core mechanism of action".

[0066] Finally, the system integrates all analysis results into a structured analysis report and outputs it in portable document format (PDF). This report includes at least the following: a list of screened bioactive components; a network topology diagram of validated core synergistic modules (visualizing multi-component-multi-target relationships); dynamic perturbation simulation curves of core signaling pathways (showing changes in the concentration of key downstream molecules over time); an enrichment map of GSEA (i.e., a "mountain map," visually illustrating the enrichment process); and all relevant quantitative scores, including module identification parameters, pathway enrichment p-values, GSEA NES, FDR, and the final comprehensive validation score V. This report provides a complete, multi-layered chain of scientific evidence for the mechanism of action of traditional Chinese medicine compound formulas, from chemical components to network modules, pathway dynamics, and experimental validation.

[0067] Please see Figure 2 This invention also provides a system for analyzing the mechanism of action of traditional Chinese medicine based on network pharmacology, used to perform the above-described methods. The functional block diagram of this system is shown below. Figure 2 As shown, the system includes: The bioactive ingredient screening module is responsible for performing the functions of step S1. Internally, it integrates a connector to a traditional Chinese medicine (TCM) ingredient database and a service that encapsulates a pre-trained graphical convolutional neural network (ADMET) prediction model. The module's input interface receives TCM compound data and outputs a structured bioactive ingredient profile. The module's internal processor performs querying, molecular graph transformation, model invocation, and filtering logic based on dual thresholds of oral bioavailability and drug-likeness.

[0068] The tissue-specific target enrichment module is responsible for performing the functions of step S2. Internally, it integrates connectors with human tissue gene expression profile databases (such as HPA) and target-compound interaction databases (such as STITCH). The module receives disease names and bioactive component profiles as input. Its core logic unit performs screening for tissue-highly expressed genes, querying for potential compound targets, and performing efficient hash set intersection operations, ultimately outputting a tissue-specific set of drugs and targets highly correlated with the disease phenotype.

[0069] The synergistic effect module discovery module is responsible for performing the functions of step S3. This module receives a bioactive component profile and a disease-drug target set as input. Internally, it contains a graph construction engine for generating a weighted compound-target bipartite graph network. The core component is an implementation of a weighted randomized block model algorithm, which performs community detection analysis on the constructed network. The module ultimately outputs a set of statistically significant tightly connected subnetwork structures identified as synergistic effect modules.

[0070] The pathway dynamics perturbation simulation module performs the functions of step S4. This module receives a synergistic module as input and connects to a biological pathway database (such as KEGG). Internally, it includes a pathway enrichment analysis unit, a KGML file parser, and an automatic generator of ordinary differential equations. A key component is a high-performance numerical solver used to simulate the dynamic response of core regulatory pathways under compound perturbation. The module's output is a quantitative, time-series-informed predictive pathway perturbation map.

[0071] The experimental evidence-based validation module is responsible for executing the functions in step S5. Its input interface receives externally uploaded experimental omics data and the predicted maps generated by the pathway dynamic perturbation simulation module. At its core is a pre-sorted gene set enrichment analysis (GSEA) algorithm engine, which performs sorting, gene set extraction, enrichment score calculation, permutation testing, and comprehensive validation score calculation. Based on the comparison of this score with a preset threshold, this module provides final confirmation or refutation of the predicted synergistic mechanism, integrates all analysis results, and generates and outputs a detailed final analysis report containing multi-dimensional visualization charts.

[0072] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for analyzing the mechanism of traditional Chinese medicine based on network pharmacology, characterized in that, The method comprises the following steps: S1: obtaining traditional Chinese medicine compound data, calling a pre-trained graph convolutional neural network pharmacokinetic property prediction model to calculate the oral bioavailability and drug-likeness of each compound, filtering according to a preset threshold, and screening to form a bioactive ingredient spectrum; S2: based on the disease name associated with the indications of the traditional Chinese medicine compound, querying the human tissue gene expression atlas database and the target-compound interaction database, extracting a set of high-expression genes specific to the diseased tissue, and performing intersection operation with a set of potential target points of the bioactive ingredient spectrum to generate a disease-drug target point set; S3: based on the bioactive ingredient spectrum and the disease-drug target point set, a weighted compound-target bipartite graph network is constructed, a weighted random block model algorithm is applied for community discovery, a sub-network structure composed of multiple compounds synergistically acting on a functionally related target point cluster is identified, and the synergistic action module is determined; S4: for the synergistic action module, the core signal pathway regulated thereby is mapped and determined, a dynamic model of the ordinary differential equation set of the core signal pathway is constructed, and the target action of the synergistic action module is taken as an initial disturbance condition for numerical solution, and a predictive pathway disturbance atlas is generated. 2.The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 1, characterized in that, The bioactive ingredient spectrum comprises compound molecular structure identifiers, oral bioavailability prediction values, drug-likeness prediction values, and information about the medicinal materials, the disease-drug target point set comprises UniProt IDs and gene symbols of protein targets, the synergistic action module comprises compound and target point sets as nodes and weighted interaction relationships as connecting edges, and the predictive pathway disturbance atlas comprises a quantitative data matrix of the concentration or activity of each molecular component in the pathway changing over time. 3.The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 1, characterized in that, The obtaining step of the bioactive ingredient spectrum is specifically: S111: receiving traditional Chinese medicine compound data composed of medicinal material names, querying a traditional Chinese medicine ingredient database, extracting all chemical ingredients corresponding to each medicinal material and their molecular structure identifiers to form an initial compound set; S112: converting the molecular structure identifiers of each compound in the initial compound set into graph representation data structures, and inputting them into a pre-trained graph convolutional neural network model to respectively output oral bioavailability prediction probability values and drug-likeness prediction probability values of each compound; S113: based on preset oral bioavailability threshold values and drug-likeness threshold values, double-filtering the prediction probability values, screening out compounds higher than the corresponding threshold values, integrating their molecular information and prediction values, and generating the bioactive ingredient spectrum.

4. The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 3, characterized in that, The obtaining step of the disease-drug target point set is specifically: S211: based on the disease name associated with the disease tissue type, querying the human tissue gene expression atlas database to screen out genes with gene expression higher than a preset expression level percentage in the diseased tissue to form a set of tissue-specific highly expressed genes; S212: based on the bioactive ingredient spectrum, querying the target-compound interaction database to obtain all potential action targets corresponding to each active ingredient, and performing set intersection operation with the set of tissue-specific highly expressed genes to obtain the disease-drug target point set.

5. The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 4, characterized in that, The obtaining step of the synergistic module is specifically as follows: S311: taking the compounds in the bioactive ingredient spectrum and the targets in the disease-drug target set as two types of nodes respectively, assigning weights to the connection edges according to the binding strength scores in the target-compound interaction database, and constructing a weighted bipartite graph network model; S312: applying a weighted random block model algorithm to the weighted bipartite graph network model for community discovery analysis, identifying and outputting a subnetwork structure with an internal weighted connection density significantly higher than an intercommunity connection density, and determining the subnetwork structure as the synergistic module. 6.The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 5, characterized in that, The obtaining step of the predictive pathway perturbation atlas is specifically as follows: S411: mapping the target set in the synergistic module to a biological pathway database, determining a core signaling pathway regulated by the module based on the distribution frequency of the targets in each pathway, and determining a core signaling pathway regulated by the module based on the distribution frequency of the targets in each pathway; S412: obtaining the topological structure of the core signaling pathway, converting each molecular component and its mutual activation or inhibition relationship in the pathway into a set of ordinary differential equations, and constructing a pathway dynamic system model; S413: quantifying the activation or inhibition effect of the compounds on the targets in the synergistic module as an initial perturbation condition of the set of ordinary differential equations, solving the set of ordinary differential equations by a numerical integration algorithm, and generating the predictive pathway perturbation atlas. 7.The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 1, characterized in that, The method further comprises a step S5: S5: calling the predictive pathway perturbation atlas, comparing the predictive pathway perturbation atlas with experimental omics data measured on a biological model for the traditional Chinese medicine compound, applying a pre-ranked gene set enrichment analysis algorithm to calculate a comprehensive verification score, and confirming a core action mechanism represented by the synergistic module; The comprehensive verification score includes a normalized enrichment score, a false discovery rate, and a verification determination result based on a compliance threshold. 8.The network pharmacology-based mechanism of action analysis method of traditional Chinese medicine according to claim 7, characterized in that, The confirmation step of the core action mechanism is specifically as follows: S511: receiving the experimental omics data, filtering the experimental omics data according to a statistical significance value, and arranging the experimental omics data in descending order according to the logarithmic value of the expression change fold to generate a pre-ordered experimental gene list; S512: based on the predictive pathway perturbation atlas, extracting downstream molecules with a concentration or activity change rate exceeding a preset perturbation threshold to form a gene set to be verified; S513: taking the pre-ordered experimental gene list and the gene set to be verified as inputs, executing a pre-ranked gene set enrichment analysis algorithm, calculating a normalized enrichment score and a false discovery rate, and generating the comprehensive verification score, and confirming the core action mechanism when the score is higher than a compliance threshold.

Images