模块化"脑网络"的虚拟损伤——枢纽 vs 连接者

# 生成一个 6 模块的随机加权网络（模块内连接密、模块间稀疏），按度（模块内枢纽）与参与系数（连接者）选损伤节点，
# 删除节点后比较全局效率与模块度（用已知的模块划分计算）的下降，并与随机删除比较。
# 这只是 Alstott 等（2009）、Gratton 等（2012）思路的拓扑简化：真实研究用结构连接组和动力学模型或病人数据。
# 解读：删除连接者对模块间整合（效率）的损害通常大于同等数量的随机节点；模块度可能反而升高（模块更"孤立"）。
import numpy as np
import networkx as nx

rng = np.random.default_rng(2)
sizes = [20] * 6
P = np.full((6, 6), 0.01); np.fill_diagonal(P, 0.35)
G = nx.stochastic_block_model(sizes, P.tolist(), seed=3)
comm = {v: G.nodes[v]["block"] for v in G}
# 指定少数"连接者"：给 6 个节点添加跨模块连接
connectors = [i * 20 for i in range(6)]
for c in connectors:
    for b in range(6):
        if b != comm[c]:
            for w in rng.choice(range(b * 20, b * 20 + 20), 2, replace=False):
                G.add_edge(c, int(w))

def participation(G, comm):
    pc = {}
    for v in G:
        k = G.degree(v)
        if k == 0:
            pc[v] = 0; continue
        counts = np.bincount([comm[w] for w in G[v]], minlength=6)
        pc[v] = 1 - np.sum((counts / k) ** 2)
    return pc

def lesion_effect(G, targets):
    H = G.copy(); H.remove_nodes_from(targets)
    parts = [{v for v in H if comm[v] == b} for b in range(6)]
    parts = [p for p in parts if p]
    return nx.global_efficiency(H), nx.community.modularity(H, parts)

pc = participation(G, comm); deg = dict(G.degree())
within = [v for v in G if v not in connectors]
hubs = sorted(within, key=lambda v: -deg[v])[:6]
conn = sorted(G, key=lambda v: -pc[v])[:6]
E0, Q0 = lesion_effect(G, [])
print("完整网络: 全局效率 %.3f  模块度 %.3f" % (E0, Q0))
eff = [(E0, Q0)]
for name, t in [("模块内枢纽(高度、低参与系数)", hubs), ("连接者(高参与系数)", conn)]:
    E, Q = lesion_effect(G, t); eff.append((E, Q)); print("%s: 效率 %.3f  模块度 %.3f" % (name, E, Q))
rand = [lesion_effect(G, list(rng.choice(G.number_of_nodes(), 6, replace=False))) for _ in range(30)]
print("随机 6 个节点: 效率 %.3f  模块度 %.3f" % tuple(np.mean(rand, axis=0)))

import matplotlib.pyplot as plt
eff = np.vstack([eff, np.mean(rand, axis=0)]); err = [0, 0, 0, np.std([r[0] for r in rand])]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
pos = nx.spring_layout(G, seed=1)
col = ["tab:red" if v in conn else ("tab:orange" if v in hubs else "0.75") for v in G]
nx.draw_networkx_edges(G, pos, ax=ax1, edge_color="0.6", alpha=0.3, width=0.5)
nx.draw_networkx_nodes(G, pos, ax=ax1, node_size=[110 if v in conn or v in hubs else 40 for v in G], node_color=col)
ax1.set_title("Connectors (red) vs within-module hubs (orange)"); ax1.axis("off")
ax2.bar(["intact", "hubs", "connectors", "random"], eff[:, 0], color=["0.4", "tab:orange", "tab:red", "0.75"], width=0.6)
ax2.errorbar(3, eff[3, 0], yerr=err[3], c="k", capsize=4, label="random: mean +/- SD (30 draws)")
ax2.set_ylim(0, E0 * 1.3); ax2.set_ylabel("global efficiency"); ax2.legend(fontsize=8, loc="upper right")
ax2.set_title("Removing 6 connectors hurts integration most")
plt.tight_layout(); plt.show()