Network Pharmacology + Molecular Docking + Molecular Dynamics: How to Do the Validation Section

Summary

Five things that make the validation section hold up

When reviewers read the validation section, they mainly check the five points below. The first three have to be settled before any computation starts.

The interpretation paper of the WFCMS Network Pharmacology Evaluation Method Guidance surveyed network pharmacology papers indexed in CNKI (up to June 2021): fewer than 30% validated their results experimentally, about one third of those with in vitro or in vivo experiments, and most studies only inferred their results from the literature. Reviewers have seen many manuscripts that rely on docking alone, and docking is no longer accepted as validation.

Full workflow

From compound to figure: input, tool, output and figure for each step

Tool versions are current as of 2026-10-10. Steps 1, 9, 10 and 12 decide whether the validation section holds up; for steps 2 to 8, reviewers mainly check thresholds and versions.

Databases

Common databases: 2026 status and practical experience

The table reflects actual access on 2026-10-10. In the methods, give three items for each database: version or URL, access date and screening criteria.

The WFCMS guidance interpretation paper examined 13 network pharmacology papers on ginseng in CNKI: even with the same database and the same thresholds, the number of included compounds differed, because thresholds were changed and compounds added from the literature were not explained. In TCMSP, 22 ginseng compounds have OB > 30% and DL > 0.18, and 17 of them have predicted targets. Keeping the raw export tables and providing them as supplementary files is the most direct way to let others reproduce your compound and target counts.

Thresholds

Screening thresholds: common settings, sources and how to report them

Most thresholds have no agreed standard. The reliable approach is to choose a value with a source, fix it before seeing the results, and state the reason in the methods.

PPI network

Target overlap and core targets: the median-degree approach and its limits

The usual approach is to harmonize compound targets and disease genes to HGNC symbols, take the overlap, submit it to STRING, compute topology in Cytoscape, and take nodes with degree at or above the median as core targets.

Core targets obtained this way are highly uniform. Diao et al. (2026) surveyed 1,038 natural-product network pharmacology studies published from October 2023 to June 2024: AKT1 appeared among the core targets in 48.8% of studies, TNF in 39.1% and EGFR in 34.8%; AKT1 appeared in 54.9% of database-only studies and in 19.6% of studies that integrated omics or experimental data. Quercetin appeared in 71.9% and 42.3% of the two groups, respectively.

One reason is that STRING text-mining evidence favors well-studied genes. The script below first checks which protein each gene symbol maps to with get_string_ids, then fetches a network from a pinned STRING API version, flags edges supported only by text mining and lists core targets by median degree. Tested on this page (2026-10-10, STRING 12.5 API, Python 3 standard library): with 10 genes that recur in network pharmacology papers (AKT1, TNF, IL6, PTGS2, ESR1, STAT3, MAPK1, CASP3, BCL2, VEGFA) at required_score=700, the network had 33 edges, 5 of them supported only by text mining; the 0.88 score for AKT1–IL6 came entirely from text mining. The 9 connected nodes had a median degree of 8, and 5 nodes tied at or above the median.

The same test showed that in STRING 12.0 and 12.5 a query for VEGFA returns COL18A1 (9606.ENSP00000352798) without any error; queries by UniProt accession P15692 or Entrez ID 7422 return nothing, while STRING 11.5 still matches VEGFA correctly. If the gene list is submitted directly to the network endpoint, VEGFA enters the network as COL18A1 (at required_score=400, COL18A1 appears with 3 edges). VEGFA is a frequent core target in network pharmacology papers, so check the mapping before submission and report genes with inconsistent mapping separately in the methods.

Ways to reduce this uniformity: check whether core targets change in the disease using transcriptomic or proteomic data; report betweenness and other metrics alongside degree; and provide the full target list as a supplement so reviewers can judge.

# Usage: python3 string_hubs.py overlap_genes.txt 700
# overlap_genes.txt: one overlapping gene symbol per line; second argument is required_score (0-1000)
import sys, csv, io, statistics, urllib.parse, urllib.request
from collections import defaultdict

API = "https://version-12-5.string-db.org/api/tsv"


def post(method, **params):
    params.update({"species": 9606, "caller_identity": "np_validation_tutorial"})
    data = urllib.parse.urlencode(params).encode()
    with urllib.request.urlopen(f"{API}/{method}", data=data) as r:
        return list(csv.DictReader(io.StringIO(r.read().decode()), delimiter="\t"))


