Modelos Biomatemáticos I. Notas 6 (parte 4) — MATERIAL EN REVISIÓN

Por Mariana Paulin

6.6 Aplicaciones de matrices en modelos biológicos

Sustitución de nucleótidos en evolución

La evolución molecular analiza cómo cambian las secuencias de ADN y proteínas a lo largo del tiempo, donde uno de los procesos más importantes es la sustitución de nucleótidos, que consiste en el reemplazo de una base nitrogenada (adenina, citosina, guanina o timina) por otra dentro de la cadena de ADN. Gracias a estas mutaciones puntuales se da buena parte de la diversidad genética y conforman la base para estudiar la divergencia entre especies, reconstruir árboles filogenéticos y comprender los patrones evolutivos a nivel molecular.
Para modelar estos cambios, se utilizan procesos estocásticos, en particular cadenas de Markov, que permiten describir la probabilidad de transición de una base a otra  en un intervalo de tiempo mediante matrices de transición. (Yang, 2014)

Modelo matricial de sustitución de nucleótidos

Para modelar cómo cambia un nucleótido a otro a lo largo del tiempo, se supone que:
• El estado futuro de una base (por ejemplo, si una A mutará a G, C o T) depende únicamente del estado actual.
• Las probabilidades de sustitución entre bases son constantes en el tiempo (modelo homogéneo).
• El sistema se representa mediante una matriz de transición de tamaño 4×4, donde las filas corresponden al estado actual (base presente) y las columnas al estado futuro (base que puede ser).

Una matriz de transición $P$ tiene la forma

$P = \begin{bmatrix} p_{AA} & p_{AC} & p_{AG} & p_{AT} \\ p_{CA} & p_{CC} & p_{CG} & p_{CT} \\ p_{GA} & p_{GC} & p_{GG} & p_{GT} \\ p_{TA} & p_{TC} & p_{TG} & p_{TT} \\ \end{bmatrix}$

donde $p_{ij}$​ representa la probabilidad de que la base $i$ cambie a la base $j$ en un tiempo determinado.

Una propiedad de estas matrices es que cada fila debe sumar 1:

$$\sum_{j} p_{ij} = 1 \quad \forall i \in \{A,C,G,T\}$$

Esto es porque las probabilidades de cambiar de la base $i$ a cualquiera de las cuatro bases (incluyendo la misma) cubren todas las posibilidades y, por lo tanto, deben sumar 1.  Las matrices que cumplen esto se llaman matrices estocásticas por filas. (Rodríguez et al., 1990)

Ejemplo
Supongamos que una población de virus de ARN muta durante la replicación y que la sustitución entre nucleótidos se representa con la siguiente matriz de transición:

$M = \begin{bmatrix} 0.85 & 0.05 & 0.05 & 0.05 \\ 0.10 & 0.80 & 0.05 & 0.05 \\ 0.10 & 0.05 & 0.80 & 0.05 \\ 0.15 & 0.05 & 0.05 & 0.75 \\ \end{bmatrix}$

De modo que la fila 1 corresponde a A, la fila 2 a C, la 3 a G y la 4 a T. Entonces, por ejemplo, un nucleótido A tiene una probabilidad de 85 % de seguir siendo A, y 5 % de mutar a C, G o T.
Además, se sabe que en la secuencia hay 1000 nucleótidos distribuidos en un estado inicial como

$\mathbf{x}(0) = \begin{bmatrix} 300 \\ 200 \\ 300 \\ 200 \end{bmatrix}$

Se desea estimar la frecuencia esperada de cada nucleótido después de una generación.

Respuesta esperada
Como las frecuencias pueden ser expresadas  como proporciones, podemos dividir entre el total

$\mathbf{x}(0) = \begin{bmatrix} 0.3 \\ 0.2 \\ 0.3 \\ 0.2 \end{bmatrix}$

Para obtener el estado después de una generación, multiplicamos

$\mathbf{x}(1) = \mathbf{x}(0) \cdot M$

$\Rightarrow \mathbf{x}(1) = \begin{bmatrix} 0.3 \\ 0.2 \\ 0.3 \\ 0.2 \end{bmatrix} \cdot \begin{bmatrix} 0.85 & 0.05 & 0.05 & 0.05 \\ 0.10 & 0.80 & 0.05 & 0.05 \\ 0.10 & 0.05 & 0.80 & 0.05 \\ 0.15 & 0.05 & 0.05 & 0.75 \\ \end{bmatrix}$

