1 Introducción

El VIH existe en la célula huésped como virus latente o en replicación activa. Muchos aspectos de este proceso se continúan comprendiendo por lo que la modelación estocástica proporciona una herramienta útil al estudiar las condiciones y los desencadenantes que hacen que el virus pase de un estado de latencia a uno de replicación viral activa.

En esta investigación, se generan trayectorias de un modelo con depósito de latencia de infección (Aavani and Allen 2019) que da cuenta del rol de las células T CD4 infectadas de reproducir partículas de virus, y del opuesto, de activar las respuestas inmunes cuando están libres de infección. La respuesta inmune CTL durante la etapa asintomática de la infección y el comportamiento de la solución con las tres etapas de infección, las etapas asintomática y final del SIDA, se ven reflejadas en el modelo.

La implementación de modelación estocástica par atacar problemas de modelación en VIH conlleva problemas asociados al desarrollo de metodologías de simulación y de estimación paramétrica bayesiana. Ante la importancia de facilitar el uso de herramientas estadísticas dentro de estudios con diversas experiencias de área se proporciona una interfaz gráfica de usuario (GUI) mientras se interactúa localmente con \(\texttt{R}\) o directamente desde la aplicación web (https://randomverse.shinyapps.io/RV_int/). La librería asociada a la aplicación se encuentra alojada en github (https://github.com/jp-rgb/randomverse).

En el área de modelación estocástica para aplicaciones en biología de sistemas de VIH existen toda una serie de librerías en \(\texttt{R}\) para atacar distintas etapas del modelado. \(\texttt{smfsb}\) (Wilkinson 2018) es la referencia principa para combinar modelación estocástica, aplicaciones biológicas y desarrollo \(\texttt{R}\). Las librerías \(\texttt{adaptivetau}\) (Johnson 2019) y \(\texttt{GillespieSSA}\) (Pineda-Krch and Cannoodt 2022) implementan el algoritmo exacto de simulación de Gillespie (método directo) así como otros métodos de aproximación. \(\texttt{GillespieSSA2}\) (Cannoodt et al. 2020) y \(\texttt{sar}\) ((Insp-Rh, n.d.)) presentan mejoras sobre \(\texttt{GillespieSSA}\) con una implementación más veloz del algoritmo de Gillespie. Las librería \(\texttt{lnar}\) (Giagos 2010) utiliza inferencia basada en verosimilitud de modelos de kinética estocástica usando una Aproximación de Ruido Lineal.

El material antes mencionado se ha desarrollado con un enfoque, por lo tanto, la aplicación unifica en una herramienta estadística las tareas de simulación e inferencia sobre los parámetros de un modelo de ODEs/SDEs sin la necesidad de emplear rutinas computacionales intensivas.

1.1 Objetivos

La modelación estocástica y la inferencia de datos para el VIH se ha vuelto cada vez más compleja. Esto puede presentar desafíos para desarrollar aplicaciones biológicas y terapéuticas, por lo que se ha vuelto díficil obtener información sobre la dinámica viral sin asistencia de una metodología estadística apropiada para administrar y procesar estos conjuntos de datos. Una respuesta para abordar las demandas estadísticas asociadas con la modelación estocástica es formar equipos en transdisciplina o colaboraciones. Idealmente los conjuntos de habilidades complementarias cubren todas las necesidades potenciales en la elaboración de modelos para VIH. En la práctica, sin embargo, ciertos análisis estadísticos se llevan a cabo con tanta frecuencia que tiene sentido crear herramientas que permitan a todas las partes del equipo realizarlas fácilmente, a pesar del esfuerzo inicial requerido para configurar la infraestructura necesaria. Los objetivos específicos de este proyecto de tesis son los siguientes:

  • Diseñar y construir una aplicación web utilizando el lenguaje de programación \(\texttt{R}\) que facilite el análisis de análisis de dpósitos de latencia en VIH. Para tal fin, se creará una aplicación que se alimenta de los parámetros de un sistema de ODEs/SDEs y permite la construcción de trayetorias de un modelo con latencia desde una perspectiva determinista o estocástica y sus visualizaciones a través de una interfaz GUI. La documentación para crear la aplicación se alojará en GitHub, lo que permitirá su fácil descargar para trabajar con el código.
  • Escribir documentación de “Guía del usuario” que detalle las funciones y el uso de la aplicación descrita en el punto anterior. Esta guía también incluirá información sobre la configuración de sistemas necesaria para ejecutar las aplicaciones.
  • Implementar el modelo de latencia del VIH desde el enfoque determinista empleando el método de Euler y desde el enfoque estocástico con el algoritmo de Gillespie, a través de las funcionalidades de la aplicación.

1.2 Alcance

La aplicación puede ayudar a profesionales, estudiantes, investigadores y académicos que se enfrenten a problemas de modelación estocástica para estudiar los depósitos de latencia del VIH y no cuenten con la suficiente infrastructura estadística para analizar dichos procesos

La documentación en la librería en \(\texttt{R}\) permite contar con casos de estudio, guías del ususario y manuales de referencia para poder explotar las capacidades de modelación estocástica en la aplicación web.

En general la aplicacion es una herramienta de libre acceso e intuitiva para estudiar modelos estocásticos para el VIH.

2 Metodología

2.1 Representación analítica

El modelo sigue el papel de las células T CD4 en tres etapas durante la infección por VIH: (1) aguda, (2) asintomática y (3) SIDA. (1) La primera etapa de la infección tiene células T CD4 infectadas y la carga viral aumenta exponencialmente. Las células T CD4 disminuyen significativamente, pero se activa la disminución del nivel de replicación del virus. Además, se produce una gran cantidad de anticuerpos que disminuye drásticamente el número de viriones en la sangre. Los síntomas duran menos de 14 días en promedio. (2) La segunda etapa es la etapa asintomática. Esta etapa es prolongada, dura de 6 a 15 años y se caracteriza por una disminución lenta de las células T CD4. En esta etapa se suelen mostrar pocos o ningún síntoma. En esta etapa, la carga viral aumenta y las células T CD4 disminuyen con el tiempo a un ritmo relativamente constante. (3) La tercera y última etapa es la etapa del SIDA. La formación de anticuerpos cae drásticamente. Eventualmente, las células T CD4 disminuyen a un nivel tan bajo que ya no pueden funcionar correctamente para activar la respuesta inmunitaria.

A través de un sistema simple de ODE se representan los roles opuestos de las células T CD4. Al descartar la respuesta inmunitaria, el modelo se simplifica a un modelo básico para el VIH-SIDA. El modelo básico incluye células objetivo sanas, células objetivo infectadas y virus o viriones libres. Varios modelos incluyen los efectos de la respuesta inmune en este modelo básico. La respuesta CTL se introduce como un sistema de ODEs de cuatro dimensiones, pues los modelos recientes de VIH se centran en diferentes aspectos de la respuesta de CTL como respuesta adaptativa del huésped a un patógeno.

En el modelo la estimulación de los CTL se produce en presencia de células T CD4 y partículas de virus infecciosos y no se debe directamente a la presencia de células infectadas ya que el genoma viral del VIH es capaz de integrarse en las células T CD4 en reposo para producir una infección latente. El virus latente que reside dentro de la célula en reposo está protegido de la respuesta inmunitaria, por lo que el tratamiento con terapia antirretroviral es ineficaz. El reservorio latente es una barrera importante para la erradicación del VIH y, por lo tanto, los agentes que revierten la latencia son un área importante de investigación para su tratamiento. El modelo básico con cuatro ODE se amplía al incluir el depósito o reservorio de latencia. Otros modelos similares se han propuesto con un reservorio latente pero sin una respuesta CTL para investigar el papel de múltiples tratamientos con fármacos antirretrovirales. También se han propuesto y estudiado modelos con un reservorio latente pero sin una respuesta CTL y modelos de infección con un reservorio latente en combinación con terapia antirretroviral.

El modelo para la dinámica del virus consta de un sistema de cinco ODE, incluidas las variables \(C\), \(I\), \(F\) y \(V\) que representan la densidad de células T CD4 no infectadas o sanas, células T CD4 infectadas, células CTL y partículas de virus libres, respectivamente, al tiempo \(t\). Una característica importante de la infección por VIH es el desarrollo de un reservorio latente (LR) de infección. Para su reproducción, el VIH integra su genoma en el ADN de células T CD4 sanas. Durante la infección algunas de las células T CD4 infectadas se convierten en células en reposo como respuesta fisiológica natural de la respuesta inmunitaria adaptativa del huésped a la infección. El genoma del virus permanece integrado hasta que la célula huésped se reactiva ante infecciones pasando al estado \(I\) y comenzando a replicar viriones. La vida media de LR es de unos 44 meses. Así, al integrar la latencia al modelo se consideran dos tipos de células infectadas, ya sea el reservorio latente \(L\) que no produce viriones pero puede infectarse activamente y las células infectadas que se reproducen activamente \(I\) . Una fracción \(η\) de células T CD4 se infecta de forma latente y la fracción restante \((1 − η)\) se infecta activamente, capaz de producir viriones. El parámetro \(d_0\) es la tasa de muerte de las células latentes y \(a_L\) es la tasa de transición de células con infección latente a células con infección activa.

El modelo ODE tiene la siguiente forma:

\[\begin{equation} \begin{aligned} \frac{dC}{dt} &= \lambda(C^*-C)-\beta CV \\ \frac{dL}{dt} &= \eta \beta CV - d_0L - a_LL \\ \frac{dI}{dt} &= (1-\eta ) \beta CV - aI + a_LL - \rho FI \\ \frac{dC}{dt} &= \frac{eCVF}{C^{*}+F} - bF \\ \frac{dV}{dt} &= aNI-kV. \end{aligned} \tag{1} \end{equation}\]

Las condiciones iniciales del sistema satisfacen \(0 < C (0) < C^∗\) , \(V (0) > 0, I (0) ≥ 0, F (0) > 0\) y \(C ( 0 ) + I ( 0 ) < C^∗\) . La constante \(λ C^∗\) es la tasa de crecimiento de las células T CD4 no infectadas. El parámetro \(β\) es la velocidad a la que las células sanas se infectan con una partícula de virus. Los parámetros \(λ, a, b\) y \(k\) son las tasas de muerte de las células T CD4 no infectadas, las células T CD4 infectadas, los CTL y las partículas de virus. El parámetro \(ρ\) es la velocidad a la que las células CTL eliminan del cuerpo las células infectadas. Además se deben controlar dos parámetros importantes previo a la estimación del sistema. El primero es el parámetro \(\eta\). El segundo es \(e\), la tasa de estimulación de los CTL. La expresión \(eCVF\) por sí sola proporciona una tasa poco realista de inducción de la respuesta inmunitaria para un número insignificante de partículas virales (\(V\) pequeña). La expresión del denominador \(C^∗ + F\) asume que la respuesta de CTL se satura a niveles altos de \(F\) y, además, el término constante \(C^∗\) limita efectivamente la cantidad de estimulación del virus en CTL, ya que la producción de virus proviene de las células \(I\) que se reproducen activamente.

La reproducción viral por parte de las células T CD4 infectadas genera \(N\) nuevas partículas virales producidas por cada célula infectada durante su vida (el promedio de vida es \(1/a\)). Por lo tanto, \(aN\) es la tasa de producción de virus por parte de una célula infectada. Generalmente, se da el caso de que la tasa de muerte de las células infectadas es mayor que la de las células sanas. Así, se hace la suposición \(a ≥ \lambda\).

Se debe contar con \(0 < η < 1\) y \(0 < a_L < ∞\) . El papel de la estimulación de las células T CD4 sobre el reservorio latente, \(L\), se da a través del parámetro \(e\).

2.2 Enfoque determinista

Cuando los sistemas en biología plantean ODEs que resultan complicadas para atacar analíticamente, las soluciones pueden ser examinadas integradndo numéricamente el ODE. El enfoque más simple y práctico es el Método de Euler de primer orden.

La idea básica de integrar numéricamente un sistema de ODEs es tomar la representación analítica dada como un vector de ODE como en la Equación (2).

\[\begin{equation} \frac{dX}{dt} = f(X), \tag{2} \end{equation}\]

donde \(X\) es un vector \(p\)-dimensional y \(f (·) : R^p → R^p\) es una función \(p-\)dimensional arbitraria (no-lineal) de \(X\).

Se emplea la definición de derivada definida en Ecuación (3)

\[\begin{equation} \frac{dX}{dt}(t) = \lim_{\Delta t \to 0} \frac{X(t + \Delta t)- X(t)}{\Delta t}, \tag{3} \end{equation}\]

con un pequeño incremento \(∆t\) Ecuación (4) es obtenida

\[\begin{equation} \frac{X(t+\Delta t)-X(t)}{\Delta t} \approx f(X(t)) \tag{4} \end{equation}\]

y reacomodando, Ecuación (5) da un método simple para calcular \(X(t + ∆t)\) a partir de \(X(t)\). Tomando como condición inicial \(X(0)\), entonces \(X(∆t), X(2∆t), X(3∆t), . . .\) se pueden calcular para obtener la dinámica completa del sistema.

\[\begin{equation} X(t + ∆t) \approx X(t) + ∆tf (X(t)). \tag{5} \end{equation}\]

La versión en shiny para el simulador determinista (en función del horizonte de tiempo) se encuentra en la pestaña \(\texttt{Deterministic}\) y su algoritmo se puede resumir como sigue:

  1. Iniciar sistema con condiciones inciciales \(t = 0\), tasas de reacción iniciales \(f(X(t))\) y números de moléculas para cada especie, \(X= (x_1 , x _2 , . . . , x_u)\).
  2. Un \(\Delta t\) de tiempo es fijada
  3. El tiempo al siguiente evento es calculado como \(t : = t + \Delta t\).
  4. EL siguiente evento es calculado como \(X = X + f(X(t)) \Delta t\).
  5. Los valores actualizados para \(X\) y \(t\) se apilan en forma de lista.
  6. Si \(t < T_{max}\), el algoritmo regrea al paso \(2\).
  7. La salida, tanto de \(X\) como de \(t\) se encuadra como un arreglo de datos una vez que el bucle para.

Para un modelo con carga viral alta se tienen valores de \(e=10\) y \(\eta=0.5\). Las realizaciones para cada especie de la Figura 1 muestran el avance de las partículas del virus libres y como comienza a descender rápidamente la cantidad de células T CD4 sanas e iniciar con la infección.

Una realizazción de simulación determinista de la dinámica de la respueta inmune del VIH por un periodo de 200 días. Parámetros e=10 y eta=.5

Figure 1: Una realizazción de simulación determinista de la dinámica de la respueta inmune del VIH por un periodo de 200 días
Parámetros e=10 y eta=.5

2.3 Simulación estocástica

Para llevar a cabo el estudio de simulción se toma un sistema con \(u\) species \(X_1 , X_2 , . . . , X_u\) y \(v\) reacciones \(R_1 , R_2 , . . . , R_v\). El número de moléculas de \(X_i\) al tiempo \(t\) es \(X_{it}\), y el estado del sistema al tiempo \(t\) is \(X_t = (X_{1t} , X_{2t} , . . . , X_{ut})^T\). El númerod e reacciones del tipo \(R_i\) en la ventana de tiempo \((0, t]\) es \(R_{it}\), y así \(R_t = (R_{1t} , . . . , R_{vt} )^T\). Entonces la Ecuación (6) actualiza el estado del sistema.

\[\begin{equation} X_t − X_0 = SR_t. \tag{6} \end{equation}\]

Aquí,

\[\begin{equation} S = (Post − Pre)^T, \tag{7} \end{equation}\]

es la \(u × v\) matriz de estequiometrías de la red del modelo. Cada constante estocástica para las tasas de las reacciones \(R_i\) se denota con \(c_i\) y su tasa o ley de reacción asociada es \(h_i (x, c_i )\), donde \(x = (x_1 , x_2 , . . . , x_u )^T\) es el estado actual del sistema. Cuando el estado del sistema es \(x\) atl tiempo \(t\), la probabilidad de que en \(R_i\) reacciones (o transiciones) ocurran en el intervalo \((t, t + dt]\) viene dado por \(h_i (x, c_i ) dt\). En la ausencia de cualesquiera otras reacciones llevándose a cabo, el tiempo para cada evento de reacción es una cantidad aleatoria exponencial denotada con \(Exp(h_i (x, c_i ))\).

Como el sistema consiste de \(v\) reacciones y la tasa de reacción del tipo de reacción \(i\) es \(h_i (x, c_i )\), la tasa de reacción para que ocurra algún tipo de reacción sigue la Ecuatión (8).

\[\begin{equation} h_0 (x, c) ≡ \sum_{i=1}^v h_i (x, c_i ). \tag{8} \end{equation}\]

Un evento de reacción \(Exp(h_0 (x, c))\) que actualiza el proceso ocurre de forma aleatoria proporcional a \(h_i (x, c_i )\) (la reacción será de tipo \(i\) con probabilidad \(h_i (x, c_i )/h_0 (x, c)\)) y es independiente al tiempo del siguiente evento. La simulación progresa usando y el tiempo al siguiente evento y el tipo de evento se actualizan. Este procedimiento de simulación de eventos discretos estándar se conoce como el algoritmo de Gillespie.

La versión en \(\texttt{shiny}\) de este algoritmo está implementada bajo la pestaña \(\texttt{Stochastic}\) y se resume como sigue:

  1. El sistema empieza con condiciones iniciales \(t = 0\), constantes de las tasas de reacción \(c_1 , c _2 , . . . , c_v\) y números inciales de moléculas para cada especie, \(x_1 , x _2 , . . . , x_u\) .
  2. Una función de riesgo \(h_i (x, c_i )\) es calculada para cada \(i = 1, 2, . . . , v\) basada en el estado actual \(x\).
  3. La función de riesgo combinada \(h_0 (x, c) ≡ \sum_{i=1}^v h_i (x, c_i )\) es obtenida.
  4. Se usa una muestra de una cantidad aleatoria \(Exp(h_0 (x, c))\) para calcular el tiempo alsiguiente evento \(t'\).
  5. El timepo se actualiza como \(t : = t + t'\) .
  6. Se usa una muestra de una cantidad aleatoria discreta con probabilidades \(h_i (x, c_i ) / h_0 (x, c), i = 1, 2, . . . , v\) para calcular el índide del tipo \(j\) de la reacción al siguiente evento.
  7. EL estado de \(x\) se actualiza de acuerdo a la reacción \(j\) como \(x : = x + S (j)\) donde \(S(j)\) denota la \(j\)th columna de la matriz de estoquiometrías \(S\).
  8. Los valores actualizados de \(x\) y \(t\) se apilan en la forma de una lista.
  9. Si \(t < T_{max}\), el algoritmo regresa al paso \(2\).
  10. La salida de \(x\) and \(t\) se enmarca en un erreglo de datos al terminar el bucle.

Session info

## R version 4.2.2 Patched (2022-11-10 r83330)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=es_MX.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_MX.UTF-8        LC_COLLATE=es_MX.UTF-8    
##  [5] LC_MONETARY=es_MX.UTF-8    LC_MESSAGES=es_MX.UTF-8   
##  [7] LC_PAPER=es_MX.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_MX.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] BiocStyle_2.24.0
## 
## loaded via a namespace (and not attached):
##  [1] bookdown_0.29       digest_0.6.29       R6_2.5.1           
##  [4] jsonlite_1.8.0      magrittr_2.0.3      evaluate_0.16      
##  [7] highr_0.9           stringi_1.7.8       cachem_1.0.6       
## [10] rlang_1.0.6         cli_3.4.0           rstudioapi_0.14    
## [13] jquerylib_0.1.4     bslib_0.4.0         rmarkdown_2.16     
## [16] tools_4.2.2         stringr_1.4.0       xfun_0.33          
## [19] yaml_2.3.5          fastmap_1.1.0       compiler_4.2.2     
## [22] BiocManager_1.30.18 htmltools_0.5.4     knitr_1.40         
## [25] sass_0.4.2

References

Aavani, Pooya, and Linda J. S. Allen. 2019. “The Role of CD4 t Cells in Immune System Activation and Viral Reproduction in a Simple Model for HIV Infection.” Applied Mathematical Modelling 75: 210–22. https://doi.org/10.1016/j.apm.2019.05.028.
Cannoodt, Robrecht, Wouter Saelens, Louise Deconinck, and Yvan Saeys. 2020. “Dyngen: A Multi-Modal Simulator for Spearheading New Single-Cell Omics Analyses.” bioRxiv, February. https://doi.org/10.1101/2020.02.06.936971.
Chang, Winston, Joe Cheng, JJ Allaire, Carson Sievert, Barret Schloerke, Yihui Xie, Jeff Allen, Jonathan McPherson, Alan Dipert, and Barbara Borges. 2022. Shiny: Web Application Framework for r. https://CRAN.R-project.org/package=shiny.
Giagos, Vasileios. 2010. “Inference for Auto-Regulatory Genetic Networks Using Diffusion Process Approximations.” PhD thesis.
Gillespie, Daniel T. 1977. “Exact Stochastic Simulation of Coupled Chemical Reactions.” The Journal of Physical Chemistry 81 (25): 2340–61. https://doi.org/10.1021/j100540a008.
Insp-Rh. n.d. “INSP-RH/Ssar: Fast r Implementation of Gillespie’s Stochastic Simulation Algorithm.” GitHub. https://github.com/INSP-RH/ssar.
Johnson, Philip. 2019. Adaptivetau: Tau-Leaping Stochastic Simulation. https://CRAN.R-project.org/package=adaptivetau.
Pineda-Krch, Mario, and Robrecht Cannoodt. 2022. GillespieSSA: Gillespie’s Stochastic Simulation Algorithm (SSA). https://CRAN.R-project.org/package=GillespieSSA.
Wilkinson, Darren. 2018. Smfsb: Stochastic Modelling for Systems Biology. https://CRAN.R-project.org/package=smfsb.