張 馳

(1 中國科學(xué)院古脊椎動物與古人類研究所，中國科學(xué)院脊椎動物演化與人類起源重點實驗室 北京 100044)

(2 中國科學(xué)院生物演化與環(huán)境卓越創(chuàng)新中心 北京 100044)

1 簡介

推斷類群的系統(tǒng)關(guān)系和分異時間是分支系統(tǒng)學(xué)分析的基礎(chǔ)。如何合理利用化石的形態(tài)和地質(zhì)年代的數(shù)據(jù)來完成此類推斷一直是一個棘手的問題。傳統(tǒng)的推斷過程往往采用分步的策略。首先單獨從形態(tài)數(shù)據(jù)出發(fā)，利用簡約法構(gòu)建類群的系統(tǒng)關(guān)系。這個關(guān)系只有拓撲結(jié)構(gòu)而沒有時間的信息，而且只代表給定時間和搜索條件下最優(yōu)的結(jié)果(即最簡約的樹或其合意樹)。然后固定這個拓撲結(jié)構(gòu)，從化石年代出發(fā)，利用最小枝長法(Laurin,2004)或等距枝長法(Brusatte et al., 2008)確定樹中內(nèi)部節(jié)點的分異時間。最小枝長法是把當(dāng)前節(jié)點的時間往前推一個百萬年作為祖先節(jié)點的時間；而等距枝長法則是把祖先節(jié)點到后代節(jié)點的時間的中點最為當(dāng)前節(jié)點的時間。至于演化速率，可以通過祖先狀態(tài)重建和每個樹枝上特征變化的次數(shù)，結(jié)合上一步推斷的時間來估計。這種分步計算的方式比較直觀，因此被用于實際數(shù)據(jù)分析(Wang and Lloyd, 2016)。然而，該策略存在諸多缺陷。首先，它每一步都忽略了推斷中的不確定性，包括樹的拓撲結(jié)構(gòu)、分異時間和祖先特征狀態(tài)；其次，每一步都只利用了一部分數(shù)據(jù)信息，如建樹時只用到形態(tài)特征，定年時只用到化石年代；再次，定年的方式很主觀，對化石數(shù)量的增減很敏感，且不適用于現(xiàn)生類群；最后，整個過程缺乏一個嚴謹?shù)慕y(tǒng)計學(xué)框架，無法對不同的模型假設(shè)進行檢驗。

近些年開發(fā)的貝葉斯支端定年法(Bayesian tip dating) (Ronquist et al., 2012;Gavryushkina et al., 2014; Zhang et al., 2016)很好地克服了上述問題。貝葉斯支端定年法把化石形態(tài)和年代數(shù)據(jù)整合在一次完整的計算過程中，能夠盡可能地利用數(shù)據(jù)信息，同時考慮了樹的拓撲結(jié)構(gòu)、分異時間、演化速率以及化石年代的不確定性。該方法通過統(tǒng)計模型來描述特征的演化、類群的生滅以及化石的采樣等過程，并借助相對成熟的貝葉斯統(tǒng)計框架和計算方法來進行參數(shù)估計和模型選擇。但該方法相對來說比較復(fù)雜，需要較多的統(tǒng)計學(xué)知識，對古生物學(xué)家來說往往難于理解和上手，而且鮮有系統(tǒng)闡述支端定年法計算過程和參數(shù)意義的文獻(Gavryushkina and Zhang, 2020)。本文逐層剖析支端定年法的計算過程，解釋其中用到的重要模型和參數(shù)意義，旨在一定程度上為古生物學(xué)家分析實際數(shù)據(jù)提供參考。

本文首先介紹描述時間樹的石化生滅過程(fossilized birth-death process)模型(Stadler,2010), 然后介紹描述特征演化速率的寬松形態(tài)鐘(relaxed clock)模型，接著介紹描述特征狀態(tài)變化的Mk模型(Lewis, 2001), 再通過貝葉斯公式把上述模型聯(lián)系起來，最后介紹估計參數(shù)后驗分布的馬氏鏈蒙特卡羅(Markov chain Monte Carlo, MCMC)算法。附錄提供了計算中生代鳥類數(shù)據(jù)(Zhang and Wang, 2019)的MrBayes命令。

2 分異時間

