Nota al lector: Si utilizas estos materiales e imágenes para su clase o encuentras algún error, ponte en contacto con Rosana Zenil-Ferguson en . Agradezco que reconozcas este trabajo.

Modelo de evolución para un caracter discreto

En este taller nos centraremos en definir un modelo de macroevolución para un caracter discreto, así como en evaluar su inferencia. Para ello, utilizaremos un tipo de modelo estocástico muy conocido llamado cadena de Markov. Las cadenas de Markov son herramientas muy populares entre los desarrolladores de modelos porque siguen un proceso aleatorio a lo largo del tiempo. En la literatura macroevolutiva, algunos de los modelos más populares incluyen el movimiento browniano (BM) y el Ornstein-Uhlenbeck (OU), que son cadenas de Markov que siguen un rasgo continuo a lo largo del tiempo. Estos dos modelos tan populares poseen propiedades matemáticas que nos permiten estimar las tasas evolutivas (cambio a lo largo del tiempo) utilizando las distancias filogenéticas.

Para estados discretos, la situación es diferente. También utilizamos cadenas de Markov, denominadas modelos Mkn, pero no podemos incluir fácilmente las distancias filogenéticas de la misma manera que los modelos BM u OU. Como veremos en las siguientes secciones, para cada rama de la filogenia debemos calcular una matriz de probabilidad, lo que dificulta su cálculo.

Para todos los ejemplos que analizaremos en este taller asumimos un árbol filogenético cuyas puntas son estados discretos de un rasgo (Fig. 1). Antes de pasar a las secciones computacionales, comenzaré con la notación formal y la enumeración de los nodos. Los nodos se enumeran siguiendo el formato Newick de R: primero las puntas, luego la raíz y, finalmente, los nodos internos, ordenados en el tiempo.


Figura 1. Árbol filogenético simple con un rasgo discreto. Los números en azul representan los números de nodo según el formato Newick en R
Figura 1. Árbol filogenético simple con un rasgo discreto. Los números en azul representan los números de nodo según el formato Newick en R

¿Qué es una Cadena de Markov (CTMC)?

Un proceso de Markov \(X_t\) es un modelo estocástico cuya propiedad es que la probabilidad de cualquier comportamiento futuro particular del proceso, cuando se conoce con exactitud su estado actual, no se ve alterada por el conocimiento adicional sobre su comportamiento pasado. En otras palabras, “el futuro sólo depende del presente y no del pasado”. El tiempo \(t\) puede ser continuo o discreto, y los valores de la variable que seguimos a través del tiempo \(X(t)\) pueden ser continuos (es decir, la masa corporal) o discretos (es decir, el número de cromosomas).

Los siguientes son tres supuestos importantes que consideraremos para este taller:

  1. En métodos filogenéticos comparativos, utilizamos \(\{X(t), t\geq 0\}\) para representar la evolución de un caracter a lo largo del tiempo \(t\), medido en millones de años.

  2. Los rasgos que modelaremos aquí son discretos, lo que significa que \(X(t)\) solo toma valores de enteros no negativos \(\mathbb{Z}^+=\{0,1,2,...\}\). Usaremos la notación \(X_i^n(t)\) para representar lo que sucede a lo largo del árbol filogenético. \(X\) es el caracter, \(t\) es el tiempo, \(i\) es el valor del rasgo y \(n\) es el número de nodo.


Figura 2. Notación para una cadena de Markov de tiempo continuo sobre un árbol filogenético al rastrear un caracter a lo largo de los nodos
Figura 2. Notación para una cadena de Markov de tiempo continuo sobre un árbol filogenético al rastrear un caracter a lo largo de los nodos


  1. La propiedad markoviana “el futuro no depende del pasado, sólo del presente” se representa matemáticamente con una probabilidad condicional \[{P}(X_i(t)|X_j(u), X_k(v))={P}(X_i(t)|X_j(u)) \textrm{ with } v < u < t\]

Figura 3. La propiedad markoviana establece que el camino rojo y el azul tienen exactamente la misma probabilidad.
Figura 3. La propiedad markoviana establece que el camino rojo y el azul tienen exactamente la misma probabilidad.


Nota adicional sobre la propiedad markoviana

La propiedad markoviana, para los árboles filogenéticos, establece que, dado un nodo, basta con centrarse en su ancestro común más reciente para calcular la probabilidad de un cambio de estado.