genes = [g.strip() for g in open(sys.argv[1]) if g.strip()]
score = int(sys.argv[2]) if len(sys.argv) > 2 else 400

# 1. Check symbol mapping first: STRING silently maps unmatched symbols to other proteins
mapped = post("get_string_ids", identifiers="\r".join(genes), echo_query=1, limit=1)
hit = {m["queryItem"]: m for m in mapped}
for g in genes:
    if g not in hit:
        print(f"WARNING not found in STRING: {g}")
    elif hit[g]["preferredName"].upper() != g.upper():
        print(f"WARNING {g} mapped to {hit[g]['preferredName']} ({hit[g]['stringId']})")
ids = [hit[g]["stringId"] for g in genes if g in hit and hit[g]["preferredName"].upper() == g.upper()]

# 2. Fetch the network using only consistently mapped STRING IDs
rows = post("network", identifiers="\r".join(ids), required_score=score)

edges = {}
for x in rows:
    key = tuple(sorted((x["preferredName_A"], x["preferredName_B"])))
    other = max(float(x[c]) for c in ("nscore", "fscore", "pscore", "ascore", "escore", "dscore"))
    edges[key] = (float(x["score"]), float(x["tscore"]), other)

tm_only = [k for k, (s, t, o) in edges.items() if o == 0]
print(f"edges={len(edges)}  textmining_only={len(tm_only)}")
for a, b in tm_only:
    print(f"  textmining only: {a}-{b} score={edges[(a, b)][0]}")

deg = defaultdict(int)
for a, b in edges:
    deg[a] += 1
    deg[b] += 1
med = statistics.median(deg.values())
print(f"nodes={len(deg)}  median_degree={med}")
for g, d in sorted(deg.items(), key=lambda kv: -kv[1]):
    print(f"{g}\t{d}\t{'>=median' if d >= med else ''}")

python3 string_hubs.py overlap_genes.txt 700

Enrichment

GO/KEGG enrichment: report the background gene set and multiple-testing correction

The two common gaps in enrichment analysis are an unreported background gene set and missing multiple-testing correction.

Wijesooriya et al. (2022) reviewed 186 articles that used enrichment analysis: 95% of over-representation analyses did not use or did not describe an appropriate background gene set, and 43% did not correct for multiple testing. In network pharmacology studies, the AGE-RAGE and PI3K-Akt pathways appear in 54.9% and 54.7% of studies, respectively (Diao et al. 2026).

A workable approach is to use the union of all compound targets and disease genes as the background and to report the background size, correction method and KEGG access date in the methods. In clusterProfiler, the ont argument of enrichGO() defaults to "MF", so omitting it returns molecular function only; enrichKEGG() downloads current KEGG data online by default, so record the access date.

Tested on this page (2026-10-10, R 4.5.3, clusterProfiler 4.18.4, org.Hs.eg.db 3.22.0, macOS arm64): with the 10 genes above as the overlap set and their 191 STRING 12.5 interaction partners as the background, the script ran successfully. With this background, 170 GO BP terms and 20 KEGG pathways were significant; without universe (the default background of all annotated genes), the counts were 1,493 and 132; omitting ont returned 77 MF terms. Under both backgrounds, the top KEGG pathway was AGE-RAGE signaling pathway in diabetic complications. The example background is for demonstration only; in a real analysis use the union of all compound targets and disease genes.

# Dependencies: BiocManager::install(c("clusterProfiler", "org.Hs.eg.db"))
library(clusterProfiler)
library(org.Hs.eg.db)

hits <- readLines("overlap_genes.txt")     # overlap genes (symbols)
bg   <- readLines("background_genes.txt")  # background: union of all compound targets and disease genes