• Componente A:
$0.3(0.85) + 0.2(0.10) + 0.3(0.10) + 0.2(0.15)$
$ = 0.255 + 0.02 + 0.03 + 0.03 = 0.335$
• Componente C:
$0.3(0.05) + 0.2(0.80) + 0.3(0.05) + 0.2(0.05)$
$ = 0.015 + 0.16 + 0.015 + 0.01 = 0.2$
• Componente G:
$0.3(0.05) + 0.2(0.05) + 0.3(0.80) + 0.2(0.05)$
$ = 0.015 + 0.01 + 0.24 + 0.01 = 0.275$
• Componente T:
$0.3(0.05) + 0.2(0.05) + 0.3(0.05) + 0.2(0.75)$
$ = 0.015 + 0.01 + 0.015 + 0.15 = 0.19$

Entonces

$\mathbf{x}(1) = \begin{bmatrix} 0.335 \\ 0.2 \\ 0.275 \\ 0.19 \end{bmatrix}$

Después de una generación de mutaciones, se espera que la proporción de nucleótidos cambie ligeramente: aumenta la proporción de A, de 30 % a 33.5 %; disminuye G, de 30 % a 27.5 %; C y T se mantienen igual o disminuyen levemente.

Ejemplo
Supongamos que estamos estudiando la evolución de una región específica del ADN en una población de bacterias. Queremos entender cómo las bases nitrogenadas cambian a lo largo de generaciones bajo un modelo simple donde todas las sustituciones son igualmente probables, excepto que la base puede permanecer igual con una probabilidad mayor.
Un modelo clásico para este caso es el modelo de Jukes-Cantor, que asume que todas las sustituciones de una base a otra tienen la misma probabilidad y el proceso es homogéneo y estacionario en el tiempo.

Supongamos que para un intervalo de tiempo $t$, la matriz de transición es:

$M = \begin{bmatrix} 0.7 & 0.1 & 0.1 & 0.1 \\ 0.1 & 0.7 & 0.1 & 0.1 \\ 0.1 & 0.1 & 0.7 & 0.1 \\ 0.1 & 0.1 & 0.1 & 0.7 \\ \end{bmatrix}$​​

Esto significa, por ejemplo, que si en la generación actual la base es A, en la siguiente generación hay un 70 % de probabilidad de que siga siendo A, y un 10 % de cambiar a C, G o T respectivamente.
Si comenzamos con una población en la que la base en una posición es A con probabilidad 1 (seguro), ¿cuál será la distribución de bases después de un intervalo de tiempo?
La distribución inicial es

$\mathbf{x}(0) = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix}$

Para obtener la distribución en la siguiente generación, multiplicamos por la matriz $M$

$\mathbf{x}(1) = \mathbf{x}(0) \cdot M = \begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \end{bmatrix} \cdot \begin{bmatrix} 0.7 & 0.1 & 0.1 & 0.1 \\ 0.1 & 0.7 & 0.1 & 0.1 \\ 0.1 & 0.1 & 0.7 & 0.1 \\ 0.1 & 0.1 & 0.1 & 0.7 \\ \end{bmatrix} = \begin{bmatrix} 0.7 \\ 0.1 \\ 0.1 \\ 0.1 \end{bmatrix}$

Por lo que después de ese intervalo, la base tiene un 70 % de seguir siendo A y un 10 % de ser cada una de las otras bases.

Ejercicio
Se tiene la siguiente matriz de transición para las sustituciones de nucleótidos en una población de virus:

$M = \begin{bmatrix} 0.6 & 0.2 & 0.1 & 0.1 \\ 0.1 & 0.7 & 0.1 & 0.1 \\ 0.2 & 0.1 & 0.6 & 0.1 \\ 0.1 & 0.2 & 0.1 & 0.6 \\ \end{bmatrix}$

Si inicialmente la base es G con probabilidad 1, ¿cuál será la distribución de bases después de un intervalo de tiempo?
Supón que el vector de estado inicial es $\mathbf{x}(0) = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix}$

Respuesta esperada
Multiplicamos

