11.4: Пристосування моделей народження-смерті до часів розгалуження

Інший підхід, який використовує більше інформації в філогенетичному дереві, включає пристосування моделей народження-смерті до розподілу часу розгалуження. Цей підхід простежується аж до Юле (1925), який вперше застосував моделі стохастичних процесів до росту філогенетичних дерев. Зовсім недавно серія робіт Раупа та його колег (Raup et al. 1973; Raup and Gould 1974; Gould et al. 1977; Raup 1985) підштовхнула сучасні підходи до кількісної макроеволюції шляхом імітації випадкових кладів, а потім демонструючи, як змінні такі клади, вирощені під простими моделями народження-смерті, можуть бути.

Більшість сучасних підходів до пристосування моделей народження-смерті до філогенетичних дерев використовують інтервали між подіями видоутворення на дереві - «час очікування» між послідовним видобутком - для оцінки параметрів моделей народження-смерті. На малюнку 10.2 показані ці терміни очікування. Часто відомості про закономірність накопичення видів у філогенетичному дереві узагальнюються ділянкою лінійно-наскрізного часу (LTT), який є сюжетом кількості родовищ у дереві проти часу (див. Рис. Як я ввів у главі 10, вісь y графіків LTT трансформується в журнал, так що очікуваний візерунок за моделлю чистого народження з постійною швидкістю є прямою лінією. Зауважте також, що графіки LTT ігнорують відносний порядок подій видоутворення. Stadler (2013b) називає моделі, що обґрунтовують такий підхід, «видообмінними» моделями - ми можемо змінити ідентичність видів в будь-який момент часу без зміни очікуваної поведінки моделі. Через це підходи до розуміння моделей народження-смерті на основі часу розгалуження відрізняються від - і доповнюють - підходів, заснованих на топології дерева, як баланс дерева.

Як обговорювалося в попередньому розділі, хоча ми часто не маємо інформації про вимерлі види в кладі, ми все ще можемо (теоретично) зробити висновок про наявність вимирання з ділянки LTT. Сигналом вимирання є надлишок молодих родовищ, що розглядається як «тяга недавнього» на наших ділянках LTT (рис. 10.10). У наступному розділі я продемонструю, як статистичні підходи можуть більш суворо зафіксувати цю закономірність.

Розділ 11.4a: Імовірність часу очікування за моделлю народження-смерть

Для того, щоб використовувати методи ML та Байєса для оцінки параметрів моделей народження-смерть за порівняльними даними, нам потрібно записати ймовірності часу очікування між подіями видоутворення в дереві. У літературі є трохи варіацій у позначеннях, тому я буду слідувати за Стадлером (2013b) та Maddison et al. (2007), серед інших, для підтримки послідовності. Будемо вважати, що клад починається в час t ₁ з пари видів. Більшість аналізів дотримуються цієї конвенції, і умова, що процес починається в момент t ₁, представляючи вузол в корені дерева. Це має сенс, оскільки ми рідко маємо інформацію про вік стовбура нашої клади. Ми також обумовимо обидві ці початкові лінії, що збереглися до наших днів, оскільки це вимога для отримання дерева з цим віком крони (наприклад, рівняння Stadler 2013a 5).

Події видоутворення та вимирання відбуваються в різний час, і процес закінчується в момент 0, коли клада має n існуючих видів - тобто ми вимірюємо час назад від сьогоднішніх днів. Вимирання призведе до видів, які не поширюються весь шлях до часу 0. Наразі ми будемо вважати, що у нас є лише дані про існуючі види. Ми будемо посилатися на філогенетичне дерево, яке показує час розгалуження, що ведуть до існуючих видів, як реконструйоване дерево (Nee et al. 1994). Для реконструйованого дерева з n видами існує n − 1 раз видоутворення, яке ми позначимо як t ₁, t ₂, t ₃,..., t _{n − 1}. Листя нашого ультраметричного дерева все закінчується в момент 0.

Зауважте, що в цьому позначенні t ₁ > t ₂ >... > t _{n − 1} > 0, тобто наші часи видоутворення вимірюються назад від кінчиків, і коли ми збільшуємо індекс, час постійно зменшується [це важливо нотаційна різниця між Stadler (2013a), що використовується тут, і Nee (1994 та інші), остання з яких враховує часові інтервали між подіями видоутворення, наприклад t ₁ − t ₂ у наших позначеннях]. Наразі ми припустимо повну вибірку; тобто всі n видів, що живуть на сьогоднішній день, представлені на дереві.

Тепер ми виведемо ймовірність спостереження множини разів видоутворення t ₁, t ₂,..., t _{n − 1} з огляду на існуючу різноманітність clade, n та параметрів нашої моделі народження-смерть λ і μ. Для цього ми дуже уважно стежимо за підходом, заснованим на диференціальних рівняннях, запроваджених Maddison et al. (2007).

Загальна ідея полягає в тому, що ми присвоїмо значення цим ймовірностям на кінчиках дерева, а потім визначимо набір правил для їх оновлення, коли ми перетікаємо назад через дерево від кінчиків до кореня. Коли ми приїдемо до кореня дерева, ми матимемо ймовірність спостерігати власне дерево з урахуванням нашої моделі - тобто ймовірність. Це ще один алгоритм обрізки.

Для початку нам потрібно відстежувати дві ймовірності: D _N (t), ймовірність того, що родовід в якийсь час t в минулому еволюціонує в існуючий клад N, як спостерігається сьогодні; і E (t), ймовірність того, що родовід у деяких час t повністю вимерне і не залишить нащадків в наші дні. (Пізніше ми переглянемо E (t), щоб він включав можливість того, що лінії мають нащадків, але жодна не була вибіркова в наших даних). Потім ми застосовуємо ці ймовірності до дерева, використовуючи три основні ідеї (рис. 11.4)

Малюнок 11.4. Загальна схема етапів розрахунку ймовірності дерева за моделлю народження-смерть. Образ, натхненний Меддісоном та ін. (2007) і створений автором, може бути використаний повторно за ліцензією CC-BY-4.0.

1. Визначаємо свої вихідні точки на кінчиках дерева.

2. Ми визначаємо, як змінюються ймовірності, визначені в (1), коли ми рухаємося назад уздовж гілок дерева.

3. Ми визначаємо, що відбувається з нашими ймовірностями у вузлах дерева.

Потім, починаючи з кінчиків дерева, пробираємося до кореня. На кожному кінчику ми маємо початкове значення як для D _N (t), так і для E (t). Ми рухаємося назад уздовж гілок дерева, оновлюючи обидві ймовірності, як ми йдемо, використовуючи крок 2. Коли дві гілки збираються разом у вузлі, ми об'єднуємо ці ймовірності за допомогою кроку 3.

Таким способом пройдемося по дереву, починаючи з кінчиків і проходячи над кожною гілкою і вузлом (рис. 11.4). Коли ми дійдемо до кореня, ми матимемо D _N (t _{r o o t}), що є повною ймовірністю, що ми хочемо.

Ви можете задатися питанням, чому нам потрібно обчислити як D _N (t), так і E (t), якщо ймовірність захоплена D _N (t) в корені. Причина полягає в тому, що ймовірність спостереження за деревом залежить від цих ймовірностей вимирання, обчислених назад у часі. Нам потрібно відстежувати E (t), щоб знати про D _N (t) і як він змінюється. Нижче ви побачите, що E (t) з'являється безпосередньо в наших диференціальних рівняннях для D _N (t).

По-перше, відправна точка. Оскільки кожен наконечник i представляє живу лінію, ми знаємо, що вона жива на сьогоднішній день - тому ми можемо визначити D _N (t) =1. Ми також знаємо, що він не згасне до включення в дерево, тому E (t) = 0. Це дає наші початкові значення для двох ймовірностей на кожному кінчику дерева (рис. 11.5).

Малюнок 11.5. Початкові точки на кінчиках дерев для обчислення ймовірності ймовірності. Образ, натхненний Меддісоном та ін. (2007) і створений автором, може бути використаний повторно за ліцензією CC-BY-4.0.

Далі, уявіть, що ми рухаємося назад вздовж деякої ділянки гілки дерева без вузлів. Ми розглянемо довільну гілку дерева. Оскільки ми йдемо назад у часі, ми почнемо з якогось вузла в дереві N, який відбувається за раз t _N, і позначимо час, що йде назад в минуле як t (рис. 11.6).

Малюнок 11.6. Оновлення D _N (t) та обчислення E (t) уздовж гілки дерева. Образ, натхненний Меддісоном та ін. (2007) і створений автором, може бути використаний повторно за ліцензією CC-BY-4.0.

Оскільки ця ділянка гілки існує в нашому дереві, ми знаємо дві речі: рід не вимерв за цей час, і якщо видоутворення відбулося, рід, який відколовся, не дожив до наших днів. Ми можемо зафіксувати ці дві можливості в диференціальному рівнянні, яке враховує, як змінюється наша загальна ймовірність протягом деякої дуже малої одиниці часу (Maddison et al. 2007).

\[ \frac{dD_N(t)}{dt} = -(\lambda + \mu) D_N(t) + 2 \lambda E(t) D_N(t) \label{11.12} \]

Тут перша частина рівняння, − (λ + μ) D _N (t), представляє ймовірність не видоутворення і не згасне, тоді як друга частина, 2 λ E (t) D _N (t), являє собою ймовірність видоутворення з подальшим остаточним вимиранням однієї з двох дочірніх родовищ. 2 у цьому рівнянні з'являється тому, що ми повинні враховувати той факт, що, слідуючи видобутку від предка до дочок А та Б, ми побачимо однакову закономірність незалежно від того, хто з двох нащадків дожив до теперішнього часу.

Нам також потрібно обчислити нашу ймовірність вимирання, що йде назад через дерево (Maddison et al. 2007):

\[ \frac{dE(t)}{dt} = \mu - (\mu + \lambda) E(t) + \lambda E(t)^2 \label{11.13} \]

Три частини цього рівняння представляють три способи, якими родовід може не досягти цього дня: або він вимерне протягом розглянутого інтервалу (μ), він переживає цей інтервал, але згасає деякий час пізніше (− (μ + λ) E (t )), або він видобуває в інтервалі, але обидва нащадки вимерли до сьогоднішнього дня (λ E (t) ²) (Maddison et al. 2007). На відміну від терміна D _N (t), ця ймовірність залежить тільки від часу, а не від топологічної структури дерева.

Ми також уточнимо, що λ > μ; можна розслабити це припущення, але це ускладнює рішення.

Ми можемо вирішити ці рівняння так, що ми зможемо оновити ймовірність переміщення назад уздовж будь-якої гілки дерева довжиною t. спочатку розв'язавши рівняння 11.13 і використовуючи нашу початкову умову E (0) = 0:

\[ E(t) = 1 - \frac{\lambda-\mu}{\lambda - (\lambda-\mu)e^{(\lambda - \mu)t}} \label{11.14} \]

Тепер ми можемо замінити цей вираз для E (t) в Equation\ ref {11.12} і вирішити

\[ D_N(t) = e^{-(\lambda - \mu)(t - t_N)} \frac{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_N})^2}{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t})^2} \cdot D_N(t_N) \label{11.15}\]

