数学推导与算法

本节推导幂律指数的最大似然估计并给出 Clauset 等（2009）的完整检验流程，写出网络基序检测的算法，求解一致型 FFL 与负自调控的 ODE，定量得到“符号敏感延迟”与“加速响应”；随后给出通量平衡分析的线性规划性质、MCL 算法、L3 链路预测的推导思路，以及差异相关检验所需的样本量。

幂律拟合：最大似然、KS 与 bootstrap

连续情形的 MLE。设度 $k\ge k_{\min}$ 的尾部服从 $p(k)=\dfrac{\alpha-1}{k_{\min}}\Big(\dfrac{k}{k_{\min}}\Big)^{-\alpha}$，尾部样本 $k_1,\dots,k_{n}$。对数似然

$$\ln\mathcal L(\alpha)=n\ln(\alpha-1)-n\ln k_{\min}-\alpha\sum_{i=1}^{n}\ln\frac{k_i}{k_{\min}}.$$

令 $\partial\ln\mathcal L/\partial\alpha=\dfrac{n}{\alpha-1}-\sum_i\ln\dfrac{k_i}{k_{\min}}=0$，得

$$\hat\alpha=1+n\Big[\sum_{i=1}^{n}\ln\frac{k_i}{k_{\min}}\Big]^{-1},\qquad \operatorname{SE}(\hat\alpha)\approx\frac{\hat\alpha-1}{\sqrt n}.$$

度是整数，严格应使用离散幂律（归一化常数为 Hurwitz zeta 函数）数值求解；Clauset 等给出的近似是把 $k_{\min}$ 换成 $k_{\min}-\tfrac12$，在 $k_{\min}$ 不太小时误差很小。

为什么不用双对数回归？ 对直方图在双对数坐标下做最小二乘，误差不满足回归假设（尾部计数少、方差大，分箱方式随意），得到的斜率有系统偏差，$R^2$ 也不能作为幂律的证据——对数正态分布同样能给出很高的 $R^2$。

算法 1  幂律检验（Clauset, Shalizi & Newman 2009）
输入：度序列 {k_i}；bootstrap 次数 B（p 值精确到 ±ε 约需 B ≈ 1/(4ε²)，如 ε=0.01 时 2500）
输出：k_min, α̂, 拟合优度 p 值，与替代分布的似然比检验结果
1. 对每个候选 k_min（度序列中的每个不同取值）：
2.     用尾部 {k_i ≥ k_min} 求 α̂(k_min)（MLE）
3.     D(k_min) ← max_k |F_emp(k) − F_fit(k)|      // 尾部经验 CDF 与拟合 CDF 的 KS 距离
4. k_min* ← argmin D；记 D* 与尾部样本量 n_tail
5. 对 b = 1..B（半参数 bootstrap）：
6.     对每个样本位置：以概率 n_tail/n 从拟合的幂律（k ≥ k_min*）抽值，
       否则从原数据中 k < k_min* 的部分有放回抽值
7.     对合成数据重复步骤 1–4（包括重新选 k_min），得 D_b
8. p ← #{D_b ≥ D*} / B；p ≤ 0.1 时拒绝幂律
9. 与对数正态、指数、截断幂律分别做对数似然比 R 与 Vuong 检验：
       R > 0 且显著 → 偏向幂律；R < 0 且显著 → 偏向替代分布；不显著 → 无法区分
复杂度：步骤 1–4 为 O(U·n)（U 为不同度值数）；总计 O(B·U·n)

两点要强调：$p>0.1$ 只说明“不能排除幂律”，并不证明幂律；对 PPI 这样只有几千个节点、尾部只有几十个点的网络，检验功效往往很低，最常见的结论是“与对数正态无法区分”。

网络基序检测

算法 2  三节点基序的富集检验
输入：有向网络 G；随机网络数 R（如 1000）；子图大小 3
输出：每类三节点连通子图的 Z 分数与经验 p 值
1. 枚举 G 中全部连通三节点诱导子图，按同构类别计数 N_real(c)，13 类
   （可用 ESU 等算法；稀疏图上约为 O(m · k_max)，m 为边数，k_max 为最大度）
2. 对 r = 1..R：
3.     用有向边交换生成随机网络 G_r：交换 (a→b),(c→d) 为 (a→d),(c→b)，
       保持每个节点的入度与出度；互惠边（a↔b）应作为整体单独交换以保持其数量