時間樹(timetree)代表類群的系統(tǒng)關(guān)系和分異時間，它的概率分布可以通過石化生滅過程(Stadler, 2010)來給出。該過程描述了從這些類群的最近共同祖先(樹根)開始，分異、滅絕、采集化石和采樣現(xiàn)生類群這一系列事件的發(fā)生，并對應(yīng)于一棵完整樹(圖1A)。但是實際數(shù)據(jù)分析中無法推斷這個完整樹，只能推斷和樣本相關(guān)的部分，即樣本樹(圖1B)。記每個樹枝的分異速率(或叫成種率)為λ, 滅絕速率為μ, 沿每個樹枝的化石采樣速率為ψ, 現(xiàn)生類群的采樣概率(或采樣比例)為ρ, 通過建立并求解一系列常微分方程，可以得到給定λ,μ,ψ,ρ時，樣本時間樹T= {τ, t}的概率分布，記為P(T|λ,μ,ψ,ρ), 其中τ代表拓撲結(jié)構(gòu)，t代表以百萬年為單位的分異時間。

圖1 石化生滅過程可能得到的時間樹Fig. 1 Example timetree generated from the fossilized birth-death process

在MrBayes軟件中，該生滅過程以樹根時間t1為起始條件(圖2), 計算時要指定t1的先驗分布。通常這個先驗比較寬泛(最大范圍從0到無窮), 不過一般能從研究的類群預(yù)估一個更精確的范圍，例如，其下界不會早于最古老化石的年代?；臅r間可以固定為具體的數(shù)值(百萬年前), 也可以用一個均勻分布給出時間的上下界。在計算時還需要提供現(xiàn)生類群大致的采樣比例(ρ)?，F(xiàn)生類群可以有兩種采樣策略，一種是均勻隨機采樣(random), 另一種是多樣化采樣(diversity) (Zhang et al., 2016), 可以根據(jù)實際數(shù)據(jù)的情況自行選擇，后者可能更符合高階元類群的采樣模式(如每個科只取一個代表的屬或者每個屬取一個代表物種)。對于分異、滅絕和化石采樣速率，程序為了設(shè)置先驗方便，重新參數(shù)化為d=λ-μ,v=μ/λ,s=ψ/ (μ+ψ)。d的默認先驗為指數(shù)分布(范圍從0到無窮),v和s的默認先驗為均勻分布(范圍從0到1,更一般地為貝塔分布)。這樣，時間樹包括分異時間等參數(shù)的先驗分布基本就確定了。

圖2 用于各個概率分布公式中參數(shù)的示例樹上化石F1和F2的特征狀態(tài)為0,現(xiàn)生類群S1和S2的特征狀態(tài)為1,內(nèi)部節(jié)點的特征狀態(tài)用x0, x1, x2表示化石的年代分別為t3 = 100 Ma, t5 = 50 Ma,樹根的時間為t1, 其他分異時間為t2和t4記t = {t1, t2, t3, t4, t5}。各個枝上的特征演化速率為r = {r1, r2, r3, r4, r5}Fig. 2 Example parameters and symbols used in the probability distributions The character states for fossils F1 and F2 are 0, for extant taxa S1 and S2 are 1, for internal nodes are x0, x1, x2. The ages of fossils F1 and F2 are t3 = 100 Ma and t5 = 50 Ma The root age is t1 and the remaining divergence times are t2 and t4. Denote t = {t1,t2, t3, t4, t5}. The evolutionary rates on the branches are r = {r1, r2, r3, r4, r5}

需特別提到的是，有些數(shù)據(jù)只包含了化石而沒有現(xiàn)生類群。這時一般假設(shè)生滅過程在未到達現(xiàn)今時間點所有類群就都滅絕了。因此現(xiàn)生類群不論采用何種采樣策略和比例(程序默認為1.0), 都沒有樣本被采集到。不過由于MrBayes軟件的限制，最年輕的樣本總被顯示為現(xiàn)生類群(時間為0), 因此需要把整棵樹的時間軸進行相應(yīng)的平移，或者說分異時間要加上年代最近的化石的年代。這只是結(jié)果顯示的問題而不是程序計算的錯誤(后續(xù)版本將修復(fù)這個問題)。