Probabilidades condicionales de las cadenas de Markov

Como se indica en la propiedad markoviana, estamos interesados en la probabilidad condicional \(P(X_i^n(t)|X_j^m(u))\) de pasar de un estado \(i\) a \(j\) en \(t-u\) unidades de tiempo. Esta probabilidad suele representarse mediante la notación \(p_{ij}(t-u)\).

Por ejemplo, en la Fig. 3, seguimos la evolución de un rasgo \(X\) que tiene estados \(0,1\). Imaginemos, por ejemplo, que este caracter corresponde a la presencia o ausencia de cola en reptiles. Nos interesa entonces calcular la probabilidad de que la cola evolucionara desde el nodo interno 5 al terminal 2 en \(10-3=7\) millones de años, dado que en el nodo ancestral 5, el estado del caracter corresponde a la ausencia de cola, es decir:

\[P(X^2(10)=1|X^5(3)=0)=P_{01}(10-3)=p_{01}(7)\]

Esta ecuación que acabo de escribir se denomina probabilidad de transición de evolucionar de 0 a 1 en 7 unidades de tiempo. En general, podríamos escribir todas las probabilidades de transición en una matriz.

\[P(t)= \left( {\begin{array}{cc} p_{00}(t) & p_{01}(t) \\ p_{10}(t) & p_{11}(t) \\ \end{array} } \right)\]

Esta matriz está organizada de tal manera que el primer subíndice de \(p\) es la fila y representa la condición o valor inicial del rasgo (también conocido como “de dónde sale”), y el segundo subíndice es el valor final del rasgo (también conocido como “de dónde llega”). La condición en la probabilidad hace una suposición importante “dado que el ancestro común más reciente no tenía cola”. Además, las filas suman 1, por ejemplo, \(p_{00}(t)+p_{01}(t)=1\), ya que, dado que el rasgo estaba en el estado 0, podría permanecer en el mismo 0 o evolucionar a 1 con probabilidad 1.

Cálculo de las probabilidades de transición

Las probabilidades de transición serían una excelente opción si pudiéramos calcularlas. Representan la probabilidad de cambio evolutivo considerando la incertidumbre. Lamentablemente, calcular estas probabilidades es muy difícil para los modelos Mkn, especialmente aquellos con un gran número de estados discretos.

Vamos a realizar un recorrido utilizando una nueva herramienta matemática para definir modelos Mkn (o cualquier modelo de Markov con rasgos discretos) llamada tasas de transición y matriz Q, y luego volveremos a identificar cómo se utilizan las tasas de transición para definir las probabilidades de transición.

Tasas de transición

Dado que no podemos calcular directamente las probabilidades, proponemos utilizar las tasas de transición. Como su nombre indica, las tasas de transición representan la cantidad de cambio a lo largo del tiempo. Equivalentes a comprender la velocidad a la que se conduce (kilómetros por hora), las tasas de transición nos indican la cantidad de cambios que esperamos de un valor de estado al siguiente en un millón de años de evolución (nota: las unidades de tiempo dependerán de la datación del árbol; si el árbol filogenético está calibrado en el tiempo, las unidades son millones de años; pero si se utiliza un árbol resultante de una estimación de reloj molecular, las unidades son tasas de cambio molecular).

Las tasas que utilizaremos son tasas de transición instantáneas, lo que significa que los cambios de un estado al siguiente ocurren en un intervalo de tiempo infinitesimalmente pequeño. En matemáticas, cuando usamos la palabra infinitesimal, significa que estamos definiendo una derivada (o límite cuando el tiempo tiende a 0), por lo que la relación formal entre una probabilidad de transición y una tasa de transición es:

\[\lim_{h\to 0} P(X(t+h)=j|X(t)=j)=\lim_{h \to 0} p_{ij}(h)=q_{ij}\]

Ahora, volvamos a nuestro ejemplo simplificado con dos estados: presencia y ausencia de colas. Para definir un modelo Mkn, solo necesitaríamos dos tasas de transición: \(q_{01}\) y \(q_{10}\). Con frecuencia, los modelos Mkn se representan mediante diagramas de círculo y flecha, como en la Fig. 4.
Figura 4. Diagrama de círculo y flecha para un modelo Mkn con dos estados y dos tasas de transición
Figura 4. Diagrama de círculo y flecha para un modelo Mkn con dos estados y dos tasas de transición