4.     计数 N_r(c)
5. Z(c) ← (N_real(c) − mean_r N_r(c)) / sd_r N_r(c)
6. p(c) ← #{N_r(c) ≥ N_real(c)} / R
7. （可选）显著性谱 SP(c) ← Z(c) / sqrt(Σ_c' Z(c')²)，用于跨网络比较
复杂度：O(R · (Q·m + 枚举开销))，Q 为每边交换次数

注意：（1）当 $\sigma_{\text{rand}}$ 很小时 $Z$ 会极大，分布也未必近似正态，应同时看经验 $p$ 值和计数比 $N_{\text{real}}/\langle N_{\text{rand}}\rangle$；（2）检验四节点基序时，零模型最好也保持三节点子图的计数，否则四节点“富集”可能只是三节点结构的副产品（Milo 等, 2002 讨论了这一点）。

一致型 1 型 FFL 的延迟：解析解

沿用前馈环中的模型：$X$ 在 $t=0$ 时被激活且 $X>K_{XY}$、$X>K_{XZ}$，初始 $Y(0)=Z(0)=0$。

Y 的动力学：$\dot Y=\beta_Y-\alpha Y$，解为 $Y(t)=Y_{\text{st}}\big(1-e^{-\alpha t}\big)$，稳态 $Y_{\text{st}}=\beta_Y/\alpha$。

Z 何时开始表达：“与”门要求 $Y>K_{YZ}$。令 $Y(T_{\text{on}})=K_{YZ}$，得

$$T_{\text{on}}=\frac1\alpha\ln\frac{1}{1-K_{YZ}/Y_{\text{st}}}.$$

$T_{\text{on}}$ 只在 $K_{YZ}<Y_{\text{st}}$ 时有限；$K_{YZ}$ 越接近 $Y_{\text{st}}$，延迟越长。之后 $Z(t)=\dfrac{\beta_Z}{\alpha}\big(1-e^{-\alpha(t-T_{\text{on}})}\big)$，$t>T_{\text{on}}$。

关闭时没有延迟：若 $X$ 在某时刻失活，$\theta(X>K_{XZ})$ 立即变为 0，$Z$ 的生成项立即消失，$Z$ 按 $e^{-\alpha t}$ 衰减，与 $Y$ 的残余量无关。于是“开”有延迟、“关”无延迟，称为符号敏感延迟。若把逻辑门换成“或”门，结论反过来：开启立即发生，关闭时 $Z$ 要等 $Y$ 衰减到 $K_{YZ}$ 以下才停止。

持续性检测：一个持续时间 $\Delta<T_{\text{on}}$ 的输入脉冲期间 $Y$ 达不到阈值，$Z$ 完全不响应。这给出了滤波截止的定量判据。注意上述推导依赖阶跃函数和相同的降解速率 $\alpha$；用 Hill 函数替换阶跃函数后，延迟变得平滑，但定性结论保持。

负自调控的上升时间

简单调控：$\dot X=\beta-\alpha X$，$X(0)=0$，解为 $X(t)=\frac{\beta}{\alpha}(1-e^{-\alpha t})$。令 $X=\frac{\beta}{2\alpha}$ 得 $T_{1/2}=\ln2/\alpha$，与 $\beta$ 无关：要想更快，只能提高降解速率，而这要以持续“边造边拆”的能量浪费为代价。

负自调控（阶跃近似）：$\dot X=\beta\,\theta(X<K)-\alpha X$。当 $X<K$ 时与简单调控相同，$X(t)=\frac{\beta}{\alpha}(1-e^{-\alpha t})$；当 $X$ 到达 $K$ 时生成停止，稳态为 $X_{\text{st}}=K$（要求 $\beta/\alpha>K$）。上升时间满足 $\frac{\beta}{\alpha}(1-e^{-\alpha T})=\frac K2$：

$$T_{1/2}^{\text{NAR}}=\frac1\alpha\ln\frac{1}{1-\alpha K/(2\beta)}\approx\frac{K}{2\beta}\quad(\beta\gg\alpha K).$$

与稳态同为 $K$ 的简单调控（$\beta'=\alpha K$，$T_{1/2}=\ln2/\alpha$）相比，比值为 $\frac{\alpha K}{2\beta\ln2}$——强启动子（大 $\beta$）配合负自调控，就能在不改变稳态的前提下大幅加快响应。用 Hill 函数 $\beta K^n/(K^n+X^n)$ 替代阶跃函数时，加速随 $n$ 增大而增强；代码示例中 $K=0.5$、$n=2$ 时 $T_{1/2}\approx0.15$，而简单调控为 0.69。