$\mathbf{x}(1) = \mathbf{x}(0) \cdot M = \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix} \cdot \begin{bmatrix} 0.6 & 0.2 & 0.1 & 0.1 \\ 0.1 & 0.7 & 0.1 & 0.1 \\ 0.2 & 0.1 & 0.6 & 0.1 \\ 0.1 & 0.2 & 0.1 & 0.6 \\ \end{bmatrix} = \begin{bmatrix} 0.2 \\ 0.1 \\ 0.6 \\ 0.1 \end{bmatrix}$

Luego, después del intervalo, la base será A con probabilidad 0.2, C con 0.1, G con 0.6, y T con 0.1.

Modelo SIR en epidemiología 

El modelo SIR es un modelo matemático que describe cómo se propaga una enfermedad infecciosa en una población. Se llama SIR porque divide a la población en tres clases: susceptibles, individuos sanos que pueden enfermar si tienen contacto con una persona infectada; infectados, individuos que actualmente tienen la enfermedad y pueden contagiarla; recuperados, individuos que ya no están enfermos y no pueden contagiar ni contagiarse.

Los cambios se modelan con ecuaciones diferenciales, que indican cómo varía cada grupo en función del tiempo.

$\begin{aligned} \frac{dS}{dt} &= -\beta S I \\ \frac{dI}{dt} &= \beta S I-\gamma I \\ \frac{dR}{dt} &= \gamma I \end{aligned}$

donde
• $\beta$ es la tasa de transmisión (contacto entre susceptibles e infectados),
• $\gamma$ es la tasa de recuperación.

Este modelo describe el cambio de una epidemia en el tiempo, siendo no lineal debido al término $\beta S I$. No obstante, se puede linealizar el sistema y representarlo en forma matricial.
Para simular paso a paso (por ejemplo, día a día), sin usar derivadas directamente, usaremos el método de aproximaciones de Eüler: dada una unidad de tiempo $\Delta t$ tenemos que

$\begin{aligned} S(t + \Delta t) &\approx S(t) + \Delta t \cdot \frac{dS}{dt}(t), \\ I(t + \Delta t) &\approx I(t) + \Delta t \cdot \frac{dI}{dt}(t), \\ R(t + \Delta t) &\approx R(t) + \Delta t \cdot \frac{dR}{dt}(t). \end{aligned}$

Sustituimos las derivadas del modelo SIR y queda que

$\begin{aligned} S(t + \Delta t) &\approx S(t)-\Delta t \cdot \beta S(t) I(t), \\ I(t + \Delta t) &\approx I(t) + \Delta t \cdot ( \beta S(t) I(t)-\gamma I(t) ), \\ R(t + \Delta t) &\approx R(t) + \Delta t \cdot \gamma I(t). \end{aligned}$

Ahora, para formular esto como una multiplicación matricial del tipo

$\mathbf{x}(t + \Delta t) = A \cdot \mathbf{x}(t),$

donde

$\mathbf{x}(t) = \begin{bmatrix} S(t) \\ I(t) \\ R(t) \end{bmatrix},$

queremos una matriz $A$ tal que al multiplicar por $\mathbf{x}(t)$ resulten las expresiones anteriores. Una manera es asumir que durante ese pequeño intervalo de tiempo $\Delta t$, el término no lineal $\beta S(t) I(t)$ se puede “incorporar” en la matriz, si mantenemos $I(t)$ como constante en ese momento para la construcción de la matriz (Suryanto, 2014).

Tomamos $\Delta t = 1$ para tener que

$\begin{aligned} S(t+1) &= (1-\beta I(t)) \, S(t), \\ I(t+1) &= \beta I(t) \, S(t) + (1-\gamma) I(t), \\ R(t+1) &= \gamma I(t) + R(t). \end{aligned}$

Esto corresponde a decir que

$\begin{bmatrix} S(t+1) \\ I(t+1) \\ R(t+1) \end{bmatrix} = \begin{bmatrix} 1-\beta I(t) & 0 & 0 \\ \beta I(t) & 1-\gamma & 0 \\ 0 & \gamma & 1 \end{bmatrix} \cdot \begin{bmatrix} S(t) \\ I(t) \\ R(t) \end{bmatrix}​$