Пам'ятайте, що t _N - глибина (виміряна з сьогоднішнього дня) вузла N (рис. 11.6).

Нарешті, нам потрібно розглянути, що відбувається, коли дві гілки збираються разом у вузлі. Оскільки існує вузол, ми знаємо, що було подія видоутворення. Розрахунки ймовірностей, що стікають по кожній гілці, множимо на ймовірність події видоутворення [Maddison et al. (2007); рис. 11.7].

Малюнок 11.7. Оновлення D _N (t) та E (t) уздовж гілки дерева. Образ, натхненний Меддісоном та ін. (2007) і створений автором, може бути використаний повторно за ліцензією CC-BY-4.0.

Отже:

\[D_{N′}(t)=D_N(t)D_M(t)λ \label{11.16}\]

Де clade N '- це клад, що складається з комбінації двох сестринських кладів N і M.

Щоб застосувати цей підхід до всього філогенетичного дерева, ми множимо Рівняння\ ref {11.15} і\ ref {11.16} по всіх гілках і вузлах дерева. Таким чином, повна ймовірність є (Maddison et al. 2007; Morlon et al. 2011):

\[ L(t_1, t_2, \dots, t_n) = \lambda^{n-1} \big[ \prod_{k = 1}^{2n-2} e^{(\lambda-\mu)(t_{k,b} - t_{k,t})} \cdot \frac{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_{k,t}})^2}{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_{k,b}})^2} \big] \label{11.17}\]