Una segunda forma, más intuitiva para los biólogos de comprender las tasas de transición es su relación con los tiempos de espera. Para las cadenas de Markov, el tiempo de espera entre la transición de un estado a otro tiene una distribución de probabilidad exponencial con media \(1/q_{01}\). Observa el ejemplo de distribución exponencial a continuación. Cuando la tasa de transición \(q_{01}=0.5\), el tiempo de espera para evolucionar de 0 a 1 es, en promedio, \(\frac{1}{0,5}=2\) millones de años, y cuando el tiempo de espera \(q_{10}=0,03\), el tiempo de espera para la transición de 1 a 0 es, en promedio, \(\frac{1}{0.03}=33.3\) millones de años.

library(ggplot2)
library(reshape2)
q_01=0.5
q_10=0.03

waiting_times <- data.frame(wait_01=rexp(1000, rate=q_01),wait_10=rexp(1000, rate=q_10))

data=melt(waiting_times)

ggplot(data,aes(x=value, fill=variable)) +
  geom_density(alpha=0.5)+
  coord_cartesian(xlim=c(0, 50))+
  geom_vline(xintercept = 1/q_01,color="red")+
  geom_vline(xintercept = 1/q_10,color="blue")+
  labs(x="waiting times",y="Exponential probability distribution")+
  theme_classic()

Matriz de transición

Ahora que comprendemos bien qué son las tasas de transición, podemos definir la siguiente herramienta matemática necesaria para definir una cadena de Markov para un caracter discreto. Esta es la matriz de transición para dos estados, llamada matriz Q.

\[Q= \left( {\begin{array}{cc} - q_{01} & q_{01} \\ q_{10} & -q_{10} \\ \end{array} } \right)\]

Cada renglón de la matriz Q suma 0 (es decir, \(q_{10}+ (-q_{10})=0\)), la razón de esto es porque la matriz Q es simplemente la derivada de la matriz de transición de probabilidad \(P\) que definimos anteriormente.

\[P(t)=e^{Qt}\] \[\frac{d P(t)}{dt}= e^{Qt}Q=P(t)Q\]

Volviendo de la matriz Q a la matriz de probabilidades de transición

Dada la matriz Q, basta con exponenciarla y multiplicarla por el tiempo t para obtener las probabilidades de transición. Este paso es complejo porque la exponencial de una matriz debe calcularse mediante aproximaciones. Los matemáticos han dedicado mucho tiempo y esfuerzo a comprender la exponencial de matrices para calcularlas numéricamente. Para rasgos con más de dos estados, es imposible obtener una solución analítica, pero en nuestro pequeño ejemplo, la matriz de probabilidad de transición se ve así.

\[P(t)=\frac{1}{q_{01}+q_{10}} \left( {\begin{array}{cc} q_{10}+q_{01}e^{-(q_{01}+q_{10})t} & q_{01}-q_{01}e^{-(q_{01}+q_{10})t} \\ q_{10}-q_{10}e^{-(q_{01}+q_{10})t} & q_{01}+q_{10}e^{-(q_{01}+q_{10})t} \\ \end{array} } \right)\]

#Example of calculating the transition probabilities from transition rates for a model with two states
# time is 1 million years. 


prob.mat<-function(q_01,q_10,t){
  exp.val<-exp(-(q_01+q_10)*t)
  probabilities<-1/(q_01+q_10)*matrix(c(q_10+q_01*exp.val, q_10-q_10*exp.val, q_01-q_01*exp.val, q_01+q_10*exp.val),nrow=2)
  return(probabilities)
}

(P.mat<- prob.mat(q_01=0.5, q_10=0.03, t=1))
##            [,1]      [,2]
## [1,] 0.61189148 0.3881085
## [2,] 0.02328651 0.9767135

Como se puede observar, este ya es un conjunto de ecuaciones muy complejo; por lo tanto, al definir modelos para un gran número de estados, siempre opta por definir una matriz Q en lugar de manejar la matriz P. En el software, la exponencial se calcula numéricamente de la siguiente manera.

#install.packages("expm")
library(expm)
(Q.mat<-matrix(c(-0.5,0.03,0.5,-0.03), nrow=2))
##       [,1]  [,2]
## [1,] -0.50  0.50
## [2,]  0.03 -0.03
(P.mat<-expm(Q.mat*1))
##            [,1]      [,2]
## [1,] 0.61189148 0.3881085
## [2,] 0.02328651 0.9767135