donde $I_0$ representa una cantidad inicial de infectados en la primera fase. Esta matriz permite ver cómo cambian los compartimentos en pequeños intervalos discretos  (Bonhoeffer & Schafroth, 2020).

Ejemplo
Supongamos que en una población cerrada de 1000 personas, 780 son susceptibles, 220 están infectadas y 0 se han recuperado; y la enfermedad tiene
• Tasa de transmisión $\beta = 0.03$
• Tasa de recuperación $\gamma = 0.1$
¿Cuántas personas habrá en cada grupo al día siguiente?

Respuesta esperada
El estado inicial es

$\mathbf{x}(0) = \begin{bmatrix} S(0) \\ I(0) \\ R(0) \end{bmatrix} = \begin{bmatrix} 0.78 \\ 0.22 \\ 0 \end{bmatrix}$​​

La matriz de transición es

$A = \begin{bmatrix} 1-\beta I(0) & 0 & 0 \\ \beta I(0) & 1-\gamma & 0 \\ 0 & \gamma & 1 \end{bmatrix}$

Sustituimos los valores: $\beta = 0.03,\quad I(0) = 0.22,\quad \gamma = 0.1$

$\Rightarrow A = \begin{bmatrix} 1-0.03(0.22) & 0 & 0 \\ 0.03(0.22) & 1-0.1 & 0 \\ 0 & 0.1 & 1 \end{bmatrix} = \begin{bmatrix} 0.9934 & 0 & 0 \\ 0.0066 & 0.9 & 0 \\ 0 & 0.1 & 1 \end{bmatrix}$

Calculamos el estado al siguiente momento

$\mathbf{x}(1) = A \cdot \mathbf{x}(0) = \begin{bmatrix} 0.9934 & 0 & 0 \\ 0.0066 & 0.9 & 0 \\ 0 & 0.1 & 1 \end{bmatrix} \cdot \begin{bmatrix} 0.78 \\ 0.22 \\ 0 \end{bmatrix} = \begin{bmatrix} 0.77766 \\ 0.20034 \\ 0.022 \end{bmatrix}$​​

Luego $\mathbf{x}(1) = \begin{bmatrix} 0.774 \\ 0.204 \\ 0.022 \end{bmatrix}$

Entonces, después de un día, aproximadamente: 6 personas dejaron de ser susceptibles, hay 16 infectados menos que al inicio y 22 personas se recuperaron.

Ejercicio
Supongamos que en una comunidad de 500 personas, 495 son susceptibles, 5 están infectadas y 0 están recuperadas. Y se tienen las tasas: $\beta = 0.001$ y $\gamma = 0.2$.
a. Escribe el vector de estado inicial
b. Escribe la matriz de transición
c. Calcula cómo estará distribuida la población al día siguiente.

Respuesta esperada
a.

$\mathbf{x}(0) = \begin{bmatrix} S(0) \\ I(0) \\ R(0) \end{bmatrix} = \begin{bmatrix} \frac{495}{500} \\ \frac{5}{500} \\ 0 \end{bmatrix} = \begin{bmatrix} 0.99 \\ 0.01 \\ 0 \end{bmatrix}$

b. Tenemos

$A = \begin{bmatrix} 1-\beta I(0) & 0 & 0 \\ \beta I(0) & 1-\gamma & 0 \\ 0 & \gamma & 1 \end{bmatrix}$​​

Sustituimos los valores: $\beta = 0.001,\quad I(0) = 0.01,\quad \gamma = 0.2$

$\Rightarrow A = \begin{bmatrix} 1-0.001 \cdot 0.01 & 0 & 0 \\ 0.001 \cdot 0.01 & 0.8 & 0 \\ 0 & 0.2 & 1 \end{bmatrix} = \begin{bmatrix} 0.99999 & 0 & 0 \\ 0.00001 & 0.8 & 0 \\ 0 & 0.2 & 1 \end{bmatrix}$​​

c.

$\mathbf{x}(1) = \begin{bmatrix} 0.99999 & 0 & 0 \\ 0.00001 & 0.8 & 0 \\ 0 & 0.2 & 1 \end{bmatrix} \cdot \begin{bmatrix} 0.99 \\ 0.01 \\ 0 \end{bmatrix} = \begin{bmatrix} 0.9809901 \\ 0.0080099 \\ 0.002 \end{bmatrix}$