除了石化生滅過程之外，MrBayes還為時間樹提供了一個均勻分布的先驗(Ronquist et al., 2012)。它沒有生滅和采樣這些參數(shù)，只依賴于樹根時間t1, 因此只需設(shè)置t1的先驗分布。均勻分布時常被認為是沒有信息的先驗，但它實際上往往帶有很強的信息，在定年這個問題上，會影響對分異時間的估計。而石化生滅過程看似參數(shù)很多，但實際上其設(shè)置可以很靈活。例如，分異、滅絕和化石采樣速率都可以隨時間變化，在不同的時間段內(nèi)各自獨立(Gavryushkina et al., 2014; Zhang et al., 2016)。這可能更符合實際生物學(xué)過程，同時還能推測凈成種速率和化石采樣速率隨時間的變化。

3 演化速率

形態(tài)特征的演化速率是指每百萬年每個特征期望的變化次數(shù)。對于給定的一段時間，演化速率越快，則特征最終期望的變化次數(shù)越多。一般給每個樹枝一個演化速率參數(shù)，記為r (圖2)。時鐘模型應(yīng)用于形態(tài)數(shù)據(jù)時被稱為形態(tài)鐘模型，類比于用在分子數(shù)據(jù)時的分子鐘模型。嚴格鐘(strict clock)模型假設(shè)演化速率在各個樹枝上都相同，這通常不適用于形態(tài)數(shù)據(jù)，實際分析時往往需使用寬松鐘(relaxed clock)模型。寬松鐘模型可以分為兩類，一類為獨立速率，另一類為自相關(guān)(autocorrelated)速率，區(qū)別在于r的概率分布P(r)不同。

獨立速率模型假設(shè)枝上的演化速率彼此獨立，它們都服從均值相同的某個概率分布。常用的概率分布包括伽馬分布(Lepage et al., 2007)和對數(shù)正態(tài)分布(Drummond et al.,2006)。分布的均值也被稱為基準速率(base rate), 反映平均的演化速率大小。分布的方差則反映演化速率在樹枝之間變化的劇烈程度：方差較小時各個速率相差不大，這意味著演化速率在整棵樹上沒有明顯的差異；而方差越大，不同樹枝上速率的差異越明顯。

自相關(guān)速率模型假設(shè)后代樹枝上的演化速率依賴于臨近祖先那枝上的速率(例如r2和r5都依賴r1,r3和r4都依賴r2)。當(dāng)前樹枝上的速率一般假設(shè)服從對數(shù)正態(tài)分布(Kishino et al., 2001; Thorne and Kishino, 2002), 其均值為臨近祖先節(jié)點的速率。同理，分布的方差也反映演化速率在樹枝之間變化的劇烈程度。

這兩類速率模型往往對分異時間的估計也有影響，這主要是因為自相關(guān)速率模型會傾向于速率的變化是漸進的，而獨立速率模型沒有這種限制，會更適應(yīng)臨近樹枝間速率變化比較劇烈的情況。對化石形態(tài)數(shù)據(jù)來說，獨立速率模型可能更適用。

默認的情況下，形態(tài)數(shù)據(jù)矩陣中所有特征都共享每個枝上的演化速率，因此，這個速率代表的是所有特征的平均情況。如果需要考慮不同特征演化速率的異質(zhì)性，就需要對特征進行分區(qū)。一般可以按照不同特征類型或不同身體部位或功能來分，每個分區(qū)內(nèi)的特征共享一組演化速率，而分區(qū)之間特征演化速率的模式是獨立的，這樣就可以推斷不同部位或功能相關(guān)特征隨時間會發(fā)生怎樣的變化(Lee, 2016; Zhang and Wang,2019)。需注意的是，分區(qū)越多，每個分區(qū)內(nèi)的特征數(shù)量就越少，因此能夠估計演化速率參數(shù)的信息就越少，會造成方差很大甚至參數(shù)個數(shù)超過特征數(shù)量導(dǎo)致無法進行參數(shù)估計。因此，要在考慮演化速率異質(zhì)性和分區(qū)數(shù)量之間做一個權(quán)衡。

4 特征狀態(tài)變化

