跳到主要内容
Sandbox Physics
演化统计物理

L4 · 第 7 章

读《既非拉动也非推动》并复核其中的漂变结论

Sandbox Physics ·

Birzu、Hallatschek 和 Korolev 认为，由增长和扩散形成的波有三类而不是两类，而中间那一类保住的遗传多样性，远少于推动波应有的水平。本章是一份对照 Lab 阅读这篇论文的导读：每个结论在哪里、靠哪个量支撑，哪些数我们复现了、哪些没有，以及原因。

一段话说清结论

拉动波里，尖端少数个体控制整道波，所以涨落随每格个体数 NN 只按对数减小。直觉上，推动波应像任何大种群一样把涨落平均掉，涨落按 1/N1/N 减小。论文发现，这只对完全推动波成立。在 1<v/vF<3/(22)1<v/v_{\mathrm F}<3/(2\sqrt2) 的半推动波里，前沿游走和多样性丢失速率都按 NN 的幂次减小，指数介于 00 和 −1-1 之间。

证据分两部分：在不同协作程度和种群密度下对踏脚石模型的模拟，以及一套解释 1/N1/N 标度在哪里失效的微扰理论。所谓普适性，即分界只取决于速度比，依据的是三种不同的增长律落在同一条曲线上。

来源：G. Birzu

先看哪里

把式 (3) 和图 1 一起读：带协作程度 BB 的增长律，以及 B=2B=2 之后速度超过 Fisher 速度。方法部分的式 (11) 给出模拟：每格先扩散、再增长、再做多项式重抽样。Lab 运行的就是这个模型。

图 5 承载关于多样性的结论。A 图是杂合度 HH 的指数衰减，B 图是衰减率 Λ\Lambda 在对数坐标下随 NN 的变化，C 图是拟合指数 αh\alpha_h 随 BB 的变化。式 (6) 和 (7) 给出两类推动波的分界。推导和闭式解在补充材料里。

H(t)=⟨2f(1−f)⟩∼e−Λt,Λ∝NαhH(t)=\bigl\langle 2f(1-f)\bigr\rangle\sim e^{-\Lambda t},\qquad \Lambda\propto N^{\alpha_h}
测量对象：前沿的平均杂合度和它的衰减率，再由衰减率得出给波分类的指数。

来源：G. Birzu

衰减率是怎样测的

补充材料里关于模拟时长的那一节很短，值得细读。杂合度先在模拟区域各格上平均，再对 10001000 次独立运行平均。对 ln⁡H\ln H 随时间的拟合，终点取在仍有至少 5050 次运行保有两种等位基因的最后时刻；起点则选在拟合优度最高处，且窗口至少 10001000 代。指数只用 N>104N>10^4 的数据拟合，因为幂律只是渐近成立。

对读者来说有两点要紧。早期非指数的衰减部分是有意排除的。在有限的 NN 范围内拟合出的指数，也未必等于渐近值。图 5B 中拉动波的有效指数约为 −0.33-0.33；理论给出对数律，对应指数为零。对前沿游走，van Saarloos 的综述记录了同一类规律：前沿扩散常数按 1/ln⁡3N1/\ln^3N 减小，并指出 Brunet 和 Derrida 要借助特殊方法才能把 NN 推到足以看清它的程度。

来源：G. Birzu · W. van Saarloos

读懂式 (5)

式 (5) 是多样性丢失速率的微扰结果。前因子 1/N1/N 就是中心极限定理的标度。积分之比给前沿上每个位置加权：分母统计祖先，权重与第 6 章的 ρ2evζ/D\rho^2e^{v\zeta/D} 相同；分子统计漂变在哪里起作用。

尖端 ρ∼e−kζ\rho\sim e^{-k\zeta}，所以分子的被积函数按 e(2v/D−3k)ζe^{(2v/D-3k)\zeta} 变化。只有 2v/D<3k2v/D<3k 时它才收敛，也就是 v/vFv/v_{\mathrm F} 超过 3/(22)3/(2\sqrt2) 时。低于这一点，积分发散，需要在低密度处截断，而截断带来的 NN 的幂次不再是 −1-1。对 B>4B>4 的协作增长律，补充材料式 (S122) 给出了积分的闭式。