Si multiplicamos por 500 para obtener los valores absolutos, tendemos que al día siguiente hay aproximadamente 491 individuos susceptibles, 5 individuos infectados y 1 recuperado.

Diagrama cinemático en conducta animal

En etología (el estudio científico del comportamiento animal), es común modelar el comportamiento de un organismo como una secuencia de estados discretos: por ejemplo, buscar alimento, alimentarse, descansar, explorar, etcétera.
Un diagrama cinemático es una representación de estas transiciones de un estado a otro, usando una matriz de transición que permita cuantificar y predecir cómo un animal cambia de comportamiento con el tiempo.
Este enfoque es análogo a un sistema de Markov en tiempo discreto, donde cada estado es una clase o comportamiento, así, cada comportamiento observado en un instante depende sólo del estado anterior y las probabilidades o tasas de transición se representan mediante una matriz cuadrada. (Metz, Dienske, de Jonge, 1983)

Representación matricial

El modelo se escribe como $\mathbf{x}(t+1) = A \cdot \mathbf{x}(t)$ donde
• $\mathbf{x}(t)$ es el vector que representa la proporción (o número absoluto) de individuos en cada estado conductual en el instante $t$,
• $A$ es la matriz de transición donde cada entrada $a_{ij}$​ representa la probabilidad de que un individuo en el estado $j$ en $t$ pase al estado $i$ en $t+1$.
En ese caso, $A$ es una matriz estocástica por columnas.

Ejemplo
Supongamos que un ave muestra tres comportamientos observables que categorizamos en estados: estado 1, buscar alimento; estado 2, comer; estado 3, descansar.
Supongamos que se realiza un registro conductual durante cada hora y se estima la siguiente matriz de transición

$A = \begin{bmatrix} 0.1 & 0.6 & 0.3 \\ 0.6 & 0.2 & 0.2 \\ 0.3 & 0.2 & 0.5 \end{bmatrix}$

Es decir, se observó que

• Desde el estado 2 (comer), el 60 % de las veces el ave pasa a buscar alimento (estado 1), el 20 % se mantiene comiendo, y el 20 % pasa a descansar.
• Desde el estado 3 (descansar), el 50 % de las veces continúa descansando, un 30 % pasa a buscar alimento, y un 20 % a comer.

Supongamos que al tiempo $t = 0$, la distribución comportamental es

$\mathbf{x}(0) = \begin{bmatrix} 0.2 \\ 0.5 \\ 0.3 \end{bmatrix}$

Es decir, el 20 % está buscando alimento, el 50 % está comiendo, y el 30 % está descansando.
Si queremos calcular la distribución en el siguiente instante:

$\mathbf{x}(1) = A \cdot \mathbf{x}(0) = \begin{bmatrix} 0.1 & 0.6 & 0.3 \\ 0.6 & 0.2 & 0.2 \\ 0.3 & 0.2 & 0.5 \end{bmatrix} \cdot \begin{bmatrix} 0.2 \\ 0.5 \\ 0.3 \end{bmatrix}$

$\begin{aligned} x_1(1) &= 0.1(0.2) + 0.6(0.5) + 0.3(0.3) = 0.02 + 0.30 + 0.09 = 0.41 \\ x_2(1) &= 0.6(0.2) + 0.2(0.5) + 0.2(0.3) = 0.12 + 0.10 + 0.06 = 0.28 \\ x_3(1) &= 0.3(0.2) + 0.2(0.5) + 0.5(0.3) = 0.06 + 0.10 + 0.15 = 0.31 \end{aligned}$

$\Rightarrow \mathbf{x}(1) = \begin{bmatrix} 0.41 \\ 0.28 \\ 0.31 \end{bmatrix}$​​

Podemos decir que después de una unidad de tiempo, el comportamiento dominante se ha desplazado hacia Buscar alimento.

Comportamiento a largo plazo

En el estudio del comportamiento animal, además de conocer cómo cambia la conducta momento a momento, también podemos entender cómo se estabiliza un patrón de comportamiento a largo plazo. Por ejemplo, podemos preguntarnos ¿cuánto tiempo pasa un animal, en promedio, comiendo o descansando si lo observamos durante muchas horas o muchos días?
Para responder esto, debemos estudiar el comportamiento asintótico del sistema, si el modelo de transición de comportamientos se considera como una cadena de Markov en tiempo discreto. Es decir, si suponemos que el comportamiento de un animal se modela por una matriz de transición $A$ que satisface