有了分異時間和演化速率，就可以計算在給定時間段t和演化速率r的情況下，形態(tài)特征從一個狀態(tài)變?yōu)榱硪粋€狀態(tài)的概率(稱為轉(zhuǎn)移概率)。這個概率由Mk模型(Lewis,2001)給出。Mk模型是描述特征狀態(tài)變化最簡單的模型，它假設(shè)狀態(tài)之間轉(zhuǎn)換的速率是相等的。這里只以兩個狀態(tài)的特征為例，用P00(r,t)表示狀態(tài)0保持不變的概率，P01(r,t)表示從0變到1的概率，P10(r,t)表示從1變到0的概率，P11(r,t)表示狀態(tài)1保持不變的概率，則P00(r,t) =P11(r,t) = 1/2 + 1/2 × exp(-2rt),P01(r,t) =P10(r,t) = 1/2 - 1/2 × exp(-2rt)。從公式中可以發(fā)現(xiàn)，時間t和速率r總是以乘積的形式出現(xiàn)，因此，在沒有化石年代信息時，兩者是不可識別的。換句話說，單單依靠形態(tài)數(shù)據(jù)來建樹，樹的枝長為時間和速率的乘積，即距離，以每個特征期望的變化次數(shù)為單位。只有同時利用化石形態(tài)和年代數(shù)據(jù)才能將分異時間和演化速率單獨估計出來。

以圖2為例，化石F1和F2的特征狀態(tài)為0, 現(xiàn)生類群S1和S2的特征狀態(tài)為1, 內(nèi)部節(jié)點的特征狀態(tài)未知用x0,x1,x2表示。F1和F2的時間分別為100和50 Ma。那么根據(jù)Mk模型，給定時間樹T= {τ, t}和速率r時，特征狀態(tài)列0011的概率為P(0011|T, r) = ∑x0∑x1∑x2Px0x1(t1-t2,r1)Px1x2(t2-t4,r2)Px21(t4,r3)Px20(t4-t5,r4)Px11(t2,r5)Px00(t1-t3,r6)。

其中西格瑪符號代表對特征在內(nèi)部節(jié)點所有可能狀態(tài)的求和。由于形態(tài)特征矩陣往往只包含可變的特征，因此這個概率還要除以所有可變狀態(tài)的概率，即P(0011 |T, r)/[1 -P(0000 |T, r) -P(1111 |T, r)]。帶有這一校正的Mk模型稱為Mkv模型(Lewis, 2001)。

假設(shè)形態(tài)矩陣中的特征都彼此獨立，那么就可以計算每一列特征在樹上的概率，再把這些概率乘起來。這個概率被稱為似然函數(shù)，表示為P(D|T, r), 其中D代表形態(tài)特征矩陣數(shù)據(jù)。

5 貝葉斯公式

在統(tǒng)計推斷時，參數(shù)都是未知的隨機變量，需要根據(jù)數(shù)據(jù)來估計它們的分布，即計算P(T, r,θ|D), 該分布稱為后驗分布，其中T= {τ, t}為時間樹，r為演化速率，θ代表其它參數(shù)(包括λ,μ,ψ等)。根據(jù)分層貝葉斯公式，可得P(T, r,θ|D) =P(D|T, r)P(r)P(T|θ)P(θ)/P(D)。等號右側(cè)分子中，第一項似然函數(shù)在第四節(jié)給出，第二項演化速率的先驗分布在第三節(jié)給出，第三項和第四項為時間樹及其參數(shù)的先驗分布在第二節(jié)給出。這樣公式分子中各項都可以計算了。分母P(D)是特征數(shù)據(jù)的概率，這需要計算對所有參數(shù)的多重積分，實際上基本無法給出解析表達式，只能通過數(shù)值算法進行近似。所以貝葉斯計算在絕大多數(shù)情況下會使用馬氏鏈蒙特卡羅算法。

6 馬氏鏈蒙特卡羅算法

馬氏鏈蒙特卡羅(MCMC)算法通過構(gòu)造馬爾可夫鏈，使其平穩(wěn)分布為要估計的后驗分布。這里為了簡化，只以一維參數(shù)的情形為示例(圖3)。實際分析中，參數(shù)一般是多維的(如τ, t, r,θ), 不過算法的思想是類似的。

MCMC采用的Metropolis-Hastings算法(Metropolis et al., 1953; Hastings, 1970)可分為如下幾步：

(1) 為θ設(shè)定任意初始值；

(2) 基于θ當(dāng)前的值，建議一個新的值θ’, 例如θ’ ～ uniform(θ- w/2,θ+ w/2);

