DEC_introduccion

Introducción al modelo DEC

Este tutorial fue traducido y modificado por Vanessa Parra a partir del tutorial “Introduction to phylogenetic biogeography with the DEC model” disponible aquí y escrito por Michael J. Landis y Sarah K. Swiston

Introducción

Muchos procesos evolutivos fundamentales, como adaptación, especiación, y extinción, ocurren en un contexto espacial. Cuando el componente histórico de dicho contexto espacial no puede observarse directamente, como suele ser el caso, se puede aplicar inferencia biogeográfica para estimar las áreas de distribución ancestrales de las especies. Este enfoque se basa en la integración de información filogenética, molecular y geográfica para modelar las distribuciones de las especies como el resultado de distintos procesos biogeográficos. Modelar adecuadamente estos procesos requiere ciertas consideraciones especiales, como la forma en que las áreas de distribución se heredan tras eventos de especiación, la influencia de eventos geológicos sobre las tasas de dispersión y los factores que afectan las tasas de dispersión y extirpación. Uno de los principales desafíos técnicos en la modelación de la evolución de las áreas de distribución consiste en cómo traducir estos procesos naturales en procesos estocásticos que sigan siendo manejables para la inferencia.

Este tutorial presenta primero una introducción general a algunos de estos modelos y, posteriormente, describe como realizar inferencia bayesiana de biogeografía histórica utilizando un modelo de Dispersión–Extinción–Cladogénesis (DEC) en RevBayes.

Descripción general del modelo de Dispersión–Extinción–Cladogénesis (DEC)