鲁棒性：阶跃近似下 $X_{\text{st}}=K$ 与 $\beta$ 无关，因此产生速率的细胞间波动（如质粒拷贝数、RNA 聚合酶量的差异）不会传到稳态水平上——这就是负自调控降低细胞间差异的定性原因。

通量平衡分析的线性规划性质

FBA 是标准线性规划：可行域 $\{v: Sv=0,\ l\le v\le u\}$ 是凸多面体，线性目标的最优值必在某个顶点上取得，可用单纯形法或内点法在多项式时间内求解，基因组尺度模型（数千个反应）几秒内即可求解。其对偶问题为每个代谢物给出一个影子价格 $y_i$（对偶变量）：它表示放松第 $i$ 个稳态约束（额外提供一个单位的代谢物 $i$）时目标值的变化率，可以用来判断哪种营养或前体限制了生长。最优解往往不唯一（多个顶点或一整条棱都最优），这就是需要通量变异性分析的原因：对每个反应 $j$ 求解 $$\min/\max\; v_j\quad\text{s.t.}\quad Sv=0,\ l\le v\le u,\ c^\top v\ge \gamma\, Z^\ast ,$$ 其中 $Z^\ast$ 为最优目标值，$\gamma\in(0,1]$ 允许一定的次优。

MCL 算法

算法 3  马尔可夫聚类（Enright 等 2002）
输入：无向（可加权）网络的邻接矩阵 A；膨胀参数 r > 1（常用 1.4–5）
输出：节点的簇划分
1. M ← A + I；把每一列归一化为和 1（列随机矩阵）
2. 重复直到 M 收敛：
3.     扩张：M ← M · M                      // 两步随机游走
4.     膨胀：M_ij ← M_ij^r；再按列归一化     // 强化强流、削弱弱流
5.     （可选）剪枝：把很小的元素置 0 以保持稀疏
6. 收敛后 M 的非零行称为“吸引子”；每一列 j 归入其值最大的吸引子行，同一吸引子的节点构成一个簇
复杂度：稠密实现每次迭代 O(N³)；带剪枝的稀疏实现约为 O(N·k²)，k 为每列保留的非零元数

扩张使流在同一稠密区域内均匀化，膨胀让“簇内”高概率的流越来越强、“簇间”低概率的流逐渐消失。$r$ 越大，膨胀越强，簇越小。MCL 每个节点只属于一个簇（不允许重叠），且簇的数量由数据与 $r$ 决定，无需预先指定。

L3 链路预测的推导思路

设两个蛋白 $X$、$Y$ 的结合依赖互补界面。若 $X$ 与 $V$ 相似（共享伙伴），而 $V$ 与 $Y$ 结合，则 $X$ 也可能与 $Y$ 结合。“$X$ 与 $V$ 共享伙伴 $U$”对应路径 $X$–$U$–$V$，再加上 $V$–$Y$，就是长度为 3 的路径 $X$–$U$–$V$–$Y$。对所有中间节点求和即 $(A^3)_{XY}$。高度节点会贡献大量路径而淹没信号，因此 Kovács 等（2019）按 $1/\sqrt{k_Uk_V}$ 加权，这等价于对矩阵 $D^{-1/2}AD^{-1/2}$ 的中间部分归一化：$p=A\,D^{-1/2}AD^{-1/2}\,A$。相比之下，共同邻居打分 $(A^2)_{XY}$ 隐含的是“相似者相连”（同配）假设，在互补结合主导的蛋白网络中并不成立。

差异相关检验的样本量

由 Fisher 变换，检测 $\rho_A$ 与 $\rho_B$ 之差所需的每组样本量（双侧水平 $\alpha_0$、功效 $1-\beta_0$）约为

$$n\approx 3+\frac{2\,(z_{1-\alpha_0/2}+z_{1-\beta_0})^2}{(\operatorname{atanh}\rho_A-\operatorname{atanh}\rho_B)^2}.$$

算例：$\rho_A=0.5$、$\rho_B=0.2$，$\operatorname{atanh}$ 之差为 $0.549-0.203=0.346$；$\alpha_0=0.05$、功效 0.8 时 $(1.96+0.84)^2=7.84$，$n\approx3+2\times7.84/0.120\approx134$。若再对 1000 个基因对做 Bonferroni 校正（$\alpha_0=5\times10^{-5}$，$z\approx4.06$），$(4.06+0.84)^2=24.0$，$n\approx3+2\times24.0/0.120\approx402$。这说明中等程度的相关变化需要数百个样本才能可靠地在多重比较下检出。