(3) 如果π(θ’) ＞π(θ), 就接受θ’; 否則，以概率α=π(θ’) /π(θ) 接受θ’;

(4) 如果θ’被接受，就更新θ=θ’; 否則，保持θ不變；

(5) 記錄θ的值，回到步驟(2)。

注意到，在計算概率α的時候，會計算后驗分布的比，這樣后驗分布分母的部分就約掉了，只剩分子部分的比。也就是說，只要能夠把分子部分寫出解析表達式，MCMC算法就可以使用來估計參數(shù)的后驗分布了。

計算結(jié)束后，就收集到一些θ的樣本。由于參數(shù)的初始狀態(tài)往往比較差，MCMC需要經(jīng)過很多代才收斂到后驗概率密度比較高的地方，因此在估計后驗分布的時候會舍棄初始的一些樣本(burn-in), 只用MCMC鏈收斂后記錄的樣本來估計后驗分布。MrBayes默認舍棄前25%的樣本。同時，MCMC鏈還要迭代足夠多次，以保證有足夠多的有效樣本來估計參數(shù)的后驗分布。一般需要有效樣本大小(ESS)大于100。

實際運算中，最好獨立運行至少兩次MCMC, 以確保兩次的結(jié)果是一致的。有時鏈長不夠，或者不同的運算卡在不同的后驗分布區(qū)域，都會導(dǎo)致估計的結(jié)果不一致。這時調(diào)整MCMC的設(shè)置或者改善模型都可能幫助MCMC算法發(fā)揮更好的效能。使用Metropolis-coupled MCMC也是有效跨越多峰分布的手段(Lakner et al., 2008)。該算法同時運行多條MCMC鏈，一條為冷鏈(cold chain), 其余為熱鏈(hot chains), 熱鏈和冷鏈之間可以相互交換。當(dāng)然MCMC樣本只從冷鏈中采集，熱鏈只是用來幫助跨越多峰的。MrBayes默認同時獨立運行兩次，每次運算使用四條鏈，其中一條為冷鏈，其余三條為熱鏈。

圖3 馬氏鏈蒙特卡羅算法在一維情形的示例參數(shù)θ的后驗分布π(θ)為估計的目標(biāo)(曲線)Fig. 3 Illustration of the MCMC algorithm for one dimensional parameter The posterior distribution π(θ) of parameter θ is the target distribution to be estimated (curve)

7 討論

本文從貝葉斯統(tǒng)計計算的角度分層剖析了支端定年法的原理和計算過程。貝葉斯后驗分布包含先驗分布和似然函數(shù)，其中先驗分布的兩個重要組成部分：分異時間和演化速率模型，在定年分析中尤為關(guān)鍵，是影響定年準確性的主要因素。

石化生滅過程作為描述類群分異、滅絕和采樣的隨機過程，具有較大的靈活性。不過該模型還有待完善之處。現(xiàn)生類群的采樣方式可以是隨機的或多樣化的，這兩種情形可能都是極端，真實的采樣模式可能介于兩者之間，或者有的支系是隨機采樣，有的支系是多樣化采樣，但目前還沒有模型能夠支持這種情況。在更好的模型被開發(fā)出來之前，可能只能調(diào)整數(shù)據(jù)來盡量符合其中一種采樣策略。這種處理方式在現(xiàn)生類群比較少或根本沒有時一般問題不大。另外，分異速率、滅絕速率和化石采樣速率可以按分段的方式隨時間變化，不過在同一時間段內(nèi)，所有樹枝都共享同一速率。對于不同支系明顯受到的選擇壓力不同或化石保存的完整性明顯不同等情況，按支系分段而非時間分段的模型(Barido-Sottani et al., 2020)可能更合適，雖然這一類模型本身也有其他限制。不論怎樣，石化生滅過程只是作為時間的先驗分布，當(dāng)數(shù)據(jù)量比較大時(包括化石在樹上的分布程度和數(shù)量以及形態(tài)特征的數(shù)量和完整度), 數(shù)據(jù)在推斷中會起主導(dǎo)作用而先驗的影響減少。但是實際情況往往比較復(fù)雜，數(shù)據(jù)量也不盡如人意，這時考察不同先驗的影響就尤為重要。