El proceso de Dispersión–Extirpación–Cladogénesis (DEC) modela la evolución de los rangos geográficos como un proceso de valores discretos Ree et al. 2005; Ree y Smith 2008(#ree2008). Existen tres componentes clave para comprender el modelo DEC: el rango como caracter, la evolución anagenética del rango y la evolución cladogenética del rango (Figura 1).

Figura 1. Esquema del comportamiento del modelo DEC. Se muestran dos eventos anagenéticos (a,b) y cinco eventos cladogenéticos (c–g) para un sistema con dos áreas. Las áreas están sombreadas cuando están habitadas por un linaje determinado y se dejan en blanco cuando están deshabitadas. El tiempo avanza de izquierda a derecha. (a) Dispersión: se añade una nueva área al rango de la especie. (b) Extirpación (o extinción local): el rango de la especie pierde un área previamente habitada. (c) Simpatría estrecha: cuando el rango ancestral contiene un área, ambos linajes hijos heredan esa área. (d) Simpatría de subconjunto: cuando el rango ancestral es amplio, un linaje hijo hereda el rango ancestral y el otro hereda solo un área. (e) Alopatría (o vicarianza): cuando el rango ancestral es amplio, un linaje hijo hereda un subconjunto de las áreas ancestrales mientras que el otro linaje hijo hereda todas las áreas ancestrales restantes. (f) Simpatría generalizada: cuando el rango ancestral está muy extendido, ambos linajes hijos heredan el rango ancestral. (g) Dispersión por salto (o especiación fundadora): un linaje hijo hereda el rango ancestral mientras que el otro hereda una nueva área desocupada.

Caracteres de rango discreto

DEC trata los rangos de taxones como datos de presencia-ausencia, es decir, el rango es el conjunto de áreas discretas (pre-definidas) donde una especie se observa. Por ejemplo, digamos que hay tres áreas, A, B y C. Si una especie está presente en las áreas A y C, entonces su rango es igual a AC, que también se puede codificar en el vector de bits de longitud 3, 101. Los vectores de bits también se pueden transformar en un estado de valor entero, por ejemplo , el número binario 101 es igual al entero 5. Ten en cuenta que hay que agregar 1 al valor entero del estado de interés para acceder a ese estado desde un objeto RevBayes, por ejemplo , el rango AC con valor entero 5 se accede en el índice 6.

Tabla 1. Ejemplo de representaciones discretas de rangos geográficos en un análisis con las áreas A, B y C.

Range	Bits	Size	Integer
∅	000	0	0
A	100	1	1
B	010	1	2
C	001	1	3
AB	110	2	4
AC	101	2	5
BC	011	2	6
ABC	111	3	7

La representación decimal de los estados de rango rara vez se utiliza en las discusiones, pero es útil tenerla en cuenta al considerar el número total de posibles distribuciones para una especie y al procesar los resultados.

Evolución anagenética del rango

En el contexto del modelo DEC, anagénesis se refiere a la evolución del rango geográfico que ocurre entre eventos de especiación, a lo largo de las ramas del linaje. Existen dos tipos de eventos anagenéticos: dispersión (Figura 1a) y extinción local o extirpación (Figura 1b).

Dado que el modelo DEC utiliza rangos geográficos de valores discretos, la anagénesis se modela mediante una cadena de Markov en tiempo continuo. Este enfoque permite calcular la probabilidad de transición de un carácter que cambia del estado i al estado j en un tiempo t mediante la exponenciación matricial

\[ P_{ij}(t) = \left[ \exp(Qt) \right]_{ij} \] donde \(Q\) es la matriz de tasas instantáneas que define las tasas de cambio entre todos los pares de caracteres, y \(P\) es la matriz de probabilidades de transición. Los índices i y j representan diferentes rangos geográficos, cada uno de los cuales se codifica como el conjunto de áreas ocupadas por la especie en un momento t. La probabilidad se integra sobre todos los escenarios posibles de transiciones de caracteres que pueden ocurrir durante el intervalo de tiempo t, siempre que la cadena comience en el rango i y termine en el rango j.

De esta manera, la matriz \(Q\) puede parametrizarse para reflejar las clases permitidas de eventos de evolución de rango con parámetros biológicamente significativos. Para tres áreas,las tasas incluidas en la matriz de tasas anagenéticas son:

Instantaneous rate matrix (\(Q\)) for three areas. States are ordered as \(\varnothing\), A, B, C, AB, AC, BC, ABC.

	∅	A	B	C	AB	AC	BC	ABC
∅	−	0	0	0	0	0	0	0
A	\(e_A\)	−	0	0	\(d_{AB}\)	\(d_{AC}\)	0	0
B	\(e_B\)	0	−	0	\(d_{BA}\)	0	\(d_{BC}\)	0
C	\(e_C\)	0	0	−	0	\(d_{CA}\)	\(d_{CB}\)	0
AB	0	\(e_B\)	\(e_A\)	0	−	0	0	\((d_{AC}+d_{BC})\)
AC	0	\(e_C\)	0	\(e_A\)	0	−	0	\((d_{AB}+d_{CB})\)
BC	0	0	\(e_C\)	\(e_B\)	0	0	−	\((d_{BA}+d_{CA})\)
ABC	0	0	0	0	\(e_C\)	\(e_B\)	\(e_A\)	−

donde \(e = (e_A, e_B, e_C)\) son las tasas de extinción local por área, y \(d = (d_{AB}, d_{AC}, d_{BC}, d_{BA}, d_{CA}, d_{CB})\) son las tasas de dispersión entre áreas. Nótese que la suma de las tasas de salida del rango nulo (∅) es cero; es decir, cualquier linaje que pierde todas las áreas de su rango permanece en ese estado de manera permanente.

Para desarollar nuestra intuición, a continuación construiremos una matriz de tasa DEC en RevBayes.

Crea un nuevo directorio en tu computadora llamado RB_biogeo_tutorial.
Navega hasta el directorio del tutorial y ejecuta el archivo binario de RevBayes.

Supongamos que tenemos tres áreas

n_areas <- 3

Primero, creamos una matriz de tasas de dispersión entre pares de áreas, con tasas \(d_{AB} = d_{AC} = \cdots = d_{CB} = 1\)

for (i in 1:n_areas) {
    for (j in 1:n_areas) {
        dr[i][j] <- 1.0
    }
}

A continuación, vamos a crear las tasas de extirpación con valores \(e_{A} = e_{B} = e_{C} = 1\)

for (i in 1:n_areas) {
    for (j in 1:n_areas) {
        er[i][j] <- 0.0
    }
    er[i][i] <- 1.0
}

Cuando la matriz de tasas de extirpación es una matriz diagonal (es decir, todas las entradas fuera de la diagonal son cero), las tasas de extirpación son mutuamente independientes, tal como se describe en Ree et al. (2005). En otros escenarios se podrían explorar modelos más complejos que penalizan áreas de distribución extensas que abarcan áreas desconectadas.

Para continuar, creamos la matriz de tasas del modelo DEC a partir de las tasas de dispersión (dr) y las tasas de extirpación (er).

Q_DEC := fnDECRateMatrix(dispersalRates=dr, extirpationRates=er)
Q_DEC


[ [ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000 ] ,
     [ 1.0000, -3.0000, 0.0000, 0.0000, 1.0000, 1.0000, 0.0000, 0.0000 ] ,
     [ 1.0000, 0.0000, -3.0000, 0.0000, 1.0000, 0.0000, 1.0000, 0.0000 ] ,
     [ 1.0000, 0.0000, 0.0000, -3.0000, 0.0000, 1.0000, 1.0000, 0.0000 ] ,
     [ 0.0000, 1.0000, 1.0000, 0.0000, -4.0000, 0.0000, 0.0000, 2.0000 ] ,
     [ 0.0000, 1.0000, 0.0000, 1.0000, 0.0000, -4.0000, 0.0000, 2.0000 ] ,
     [ 0.0000, 0.0000, 1.0000, 1.0000, 0.0000, 0.0000, -4.0000, 2.0000 ] ,
     [ 0.0000, 0.0000, 0.0000, 0.0000, 1.0000, 1.0000, 1.0000, -3.0000 ] ]

Calculamos las probabilidades de transición anagenética para una rama de longitud 0.2.

tp_DEC <- Q_DEC.getTransitionProbabilities(rate=0.2)
tp_DEC

[ [ 1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000 ], [ 0.000, 0.673, 0.013, 0.013, 0.123, 0.123, 0.005, 0.050 ], [ 0.000, 0.013, 0.673, 0.013,
   0.123, 0.005, 0.123, 0.050 ], [ 0.000, 0.013, 0.013, 0.673, 0.005, 0.123, 0.123, 0.050 ], [ 0.000, 0.107, 0.107, 0.004, 0.502, 0.031, 0.031, 0.218 ], [
   0.000, 0.107, 0.004, 0.107, 0.031, 0.502, 0.031, 0.218 ], [ 0.000, 0.004, 0.107, 0.107, 0.031, 0.031, 0.502, 0.218 ], [ 0.000, 0.021, 0.021, 0.021, 0.107,
   0.107, 0.107, 0.616 ] ]

Observa cómo la estructura de la matriz de tasas se refleja directamente en la matriz de probabilidades de transición. Por ejemplo, los rangos que están separados por múltiples eventos de dispersión y extirpación presentan las probabilidades más bajas: la transición desde el rango A hasta BC requiere un mínimo de tres eventos y tiene una probabilidad de 0.005.

Asimismo, observa que la probabilidad de entrar o salir del rango nulo es cero. De manera predeterminada, RevBayes condiciona el proceso de evolución anagenética del rango a no entrar nunca en el rango nulo al calcular las probabilidades de transición (nullRange = "CondSurv"). Esto permite que el modelo utilice las mismas probabilidades de transición tanto para la simulación como para la inferencia. Massana et al. (2015) señalaron por primera vez que el rango nulo, un estado absorbente no observado, da lugar a estimaciones anómalas de las tasas de extirpación y del tamaño del rango. Como solución, propusieron eliminar el rango nulo del espacio de estados, lo cual puede implementarse mediante la opción nullRange = "Exclude". La opción nullRange = "Include" no aplica ningún tratamiento especial al rango nulo y produce las probabilidades originales descritas por Ree et al. (2005).

Evolución cladogenética del rango

El componente cladogenético del modelo DEC describe los cambios evolutivos que acompañan a los eventos de especiación (Figuras 1c–g). En el contexto de la evolución del rango geográfico, las especies hijas no necesariamente heredan el rango ancestral de manera idéntica. Para cada nodo interno del árbol reconstruido, puede ocurrir uno de varios tipos de eventos cladogenéticos, algunos de los cuales se describen a continuación.

Comenzando por el caso más simple, supongamos que el rango de una especie es \(A\) en el momento previo a un evento de especiación en un nodo interno. Dado que el rango ancestral de la especie es de tamaño 1, ambos línajes hijos necesariamente heredan el rango ancestral (\(A\)). En el lenguaje de DEC, este proceso se denomina simpatría estrecha (narrow sympatry; Figura 1c).

Ahora, supongamos que el rango ancestral es \(ABC\). Bajo un escenario de simpatría de subconjuntos (subset sympatry), uno de los linajes hijos hereda de manera idéntica el rango ancestral (\(ABC\)), mientras que el otro hereda únicamente una sola área, es decir, \(A\), \(B\) o \(C\) (Figura 1d).

En el caso de cladogénesis alopátrica, el rango ancestral se divide entre las dos líneas hijas, por ejemplo, una puede heredar \(AB\) y la otra \(C\) (Figura 1e). En la cladogénesis simpátrica de rangos amplios (widespread sympatric cladogenesis), ambas linajes hijos heredan el rango ancestral completo (\(ABC\)) (Figura 1f).

Finalmente, si el rango ancestral es \(A\), la cladogénesis por dispersión de salto (jump dispersal cladogenesis) da lugar a que un linaje hijo herede el rango ancestral (\(A\)), mientras que la otra hereda un área previamente no ocupada, como \(B\) o \(C\) (Figura 1g). Para una revisión detallada de los distintos tipos de transiciones cladogenéticas descritas en la literatura, véase Matzke (2012).

A continuación, construímos la matriz de probabilidad de eventos cladogenéticos:

clado_event_types = [ "s", "a" ]
clado_event_probs <- simplex( 1, 1 )
P_DEC := fnDECCladoProbs(eventProbs=clado_event_probs,
                         eventTypes=clado_event_types,
                         numCharacters=n_areas)

clado_event_types define qué tipos de eventos cladogenéticos se utilizan en el modelo. Los códigos "a" y "s" representan alopatría y simpatría de subconjuntos, respectivamente, tal como se describen en Ree et al. (2005). Otros tipos de eventos cladogenéticos incluyen la dispersión por salto "j" Matzke (2014) y la simpatría completa "f" Landis et al. 2013.

La matriz de probabilidad de eventos cladogenéticos asume que los vectores eventProbs y eventTypes comparten el mismo orden, de modo que cada probabilidad se asocia correctamente con su tipo de evento correspondiente.

A continuación, imprimimos las probabilidades de transición cladogenética.

P_DEC

       [
         ( 1 -> 1, 1 ) = 1.0000,
         ( 2 -> 2, 2 ) = 1.0000,
         ( 3 -> 3, 3 ) = 1.0000,
         ...
         ( 7 -> 7, 1 ) = 0.0833,
         ( 7 -> 7, 2 ) = 0.0833,
         ( 7 -> 7, 3 ) = 0.0833
       ]

La matriz de probabilidad cladogenética se vuelve muy dispersa cuando tenemos un gran número de áreas, por lo que únicamente se muestran los valores distintos de cero. Cada fila reporta un triplete de estados (el estado ancestral y los dos estados hijos) junto con la probabilidad asociada a ese evento. Como se trata de probabilidades dado un evento de cladogénesis, la suma de las probabilidades de todos los resultados cladogenéticos posibles, para cierto estado ancestral, es igual a uno.

Para tener en cuenta

Las probabilidades de cambio anagenético a lo largo de los linajes deben considerar todas las combinaciones posibles de estados iniciales y finales. Para 3 áreas, existen 8 estados posibles, y por lo tanto, \(8×8=64\) probabilidades correspondientes a pares de estados. En el caso de los cambios cladogenéticos, se requieren probabilidades de transición para todas las combinaciones posibles del estado antes de la cladogénesis, después de la cladogénesis para el linaje izquierdo y después de la cladogénesis para el linaje derecho. Como en el caso anterior, para tres áreas, hay \(8×8×8=512\) probabilidades cladogenéticas.

por supuesto, este modelo puede especificarse para más de tres áreas. Sin embargo, resulta ilustrativo considerar cómo crece el tamaño de la matriz \(Q\) a medida que aumenta el número de áreas \(N\). Para tres áreas, la matriz \(Q\) tiene dimensiones \(8×8\). Para diez áreas, \(Q\) es de tamaño \(2^{10} \times 2^{10} = 1024 \times 1024\), lo que se aproxima al límite práctico de tamaño para matrices que pueden ser exponenciadas en un tiempo razonable. Para veinte áreas, \(Q\) tendría dimensiones de \(2^{20} \times 2^{20} \approx 10^{6} \times 10^{6}\), haciendo inviable la exponenciación matricial. Por lo tanto, la selección de las áreas discretas en un análisis DEC debe realizarse considerando cuidadosamente los objetivos del análisis y el nivel de complejidad que se desea capturar.