$\mathbf{x}(t+1) = A \cdot \mathbf{x}(t)$

entonces, bajo ciertas condiciones, existe un vector estacionario $\mathbf{v}$ tal que

$$A \cdot \mathbf{v} = \mathbf{v}, \quad \text{con} \quad \sum_{i} v_i = 1$$

Este vector $\mathbf{v}$ representa una distribución invariante, o sea que si el comportamiento del animal en el tiempo $t$ está dado por $\mathbf{x}(t) = \mathbf{v}$, entonces en todos los tiempos futuros se mantendrá igual: $\mathbf{x}(t+1) = \mathbf{v}, \mathbf{x}(t+2) = \mathbf{v}, \ldots$ (Schinazi, 1999).
Entonces, si lo vemos desde el punto de vista etológico, el vector estacionario representará la proporción de tiempo que un animal pasará en promedio en cada estado conductual si se lo observa por un periodo suficientemente largo. Esto no implica que el animal esté congelado en esos estados, sino que, a lo largo del tiempo su comportamiento oscilará entre ellos con una proporción promedio estable.
Ahora, como se mencionó anteriormente, se requieren ciertas condiciones para que exista y se alcance esta distribución estable. No todas las matrices de transición garantizan que exista un comportamiento estacionario alcanzable desde cualquier estado inicial. Para que sí exista y se cumpla

$$\lim_{t \to \infty} \mathbf{x}(t) = \mathbf{v}$$

es necesario que la matriz $A$ sea regular, lo cual en términos prácticos significa que alguna potencia de $A$ tiene todas sus entradas estrictamente positivas. Esto implica que es posible ir de cualquier estado a cualquier otro, y que no hay estados absorbentes, o sea de los que ya no se puede salir.

Para encontrar el vector estacionario $\mathbf{v}$, se resuelve el sistema

$$A \cdot \mathbf{v} = \mathbf{v}, \quad \sum_{i} v_i = 1$$

Que es lo mismo que resolver

$$(A-I)\mathbf{v} = 0, \quad \text{con} \quad \sum_{i} v_i = 1$$

Este sistema suele resolverse algebraicamente o numéricamente. En la práctica, también puede aproximarse aplicando repetidamente la matriz a un estado inicial arbitrario (iteraciones sucesivas de x(t+1)=A⋅x(t)\mathbf{x}(t+1) = A \cdot \mathbf{x}(t)x(t+1)=A⋅x(t)) hasta que se estabilicen los valores.

Ejercicio
Una tortuga exhibe tres estados de comportamiento: explorar (1), comer (2) y dormir (3). Supongamos que con observaciones se obtiene la matriz de transición

$A = \begin{bmatrix} 0.3 & 0.1 & 0.2 \\ 0.4 & 0.6 & 0.3 \\ 0.3 & 0.3 & 0.5 \end{bmatrix}$​​

a. Verifica si $A$ es una matriz estocástica por columnas.
b. Si la distribución inicial es $\mathbf{x}(0) = \begin{bmatrix} 0.3 \\ 0.5 \\ 0.2 \end{bmatrix}$​​, encuentra $\mathbf{x}(1)$.
c. Calcula la distribución estacionaria $\mathbf{v}$.

Respuesta esperada
a. Sabemos que una matriz es estocástica por columnas si cada columna suma 1.
Columna 1: $0.3 + 0.4 + 0.3 = 1$
Columna 2: $0.1 + 0.6 + 0.3 = 1$
Columna 3: $0.2 + 0.3 + 0.5 = 1$
Entonces sí, $A$ es estocástica por columnas.

b. 

$\mathbf{x}(1) = A \cdot \mathbf{x}(0) = \begin{bmatrix} 0.3 & 0.1 & 0.2 \\ 0.4 & 0.6 & 0.3 \\ 0.3 & 0.3 & 0.5 \\ \end{bmatrix} \cdot \begin{bmatrix} 0.3 \\ 0.5 \\ 0.2 \end{bmatrix}$​​