Λ=1N ∫−∞∞γf ρ3(ζ) e2vζ/D dζ[∫−∞∞ρ2(ζ) evζ/D dζ]2\Lambda=\frac{1}{N}\,\frac{\displaystyle\int_{-\infty}^{\infty}\gamma_f\,\rho^3(\zeta)\,e^{2v\zeta/D}\,d\zeta}{\Bigl[\displaystyle\int_{-\infty}^{\infty}\rho^2(\zeta)\,e^{v\zeta/D}\,d\zeta\Bigr]^2}
论文式 (5)：多样性丢失速率，ρ\rho 是随前沿运动坐标中的归一化密度分布。
Λ=γf4πNg0B2D (B+4)tan⁡2πB(B>4)\Lambda=\frac{\gamma_f}{4\pi N}\sqrt{\frac{g_0B}{2D}}\,(B+4)\tan\frac{2\pi}{B}\qquad(B>4)
补充材料式 (S122)：协作增长律下的同一速率，适用于完全推动波。

来源：G. Birzu

完整例题：复现一个衰减率

题目。取 B=10B=10、g0=0.01g_0=0.01、D=0.125D=0.125、γf=1\gamma_f=1、N=300N=300。预测 Λ\Lambda，并与 Lab 模型比较。

闭式。0.01×10/0.25=0.632\sqrt{0.01\times10/0.25}=0.632，B+4=14B+4=14，tan⁡(2π/10)=0.727\tan(2\pi/10)=0.727。乘积为 6.436.43，除以 4π×300=37704\pi\times300=3770，得 Λ=1.71×10−3\Lambda=1.71\times10^{-3}（每代）。改用格点前沿的分布计算式 (5)，得 1.70×10−31.70\times10^{-3}。

模拟。320320 次固定种子的运行，先预热 15001500 代，再把两种等位基因按 5050:5050 重新混合，杂合度在前沿后方十格到前方三格之间平均。在第 10001000 到 21002100 代之间对 ln⁡H\ln H 做直线拟合，得 Λ=(1.68±0.20)×10−3\Lambda=(1.68\pm0.20)\times10^{-3}，误差用对运行的刀切法估计。与微扰值之比为 0.990.99。

第二个点。N=1000N=1000 时测得 (4.3±0.7)×10−4(4.3\pm0.7)\times10^{-4}，预测 5.1×10−45.1\times10^{-4}，低 1.11.1 个标准误。两个规模之间速率下降了 3.9±0.83.9\pm0.8 倍，严格的 1/N1/N 律给出 3.333.33。

结论。对完全推动波，Lab 模型在几百次运行的误差范围内，复现了 Λ\Lambda 的大小和它的 1/N1/N 标度。

来源：G. Birzu

这次复现没有说明什么

拉动波在 Lab 中测得：N=100N=100 时 Λ=(1.96±0.21)×10−3\Lambda=(1.96\pm0.21)\times10^{-3}，N=1000N=1000 时 (0.92±0.11)×10−3(0.92\pm0.11)\times10^{-3}。比值 2.1±0.32.1\pm0.3 远小于 1/N1/N 律的 1010 倍，也低于渐近 ln⁡3N\ln^3N 律的 3.43.4。换算成有效指数是 −0.33±0.07-0.33\pm0.07，与图 5B 的值接近；不过我们的 NN 比论文的低，这种一致可能部分是巧合。

Lab 的做法与论文有几处不同：运行几百次而不是一千次；拟合窗口截止于剩 3030 次存活运行而不是 5050 次；杂合度只在前沿附近十四格上平均，而不是整个区域；NN 远低于 10410^4。我们没有测半推动波的指数、前沿扩散常数，也没有检验图 6 中三种增长律的重合。这些才是论文的核心新结果，要检验它们需要规模大得多的运行。Birzu、Matin、Hallatschek 和 Korolev 的后续工作报告，扩散随密度增加时，前沿的漂变会被压低，所以波属于哪一类，既取决于个体怎样增长，也取决于它们怎样移动。Lab 中的扩散是固定的。

来源：G. Birzu · G. Birzu

到实验里试一试

看中心极限标度在哪里成立，又在哪里失效

  1. 选“完全推动波：后方也有祖先”，看多样性读数下降。预测容纳量从一千升到约三千时，它会下降得慢多少。
  2. 做这个改动并比较。再选“拉动波：尖端胜出”，做同样的改动，预测变慢的幅度是否一样大。
  3. 把协作程度设为三。用前沿尖端的情况解释，为什么论文要靠截断来估计这种波中的漂变。

