{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1968217e",
   "metadata": {},
   "source": [
    "# Clase 2 — Análisis transitorio de CTMC\n",
    "# Ráfagas transcripcionales de un gen (sistema univariado con estado absorbente)\n",
    "\n",
    "**Fuente.** Radulescu, O. et al. (2024). *\"Identifying Markov chain models from time-to-event data: an algebraic approach\"*.\n",
    "\n",
    "**Contexto del paper.** El paper modela la expresión génica como una CTMC: el promotor de un gen alterna entre varios **estados no observables** (configuraciones de encendido/apagado) hasta alcanzar un estado que dispara la transcripción — un evento observable (\"ráfaga\" o *burst* transcripcional). El tiempo entre observaciones sigue una **distribución phase-type**, que es justamente la distribución del tiempo hasta la absorción en la CTMC subyacente (su modelo $M9$, Figura 2). Hoy recorremos ese modelo paso a paso, con foco en las **probabilidades transientes** $\\pi(t)$: ¿qué tan probable es que el gen ya haya producido una ráfaga hacia el instante $t$?\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c71de1cb",
   "metadata": {},
   "source": [
    "## 1. Descripción del sistema\n",
    "\n",
    "El promotor de un gen puede estar en una de dos configuraciones \"cerradas\" (inactivas) o en una configuración \"abierta\" (poised), desde la cual puede dispararse una ráfaga de transcripción. Una vez ocurre la ráfaga, se considera el evento observado (para el análisis del tiempo hasta el evento, tratamos ese estado como absorbente).\n",
    "\n",
    "**Variable de estado.** $X(t)$: configuración del promotor del gen en el instante $t$ (en minutos).\n",
    "\n",
    "## 2. Espacio de estados\n",
    "\n",
    "$$S=\\{1,2,3,4\\}$$\n",
    "- Estado 1: promotor cerrado, configuración A.\n",
    "- Estado 2: promotor cerrado, configuración B.\n",
    "- Estado 3: promotor abierto / poised (puede disparar la ráfaga o volver a cerrarse).\n",
    "- Estado 4: ráfaga transcripcional (evento observable, **estado absorbente** para este análisis).\n",
    "\n",
    "## 3. Datos / parámetros\n",
    "\n",
    "Tasas ilustrativas (del orden típico de la cinética de activación de promotores, en 1/minuto):\n",
    "- $k_1=0.05$: tasa de apertura desde el estado 1 hacia el estado 3.\n",
    "- $k_2=0.08$: tasa de apertura desde el estado 2 hacia el estado 3.\n",
    "- $k_3=0.02$: tasa de cierre desde el estado 3 hacia el estado 1.\n",
    "- $k_4=0.03$: tasa de cierre desde el estado 3 hacia el estado 2.\n",
    "- $k_5=0.5$: tasa de disparo de la ráfaga desde el estado 3 (mucho más rápida que las tasas de apertura/cierre)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4aeb3b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# importar librerías\n",
    "import numpy as np\n",
    "from jmarkov.ctmc import ctmc\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb8837b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 1: define las tasas k1..k5\n",
    "k1, k2, k3, k4, k5 = None, None, None, None, None\n",
    "print(f'k1={k1}, k2={k2}, k3={k3}, k4={k4}, k5={k5}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d99b2a12",
   "metadata": {},
   "source": [
    "## 4. Vector de estados\n",
    "\n",
    "$$\\text{estados}=[1,2,3,4]$$\n",
    "\n",
    "## 5. Vector de condiciones iniciales\n",
    "\n",
    "El gen arranca con el promotor cerrado en la configuración 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d8753f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 2: vector de condiciones iniciales alpha\n",
    "alpha = np.array([])\n",
    "print(f'alpha={alpha}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ef77d3f",
   "metadata": {},
   "source": [
    "## 6. Gráfico de tasas\n",
    "\n",
    "A diferencia de una cadena puramente secuencial, aquí el estado 3 tiene **tres** salidas posibles (puede volver a 1, volver a 2, o disparar la ráfaga hacia 4):\n",
    "\n",
    "```\n",
    "   1  ---k1--->  3  ---k3--->  1   (regresa)\n",
    "                 |\n",
    "   2  ---k2--->  3  ---k4--->  2   (regresa)\n",
    "                 |\n",
    "                 +---k5--->  4  (absorbente: rafaga)\n",
    "```\n",
    "\n",
    "## 7. Tasas de transición\n",
    "\n",
    "| Desde | Hacia | Tasa |\n",
    "|---|---|---|\n",
    "| 1 | 3 | $k_1$ |\n",
    "| 2 | 3 | $k_2$ |\n",
    "| 3 | 1 | $k_3$ |\n",
    "| 3 | 2 | $k_4$ |\n",
    "| 3 | 4 | $k_5$ |\n",
    "\n",
    "## 8. Matriz generadora\n",
    "\n",
    "Esta es exactamente la estructura del modelo $M9$ del paper (su Ecuación 2):\n",
    "\n",
    "$$Q=\\begin{pmatrix}-\\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c7ae7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 3: construye la matriz generadora Q (4x4)\n",
    "Q = np.array([\n",
    "   \n",
    "])\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q.sum(axis=1), 0), 'cada fila de Q debe sumar 0'\n",
    "print('OK: Q es una matriz generadora válida (filas suman 0).')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dbb7d87",
   "metadata": {},
   "source": [
    "## 9. Tasa total de salida de cada estado\n",
    "\n",
    "$q_i=-Q_{ii}$. El estado 3 tiene la tasa de salida más alta, porque de ahí puede ocurrir cualquiera de tres transiciones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d71a738",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 4: tasa total de salida de cada estado\n",
    "q = None  # usa -np.diag(Q)\n",
    "for s, qi in zip([1, 2, 3, 4], q):\n",
    "    print(f'estado {s}: q = {qi:.4f} por minuto')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80af290b",
   "metadata": {},
   "source": [
    "## 10. Tiempos de permanencia\n",
    "\n",
    "$E[T_i]=1/q_i$. El estado 3 es, en promedio, el más breve (por su alta tasa total de salida)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f7f8171",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 5: tiempo medio de permanencia en los 3 estados transientes\n",
    "for s, qi in zip([1, 2, 3], q[:3]):\n",
    "    print(f'estado {s}: E[T] = {None} minutos')  # usa 1/qi\n",
    "print('estado 4: E[T] = infinito (absorbente)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90eff78d",
   "metadata": {},
   "source": [
    "## 11. Probabilidades de la cadena embebida\n",
    "\n",
    "Desde 1 y desde 2 solo hay un destino posible ($P=1$ hacia el estado 3). Desde el estado 3, en cambio, hay tres posibles destinos, con probabilidades proporcionales a sus tasas:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9350e263",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 6: construye la matriz de la cadena embebida P\n",
    "n = 4\n",
    "P = np.zeros((n, n))\n",
    "for i in range(n):\n",
    "    if q[i] > 0:\n",
    "        for j in range(n):\n",
    "            if i != j:\n",
    "                P[i, j] = None  #TODO \n",
    "print(f'P (cadena embebida)=\\n{np.round(P,4)}')\n",
    "print(f'P_3,4 (prob. de disparar la rafaga al salir de 3) = {P[2,3]:.4f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "606b081c",
   "metadata": {},
   "source": [
    "## 12. Probabilidades transientes $\\pi(t)$\n",
    "\n",
    "Con `jmarkov`, evaluamos $\\pi(t)=\\alpha\\,e^{Qt}$ para varios tiempos, sin resolver la ecuación diferencial a mano."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38acbe92",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 7: construye el modelo y evalúa pi(t) para varios tiempos\n",
    "mc =  #CONTRUYA LA CTMC con la matriz generadora Q\n",
    "tiempos = [1, 5, 10, 20, 40, 80]\n",
    "for t in tiempos:\n",
    "    pi_t = None  # calcule las probabilidad transientes\n",
    "    print(f'pi({t} min)={np.round(pi_t,4)}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b73d1ac6",
   "metadata": {},
   "source": [
    "## 13. Función de supervivencia (aún sin ráfaga)\n",
    "\n",
    "Esta es exactamente $S(t)=\\mathbb{P}[T_s>t]$ del paper (su Ecuación 1): la probabilidad de que **todavía no** haya ocurrido la ráfaga en el instante $t$.\n",
    "\n",
    "$$S(t)=1-\\pi_4(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927362ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 8: calcula S(t) para cada uno de los mismos tiempos\n",
    "for t in tiempos:\n",
    "    pi_t = ... #calcule las pribabilidades transientes\n",
    "    S_t = None\n",
    "    print(f'S({t} min) = {S_t:.4f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "682b7573",
   "metadata": {},
   "source": [
    "## 14. Supervivencia en el tiempo (gráfico)\n",
    "\n",
    "Evalúa $S(t)$ en una malla fina y grafica cómo decae — esta curva es la distribución phase-type del paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c3657a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 9: evalúa S(t) en una malla fina y grafica\n",
    "t_fino = np.linspace(0, 80, 200)\n",
    "S_fino = None  # arreglo con S(t) para cada t en t_fino\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4.5))\n",
    "ax.plot(t_fino, S_fino, color='#4C72B0', lw=2)\n",
    "ax.axhline(0.5, color='#C44E52', ls='--', lw=1, label='S(t) = 0.5')\n",
    "ax.set_xlabel('t (minutos)')\n",
    "ax.set_ylabel('S(t)')\n",
    "ax.set_title('Probabilidad de no haber disparado la ráfaga vs. tiempo')\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09a47e38",
   "metadata": {},
   "source": [
    "## 15. Distribución estacionaria y medidas de desempeño\n",
    "\n",
    "Como el estado 4 es absorbente, $\\pi(t)\\to(0,0,0,1)$. El tiempo medio hasta la primera ráfaga (tiempo medio hasta la absorción) se obtiene resolviendo $-Q_{TT}\\,\\mathbf{m}=\\mathbf{1}$ con $Q_{TT}$ la submatriz de los 3 estados transientes — a diferencia del caso de la fachada (Clase 2, Caso B anterior), aquí el sistema **no** es una simple cadena secuencial: el estado 3 puede regresar a 1 o 2 antes de disparar la ráfaga, así que el tiempo medio no es simplemente una suma de $1/q_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67c722e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 10: distribución estacionaria y tiempo medio hasta la primera ráfaga\n",
    "ss = None  # usa mc.steady_state()\n",
    "print(f'pi_estacionaria = {np.round(ss,4)}')\n",
    "\n",
    "Q_TT = Q[:3, :3]\n",
    "m = None  # usa np.linalg.solve(-Q_TT, np.ones(3))\n",
    "print(f'tiempo medio hasta la primera rafaga, partiendo de cada estado: {np.round(m,2)} minutos')\n",
    "print(f'tiempo medio hasta la primera rafaga (arrancando en estado 1) = {m[0]:.2f} minutos')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90ec2ffd",
   "metadata": {},
   "source": [
    "## 16. Costeo: costo esperado a partir de probabilidades\n",
    "\n",
    "Supongamos que mantener el promotor en cada configuración tiene un **costo metabólico** por minuto para la célula (energía de mantener la cromatina en cada estado), y que el disparo de la ráfaga (estado 4) no acumula más costo (ya se considera el evento observado):\n",
    "\n",
    "| Estado | Costo $c_i$ (por minuto) |\n",
    "|---|---|\n",
    "| 1 (cerrado A) | 1.0 |\n",
    "| 2 (cerrado B) | 1.5 |\n",
    "| 3 (abierto/poised) | 3.0 |\n",
    "| 4 (ráfaga, absorbente) | 0 |\n",
    "\n",
    "**Costo esperado instantáneo (transitorio), usando $\\pi(t)$:**\n",
    "\n",
    "$$E[C(t)]=\\sum_{i=1}^{4} c_i\\,\\pi_i(t) = \\pi(t)\\cdot\\mathbf{c}$$\n",
    "\n",
    "**Costo esperado total hasta la primera ráfaga.** Esta es la generalización natural del tiempo medio hasta la absorción (sección 15): en vez de resolver $-Q_{TT}\\mathbf{m}=\\mathbf{1}$, se resuelve con el vector de costos $\\mathbf{c}_T=(c_1,c_2,c_3)$ en el lado derecho:\n",
    "\n",
    "$$-Q_{TT}\\,\\mathbf{m}_c=\\mathbf{c}_T \\quad\\Longrightarrow\\quad \\mathbf{m}_c=(-Q_{TT})^{-1}\\mathbf{c}_T$$\n",
    "\n",
    "$m_c(i)$ es el costo metabólico total esperado, partiendo del estado $i$, hasta que ocurre la ráfaga."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33c8d2ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 11: costo esperado instantaneo y costo esperado total hasta la primera rafaga\n",
    "c = np.array([1.0, 1.5, 3.0, 0.0])\n",
    "\n",
    "for t in [1, 5, 10, 20]:\n",
    "    pi_t = ...   # calcule las pribabilidades transientes\n",
    "    E_C_t = None  # producto punto pi_t @ c\n",
    "    print(f'E[C({t} min))] = {E_C_t:.4f}')\n",
    "\n",
    "c_T = c[:3]\n",
    "m_c = None  # usa np.linalg.solve(-Q_TT, c_T)\n",
    "print(f'\\ncosto esperado total hasta la primera rafaga, partiendo de cada estado: {np.round(m_c,3)}')\n",
    "print(f'costo esperado total (arrancando en estado 1) = {m_c[0]:.3f}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert m_c[0] > 0\n",
    "print('OK: el costo esperado total hasta la rafaga es positivo.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1ea2814",
   "metadata": {},
   "source": [
    "## Conclusión\n",
    "\n",
    "_¿Por qué el tiempo medio hasta la primera ráfaga no es simplemente $1/q_1+1/q_3$ (como sí lo era en el caso de la fachada)? ¿Qué papel juega la posibilidad de que el promotor \"regrese\" del estado 3 a los estados 1 o 2 antes de disparar la ráfaga?_"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}