$\mathbf{x_1}(1) = 0.3(0.3) + 0.1(0.5) + 0.2(0.2) = 0.09 + 0.05 + 0.04 = 0.18$
$\mathbf{x_2}(1) = 0.4(0.3) + 0.6(0.5) + 0.3(0.2) = 0.12 + 0.30 + 0.06 = 0.48$
$\mathbf{x_3}(1) = 0.3(0.3) + 0.3(0.5) + 0.5(0.2) = 0.09 + 0.15 + 0.10 = 0.34$
$\Rightarrow \mathbf{x}(1) = \begin{bmatrix} 0.18 \\ 0.48 \\ 0.34 \end{bmatrix}$

Esto quiere decir que, después de un intervalo de tiempo, la mayoría de las tortugas están comiendo (48 %), seguidas de dormir (34 %) y explorar (18 %).

c. Comenzamos calculando $A-I$

$A = \begin{bmatrix} 0.3 & 0.1 & 0.2 \\ 0.4 & 0.6 & 0.3 \\ 0.3 & 0.3 & 0.5 \end{bmatrix}, \quad I = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$

$A-I = \begin{bmatrix} -0.7 & 0.1 & 0.2 \\ 0.4 & -0.4 & 0.3 \\ 0.3 & 0.3 & -0.5 \end{bmatrix}$

Entonces buscamos $\mathbf{v}$ tal que

$(A-I)\mathbf{v} = 0, \quad \text{con } v_1 + v_2 + v_3 = 1$

Escribimos el sistema como ecuaciones para resolverlo

$\begin{cases} -0.7v_1 + 0.1v_2 + 0.2v_3 = 0 \quad \text{(ec1)} \\ 0.4v_1-0.4v_2 + 0.3v_3 = 0 \quad \text{(ec2)} \\ 0.3v_1 + 0.3v_2-0.5v_3 = 0 \quad \text{(ec3)} \\ v_1 + v_2 + v_3 = 1 \quad \text{(ec4)} \end{cases}$

Usamos ec1 y ec2 para resolver por eliminación:
multiplicamos ec1 por 4 y ec2 por 7 para igualar los coeficientes de $v_1$​, obtenemos
ec1: $-2.8v_1 + 0.4v_2 + 0.8v_3 = 0$
ec2: $2.8v_1-2.8v_2 + 2.1v_3 = 0$

Sumamos ambas ecuaciones
$(-2.8 + 2.8)v_1 + (0.4-2.8)v_2 + (0.8 + 2.1)v_3 = 0$
$\Rightarrow -2.4v_2 + 2.9v_3 = 0$

entonces

$$2.9v_3 = 2.4v_2 \Rightarrow v_3 = \frac{2.4}{2.9}v_2 \approx 0.8276v_2$$

Sustituimos $v_3$​ en términos de $v_2$ en ec4 para obtener todos los valores, si tenemos inicialmente ec4: $v_1 + v_2 + v_3 = 1$, entonces
$v_1 + v_2 + 0.8276v_2 = 1$
$v_1 + 1.8276v_2 = 1$
$\Rightarrow v_1 = 1-1.8276v_2$
Este sistema tiene infinitas soluciones, por lo que elegimos a $v_2$ como variable libre y la igualamos a $t$, de modo que

$v_2 = t \Rightarrow v_3 = 0.8276t \Rightarrow v_1 = 1-t-0.8276t = 1-1.8276t$

luego, la familia de soluciones es

$\mathbf{v}(t) = \begin{bmatrix} 1-1.8276t \\ t \\ 0.8276t \end{bmatrix}$​​

Sabemos que para que $\mathbf{v}$ sea una distribución de probabilidad válida, se deben cumplir las condiciones:
• Los valores deben sumar 1 (se cumple por construcción).
• Todos los valores deben ser mayores o iguales que cero.

Observemos que $v_1 = 1-1.8276t$ podría volverse negativa si $t$ es muy grande, entonces tenemos la condición

$$1-1.8276t \geq 0 \Rightarrow t \leq \frac{1}{1.8276} \approx 0.547$$

Entonces, cualquier valor de $t$ en el intervalo $[0, 0.547]$ generará una distribución válida. Dentro del intervalo válido, tomamos arbitrariamente $t=0.4$, de modo que

$t = 0.4 \Rightarrow v_2 = 0.4, \quad v_3 = 0.8276 \cdot 0.4 \approx 0.331, \quad v_1 = 1-0.4-0.331 = 0.269$