演化速率的先驗，即形態(tài)鐘模型，也會和時間相互作用，從而影響對最終分異時間的估計。這種影響在化石較少或化石在樹上分布很不均時尤為明顯。這主要是因為化石形態(tài)數(shù)據(jù)只提供了距離的信息(每個特征期望的變化次數(shù)), 其為時間和速率的乘積(見第四節(jié))。當(dāng)缺少化石時，就沒辦法準確提供時間的信息，那么對于同樣的距離，可以是很長時間速率很慢，也可以是很短時間但速率很快，具體是怎樣就只能取決于時間和演化速率的先驗了。對于某些支系時間估計得明顯偏大或偏小，但又沒有化石來校正的情況，可以通過添加內(nèi)部節(jié)點的校正分布來得到更合理的時間估計(O’Reilly and Donoghue, 2016)。在完全沒有化石只有現(xiàn)生類群時，節(jié)點定年法(另一類型定年方法)(Yang and Rannala, 2006)就是通過使用內(nèi)部節(jié)點的校正分布來完成的。

描述形態(tài)特征狀態(tài)變化的模型也有很大的改進空間，其中涉及更多的建模和隨機模擬等工作，不是本研究的重點，這里只簡單討論一下Mk模型對定年可能的影響。Mk模型假設(shè)特征各個狀態(tài)間轉(zhuǎn)變的速率都是相等的。這種轉(zhuǎn)變有無序和有序之分(只對三個及以上狀態(tài)的特征)。無序是指特征可以直接從一個狀態(tài)變?yōu)槿我馄渌麪顟B(tài)(如從0直接變?yōu)?), 而有序是指特征只能在臨近狀態(tài)間直接變化(如從0到1, 從1到2, 再從2到3)。顯然，有序需要更多次變化(也就是更長的距離)才能從當(dāng)前狀態(tài)變?yōu)椴幌噜彽臓顟B(tài)，因此對有序的特征假設(shè)無序的變化會低估距離。更復(fù)雜的情況是，各個狀態(tài)間轉(zhuǎn)變的速率未必相等甚至相差很多，極端情況像Dollo特征甚至是不可逆的。這時使用Mk模型也會造成距離估計的偏差。在計算時還假設(shè)不同特征之間都是獨立的，如果有些特征有較強的相關(guān)性，則會導(dǎo)致演化距離的高估。前面提到，距離是時間和速率的乘積，因此在化石(時間信息)很豐富的理想情況下，距離的偏差會主要體現(xiàn)到演化速率上而對分異時間影響較小。但是分析實際數(shù)據(jù)時會更復(fù)雜一些，要具體問題具體分析。相關(guān)的工作還較少(Klopfstein et al., 2019), 需要更多后續(xù)研究來更詳細地考察。

最后提到貝葉斯計算使用的MCMC算法。該算法的策略與簡約法和似然法有明顯不同。簡約法尋找的是簡約樹長最小的樹，似然法尋求的是似然值最大時參數(shù)的估計(即最大似然樹)。因此在設(shè)計樹的搜索方法時，只要盡可能快速地找到最優(yōu)的樹就可以了。MCMC算法是為了估計參數(shù)的后驗分布，這不僅僅是一個點，而是參數(shù)的空間。因此MCMC算法的效率涉及收斂(convergence)和混合(mixing)兩部分。收斂是指MCMC達到分布概率密度高的區(qū)域，混合是指MCMC能夠按概率分布進行取樣。提高收斂速度相對容易，可以通過如簡約樹長向?qū)У姆绞絹砜焖僬业礁怕蚀蟮臉?Zhang et al.,2020)。提高混合則更困難，需要設(shè)計更好的建議(proposal)方法，這是貝葉斯計算的重點也是難點。

總之，貝葉斯支端定年法作為整合的分析方法，能夠結(jié)合化石形態(tài)和年代數(shù)據(jù)以及現(xiàn)生類群的形態(tài)和分子數(shù)據(jù)來推斷類群的系統(tǒng)關(guān)系，分異時間和演化速率，同時考慮了樹的拓撲結(jié)構(gòu)、分異時間、演化速率以及化石年代的不確定性。但該方法仍處于發(fā)展初期，模型和算法的諸多方面還亟待完善，因此還有很多工作需要做。

Supplementary material can be found on the website of Vertebrata PalAsiatica (http://www.vertpala.ac.cn/EN/2096-9899/home.shtml) in Vol. 59, Issue 4.