hits_id <- bitr(hits, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
bg_id   <- bitr(bg,   fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

ego <- enrichGO(gene = hits_id, universe = bg_id, OrgDb = org.Hs.eg.db,
                keyType = "ENTREZID", ont = "BP",  # ont defaults to "MF"
                pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2,
                readable = TRUE)
ekegg <- enrichKEGG(gene = hits_id, universe = bg_id, organism = "hsa",
                    keyType = "kegg", pAdjustMethod = "BH",
                    pvalueCutoff = 0.05, qvalueCutoff = 0.2)

write.csv(as.data.frame(ego),   "go_bp.csv", row.names = FALSE)
write.csv(as.data.frame(ekegg), "kegg.csv",  row.names = FALSE)
writeLines(c(paste("hits", length(hits_id)),
             paste("background", length(bg_id)),
             paste("KEGG accessed", Sys.Date())), "enrichment_meta.txt")

Molecular docking

Docking validation: where the cutoffs come from, positive controls and redocking

In a network pharmacology paper, docking provides structural plausibility for the core compounds and core targets. A docking score does not prove binding.

Che and Zhang (2025) reviewed 35 network pharmacology papers published in 2024: about 15 used −5.0 kcal/mol as an affinity cutoff and the others used values from −1.2 to −7.0, none with a traceable source; scores from different docking programs are not comparable; none of the 35 papers used a positive control; 4 explicitly used blind docking, and in another 13 different ligands docked to different regions of the protein, suggesting blind docking. In their test on 285 complexes from CASF-2016, Vina scores from site-specific docking correlated with experimental affinity at r = 0.604, against 0.387 for blind docking; for blind docking, the median RMSD was 3.37 Å at exhaustiveness 8 and 2.21 Å at 64, at about 8 times the run time.

Receptor and ligand preparation, box setup and complete commands are covered in /learn/autodock-vina-docking.

1. 01

   Choose the structure

   Prefer human crystal structures with a co-crystallized ligand and a complete binding pocket, and report the PDB ID and resolution. If only a predicted structure is available, check the pLDDT of the pocket residues first (see /learn/alphafold-results).

2. 02

   Define the site

   Define the docking box from the co-crystallized ligand position or from functional residues reported in the literature. If the site is unknown, run pocket prediction first and cross-check with a second docking program.

3. 03

   Redock

   Remove the co-crystallized ligand and redock it with the same settings. A heavy-atom RMSD ≤ 2 Å between the top pose and the crystal pose shows that the settings reproduce the known binding mode (PoseBusters uses the same criterion).

4. 04

   Positive control

   Dock a known inhibitor of the target, or the co-crystallized ligand, in the same box with the same settings. Compare compound scores with the positive control, and take only compounds that score close to or better than it into MD.

5. 05

   Report the settings

   Report the box center and size, exhaustiveness, random seed and Vina version, with 2D/3D interaction diagrams and the redocking overlay.

Molecular dynamics

Molecular dynamics validation: replicas, metrics and MM/PBSA

Molecular docking and molecular dynamics are two successive steps: docking searches ligand poses on a nearly rigid receptor and scores them, and molecular dynamics tests whether a pose persists over time in explicit solvent. If the docking site is wrong, MD cannot correct it.

There is no agreed simulation length. The MD reliability checklist published in Communications Biology in 2023 asks for evidence of convergence and at least 3 independent simulations per condition; it does not set a fixed number of nanoseconds. A common practical way to judge equilibration: least-squares fit each frame on the protein backbone to a reference structure, compute the RMSD, apply a running average to the curve, and treat the structure as equilibrated only from the point where the curve no longer rises or falls overall.

Known limitations of MM/PBSA (Genheden and Ryde 2015): the energies are usually too large; ranking is often unreliable when ligand affinities differ by less than 12 kJ/mol (about 2.9 kcal/mol); r² against experiment is about 0.3 across PDBbind and 0 to 0.8 for individual proteins; the entropy term has the largest statistical uncertainty and does not improve results in large tests, so it is often omitted; results are often best with a solute dielectric constant of 2 to 4; and the single-trajectory approach is usually more accurate than the three-trajectory approach. State in the paper whether entropy was computed and which dielectric constant was used.

The current gmx_MMPBSA version is 1.7.0. When copying older tutorials, note that extdiel in &gb has defaulted to 78.5 since 1.5.0 (80.0 before); with Interaction Entropy, report σIE, and treat results as unreliable when σIE exceeds about 3.6 kcal/mol. System setup, mdp parameters and commands are covered in /learn/gromacs-protein-ligand, and computing and plotting RMSD, RMSF, hydrogen bonds and SASA in /learn/md-trajectory-analysis.

• Run at least 3 independent simulations per complex (different initial-velocity random seeds) and report the mean and standard deviation.
• Show evidence that the analyzed properties have equilibrated, and state how the equilibration and analysis segments are divided.
• Run at least one simulation from a different starting conformation (for example the second-ranked docking pose) to show that conclusions do not depend on the starting structure.
• Provide a system setup table: box dimensions, total atoms, number of water molecules, salt concentration, force field and water model, protonation states.
• Report software versions and provide the initial coordinates, input parameter files and final coordinates.
• Simulate the positive-control complex with the same protocol.

Reviewer comments

Common reviewer questions and how to address them

The questions below come from journal review requirements and methodological review papers. The responses can be completed before submission.

Practice in China

Converting TCMSP target names to gene symbols and building Cytoscape networks: forum practice and a test

The target table of the old TCMSP gives protein names (the target_name column), which must be converted to gene symbols before intersection and enrichment. The usual approach in Chinese experience posts is to download human reviewed UniProt entries and match protein names to gene names (Excel VLOOKUP or an R script). Below are test results for this step, the cases that need manual handling, and common problems when importing into Cytoscape.

# Usage: python3 tcmsp2symbol.py tcmsp_targets.txt > target_symbol.tsv
# tcmsp_targets.txt: one TCMSP target name per line (target_name column, deduplicated)
# Python 3 standard library only; downloads human reviewed UniProt entries (about 20,000 rows) at run time
import csv, io, re, sys, urllib.request

URL = ("https://rest.uniprot.org/uniprotkb/stream?query=organism_id:9606+AND+reviewed:true"
       "&fields=accession,protein_name,gene_primary&format=tsv")
rec, alt = {}, {}
with urllib.request.urlopen(URL, timeout=600) as r:
    for row in csv.DictReader(io.StringIO(r.read().decode()), delimiter="\t"):
        gene = row["Gene Names (primary)"]
        if not gene:
            continue
        pn = row["Protein names"].split(" [Cleaved")[0]   # drop the cleaved-chain description
        rec.setdefault(re.sub(r"\s*\(.*$", "", pn).strip().lower(), set()).add(gene)  # recommended name
        for name in re.findall(r"\(([^()]*)\)", pn):      # alternative names in parentheses
            alt.setdefault(name.strip().lower(), set()).add(gene)

names = [l.strip() for l in open(sys.argv[1], encoding="utf-8") if l.strip()]
n_rec = n_alt = 0
print("target_name\tsymbols\tmatch")
for n in dict.fromkeys(names):
    key = n.lower()
    if key in rec:
        genes, how = rec[key], "recommended"; n_rec += 1
    elif key in alt:
        genes, how = alt[key], "alternative"; n_alt += 1   # alternative-name match: check each one by hand
    else:
        genes, how = set(), "none"                          # unmatched: look up in UniProt by hand
    flag = "ambiguous" if len(genes) > 1 else how
    print(f"{n}\t{'/'.join(sorted(genes))}\t{flag}")
print(f"targets={len(dict.fromkeys(names))} recommended={n_rec} alternative={n_alt} "
      f"unmatched={len(dict.fromkeys(names)) - n_rec - n_alt}", file=sys.stderr)

Hand off to an agent

Hand this workflow to an agent

The computational steps are fixed, produce many files and run for a long time, so they suit a science agent. Judgment calls and experiments remain yours.

Example one-sentence instruction: "Study Scutellaria baicalensis for ulcerative colitis: collect compounds and targets from HERB 2.0 and BATMAN-TCM 2.0, screen with OB ≥ 30%, DL ≥ 0.18 and the SwissTargetPrediction top 15; take disease genes from OMIM, DisGeNET curated and GEO differential genes; build a STRING 12.5 network and flag edges supported only by text mining; run GO/KEGG enrichment with the union of all compound targets and disease genes as background; do site-specific docking, redocking and positive-control docking for the top 5 core targets; run 3 replicas of 100 ns for each of the 3 best-scoring complexes and the positive control, compute MM/GBSA and make the figures."

1. 01

   Retrieval and records

   Queries each database, saves raw export tables and records versions and access dates.

2. 02

   Network and enrichment

   Runs Python and R scripts for the overlap, PPI, core targets and enrichment, and writes tables and figures.

3. 03

   Docking

   Downloads PDB structures, prepares receptors and ligands, runs redocking, compound docking and positive-control docking, and outputs the score table and interaction diagrams.

4. 04

   Molecular dynamics

   Builds systems in an environment with GROMACS, AmberTools and ACPYPE preinstalled, rents a GPU as needed for the replicas, and completes trajectory analysis and gmx_MMPBSA calculations.

5. 05

   Adversarial review

   Checks whether the task is actually complete, for example whether redocking RMSD is ≤ 2 Å, whether every complex has all its replicas, and whether figures show the mean and standard deviation.

• The workspace keeps: compound tables, target tables, network TSV files, enrichment results, docking logs and poses, MD input files and analysis results, and all scripts and figures.
• You still need to check: whether the core compounds are actually present in the herb and in what amount, whether the thresholds fit your study design, and whether the PDB structures and positive controls are appropriate.
• You design and run the experimental validation yourself.