Тут n - кількість підказок у дереві (зверніть увагу, що оригінальна похідна в Maddison використовує n як кількість вузлів, але я змінив його для узгодженості з рештою книги).

Добуток у Equation\ ref {11.17} взято на всі 2 n − 2 гілки дерева. Кожна гілка k має два рази пов'язані з ним, один до основи дерева, t _{k, b}, і один до кінчиків, t _{k, t}.

Більшість методів підгонки моделей народження-смерті до дерев обумовлюють існування дерева - тобто обумовлюючи те, що весь процес не вимерв до наших днів, а подія видоутворення з кореневого вузла призвела до двох збережених родовищ. Для цього обумовлення ділимо рівняння\ ref {11.17} на λ [1 − E (t _{r o o t})] ² (Morlon et al. 2011; Stadler 2013a).

Крім того, ймовірності для часу очікування народження-смерті, наприклад, у вихідному похідному від Nee, включають додатковий термін (n − 1)!. Це тому, що є (n − 1)! можливі топології для будь-якої множини n − 1 часу очікування, всі однаково вірогідні. Оскільки цей термін є постійним для заданого розміру дерева n, то його залишення не має ніякого впливу на відносні ймовірності різних значень параметрів - але про цей множник необхідно знати, якщо порівнювати ймовірності різних моделей для вибору моделі, або порівнюючи вихід різні програми (Стадлер 2013а).