练习与解答

  1. 用补充材料式 (S122) 和 Lab 参数，算 B=10B=10、N=1000N=1000 时的 Λ\Lambda。

    展开解答

    速率按 1/N1/N 标度，所以是 1.706×10−3×300/1000=5.12×10−41.706\times10^{-3}\times300/1000=5.12\times10^{-4}（每代）。Lab 测得 (4.3±0.7)×10−4(4.3\pm0.7)\times10^{-4}。

  2. 补充材料给出拉动前沿的过渡为 v(t)=vF(1−3/(4r0t))v(t)=v_{\mathrm F}\bigl(1-3/(4r_0t)\bigr)。证明对连续方程它与第 6 章的 v∗−3/(2λ∗t)v^*-3/(2\lambda^*t) 相同。

    展开解答

    连续方程中 v∗=vF=2Dr0v^*=v_{\mathrm F}=2\sqrt{Dr_0}，λ∗=r0/D\lambda^*=\sqrt{r_0/D}。于是 vF⋅3/(4r0t)=3D/(2r0 t)=3/(2λ∗t)v_{\mathrm F}\cdot3/(4r_0t)=3\sqrt{D}/(2\sqrt{r_0}\,t)=3/(2\lambda^*t)，两式一致。

  3. NN 增大十倍时某速率下降为 1/2.121/2.12。对应的有效指数 αh\alpha_h 是多少？在 N=100N=100 到 10001000 之间，ln⁡3N\ln^3N 律会给出多少？

    展开解答

    αh=−log⁡102.12=−0.33\alpha_h=-\log_{10}2.12=-0.33。ln⁡3N\ln^3N 律给出比值 (ln⁡1000/ln⁡100)3=1.53=3.4(\ln1000/\ln100)^3=1.5^3=3.4，有效指数 −0.53-0.53，尽管它的渐近指数是零。只跨一个数量级拟合出的有效指数，几乎说明不了渐近规律。

来源与复算范围

  1. G. Birzu, O. Hallatschek and K. S. Korolev, Neither pulled nor pushed: genetic drift and front wandering uncover a new class of reaction-diffusion waves, arXiv:1709.01601v2 (2018)

    第 2 版，全文。式 (2)：Fisher 速度。式 (3)：协作增长律。式 (5)：多样性丢失速率的微扰结果。式 (6)—(7)：前沿尾部的衰减，以及区分半推动波与完全推动波的速度比 3/(22)3/(2\sqrt2)。图 5：三类波的多样性丢失。方法部分式 (11)：模拟模型。补充材料式 (S28)：祖先权重；式 (S122)：丢失速率的闭式；“模拟时长与数据采集”一节，式 (S149)—(S150)。

  2. G. Birzu, S. Matin, O. Hallatschek and K. S. Korolev, Genetic drift in range expansions is very sensitive to density feedback in dispersal and growth, arXiv:1903.11627v1 (2019)

    第 1 版，全文。摘要：扩散随种群密度增加时，前沿的遗传漂变被大幅压低；决定多样性丢失属于哪种情形的是一个比值，即扩张速度与前沿处扩散率和增长率几何平均之比。

  3. W. van Saarloos, Front propagation into unstable states, arXiv:cond-mat/0308540v2 (2003), Physics Reports

    第 2 版，全文。第 2.2 节，式 (24)—(27)：由色散关系求线性扩展速度 v∗v^* 与陡度 λ∗\lambda^*。第 2.9.1 节，式 (78)：拉动前沿普适的缓慢逼近 v(t)=v∗−3/(2λ∗t)v(t)=v^*-3/(2\lambda^*t)。第 7.1 节，式 (241)—(246)：粒子模型中拉动速度的截断修正，正比于 1/ln⁡2N1/\ln^2N；式 (248)：前沿扩散常数按 1/ln⁡3N1/\ln^3N 减小。

反馈这个页面有错误、看不明白，或者缺了什么？告诉我们
反馈类型

提交时会附带本页地址、Lab ID、网站语言和设备类型。请不要填写个人信息。反馈最多保存 12 个月。隐私政策