Obtenemos el vector

$\mathbf{v} \approx \begin{bmatrix} 0.269 \\ 0.4 \\ 0.331 \end{bmatrix}$​​

Este vector ya suma 1 (aproximadamente) y todos los valores son positivos. Por tanto, es una distribución estacionaria válida.
Esto quiere decir que, a largo plazo, la tortuga pasará aproximadamente 26.9 % del tiempo explorando, 40.0 % del tiempo comiendo, y 33.1 % del tiempo durmiendo.

Ejercicio
Supón que un insecto muestra cuatro estados de comportamientos: buscar alimento (1), comer (2), escapar (3), y reposar (4). Y que la matriz de transición es

$A = \begin{bmatrix} 0.1 & 0.2 & 0.4 & 0.1 \\ 0.7 & 0.3 & 0.1 & 0.2 \\ 0.1 & 0.2 & 0.4 & 0.2 \\ 0.1 & 0.3 & 0.1 & 0.5 \\ \end{bmatrix}$

a. Responde: ¿cuál es la interpretación de $a_{23} = 0.1$?
b. Comprueba que $A$ es estocástica por columnas.
c. Si $\mathbf{x}(0) = \begin{bmatrix} 0.25 \\ 0.25 \\ 0.25 \\ 0.25 \end{bmatrix}$​​, encuentra $\mathbf{x}(1)$.
d. Responde: ¿cuál es la interpretación biológica de una distribución estacionaria en este contexto?

Respuesta esperada
a. $a_{23} = 0.1$ indica que hay una probabilidad del 10 % de que un insecto que está en el estado 3, escapar, pase al estado 2, comer, en el siguiente intervalo de tiempo.
b. 
Columna 1: $0.1 + 0.7 + 0.1 + 0.1 = 1$
Columna 2: $0.2 + 0.3 + 0.2 + 0.3 = 1$
Columna 3: $0.4 + 0.1 + 0.4 + 0.1 = 1$
Columna 4: $0.1 + 0.2 + 0.2 + 0.5 = 1$
Luego, $A$ es estocástica por columnas.
c. 

$\mathbf{x}(1) = A \cdot \mathbf{x}(0) = \begin{bmatrix} 0.1 & 0.2 & 0.4 & 0.1 \\ 0.7 & 0.3 & 0.1 & 0.2 \\ 0.1 & 0.2 & 0.4 & 0.2 \\ 0.1 & 0.3 & 0.1 & 0.5 \end{bmatrix} \cdot \begin{bmatrix} 0.25 \\ 0.25 \\ 0.25 \\ 0.25 \end{bmatrix}$

$\mathbf{x_1}(1) = 0.1(0.25) + 0.2(0.25) + 0.4(0.25) + 0.1(0.25) = 0.025 + 0.05 + 0.10 + 0.025 = 0.20$
$\mathbf{x_2}(1) = 0.7(0.25) + 0.3(0.25) + 0.1(0.25) + 0.2(0.25) = 0.175 + 0.075 + 0.025 + 0.05 = 0.325$
$\mathbf{x_3}(1) = 0.1(0.25) + 0.2(0.25) + 0.4(0.25) + 0.2(0.25) = 0.025 + 0.05 + 0.10 + 0.05 = 0.225$
$\mathbf{x_4}(1) = 0.1(0.25) + 0.3(0.25) + 0.1(0.25) + 0.5(0.25) = 0.025 + 0.075 + 0.025 + 0.125 = 0.25$

$\Rightarrow \mathbf{x}(1) = \begin{bmatrix} 0.20 \\ 0.325 \\ 0.225 \\ 0.25 \end{bmatrix}$​​

Esto quiere decir que en este intervalo, el comportamiento más frecuente será comer (32.5 %), seguido de reposar (25 %), escapar (22.5 %) y buscar alimento (20 %).

d. Sabemos que una distribución estacionaria representa la proporción estable de tiempo que el insecto pasa en cada estado de comportamiento a largo plazo, independientemente del estado inicial. Es decir que se trata del comportamiento promedio que se espera del insecto si se observa durante un periodo largo, donde el patrón de cambio entre estados ya no varía tanto.

[Faltan las referencias]

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.