Враховуючи ці два фактори, повна ймовірність становить:

\[ L(\tau) = (n-1)! \frac{\lambda^{n-2} \big[ \prod_{k = 1}^{2n-2} e^{(\lambda-\mu)(t_{k,b} - t_{k,t})} \cdot \frac{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_{k,t}})^2}{(\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_{k,b}})^2} \big]}{ [1-E(t_{root})]^2} \label{11.18}\]

де:

\[ E(t_{root}) = 1 - \frac{\lambda-\mu}{\lambda - (\lambda-\mu)e^{(\lambda - \mu)t_{root}}}\label{11.19}\]

Розділ 11.4b: Використання максимальної ймовірності для відповідності моделі народження-смерті

З огляду на рівняння 11.19 для ймовірності, ми можемо оцінити рівень народжуваності та смертності, використовуючи як ML, так і байєсівський підходи. Для оцінки ML максимізуємо рівняння 11.19 над λ і μ. Для моделі чистого народження ми можемо встановити μ = 0, а максимальну оцінку ймовірності λ можна обчислити аналітично як:

\[ \lambda= \frac{n-2}{s_{branch}} \label{11.20} \]

де s _{b r a n c h} - сума довжин гілок в дереві,

\[ s_{branch} = \sum_{i=1}^{n-1} t_i + t_{n-1} \label{11.21} \]

Рівняння\ ref {11.21} також називається оцінкою Кендала-Морана швидкості видоутворення (Nee 2006).

Для моделі народження-смерть ми можемо використовувати числові методи для максимізації ймовірності над λ та μ.

Наприклад, ми можемо використовувати ML, щоб пристосувати модель народження-смерті до дерева Lupinus (Drummond et al. 2012), яке налічує 137 видів верхівки та загальний вік 16,6 мільйона років. Роблячи так, отримаємо оцінки параметрів ML λ = 0,46 і μ = 0,20, при лог-ймовірності l n L _{b d} = 262,3. Порівняйте це з моделлю чистого народження на тому ж дереві, яка дає λ = 0,35 і l n L _{p b} = 260,4. Можна порівняти придатність цих двох моделей за допомогою балів AIC: A I C _{b d} = −520.6 та A I C _{p b} = −518,8, тому модель народження-смерть має кращий (нижчий) показник АПК, але менше ніж на дві одиниці АПК. Випробування коефіцієнта правдоподібності, яке дає Δ = 3,7 і P = 0,054. Іншими словами, ми оцінюємо ненульову швидкість вимирання в кладі, але докази, що підтверджують цю модель над моделлю чистого народження, не особливо сильні. Навіть якщо цей вибір моделі трохи неоднозначний, пам'ятайте, що ми також оцінили параметри, використовуючи всю інформацію, яку ми маємо в час очікування філогенетичного дерева.

Розділ 11.4c: Використання байєсівських MCMC для відповідності моделі народження-смерті

Ми також можемо оцінити рівень народжуваності та смертності за допомогою байєсійського MCMC. Ми можемо використовувати точно описаний вище метод для віку та різноманітності кладів, але замінити Equation\ ref {11.11} для ймовірності, таким чином, використовуючи час очікування, отриманий з філогенетичного дерева, для оцінки параметрів моделі.

Застосовуючи це до люпин з тими ж попередніми ознаками, що і раніше, отримаємо задні розподіли, показані на малюнку 11.5. Середнє значення заднього для кожного параметра дорівнює λ = 0,48 і μ = 0,23, досить близьке до оцінок МЛ за цими параметрами.

Малюнок 11.8. Задній розподіл для b і d для люпінуса. Дані з (Drummond et al. 2012), образ, натхненний Меддісоном та ін. (2007) і створений автором, може бути використаний повторно за ліцензією CC-BY-4.0.